Télécharger dneigh.eso

Retour à la liste

Numérotation des lignes :

  1. C DNEIGH SOURCE BP208322 15/10/13 21:15:46 8670
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dneigh
  6. c
  7. c\Description:
  8. c Compute the eigenvalues of the current upper Hessenberg matrix
  9. c and the corresponding Ritz estimates given the current residual norm.
  10. c
  11. c\Usage:
  12. c call dneigh
  13. c ( RNORM, N, H, LDH, RITZR, RITZI, BOUNDS, Q, LDQ, WORKL, IERR )
  14. c
  15. c\Arguments
  16. c RNORM Double precision scalar. (INPUT)
  17. c Residual norm corresponding to the current upper Hessenberg
  18. c matrix H.
  19. c
  20. c N Integer. (INPUT)
  21. c Size of the matrix H.
  22. c
  23. c H REAL*8 N by N array. (INPUT)
  24. c H contains the current upper Hessenberg matrix.
  25. c
  26. c LDH Integer. (INPUT)
  27. c Leading dimension of H exactly as declared in the calling
  28. c program.
  29. c
  30. c RITZR, Double precision arrays of length N. (OUTPUT)
  31. c RITZI On output, RITZR(1:N) (resp. RITZI(1:N)) contains the real
  32. c (respectively imaginary) parts of the eigenvalues of H.
  33. c
  34. c BOUNDS Double precision array of length N. (OUTPUT)
  35. c c On output, BOUNDS contains the Ritz estimates associated with
  36. c the eigenvalues RITZR and RITZI. This is equal to RNORM
  37. c times the last components of the eigenvectors corresponding
  38. c to the eigenvalues in RITZR and RITZI.
  39. c
  40. c Q REAL*8 N by N array. (WORKSPACE)
  41. c Workspace needed to store the eigenvectors of H.
  42. c
  43. c LDQ Integer. (INPUT)
  44. c Leading dimension of Q exactly as declared in the calling
  45. c program.
  46. c
  47. c WORKL Double precision work array of length N**2 + 3*N. (WORKSPACE)
  48. c Private (replicated) array on each PE or array allocated on
  49. c the front end. This is needed to keep the full Schur form
  50. c of H and also in the calculation of the eigenvectors of H.
  51. c
  52. c IERR Integer. (OUTPUT)
  53. c Error exit flag from dlaqrb or dtrevc.
  54. c
  55. c\EndDoc
  56. c
  57. c-----------------------------------------------------------------------
  58. c
  59. c\BeginLib
  60. c
  61. c\Local variables:
  62. c xxxxxx real
  63. c
  64. c\Routines called:
  65. c dlaqrb ARPACK routine to compute the real Schur form of an
  66. c upper Hessenberg matrix and last row of the Schur vectors.
  67. c arscnd ARPACK utility routine for timing.
  68. c dmout ARPACK utility routine that prints matrices
  69. c dvout ARPACK utility routine that prints vectors.
  70. c dlacpy LAPACK matrix copy routine.
  71. c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
  72. c dtrevc LAPACK routine to compute the eigenvectors of a matrix
  73. c in upper quasi-triangular form
  74. c dgemv Level 2 BLAS routine for matrix vector multiplication.
  75. c dcopy Level 1 BLAS that copies one vector to another .
  76. c dnrm2 Level 1 BLAS that computes the norm of a vector.
  77. c dscal Level 1 BLAS that scales a vector.
  78. c
  79. c
  80. c\Author
  81. c Danny Sorensen Phuong Vu
  82. c Richard Lehoucq CRPC / Rice University
  83. c Dept. of Computational & Houston, Texas
  84. c Applied Mathematics
  85. c Rice University
  86. c Houston, Texas
  87. c
  88. c\Revision history:
  89. c xx/xx/92: Version ' 2.1'
  90. c
  91. c\SCCS Information: @(#)
  92. c FILE: neigh.F SID: 2.3 DATE OF SID: 4/20/96 RELEASE: 2
  93. c
  94. c\Remarks
  95. c None
  96. c
  97. c\EndLib
  98. c
  99. c-----------------------------------------------------------------------
  100. c
  101. subroutine dneigh (rnorm, n, h, ldh, ritzr, ritzi, bounds,
  102. & q, ldq, workl, ierr)
  103. c
  104. c %----------------------------------------------------%
  105. c | Include files for debugging and timing information |
  106. -INC TARTRAK
  107. c %----------------------------------------------------%
  108. c
  109. c
  110. c %------------------%
  111. c | Scalar Arguments |
  112. c %------------------%
  113. c
  114. integer ierr, n, ldh, ldq
  115. REAL*8
  116. & rnorm
  117. c
  118. c %-----------------%
  119. c | Array Arguments |
  120. c %-----------------%
  121. c
  122. REAL*8
  123. & bounds(n), h(ldh,n), q(ldq,n), ritzi(n), ritzr(n),
  124. & workl(n*(n+3))
  125. c
  126. c %------------%
  127. c | Parameters |
  128. c %------------%
  129. c
  130. REAL*8
  131. & one, zero
  132. parameter (one = 1.0D+0, zero = 0.0D+0)
  133. c
  134. c %------------------------%
  135. c | Local Scalars & Arrays |
  136. c %------------------------%
  137. c
  138. logical select(1)
  139. integer i, iconj, msglvl
  140. REAL*8
  141. & temp, vl(1)
  142. c
  143. c %----------------------%
  144. c | External Subroutines |
  145. c %----------------------%
  146. c
  147. external dcopy, dlacpy, dlaqrb, dtrevc, dvout, arscnd
  148. c
  149. c %--------------------%
  150. c | External Functions |
  151. c %--------------------%
  152. c
  153. REAL*8
  154. external dlapy2, dnrm2
  155. c
  156. c %---------------------%
  157. **c | Intrinsic Functions |
  158. **c %---------------------%
  159. **c
  160. ** intrinsic abs
  161. **c
  162. **c %-----------------------%
  163. **c | Executable Statements |
  164. c %-----------------------%
  165. c
  166. c
  167. c %-------------------------------%
  168. c | Initialize timing statistics |
  169. c | & message level for debugging |
  170. c %-------------------------------%
  171. c
  172. * call arscnd (t0)
  173.  
  174.  
  175. msglvl = mneigh
  176. c
  177. if (msglvl .gt. 2) then
  178. call dmout (logfil, n, n, h, ldh, ndigit,
  179. & '_neigh: Entering upper Hessenberg matrix H ')
  180. end if
  181. c
  182. c %-----------------------------------------------------------%
  183. c | 1. Compute the eigenvalues, the last components of the |
  184. c | corresponding Schur vectors and the full Schur form T |
  185. c | of the current upper Hessenberg matrix H. |
  186. c | dlaqrb returns the full Schur form of H in WORKL(1:N**2) |
  187. c | and the last components of the Schur vectors in BOUNDS. |
  188. c %-----------------------------------------------------------%
  189. c
  190. call dlacpy ('All', n, n, h, ldh, workl, n)
  191.  
  192. do 5 j = 1, n-1
  193. bounds(j) = zero
  194. 5 continue
  195. bounds(n) = one
  196.  
  197. call dlahqr(.true., .true., n, 1, n, workl, n, ritzr, ritzi, 1,
  198. & 1, bounds, 1, ierr)
  199.  
  200. c call dlaqrb (.true., n, 1, n, workl, n, ritzr, ritzi, bounds,
  201. c & ierr)
  202.  
  203.  
  204. if (ierr .ne. 0) go to 9000
  205. c
  206. if (msglvl .gt. 1) then
  207. call dvout (logfil, n, bounds, ndigit,
  208. & '_neigh: last row of the Schur matrix for H')
  209. end if
  210. c
  211. c %-----------------------------------------------------------%
  212. c | 2. Compute the eigenvectors of the full Schur form T and |
  213. c | apply the last components of the Schur vectors to get |
  214. c | the last components of the corresponding eigenvectors. |
  215. c | Remember that if the i-th and (i+1)-st eigenvalues are |
  216. c | complex conjugate pairs, then the real & imaginary part |
  217. c | of the eigenvector components are split across adjacent |
  218. c | columns of Q. |
  219. c %-----------------------------------------------------------%
  220. c
  221. call dtrevc ('R', 'A', select, n, workl, n, vl, n, q, ldq,
  222. & n, n, workl(n*n+1), ierr)
  223. c
  224. if (ierr .ne. 0) go to 9000
  225. c
  226. c %------------------------------------------------%
  227. c | Scale the returning eigenvectors so that their |
  228. c | euclidean norms are all one. LAPACK subroutine |
  229. c | dtrevc returns each eigenvector normalized so |
  230. c | that the element of largest magnitude has |
  231. c | magnitude 1; here the magnitude of a complex |
  232. c | number (x,y) is taken to be |x| + |y|. |
  233. c %------------------------------------------------%
  234. c
  235. iconj = 0
  236. do 10 i=1, n
  237. if ( abs( ritzi(i) ) .le. zero ) then
  238. c
  239. c %----------------------%
  240. c | Real eigenvalue case |
  241. c %----------------------%
  242. c
  243. temp = dnrm2( n, q(1,i), 1 )
  244. call dscal ( n, one / temp, q(1,i), 1 )
  245. else
  246. c
  247. c %-------------------------------------------%
  248. c | Complex conjugate pair case. Note that |
  249. c | since the real and imaginary part of |
  250. c | the eigenvector are stored in consecutive |
  251. c | columns, we further normalize by the |
  252. c | square root of two. |
  253. c %-------------------------------------------%
  254. c
  255. if (iconj .eq. 0) then
  256. temp = dlapy2( dnrm2( n, q(1,i), 1 ),
  257. & dnrm2( n, q(1,i+1), 1 ) )
  258. call dscal ( n, one / temp, q(1,i), 1 )
  259. call dscal ( n, one / temp, q(1,i+1), 1 )
  260. iconj = 1
  261. else
  262. iconj = 0
  263. end if
  264. end if
  265. 10 continue
  266. c
  267. call dgemv ('T', n, n, one, q, ldq, bounds, 1, zero, workl, 1)
  268. c
  269. if (msglvl .gt. 1) then
  270. call dvout (logfil, n, workl, ndigit,
  271. & '_neigh: Last row of the eigenvector matrix for H')
  272. end if
  273. c
  274. c %----------------------------%
  275. c | Compute the Ritz estimates |
  276. c %----------------------------%
  277. c
  278. iconj = 0
  279. do 20 i = 1, n
  280. if ( abs( ritzi(i) ) .le. zero ) then
  281. c
  282. c %----------------------%
  283. c | Real eigenvalue case |
  284. c %----------------------%
  285. c
  286. bounds(i) = rnorm * abs( workl(i) )
  287. else
  288. c
  289. c %-------------------------------------------%
  290. c | Complex conjugate pair case. Note that |
  291. c | since the real and imaginary part of |
  292. c | the eigenvector are stored in consecutive |
  293. c | columns, we need to take the magnitude |
  294. c | of the last components of the two vectors |
  295. c %-------------------------------------------%
  296. c
  297. if (iconj .eq. 0) then
  298. bounds(i) = rnorm * dlapy2( workl(i), workl(i+1) )
  299. bounds(i+1) = bounds(i)
  300. iconj = 1
  301. else
  302. iconj = 0
  303. end if
  304. end if
  305. 20 continue
  306. c
  307. if (msglvl .gt. 2) then
  308. call dvout (logfil, n, ritzr, ndigit,
  309. & '_neigh: Real part of the eigenvalues of H')
  310. call dvout (logfil, n, ritzi, ndigit,
  311. & '_neigh: Imaginary part of the eigenvalues of H')
  312. call dvout (logfil, n, bounds, ndigit,
  313. & '_neigh: Ritz estimates for the eigenvalues of H')
  314. end if
  315. c
  316. * call arscnd (t1)
  317. tneigh = tneigh + (t1 - t0)
  318. c
  319. 9000 continue
  320. return
  321. c
  322. c %---------------%
  323. c | End of dneigh |
  324. c %---------------%
  325. c
  326. end
  327.  
  328.  

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