1 : C PROCH3 SOURCE BP208322 10/06/29 21:15:01 6698 2 : C SUBROUTINE PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE, 3 : C INSYM, MTAB3, I ) 4 : ************************************************************************ 5 : * 6 : * P R O C H 3 7 : * ----------- 8 : * 9 : * FONCTION: 10 : * --------- 11 : * 12 : * RECHERCHE DES NBMOD MODES PROPRES LES PLUS PROCHES DE FREQ 13 : * ET CONSTRUIT L'OBJET 'SOLUTION' CORRESPONDANT. 14 : * 15 : * MODE D'APPEL: 16 : * ------------- 17 : * 17 : * CALL PROCH3( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE, INSYM ) 18 : * 19 : * PARAMETRES: (E)=ENTREE (S)=SORTIE 20 : * ----------- 21 : * 22 : * FREQ REEL DP (E) FREQUENCE A APPROCHER PAR UNE FREQUENCE 23 : * PROPRE. 24 : * NBMOD ENTIER (E) ORDRE DE MULTIPLICITE DE LA FREQ PROPRE 25 : * IPRIGI ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT 26 : * LA MATRICE DE RIGIDITE. 27 : * IPMASS ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' CONTENANT 28 : * LA MATRICE MASSE. 29 : * INF0 ENTIER (E) NOMBRE DE TERMES DIAGONAUX NEGATIFS DE LA 30 : * 'RIGIDITE' "K" DECOMPOSEE EN L.D.LT . 31 : * CE NOMBRE N'EST PAS NUL A CAUSE DE LA FACON 32 : * DONT SONT INTRODUITS LES BLOCAGES DES 33 : * D.D.L. (MULTIPLICATEURS DE LAGRANGE "LX"). 34 : * IPMODE ENTIER (S) POINTEUR DE L'OBJET 'SOLUTION' CONTENANT 35 : * LE MODE PROPRE TROUVE. 36 : * LIMAGE BOOLEEN (E) VRAI SI ON SOUHAITE DE FREQ. NEGATIVES. 37 : * 38 : * LEXIQUE: (ORDRE ALPHABETIQUE) 39 : * -------- 40 : * 41 : * IPKW2M ENTIER POINTEUR DE LA 'RIGIDITE' "DECALEE" K - W2.M 42 : * IPLVER LISTCHPO LISTE DE CHPO A ITERER. ILS FORMENT UN SOUS- 43 : * ESPACE QUI CONVERGE VERS LE SOUS-ESPACE PROPRE 44 : * DOMINANT DE DIMENSION NBMOD. SORTIE: PARTIE REELLE 45 : * DES CHPOINTS SOLUTIONS 46 : * IPLVEI LISTCHPO LISTE DE CHPO.SORTIE: PARTIE IMAGINAIRE 47 : * DES CHPOINTS SOLUTIONS 48 : * IPVECP ENTIER POINTEUR DU 'CHPOINT' QUI CONTIENT DES NOMBRES 49 : * ALEATOIRES, PUIS UN VECTEUR PROPRE. 50 : * OMEGA2 REEL DP PULSATION PROPRE TROUVEE AU CARRE. 51 : * ( NE SERT A RIEN ! ) 52 : * PROPRE REEL DP TABLEAU CONTENANT DES CARACTERISTIQUES DU 53 : * MODE PROPRE CALCULE: 54 : * PROPRE(1) = FREQ. PROPRE , 55 : * PROPRE(2) = MASSE GEN. 56 : * PROPRE(3)ET(4) ET(5) DEPL.GEN. SELON X,Y,Z 57 : * W2 REEL DP PULSATION AU CARRE A APPROCHER. 58 : * 59 : * MODE DE FONCTIONNEMENT: 60 : * ----------------------- 61 : * 62 : * LE CALCUL D'UN VECTEUR PROPRE SE FAIT PAR LA METHODE DE 63 : * L'ITERATION D'UN SOUS-ESPACE, AVEC DECALAGE INITIAL ("SHIFTING") 64 : * DE W2 65 : * 66 : * AUTEUR, DATE DE CREATION: 67 : * ------------------------- 68 : * 69 : * PASCAL MANIGOT 16 OCTOBRE 1984 70 : * 71 : * MODIF: Benoit PRABEL MARS 2009 72 : * 73 : * LANGAGE: 74 : * -------- 75 : * 76 : * ESOPE 77 : * 78 : ************************************************************************ 79 : 80 : SUBROUTINE PROCH3 ( FREQ,NBMOD,IPRIGI,IPMASS,INF0,IPMODE,LIMAGE 81 : & ,INSYM, MTAB3, I) 82 : 83 : IMPLICIT INTEGER(I-N) 84 : IMPLICIT REAL*8 (A-H,O-Z) 85 : 86 : -INC CCREEL 87 : -INC CCOPTIO 88 : -INC SMLCHPO 89 : -INC SMLREEL 90 : -INC SMCHPOI 91 : 92 : ****** 93 : * -- CONSTANTES -- 94 : *** 95 : PARAMETER (LPROPR = 10) 96 : PARAMETER (DEUXPI = (2.D0*XPI)) 97 : 98 : ****** 99 : * -- ARGUMENTS -- 100 : *** 101 : REAL*8 FREQ 102 : INTEGER NBMOD, IPRIGI, IPMASS, INF0, IPMODE 103 : LOGICAL LIMAGE 104 : CHARACTER*4 MOTCLE 105 : 106 : ****** 107 : * -- VARIABLES LOCALES -- 108 : *** 109 : POINTEUR IPLVER.MLCHPO, IPLVE1.MLCHPO, IPLVEI.MLCHPO 110 : POINTEUR IPLVAR.MLREEL, IPLVA1.MLREEL, IPLVAI.MLREEL 111 : INTEGER NBVEC 112 : REAL*8 OMEGA2, PROPRE(LPROPR) 113 : 114 : ****************************************************************** 115 : * INITIALISATION * 116 : ****************************************************************** 117 : 118 : ****** 119 : * -- CREATION DE (K-W2M) -- 120 : *** 121 : W2 = (FREQ * DEUXPI) ** 2 122 : * WRITE(6,*)'DECALAGE ', W2 123 : IF( LIMAGE ) THEN 124 : W2 = SIGN(W2,FREQ) 125 : ENDIF 125 : CALL DECALE ( IPRIGI, IPMASS, W2, IPKW2M ) 126 : IF ( IERR .NE. 0 ) RETURN 127 : 128 : ****** 129 : * -- INITIALISATION DES ITERATIONS: CREATION D'UNE LISTE 130 : * DE NBVEC 'CHPOINT' ALEATOIRE -- 131 : *** 132 : 133 : * bp (03.2009) : si on souhaite eventuellement fournir 1 CHPOINT 133 : call LIROBJ('LISTCHPO',MLCHP1,0,IRET1) 134 : if(IRET1 .ne. 0) then 135 : segact,MLCHP1 136 : NCHPO1 = MLCHP1.ICHPOI(/1) 137 : IPVECP = 0 137 : * call COPIE2(MLCHP1.ICHPOI(1),IPVECP) 138 : ICHP1 = MLCHP1.ICHPOI(1) 138 : call COPIE2(ICHP1,IPVECP) 139 : * A Faire : vérifier l'accord entre chpoint et matrices...? 140 : else 140 : CALL ALEAT1 (IPKW2M, IPVECP ) 141 : endif 142 : IF ( IERR .NE. 0 ) RETURN 143 : N1 = 1 144 : SEGINI ,IPLVER 145 : IPLVER.ICHPOI(1) = IPVECP 146 : ** on cherce le nombre de ddl sans relations 147 : nddl = 0 148 : ndlx = 0 149 : mchpoi = ipvecp 150 : segact mchpoi 151 : nsous = ipchp(/1) 152 : do 81 isous = 1, nsous 153 : msoupo = ipchp(isous) 154 : segact msoupo 155 : mpoval = ipoval 156 : segact mpoval 157 : nn = vpocha(/1) 158 : nc1 = vpocha(/2) 159 : nddl = nddl + nn*nc1 160 : do 95 inc = 1,nc1 161 : if (nocomp(inc).eq. 'LX ') ndlx = ndlx + nn 162 : 95 continue 163 : segdes msoupo, mpoval 164 : 81 continue 165 : segdes mchpoi 166 : 167 : nddl = nddl - ndlx - (ndlx/2) 168 : nbvec = max(2*nbmod, 8) 169 : if (nbvec .gt. nddl) nbvec = nddl 170 : 171 : * on continue le tirage au sort (nous aurions pu passer par Lanzos!) 172 : DO 100 IB100 = 1, NBVEC-1 173 : if ( (IRET1 .ne. 0) .and. (IB100 .lt. NCHPO1) ) then 174 : IPVECP = 0 174 : * call COPIE2(MLCHP1.ICHPOI(1 + IB100),IPVECP) 175 : ICHP1 = MLCHP1.ICHPOI(1 + IB100) 175 : call COPIE2(ICHP1,IPVECP) 176 : else 176 : CALL ALEAT1 (IPKW2M, IPVECP ) 177 : endif 178 : IF ( IERR .NE. 0 ) RETURN 179 : IPLVER.ICHPOI(**) = IPVECP 180 : 100 CONTINUE 181 : SEGDES ,IPLVER 182 : if(IRET1 .ne. 0) SEGDES,MLCHP1 183 : 184 : * cette orthonormalisation n'est pas indispensable, dans la plupart 185 : * des cas, 186 : * mais, on l'a gardée pour tester la singularité de la masse 187 : * et mettre à jour nbvec si necessaire. On n'a besoin de ce test 188 : * qu'au debut, car si l'on a eliminé les eventuels vecteurs du noyau 189 : * de la masse ils ne peuvent pas ressurgir lors du processus iteratif 190 : * car ils correspondent à des frequences infinies 191 : * IF ( INSYM .EQ. 0) THEN 191 : CALL GRAAMO ( IPLVER, IPMASS ) 192 : IF ( IERR .NE. 0 ) RETURN 193 : segact IPLVER 194 : nbvec = IPLVER.ICHPOI(/1) 195 : segdes iplver 196 : if (nbmod .gt. nbvec) nbmod = nbvec 197 : * END IF 198 : 199 : ****************************************************************** 200 : * RESOLUTION * 201 : ****************************************************************** 202 : 203 : ****** 204 : * -- RESOLUTION PAR ITERATION DU SOUS-ESPACE IPLVER 205 : * ET NORMALISATION DES MODES PROPRES -- 206 : *** 207 : *on differentie NBMOD=dim du sous espace initialement recherché et 208 : *NBMOD2 qui peut etre modifié par sespac si dernier vecteur complexe 209 : NBMOD2 = NBMOD 210 : IF (INSYM .EQ. 0) THEN 210 : CALL SESPA ( IPLVAR, IPLVER, NBMOD, IPKW2M, IPMASS ) 211 : IF ( IERR .NE. 0 ) RETURN 212 : segact iplver 213 : nblgn = IPLVER.ICHPOI(/1) 214 : ICOMP = 0 215 : mlchpo = iplver 216 : segact mlchpo*mod 217 : do 101 ibmod = 1,nbmod 218 : ix = ichpoi(ibmod) 218 : CALL MOTS1( IPLMOT, MOTCLE ) 218 : call norma1 (ix,iplmot,motcle,ixx) 218 : call dtchpo(ix) 219 : ichpoi(ibmod) = ixx 220 : 101 continue 221 : segdes mlchpo 222 : 223 : ELSE 223 : CALL SESPAC(IPLVAR,IPLVAI,IPLVER,NBMOD2,IPKW2M,IPMASS) 224 : ICOMP = 0 225 : * 226 : END IF 227 : 228 : 229 : ****************************************************************** 230 : * CREATION DE L'OBJET REPRESENTANT LE MODE MULTIPLE * 231 : ****************************************************************** 232 : 233 : ****** 234 : * -- TRANSLATION DES FREQUENCES PROPRES -- CLASSEMENT 235 : * DES MODES PROPRES -- 236 : *** 237 : SEGACT ,IPLVAR*MOD 238 : * IF (ICOMP .EQ. 0) THEN 239 : *---- Cas d'un appel a sespa (matrices symetriques => mode Reel) 240 : IF (INSYM .EQ. 0) THEN 241 : DO 200 IB200 =1, NBMOD 242 : FREQPR = IPLVAR.PROG(IB200) 242 : CALL W2FREQ ( FREQPR, W2, OMEGA2, FREQP1, LIMAGE ) 243 : IF ( IERR .NE. 0 ) RETURN 244 : IPLVAR.PROG(IB200) = FREQP1 245 : 200 CONTINUE 246 : 247 : *---- Cas d'un appel a sespaC (matrices nonsymetriques => mode Reel ou Complexe) 248 : ELSE 249 : SEGACT IPLVAI*MOD 250 : * DO 400 IB400 = 1, NBMOD2 251 : DO 400 IB400 = 1, NBMOD2 252 : FREQPR = IPLVAR.PROG(IB400) 253 : FREQPI = IPLVAI.PROG(IB400) 253 : CALL W2FRQC ( FREQPR, FREQPI, W2,XRW2,XIW2, FREQP1, FREQP2) 254 : IF ( IERR .NE. 0 ) RETURN 255 : * IPLVAR.PROG(IB400) = FREQP1 256 : * IPLVAI.PROG(IB400) = FREQP2 257 : * on se débarasse des erreurs d'arrondis de W2FRQC 258 : if (FREQPI .eq. 0.) then 259 : if ( (FREQPR + W2) .lt. 0.) then 260 : IPLVAR.PROG(IB400)= 0.D0 261 : IPLVAI.PROG(IB400)= FREQP2 262 : else 263 : IPLVAR.PROG(IB400)= FREQP1 264 : IPLVAI.PROG(IB400)= 0.D0 265 : endif 266 : else 267 : IPLVAR.PROG(IB400)= FREQP1 268 : IPLVAI.PROG(IB400)= FREQP2 269 : endif 270 : * WRITE(6,*) IB400,' FREQUENCE =',FREQP1,'+i',FREQP2, 271 : * & ' soit',(IPLVAR.PROG(IB400)),'+i',(IPLVAI.PROG(IB400)) 272 : 400 CONTINUE 273 : END IF 274 : 275 : ****** 276 : * -- ON GARDE LES NBMOD2 PREMIERS ELEMENTS PROPRES -- 277 : * -- ET ON DETRUIT LES AUTRES -- 278 : *** 279 : 280 : * -- LES VALEURS PROPRES -- 281 : 282 : * JG = NBMOD 283 : JG = NBMOD2 284 : SEGADJ IPLVAR 285 : SEGDES IPLVAR 286 : 287 : IF (INSYM .NE. 0) THEN 288 : * JG = NBMOD 289 : JG = NBMOD2 290 : SEGADJ IPLVAI 291 : SEGDES IPLVAI 292 : END IF 293 : 294 : * -- LES MODES PROPRES -- 295 : SEGACT IPLVER 296 : * si le dernier mode est complexe, alors il nous faut le (nbmod + 1)^ieme 297 : * qui correspond a la partie imaginaire indissociable de la partie reelle, 298 : * ce qui justifie l'introduction de nbmod2 299 : DO 550 IB550 = (NBMOD2 + 1), NBVEC 300 : IPCHPO = IPLVER.ICHPOI(IB550) 300 : CALL DTCHPO ( IPCHPO ) 301 : 550 CONTINUE 302 : N1 = NBMOD2 303 : SEGADJ IPLVER 304 : SEGDES IPLVER 305 : 306 : 307 : ****** 308 : * -- CLASSEMENT DES MODES PROPRES -- 309 : *** 310 : IF (INSYM .eq. 0) THEN 310 : CALL ORDVEC ( IPLVAR, IPLVER ) 311 : END IF 312 : * rem: mode complexes deja classés par module de la val propre croissant 313 : 314 : 315 : ****** 316 : * -- AFFICHAGE -- 317 : *** 318 : IF ( IIMPI .EQ. 2 ) THEN 319 : WRITE ( IOIMP, 1000 ) FREQ 320 : 1000 FORMAT( /8X, 'SOLUTION POUR LA FREQUENCE:', E12.6, 321 : 1 /8X, '---------------------------', / ) 322 : ENDIF 323 : 324 : ****** 325 : * -- CONSTRUCTION DE LA SOLUTION -- 326 : *** 327 : IF (INSYM .EQ. 0) THEN 327 : * CALL CRSOLU (FREQ,IPLVAR,IPLVER,NBMOD,IPKW2M,IPMASS 328 : * & ,IPMODE, INF0) 328 : CALL CRSOLU (W2,IPLVAR,IPLVER,NBMOD,IPKW2M,IPMASS 329 : & ,IPMODE, INF0) 330 : ELSE 331 : * WRITE(6,*)'APPEL CRSOL1 ICOMP= ',ICOMP 331 : * CALL CRSOL1(FREQ,IPLVAR,IPLVAI,IPLVER,IPLVEI,NBMOD, 332 : * & IPKW2M,IPMASS,MTAB3, ICOMP, I) 332 : CALL CRSOL1(W2,IPLVAR,IPLVAI,IPLVER,NBMOD,NBMOD2, 333 : & IPKW2M,IPMASS,MTAB3,I,INF0) 334 : END IF 335 : 336 : 337 : ****** 338 : * -- SUPPRESSION DES OBJETS DE TRAVAIL -- 339 : *** 339 : CALL DTRIGI ( IPKW2M ) 340 : IF ( IERR .NE. 0 ) RETURN 340 : CALL DTLCHP ( IPLVER ) 341 : IF ( IERR .NE. 0 ) RETURN 341 : CALL DTLREE ( IPLVAR ) 342 : IF ( IERR .NE. 0 ) RETURN 343 : IF (INSYM .NE. 0) THEN 343 : * CALL DTLCHP ( IPLVEI) 343 : CALL DTLREE ( IPLVAI ) 344 : END IF 345 : 346 : RETURN 347 : END 348 : 349 : 350 : 351 : 352 : 353 : 354 : 355 : 356 : 357 : 358 : 359 : 360 : 361 : 362 : 363 :
© Cast3M 2003 - All rights reserved.
Disclaimer