Télécharger la plaquette Cast3M


fr en

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