Télécharger lanczo.eso

Retour à la liste

Numérotation des lignes :

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

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