1 : C LANCZO SOURCE PV 11/02/18 21:15:06 6875 2 : SUBROUTINE LANCZO (IPRIGI,IPMASS,Nbmode,NMOD,nvecma,Nmodopt, 3 : $ xktale,iflu,nbonve,iplval,iphi,IPERM,icycle,Northo) 4 : 5 : ************************************************************************ 6 : * 7 : * L A N C Z O 8 : * ----------- 9 : * 10 : * FONCTION: 11 : * --------- 12 : * 13 : * APPELE PAR LE SOUS-PROGRAMME "SIMUL1". 14 : * DETERMINE D'UNE SUITE DE 'CHPOINTS' X K-ORTHOGONAUX 15 : * PAR L'ALGORITHME DE LANCZOS 16 : * ET RESOLUTION DU PROBLEME REDUIT TRIDIAGONAL ASSOCIE 17 : * 18 : * MODE D'APPEL: 19 : * ------------- 20 : * 20 : * CALL LANCZO (IPRIGI,IPMASS,Nbmode,NMOD,nvecma,Nmodopt, 21 : * $ xktale,iflu,nbonve,iplval,iphi,IPERM,icycle,ortho) 22 : * 23 : * PARAMETRES : (E)=ENTREE (S)=SORTIE 24 : * ----------- 25 : * 26 : * IPRIGI ENTIER (E) POINTEUR SUR LA MATRICE DE RIGIDITE K-W^2m 27 : * (OBJET 'RIGIDITE') EVENTUELLEMENT DECALEE. 28 : * IPMASS ENTIER (E) POINTEUR SUR LA MATRICE MASSE M 29 : * (OBJET 'RIGIDITE'). 30 : * Nbmode ENTIER (E) NOMBRE DE MODES PROPRES DEMANDES PAR 31 : * L'UTILISATEUR. 32 : * NMOD ENTIER (E) NOMBRE DE MODES VISES POUR CE CYCLE 33 : * nvecma ENTIER (E) NOMBRE DE VECTEURS X MAXIMUM 34 : * Nmodopt ENTIER (E/S) NOMBRE DE MODES OPTIMUM % tcpu 35 : * xktale REEL (E) temps pour decomposer la matrice K'=LDLT 36 : * (S) temps passé dans Lanczos 37 : * iflu ENTIER (E) =1 si terme en PI (meca flu) 38 : * nbonve ENTIER (S) NOMBRE DE MODES ESTIMES CONVERGES 39 : * iplval ENTIER (S) POINTEUR VERS LES VALEURS PROPRES CONVERGEES 40 : * iphi ENTIER (S) POINTEUR VERS LES VECTEURS PROPRES CONVERGES 41 : * IPERM ENTIER (S) POINTEUR VERS LE TABLEAU DES PERMUTATIONS SI 42 : * LES MODES (iplval et iphi) SONT RANGES DANS 43 : * UN ORDRE DIFFERENT DE 1 2 3 ... nbonve 44 : * icycle ENTIER (E) NUMERO DU CYCLE EN COURS 45 : * Northo ENTIER (E) PAR RAPPORT A COMBIEN de MODES FAUT-IL 46 : * M-ORTHOGONALISER LES VECTEURS X^i 47 : * 48 : * 49 : * MODE DE FONCTIONNEMENT: 50 : * ----------------------- 51 : * 52 : * explications donnes dans le corps du programme 53 : * pour + de detail, voir [BATHE] page 952 (ou autre ref.) 54 : * 55 : * DESCRIPTION DES VARIABLES : cf. commentaires au cours du programme 56 : * -------------------------- 57 : * 58 : * CREATION, MODIFICATIONS : 59 : * ------------------------ 60 : * 61 : * CREATION : PASCAL MANIGOT 22 AVRIL 1985 62 : * REFONTE : THIERRY CHARRAS, BENOIT PRABEL 2010 63 : * 64 : * LANGAGE : ESOPE (FORTRAN77) 65 : * 66 : ************************************************************************ 67 : * 68 : IMPLICIT INTEGER(I-N) 69 : IMPLICIT REAL*8 (A-H,O-Z) 70 : 71 : -INC CCREEL 72 : -INC CCOPTIO 73 : -INC CCHAMP 74 : -INC SMLREEL 75 : -INC SMVECTD 76 : -INC SMLENTI 77 : -INC SMRIGID 78 : -INC SMCHPOI 79 : pointeur mvect4.mvectd,mvect5.mvectd,mvect6.mvectd,mvect7.mvectd 80 : pointeur mlent4.mlenti,mlent5.mlenti,mlent6.mlenti,mlent7.mlenti 81 : pointeur mlen10.mlenti,mlre10.mlreel,mlen11.mlenti,mlre11.mlreel 82 : segment idemen(0) 83 : * 84 : * REGOUPEMENT DE RENSEIGNEMENTS POUR "LANCZO" et "SIMUL.": 85 : COMMON/CLANCZ/ IPV1,IPMV1,IPLMOX,IPLMOY 86 : COMMON/CSIMUL/ W2SHIF,XPREC21,mvecri,IPVECX,IPVECZ,IPMZ,IPFREZ 87 : * 88 : REAL*8 W2SHIF,XPREC21,XORTHO,xKTKx2 89 : PARAMETER (DEUXPI = (2.D0*XPI)) 90 : * 91 : LOGICAL CONV,ortho,quit1,PRESQ 92 : POINTEUR IPLA.MLREEL,IPLB.MLREEL,iplc.mlreel 93 : POINTEUR iplas.mlreel,iplbs.mlreel 94 : 95 : * 96 : * ----quelques initialisations ------- 97 : * 98 : xktlan = 0.D0 99 : yktort = 0.D0 100 : * -on teste l efficacite (temps CPU) toutes les itest iterations 101 : itest = 20 102 : itest0 = itest 103 : Nvecopt = nvecma 104 : * 105 : * -NOMBRV = nbre de vecteurs estimes necessaires pour le calcul des 106 : * NMOD modes demandes 107 : NOMBRV = (2*NMOD) + 8 108 : NOMBRV = MIN(NOMBRV,nvecma) 109 : * nvectm = borne sup de ce nombre 110 : nvectm = 2*nombrv 111 : nvectm = MIN(nvectm,nvecma) 112 : 113 : * -on M-orthogonalise par % aux Northo modes les + proches 114 : c Northo = min(NMOD,(nvecma/3)) 115 : ortho = (icycle.ne.1).and.(Northo.gt.0) 116 : 117 : * -on va tester la convergence apres le calcul de NMBRV1 vecteurs 118 : * (cette valeur est augmentée de 10% a chaque fois qu'on n'a pas convergée) 119 : NMBRV1 = NOMBRV*0.8 120 : NMBRV1 = min(NMBRV1,nvecma) 121 : if (iimpi.ge.6) then 122 : write(6,*) 'lanczo:' 123 : write(6,*) ' nbre vect estimés necessaires',NOMBRV,' - ',nvectm 124 : write(6,*) ' pour le calcul de ',NMOD,' modes' 125 : write(6,*) ' on teste la convergence tous les ',NMBRV1,' itérés' 126 : write(6,*) ' on orthogonalise par % a ',Northo,' modes' 127 : endif 128 : * -nombre de vecteurs ayant convergés /max atteint 129 : nbonve=0 130 : nbonmax=0 131 : * -detecte que la re-orthogonalisation a quasi-annulé le vecteur courant 132 : c RLIM=1.D-3 133 : RLIM=sqrt(XPREC21) 134 : quit1 = .false. 135 : iquit1=0 136 : c on considere qu il faut reorthognaliser si le produit scalaire 137 : * est > XORTHO (avec **1.6 XORTHO=3E-13 si xprec21=1.5E-8) 138 : c XORTHO = 1.D-12 139 : c XORTHO = XPREC21**1.5 ne suffit pas 140 : XORTHO = XPREC21**1.6 141 : 142 : * -initialisation de listentier vers chpoint/vecteur X et KX 143 : jg = nombrv 144 : segini mlent3,mlent4 145 : * -initialisation de listreel alpha et beta et sign(xtKx) 146 : segini ipla,iplb,iplc,iplas,iplbs 147 : s1=1.d0 148 : s0=1.d0 149 : 150 : * -passage chpoint->vecteur: y^0 = IPV1 -> ivev1 , r^1 = M*y^0 = IPMV1 -> ivemv1 151 : mrigid=iprigi 152 : segact,mrigid 153 : icholi=ichole 154 : mvecri=ivecri 155 : segdes mrigid 155 : call chv3(icholi,IPV1,ivev1,1) 156 : mvect1=ivev1 157 : segact mvect1*mod 158 : inc = mvect1.vectbb(/1) 159 : 160 : 161 : * -eventuellement, initialisation du vecteur de départ y^0 162 : * dans un s.e. M-orthogonal aux modes deja convergés 163 : if (ortho) then 164 : * -recup 165 : mlent6 = IPVECZ 166 : mlent7 = IPMZ 167 : mlree1 = IPFREZ 168 : c write(6,*) 'IPFREZ,IPMZ,IPVECZ=',IPFREZ,IPMZ,IPVECZ 169 : if (mlree1.gt.0) then 170 : * -on cherche les modes deja convergés les + proches du shift actuel 171 : FRSHIF = (sqrt(W2SHIF)) / DEUXPI 172 : segact,mlree1 173 : JG1=mlree1.prog(/1) 174 : Northo=min(Northo,JG1) 175 : JGsave=JG 176 : JG=JG1 177 : segini,mlen10,mlen11 178 : segini,mlre10,mlre11 179 : JG=JGsave 180 : do j1=1,JG1 181 : mlen10.lect(j1)=j1 182 : FRj1 = mlree1.prog(j1) 183 : mlre10.prog(j1)=abs(FRSHIF-FRj1) 184 : enddo 184 : call triflo(mlre10.prog(1),mlre11.prog(1), 185 : $ mlen10.lect(1),mlen11.lect(1),jg1) 186 : * mlen10(j1)= j10 siginifie que j10 est le j1ieme + proche mode 187 : segsup,mlre10,mlre11,mlen11 188 : if (iimpi.ge.6) then 189 : write(6,*)'Lanczos: M-orthog de y^0 par rapport a Z ' 190 : write(6,*)'i=',(mlen10.lect(iou),iou=1,Northo) 191 : write(6,*)'f=',(mlree1.prog(mlen10.lect(iou)),iou=1,Northo) 192 : endif 193 : segdes,mlree1 194 : * -M-orthogonalisation % aux modes deja convergés 194 : call dtchpo(ipv1) 194 : call dtchpo(ipmv1) 195 : segact,mlent6,mlent7 196 : c JG0 = mlent6.LECT(/1) 197 : c do jpo=1,JG0 198 : do j1=1,Northo 199 : jpo=mlen10.lect(j1) 200 : * on recupere x^j et r^j = M*z^j / (z^j*M*z^j) 201 : mvect6= mlent6.lect(jpo) 202 : mvect7= mlent7.lect(jpo) 203 : segact mvect6,mvect7 204 : * on calcule xtru = y^i * r^j 205 : xtru = 0.D0 206 : do iou=1,inc 207 : xtru= xtru + (mvect7.vectbb(iou)*mvect1.vectbb(iou)) 208 : enddo 209 : c if(abs(xtru).lt.1.D-12) go to 88 210 : if(abs(xtru).lt.XORTHO) go to 88 211 : * y^i = y^i - xtru * z^{j} et idem pour Ky^i 212 : do iou=1,inc 213 : mvect1.vectbb(iou) = mvect1.vectbb(iou) 214 : $ - (xtru*mvect6.vectbb(iou)) 215 : enddo 216 : 88 continue 217 : segdes mvect6,mvect7 218 : c mis en commentaire car on reothogonalise chaque vecteur x 219 : c finalement, on ne reorthogonalise que le vecteur y^0 220 : enddo 221 : segdes,mlent6,mlent7 222 : segsup,mlen10 222 : call crech2(IPV1,mvect1,mvecri,1) 222 : CALL MUCPRI(IPV1,IPMASS,IPMV1) 223 : endif 224 : endif 225 : * creation du second membre initial 225 : call chv2(icholi,IPMV1,ivemv1,1) 226 : mvect2=ivemv1 227 : segact mvect2*mod 228 : 229 : c menage 229 : call dtchpo(IPV1) 230 : segsup,mvect1,mvect2 231 : 232 : * -preparation a la resolution 233 : segini,idemen 234 : idemen(**)=0 235 : NOID=0 236 : NOEN=1 237 : * 238 : if(iimpi.ge.6) 239 : $ write(6,*) '-------------------------------------------------' 240 : * 241 : *-----Boucle iterative ----------------------------------------------* 242 : * 243 : DO 100 IB100=1,(nvecma+1) 244 : 245 : * -initialisation du temps cpu 246 : if (IB100.eq.itest) then 246 : call gibtem(xkt) 247 : xktlan = xktlan + xkt 248 : endif 249 : 250 : * -numero des vecteurs en jeu 251 : NUMV = IB100 252 : NUMV1 = NUMV - 1 253 : if(iimpi.ge.6) 254 : $ write(6,*) ' lanczo: -------i=',NUMV,'/' ,(nvecma+1),'--',rnm 255 : 256 : * -MISE A ZERO DES TERMES EN PI POUR LES ELTS LIQUIDES 257 : ipmv11 = ipmv1 257 : CALL ANCHPO(ipmv11,NOMDU(15),ipmv1) 257 : CALL DTCHPO(ipmv11) 258 : 259 : * -resolution de K*y^i = M*x^i (ou = K*y^1=My^0 si i=1) 260 : IDEMEN(1)=ipmv1 260 : CALL RESOU1 (IPRIGI,IDEMEN,NOID,NOEN,1.D-18) 261 : if (iflu.eq.1) then 262 : mrigid=iprigi 263 : segact mrigid 264 : icholi=ichole 265 : mvecri=ivecri 266 : segdes mrigid 267 : endif 268 : SEGACT IDEMEN*mod 269 : IPV2 = IDEMEN(1) 269 : call dtchpo(ipmv1) 270 : 271 : * y^i: passage chpo -> vecteur 271 : call chv3(icholi,ipv2,mvect3,1) 271 : call dtchpo(ipv2) 272 : segact mvect3*mod 273 : 274 : * -Application des coef beta et alpha 275 : * y^i = y^i - beta_{i-1} * x^{i-1} 276 : if (numv.gt.2) then 277 : mvect2=mlent3.lect(numv-2) 278 : b1s0=b1*s0 279 : do iou=1,inc 280 : mvect3.vectbb(iou) = mvect3.vectbb(iou) 281 : $ - (b1s0 * mvect2.vectbb(iou)) 282 : enddo 283 : endif 284 : 285 : * y^i = y^i - alpha_{i} * x^{i} 286 : if (numv.gt.1) then 287 : mvect1 = mlent3.lect(numv1) 288 : mvect4 = mlent4.lect(numv1) 289 : a1s1=a1*s1 290 : do iou=1,inc 291 : mvect3.vectbb(iou) = mvect3.vectbb(iou) 292 : $ - (a1s1 * mvect1.vectbb(iou)) 293 : enddo 294 : endif 295 : s0=s1 296 : 297 : * -calcul de K * y^i 297 : call crech2(ipv1,mvect3,mvecri,1) 297 : call mucpri(ipv1,IPRIGI,ipkv1) 297 : call chv2 (icholi,ipkv1,mvect4,1) 298 : segact,mvect3*mod 299 : segact,mvect4 300 : * comme on n'utilise que mvect3 et mvect4 , on detruit: 300 : call dtchpo(ipv1) 300 : call dtchpo(ipkv1) 301 : 302 : * -calcul de coe = y^i * K * y^i 303 : coe=0.D0 304 : do iou=1,inc 305 : coe= coe + mvect4.vectbb(iou)*mvect3.vectbb(iou) 306 : enddo 307 : segsup,mvect4 308 : 309 : unval = sqrt(abs(1.d0/coe)) 310 : s1 = abs(1.d0/coe)*coe 311 : xtrx = sign(1.d0,coe) 312 : xtrx1 = xtrx 313 : 314 : * -on commence par K-normaliser y^i (est-ce vraiment utile ?) 315 : c if ((coe.lt.XORTHO).and.(numv.ne.1)) then 316 : c NV = numv 317 : c WRITE(IOIMP,2000) NV,coe 318 : c 2000 FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE', 319 : c C /3X,'DIFFICULTE POUR NORMALISER LE ',I5, 320 : c C ' IEME VECTEUR ',/3X,'NORME PAR', 321 : c C ' RAPPORT A LA NORME INITIALE ',E12.5) 322 : c * on force l algo a faire les tests de convergence 323 : c * si on a trouvé des modes -> on quitte lanczos 324 : c * sinon on redémarre avec un nouveau vecteur (permis 1 fois) 325 : c quit1 = .true. 326 : c iquit1 = iquit1 + 1 327 : c endif 328 : do iou=1,inc 329 : mvect3.vectbb(iou) = mvect3.vectbb(iou)*unval 330 : enddo 331 : 332 : c xktre1 = TCPU pour resoudre + normaliser (lineaire avec numv) 333 : if (IB100.eq.itest) then 333 : call gibtem(xktre1) 334 : xktlan = xktlan + xktre1 335 : endif 336 : 337 : i08 = 0 338 : nort = 0 339 : 08 continue 340 : i08 = i08 + 1 341 : xtrum = 0.D0 342 : xtrum2 = 0.D0 343 : 344 : * -K-orthogonalisation % aux précedents (si ce n est pas le 1er passage) 345 : if (numv.gt.1) then 346 : xtrx1 = xtrx 347 : do ipo=1,numv1 348 : * on recupere x^j , r^j = K*x^j et s^j = sign(x^j*K*x^j) 349 : mvect1 = mlent3.lect(ipo) 350 : mvect2 = mlent4.lect(ipo) 351 : utru = iplc.prog(ipo) 352 : * on calcule xtru = y^i * r^j = y^i*K*x^j 353 : xtru = 0.D0 354 : do iou=1,inc 355 : xtru = xtru + (mvect2.vectbb(iou)*mvect3.vectbb(iou)) 356 : enddo 357 : xxtru = -1.D0 *xtru/utru 358 : xtrum = xtrum + abs(xxtru) 359 : xtrum2 = max(xtrum2,(abs(xxtru))) 360 : * re-K-orthogonalisation selective (=si correction xxtru > 5.D-12) 361 : c if(abs(xxtru).lt.1.D-12) go to 18 362 : if(abs(xxtru).lt.XORTHO) go to 18 363 : nort = nort + 1 364 : * y^i = y^i - (xtru/utru)*x^{j} 365 : do iou=1,inc 366 : mvect3.vectbb(iou) = mvect3.vectbb(iou) 367 : $ + (xxtru*mvect1.vectbb(iou)) 368 : enddo 369 : xtrx1=xtrx1-((xtru*xtru)/utru) 370 : 18 continue 371 : enddo 372 : endif 373 : 374 : if(iimpi.ge.7) write(6,*) xtrum,xtrum2 375 : c c if(xtrum.ge.XORTHO) goto 08 => explose 376 : c if(xtrum.ge.XPREC21) goto 08 377 : c c if(xtrum.ge.rlim) goto 08 378 : c c if(xtrum.ge.0.707) goto 08 379 : c c bp: l'algo redémarrerait-il avec un nouveau vecteur? 380 : c bp29/10/2010: avec XPREC21 comme precision, on constate qu'1 ou 2 boucle 381 : c d'orthogonalisation suffit, mais pour eviter de tourner en rond, 382 : c on limite a 2 passages 383 : if (xtrum.ge.XPREC21) then 384 : if(i08.lt.2) goto 08 385 : WRITE(IOIMP,1088) numv,xtrum 386 : 1088 FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE', 387 : C /3X,'DIFFICULTE POUR ORTHOGONALISER LE',I5,' IEME VECTEUR', 388 : C /3X,'SOMME DES PRODUITS SCALAIRE ',E12.5) 389 : * on force l algo a faire les tests de convergence 390 : * si on a trouvé des modes -> on quitte lanczos 391 : * sinon on redémarre avec un nouveau vecteur (permis 1 fois) 392 : quit1 = .true. 393 : iquit1 = iquit1 + 1 394 : endif 395 : 396 : 397 : c * -M-orthogonalisation % aux modes deja convergés 398 : c if (ortho) then 399 : c c JG0 = mlent6.LECT(/1) 400 : c c do jpo=1,JG0 401 : c do j1=1,Northo 402 : c jpo=mlen10.lect(j1) 403 : c * on recupere x^j et r^j = M*z^j / (z^j*M*z^j) 404 : c mvect6= mlent6.lect(jpo) 405 : c mvect7= mlent7.lect(jpo) 406 : c c deja fait segact mvect6,mvect7 407 : c * on calcule xtru = y^i * r^j 408 : c xtru = 0.D0 409 : c do iou=1,inc 410 : c xtru= xtru + (mvect7.vectbb(iou)*mvect3.vectbb(iou)) 411 : c enddo 412 : c c if(abs(xtru).lt.1.D-12) go to 28 413 : c if(abs(xtru).lt.XORTHO) go to 28 414 : c * y^i = y^i - xtru * z^{j} et idem pour Ky^i 415 : c do iou=1,inc 416 : c mvect3.vectbb(iou) = mvect3.vectbb(iou) 417 : c $ - (xtru*mvect6.vectbb(iou)) 418 : c enddo 419 : c 28 continue 420 : c c segdes mvect6,mvect7 421 : c c mis en commentaire car on reothogonalise chaque vecteur x 422 : c enddo 423 : c c segdes,mlent6,mlent7 424 : c endif 425 : c 426 : c bp 27/10/2010 : finalement, on ne reorthogonalise que le vecteur y^0 427 : c car semble etre trop contraignant pour lanczo 428 : 429 : 430 : c xktort = TCPU pour orthogonaliser (proport à numv*numv1) 431 : if (IB100.eq.itest) then 431 : call gibtem(xktort) 432 : xktlan = xktlan + xktort 433 : endif 434 : 435 : * -re-calcul de r^{i+1} = K*y^i 435 : call crech2(ipv1,mvect3,mvecri,1) 435 : call mucpri(ipv1,IPRIGI,ipkv1) 435 : call chv2 (icholi,ipkv1,mvect4,1) 436 : segact,mvect3*mod 437 : segact mvect4*mod 438 : 439 : * calcul de xtrx2 = y^i * K * y^i apres re-orthog 440 : xtrx2 = 0.D0 441 : do iou=1,inc 442 : xtrx2 = xtrx2 + (mvect3.vectbb(iou)*mvect4.vectbb(iou)) 443 : enddo 444 : 445 : * -on verifie que l orthogonalisation ne change pas trop le vecteur 446 : c (ce qui devrait etre vrai en arithmétique exacte pour les vecteurs non redémarrés) 447 : rnm = abs(xtrx2/xtrx) 448 : if ((rnm.lt.rlim).and.(numv.ne.1)) then 449 : NV = numv 450 : WRITE(IOIMP,1000) NV,RNM 451 : 1000 FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE', 452 : C /3X,'DIFFICULTE POUR ORTHOGONALISER LE ',I5, 453 : C ' IEME VECTEUR ',/3X,'NORME APRES ORTHOGONALISATION PAR', 454 : C ' RAPPORT A LA NORME INITIALE ',E12.5) 455 : * on force l algo a faire les tests de convergence 456 : * si on a trouvé des modes -> on quitte lanczos 457 : * sinon on redémarre avec un nouveau vecteur (permis 1 fois) 458 : quit1 = .true. 459 : iquit1 = iquit1 + 1 460 : endif 461 : 462 : * -K-normalisation : 463 : coef=1.D0/sqrt(abs(xtrx2)) 463 : call mulpar(mvect3,mvect4,IPV1,IPKV1,coef) 463 : call DTCHPO(IPKV1) 464 : 465 : s1=sign(1.D0,xtrx2) 466 : 467 : * -on recrée les chpoints x^{i+1} et K*x^{i+1} à partir des vecteurs modifiés 468 : c 469 : c -evaluation du TCPU + optimisation du nbre de vecteurs utilisés 470 : if (IB100.eq.itest) then 471 : c xktres = TCPU pour resoudre + normaliser (lineaire avec numv) 471 : call gibtem(xktre2) 472 : xktlan = xktlan + xktre2 473 : xktres = xktre0 + xktre1 + xktre2 474 : itest = int(1.2*itest) + 1 475 : 476 : if (IB100.eq.itest0) then 477 : yktort = xktort / numv1 478 : else 479 : yktort = 0.5*(yktort + (xktort / numv1)) 480 : endif 481 : 482 : c si xktort pas tres representatif, on le remplace 483 : if(xktort.le.1.) xktort = numv1*yktort 484 : c 485 : c si orthog. ne coute pas cher, on la retestera + tard 486 : if (xktort.lt.0.6.or.IB100.eq.itest0) then 487 : * 488 : if (iimpi.ge.7) then 489 : WRITE(6,1008) IB100,xktale,xktres,xktort,nort,xopt,Nmodopt 490 : 1008 FORMAT('IB100,xktale,xktres,xktort,nort,xopt,Nmodopt =', 491 : C I5,' ',E10.5,' ',E10.5,' ',E10.5,' ',I5,' ',E10.5,' ',I5) 492 : endif 493 : 494 : c sinon calcul du nombre de vecteurs optimal %TCPU (=nopt) 495 : else 496 : c bp : on va moyenner yktort pour le rendre moins sensible et variable 497 : xkt0 = (1.*xktale) + (8.*xktres) + (32.*yktort) 498 : xkt2 = (4.*yktort) + (2.*xktort/real(nort+numv1)) 499 : * xopt = sqrt(xkt0/xkt2) -> prise en compte du rebouchage de trou dans 50% des cas (*1.5) 500 : xopt = 1.5*sqrt(xkt0/xkt2) 501 : Nmodopt = INT(xopt) + 1 502 : Nvecopt = (2*Nmodopt) + 8 503 : 504 : if (iimpi.ge.7) then 505 : WRITE(6,1009) IB100,xktale,xktres,xktort,nort,xopt,Nmodopt 506 : 1009 FORMAT('IB100,xktale,xktres,xktort,nort,xopt,Nmodopt =', 507 : C I5,' ',E10.5,' ',E10.5,' ',E10.5,' ',I5,' ',E10.5,' ',I5) 508 : endif 509 : 510 : c si on pensait ne faire qu'1 cycle et que l on veut moins 511 : c de modes que l optimum, alors on n a pas interet a cycler 512 : if (Nbmode.eq.Nmod.and.Nmodopt.gt.int(Nbmode*0.8)) then 513 : Nmodopt = Nmod 514 : Nvecopt = NOMBRV 515 : c sinon on propose de cycler avec les valeurs optimisées: 516 : else 517 : Nmodopt = min(Nmodopt,Nvecma) 518 : cbp: l'ajustement du nombre de modes calculés ne concerne que le 1er cycle* 519 : c sinon pb lorsque l on bouche les trous.... 520 : if (icycle.eq.1) then 521 : c il n y a pas d interet a calculer + de mode que demandé 522 : Nmod = min(Nmod,Nmodopt) 523 : Nvecopt= min(Nvecopt,Nvecma) 524 : NOMBRV = Nvecopt 525 : NMBRV1 = int(0.8*NOMBRV) 526 : if (iimpi.ge.6) then 527 : write(6,*) ' nbre vect estimés necessaires',NOMBRV 528 : write(6,*) ' pour le calcul de ',NMOD,' modes' 529 : write(6,*) ' on teste la convergence tous les ',NMBRV1,' itérés' 530 : endif 531 : endif 532 : endif 533 : 534 : endif 535 : endif 536 : 537 : c initialisation du TCPU 538 : if (IB100.eq.(itest-1)) then 538 : call gibtem(xkt) 539 : xktlan = xktlan + xkt 540 : endif 541 : 542 : * on evite les erreurs d'arrondi pour s^{i+1} 543 : s1=sign(1.d0,s1) 544 : 545 : * on stocke x^{i+1} et r^{i+1} = K*x^{i+1} 546 : if (numv.ge.mlent3.lect(/1)) then 547 : jg= (numv+1) *1.2 548 : segadj mlent3,mlent4,ipla,iplb,iplc,iplas,iplbs 549 : endif 550 : mlent3.lect(numv)=mvect3 551 : mlent4.lect(numv)=mvect4 552 : iplc.prog(numv)=s1 553 : 554 : * calcul de a1 (alpha_i) et b1 (beta_i) 554 : call mucpri(ipv1,ipmass,ipmv1) 554 : call DTCHPO(IPV1) 554 : call chv2(icholi,ipmv1,mvect5,1) 555 : segact,mvect5 556 : a1=0.d0 557 : b1=0.d0 558 : * alpha_i = y^i*K*x^i = x^i*M*K^-1*K*x^i = x^i*M*x^i 559 : do iou=1,inc 560 : a1=a1 + mvect3.vectbb(iou)*mvect5.vectbb(iou) 561 : enddo 562 : ipla.prog(numv)=a1 563 : * alpha_i * s_i = (x^i * M * x^i) * sign(y^i * K * y^i) 564 : * on ajoute le signe s_i car x^i = yb/ |yb*K*yb|^0.5 565 : * donc x^i*K*x^i = yb*K*yb / |yb*K*yb|^0.5 = sign(yb*K*yb) 566 : dp = a1*s1 567 : a1s1=dp 568 : iplas.prog(numv)=dp 569 : * beta_i = y^i*K*x^{i-1} = x^{i}*M*x^{i-1} 570 : if (numv.gt.1) then 571 : mvect1 = mlent3.lect(numv1) 572 : do iou=1,inc 573 : b1=b1 + mvect1.vectbb(iou)*mvect5.vectbb(iou) 574 : enddo 575 : dp=b1*b1 *s0*s1 576 : iplbs.prog(numv1)=dp 577 : iplb.prog(numv1)=b1 578 : endif 579 : segsup,mvect5 580 : 581 : 582 : c xktres = TCPU pour resoudre + normaliser (lineaire avec numv) 583 : if (IB100.eq.(itest-1)) then 583 : call gibtem(xktre0) 584 : xktlan = xktlan + xktre0 585 : endif 586 : 587 : * -Resolution du probleme projete : [T]*Phi = 1/lambda*[D]*Phi 588 : * avec D =X*K*X = diag(s_i) = [s_1 0 ... ; 0 s_2 0 ... ; ...] 589 : * et T = X*M*X = [a_1 b_1 0 0 ... ; b_1 a_2 b_2 0 ... ; ...] 590 : IF ((quit1).or.(NUMV .GE. nmbrv1)) THEN 591 : 592 : c * ajout bp 22/10/2010 593 : c * calcul de xKTKx2 = ||r^{i+1}||_2^2 = ||K*x^{i+1}||_2^2 594 : c xKTKx2 = 0.D0 595 : c do iou=1,inc 596 : c xKTKx2 = xKTKx2 + mvect4.vectbb(iou)*mvect4.vectbb(iou) 597 : c enddo 598 : c write(6,*) 'lanczo : ||K*x^{i+1}||_2^2=',xKTKx2 599 : c * fin ajout bp 22/10/2010 600 : 600 : call simu21(ipla,iplb,iplc,NMOD,conv,nbonve,numv1, 601 : $ iplval,iphi,IPERM,XPREC21) 602 : 603 : if (iimpi.ge.6) then 604 : write(6,*) 'lanczo: CONV=',CONV,' car' 605 : write(6,*) 'pour ',numv,'vecteurs /',nvectm,'maxi' 606 : write(6,*) 'on a ',nbonve,'valeurs ok /',NMOD 607 : endif 608 : 609 : c c si le nombre de bon modes diminue trop, 610 : c c c'est qu'on aurait (probablement) du s'arreter plus tot ... 611 : c if (nbonve.ge.nbonmax) then 612 : c nbonmax= nbonve 613 : c numbon = numv1 614 : c endif 615 : c if (nbonmax.gt.0) then 616 : c nbonmin=0.7*nbonmax 617 : c nbonmin=max(2,nbonmin) 618 : c if (nbonve.le.nbonmin) then 618 : c call simu21(ipla,iplb,iplc,NMOD,conv,nbonve,numbon, 619 : c $ iplval,iphi,IPERM,XPREC21) 620 : c goto 110 621 : c endif 622 : c endif 623 : c bp (18/112010): cette strategie pose pb lors des permiers appels a 624 : c simu21. on doit pouvoir eviter de rencontrer le cas du nombre de bon 625 : c modes qui diminue autrement (grace au nbre d inconnues independantes) 626 : 627 : NOMBRV = NUMV 628 : * on a atteint nvectm vecteurs sans converger: peut on augmenter nvectm? 629 : if(numv.ge.nvectm.and.nbonve.lt.NMOD) 630 : $ nvectm=min(int(1.5*nvectm),nvecma) 631 : 632 : * on n a pas trouvé 1 seul mode: vraiment pas de chance (ou mauvais shift?) 633 : * => on essaie de redémarrer avec un nouveau vecteur 634 : c if (nbonve.lt.1.and.quit1.and.iquit1.lt.2) then 635 : if (nbonve.lt.1.and.quit1.and.iquit1.lt.2) then 635 : call aleat1 (IPRIGI,IPV1) 635 : CALL MUCPRI (IPV1,IPMASS,IPMV1) 636 : c ajout bp 17.12.2010 637 : write(6,*) 'lanczo: appel a initfl' 637 : CALL INITFL (IPRIGI,IPMASS,IPMV1,IPV1,IFLU) 638 : endif 639 : 640 : c on a usé la 2eme chance : on ne doit plus iterer 641 : if(quit1.and.iquit1.ge.2) goto 110 642 : 643 : * on quitte le processus lorsqu on a convergé ou atteint nvectm vecteurs 644 : * convergence? 645 : CONV = (CONV.or.(nbonve.ge.NMOD)) 646 : 647 : * si on a convergé sur NMOD et on a quasiment atteint le nombre total de 648 : * modes Nbmode, alors on n a pas interet a cycler 649 : PRESQ = (nbonve.ge.int(0.8*Nbmode)).and.(nbonve.lt.Nbmode) 650 : if(CONV.and.PRESQ.and.numv.lt.(0.8*nvectm)) goto 99 651 : 652 : c on a cvgé pour ce cycle ou atteint de nombre max de vecteurs : 653 : c on ne peut plus iterer 654 : c if(quit1.or.CONV.or.numv.ge.nvectm) goto 110 655 : if(CONV.or.numv.ge.nvectm) goto 110 656 : 657 : 99 continue 658 : NMBRV1 = NUMV + (int(0.1*NMBRV1) + 1) 659 : 660 : ENDIF 661 : 662 : 100 continue 663 : *-----fin de la Boucle iterative --------------------------------------* 664 : * 665 : 110 continue 666 : 667 : if (iimpi.ge.6) then 668 : WRITE(6,*) '----------------- Lanczos -----------------' 669 : WRITE(IOIMP,1990) NUMV,nvecma,NBonve,NMOD 670 : 1990 FORMAT(2X,'NOMBRE DE VECTEURS DE BASE',I5,' /',I5,' MAX' 671 : C ,/2X,'NOMBRE DE MODES TROUVES',I5,' /',I5,' DEMANDES') 672 : endif 673 : 674 : * 675 : *----- MENAGE --------------------------------------* 676 : * 677 : c * -on desactive x^j et r^j = M*z^j / (z^j*M*z^j) 678 : c if (ortho) then 679 : c do j1=1,Northo 680 : c jpo=mlen10.lect(j1) 681 : c mvect6= mlent6.lect(jpo) 682 : c mvect7= mlent7.lect(jpo) 683 : c segdes,mvect6,mvect7 684 : c enddo 685 : c segdes,mlent6,mlent7 686 : c segsup,mlen10 687 : c endif 688 : 689 : * -on ne garde que les listentier vers les chpoints/vecteurs X et MX 690 : jg=numv 691 : segadj,mlent3 692 : * on detruit ici les alpha*s(IPTTA=iplas), beta*s(ipttb=iplbs) 692 : CALL DTLREE(iplas) 692 : CALL DTLREE(iplbs) 692 : CALL DTLREE(ipla) 692 : CALL DTLREE(iplb) 692 : CALL DTLREE(iplc) 693 : 694 : * -on supprime les vecteurs MX et KX 695 : do iou=1,mlent4.lect(/1) 696 : mvect1=mlent4.lect(iou) 697 : segsup,mvect1 698 : enddo 699 : segsup,mlent4 700 : 701 : * on passe mlent3 a travers IPVECX 702 : IPVECX=mlent3 703 : 704 : segsup,idemen 705 : 705 : call gibtem(xkt) 706 : xktlan = xktlan + xkt 707 : if(iimpi.ge.6) write(6,*) 'temps passe dans Lanczos=',xktlan 708 : xktale = xktlan 709 : 710 : RETURN 711 : END 712 : 713 : 714 : 715 : 716 : 717 : 718 : 719 : 720 : 721 : 722 : 723 :
© Cast3M 2003 - All rights reserved.
Disclaimer