1 : C SIMU21 SOURCE FANDEUR 10/12/14 21:19:41 6806 2 : 3 : SUBROUTINE SIMU21(IPLA,IPLB,iplc,NMOD,CONV,nbonve,numv1, 4 : $ iplval,iphi,IPERM,XPREC21) 5 : 6 : ********************************************************************* 7 : * 8 : * S I M U 2 1 9 : * ----------- 10 : * 11 : * FONCTION: 12 : * --------- 13 : * Passage au petit probleme (calcul de T et de D), 14 : * Resolution du petit probleme [T-(1/lambda)D]*Phi=0 15 : * Test de convergence 16 : * 17 : * PARAMETRES : (E)=ENTREE (S)=SORTIE 18 : * ----------- 19 : * IPLA LISTREEL (E) alpha 20 : * IPLB LISTREEL (E) beta 21 : * iplc LISTREEL (E) sign(d) 22 : * NMOD ENTIER (E) nombre de modes demandés pour ce cycle 23 : * CONV LOGIQUE (S) = vrai si on a trouvé nbonve > NMOD 24 : * nbonve ENTIER (S) nombre de modes ayant convergés 25 : * numv1 ENTIER (E) dimension du probleme reduit 26 : * iplval ENTIER (S) POINTEUR VERS LES VALEURS PROPRES CONVERGEES 27 : * iphi ENTIER (S) POINTEUR VERS LES VECTEURS PROPRES CONVERGES 28 : * IPERM ENTIER (S) POINTEUR VERS LE TABLEAU DES PERMUTATIONS SI 29 : * LES MODES (iplval et iphi) SONT RANGES DANS 30 : * UN ORDRE DIFFERENT DE 1 2 3 ... nbonve 31 : * XPREC21 FLOTTANT (E) PRECISION SUR LES VALEURS PROPRES 32 : * 33 : ********************************************************************* 34 : 35 : IMPLICIT INTEGER(I-N) 36 : IMPLICIT REAL*8 (A-H,O-Z) 37 : 38 : -INC CCOPTIO 39 : -INC SMLREEL 40 : -INC SMLENTI 41 : -INC SMLCHPO 42 : 43 : *=======================================================* 44 : * TYPE MATRIX et DECLARATIONS 45 : *=======================================================* 46 : 47 : SEGMENT MATRIX 48 : REAL*8 A(N,N) 49 : ENDSEGMENT 50 : 51 : POINTEUR IPRIG1.MATRIX, IPMAS1.MATRIX, IPHI.MATRIX 52 : POINTEUR IPLVAL.MLREEL ,IPLERR.MLREEL 53 : POINTEUR IPLA.MLREEL,IPLB.MLREEL,iplc.mlreel 54 : POINTEUR IPERM.MLENTI 55 : LOGICAL CONV,CONV2 56 : REAL*8 XPREC21 57 : 58 : mlree1=ipla 59 : mlree2=iplb 60 : ildim=numv1 61 : JG=ILDIM 62 : SEGINI,IPLERR,IPERM 63 : N=ILDIM 64 : SEGINI,iprig1,ipmas1 65 : do ibr=1,ildim 66 : iprig1.a(ibr,ibr)=iplc.prog(ibr) 67 : enddo 68 : do ibr=1,ildim-1 69 : ipmas1.a(ibr,ibr)=ipla.prog(ibr) 70 : ipmas1.a(ibr,ibr+1)=iplb.prog(ibr) 71 : ipmas1.a(ibr+1,ibr)=iplb.prog(ibr) 72 : enddo 73 : ipmas1.a(ildim,ildim)=ipla.prog(ildim) 74 : * 75 : c write(6,*) 'simu21: appel sespa3 pour resoudre [T-(1/l)D]Phi=0' 76 : c write(6,*) ' D=',(iplc.prog(iuy),iuy=1,ildim) 77 : c write(6,*) ' alpha=',(ipla.prog(iuy),iuy=1,ildim) 78 : c write(6,*) ' beta=',(iplb.prog(iuy),iuy=1,ildim) 79 : 80 : *=======================================================* 81 : * Resolution du systeme reduit tridiagonal 82 : *=======================================================* 82 : CALL SESPA3 ( IPMAS1, IPRIG1, IPHI, IPLVAL) 83 : IF(IERR .NE. 0) RETURN 84 : 85 : 86 : *=======================================================* 87 : * Test d'arret basé sur un estimateur du residu qui est egalement lié 88 : * a l'erreur relative des valeurs propres 89 : *=======================================================* 90 : * 91 : * REMARQUE IMPORTANTE : 92 : * bp 25/10/2010 : on ne parvient pas a demontrer (ni analytiquement et 93 : * ni numeriquement) ce lien avec l algo de Lanczo.eso tel quil est ecrit 94 : * (utilisation de la K-norme) alors qu il est evident avec la M-norme 95 : * cclusion : il faudrait passer a la M-norme un jour... 96 : * pour l instant, on se contente d etre + severe sur la precision demandee 97 : XPREC0 = XPREC21 98 : XPREC21 = XPREC21**1.5 99 : 100 : XBETA=IPLB.PROG(ILDIM) 101 : DO 100 IB100 = 1 , ILDIM 102 : XLAMB=IPLVAL.PROG(IB100) 103 : XPHIN= IPHI.A( ILDIM, IB100 ) 104 : * erreur absolue 105 : * XERR = ABS( ( XBETA * XPHIN ) ) 106 : * erreur relative (sur XLAMB=1/lambda) 107 : XERR = ABS( ( XBETA * XPHIN )/ XLAMB ) 108 : IPLERR.PROG( IB100 ) = XERR 109 : 100 CONTINUE 110 : 111 : * lambda = 1 / (1/ lambda) 112 : DO 200 IB200=1,ILDIM 113 : iplval.prog(ib200)=1.D0/IPLVAL.PROG(IB200) 114 : 200 CONTINUE 115 : c if (iimpi.ge.6) then 116 : c write(6,*) 'simu21: iplval=',(iplval.prog(iou),iou=1,ILDIM) 117 : c write(6,*) 'simu21: IPLERR=',(IPLERR.prog(iou),iou=1,ILDIM) 118 : c endif 119 : 120 : 121 : *=========================================================* 122 : * On regarde si les NMODs premieres valeurs ont convergé 123 : *=========================================================* 124 : 125 : jg=ildim 126 : segini,mlree1,mlree2 127 : 128 : c----- Stratégie 1 : On trie par module croissant --------- 128 : call trifla(iplval.prog(1),mlree1.prog(1),iplerr.prog(1), 129 : $ mlree2.prog(1),jg,IPERM.lect(1)) 130 : 131 : 132 : c----- Stratégie 2 : On trie par erreur croissante --------- 133 : c Possible car tri des vp fait par simul7 + tard 134 : c On obtient + de modes par cycle au prix d un + grand nombre de trous 135 : c => implique une modification de strate(gie) pour aller moins loin 136 : c et pour controler le spectre qui est tres diffus... 136 : c call trifla(iplerr.prog(1),mlree2.prog(1), 137 : c $ iplval.prog(1),mlree1.prog(1),jg,IPERM.lect(1)) 138 : cbp 10/2010: finalement abandonné car difficile a mettre en oeuvre... 139 : 140 : segsup,mlree1,mlree2 141 : IF (IERR .NE. 0) RETURN 142 : if (iimpi.ge.6) then 143 : write(6,*) 'simu21: iplval=',(iplval.prog(iou),iou=1,ILDIM) 144 : write(6,*) 'simu21: IPLERR=',(IPLERR.prog(iou),iou=1,ILDIM) 145 : endif 146 : 147 : c convergence sur l'erreur relative 148 : CONV = .TRUE. 149 : nbonve0=nbonve 150 : nbonve=NMOD 151 : 152 : 153 : c------boucle sur les valeur propres 154 : DO 300 IB300 = 1,ildim 155 : 156 : c -erreur relative < tolerance? 157 : EPS1 = IPLERR.PROG( IB300 ) 158 : EPS=ABS(EPS1) 159 : IF ( EPS .GT. XPREC21 ) THEN 160 : if (IB300.eq.1) then 161 : nbonve=ib300-1 162 : if(ib300.le.NMOD+1) CONV=.FALSE. 163 : GOTO 310 164 : endif 165 : c pour eviter de redemarrer a cause d un mode multiple, 166 : c puisque cet estimateur d erreur est assez mechant, 167 : c on va tolerer une erreur relative plus grande 168 : EVP1 = (iplval.prog(IB300)) / (iplval.prog(IB300-1)) 169 : EVP1 = ABS (EVP1 - 1.D0) 170 : if (EVP1.LT.XPREC21.AND.EPS.LT.(2.D2*XPREC21)) then 171 : if (iimpi.ge.6) then 172 : write(6,108) IB300,EPS 173 : 108 FORMAT(2X,'on sauve le candidat',I5, 174 : $ ' car mode multiple et erreur acceptable=',E10.5) 175 : endif 176 : elseif(EVP1.LT.(2.D2*XPREC21).AND.IB300.gt.2) then 177 : if (iimpi.ge.6) then 178 : write(6,109) IB300,EPS 179 : 109 FORMAT(2X,'on annule le candidat',I5,' car ', 180 : $ 'probable mode multiple et erreur inacceptable=',E10.5) 181 : endif 182 : nbonve=ib300-2 183 : if(ib300.le.NMOD+1) CONV=.FALSE. 184 : GOTO 310 185 : else 186 : * a priori ce n est pas un mode multiple et erreur inacceptable 187 : nbonve=ib300-1 188 : if(ib300.le.NMOD+1) CONV=.FALSE. 189 : GOTO 310 190 : endif 191 : ENDIF 192 : 300 CONTINUE 193 : c------fin de boucle sur les valeur propres 194 : 310 CONTINUE 195 : 196 : * on remet la valeur initiale 197 : XPREC21 = XPREC0 198 : *=========================================================* 199 : * destruction des objets de travail 200 : *=========================================================* 201 : mlreel=IPLERR 202 : segsup,mlreel 203 : mlree1=iplval 204 : mlree2=IPERM 205 : SEGDES,mlree1,mlree2 206 : SEGSUP,IPRIG1,IPMAS1 207 : 208 : RETURN 209 : END 210 : 211 :
© Cast3M 2003 - All rights reserved.
Disclaimer