Télécharger la plaquette Cast3M


fr en

Télécharger simul3.eso

Retour à la liste

Numérotation des lignes : 

simul3

  1. C SIMUL3    SOURCE    PV        17/12/20    21:15:40     9674           
  2.       SUBROUTINE SIMUL3 (IPKW2M,IPMASS,IPVALP,IPHI,IPERM,NMOD,nbonve,
  3.      $                   Nramax,IPSOLU,LIMAGE)
  4. c      $                   nbonZ,nrate,IPSOLU,LIMAGE)
  5. ***********************************************************************
  6. *
  7. *                             S I M U L 3
  8. *                             -----------
  9. *
  10. * FONCTION:
  11. * ---------
  12. *
  13. *     APPELE PAR LE SOUS-PROGRAMME "SIMUL1".
  14. *     CALCUL DES VECTEURS PROPRES.
  15. *     (VIBRATIONS - OPTION SIMULTANE)
  16. *
  17. *
  18. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE
  19. * -----------
  20. *
  21. *     IPKW2M  (E)  POINTEUR DE L'OBJET 'RIGIDITE' K-w^2*M
  22. *     IPMASS  (E)  POINTEUR DE L'OBJET 'RIGIDITE' M
  23. *     IPVALP  (E)  POINTEUR DE L'OBJET 'LISTREEL' CONTENANT LA SUITE
  24. *                  DE VALEURS PROPRES DU PETIT PROBLEME AUX VALEURS
  25. *                  PROPRES.
  26. *     IPHI    (E)  POINTEUR DE L'OBJET MATRIX des coef des vect propres
  27. *     IPERM   (E)  LISTENTIER pour ordonner les vect p comme les val p
  28. *     NMOD    (E)  NOMBRE DE MODES PROPRES DEMANDE.
  29. *     nbonve  (E)  NOMBRE DE MODES ESTIMES CONVERGES PAR SIMU21
  30. *     IPSOLU  (S)  POINTEUR DE L'OBJET 'SOLUTION' CONTENANT LA SUITE
  31. *                  DE MODES PROPRES SOLUTIONS.
  32. *     LIMAGE  (E)  LOGICAL =vrai si valeurs propres negatives
  33. *                  (matrice K non definie positive par ex.)
  34. *
  35. *     + PARAMETRES PASSES PAR LE COMMON "CSIMUL".
  36. *
  37. * DESCRIPTION DES VARIABLES : cf. commentaires au cours du programme
  38. * --------------------------
  39. *
  40. * MODE DE FONCTIONNEMENT:
  41. * -----------------------
  42. *     -DEDUCTION DES FREQUENCES PROPRES DU PROBLEME PHYSIQUE,
  43. *     -RECOMBINAISON DES VECTEURS PROPRES Z=X*PHI,
  44. *     -CALCUL DU RESIDU,
  45. *      si(cvgce): ORTHOGONALISATION % MODES DEJA CONVERGES
  46. *      sinon:     2 ITERATIONS INVERSES + ORTHOGONALISATION
  47. *     -CREATION OU COMPLETION DE L'OBJET "SOLUTION" IPSOLU via SIMUL6
  48. *
  49. * LANGAGE : ESOPE (FORTRAN77)
  50. *
  51. ************************************************************************
  52. *
  53.       IMPLICIT INTEGER(I-N)
  54.       IMPLICIT REAL*8(A-H,O-Z)
  55. -INC CCREEL
  56.  
  57. -INC PPARAM
  58. -INC CCOPTIO
  59. -INC SMLREEL
  60. -INC SMLENTI
  61. -INC SMVECTD
  62. -INC SMRIGID
  63. *
  64. *     REGOUPEMENT DE RENSEIGNEMENTS POUR "LANCZO" et "SIMUL.":
  65.       COMMON/CSIMUL/ W2SHIF,XPREC21,mvecri,IPVECX,IPVECZ,IPMZ,IPFREZ
  66.       REAL*8 W2SHIF,ZTMZ,XPREC21,XORTHO
  67.       REAL*8 DEUXPI
  68.       PARAMETER (DEUXPI = (2.D0*XPI))
  69.  
  70.       SEGMENT MATRIX
  71.        REAL*8 A(N,N)
  72.       ENDSEGMENT
  73.       POINTEUR IPHI.MATRIX
  74.       pointeur mlent5.mlenti,mlent6.mlenti,mlent7.mlenti
  75.       pointeur mvect5.MVECTD,mvect6.MVECTD,mvect7.MVECTD
  76.       pointeur mvect4.MVECTD
  77.  
  78.       PARAMETER (LPROPR = 5)
  79.       LOGICAL CONVRG
  80.       REAL*8 PROPRE(LPROPR)
  81. *
  82. *
  83.       REAL*8 VALPP,FREP,W2,W2CALC
  84. *
  85.  
  86. *
  87. ***********************************************************************
  88. *
  89.       LOGICAL LIMAGE
  90. *
  91. *---------------------------------------------------------------------*
  92. *     -- INITIALISATIONS --
  93. *---------------------------------------------------------------------*
  94. *
  95. c inutile car fait par simul1 lorsqu on shifte:     IPSOLU = 0
  96. *
  97.       MLREEL=IPVALP
  98.       SEGACT,MLREEL
  99.       NMODMAX= PROG(/1)
  100.       MLENTI=IPERM
  101.       SEGACT,MLENTI
  102.       SEGACT,IPHI
  103.       N=IPHI.A(/1)
  104. *     recup des vecteurs de Lanczos de ce cycle
  105.       mlent3 = IPVECX
  106.       NVEC3  = mlent3.lect(/1)
  107.       mvect3 = mlent3.lect(1)
  108.       inc = mvect3.VECTBB(/1)
  109.  
  110. *     recup des modes deja convergés aux précédents cycles (=IPVECZ)
  111.       if (IPVECZ.ne.0) then
  112.         mlent5 = IPVECZ
  113.         mlent6 = IPMZ
  114.         mlree1 = IPFREZ
  115.         segact,mlent5*mod
  116.         JG0 = mlent5.LECT(/1)
  117. c         JG = JG0 + nbonve
  118.         JG = JG0 + NMODMAX
  119.         segadj,mlent5
  120.         segadj,mlent6
  121.         segadj,mlree1
  122.         segact,mlent5*mod
  123.         segact,mlent6*mod
  124.         segact,mlree1*mod
  125. *     creation des modes convergés (=IPVECZ)
  126.       else
  127.         JG0 = 0
  128. c         JG = nbonve
  129.         JG = NMODMAX
  130.         segini,mlent5
  131.         IPVECZ = mlent5
  132.         segini,mlent6
  133.         IPMZ = mlent6
  134.         segini,mlree1
  135.         IPFREZ = mlree1
  136.       endif
  137.  
  138. c     on considere qu il faut reorthognaliser si le produit scalaire
  139. c     est > XORTHO (avec **1.6 XORTHO=3E-13 si xprec21=1.5E-8)
  140.       XORTHO = XPREC21**1.6
  141. c c     on ne fait meme pas d iterations inverses si residu trop grand
  142. c       XRMAX = sqrt(XPREC21)
  143. c => mauvaise idée pour dyne_devo_10.dgibi
  144. c     limite pour detecter la presence dans la solution
  145. c       RLIM = 1.D-5
  146.       RLIM = sqrt(XPREC21)
  147.  
  148.       IF (IIMPI .ge. 6) THEN
  149.       WRITE(IOIMP,*) '------------------- simul3 ---------------------'
  150.       WRITE(IOIMP,*) 'candidat| mode    |    freq    |  residu relatif'
  151.       ENDIF
  152.  
  153. c       -recup de icholi pour le passage vecteur <-> chpoint
  154.          mrigid = IPKW2M
  155.          segact,mrigid
  156.          icholi = ichole
  157.          segdes,mrigid
  158. c       call gibtem(xkt30)
  159. c       WRITE(IOIMP,*) '=> simul3        xkt30=',xkt30
  160.  
  161. c     shift utilise pour ce cycle
  162.       FRSHIF = (sqrt(abs(W2SHIF))) / DEUXPI
  163. *
  164. *---------------------------------------------------------------------*
  165. *     -- Boucle sur les candidats a devenir des modes --
  166. *---------------------------------------------------------------------*
  167. *
  168. c     Au min, On teste les NMOD demandés ou nbonve supposé convergés
  169.       NMODIN = MAX(nbonve,NMOD)
  170. c     Au max, On teste tous les candidats (tant que ca converge)
  171. c       N105 = MIN(nbonve,NMODMAX)
  172.       N105 = NMODMAX
  173. c     Nombre de modes ratés, convergés
  174. c       NRATE = 0
  175.       IEME = 0
  176.  
  177.       DO 105 IB105 = 1, N105
  178. *
  179.          IEME = IEME + 1
  180. *
  181. *------- RETOUR AU PROBLEME PHYSIQUE INITIAL -------------------------*
  182. *
  183. *       -IEME FREQUENCE PROPRE:
  184.          VALPP = PROG(IB105)
  185.          W2CALC = VALPP
  186.          CALL W2FREQ(W2CALC,W2SHIF, W2,FREP,LIMAGE)
  187.          IF(IERR .NE. 0) RETURN
  188. *
  189. *       -IEME VECTEUR PROPRE:
  190.          segini,mvect1
  191.          IEME2 = LECT(IB105)
  192.          DO jvec=1,N
  193.             XPHI = IPHI.A(jvec,IEME2)
  194.             mvect3 = mlent3.lect(jvec)
  195.             do i=1,inc
  196.                mvect1.VECTBB(i) = mvect1.VECTBB(i)
  197.      $                   + ((mvect3.VECTBB(i)) * XPHI)
  198.             enddo
  199.          ENDDO
  200.  
  201. *       -creation du chpoint
  202.          call crech2(IPVECP,mvect1,mvecri,1)
  203.          segact,mvect1*mod
  204. c       call gibtem(xkt30z)
  205. c       WRITE(IOIMP,*) '=> simul3  z=x*p xkt30=',xkt30z
  206.  
  207.          IF (IIMPI .EQ. 369) THEN
  208.             WRITE (IOIMP,*) 'simul3: VECTEUR PROPRE PHYSIQUE_',IEME
  209.             CALL ECCHPO (IPVECP,0)
  210.          ENDIF
  211.  
  212. *        on suppose qu on a convergé
  213.          CONVRG = .true.
  214.  
  215. *
  216. *------- CALCUL DU RESIDU --------------------------------------------*
  217. *
  218. *lorsqu on n est pas sur de la convergence on teste directement le residu
  219.  
  220. *parfois, il est necessaire d iterer un peu pour "purger" la solution
  221.          N101 = 30
  222. c          N101 = 2
  223.          I101 = 0
  224.   101    CONTINUE
  225.          I101 = I101 + 1
  226. *
  227. c       -calcul de [K-w^2M]*Z
  228.          call MUCPRI(IPVECP,IPKW2M,IPKV1)
  229.          call chv2(icholi,IPKV1,mvect4,1)
  230.          call dtchpo(IPKV1)
  231. c       -calcul de MZ
  232.          call MUCPRI(IPVECP,IPMASS,IPMV1)
  233.          call chv2(icholi,IPMV1,mvect2,1)
  234.  
  235. *       -calcul de (Z^T.M.Z)
  236.          ZTMZ = 0.D0
  237.          do i=1,inc
  238.             ZTMZ = ZTMZ + (mvect2.vectbb(i)*mvect1.vectbb(i))
  239.          enddo
  240.  
  241. c       -calcul du coef de Rayleigh ZRAY1 = Z^T K Z / Z^T M Z
  242. c       (est notamment utile lorsqu'on a fait appel a itinv)
  243.          IF (.NOT.CONVRG) THEN
  244.            ZNUM1 = 0.D0
  245.            ZDEN1 = ZTMZ
  246. c            ZDEN1 = 0.D0
  247.            do i=1,inc
  248.               ZNUM1 = ZNUM1 + (mvect1.VECTBB(i)*mvect4.VECTBB(i))
  249. c               ZDEN1 = ZDEN1 + (mvect1.VECTBB(i)*mvect2.VECTBB(i))
  250.            enddo
  251.            ZRAY1 = ZNUM1 / ZDEN1
  252.            VALPP = ZRAY1
  253.            W2CALC = VALPP
  254.            CALL W2FREQ(W2CALC,W2SHIF, W2,FREP,LIMAGE)
  255.          ENDIF
  256.  
  257. c        calcul du residu = Kz - lambda*Mz et de sa norme au carré
  258.          VALPM = -1.D0*VALPP
  259.          XNUM1 = 0.D0
  260.          do i=1,inc
  261.             XRES1 = mvect4.VECTBB(i) + (VALPM * mvect2.VECTBB(i))
  262.             XNUM1 = XNUM1 + (XRES1*XRES1)
  263.          enddo
  264.  
  265. c        calcul d une ref |K*z|_2^2 ou |M*z|_2^2
  266.          XDEN1=0.D0
  267.          XDEN2=0.D0
  268.          do i=1,inc
  269.             XDEN1 = XDEN1 + (mvect4.VECTBB(i)*mvect4.VECTBB(i))
  270.             XDEN2 = XDEN2 + (mvect2.VECTBB(i)*mvect2.VECTBB(i))
  271.          enddo
  272.          segsup,mvect4
  273. c          XRES1 = XNUM1 / XDEN1
  274. c          if(XRES1.gt.1.D-5)
  275.          XDEN1m=max(XDEN1,XDEN2)
  276.          XRES1 = XNUM1 / XDEN1m
  277.  
  278.          IF (IIMPI.ge.6) then
  279.           WRITE(IOIMP,109) IB105,IEME,frep,XRES1,XNUM1,XDEN1,XDEN2
  280.  109      FORMAT(1X,I5,3X,I5,3X,F12.3,5X,E10.5,
  281.      &            5X,E10.5,1X,E10.5,1X,E10.5)
  282.          endif
  283. *
  284. *------- DIFFERENT CAS DE FIGURES ------------------------------------*
  285. *
  286.  
  287. *------- a. on a convergé: le residu relatif < tolerance -------------*
  288. cbp13.01.2010     if (XRES1.le.1.D-4) then -> pas assez contraignant ?
  289. cbp25.01.2010     if (XRES1.le.1.D-6) then -> pas suffisant pour K-orthog
  290. c          if (XRES1.le.1.D-8)  then
  291.          if (XRES1.le.XPREC21)  then
  292.  
  293.                 CONVRG = .true.
  294. c       call gibtem(xkt31)
  295. c       WRITE(IOIMP,*) '=> simul3 CONVRG xkt31=',xkt31
  296.  
  297.  
  298. *------- b. le residu relatif > tolerance  ---------------------------*
  299. *           et on avait convergé en lambda
  300. c           =>  on essaie de converger en residu = "Purge"
  301. c          elseif((IEME.le.nbonve).and.(I101.le.N101)) then
  302.          elseif(I101.le.N101) then
  303. c          elseif((I101.le.N101).and.(XRES1.le.XRMAX)) then
  304.  
  305. *            -puisque X et MX vont changer ...
  306.               segsup,mvect1,mvect2
  307.  
  308. *            -iterations inverses avec z en vecteur initial
  309.               CALL ITINV (IPKW2M,IPMASS,IPVECP,PROPRE,CONVRG,1,100
  310.      &                    ,1.D-5,1.D0,IPMV1)
  311. * bp: avec des petites matrices (3x3 par ex), ITINV ne peut rien si on
  312. * ne shifte pas => il est urgent d'améliorer l'estimation des vp faite
  313. * par sespa3 dans simu21 (utilisation du QR par ex.)
  314. c               if(IIMPI.ge.6)
  315. c      &           WRITE(IOIMP,*) ' simul3: Residu=',XRES1,
  316. c      &          ' appel',I101,'a itinv cvg?',CONVRG,' vers',(PROPRE(1))
  317.               if(IERR .NE. 0) RETURN
  318.  
  319.               call chv3(icholi,IPVECP,mvect1,1)
  320.               call dtchpo(IPVECP)
  321.  
  322. *            -M-orthogonalisation % ...
  323.               I08 = 0
  324.    08         continue
  325.               I08 = I08 + 1
  326.               xtrum = 0.D0
  327.  
  328. c *             ... aux modes des cycles précedents convergés...
  329. c               if (JG0.gt.1)  then
  330. c                  if(IIMPI.ge.7)
  331. c      &           WRITE(IOIMP,*) ' simul3: M-orthog % cycles precedents'
  332. c                  do 09 jpo=1,JG0
  333. c                    FRjpo = mlree1.prog(jpo)
  334. c                    FRDIST = 1.001D0*(FREP-FRSHIF)
  335. c                    FRDISj = (FRjpo-FRSHIF)
  336. c                    FRSIGN = FRDIST*FRDISj
  337. c *                  ...proches du modes testés!
  338. c                    if((FRSIGN.lt.0.D0).or.(FRDIST.lt.FRDISj)) goto 09
  339. c *                  on recupere  z^j  et  r^j = M*z^j / (z^j*M*z^j)
  340. c                    mvect6= mlent5.lect(jpo)
  341. c                    mvect7= mlent6.lect(jpo)
  342. c                    segact mvect6,mvect7
  343. c *                  on calcule xtru = y^i * r^j
  344. c                    xtru = 0.D0
  345. c                    do iou=1,inc
  346. c                       xtru=xtru+(mvect1.vectbb(iou)*mvect7.vectbb(iou))
  347. c                    enddo
  348. c                    xtrum = xtrum + abs(xtru)
  349. c                    if (abs(xtru).ge.XORTHO) then
  350. c *                    y^i = y^i - xtru * z^{j}
  351. c                      do iou=1,inc
  352. c                        mvect1.vectbb(iou) = mvect1.vectbb(iou)
  353. c      $                              - (xtru*mvect6.vectbb(iou))
  354. c                      enddo
  355. c                    endif
  356. c                    segdes mvect6,mvect7
  357. c    09            continue
  358. c               endif
  359.  
  360. *             ... aux modes de ce cycle convergés
  361.               if (IEME.gt.1)  then
  362.                  if(IIMPI.ge.7)
  363.      &           WRITE(IOIMP,*) ' simul3: M-orthog % modes de ce cycle'
  364.      &                   ,(JG0+1),' a ',(JG0+IEME-1)
  365.                  do jpo=(JG0+1),(JG0+IEME-1)
  366. *                  on recupere  z^j  et  r^j = M*z^j / (z^j*M*z^j)
  367.                    mvect6= mlent5.lect(jpo)
  368.                    mvect7= mlent6.lect(jpo)
  369. cbp inutile car jpo>JG0    segact mvect6,mvect7
  370. *                  on calcule xtru = y^i * r^j
  371.                    xtru = 0.D0
  372.                    do iou=1,inc
  373.                       xtru=xtru+(mvect1.vectbb(iou)*mvect7.vectbb(iou))
  374.                    enddo
  375.                    xtrum = xtrum + abs(xtru)
  376.                    if (abs(xtru).ge.XORTHO) then
  377. *                    y^i = y^i - xtru * z^{j}
  378.                      do iou=1,inc
  379.                        mvect1.vectbb(iou) = mvect1.vectbb(iou)
  380.      $                              - (xtru*mvect6.vectbb(iou))
  381.                      enddo
  382.                    endif
  383. cbp inutile car jpo>JG0    segdes mvect6,mvect7
  384.                  enddo
  385.               endif
  386.  
  387. *             a t on besoin de re-orthogonaliser ?
  388.               if (xtrum.ge.0.707D0)  then
  389. *               2 re-orthogonalisation c'est suffisant !
  390.                 if (I08.ge.2) then
  391.                    call crech2(IPVECP,mvect1,mvecri,1)
  392.                    segact,mvect1*mod
  393.                    CONVRG = .false.
  394.                    goto 101
  395.                 endif
  396.                 goto 08
  397.               endif
  398. *            -fin de la M-orthogonalisation
  399.  
  400.               call crech2(IPVECP,mvect1,mvecri,1)
  401.               segact,mvect1*mod
  402.  
  403.               CONVRG = .false.
  404. c              call DTCHPO(IPMV1)
  405.               goto 101
  406.  
  407. *------- c. le residu relatif > tolerance  ---------------------------*
  408. *           et on N'avait PAS convergé en lambda
  409. c           =>  on S'ARRETE LA
  410.          else
  411.  
  412.                CONVRG = .false.
  413.  
  414.          endif
  415.  
  416. *------- Non convergence on fait le menage et on retourne ------------*
  417.          IF (.NOT.CONVRG) THEN
  418.                IEME = IEME - 1
  419. *             -On detruit vecteurs, chpoints
  420.                segsup,mvect1,mvect2
  421.                call DTCHPO(IPVECP)
  422.                call DTCHPO(IPMV1)
  423.                goto 905
  424. c c *              si on a testé plus de NMODIN candidats, on sort
  425. c c                if(IB105.ge.NMODIN) goto 905
  426. c c c *              si on teste un candidat appartenant aux NMOD 1ers modes,
  427. c c c *              on le compte si on l a raté
  428. c c c                if((nbonZ+IB105).le.NMOD) NRATE=NRATE+1
  429. c c                goto 105
  430. c *              si on a testé plus de NMODIN candidats, on sort
  431. c                if(IB105.ge.NMODIN) goto 905
  432. c                NRATE=NRATE+1
  433. c *              si on a raté trop de modes, on sors aussi
  434. c                if(NRATE.gt.Nramax) goto 905
  435. c                goto 105
  436.          ENDIF
  437.  
  438. *------- Convergence : on fait aussi le menage  ----------*
  439.          call DTCHPO(IPMV1)
  440. *
  441. *       -CE MODE A T IL DEJA ETE CALCULE ?
  442.          modif205 = 0
  443.    09    continue
  444.          xtrum = 0.D0
  445.          FRREF = max(FRSHIF,abs(FREP))
  446.          FRREF = max(FRREF,1.D0)
  447.          NB205 = JG0 + IEME - 1
  448.          do 205 jb205=1,NB205
  449.             FR205 = mlree1.prog(jb205)
  450.             ECARFR = abs(FREP-FR205)
  451.             ECAREL = ECARFR / FRREF
  452. *           on a trouvé un autre mode a la meme freq => on orthogonalise
  453. c bp 10/12/2010 : quand les shifts de 2 cycles sont tres differents,
  454. c                 on peut avoir de plus gros ecart que prévu....
  455.             if ((ECAREL.lt.1.D-2).or.(ECARFR.lt.1.D0))  then
  456. *              on recupere  z^j  et  r^j = M*z^j / (z^j*M*z^j)
  457.                mvect6= mlent5.lect(jb205)
  458.                mvect7= mlent6.lect(jb205)
  459.                if(jb205.le.JG0)  segact,mvect6,mvect7
  460. *              on calcule xtru = y^i * r^j
  461.                xtru = 0.D0
  462.                do iou=1,inc
  463.                   xtru=xtru+(mvect1.vectbb(iou)*mvect7.vectbb(iou))
  464.                enddo
  465.                xtrum = xtrum + abs(xtru)
  466.                if (abs(xtru).ge.XORTHO) then
  467.                  modif205 = modif205 + 1
  468.                  if(IIMPI.ge.6) WRITE(IOIMP,*)
  469.      &           'ORTHOGONALISATION PAR RAPPORT AU MODE',jb205
  470. *                y^i = y^i - xtru * z^{j}
  471.                  do iou=1,inc
  472.                    mvect1.vectbb(iou) = mvect1.vectbb(iou)
  473.      $                          - (xtru*mvect6.vectbb(iou))
  474.                  enddo
  475.                endif
  476.                if(jb205.le.JG0)  segdes,mvect6,mvect7
  477.             endif
  478.   205    continue
  479. *        si grosse modif, on reorthogonalise en 2 passes
  480. *         if(xtrum.ge.0.707) goto 09   aucun sens car z non-normalisé
  481. c       call gibtem(xkt32)
  482. c       WRITE(IOIMP,*) '=> simul3 orthog xkt32=',xkt32,modif205
  483.  
  484. *       -si on a modifie le candidat on verifie qu il en reste qq chose
  485.          if (modif205.gt.0) then
  486.            segsup,mvect2
  487. *          calcul de M.Z
  488.            call DTCHPO(IPVECP)
  489.            call crech2(IPVECP,mvect1,mvecri,1)
  490.            call MUCPRI(IPVECP,IPMASS,IPMV1)
  491.            call chv2(icholi,IPMV1,mvect2,1)
  492.            call dtchpo(IPMV1)
  493.            segact,mvect1,mvect2*mod
  494. *          calcul de Z^T.M.Z
  495.            ZTMZ2 = 0.D0
  496.            do i=1,inc
  497.               ZTMZ2 = ZTMZ2 + (mvect2.vectbb(i)*mvect1.vectbb(i))
  498.            enddo
  499.          else
  500.            ZTMZ2 = ZTMZ
  501.          endif
  502.  
  503.          RNM = abs(ZTMZ2 / ZTMZ)
  504.          if (RNM.lt.RLIM) then
  505. *           le mode candidat est deja présent dans la solution
  506. *           (il en reste moins de 1/10000ieme pour xprec21=1E-8)
  507. *           => on l oublie et on passe au suivant
  508.             if(IIMPI.GE.6)
  509.      &      WRITE(IOIMP,*) 'MODE DEJA PRESENT DANS LA SOLUTION',RNM
  510.             IEME = IEME - 1
  511.             goto 105
  512.          endif
  513.  
  514. *       -le mode candidat a survecu aux differents tests
  515. *        => on l'ajoute a la solution
  516.  
  517. *       -ON STOCKE  Z  et  M.Z/(Z^T.M.Z)  POUR L ORTHOGONALISATION
  518.          coe2 = 1.D0 / ZTMZ2
  519.          do i=1,inc
  520.            mvect2.vectbb(i) = coe2 * (mvect2.vectbb(i))
  521.          enddo
  522. *       -ON STOCKE  FREP,  Z  et  M.Z/(Z^T.M.Z)  POUR L ORTHOGONALISATION
  523.          mlree1.prog(JG0+IEME) = FREP
  524.          mlent5.lect(JG0+IEME) = mvect1
  525.          mlent6.lect(JG0+IEME) = mvect2
  526.  
  527. *
  528. *------- CREATION OU COMPLETION DE L'OBJET "SOLUTION" ----------------*
  529. * solution du cycle seulement! le reste est dans idist (cf strate.eso)
  530.          CALL SIMUL6 (FREP,IPVECP,IPSOLU)
  531.          IF (IERR .NE. 0) RETURN
  532.  
  533. *
  534.   105 CONTINUE
  535. *     fin de la boucle sur les modes
  536.  
  537.   905 CONTINUE
  538. *
  539. *---- finalement on trouvé nbonve modes de residu < XPREC21 ----------*
  540.       nbonve = IEME
  541.       JG = JG0 + nbonve
  542.       segadj,mlent5
  543.       segadj,mlent6
  544.       segadj,mlree1
  545.       IF (IIMPI.EQ.7) then
  546.         WRITE(IOIMP,*) 'simul3: finalement '
  547.         WRITE(IOIMP,*) ' Modes OK pour ce cycle, au total',nbonve,JG
  548.       ENDIF
  549.  
  550. c     on desactive les vecteurs de ce cycle
  551.       if (IEME.ge.1) then
  552.       do j=(JG0+1),JG
  553.          mvect1=mlent5.lect(j)
  554.          mvect2=mlent6.lect(j)
  555.          segdes,mvect1,mvect2
  556.       enddo
  557.       endif
  558.  
  559.       segdes,mlent5,mlent6
  560.       segdes,mlree1
  561.  
  562. *
  563. *----- MENAGE  --------------------------------------*
  564. *
  565.       SEGSUP,MLREEL,MLENTI,IPHI
  566. *
  567. *
  568.       END
  569.  
  570.  
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.  

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