Télécharger qzredu.eso

Retour à la liste

Numérotation des lignes :

  1. C QZREDU SOURCE BP208322 12/10/03 21:15:29 7517
  2. SUBROUTINE QZREDU (MATA, MATB, EPS1, CALV, MATZ, CERR)
  3. c SUBROUTINE QZREDU (XMATA, XMATB, EPS1, CALV, XMATZ,NM,N, CERR)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. *
  7. *************************************************************************
  8. * *
  9. * MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE *
  10. * D'UNE MATRICE DE HESSENBERG SUPERIEURE (MATA) *
  11. * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB *
  12. * ( DEUXIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) *
  13. * _____________________________________________________________________ *
  14. * *
  15. * *
  16. * DATE : le 27 Mars 1995 *
  17. * AUTEURS : J.ANTUNES , Nicolas BENECH *
  18. * _____________________________________________________________________ *
  19. * *
  20. * MODULE(S) APPELANT(S) : VIBRAC *
  21. * *
  22. * MODULE(S) APPELE(S) : *
  23. * _____________________________________________________________________ *
  24. * *
  25. * EN ENTREE : *
  26. * - MATA : matrice reelle Hessenberg superieure (XMATRI) *
  27. * - MATB : matrice reelle triangulaire superieure (XMATRI) *
  28. * - EPS1 : tolerance pour les elements negligeables *
  29. * si EPS1 = 0. ou EPS1 < 0., un element est neglige *
  30. * s'il est inferieur a l'erreur d'arrondi multiplie *
  31. * par la norme de la matrice *
  32. * si EPS1 > 0., un element est neglige s'il est inferieur *
  33. * a EPS1 multiplie par la norme de la matrice *
  34. * (plus rapide mais moins fiable) *
  35. * - CALV : indique si l'on calcule ou non les vecteurs propres *
  36. * *
  37. * EN SORTIE : *
  38. * - MATA : matrice quasi-triangulaire superieure *
  39. * - MATB : matrice sous forme triangulaire superieure *
  40. * - MATZ : produit des transformations effectuees (si CALV=.TRUE.)*
  41. * - CERR : code des erreurs *
  42. * CERR = 0 : retour normal *
  43. * CERR = J : plus de 30*N iterations pour le calcul *
  44. * de la valeur propre J *
  45. * *
  46. *************************************************************************
  47. *
  48. -INC CCOPTIO
  49. -INC SMRIGID
  50. C
  51. INTEGER I , J , K , L , EN , K1 , K2 , LD , LL , L1 , NA ,
  52. & ISH , ITN , ITS , KM1 , LM1 , ENM2 , CERR , LOR1 ,
  53. & ENORN
  54. C
  55. LOGICAL CALV , NOTLAS
  56. C-----------------------------------------------------------------------
  57. POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  58. c REAL*8 XMATA(NM,N),XMATB(NM,N),XMATZ(NM,N)
  59. C
  60. c SEGACT , MATA*mod
  61. c SEGACT , MATB*mod
  62. c SEGACT , MATZ*mod
  63. NM=MATA.RE(/1)
  64. N=MATA.RE(/2)
  65. C-----------------------------------------------------------------------
  66. C
  67. if(iimpi.ge.6) then
  68. write(ioimp,*) 'A='
  69. do iou=1,N
  70. write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N)
  71. enddo
  72. write(ioimp,*) 'B='
  73. do iou=1,N
  74. write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N)
  75. enddo
  76. endif
  77. C***********************************************************************
  78. C***********************************************************************
  79. C***********************************************************************
  80. C
  81. C=======================================================================
  82. C========================= DEBUT DES CALCULS ===========================
  83. C=======================================================================
  84. C
  85. CERR = 0
  86. C
  87. C-----------------------------------------------------------------------
  88. C---------------------- CALCUL DE EPSA ET EPSB ----------------------
  89. C-----------------------------------------------------------------------
  90. C
  91. ANORM = 0.0D0
  92. BNORM = 0.0D0
  93. C
  94. DO 30 I = 1, N
  95. ANI = 0.0D0
  96. IF ( I .NE. 1 ) ANI = ABS(MATA.RE(I,I-1,1))
  97. BNI = 0.0D0
  98. C
  99. DO 20 J = I, N
  100. ANI = ANI + ABS(MATA.RE(I,J,1))
  101. BNI = BNI + ABS(MATB.RE(I,J,1))
  102. 20 CONTINUE
  103. C
  104. IF ( ANI .GT. ANORM ) ANORM = ANI
  105. IF ( BNI .GT. BNORM ) BNORM = BNI
  106. 30 CONTINUE
  107. C
  108. IF ( ANORM .EQ. 0.0D0 ) ANORM = 1.0D0
  109. IF ( BNORM .EQ. 0.0D0 ) BNORM = 1.0D0
  110. EP = EPS1
  111. IF ( EP .GT. 0.0D0 ) GOTO 50
  112. C
  113. EP = EPSLON(1.0D0)
  114. IF (EP .GT. 0.0D0) GO TO 50
  115. C ********** COMPUTE ROUNDOFF LEVEL IF EPSLON(1.) IS ZERO **********
  116. EP = 1.0D0
  117. 40 EP = EP / 2.0
  118. IF (1.0D0 + EP .GT. 1.0D0) GO TO 40
  119. 50 EPSA = EP * ANORM
  120. EPSB = EP * BNORM
  121. C ********** REDUCE A TO QUASI-TRIANGULAR FORM, WHILE
  122. C KEEPING B TRIANGULAR **********
  123. LOR1 = 1
  124. ENORN = N
  125. EN = N
  126. c ITN = 30*N
  127. c bp 2012-10-02 : on augmente le nombre d'iterations max
  128. ITN = 60*N
  129. C
  130. C-----------------------------------------------------------------------
  131. C------------------------- DEBUT DU PAS QZ ---------------------------
  132. C-----------------------------------------------------------------------
  133. C
  134. 60 IF ( EN .LE. 2 ) GOTO 1001
  135. IF ( .NOT. CALV ) ENORN = EN
  136. ITS = 0
  137. NA = EN - 1
  138. ENM2 = NA - 1
  139. 70 ISH = 2
  140. C
  141. C-----------------------------------------------------------------------
  142. C------------------- VERIFICATION DE LA CONVERGENCE --------------------
  143. C-----------------------------------------------------------------------
  144. C
  145. DO 80 LL = 1, EN
  146. LM1 = EN - LL
  147. L = LM1 + 1
  148. IF ( L .EQ. 1 ) GOTO 95
  149. IF ( ABS(MATA.RE(L,LM1,1)) .LE. EPSA ) GOTO 90
  150. 80 CONTINUE
  151. C
  152. 90 MATA.RE(L,LM1,1) = 0.0D0
  153. IF ( L .LT. NA ) GOTO 95
  154. C ********** 1-BY-1 OR 2-BY-2 BLOCK ISOLATED **********
  155. EN = LM1
  156. GOTO 60
  157. C ********** CHECK FOR SMALL TOP OF B **********
  158. 95 LD = L
  159. 100 L1 = L + 1
  160. C
  161. B11 = MATB.RE(L,L,1)
  162. IF ( ABS(B11) .GT. EPSB ) GOTO 120
  163. MATB.RE(L,L,1) = 0.0D0
  164. C
  165. S = ABS(MATA.RE(L,L,1)) + ABS(MATA.RE(L1,L,1))
  166. U1 = MATA.RE(L,L,1) / S
  167. U2 = MATA.RE(L1,L,1) / S
  168. C
  169. C
  170. R = SIGN(SQRT(U1*U1+U2*U2),U1)
  171. V1 = -(U1 + R) / R
  172. V2 = -U2 / R
  173. U2 = V2 / V1
  174. C
  175. DO 110 J = L, ENORN
  176. T = MATA.RE(L,J,1) + U2 * MATA.RE(L1,J,1)
  177. MATA.RE(L,J,1) = MATA.RE(L,J,1) + T * V1
  178. MATA.RE(L1,J,1) = MATA.RE(L1,J,1) + T * V2
  179. T = MATB.RE(L,J,1) + U2 * MATB.RE(L1,J,1)
  180. MATB.RE(L,J,1) = MATB.RE(L,J,1) + T * V1
  181. MATB.RE(L1,J,1) = MATB.RE(L1,J,1) + T * V2
  182. 110 CONTINUE
  183. C
  184. IF ( L .NE. 1 ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1)
  185. LM1 = L
  186. L = L1
  187. GOTO 90
  188. 120 CONTINUE
  189. A11 = MATA.RE(L,L,1) / B11
  190. A21 = MATA.RE(L1,L,1) / B11
  191. IF ( ISH .EQ. 1 ) GOTO 141
  192. C
  193. C-----------------------------------------------------------------------
  194. C---------------------- STRATEGIE DE L ITERATION -----------------------
  195. C-----------------------------------------------------------------------
  196. C
  197. C ********** ITERATION STRATEGY **********
  198. IF ( ITN .EQ. 0 ) GOTO 1000
  199. c bp 2012-10-02 : on change le moment d'appel au "ad hoc shift" dont on
  200. c ne trouve pas de justification mathematique mais on la cale
  201. c numeriquement sur un cas test problematique : brasero_tube3D.dgibi
  202. c l'ideal serait d'avoir un critere qui permet de declencher de shift
  203. c (ou bien de passer a linpack et DHGEQZ)
  204. c IF ( ITS .EQ. 10 ) GOTO 155
  205. IF ( ITS .EQ. 10*N) GOTO 155
  206.  
  207. C
  208. C ********** DETERMINE TYPE OF SHIFT **********
  209. B22 = MATB.RE(L1,L1,1)
  210. IF ( ABS(B22) .LT. EPSB ) B22 = EPSB
  211. B33 = MATB.RE(NA,NA,1)
  212. IF ( ABS(B33) .LT. EPSB ) B33 = EPSB
  213. B44 = MATB.RE(EN,EN,1)
  214. IF ( ABS(B44) .LT. EPSB ) B44 = EPSB
  215. A33 = MATA.RE(NA,NA,1) / B33
  216. A34 = MATA.RE(NA,EN,1) / B44
  217. A43 = MATA.RE(EN,NA,1) / B33
  218. A44 = MATA.RE(EN,EN,1) / B44
  219. B34 = MATB.RE(NA,EN,1) / B44
  220. T = 0.5D0 * (A43 * B34 - A33 - A44)
  221. R = T * T + A34 * A43 - A33 * A44
  222. IF ( R .LT. 0.0D0 ) GOTO 150
  223. C
  224. C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A **********
  225. ISH = 1
  226. R = SQRT(R)
  227. SH = -T + R
  228. S = -T - R
  229. IF ( ABS(S-A44) .LT. ABS(SH-A44) ) SH = S
  230. C
  231. C ********** LOOK FOR TWO CONSECUTIVE SMALL
  232. C SUB-DIAGONAL ELEMENTS OF A.
  233. C FOR L=EN-2 STEP -1 UNTIL LD DO -- **********
  234. DO 130 LL = LD, ENM2
  235. L = ENM2 + LD - LL
  236. IF ( L .EQ. LD ) GOTO 140
  237. LM1 = L - 1
  238. L1 = L + 1
  239. T = MATA.RE(L,L,1)
  240. IF ( ABS(MATB.RE(L,L,1)) .GT. EPSB )
  241. & T = T - SH * MATB.RE(L,L,1)
  242. IF (ABS(MATA.RE(L,LM1,1)).LE.ABS(T/MATA.RE(L1,L,1))*EPSA)
  243. & GOTO 100
  244. 130 CONTINUE
  245. 140 continue
  246. C
  247. 141 A1 = A11 - SH
  248. A2 = A21
  249. IF ( L .NE. LD ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1)
  250. GOTO 160
  251. C
  252. C ********** DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A **********
  253. 150 continue
  254. A12 = MATA.RE(L,L1,1) / B22
  255. A22 = MATA.RE(L1,L1,1) / B22
  256. B12 = MATB.RE(L,L1,1) / B22
  257. A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11)
  258. & / A21 + A12 - A11 * B12
  259. A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11)
  260. & + A43 * B34
  261. A3 = MATA.RE(L1+1,L1,1) / B22
  262. GOTO 160
  263. C
  264. C ********** AD HOC SHIFT **********
  265. 155 A1 = 0.0D0
  266. A2 = 1.0D0
  267. A3 = 1.1605D0
  268. 160 ITS = ITS + 1
  269. ITN = ITN - 1
  270. if(iimpi.ge.6) write(*,*)'iteration ',ITS,' A1,A2,A3=',A1,A2,A3
  271. IF ( .NOT. CALV ) LOR1 = LD
  272. C
  273. C-----------------------------------------------------------------------
  274. C------------------------- BOUCLE PRINCIPALE ---------------------------
  275. C-----------------------------------------------------------------------
  276. C
  277. DO 260 K = L, NA
  278. NOTLAS = K .NE. NA .AND. ISH .EQ. 2
  279. K1 = K + 1
  280. K2 = K + 2
  281. KM1 = MAX0(K-1,L)
  282. LL = MIN0(EN,K1+ISH)
  283. IF ( NOTLAS ) GOTO 190
  284. C
  285. C ********** ZERO A(K+1,K-1) **********
  286. IF ( K .EQ. L ) GOTO 170
  287. A1 = MATA.RE(K,KM1,1)
  288. A2 = MATA.RE(K1,KM1,1)
  289. 170 S = ABS(A1) + ABS(A2)
  290. IF ( S .EQ. 0.0D0 ) GOTO 70
  291. U1 = A1 / S
  292. U2 = A2 / S
  293. R = SIGN(SQRT(U1*U1+U2*U2),U1)
  294. V1 = -(U1 + R) / R
  295. V2 = -U2 / R
  296. U2 = V2 / V1
  297. C
  298. DO 180 J = KM1, ENORN
  299. T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1)
  300. MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1
  301. MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2
  302. T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1)
  303. MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1
  304. MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2
  305. 180 CONTINUE
  306. C
  307. IF ( K .NE. L ) MATA.RE(K1,KM1,1) = 0.0D0
  308. GOTO 240
  309. C
  310. C ********** ZERO A(K+1,K-1) AND A(K+2,K-1) **********
  311. 190 IF ( K .EQ. L ) GOTO 200
  312. A1 = MATA.RE(K,KM1,1)
  313. A2 = MATA.RE(K1,KM1,1)
  314. A3 = MATA.RE(K2,KM1,1)
  315. 200 S = ABS(A1) + ABS(A2) + ABS(A3)
  316. IF ( S .EQ. 0.0D0 ) GOTO 260
  317. U1 = A1 / S
  318. U2 = A2 / S
  319. U3 = A3 / S
  320. R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
  321. V1 = -(U1 + R) / R
  322. V2 = -U2 / R
  323. V3 = -U3 / R
  324. U2 = V2 / V1
  325. U3 = V3 / V1
  326. C
  327. DO 210 J = KM1, ENORN
  328. T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1)
  329. & + U3 * MATA.RE(K2,J,1)
  330. MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1
  331. MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2
  332. MATA.RE(K2,J,1) = MATA.RE(K2,J,1) + T * V3
  333. T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1)
  334. & + U3 * MATB.RE(K2,J,1)
  335. MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1
  336. MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2
  337. MATB.RE(K2,J,1) = MATB.RE(K2,J,1) + T * V3
  338. 210 CONTINUE
  339. C
  340. IF ( K .EQ. L ) GOTO 220
  341. MATA.RE(K1,KM1,1) = 0.0D0
  342. MATA.RE(K2,KM1,1) = 0.0D0
  343. C
  344. C ********** ZERO B(K+2,K+1) AND B(K+2,K) **********
  345. 220 continue
  346. S = ABS(MATB.RE(K2,K2,1)) + ABS(MATB.RE(K2,K1,1))
  347. & + ABS(MATB.RE(K2,K,1))
  348. IF ( S .EQ. 0.0D0 ) GOTO 240
  349. U1 = MATB.RE(K2,K2,1) / S
  350. U2 = MATB.RE(K2,K1,1) / S
  351. U3 = MATB.RE(K2,K,1) / S
  352. R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1)
  353. V1 = -(U1 + R) / R
  354. V2 = -U2 / R
  355. V3 = -U3 / R
  356. U2 = V2 / V1
  357. U3 = V3 / V1
  358. C
  359. DO 230 I = LOR1, LL
  360. T = MATA.RE(I,K2,1) + U2 * MATA.RE(I,K1,1)
  361. & + U3 * MATA.RE(I,K,1)
  362. MATA.RE(I,K2,1) = MATA.RE(I,K2,1) + T * V1
  363. MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V2
  364. MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V3
  365. T = MATB.RE(I,K2,1) + U2 * MATB.RE(I,K1,1)
  366. & + U3 * MATB.RE(I,K,1)
  367. MATB.RE(I,K2,1) = MATB.RE(I,K2,1) + T * V1
  368. MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V2
  369. MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V3
  370. 230 CONTINUE
  371. C
  372. MATB.RE(K2,K,1) = 0.0D0
  373. MATB.RE(K2,K1,1) = 0.0D0
  374. IF ( .NOT. CALV ) GOTO 240
  375. C
  376. DO 235 I = 1, N
  377. T = MATZ.RE(I,K2,1) + U2 * MATZ.RE(I,K1,1)
  378. & + U3 * MATZ.RE(I,K,1)
  379. MATZ.RE(I,K2,1) = MATZ.RE(I,K2,1) + T * V1
  380. MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V2
  381. MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V3
  382. 235 CONTINUE
  383. C
  384. C ********** ZERO B(K+1,K) **********
  385. 240 CONTINUE
  386. S = ABS(MATB.RE(K1,K1,1)) + ABS(MATB.RE(K1,K,1))
  387. IF ( S .EQ. 0.0D0 ) GOTO 260
  388. U1 = MATB.RE(K1,K1,1) / S
  389. U2 = MATB.RE(K1,K,1) / S
  390. R = SIGN(SQRT(U1*U1+U2*U2),U1)
  391. V1 = -(U1 + R) / R
  392. V2 = -U2 / R
  393. U2 = V2 / V1
  394. C
  395. DO 250 I = LOR1, LL
  396. T = MATA.RE(I,K1,1) + U2 * MATA.RE(I,K,1)
  397. MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V1
  398. MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V2
  399. T = MATB.RE(I,K1,1) + U2 * MATB.RE(I,K,1)
  400. MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V1
  401. MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V2
  402. 250 CONTINUE
  403. C
  404. MATB.RE(K1,K,1) = 0.0D0
  405. C
  406. IF ( .NOT. CALV ) GOTO 260
  407. C
  408. DO 255 I = 1, N
  409. T = MATZ.RE(I,K1,1) + U2 * MATZ.RE(I,K,1)
  410. MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V1
  411. MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V2
  412. 255 CONTINUE
  413. C
  414. 260 CONTINUE
  415. C
  416. GOTO 70
  417. C
  418. C-----------------------------------------------------------------------
  419. C---- INDICATION D ERREUR SI PAS DE CONVERGENCE EN 30*N ITERATIONS -----
  420. C-----------------------------------------------------------------------
  421. C
  422. 1000 CERR = EN
  423. C
  424. C-----------------------------------------------------------------------
  425. C------ ON SAUVE EPSB POUR UTILISATION DANS QZVALP ET QZVECP ------
  426. C-----------------------------------------------------------------------
  427. C
  428. 1001 IF ( N .LE. 1 ) GOTO 1002
  429. MATB.RE(N,1,1) = EPSB
  430. c SEGDES , MATA,MATB,MATZ
  431. C
  432. C=======================================================================
  433. C========================== FIN DES CALCULS ============================
  434. C=======================================================================
  435. C
  436. if(iimpi.ge.4) then
  437. if(CERR.eq.0) then
  438. write(ioimp,*)'convergence en',ITS,' iterations ITN=',ITN
  439. else
  440. write(ioimp,*)'NON convergence en',ITS,' iterations ITN=',ITN
  441. write(ioimp,*) 'A='
  442. do iou=1,N
  443. write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N)
  444. enddo
  445. write(ioimp,*) 'B='
  446. do iou=1,N
  447. write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N)
  448. enddo
  449. endif
  450. endif
  451.  
  452. 1002 RETURN
  453. END
  454.  
  455.  
  456.  
  457.  
  458.  

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