1 : C ITINVC SOURCE BP208322 09/03/20 21:15:10 6331
2 : SUBROUTINE ITINVC (IPA,IPB,IPRX,IPIX,PROPRE,CONVRG,ITERMX,IPMX)
3 : IMPLICIT INTEGER(I-N)
4 : IMPLICIT REAL*8 (A-H,O-Z)
5 : ************************************************************************
6 : *
7 : * I T I N V C
8 : * -----------
9 : *
10 : * FONCTION:
11 : * ---------
12 : *
13 : * RESOUDRE, PAR ITERATIONS INVERSES, UN SYSTEME D'EQUATIONS:
14 : * |A|.(X) = V.|B|.(X)
15 : * |A| ET |B| ETANT 2 'RIGIDITE',
16 : * (X) UN 'CHPOINT' A DETERMINER ET
17 : * V UN 'FLOTTANT' EGALEMENT A DETERMINER.
18 : *
19 : * ("ITINVC" VAUT POUR IT-ERATIONS INV-ERSES C-OMPLEXES)
20 : *
21 : * MODE D'APPEL:
22 : * -------------
23 : *
23 : * CALL ITINVC (IPA,IPB,IPRX,IPIX,PROPRE,CONVRG,ITERMX,IPMX)
24 : *
25 : * PARAMETRES: (E)=ENTREE (S)=SORTIE
26 : * -----------
27 : *
28 : * IPA ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' |A|.
29 : * IPB ENTIER (E) POINTEUR DE L'OBJET 'RIGIDITE' |B|.
30 : * IPRX ENTIER (E) POINTEUR DE L'OBJET 'CHPOINT' DE DEPART.
31 : * (S) PARTIE REELLE DEL'OBJET 'CHPOINT' SOLUTION.
32 : * IPIX ENTIER (S) PARTIE IMAGINAIRE DE L'OBJET SOLUTION
33 : * PROPRE REEL DP (S) TABLEAU CONTENANT DES CARACTERISTIQUES DU
34 : * MODE PROPRE CALCULE. ACTUELLEMENT,
35 : * PROPRE(1) = PARTIE REELLE DU MODE PRORPE,
36 : * PROPRE(2) = (X)T.|B|.(X) , (X) 'CHPOINT'
37 : * SOLUTION,
38 : * PROPRE(3,4,5) DEPL.GEN. REELS SELON X,Y,Z
39 : * PROPRE(6)= PARTIE IMAGINAIRE DU MODE PROPRE
40 : * PROPRE(7)=PARTIE IM. DE XT.|B|.X
41 : * PROPRE(8,9,10) PARTIE IM DES DEP. GEN.
42 : * CONVRG LOGIQUE (S) INDIQUE PAR .TRUE. OU .FALSE. SI LA
43 : * CONVERGENCE A EU LIEU OU NON.
44 : * ITERMX ENTIER (E) NOMBRE MAXIMUM D'ITERATIONS PERMIS.
45 : *
46 : * LEXIQUE: (ORDRE ALPHABETIQUE)
47 : * --------
48 : *
49 : * DIFREL REEL SP VOIR LE S.P. "ITINV1".
50 : * IACCEL ENTIER NOMBRE D'ITERATIONS CONSECUTIVES EFFECTUEES
51 : * SANS ACCELERATION DE CONVERGENCE.
52 : * IPX0 ENTIER VOIR LE S.P. "ITINV1".
53 : * IPX1 ENTIER VOIR LE S.P. "ITINV1".
54 : * IPX2 ENTIER VOIR LE S.P. "ITINV1".
55 : * NBITER ENTIER NOMBRE D'ITERATIONS EFFECTUEES.
56 : * NUMXBX ENTIER NUMERO DE LA DERNIERE ITERATION OU L'ON A
57 : * CALCULE "XT.B.X" POUR 2 'CHPOINT' ITERES
58 : * CONSECUTIFS.
59 : * VALPP REEL DP VALEUR PROPRE REELLE ASSOCIEE AU 'CHPOINT' SOLUTION.
60 : *
61 : * MODE DE FONCTIONNEMENT:
62 : * -----------------------
63 : *
64 : * METHODE DES ITERATIONS INVERSES:
65 : *
66 : * LA SUITE "(X)I" TELLE QUE:
67 : * |A| . (X)I+1 = |B| . (X)I
68 : * TEND VERS LA (OU UNE DES) SOLUTION(S) DE:
69 : * |A| . (X) = V . |B| . (X)
70 : * CORRESPONDANT AU PLUS PETIT V SOLUTION (EN VALEUR ABSOLUE) SOUS
71 : * RESERVE QUE LE (X)1 DE DEPART N'EST PAS B-ORTHOGONAL AU (X)
72 : * SOLUTION.
73 : *
74 : * SOUS-PROGRAMMES APPELES:
75 : * ------------------------
76 : *
77 : * DTCHPO, ITINV1, XTMX, YTX1, RAYLEI, (?), VRFMOD,DEPGE2 ,DTCHPM
78 : *
79 : * AUTEUR, DATE DE MODIFICATION:
80 : * -------------------------
81 : *
82 : * C. LE BIDEAU JUILLET 2001
83 : * Benoit Prabel mars 2009
84 : *
85 : * LANGAGE:
86 : * --------
87 : *
88 : * FORTRAN77 + ESOPE
89 : *
90 : ************************************************************************
91 : *
92 : -INC CCOPTIO
93 : -INC SMLMOTS
94 : -INC SMCHPOI
95 : -INC CCHAMP
96 : *
97 : REAL*8 PROPRE(*), XRVP, XIVP
98 : REAL*8 TRNW1(100),TREPS(100)
99 : *
100 : COMMON/CITINV/ NBITER,IACCEL,NUMAC,IPX2,IPX0,IPX1,IPBX1,
101 : C IBBX1,IBBX2,ITPRO,DIFREL
102 : *
103 : LOGICAL CONVRG
104 : *
105 : PARAMETER (INFINI = 9999)
106 : c precision sur la norme 2
107 : * (rem: pour matrices sym, la precision est de 1.D-5 sur la norme infinie => on est donc bcp + exigeant)
108 : PARAMETER (XPREC1 = 1.D-8)
109 : PARAMETER (XPREC2 = 1.D-8)
110 : *
111 : IF (IIMPI .EQ. 747) THEN
111 : CALL GIBTEM(XKT)
112 : INTERR(1)=XKT
112 : CALL ERREUR(-259)
113 : ENDIF
114 : *
115 : * PREPARATION DES ITERATIONS:
116 : IBBX1=0
117 : IBBX2 = IPMX
118 : IPX2 = IPRX
119 : NUMAC = 100
120 : * bp: on n'accelere pas la convergence pour l instant
121 : * car avec des matrices non-sym, la direction du vecteur tourne
122 : * => il faudra adapter l'acceleration au cas complexe
123 : * car parfois cvge tres lente...
124 : NBITER = 0
125 : IACCEL = 0
126 : NUMXBX = -10
127 : X1BX1 = 1.D10
128 : IPLMOX=0
129 : IPLMOY=0
130 : RNW1=1.
131 : REPS=1.
132 : DIFREL = 1.E10
133 : CONVRG = .false.
134 : C
135 : C PREPARATION DES TABLEAUX DONNANT LA CORRESPONDANCE DES NOMS
136 : C D INCONNUE DANS X ET MX STOCKE DANS UN LIST MOT
137 : C
137 : CALL CORRSP(ipa,IPRX,IPMX,IPLMOX,IPLMOY)
138 : C
139 : * -- DEBUT DES ITERATIONS INVERSES -------------------
140 : *
141 : DO 5000 I = 1,ITERMX
142 :
142 : IF (IBBX1.NE.0) CALL DTCHPO(IBBX1)
143 : IBBX1 = IBBX2
144 :
144 : IF (I.GE.2) CALL COPIE2 (IPX1,IPREC)
145 : *** réalisation d'une iteration inverse (avec une eventuelle acceleration)
145 : CALL ITINV1 (IPA,IPB)
146 : IF (IERR .NE. 0) RETURN
147 :
148 : *** on ne teste la convergence qu'apres 3 iterations
149 : IF (I .GT. 3) THEN
150 : *
151 : * ON RECUPERE DEUX VECTEURS IPREC ET IPX1
152 : * 1ER VECTEUR DE BASE
152 : CALL XTY1 (IPREC,IPREC,IPLMOX,IPLMOX,RNNEC)
153 : RNEC = SQRT( ABS(RNNEC))
153 : CALL MUCHPO (IPREC, RNEC, IV, -1)
154 : *
155 : * 2EME VECTEUR DE BASE
155 : CALL XTY1 (IPX1, IV, IPLMOX, IPLMOX, RNX1)
155 : CALL ADCHPO (IPX1, IV, IW1, 1.D0, -RNX1)
155 : CALL XTY1 (IW1, IW1, IPLMOX, IPLMOX, RNNW1)
156 : RNW1 = SQRT( ABS(RNNW1))
157 : TRNW1(I)=RNW1
158 : *
159 : * SI 2E VECTEUR DE BASE NUL : 1 VP REELLE simple
160 : * IF ((I .GE. 3) .AND. (RNW1 .LT. XPREC1)) THEN
161 : IF ((I .GE. 3) .AND. (RNW1 .LT. (XPREC1*RNEC))) THEN
162 : CONVRG = .TRUE.
163 : GOTO 3001
164 : END IF
164 : CALL MUCHPO (IW1, RNW1, IW, -1)
165 :
166 : * CRITERE DE CONVERGENCE vers 1 paire de vp complexes conjuguées
166 : CALL XTY1 (IPX2, IV, IPLMOX, IPLMOX, P1)
166 : CALL XTY1 (IPX2, IW, IPLMOX, IPLMOX, P2)
167 : IF (IERR .NE. 0 ) RETURN
167 : CALL ADCHPO (IV, IW, IPX4, P1, P2)
167 : CALL ADCHPO (IPX2, IPX4, IEPS, 1.D0, -1.D0)
167 : CALL XTY1 (IEPS, IEPS, IPLMOX, IPLMOX, RNEPS)
168 : IF (IERR .NE. 0 ) RETURN
169 : REPS = SQRT(ABS(RNEPS))
170 : TREPS(I)=REPS
171 : *
172 : * SI IPX2 appartient au s.e. engendré par (IV,IW) : RETOUR VALEURS PROPRES
173 : * IF (RNEPS .LT. 1.E-6) THEN
174 : * IF (REPS .LT. (XPREC2)) THEN
175 : IF (REPS .LT. (XPREC1*RNEC)) THEN
175 : CALL RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA,
176 : & IPRX, IPIX, ICVG)
177 : * IF ( ABS(DELTA) .GE. 1.D-14) THEN
178 : IF ( ABS(DELTA) .GE. 1.D-10) THEN
179 : IF (DELTA .LT. 0.D0) THEN
180 : * cas d une paire de valeur propre Complexe
181 : CONVRG = .TRUE.
182 : GOTO 4000
183 : ELSE
184 : * cas d une valeur propre Réelle simple
185 : ** chgt de stratégie concernant l'acceleration de cvgce
186 : * if(NUMAC .ne. 5) then
187 : * NUMAC = 5
188 : * IACCEL = NUMAC - 1
189 : * endif
190 : * GOTO 4900
191 : * cas d une valeur propre Réelle simple pas encore assez bien calculée!
192 : if(ICVG .eq. 0) goto 4900
193 : CONVRG = .TRUE.
194 : * cas de 2 valeurs propre Réelle Double
195 : if(ICVG .eq. 2) GOTO 3002
196 : * cas de 2 valeurs propre Réelle simples de meme module!
197 : if(ICVG .eq. 3) GOTO 3003
198 : END IF
199 : ELSE
200 : * cas d une valeur propre Réelle double
201 : CONVRG = .TRUE.
202 : GOTO 3002
203 : END IF
204 : END IF
205 :
206 :
207 : 4900 CONTINUE
208 : * si l'acceleration de cvgce n'ameliore pas le taux de cvgce, alors on arrete
209 : * IF((NUMAC .ne. 100) .and. (IACCEL .EQ. 0)) then
210 : * XTAU2 = (LOG(TREPS(I))) - (LOG(TREPS(I-1)))
211 : * XTAU1 = (LOG(TREPS(I-1))) - (LOG(TREPS(I-2)))
212 : * if(XTAU2 .ge. XTAU1) then
213 : * NUMAC = 100
214 : * endif
215 : * ENDIF
216 : * segdes,MPOVAL,MSOUPO,MCHPOI
217 :
218 :
219 : * si on a atteint le nbre maxi d iterations on retourne ce que l 'on a
220 : IF (NBITER .GE. ITERMX) THEN
221 : CONVRG = .FALSE.
222 : * GOTO 302
222 : CALL RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA,
223 : & IPRX, IPIX, ICVG)
224 : IF ( ABS(DELTA) .GE. 1.D-10) THEN
225 : IF (DELTA .LT. 0.D0) THEN
226 : GOTO 4000
227 : ELSE
228 : if(XIVP .ne. 0.) GOTO 3002
229 : GOTO 3001
230 : END IF
231 : ELSE
232 : GOTO 3002
233 : END IF
234 :
235 : END IF
236 :
237 : NBITER = I
238 :
239 : *** cas ou I < 3 iterations
240 : ELSE
241 : TRNW1(I)=RNW1
242 : TREPS(I)=REPS
243 :
244 : END IF
245 : *** fin des tests de convergence réalisés apres 3 iterations
246 : *
247 : *
248 : 5000 CONTINUE
249 : * write(*,*)'-- FIN DE LA BOUCLE DES ITERATIONS INVERSES --------'
249 : CALL DTCHPO(IPX1)
250 : GOTO 4000
251 :
252 :
253 : 3001 CONTINUE
254 : * cas d une valeur propre Réelle (simple?) '
255 : * ----------------------------
256 : TRNW1(I)=RNW1
257 : TREPS(I)=1.
258 : * -on n'est pas passé par Rayleigh : il faut faire le travail ici
259 : * normalisation du vecteur propre IPX2
259 : CALL MOTS1(IPLMOT,MOTCLE)
259 : call NORMA1(IPX2,IPLMOT,MOTCLE,IPX2)
260 : * vecteur propre
261 : IPRX = IPX2
262 : IPIX = 0
263 : * produits
263 : CALL MUCPRI(IPX2, IPA, IPAX2)
263 : CALL XTY1(IPX2, IPAX2, IPLMOX, IPLMOY, X2AX2)
264 : IF (IERR .NE. 0 ) RETURN
264 : CALL MUCPRI(IPX2, IPB, IPBX2)
264 : CALL XTY1(IPX2, IPBX2, IPLMOX, IPLMOY, X2BX2)
265 : IF (IERR .NE. 0 ) RETURN
266 : * valeur propre (simple)
267 : XRVP = X2AX2 / X2BX2
268 : XIVP = 0.D0
269 : PROPRE(1) = XRVP
270 : PROPRE(6) = XIVP
271 : * MASSES GENERALISEES
272 : PROPRE(2) = X2BX2
273 : PROPRE(7) = 0.D0
274 : C INTRODUCTION DES COEF. PI OU 2PI EVENTUELS + calcul DEPGEN
274 : CALL MASGEN(IPX2,PROPRE)
274 : CALL DEPGEN( IPB, IPX2,PROPRE,IBBX2,IPLMOX,IPLMOY)
275 : * write(6,*) 'itinvc trouve lambda=',XRVP
276 : GOTO 302
277 :
278 : 3002 CONTINUE
279 : * cas d une valeur propre Réelle double
280 : * (calculée par Rayleigh) ----
281 : * valeur propre (double)
282 : PROPRE(1) = XRVP
283 : PROPRE(6) = XIVP
284 : * masse généralisée (2: 1er vecteur reel , 7: 2eme vecteur reel)
284 : CALL MUCPRI(IPRX, IPB, IPBRX)
284 : CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
284 : CALL MUCPRI(IPIX, IPB, IPBIX)
284 : CALL XTY1(IPIX, IPBIX, IPLMOX, IPLMOY, XIBIX)
285 : PROPRE(2) = XRBRX
286 : PROPRE(7) = XIBIX
287 : * write(6,*) 'itinvc trouve lambda=',XRVP,' et ',XIVP
288 : GOTO 302
289 : * ...-> reste a prevoir ce cas dans crebas.eso...
290 :
291 : 3003 CONTINUE
292 : * cas de 2 valeur propre Réelle simple de meme module
293 : * (calculée par Rayleigh) ----
294 : PROPRE(1) = XRVP
295 : PROPRE(6) = 0.
296 : * masse généralisée (2: 1er vecteur reel , 7: 2eme vecteur reel)
296 : CALL MUCPRI(IPRX, IPB, IPBRX)
296 : CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
297 : PROPRE(2) = XRBRX
298 : PROPRE(7) = 0.
299 : C INTRODUCTION DES COEF. PI OU 2PI EVENTUELS + calcul DEPGEN
299 : CALL MASGEN(IPRX,PROPRE)
299 : * CALL DEPGEN( IPB, IPRX,PROPRE,IBBX2,IPLMOX,IPLMOY)
300 : * write(6,*) 'itinvc trouve lambda=',XRVP,' (et ',(-1.*XRVP),')'
301 : GOTO 302
302 :
303 : 4000 CONTINUE
304 : * cas d une paire de valeur propre Complexe'
305 : * (calculée par Rayleigh) ----
306 : * valeur propre lamdba
307 : PROPRE(1) = XRVP
308 : PROPRE(6) = XIVP
309 : * masse complexe générlaisée
309 : CALL MUCPRI(IPRX, IPB, IPBRX)
309 : CALL XTY1(IPRX, IPBRX, IPLMOX, IPLMOY, XRBRX)
309 : CALL XTY1(IPIX, IPBRX, IPLMOX, IPLMOY, XIBRX)
309 : CALL MUCPRI(IPIX, IPB, IPBIX)
309 : CALL XTY1(IPRX, IPBIX, IPLMOX, IPLMOY, XRBIX)
309 : CALL XTY1(IPIX, IPBIX, IPLMOX, IPLMOY, XIBIX)
310 : PROPRE(2) = XRBRX - XIBIX
311 : PROPRE(7) = XIBRX + XRBIX
312 : * write(*,*) 'itinvc trouve lambda=',XRVP,' + i',XIVP
313 : GOTO 302
314 :
315 : *
316 : * IMPRESSIONS
317 : *
318 : 302 CONTINUE
319 : *
320 : IF (IIMPI.EQ.2) WRITE (IOIMP,2000) NBITER
321 : 2000 FORMAT (//,1X,I3,' ITERATIONS INVERSES ONT ETE EFFECTUEES.'///)
322 : IF (IIMPI.EQ.30) WRITE(IOIMP,1000) (PROPRE(I),I=1,10)
323 : 1000 FORMAT(/10X,'SBR ITINVC',/10X,5(E12.5,1X))
324 : *
325 : * write(6,*) 'i Residu1 Residu2'
326 : * do ii=1,I
327 : * write(6,*) ii,(TRNW1(ii)),(TREPS(ii))
328 : * enddo
329 :
330 : *
331 : * SUPPRESSION DES 'CHPOINT' DE TRAVAIL:
332 : IF ( (IACCEL + 1) .EQ. NUMAC) THEN
332 : CALL DTCHPO (IPX0)
333 : END IF
333 : CALL DTCHPO (IBBX2)
334 : MLMOTS =IPLMOX
335 : MLMOT1 =IPLMOY
336 : SEGSUP MLMOTS,MLMOT1
337 : *
338 : IF (IIMPI .EQ. 30) THEN
338 : CALL GIBTEM(XKT)
339 : INTERR(1)=XKT
339 : CALL ERREUR(-259)
339 : CALL GIBTEM(XKT)
340 : INTERR(1)=XKT
340 : CALL ERREUR(-259)
341 : END IF
342 : *
343 : END
344 :
345 :
346 :
© Cast3M 2003 - All rights reserved.
Disclaimer