Télécharger la plaquette Cast3M


fr en

Télécharger sespac.eso

Retour à la liste

Numérotation des lignes : 

sespac

  1. C SESPAC    SOURCE    CB215821  20/11/25    13:39:42     10792          
  2. C     SUBROUTINE SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS )
  3. ************************************************************************
  4. *
  5. *                                SESPAC
  6. *                             -----------
  7. *
  8. * FONCTION:
  9. * ---------
  10. *
  11. *     CONSTRUIT NBMOD ELEMENTS PROPRES EN ITERANT LE SOUS-ESPACE
  12. *     IPLVEC, JUSQU'A LA CONVERGENCE DE CELUI-CI. OBTENTION DE VALEURS PROPRES
  13. *     RELLES OU COMPLEXES.
  14. *
  15. * MODE D'APPEL:
  16. * -------------
  17. *
  18. *     CALL SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS )
  19. *
  20. * PARAMETRES:   (E)=ENTREE   (S)=SORTIE
  21. * -----------
  22. *
  23. *     IPLVAL  ENTIER    (S)  POINTEUR DE L'OBJET 'LISTREEL' CONTENANT
  24. *                            LA SUITE DE VALEURS PROPRES OBTENUES.
  25. *     IPLVEC  ENTIER    (E)  POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  26. *                            LE SOUS-ESPACE INITIAL.
  27. *     IPLVEC  ENTIER    (S)  POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT
  28. *                            LE SOUS-ESPACE FINAL. EN ORTHONORMALISANT
  29. *                            LES 'CHPOINT' DE CET ESPACE ON OBTIENT LES
  30. *                            VECTEURS PROPRES RECHERCHES.
  31. *     NBMOD   ENTIER    (E)  NOMBRE DE VECTEURS RECHERCHES. IPLVEC
  32. *                            CONTIENT PLUS QUE NBMOD 'CHPO', CAR CECI
  33. *                            PERMET DE CONVERGER PLUS RAPIDEMENT.
  34. *
  35. *     IPKW2M  ENTIER    (E)  POINTEUR SUR L'OBJET 'RIGIDITE' K
  36. *     IPMASS  ENTIER    (E)  POINTEUR SUR L'OBJET 'RIGIDITE' M
  37. *
  38. *
  39. * AUTEUR, DATE DE CREATION:
  40. * -------------------------
  41. *
  42. *   C. LE BIDEAU  09 / 2001  ( FORTRAN + ESOPE )
  43. *   MODIF     Benoit Prabel    Mars 2009
  44. *
  45. ***********************************************************************
  46.  
  47. *      SUBROUTINE SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS )
  48.       SUBROUTINE SESPAC(IPLVAR,IPLVAI,IPLVEC,NBMOD,IPKW2M,IPMASS)
  49.  
  50.       IMPLICIT INTEGER(I-N)
  51.       IMPLICIT REAL*8 (A-H,O-Z)
  52.  
  53.       DATA XZERO / 0.0D0 /
  54.  
  55. -INC PPARAM
  56. -INC CCOPTIO
  57. -INC  SMLCHPO
  58. -INC  SMLREEL
  59. -INC  SMCHPOI
  60. -INC  SMRIGID
  61.  
  62. ******
  63. *     -- CONSTANTES --
  64. ***
  65. *   nbre d'iteration de la boucle complete et de celle sur pb global (=puissance)
  66.       PARAMETER ( ITERMX0 = 10 )
  67.       PARAMETER ( ITERMX1 = 3 )
  68. *      PARAMETER ( ITERMX0 = 3 )
  69. *      PARAMETER ( ITERMX1 = 2 )
  70. *      PARAMETER ( ITERMX1 = 1 )
  71.  
  72. ******
  73. *     -- ARGUMENTS --
  74. ***
  75.       POINTEUR IPLVEC.MLCHPO
  76. *      IPLVER.MLCHPO,IPLVEI.MLCHPO
  77.       POINTEUR IPLVAR.MLREEL,IPLVAI.MLREEL
  78.       INTEGER  IPKW2M, IPMASS, NBMOD
  79.  
  80. ******
  81. *     -- VARIABLES LOCALES --
  82. ***
  83.  
  84.       POINTEUR IPLCH1.MLCHPO,IPLCH2.MLCHPO
  85.       POINTEUR ILVA1R.MLREEL,ILVA1I.MLREEL,ILVA2R.MLREEL,ILVA2I.MLREEL
  86. cbp       POINTEUR ILAMBR.XMATRI,ILAMBI.XMATRI,ILAMBD.XMATRI,IPZ.XMATRI
  87.       POINTEUR IPZ.XMATRI
  88.       POINTEUR ILAMBR.MLREEL,ILAMBI.MLREEL,ILAMBD.MLREEL
  89.       INTEGER ILDIM,NERR
  90.       LOGICAL CALCV,CVGVAL,CVGCHP
  91.  
  92. ****************************************************************
  93. *                    INITIALISATIONS                           *
  94. ****************************************************************
  95.  
  96. *     ouverture des vecteurs itérés
  97.       SEGACT IPLVEC*MOD
  98.       ILDIM = IPLVEC.ICHPOI( /1 )
  99.  
  100. *     recherche des listmots inconnue (IPLMOX) et duale (IPLMOY)
  101.       IPCHPO = IPLVEC.ICHPOI(1)
  102.       CALL MUCPRI (IPCHPO, IPKW2M, IPCX)
  103.       CALL CORRSP (IPKW2M, IPCHPO, IPCX, IPLMOX, IPLMOY)
  104.       CALL DTCHPO (IPCX)
  105. *      CALL SESBAS (IPLVEC, IPKW2M, IPMASS, IPLMOX, IPLMOY)
  106.  
  107. *     creation des valeurs propres (initialisées à 0)
  108.       JG = ILDIM
  109.       SEGINI,ILVA1R,ILVA1I,ILVA2R,ILVA2I
  110.       DO JVEC = 1,ILDIM
  111. *         ILVA1.PROG(JVEC) = 0.D0
  112.          ILVA2R.PROG(JVEC) = 0.D0
  113.          ILVA2I.PROG(JVEC) = 0.D0
  114.       ENDDO
  115.  
  116. *     creation des vecteurs propres (initialisées à IPLVEC)
  117.       N1 = ILDIM
  118.       SEGINI,IPLCH1
  119.       DO JVEC = 1, ILDIM
  120.          IPCHPO = IPLVEC.ICHPOI(JVEC)
  121.          CALL COPIE2 ( IPCHPO, IPCH1 )
  122.          IF ( IERR .NE. 0 ) RETURN
  123.          IPLCH1.ICHPOI(JVEC)  = IPCH1
  124.       ENDDO
  125.  
  126. *     impressions
  127. *      do JVEC=1,ILDIM
  128. *        MCHPOI = IPLVEC.ICHPOI(JVEC)
  129. *        segact,MCHPOI
  130. *        MSOUPO=IPCHP(1)
  131. *        segact,MSOUPO
  132. *        MPOVAL=IPOVAL
  133. *        segact,MPOVAL
  134. *        write(*,*) ' X^0_',JVEC,'=',VPOCHA(1,1),VPOCHA(1,2),VPOCHA(1,3)
  135. *        segdes,MPOVAL,MSOUPO,MCHPOI
  136. *      enddo
  137.  
  138.       SEGDES,IPLVEC,IPLCH1
  139.       IPLCH2 = IPLVEC
  140.       IPLVEC = 0
  141.  
  142. ****************************************************************
  143. *  ITERMX0  ITERATIONS associant iterations GLOBALES et sur systeme reduit
  144. *
  145. *  REPETER JUSQU'A: *  CONVERGER
  146. *                   *  DEPASSER ITERMX0 ITERATIONS
  147. ****************************************************************
  148.  
  149.       DO 100 IB100 = 1,ITERMX0
  150. *      write(*,*) 'SESPAC.eso : ITERATION GLOBALE #',IB100,'------------'
  151.  
  152.  
  153. *------ MISE A JOUR DE IPLVA1 -----------------------------------
  154.          CALL DTLREE ( ILVA1R )
  155.          CALL DTLREE ( ILVA1I )
  156.          IF ( IERR .NE. 0 ) RETURN
  157.          ILVA1R = ILVA2R
  158.          ILVA1I = ILVA2I
  159.  
  160. *------ MISE A JOUR DE IPLCH1 -----------------------------------
  161. *        On commence par le detruire'
  162.          SEGACT,IPLCH1
  163.          DO JVEC = 1, ILDIM
  164.            IPCHPO = IPLCH1.ICHPOI(JVEC)
  165.             CALL DTCHPO(IPCHPO)
  166.             IF ( IERR .NE. 0 ) RETURN
  167.          ENDDO
  168.          SEGDES,IPLCH1
  169.          CALL DTLCHP(IPLCH1)
  170.          IF ( IERR .NE. 0 ) RETURN
  171.  
  172. *        et on recopie IPLCH2 dans IPLCH1 '
  173.          SEGACT ,IPLCH2
  174.          N1 = 0
  175.          SEGINI ,IPLCH1
  176.          DO JVEC = 1, ILDIM
  177.             IPCHPO = IPLCH2.ICHPOI(JVEC)
  178.             CALL COPIE2 ( IPCHPO, IPCHP1 )
  179.             IF ( IERR .NE. 0 ) RETURN
  180.             IPLCH1.ICHPOI(**) = IPCHP1
  181.          ENDDO
  182.          SEGDES ,IPLCH1, IPLCH2
  183.  
  184.  
  185. *------- ITERMX1 ITERATIONs DU SOUS-ESPACE IPLVEC --------------------
  186.          DO 200 IB200=1,ITERMX1
  187. *            write(*,*) 'i0,i1=',IB100,IB200,' -> on va dans SESPC1'
  188.             CALL SESPC1 ( IPLCH2, IPKW2M, IPMASS )
  189.             IF ( IERR .NE. 0 ) RETURN
  190.  200     CONTINUE
  191.  
  192.  
  193. *------- PROJECTION DE K ET DE M SUR CE SOUS-ESPACE ------------------
  194. *         write(*,*) 'appel a SESPC2'
  195.          CALL SESPC2 ( IPLCH2, IPKW2M, IPARED )
  196.          CALL SESPC2 ( IPLCH2, IPMASS, IPBRED )
  197.          IF ( IERR .NE. 0 ) RETURN
  198.  
  199.  
  200. *         write(*,*) 'appel aux routines du QZ'
  201. *------- RESOLUTION DU PROBLEME REDUIT PAR LA METHODE DU QZ ---------
  202. * 1./ MISE SOUS FORME DE HESSENBERG SUPERIEURE DE A
  203. *     ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE B
  204.       CALCV = .TRUE.
  205.       CALL QZHESS(IPARED,IPBRED,CALCV,IPZ)
  206.       IF (IERR .NE. 0) RETURN
  207.  
  208. * 2./ MISE SOUS FORME QUASI-TRIANGULAIRE SUP D'UNE MATRICE DE HESSENBERG SUP A
  209. *     ET   TRIANGULARISATION SUPERIEURE SIMULTANEE de B
  210. *      CALL QZREDU(IPARED,IPBRED,0.0D0,CALCV,IPZ,NERR)
  211.       CALL QZREDU(IPARED,IPBRED,XZERO,CALCV,IPZ,NERR)
  212. *      XPREC2 = EPSLON(1.0D0)
  213. *      XPREC1 = XPREC2**0.5
  214. *      write(*,*) 'XPREC=',XPREC2,XPREC1
  215. *      CALL QZREDU(IPARED,IPBRED,XPREC1,CALCV,IPZ,NERR)
  216.       IF (IERR .NE. 0) RETURN
  217.       IF (NERR .NE. 0) WRITE (*,*) 'Mode ',NERR,' : trop d iterations !'
  218.  
  219. * 3./ FIN DE MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE DE A
  220. *     ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE B
  221. *     et EXTRACTION DES VALEURS PROPRES
  222.       CALL QZVALP(IPARED,IPBRED,ILAMBR,ILAMBI,ILAMBD,CALCV,IPZ)
  223.       IF (IERR.NE.0) RETURN
  224.  
  225. * 4./ CALCUL DES VECTEURS PROPRES COMPLEXES
  226.       CALL QZVECP(IPARED,IPBRED,ILAMBR,ILAMBI,ILAMBD,IPZ)
  227.       IF (IERR.NE.0) RETURN
  228.  
  229. *         write(*,*) 'appel a SESPC4'
  230. *------- CALCUL D'UNE APPROX. DES VECTEURS PROPRES
  231. *        A PARTIR DE LEURS PROJECTIONS SUR LE SOUS-ESPACE ---------
  232.          CALL SESPC4 ( IPLCH2, IPZ, ILAMBI )
  233.          IF ( IERR .NE. 0 ) RETURN
  234.  
  235.  
  236. *------- ECRITURE DES VALEURS PROPRES CORRESPONDANTES -------------
  237. c          segact,ILAMBR,ILAMBI,ILAMBD
  238.          JG = ILDIM
  239.          segini,ILVA2R,ILVA2I
  240.          DO JVEC = 1,ILDIM
  241.           XLAMBR = ILAMBR.PROG(JVEC)
  242.           XLAMBI = ILAMBI.PROG(JVEC)
  243.           XLAMBD = ILAMBD.PROG(JVEC)
  244.           ILVA2R.PROG(JVEC)= (XLAMBR) / (XLAMBD)
  245. *         la precision sur la determination des vp étant de 1.0D-6 on met a 0 si <1.D-5
  246. *          if(XLAMBI .lt. (1.D-5*XLAMBD))    XLAMBI = 0.D0
  247.           ILVA2I.PROG(JVEC)= (XLAMBI) / (XLAMBD)
  248. *          write(*,*) 'Lambda_',JVEC,'=',(ILVA2R.PROG(JVEC)),' + i',
  249. *     1    ILVA2I.PROG(JVEC)
  250.          ENDDO
  251.          segsup,ILAMBR,ILAMBI,ILAMBD
  252.  
  253.  
  254. *------- ON TESTE LA CONVERGENCE -----------------------------------
  255. *         write(*,*) 'appel a SESPC5'
  256.          NBMOD2 = NBMOD
  257.          CALL SESPC5 ( ILVA1R,ILVA1I, ILVA2R,ILVA2I, IPLCH2,
  258.      1                 IPKW2M,IPMASS, CVGVAL,CVGCHP, NBMOD2 )
  259. *         write(*,*) 'CVGVAL,CVGCHP=',CVGVAL,CVGCHP
  260.          IF ( IERR .NE. 0 ) RETURN
  261.  
  262. *        -- SI ON A CONVERGE, C'EST FINI !   --
  263. *         IF ( CVGVAL .and. CVGCHP ) THEN
  264.          IF ( CVGCHP ) THEN
  265. *            IF ( IIMPI .EQ. 2 ) THEN
  266.                WRITE ( IOIMP, 1000 ) IB100
  267.  1000          FORMAT( /1X, 'On a effectue ',I2,' iterations.', /)
  268. *            ENDIF
  269.             GOTO 110
  270.          ENDIF
  271.  
  272. *        -- SI NON, PAS DE CONVERGE, MAIS ON RENVOIE LA SOLUTION !
  273.          IF ( IB100 .EQ. ITERMX0 ) THEN
  274.             WRITE ( IOIMP, 2000 ) ITERMX0
  275.  2000       FORMAT( /1X, 'Pas de convergence apres ',I2,' iterations.',
  276.      1              /1X, 'La solution est quand meme renvoyee.',
  277.      2              /1X, 'L''execution continue ...', / )
  278.  
  279.          ENDIF
  280.  100  CONTINUE
  281.  110  CONTINUE
  282.  
  283.  
  284.  
  285. ** estimation d'une borne superieure de l'erreur sur les valeurs propres
  286. ** (c.f. Argyris, Wilkinson)
  287. *... a revoir ....
  288.  
  289. ****************************************************************
  290. *                  NETTOYAGE DE LA MEMOIRE                     *
  291. ****************************************************************
  292.  
  293. ******
  294. *     -- ON DETRUIT IPLVA1 ET IPLCH1
  295. ***
  296.       CALL DTLREE ( ILVA1R )
  297.       CALL DTLREE ( ILVA1I )
  298.       IF ( IERR .NE. 0 ) RETURN
  299.  
  300.       SEGACT ,IPLCH1
  301.       DO 400 IB400 = 1, ILDIM
  302.          IPCHPO = IPLCH1.ICHPOI(IB400)
  303.          CALL DTCHPO (IPCHPO)
  304.          IF ( IERR .NE. 0 ) RETURN
  305.  400  CONTINUE
  306.       SEGDES ,IPLCH1
  307.       CALL DTLCHP ( IPLCH1 )
  308.       IF ( IERR .NE. 0 ) RETURN
  309.  
  310. ******
  311. *     -- ON RENVOIE LES VALEURS ET VECTEURS PROPRES --
  312. ***
  313.       IPLVAR = ILVA2R
  314.       IPLVAI = ILVA2I
  315.       IPLVEC = IPLCH2
  316. *     on modifie le nbmod proposé si le dernier mode est complexe
  317. *     et que l'on a besoin du nbmod+1^ieme vecteur
  318. *      write(*,*) 'sespac: NBMOD,NBMOD2=',NBMOD,NBMOD2
  319.       NBMOD = NBMOD2
  320.  
  321.       RETURN
  322.       END
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  

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