Télécharger dneigh.eso

Retour à la liste

Numérotation des lignes :

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

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