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