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    MB234859  26/06/25    21:15:16     12580          
  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.          INSYM=0
  277.          CALL RESOU1 (IPRIGI,IDEMEN,NOID,NOEN,xspetl,0,1,INSYM,0)
  278.          if (iflu.eq.1) then
  279.            mrigid=iprigi
  280.            segact mrigid
  281.            icholi=ichole
  282.            mvecri=ivecri
  283.            segdes mrigid
  284.          endif
  285.          SEGACT IDEMEN*mod
  286.          IPV2 = IDEMEN(1)
  287.          call dtchpo(ipmv1)
  288.  
  289. *        y^i: passage chpo -> vecteur
  290.          call chv3(icholi,ipv2,mvect3,1)
  291.          call dtchpo(ipv2)
  292.          segact mvect3*mod
  293.  
  294. *       -Application des coef beta et alpha
  295. *        y^i = y^i - beta_{i-1} * x^{i-1}
  296.          if (numv.gt.2) then
  297.            mvect2=mlent3.lect(numv-2)
  298.            b1s0=b1*s0
  299.            do iou=1,inc
  300.               mvect3.vectbb(iou) = mvect3.vectbb(iou)
  301.      $                           - (b1s0 * mvect2.vectbb(iou))
  302.            enddo
  303.          endif
  304.  
  305. *        y^i = y^i - alpha_{i} * x^{i}
  306.          if (numv.gt.1) then
  307.            mvect1 = mlent3.lect(numv1)
  308.            mvect4 = mlent4.lect(numv1)
  309.            a1s1=a1*s1
  310.            do iou=1,inc
  311.               mvect3.vectbb(iou) = mvect3.vectbb(iou)
  312.      $                           - (a1s1 * mvect1.vectbb(iou))
  313.            enddo
  314.          endif
  315.          s0=s1
  316.  
  317. *       -calcul de K * y^i
  318.          call crech2(ipv1,mvect3,mvecri,1)
  319.          call mucpri(ipv1,IPRIGI,ipkv1)
  320.          call chv2 (icholi,ipkv1,mvect4,1)
  321.          segact,mvect3*mod
  322.          segact,mvect4
  323. *        comme on n'utilise que mvect3 et mvect4 , on detruit:
  324.          call dtchpo(ipv1)
  325.          call dtchpo(ipkv1)
  326.  
  327. *       -calcul de coe = y^i * K * y^i
  328.          coe=0.D0
  329.          do iou=1,inc
  330.             coe= coe + mvect4.vectbb(iou)*mvect3.vectbb(iou)
  331.          enddo
  332.          segsup,mvect4
  333.  
  334.          unval = sqrt(abs(1.d0/coe))
  335.          s1    = abs(1.d0/coe)*coe
  336.          xtrx  = sign(1.d0,coe)
  337.          xtrx1 = xtrx
  338.  
  339. *       -on commence par K-normaliser y^i  (est-ce vraiment utile ?)
  340. c          if ((coe.lt.XORTHO).and.(numv.ne.1)) then
  341. c             NV = numv
  342. c             WRITE(IOIMP,2000) NV,coe
  343. c 2000        FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE',
  344. c      C      /3X,'DIFFICULTE POUR NORMALISER  LE ',I5,
  345. c      C      ' IEME VECTEUR ',/3X,'NORME PAR',
  346. c      C      ' RAPPORT A LA NORME INITIALE ',E12.5)
  347. c *           on force l algo a faire les tests de convergence
  348. c *           si on a trouvé des modes -> on quitte lanczos
  349. c *           sinon on redémarre avec un nouveau vecteur (permis 1 fois)
  350. c             quit1 = .true.
  351. c             iquit1 = iquit1 + 1
  352. c          endif
  353.          do iou=1,inc
  354.             mvect3.vectbb(iou) = mvect3.vectbb(iou)*unval
  355.          enddo
  356.  
  357. c        xktre1 = TCPU pour resoudre + normaliser (lineaire avec numv)
  358.          if (IB100.eq.itest) then
  359.             call gibtem(xktre1)
  360.             xktlan = xktlan + xktre1
  361.          endif
  362.  
  363.          i08 = 0
  364.          nort = 0
  365.    08    continue
  366.          i08 = i08 + 1
  367.          xtrum = 0.D0
  368.          xtrum2 = 0.D0
  369.  
  370. *       -K-orthogonalisation % aux précedents (si ce n est pas le 1er passage)
  371.          if (numv.gt.1) then
  372.            xtrx1 = xtrx
  373.            do ipo=1,numv1
  374. *            on recupere  x^j  ,  r^j = K*x^j  et s^j = sign(x^j*K*x^j)
  375.              mvect1 = mlent3.lect(ipo)
  376.              mvect2 = mlent4.lect(ipo)
  377.              utru   = iplc.prog(ipo)
  378. *            on calcule xtru = y^i * r^j = y^i*K*x^j
  379.              xtru = 0.D0
  380.              do iou=1,inc
  381.                 xtru = xtru + (mvect2.vectbb(iou)*mvect3.vectbb(iou))
  382.              enddo
  383.              xxtru = -1.D0 *xtru/utru
  384.              xtrum = xtrum + abs(xxtru)
  385.              xtrum2 = max(xtrum2,(abs(xxtru)))
  386. *            re-K-orthogonalisation selective (=si correction xxtru > 5.D-12)
  387. c              if(abs(xxtru).lt.1.D-12) go to 18
  388.              if(abs(xxtru).lt.XORTHO) go to 18
  389.              nort = nort + 1
  390. *            y^i = y^i - (xtru/utru)*x^{j}
  391.              do iou=1,inc
  392.                 mvect3.vectbb(iou) = mvect3.vectbb(iou)
  393.      $                             + (xxtru*mvect1.vectbb(iou))
  394.              enddo
  395.              xtrx1=xtrx1-((xtru*xtru)/utru)
  396.    18        continue
  397.           enddo
  398.          endif
  399.  
  400.          if(iimpi.ge.7) write(6,*) xtrum,xtrum2
  401. c c          if(xtrum.ge.XORTHO) goto 08 => explose
  402. c          if(xtrum.ge.XPREC21) goto 08
  403. c c          if(xtrum.ge.rlim) goto 08
  404. c c          if(xtrum.ge.0.707) goto 08
  405. c c        bp: l'algo redémarrerait-il avec un nouveau vecteur?
  406. c bp29/10/2010: avec XPREC21 comme precision, on constate qu'1 ou 2 boucle
  407. c d'orthogonalisation suffit, mais pour eviter de tourner en rond,
  408. c on limite a 2 passages
  409.          if (xtrum.ge.XPREC21) then
  410.             if(i08.lt.2)  goto 08
  411.             WRITE(IOIMP,1088) numv,xtrum
  412. 1088        FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE',
  413.      C      /3X,'DIFFICULTE POUR ORTHOGONALISER LE',I5,' IEME VECTEUR',
  414.      C      /3X,'SOMME DES PRODUITS SCALAIRE ',E12.5)
  415. *           on force l algo a faire les tests de convergence
  416. *           si on a trouvé des modes -> on quitte lanczos
  417. *           sinon on redémarre avec un nouveau vecteur (permis 1 fois)
  418.             quit1 = .true.
  419.             iquit1 = iquit1 + 1
  420.          endif
  421.  
  422.  
  423. c *       -M-orthogonalisation % aux modes deja convergés
  424. c          if (ortho) then
  425. c c           JG0 = mlent6.LECT(/1)
  426. c c           do jpo=1,JG0
  427. c           do j1=1,Northo
  428. c              jpo=mlen10.lect(j1)
  429. c *            on recupere  x^j  et  r^j = M*z^j / (z^j*M*z^j)
  430. c              mvect6= mlent6.lect(jpo)
  431. c              mvect7= mlent7.lect(jpo)
  432. c c deja fait             segact mvect6,mvect7
  433. c *            on calcule xtru = y^i * r^j
  434. c              xtru = 0.D0
  435. c              do iou=1,inc
  436. c                 xtru= xtru + (mvect7.vectbb(iou)*mvect3.vectbb(iou))
  437. c              enddo
  438. c c              if(abs(xtru).lt.1.D-12) go to 28
  439. c              if(abs(xtru).lt.XORTHO) go to 28
  440. c *            y^i = y^i - xtru * z^{j}  et idem pour Ky^i
  441. c              do iou=1,inc
  442. c                 mvect3.vectbb(iou) = mvect3.vectbb(iou)
  443. c      $                             - (xtru*mvect6.vectbb(iou))
  444. c              enddo
  445. c    28        continue
  446. c c              segdes mvect6,mvect7
  447. c c            mis en commentaire car on reothogonalise chaque vecteur x
  448. c           enddo
  449. c c           segdes,mlent6,mlent7
  450. c          endif
  451. c
  452. c bp 27/10/2010 : finalement, on ne reorthogonalise que le vecteur y^0
  453. c                 car semble etre trop contraignant pour lanczo
  454.  
  455.  
  456. c        xktort = TCPU pour orthogonaliser (proport à numv*numv1)
  457.          if (IB100.eq.itest) then
  458.             call gibtem(xktort)
  459.             xktlan = xktlan + xktort
  460.          endif
  461.  
  462. *       -re-calcul de r^{i+1} = K*y^i
  463.          call crech2(ipv1,mvect3,mvecri,1)
  464.          call mucpri(ipv1,IPRIGI,ipkv1)
  465.          call chv2 (icholi,ipkv1,mvect4,1)
  466.          segact,mvect3*mod
  467.          segact mvect4*mod
  468.  
  469. *         calcul de xtrx2 = y^i * K * y^i apres re-orthog
  470.           xtrx2 = 0.D0
  471.           do iou=1,inc
  472.             xtrx2 = xtrx2 + (mvect3.vectbb(iou)*mvect4.vectbb(iou))
  473.           enddo
  474.  
  475. *        -on verifie que l orthogonalisation ne change pas trop le vecteur
  476. c         (ce qui devrait etre vrai en arithmétique exacte pour les vecteurs non redémarrés)
  477.           rnm = abs(xtrx2/xtrx)
  478.           if ((rnm.lt.rlim).and.(numv.ne.1)) then
  479.             NV = numv
  480.             WRITE(IOIMP,1000) NV,RNM
  481. 1000        FORMAT(/3X,'OPERATEUR VIBRA OPTION SIMULTANE',
  482.      C      /3X,'DIFFICULTE POUR ORTHOGONALISER LE ',I5,
  483.      C      ' IEME VECTEUR ',/3X,'NORME APRES ORTHOGONALISATION PAR',
  484.      C      ' RAPPORT A LA NORME INITIALE ',E12.5)
  485. *           on force l algo a faire les tests de convergence
  486. *           si on a trouvé des modes -> on quitte lanczos
  487. *           sinon on redémarre avec un nouveau vecteur (permis 1 fois)
  488.             quit1 = .true.
  489.             iquit1 = iquit1 + 1
  490.           endif
  491.  
  492. *        -K-normalisation :
  493.           coef=1.D0/sqrt(abs(xtrx2))
  494.           call mulpar(mvect3,mvect4,IPV1,IPKV1,coef)
  495.           call DTCHPO(IPKV1)
  496.  
  497.           s1=sign(1.D0,xtrx2)
  498.  
  499. *       -on recrée les chpoints x^{i+1} et K*x^{i+1} à partir des vecteurs modifiés
  500. c
  501. c        -evaluation du TCPU + optimisation du nbre de vecteurs utilisés
  502.           if (IB100.eq.itest) then
  503. c           xktres = TCPU pour resoudre + normaliser (lineaire avec numv)
  504.             call gibtem(xktre2)
  505.             xktlan = xktlan + xktre2
  506.             xktres = xktre0 + xktre1 + xktre2
  507.             itest = int(1.2*itest) + 1
  508.  
  509.             if (IB100.eq.itest0) then
  510.               yktort = xktort / numv1
  511.             else
  512.               yktort = 0.5*(yktort + (xktort / numv1))
  513.             endif
  514.  
  515. c           si xktort pas tres representatif, on le remplace
  516.             if(xktort.le.1.)   xktort = numv1*yktort
  517. c
  518. c           si orthog. ne coute pas cher, on la retestera + tard
  519.             if (xktort.lt.0.6.or.IB100.eq.itest0) then
  520. *
  521.              if (iimpi.ge.7) then
  522.              WRITE(6,1008) IB100,xktale,xktres,xktort,nort,xopt,Nmodopt
  523. 1008         FORMAT('IB100,xktale,xktres,xktort,nort,xopt,Nmodopt =',
  524.      C       I5,' ',E10.5,' ',E10.5,' ',E10.5,' ',I5,' ',E10.5,' ',I5)
  525.              endif
  526.  
  527. c           sinon calcul du nombre de vecteurs optimal %TCPU (=nopt)
  528.             else
  529. c bp : on va moyenner yktort pour le rendre moins sensible et variable
  530.               xkt0 = (1.*xktale) + (8.*xktres) + (32.*yktort)
  531.               xkt2 = (4.*yktort) + (2.*xktort/real(nort+numv1))
  532. *      xopt = sqrt(xkt0/xkt2) -> prise en compte du rebouchage de trou dans 50% des cas (*1.5)
  533.               xopt = 1.5*sqrt(xkt0/xkt2)
  534.               Nmodopt = INT(xopt) + 1
  535.               Nvecopt = (2*Nmodopt) +  8
  536.  
  537.               if (iimpi.ge.7) then
  538.              WRITE(6,1009) IB100,xktale,xktres,xktort,nort,xopt,Nmodopt
  539. 1009         FORMAT('IB100,xktale,xktres,xktort,nort,xopt,Nmodopt =',
  540.      C       I5,' ',E10.5,' ',E10.5,' ',E10.5,' ',I5,' ',E10.5,' ',I5)
  541.               endif
  542.  
  543. c             si on pensait ne faire qu'1 cycle et que l on veut moins
  544. c             de modes que l optimum, alors on n a pas interet a cycler
  545.               if (Nbmode.eq.Nmod.and.Nmodopt.gt.int(Nbmode*0.8)) then
  546.                 Nmodopt = Nmod
  547.                 Nvecopt = NOMBRV
  548. c             sinon on propose de cycler avec les valeurs optimisées:
  549.               else
  550.                 Nmodopt = min(Nmodopt,Nvecma)
  551. cbp: l'ajustement du nombre de modes calculés ne concerne que le 1er cycle*
  552. c    sinon pb lorsque l on bouche les trous....
  553.                 if (icycle.eq.1) then
  554. c               il n y a pas d interet a calculer + de mode que demandé
  555.                 Nmod = min(Nmod,Nmodopt)
  556.                 Nvecopt= min(Nvecopt,Nvecma)
  557.                 NOMBRV = Nvecopt
  558.                 NMBRV1 = int(0.8*NOMBRV)
  559.       if (iimpi.ge.6) then
  560.        write(6,*) ' nbre vect estimés necessaires',NOMBRV
  561.        write(6,*) ' pour le calcul de ',NMOD,' modes'
  562.        write(6,*) ' on teste la convergence tous les ',NMBRV1,' itérés'
  563.       endif
  564.                 endif
  565.               endif
  566.  
  567.             endif
  568.           endif
  569.  
  570. c        initialisation du TCPU
  571.          if (IB100.eq.(itest-1)) then
  572.             call gibtem(xkt)
  573.             xktlan = xktlan + xkt
  574.           endif
  575.  
  576. *        on evite les erreurs d'arrondi pour s^{i+1}
  577.          s1=sign(1.d0,s1)
  578.  
  579. *        on stocke x^{i+1} et  r^{i+1} = K*x^{i+1}
  580.          if (numv.ge.mlent3.lect(/1)) then
  581.            jg= (numv+1) *1.2
  582.            segadj mlent3,mlent4,ipla,iplb,iplc,iplas,iplbs
  583.          endif
  584.          mlent3.lect(numv)=mvect3
  585.          mlent4.lect(numv)=mvect4
  586.          iplc.prog(numv)=s1
  587.  
  588. *  calcul de a1 (alpha_i) et b1 (beta_i)
  589.          call mucpri(ipv1,ipmass,ipmv1)
  590.          call DTCHPO(IPV1)
  591.          call chv2(icholi,ipmv1,mvect5,1)
  592.          segact,mvect5
  593.          a1=0.d0
  594.          b1=0.d0
  595. *        alpha_i = y^i*K*x^i = x^i*M*K^-1*K*x^i = x^i*M*x^i
  596.          do iou=1,inc
  597.             a1=a1 + mvect3.vectbb(iou)*mvect5.vectbb(iou)
  598.          enddo
  599.          ipla.prog(numv)=a1
  600. *        alpha_i * s_i = (x^i * M * x^i) * sign(y^i * K * y^i)
  601. *        on ajoute le signe s_i car x^i = yb/ |yb*K*yb|^0.5
  602. *        donc x^i*K*x^i = yb*K*yb / |yb*K*yb|^0.5 = sign(yb*K*yb)
  603.          dp = a1*s1
  604.          a1s1=dp
  605.          iplas.prog(numv)=dp
  606. *        beta_i = y^i*K*x^{i-1} = x^{i}*M*x^{i-1}
  607.          if (numv.gt.1) then
  608.            mvect1 = mlent3.lect(numv1)
  609.            do iou=1,inc
  610.              b1=b1 + mvect1.vectbb(iou)*mvect5.vectbb(iou)
  611.            enddo
  612.            dp=b1*b1 *s0*s1
  613.            iplbs.prog(numv1)=dp
  614.            iplb.prog(numv1)=b1
  615.          endif
  616.          segsup,mvect5
  617.  
  618.  
  619. c        xktres = TCPU pour resoudre + normaliser (lineaire avec numv)
  620.          if (IB100.eq.(itest-1)) then
  621.             call gibtem(xktre0)
  622.             xktlan = xktlan + xktre0
  623.          endif
  624.  
  625. *       -Resolution du probleme projete : [T]*Phi = 1/lambda*[D]*Phi
  626. *        avec D =X*K*X = diag(s_i)  = [s_1 0 ... ; 0 s_2 0 ... ; ...]
  627. *        et T = X*M*X = [a_1 b_1 0  0 ... ; b_1 a_2 b_2 0 ... ; ...]
  628.          IF ((quit1).or.(NUMV .GE. nmbrv1)) THEN
  629.  
  630. c *        ajout bp 22/10/2010
  631. c *        calcul de xKTKx2 = ||r^{i+1}||_2^2 = ||K*x^{i+1}||_2^2
  632. c          xKTKx2 = 0.D0
  633. c          do iou=1,inc
  634. c             xKTKx2 = xKTKx2 + mvect4.vectbb(iou)*mvect4.vectbb(iou)
  635. c          enddo
  636. c          write(6,*) 'lanczo : ||K*x^{i+1}||_2^2=',xKTKx2
  637. c *        fin ajout bp 22/10/2010
  638.  
  639.            call simu21(ipla,iplb,iplc,NMOD,conv,nbonve,numv1,
  640.      $                 iplval,iphi,IPERM,XPREC21)
  641.  
  642.            if (iimpi.ge.6) then
  643.              write(6,*) 'lanczo: CONV=',CONV,'  car'
  644.              write(6,*) 'pour ',numv,'vecteurs /',nvectm,'maxi'
  645.              write(6,*) 'on a ',nbonve,'valeurs ok /',NMOD
  646.            endif
  647.  
  648. c c          si le nombre de bon modes diminue trop,
  649. c c          c'est qu'on aurait (probablement) du s'arreter plus tot ...
  650. c            if (nbonve.ge.nbonmax) then
  651. c               nbonmax= nbonve
  652. c               numbon = numv1
  653. c            endif
  654. c            if (nbonmax.gt.0) then
  655. c              nbonmin=0.7*nbonmax
  656. c              nbonmin=max(2,nbonmin)
  657. c              if (nbonve.le.nbonmin) then
  658. c                call simu21(ipla,iplb,iplc,NMOD,conv,nbonve,numbon,
  659. c      $                     iplval,iphi,IPERM,XPREC21)
  660. c                goto 110
  661. c              endif
  662. c            endif
  663. c bp (18/112010): cette strategie pose pb lors des permiers appels a
  664. c simu21. on doit pouvoir eviter de rencontrer le cas du nombre de bon
  665. c modes qui diminue autrement (grace au nbre d inconnues independantes)
  666.  
  667.            NOMBRV = NUMV
  668. * on a atteint nvectm vecteurs sans converger: peut on augmenter nvectm?
  669.            if(numv.ge.nvectm.and.nbonve.lt.NMOD)
  670.      $           nvectm=min(int(1.5*nvectm),nvecma)
  671.  
  672. * on n a pas trouvé 1 seul mode: vraiment pas de chance (ou mauvais shift?)
  673. * => on essaie de redémarrer avec un nouveau vecteur
  674. c            if (nbonve.lt.1.and.quit1.and.iquit1.lt.2) then
  675.            if (nbonve.lt.1.and.quit1.and.iquit1.lt.2) then
  676.               call aleat1 (IPRIGI,IPV1)
  677.               CALL MUCPRI (IPV1,IPMASS,IPMV1)
  678. c             ajout bp 17.12.2010
  679.               write(6,*) 'lanczo: appel a initfl'
  680.               CALL INITFL (IPRIGI,IPMASS,IPMV1,IPV1,IFLU)
  681.            endif
  682.  
  683. c          on a usé la 2eme chance : on ne doit plus iterer
  684.            if(quit1.and.iquit1.ge.2)   goto 110
  685.  
  686. * on quitte le processus lorsqu on a convergé ou atteint nvectm vecteurs
  687. *          convergence?
  688.            CONV = (CONV.or.(nbonve.ge.NMOD))
  689.  
  690. * si on a convergé sur NMOD et on a quasiment atteint le nombre total de
  691. * modes Nbmode, alors on n a pas interet a cycler
  692.            PRESQ = (nbonve.ge.int(0.8*Nbmode)).and.(nbonve.lt.Nbmode)
  693.            if(CONV.and.PRESQ.and.numv.lt.(0.8*nvectm))  goto 99
  694.  
  695. c          on a cvgé pour ce cycle ou atteint de nombre max de vecteurs :
  696. c          on ne peut plus iterer
  697. c            if(quit1.or.CONV.or.numv.ge.nvectm)   goto 110
  698.            if(CONV.or.numv.ge.nvectm)   goto 110
  699.  
  700.    99      continue
  701.            NMBRV1 = NUMV + (int(0.1*NMBRV1) + 1)
  702.  
  703.          ENDIF
  704.  
  705.   100 continue
  706. *-----fin de la Boucle iterative --------------------------------------*
  707. *
  708.   110 continue
  709.  
  710.       if (iimpi.ge.6) then
  711.         WRITE(6,*) '----------------- Lanczos -----------------'
  712.         WRITE(IOIMP,1990) NUMV,nvecma,NBonve,NMOD
  713. 1990    FORMAT(2X,'NOMBRE DE VECTEURS DE BASE',I5,' /',I5,' MAX'
  714.      C  ,/2X,'NOMBRE DE MODES TROUVES',I5,' /',I5,' DEMANDES')
  715.       endif
  716.  
  717. *
  718. *----- MENAGE  --------------------------------------*
  719. *
  720. c *    -on desactive  x^j  et  r^j = M*z^j / (z^j*M*z^j)
  721. c       if (ortho) then
  722. c        do j1=1,Northo
  723. c           jpo=mlen10.lect(j1)
  724. c           mvect6= mlent6.lect(jpo)
  725. c           mvect7= mlent7.lect(jpo)
  726. c           segdes,mvect6,mvect7
  727. c        enddo
  728. c        segdes,mlent6,mlent7
  729. c        segsup,mlen10
  730. c       endif
  731.  
  732. *    -on ne garde que les listentier vers les chpoints/vecteurs X et MX
  733.       jg=numv
  734.       segadj,mlent3
  735. *     on detruit ici les alpha*s(IPTTA=iplas), beta*s(ipttb=iplbs)
  736.       CALL DTLREE(iplas)
  737.       CALL DTLREE(iplbs)
  738.       CALL DTLREE(ipla)
  739.       CALL DTLREE(iplb)
  740.       CALL DTLREE(iplc)
  741.  
  742. *    -on supprime les vecteurs MX et KX
  743.       do iou=1,mlent4.lect(/1)
  744.        mvect1=mlent4.lect(iou)
  745.        segsup,mvect1
  746.       enddo
  747.       segsup,mlent4
  748.  
  749. *     on passe mlent3 a travers IPVECX
  750.       IPVECX=mlent3
  751.  
  752.       segsup,idemen
  753.  
  754.       call gibtem(xkt)
  755.       xktlan = xktlan + xkt
  756.       if(iimpi.ge.6)  write(6,*) 'temps passe dans Lanczos=',xktlan
  757.       xktale = xktlan
  758.  
  759.       RETURN
  760.       END
  761.  
  762.  
  763.  
  764.  
  765.  
  766.  
  767.  
  768.  
  769.  
  770.  
  771.  
  772.  
  773.  
  774.  
  775.  
  776.  
  777.  
  778.  
  779.  
  780.  
  781.  
  782.  

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