Esope source files

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