Télécharger lanczo.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales