C SESPAC SOURCE FANDEUR 13/01/29 21:16:26 7683 C SUBROUTINE SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS ) ************************************************************************ * * SESPAC * ----------- * * FONCTION: * --------- * * CONSTRUIT NBMOD ELEMENTS PROPRES EN ITERANT LE SOUS-ESPACE * IPLVEC, JUSQU'A LA CONVERGENCE DE CELUI-CI. OBTENTION DE VALEURS PROPRES * RELLES OU COMPLEXES. * * MODE D'APPEL: * ------------- * * CALL SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS ) * * PARAMETRES: (E)=ENTREE (S)=SORTIE * ----------- * * IPLVAL ENTIER (S) POINTEUR DE L'OBJET 'LISTREEL' CONTENANT * LA SUITE DE VALEURS PROPRES OBTENUES. * IPLVEC ENTIER (E) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT * LE SOUS-ESPACE INITIAL. * IPLVEC ENTIER (S) POINTEUR DE L'OBJET 'LISTCHPO' CONTENANT * LE SOUS-ESPACE FINAL. EN ORTHONORMALISANT * LES 'CHPOINT' DE CET ESPACE ON OBTIENT LES * VECTEURS PROPRES RECHERCHES. * NBMOD ENTIER (E) NOMBRE DE VECTEURS RECHERCHES. IPLVEC * CONTIENT PLUS QUE NBMOD 'CHPO', CAR CECI * PERMET DE CONVERGER PLUS RAPIDEMENT. * * IPKW2M ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' K * IPMASS ENTIER (E) POINTEUR SUR L'OBJET 'RIGIDITE' M * * * AUTEUR, DATE DE CREATION: * ------------------------- * * C. LE BIDEAU 09 / 2001 ( FORTRAN + ESOPE ) * MODIF Benoit Prabel Mars 2009 * *********************************************************************** * SUBROUTINE SESPAC ( IPLVAL, IPLVEC, NBMOD, IPKW2M, IPMASS ) SUBROUTINE SESPAC(IPLVAR,IPLVAI,IPLVEC,NBMOD,IPKW2M,IPMASS) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) DATA XZERO / 0.0D0 / -INC CCOPTIO -INC SMLCHPO -INC SMLREEL -INC SMCHPOI -INC SMRIGID ****** * -- CONSTANTES -- *** * nbre d'iteration de la boucle complete et de celle sur pb global (=puissance) PARAMETER ( ITERMX0 = 10 ) PARAMETER ( ITERMX1 = 3 ) * PARAMETER ( ITERMX0 = 3 ) * PARAMETER ( ITERMX1 = 2 ) * PARAMETER ( ITERMX1 = 1 ) ****** * -- ARGUMENTS -- *** POINTEUR IPLVEC.MLCHPO * IPLVER.MLCHPO,IPLVEI.MLCHPO POINTEUR IPLVAR.MLREEL,IPLVAI.MLREEL INTEGER IPKW2M, IPMASS, NBMOD ****** * -- VARIABLES LOCALES -- *** POINTEUR IPLCH1.MLCHPO,IPLCH2.MLCHPO POINTEUR ILVA1R.MLREEL,ILVA1I.MLREEL,ILVA2R.MLREEL,ILVA2I.MLREEL POINTEUR ILAMBR.XMATRI,ILAMBI.XMATRI,ILAMBD.XMATRI,IPZ.XMATRI INTEGER ILDIM,NERR LOGICAL CALCV,CVGVAL,CVGCHP **************************************************************** * INITIALISATIONS * **************************************************************** * ouverture des vecteurs itérés SEGACT IPLVEC*MOD ILDIM = IPLVEC.ICHPOI( /1 ) * recherche des listmots inconnue (IPLMOX) et duale (IPLMOY) IPCHPO = IPLVEC.ICHPOI(1) CALL MUCPRI (IPCHPO, IPKW2M, IPCX) CALL CORRSP (IPKW2M, IPCHPO, IPCX, IPLMOX, IPLMOY) CALL DTCHPO (IPCX) * CALL SESBAS (IPLVEC, IPKW2M, IPMASS, IPLMOX, IPLMOY) * creation des valeurs propres (initialisées à 0) JG = ILDIM SEGINI,ILVA1R,ILVA1I,ILVA2R,ILVA2I DO JVEC = 1,ILDIM * ILVA1.PROG(JVEC) = 0.D0 ILVA2R.PROG(JVEC) = 0.D0 ILVA2I.PROG(JVEC) = 0.D0 ENDDO * creation des vecteurs propres (initialisées à IPLVEC) N1 = ILDIM SEGINI,IPLCH1 DO JVEC = 1, ILDIM IPCHPO = IPLVEC.ICHPOI(JVEC) CALL COPIE2 ( IPCHPO, IPCH1 ) IF ( IERR .NE. 0 ) RETURN IPLCH1.ICHPOI(JVEC) = IPCH1 ENDDO * impressions * do JVEC=1,ILDIM * MCHPOI = IPLVEC.ICHPOI(JVEC) * segact,MCHPOI * MSOUPO=IPCHP(1) * segact,MSOUPO * MPOVAL=IPOVAL * segact,MPOVAL * write(*,*) ' X^0_',JVEC,'=',VPOCHA(1,1),VPOCHA(1,2),VPOCHA(1,3) * segdes,MPOVAL,MSOUPO,MCHPOI * enddo SEGDES,IPLVEC,IPLCH1 IPLCH2 = IPLVEC IPLVEC = 0 **************************************************************** * ITERMX0 ITERATIONS associant iterations GLOBALES et sur systeme reduit * * REPETER JUSQU'A: * CONVERGER * * DEPASSER ITERMX0 ITERATIONS **************************************************************** DO 100 IB100 = 1,ITERMX0 * write(*,*) 'SESPAC.eso : ITERATION GLOBALE #',IB100,'------------' *------ MISE A JOUR DE IPLVA1 ----------------------------------- CALL DTLREE ( ILVA1R ) CALL DTLREE ( ILVA1I ) IF ( IERR .NE. 0 ) RETURN ILVA1R = ILVA2R ILVA1I = ILVA2I *------ MISE A JOUR DE IPLCH1 ----------------------------------- * On commence par le detruire' SEGACT,IPLCH1 DO JVEC = 1, ILDIM IPCHPO = IPLCH1.ICHPOI(JVEC) CALL DTCHPO(IPCHPO) IF ( IERR .NE. 0 ) RETURN ENDDO SEGDES,IPLCH1 CALL DTLCHP(IPLCH1) IF ( IERR .NE. 0 ) RETURN * et on recopie IPLCH2 dans IPLCH1 ' SEGACT ,IPLCH2 N1 = 0 SEGINI ,IPLCH1 DO JVEC = 1, ILDIM IPCHPO = IPLCH2.ICHPOI(JVEC) CALL COPIE2 ( IPCHPO, IPCHP1 ) IF ( IERR .NE. 0 ) RETURN IPLCH1.ICHPOI(**) = IPCHP1 ENDDO SEGDES ,IPLCH1, IPLCH2 *------- ITERMX1 ITERATIONs DU SOUS-ESPACE IPLVEC -------------------- DO 200 IB200=1,ITERMX1 * write(*,*) 'i0,i1=',IB100,IB200,' -> on va dans SESPC1' CALL SESPC1 ( IPLCH2, IPKW2M, IPMASS ) IF ( IERR .NE. 0 ) RETURN 200 CONTINUE *------- PROJECTION DE K ET DE M SUR CE SOUS-ESPACE ------------------ * write(*,*) 'appel a SESPC2' CALL SESPC2 ( IPLCH2, IPKW2M, IPARED ) CALL SESPC2 ( IPLCH2, IPMASS, IPBRED ) IF ( IERR .NE. 0 ) RETURN * write(*,*) 'appel aux routines du QZ' *------- RESOLUTION DU PROBLEME REDUIT PAR LA METHODE DU QZ --------- * 1./ MISE SOUS FORME DE HESSENBERG SUPERIEURE DE A * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE B CALCV = .TRUE. CALL QZHESS(IPARED,IPBRED,CALCV,IPZ) IF (IERR .NE. 0) RETURN * 2./ MISE SOUS FORME QUASI-TRIANGULAIRE SUP D'UNE MATRICE DE HESSENBERG SUP A * ET TRIANGULARISATION SUPERIEURE SIMULTANEE de B * CALL QZREDU(IPARED,IPBRED,0.0D0,CALCV,IPZ,NERR) CALL QZREDU(IPARED,IPBRED,XZERO,CALCV,IPZ,NERR) * XPREC2 = EPSLON(1.0D0) * XPREC1 = XPREC2**0.5 * write(*,*) 'XPREC=',XPREC2,XPREC1 * CALL QZREDU(IPARED,IPBRED,XPREC1,CALCV,IPZ,NERR) IF (IERR .NE. 0) RETURN IF (NERR .NE. 0) WRITE (*,*) 'Mode ',NERR,' : trop d iterations !' * 3./ FIN DE MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE DE A * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE B * et EXTRACTION DES VALEURS PROPRES CALL QZVALP(IPARED,IPBRED,ILAMBR,ILAMBI,ILAMBD,CALCV,IPZ) IF (IERR.NE.0) RETURN * 4./ CALCUL DES VECTEURS PROPRES COMPLEXES CALL QZVECP(IPARED,IPBRED,ILAMBR,ILAMBI,ILAMBD,IPZ) IF (IERR.NE.0) RETURN * write(*,*) 'appel a SESPC4' *------- CALCUL D'UNE APPROX. DES VECTEURS PROPRES * A PARTIR DE LEURS PROJECTIONS SUR LE SOUS-ESPACE --------- CALL SESPC4 ( IPLCH2, IPZ, ILAMBI ) IF ( IERR .NE. 0 ) RETURN *------- ECRITURE DES VALEURS PROPRES CORRESPONDANTES ------------- segact,ILAMBR,ILAMBI,ILAMBD JG = ILDIM segini,ILVA2R,ILVA2I DO JVEC = 1,ILDIM XLAMBR = ILAMBR.RE(1,JVEC) XLAMBI = ILAMBI.RE(1,JVEC) XLAMBD = ILAMBD.RE(1,JVEC) ILVA2R.PROG(JVEC)= (XLAMBR) / (XLAMBD) * la precision sur la determination des vp étant de 1.0D-6 on met a 0 si <1.D-5 * if(XLAMBI .lt. (1.D-5*XLAMBD)) XLAMBI = 0.D0 ILVA2I.PROG(JVEC)= (XLAMBI) / (XLAMBD) * write(*,*) 'Lambda_',JVEC,'=',(ILVA2R.PROG(JVEC)),' + i', * 1 ILVA2I.PROG(JVEC) ENDDO segsup,ILAMBR,ILAMBI,ILAMBD *------- ON TESTE LA CONVERGENCE ----------------------------------- * write(*,*) 'appel a SESPC5' NBMOD2 = NBMOD CALL SESPC5 ( ILVA1R,ILVA1I, ILVA2R,ILVA2I, IPLCH2, 1 IPKW2M,IPMASS, CVGVAL,CVGCHP, NBMOD2 ) * write(*,*) 'CVGVAL,CVGCHP=',CVGVAL,CVGCHP IF ( IERR .NE. 0 ) RETURN * -- SI ON A CONVERGE, C'EST FINI ! -- * IF ( CVGVAL .and. CVGCHP ) THEN IF ( CVGCHP ) THEN * IF ( IIMPI .EQ. 2 ) THEN WRITE ( IOIMP, 1000 ) IB100 1000 FORMAT( /1X, 'On a effectue ',I2,' iterations.', /) * ENDIF GOTO 110 ENDIF * -- SI NON, PAS DE CONVERGE, MAIS ON RENVOIE LA SOLUTION ! IF ( IB100 .EQ. ITERMX0 ) THEN WRITE ( IOIMP, 2000 ) ITERMX0 2000 FORMAT( /1X, 'Pas de convergence apres ',I2,' iterations.', 1 /1X, 'La solution est quand meme renvoyee.', 2 /1X, 'L''execution continue ...', / ) ENDIF 100 CONTINUE 110 CONTINUE ** estimation d'une borne superieure de l'erreur sur les valeurs propres ** (c.f. Argyris, Wilkinson) *... a revoir .... **************************************************************** * NETTOYAGE DE LA MEMOIRE * **************************************************************** ****** * -- ON DETRUIT IPLVA1 ET IPLCH1 *** CALL DTLREE ( ILVA1R ) CALL DTLREE ( ILVA1I ) IF ( IERR .NE. 0 ) RETURN SEGACT ,IPLCH1 DO 400 IB400 = 1, ILDIM IPCHPO = IPLCH1.ICHPOI(IB400) CALL DTCHPO (IPCHPO) IF ( IERR .NE. 0 ) RETURN 400 CONTINUE SEGDES ,IPLCH1 CALL DTLCHP ( IPLCH1 ) IF ( IERR .NE. 0 ) RETURN ****** * -- ON RENVOIE LES VALEURS ET VECTEURS PROPRES -- *** IPLVAR = ILVA2R IPLVAI = ILVA2I IPLVEC = IPLCH2 * on modifie le nbmod proposé si le dernier mode est complexe * et que l'on a besoin du nbmod+1^ieme vecteur * write(*,*) 'sespac: NBMOD,NBMOD2=',NBMOD,NBMOD2 NBMOD = NBMOD2 RETURN END