C EXDIAR SOURCE CB215821 22/07/20 15:39:40 11411 SUBROUTINE EXDIAR(MRIGID,MCHPOI) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C************************************************************************* C Operateur EXDIAR C C OBJET : Extrait la diagonale d'une matrice. C C*********************************************************************** C HISTORIQUE : 06/06/1 : Première version C C*********************************************************************** -INC PPARAM -INC CCOPTIO -INC CCHAMP -INC SMCHPOI -INC SMRIGID -INC SMELEME -INC SMCOORD -INC TMTRAV * segment trav1 character*4 compd(nbincd),compdd(nbincd) endsegment segment trav2 integer cor1p(nligrp),cor1d(nligrd) endsegment * * Pour chaque segment DESCR, on a un segment CORES2 où * CORDIA donne les couples (NLIGRP,NLIGRD) qui sont des * termes diagonaux avec les cas particuliers suivant : * CORES2=0 : il n'y a pas de termes diagonaux dans ces matrices * élémentaires * CORES2=-1 : les termes diagonaux de la matrices élémentaire (carrée) * sont des termes diagonaux de la matrice globale * on a aussi a segment CORES3 où COMTRA indique * le nom d'inconnue primal+numéro d'harmonique suivant * le repérage induit par MTRA1 et MTRA2 * SEGMENT CORES1 INTEGER ITRAV2(NRIGEL) INTEGER IPCOR2(NRIGEL) INTEGER IPCOR3(NRIGEL) ENDSEGMENT SEGMENT CORES2 INTEGER CORDIA(2,NLIGP) ENDSEGMENT SEGMENT CORES3 INTEGER COMTRA(NLIGP) ENDSEGMENT * SEGMENT ICPR(nbpts) SEGMENT MTRA1 * CHARACTER*4 ICOMP(0) INTEGER ICOMP(0) ENDSEGMENT SEGMENT MTRA2 INTEGER MHAR(0) ENDSEGMENT logical descar * * Construisons d'abord l'ensemble des noms d'inconnues * primales, puis les duales associées (segment trav1) * C WRITE(IOIMP,*) 'Coucou exdiar' SEGACT MRIGID nbincd=0 segini trav1 NRIGEL=COERIG(/1) DO IRIG=1,NRIGEL DESCR=IRIGEL(3,IRIG) SEGACT DESCR do 45 iligrp=1,lisinc(/2) do 55 inc=1,nbincd if (compd(inc).eq.lisinc(iligrp)) goto 45 55 continue nbincd=nbincd+1 segadj trav1 compd(nbincd)=lisinc(iligrp) 45 continue SEGDES DESCR ENDDO do iincd=1,nbincd CALL PLACE(NOMDD,LNOMDD,idx,compd(iincd)) IF (idx.EQ.0) THEN compdd(iincd)=compd(iincd) ELSE compdd(iincd)=NOMDU(idx) ENDIF ENDDO C WRITE(IOIMP,*) ' 1) compd' C WRITE (IOIMP,2019) (compd(I),I=1,compd(/2)) C WRITE(IOIMP,*) ' 2) compdd' C WRITE (IOIMP,2019) (compdd(I),I=1,compdd(/2)) * * Maintenant, on construit les CORES1 et CORES2 qui nous * permettent de savoir quels sont les termes diagonaux * dans les sous-matrices * SEGINI CORES1 DO IRIG=1,NRIGEL DESCR=IRIGEL(3,IRIG) SEGACT DESCR * Tout d'abord le segment trav2 qui permet d'avoir * la correspondance sur les noms d'inconnues NLIGRP=LISINC(/2) NLIGRD=LISDUA(/2) segini trav2 do 410 iligrp=1,nligrp do 420 incd=1,nbincd if (compd(incd).eq.lisinc(iligrp)) then cor1p(iligrp)=incd goto 410 endif 420 continue 410 continue do 510 iligrd=1,nligrd do 520 incd=1,nbincd if (compdd(incd).eq.lisdua(iligrd)) then cor1d(iligrd)=incd goto 510 endif 520 continue * Si on a pas trouvé dans la liste des duales, on regarde * dans la liste des primales. * Ceci est pratique en mécaflu où il arrive qu'on nomme les duales * pareil que les primales (ex. 'UX' 'UX') do 525 incd=1,nbincd if (compd(incd).eq.lisdua(iligrd)) then cor1d(iligrd)=incd goto 510 endif 525 continue 510 continue * Voyons si le descripteur est carré, ce qui est un * cas courant qui simplifie notablement les choses descar=.false. if (nligrp.ne.nligrd) goto 600 do iligr=1,nligrp if (cor1p(iligr).ne.cor1d(iligr)) goto 600 if (noelep(iligr).ne.noeled(iligr)) goto 600 enddo descar=.true. 600 continue * Si oui, tout est simple, sinon il faut rechercher * une correspondance primale-duale if (descar) then ipcor2(irig)=-1 else nligp=nligrp segini cores2 nligp=0 do iligrp=1,nligrp do iligrd=1,nligrd if (cor1p(iligrp).eq.cor1d(iligrd).and. $ noelep(iligrp).eq.noeled(iligrd)) then nligp=nligp+1 cordia(1,nligp)=iligrp cordia(2,nligp)=iligrd endif enddo enddo if (nligp.eq.0) then segsup cores2 ipcor2(irig)=0 else segadj cores2 ipcor2(irig)=cores2 endif endif itrav2(irig)=trav2 C Write(ioimp,*) 'irig=',irig C Write(ioimp,*) 'descar=',descar C if (ipcor2.gt.0) then C write(ioimp,*) 'cordia(1, )' C WRITE (IOIMP,2020) (cordia(1,I),I=1,cordia(/2)) C write(ioimp,*) 'cordia(2, )' C WRITE (IOIMP,2020) (cordia(2,I),I=1,cordia(/2)) C else C WRITE(IOIMP,*) 'ipcor2=',ipcor2 C endif SEGDES DESCR ENDDO * * Maintenant, on peut construire MTRA1, MTRA2 * (repérage des couples noms d'inconnues primales, * harmoniques) ainsi que CORES3 * et ICPR (repérage global->local des noeuds) * ce qui permettra dans une deuxième passe de remplir * le segment TMTRAV qui sera transformé en CHPOINT * SEGINI ICPR NNNOE=0 SEGINI MTRA1,MTRA2 NNIN=0 DO 100 IRIG=1,NRIGEL DESCR=IRIGEL(3,IRIG) trav2=itrav2(irig) MELEME=IRIGEL(1,IRIG) cores2=ipcor2(irig) if (cores2.eq.0) goto 99 SEGACT DESCR SEGACT MELEME if (cores2.eq.-1) then nligp=noelep(/1) else nligp=cordia(/2) endif segini cores3 do 110 iligp=1,nligp if (cores2.eq.-1) then iligrp=iligp else iligrp=cordia(1,iligp) endif do ic=1,icomp(/1) if (icomp(ic).eq.cor1p(iligrp).and. $ mhar(ic).eq.irigel(5,irig)) then comtra(iligp)=ic goto 110 endif enddo icomp(**)=cor1p(iligrp) mhar(**)=irigel(5,irig) nnin=nnin+1 comtra(iligp)=nnin 110 continue C WRITE(IOIMP,*) ' 3) comtra' C WRITE (IOIMP,2020) (comtra(I),I=1,comtra(/1)) ipcor3(irig)=cores3 * Dans le cas cores2=-1 * tous les noeuds du meleme sont référencés * sinon, on regarde ceux qui le sont if (cores2.eq.-1) then nno=num(/1) else nno=nligp endif do iel=1,num(/2) do ino=1,nno if (cores2.eq.-1) then ilno=ino else ilno=noelep(cordia(1,ino)) endif ipt=num(ilno,iel) if (icpr(ipt).eq.0) then nnnoe=nnnoe+1 icpr(ipt)=nnnoe endif enddo enddo 99 continue SEGDES MELEME SEGDES DESCR segsup trav2 100 CONTINUE C WRITE(IOIMP,*) 'nnin=',nnin C WRITE(IOIMP,*) 'nnin2=',icomp(/1) C WRITE(IOIMP,*) 'nnin3=',mhar(/1) C WRITE(IOIMP,*) ' 1) icomp' C WRITE (IOIMP,2020) (icomp(I),I=1,icomp(/1)) C WRITE(IOIMP,*) ' 2) mhar' C WRITE (IOIMP,2020) (mhar(I),I=1,mhar(/1)) C WRITE(IOIMP,*) 'nnnoe=',nnnoe C WRITE(IOIMP,*) ' 4) icpr' C WRITE (IOIMP,2020) (icpr(I),I=1,icpr(/1)) * * Maintenant, on peut remplir MTRAV * SEGINI MTRAV DO ININ=1,NNIN INCO(ININ)=compd(ICOMP(ININ)) NHAR(ININ)=MHAR(ININ) ENDDO SEGSUP MTRA1,MTRA2 segsup trav1 C xt=0.D0 DO 1100 IRIG=1,NRIGEL MELEME=IRIGEL(1,IRIG) DESCR=IRIGEL(3,IRIG) XMATRI=IRIGEL(4,IRIG) cores2=ipcor2(irig) cores3=ipcor3(irig) if (cores2.eq.0) goto 1099 SEGACT MELEME SEGACT DESCR SEGACT XMATRI nel=num(/2) if (cores2.eq.-1) then nligp=noelep(/1) else nligp=cordia(/2) endif do iel=1,nel do iligp=1,nligp if (cores2.eq.-1) then iligrp=iligp iligrd=iligrp else iligrp=cordia(1,iligp) iligrd=cordia(2,iligp) endif ININ=comtra(iligp) ino=noelep(iligrp) ipt=num(ino,iel) INNOE=icpr(ipt) C if (ipt.eq.2.and.lisinc(iligrp).eq.'UX ') then C write(ioimp,*) 'irig=',irig C write(ioimp,*) 'iel=',iel C write(ioimp,*) 'iligrp=',iligrp C write(ioimp,*) 'iligrd=',iligrd C xt=xt+RE(iligrd,iligrp,iel) C WRITE(IOIMP,*) 'RE=',RE(iligrd,iligrp,iel) C WRITE(IOIMP,*) 'xt=',xt C WRITE(IOIMP,*) 'inin=',inin C WRITE(IOIMP,*) 'innoe=',innoe C endif IGEO(INNOE)=ipt IBIN(ININ,INNOE)=1 BB(ININ,INNOE)=BB(ININ,INNOE)+ $ RE(iligrd,iligrp,iel)*COERIG(IRIG) enddo enddo if (cores2.ne.-1) segsup cores2 segsup cores3 SEGDES XMATRI SEGDES DESCR SEGDES MELEME 1099 continue 1100 CONTINUE SEGSUP ICPR SEGSUP CORES1 SEGDES MRIGID * * Création du CHPOINT résultat * CALL CRECHP(MTRAV,MCHPOI) $ SEGSUP MTRAV RETURN 2019 FORMAT (20(2X,A4) ) 2020 FORMAT (20(2X,I4) ) 2021 FORMAT (20(2X,L4) ) 2022 FORMAT(10(1X,1PG12.5)) * * End of EXDIAG * END