Télécharger lanczo.eso

Retour à la liste

Numérotation des lignes :

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

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