C MONDES SOURCE MB234859 26/06/10 21:15:48 12569 SUBROUTINE MONDES(MMATRX,MVECTX,NOEN,ISOUCI,lagdua) C C **** EXECUTE LA SOLUTION X DE (Lt D L) X=F C CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB CMB CMB Plutot la solution de L.D.Lt ou L.D.Mt (cas non symétrique) CMB Elle devrait dons s'appeller DESMON et non MONDES. CMB CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (LPCRAY=10000) INTEGER OOOVAL C ... Variable décrivant l'EXIStence du Triangle Supérieur différent C de l'inférieur (cas non symétrique) ... LOGICAL EXISTS -INC SMMATRI -INC SMELEME -INC SMVECTD POINTEUR MVECT4.MVECTD -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SILICRE -INC CCASSIS -INC TMTRAV SEGMENT ITTRV(MULRE) SEGMENT,ITRAA(NENSL) SEGMENT/BID/(IBIDON(IIMAX)) logical toutmem mvect4 = 0 mtrav=0 toutmem=.true. DNORMA=0.D0 DNORMB=0.D0 crit=xpetit crit2 = 0.D0 call oooprl(1) nbthr=nbthrs CALL OOOMRU(1) IVALMA=0 IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX 1000 FORMAT(' SUBROUTINE MONDES : POINTEUR DE LA MATRICE=',I5, 1 ' POINTEUR DU VECTEUR=',I5) C C **** ACTIVATION DES SEGMENTS C MMATRI=MMATRX SEGACT,MMATRI*MOD * BEC00SE OPTIMISEUR ITRAA=MMATRI NENSL=0 IF(NENS.NE.0) THEN NENSL=NENS SEGINI,ITRAA ENDIF NE1=NENSL EXISTS=.FALSE. IF(IILIGS.NE.0) THEN MILIGN=IILIGS SEGACT MILIGN IF(ILIGN(/1).GT.0) EXISTS=.TRUE. ENDIF MILIGN=IILIGN SEGACT,MILIGN INC=IPNO(/1) nbthr=min(nbthrs,inc/1000+1) MVECTD=MVECTX SEGACT MVECTD*MOD MDNOR=IDNORM SEGACT MDNOR IF(IDNORD.GT.0) THEN MDNO1=IDNORD SEGACT MDNO1 ELSE MDNO1=MDNOR ENDIF MDIAG=IDIAG C ... MULRE = nombre de seconds membres ... MULRE = VECTBB(/1)/INC C ... MUINC ne servira que comme une borne sur les indices de boucles ... MUINC = ( MULRE-1)*INC SEGINI ITTRV inc=vectbb(/1) * mvect1 sauve les forces * mvect2 sauve le deplacement * mvect3 vecteur de travail * mvectd evolue: force a resoudre (residu) puis increment de deplacement if (mulre.eq.1.and..not.exists) segini mvect1,mvect2 INC=IPNO(/1) C ... Multiplication du second membre par les DNOR(.) (à cause du C conditionnement des matrices de blocages) ... DNORMA=0.D0 DO 450 K=0,MUINC,INC DO 45 I=1,INC VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I) ** DNORMA=MAX(DNORMA,ABS(VECTBB(I+K))) 45 CONTINUE 450 CONTINUE ** DNORMA=DNORMA*xzprec*xzprec DNORMA=xpetit/xzprec C ... Détection du premier terme important du second membre ... II=0 DO 451 K=0,MUINC,INC II=II+1 DO 452 I=1,INC IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452 C ... Le numéro du noeud concerné par le premier terme important va à ITTRV ... ITTRV(II)=IPNO(I) GO TO 451 452 CONTINUE 451 CONTINUE C ... NNOE = nombre de noeuds concernés ... NNOE=ILIGN(/1) IDEP=NNOE C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ... DO 453 I=1,MULRE IF(ITTRV(I).LT.IDEP) THEN IF(ITTRV(I).NE.0) THEN IDEP=ITTRV(I) ELSE IDEP=1 GO TO 4530 ENDIF ENDIF 453 CONTINUE 4530 CONTINUE SEGDES MDNOR * * debut de la boucle d'amelioration de la precision si necessaire * * sauver le champs de forces initial if(mulre.eq.1.and..not.exists) then do i=1,inc mvect1.vectbb(i)=vectbb(i) enddo icorrec=0 endif 10000 continue IAA=0 SEGACT,MDIAG C C **** DESCENTE: ON RESOU L*C=B. EN FAIT ON STOCKE C DANS B C LTRK=MAX(1+LPCRAY,OOOVAL(1,4)) IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK iimax=oooval(1,1)/10 * test si la place disponible permet de tout mettre en memoire xplds=oooval(1,1)-oooval(2,3) if (real(nnoe)*ijmax.lt.xplds) iimax=0 SEGINI /err=3500/ BID goto 3501 3500 continue * place insuffisante en memoire centrale. On en fait * desactivation de la matrice en partant du debut do ivs=1,nnoe lign=ilign(ivs) segdes lign enddo segini bid 3501 continue C ... IDEP : on commence par ce noeud car le second membre est nul pour C tous les precedents ... * activation de tous les derniers noeuds car ils sont deja en memoire normalement do 2001 ivs=nnoe,1,-1 lign=ilign(ivs) segact/err=2010/ lign 2001 continue 2010 continue ilmin=ivs+1 DO 610 IVS=IDEP,NNOE LIGN=ILIGN(IVS) SEGACT /ERR=50/ LIGN IVALMA=IVALMA+VAL(/1) *pvpv IF (IVALMA.GT.NGMAXY) GOTO 50 NA=IMMM(/1) IPRELL=IPREL DO 611 IA=1,NA C ... Si l'inconnue présente un mode d'ensemble ... IF(IMMM(IA).NE.0) THEN C ... On incrémente le compteur et IAA=IAA+1 C ... on met le N° de la ligne dans ITRAA ... ITRAA(IAA)=IPRELL+IA-1 ENDIF 611 CONTINUE segact lign CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1), > NA,IPREL,MULRE,INC,IVS,LLOM,dnorma) 610 CONTINUE C ... Si on n'a pas eu de pb de mémoire on passe par là ... SEGSUP BID C ... ILMAX vaut le dernier noeud qui a pu être traité ... ILMAX=NNOE C ... On va à la division par la diagonale ... GOTO 55 50 CONTINUE C ... Si on est là, c'est à cause des pb de mémoire ... SEGSUP BID toutmem=.false. C ... IVS est le N° du LIGN qui n'a pas pu être traité ... C ... ILMAX = N° du dernier traité ... ILMAX=IVS-1 C ... IILMAX = N° du premier à traiter ... IILMAX=IVS DO 54 IVS=IILMAX,NNOE 58 CONTINUE LIGN=ILIGN(IVS) C ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ... SEGACT /ERR=56/ LIGN GOTO 57 56 CONTINUE * EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE toutmem=.false. C ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX, LIGN=ILIGN(ILMAX) if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU) C ... puis ce dernier est décrémenté ... ILMAX=ILMAX-1 IF (ILMAX.EQ.1) CALL ERREUR(5) C ... et on recommence ... GOTO 58 C ... On est là si tout s'est bien passé avec SEGACT ... 57 CONTINUE NA=IMMM(/1) IPRELL=IPREL DO 612 IA=1,NA IF(IMMM(IA).EQ.0) GOTO 612 IAA=IAA+1 ITRAA(IAA)=IPRELL+IA-1 612 CONTINUE CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1), > NA,IPREL,MULRE,INC,IVS,LLOM,dnorma) if (ivs.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU) 54 CONTINUE 55 CONTINUE C C ... À cet endroit la descente est finie. L'état des LIGN est le suivant : C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ... C C **** DIVISION PAR LE TERME DIAGONAL **** C dnorma=0.d0 idenorm=0 DO 120 K=0,MUINC,INC DO 12 I=1,INC J=I+K if (2*abs(diag(i)).gt.abs(diag(i))) goto 122 idenorm=1 122 continue VECTBB(J)=VECTBB(J)*DIAG(I) ** dnorma=max(dnorma,abs(vectbb(j))) 12 CONTINUE 120 CONTINUE ** dnorma=dnorma*xzprec*xzprec dnorma=xpetit/xzprec SEGDES,MDIAG C C **** MONTEE **** C C ... Si la matrice est asymétrique ... IF(EXISTS) THEN C ... On commence par désactiver les LIGN qui sont restés actifs ... DO 2000 I=IDEP,ILMAX LIGN=ILIGN(I) if (i.lt.ilmin) seGDES,LIGN*(NOMOD,MRU) 2000 CONTINUE C ... Puis on désactive MILIGN ... SEGDES,MILIGN C ... On change de MILIGN ... MILIGN=IILIGS SEGACT,MILIGN C ... Et on active des LIGN pointés par ce dernier ... ILMAX=IDEP-1 DO 2100 I=IDEP,NNOE LIGN=ILIGN(I) SEGACT /ERR=2110/ LIGN ILMAX = I 2100 CONTINUE C ... Et on passe à la montée ... GOTO 3000 C ... Si on n'a même pas activé le premier, alors ADIOS ... 2110 IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5) toutmem=.false. 3000 CONTINUE ENDIF J=NNOE C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ... 70 CONTINUE IF (J.EQ.ILMAX) GOTO 72 LIGN=ILIGN(J) SEGACT /ERR=68/ LIGN GO TO 67 C ... Si on a des pb avec activation, on désactive et on diminue ILMAX ... 68 CONTINUE toutmem=.false. LIGN = ILIGN( ILMAX) if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU) ILMAX=ILMAX-1 IF(ILMAX.EQ.1) CALL ERREUR (5) C ... puis on recommence la tentative ... GO TO 70 C ... On a réussi à activer ... 67 CONTINUE NA=IMMM(/1) IPRELL=IPREL CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), > NA,IPREL,MULRE,INC,dnorma) if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU) J = J-1 GO TO 70 72 CONTINUE CC FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!! CC FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1 C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ... DO 13 J=ILMAX,IDEP,-1 LIGN=ILIGN(J) NA=IMMM(/1) IPRELL=IPREL CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), > NA,IPREL,MULRE,INC,dnorma) if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU) 13 CONTINUE C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ... DO 113 J=IDEP-1,1,-1 LIGN=ILIGN(J) SEGACT LIGN NA=IMMM(/1) IPRELL=IPREL DO 1131 ILM=1,NA IF(IMMM(ILM).EQ.0) GO TO 1131 IAA=IAA+1 ITRAA(IAA)=IPRELL+ILM-1 1131 CONTINUE CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1), > NA,IPREL,MULRE,INC,dnorma) if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU) 113 CONTINUE if (idenorm.eq.1) then call erreur(1049) do k=0,muinc,inc do i=1,inc vectbb(i+k)=0.d0 enddo enddo endif * desactivation generale do ivs=1,nnoe lign=ilign(ivs) segdes lign enddo * * on verifie la precision du resultat si pas de noen ou pas de modes d'ensemble ** write(6,*) ' noen nens ',noen,nens if (jlicre.ne.0) then if (noen.ne.0.or.nens.eq.0) then if(mulre.eq.1.and..not.exists) then * U+dU if (icorrec.eq.0) then do i=1,inc mvect2.vectbb(i)=vectbb(i) enddo else do i=1,inc mvect2.vectbb(i)=mvect2.vectbb(i)+vectbb(i) enddo endif ilicre=jlicre segact ilicre ligcre=ligcrp segact ligcre * K*U * sauver l'ancien residu if (mvect4.ne.0) segsup mvect4 if (icorrec.ne.0) mvect4 = mvect3 segini mvect3 call graco7(ilicre,mvect2,mvect3,inc,nbthr,0) * F-K*U crit=xpetit/xzprec crit2 = 0.D0 do i=1,inc crit=max(crit,abs(mvect1.vectbb(i)),abs(mvect3.vectbb(i))) mvect3.vectbb(i)=mvect1.vectbb(i)-mvect3.vectbb(i) ** write(6,*) ' mvect1 mvect3 ',mvect1.vectbb(i),mvect3.vectbb(i) crit2=max(crit2,abs(mvect3.vectbb(i))) enddo * calcul lambda delta u xlbd = 1. if (icorrec.ne.0) then xd = 0. xn = 0. do i=1,inc r34 =mvect3.vectbb(i)-mvect4.vectbb(i) xn = xn + mvect4.vectbb(i) * r34 xd = xd + r34 * r34 enddo xd = max(xd,abs(xn)*xpetit,xpetit) xlbd = -xn / xd ** write(6,*) ' lambda ',xlbd if (xlbd.lt.-1..or.xlbd.gt.2.) xlbd = 0.5 endif * test if (icorrec.le.8) crite = crit * sqrt(xzprec*xszpre) * critere lache apres les premieres ite©rations. demande une reflexion approfondie. if (icorrec.ge.13) crite = crit * 1d-4 if(icorrec.gt.12.and.nensl.ne.0) crite=xgrand*xzprec ** write(6,*) ' crite crit2',crite,crit2,icorrec if (crit2.gt.crite) then * il faut iterer do i2=1,inc vectbb(i2)=mvect3.vectbb(i2)* xlbd enddo ** segsup mvect3 icorrec=icorrec+1 if (icorrec.lt.15.and.toutmem) goto 10000 endif * verif correcte segsup mvect1,mvect3 if (mvect4.ne.0) segsup mvect4 do i=1,inc vectbb(i)=mvect2.vectbb(i) enddo segsup mvect3 reaerr(1)=crit2/crit interr(1)=icorrec interr(2)=nensl *** isouci = 0 if (icorrec.ge.15.or.crit2/crit.gt.1d-7) then if (isouci.eq.0) then call erreur(1128) else call soucis(1128) call erreur(1129) endif endif if (icorrec.ge.1) then * l'erreur 1129 est informative donc pas de souci ** call erreur(1129) endif endif endif endif ** write(6,*) 'icorrec',icorrec C C **** ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE C C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ... NENS=IAA ** write(6,*) ' mondes nens nensl ',nens,nensl NBEN=0 MINCPO=IINCPO MIMIK=IIMIK SEGACT,MINCPO,MIMIK MELEME=IGEOMA SEGACT,MELEME NNOE=INCPO(/2) IINC1=INCPO(/1) * Pour creation chpoint des modes d'ensemble actifs IF (NOEN.EQ.0) then NNIN=IMIK(/2) NNNOE=NUM(/2) SEGINI MTRAV DO I=1,IMIK(/2) INCO(I)=IMIK(I) ENDDO ENDIF ** write(6,*) 'mondes creation mtrav:',mtrav,noen IF(NENS.EQ.0) GO TO 26 C ... Boucle sur les seconds membres ... DO 300 KI=0,MUINC,INC C ... XMA et XMAL seront le maximum des val. abs. des termes C du résultat partiel (avant la multiplication par MDNOR) C N° KI+1 multiplié par 1.e-10 ... XMA=xpetit/xzprec XMAL=xpetit/xzprec DO 30 KK=1,INC I=KK+KI if (ittr(kk).eq.0) then XMA=MAX(XMA,ABS(VECTBB(I))) else XMAL=MAX(XMAL,ABS(VECTBB(I))) endif 30 CONTINUE XMA=XMA*1.d-10 XMAL=XMAL*max(1d-10,crit2/crit*1d5) xmal=max(xma*1d-2,xmal) * write (6,*) ' mondes xma xmal ',xma,xmal C ... Boucle sur les modes d'ensemble ... ** write (6,*) ' nombre de modes d ensemble',nens iwrite = 0 DO 21 IA=1,NENS I1=ITRAA(IA) J=IPNO(I1) DO 22 K=1,IINC1 IF(INCPO(K,J).EQ.I1) GO TO 23 22 CONTINUE CALL ERREUR(5) call oooprl(0) RETURN 23 CONTINUE C ... Si ce n'est pas un multiplicateur, le déplacement doit être C < à XMA, sinon ERREUR ... * write (6,*) ' mondes i1 ittr vect xma ', * > i1,ittr(i1),vectbb(i1+ki),xma IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA ) then vectbb(i1+ki)=0.d0 GO TO 20 endif C ... Si c'est un multiplicateur, le multiplicateur doit être C < à XMAL sinon ERREUR ... i2=1 if (ittr(i1).ne.0) then i2=ittr(i1) vecs=vectbb(i1+ki)+vectbb(i2+ki) if (abs(vecs).le.xmal) then vectbb(i1+ki)=0.d0 vectbb(i2+ki)=0.d0 goto 20 endif endif * write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki), * > vectbb(i2+ki),k,imik(k) C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme C ... Si option NOEN on accumule les pts dans un meleme IF(NOEN.EQ.0) THEN IBIN(k,J)=1 IGEO(J)=num(1,j) vv = vectbb(i1+ki) bb(k,j)=vectbb(i1+ki) * mettre aussi le second multiplicateur de lagrange pour l'appel a dbbcf if (ittr(i1).ne.0) then jd = ipno(i2) IBIN(k,Jd)=1 IGEO(Jd)=num(1,jd) vv = vectbb(i1+ki) bb(k,jd)=vectbb(i2+ki) endif ENDIF NBEN=NBEN+1 MOTERR(1:4)=IMIK(K) INTERR(1)=NUM(1,J) IF(NOEN.EQ.0.OR.ISOUCI.EQ.1) THEN call soucis(149) GO TO 21 ENDIF CALL ERREUR(149) call oooprl(0) RETURN 20 CONTINUE C ... Messages d'information ... JJK=NUM(1,J) KIK=KI/INC +1 * on n'ecrit qu'une seule fois le message indetermination IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then if (imik(k).ne.'LX ') > WRITE(IOIMP,41) JJK,IMIK(K) iwrite=iwrite+1 endif IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then WRITE(IOIMP,42) JJK,IMIK(K) endif IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then WRITE(IOIMP,410)KIK,JJK,IMIK(K) iwrite=iwrite+1 endif IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1) & WRITE(IOIMP,420)KIK,JJK,IMIK(K) 21 CONTINUE 300 CONTINUE 41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ', * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ', * 'LA SUSDITE dans mondes') 42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ', * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ', * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes') 410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD ' *,I6,' INCONNUE ', * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ', * 'LA SUSDITE dans mondes') 420 FORMAT(' VECTEUR NUMERO ',I3,/, * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ', * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ', * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes') SEGDES,MELEME SEGDES,MINCPO SEGDES,MIMIK 26 CONTINUE if (noen.eq.0) then * test si nan ou inf dans le resultat inan = 0 DO 500 KI=0,MUINC,INC DO 501 KK=1,INC i = kk + ki if (abs(vectbb(i)).lt.xgrand) goto 501 inan = inan + 1 vectbb(i)=xgrand 501 continue 500 continue if (inan.ne.0) then nben = -inan call soucis(nben) endif endif * indiquer si besoin le nombre de modes d'ensemble excites * et le maillage des noeuds freinés * et le champoint des modes excites IF (NOEN.EQ.0) then * if (nben.ne.0) write (6,*) ' mondes nben ',nben CALL ECRENT(NBEN) * pour permettre a crechp de prendre le persistent lock call oooprl(0) ** write(6,*) 'mondes appeler crechp avec mtrav',mtrav CALL CRECHP(MTRAV,MCHPO4) call oooprl(1) if(mtrav.ne.0) SEGSUP MTRAV * dedualisation des multiplicateurs if (lagdua.ne.0) call dbbcf(mchpo4,lagdua) call ecrobj('CHPOINT ',mchpo4) endif MDNOR=IDNORM SEGACT,MDNOR DO 350 K=0,MUINC,INC DO 35 I = 1, INC VECTBB(I+K)=VECTBB(I+K)*DNOR(I) 35 CONTINUE * on force l'egalite des multiplicateurs de lagrange DO 36 I = 1, INC if (ITTR(I).ne.0) then vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))*0.5D0 vectbb(ittr(i)+k)=vectbb(i+k) endif 36 CONTINUE 350 CONTINUE SEGDES,MDNOR C ... On ne désactive MDNO1 que dans le cas où il est C différent de MDNOR ... IF(IDNORD.GT.0) THEN SEGDES,MDNO1 ENDIF SEGDES MMATRI SEGDES,MILIGN SEGDES,MVECTD SEGSUP ITTRV IF(NENSL.NE.0) SEGSUP,ITRAA IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD 1001 FORMAT(' SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5) CALL OOOMRU(0) call oooprl(0) RETURN END