mondes
C MONDES SOURCE PV090527 24/06/07 21:15:02 11935 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 -INC PPARAM -INC CCOPTIO -INC CCREEL -INC SILICRE -INC CCASSIS SEGMENT ITTRV(MULRE) SEGMENT,ITRAA(NENSL) SEGMENT/BID/(IBIDON(IIMAX)) logical toutmem toutmem=.true. DNORMA=0.D0 DNORMB=0.D0 crit=xpetit call oooprl(1) nbthr=nbthrs if (noen.eq.0) then nbelem=0 nbnn =1 nbsous=0 nbref =0 segini,ipt8 ipt8.itypel=1 endif 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 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 2001 continue 2010 continue ilmin=ivs+1 DO 610 IVS=IDEP,NNOE 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 > 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 C ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ... 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, C ... puis ce dernier est décrémenté ... ILMAX=ILMAX-1 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 > NA,IPREL,MULRE,INC,IVS,LLOM,dnorma) 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 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 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 ... 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 GO TO 67 C ... Si on a des pb avec activation, on désactive et on diminue ILMAX ... 68 CONTINUE toutmem=.false. ILMAX=ILMAX-1 C ... puis on recommence la tentative ... GO TO 70 C ... On a réussi à activer ... 67 CONTINUE NA=IMMM(/1) IPRELL=IPREL > NA,IPREL,MULRE,INC,dnorma) 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 NA=IMMM(/1) IPRELL=IPREL > NA,IPREL,MULRE,INC,dnorma) 13 CONTINUE C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ... DO 113 J=IDEP-1,1,-1 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 > NA,IPREL,MULRE,INC,dnorma) 113 CONTINUE if (idenorm.eq.1) then do k=0,muinc,inc do i=1,inc vectbb(i+k)=0.d0 enddo enddo endif * desactivation generale do ivs=ilmin,nnoe 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 (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 segini mvect3 * F-K*U do i=1,inc mvect3.vectbb(i)=mvect1.vectbb(i)-mvect3.vectbb(i) ** write(6,*) ' mvect1 mvect3 ',mvect1.vectbb(i),mvect3.vectbb(i) enddo * test * critere lache apres les premieres ite©rations. demande une reflexion approfondie. if(icorrec.gt.7.and.nensl.ne.0) crite=xgrand*xzprec ** write(6,*) ' crite crit2',crite,crit2,inc * il faut iterer enddo segsup mvect3 icorrec=icorrec+1 if (icorrec.lt.10.and.toutmem) goto 10000 endif * verif correcte segsup mvect1,mvect3 do i=1,inc vectbb(i)=mvect2.vectbb(i) enddo segsup mvect3 interr(1)=icorrec interr(2)=nensl if (isouci.eq.0) then else call soucis(1128) endif endif if (icorrec.ge.1) then * l'erreur 1129 est informative donc pas de souci ** call erreur(1129) 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 IF(NENS.EQ.0) GO TO 26 MINCPO=IINCPO MIMIK=IIMIK SEGACT,MINCPO,MIMIK MELEME=IGEOMA SEGACT,MELEME NNOE=INCPO(/2) IINC1=INCPO(/1) 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=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 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 ... if (ittr(i1).ne.0) then if (abs(vecs).le.xmal) then vectbb(i1+ki)=0.d0 goto 20 endif endif ** write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki), ** > vectbb(i2+ki),xmal 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 nbelem=ipt8.num(/2)+1 segadj ipt8 ipt8.num(1,nbelem)=num(1,j) 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 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 IF (NOEN.EQ.0) then * write (6,*) ' mondes nben ',nben segdes ipt8 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 * write (6,*) ' dans mondes ',i,ittr(i) if (vectbb(i+k).eq.0.d0.or.vectbb(ittr(i)+k).eq.0.d0) then * write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k) vectbb(i+k)=0.d0 vectbb(ittr(i)+k)=0.d0 goto 36 endif vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))/2 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
© Cast3M 2003 - Tous droits réservés.
Mentions légales