Télécharger dgetv0.eso

Retour à la liste

Numérotation des lignes :

  1. C DGETV0 SOURCE GF238795 18/02/01 21:15:13 9724
  2. c-----------------------------------------------------------------------
  3. c\BeginDoc
  4. c
  5. c\Name: dgetv0
  6. c
  7. c\Description:
  8. c Generate a random initial residual vector for the Arnoldi process.
  9. c Force the residual vector to be in the range of the operator OP.
  10. c
  11. c\Usage:
  12. c call dgetv0
  13. c ( IDO, BMAT, ITRY, INITV, N, J, V, LDV, RESID, RNORM,
  14. c IPNTR, WORKD, IERR )
  15. c
  16. c\Arguments
  17. c IDO Integer. (INPUT/OUTPUT)
  18. c Reverse communication flag. IDO must be zero on the first
  19. c call to dgetv0.
  20. c -------------------------------------------------------------
  21. c IDO = 0: first call to the reverse communication interface
  22. c IDO = -1: compute Y = OP * X where
  23. c IPNTR(1) is the pointer into WORKD for X,
  24. c IPNTR(2) is the pointer into WORKD for Y.
  25. c This is for the initialization phase to force the
  26. c starting vector into the range of OP.
  27. c IDO = 2: compute Y = B * X where
  28. c IPNTR(1) is the pointer into WORKD for X,
  29. c IPNTR(2) is the pointer into WORKD for Y.
  30. c IDO = 99: done
  31. c -------------------------------------------------------------
  32. c
  33. c BMAT Character*1. (INPUT)
  34. c BMAT specifies the type of the matrix B in the (generalized)
  35. c eigenvalue problem A*x = lambda*B*x.
  36. c B = 'I' -> standard eigenvalue problem A*x = lambda*x
  37. c B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
  38. c
  39. c ITRY Integer. (INPUT)
  40. c ITRY counts the number of times that dgetv0 is called.
  41. c It should be set to 1 on the initial call to dgetv0.
  42. c
  43. c INITV Logical variable. (INPUT)
  44. c .TRUE. => the initial residual vector is given in RESID.
  45. c .FALSE. => generate a random initial residual vector.
  46. c
  47. c N Integer. (INPUT)
  48. c Dimension of the problem.
  49. c
  50. c J Integer. (INPUT)
  51. c Index of the residual vector to be generated, with respect to
  52. c the Arnoldi process. J > 1 in case of a "restart".
  53. c
  54. c V REAL*8 N by J array. (INPUT)
  55. c The first J-1 columns of V contain the current Arnoldi basis
  56. c if this is a "restart".
  57. c
  58. c LDV Integer. (INPUT)
  59. c Leading dimension of V exactly as declared in the calling
  60. c program.
  61. c
  62. c RESID Double precision array of length N. (INPUT/OUTPUT)
  63. c Initial residual vector to be generated. If RESID is
  64. c provided, force RESID into the range of the operator OP.
  65. c
  66. c RNORM Double precision scalar. (OUTPUT)
  67. c B-norm of the generated residual.
  68. c
  69. c IPNTR Integer array of length 3. (OUTPUT)
  70. c
  71. c WORKD Double precision work array of length 2*N. (REVERSE COMMUNICATION).
  72. c On exit, WORK(1:N) = B*RESID to be used in SSAITR.
  73. c
  74. c IERR Integer. (OUTPUT)
  75. c = 0: Normal exit.
  76. c = -1: Cannot generate a nontrivial restarted residual vector
  77. c in the range of the operator OP.
  78. c
  79. c\EndDoc
  80. c
  81. c-----------------------------------------------------------------------
  82. c
  83. c\BeginLib
  84. c
  85. c\Local variables:
  86. c xxxxxx real
  87. c
  88. c\References:
  89. c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
  90. c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
  91. c pp 357-385.
  92. c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
  93. c Restarted Arnoldi Iteration", Rice University Technical Report
  94. c TR95-13, Department of Computational and Applied Mathematics.
  95. c
  96. c\Routines called:
  97. c arscnd ARPACK utility routine for timing.
  98. c dvout ARPACK utility routine for vector output.
  99. c dlarnv LAPACK routine for generating a random vector.
  100. c dgemv Level 2 BLAS routine for matrix vector multiplication.
  101. c dcopy Level 1 BLAS that copies one vector to another.
  102. c ddot Level 1 BLAS that computes the scalar product of two vectors.
  103. c dnrm2 Level 1 BLAS that computes the norm of a vector.
  104. c
  105. c\Author
  106. c Danny Sorensen Phuong Vu
  107. c Richard Lehoucq CRPC / Rice University
  108. c Dept. of Computational & Houston, Texas
  109. c Applied Mathematics
  110. c Rice University
  111. c Houston, Texas
  112. c
  113. c\SCCS Information: @(#)
  114. c FILE: getv0.F SID: 2.7 DATE OF SID: 04/07/99 RELEASE: 2
  115. c
  116. c\EndLib
  117. c
  118. c-----------------------------------------------------------------------
  119. c
  120. subroutine dgetv0
  121. & ( ido, bmat, itry, initv, n, j, v, ldv, resid, rnorm,
  122. & ipntr, workd, ierr )
  123. c
  124. c %----------------------------------------------------%
  125. c | Include files for debugging and timing information |
  126. -INC TARTRAK
  127. c %----------------------------------------------------%
  128. c
  129. c
  130. c %------------------%
  131. c | Scalar Arguments |
  132. c %------------------%
  133. c
  134. REAL*8 T0
  135. REAL*8 T1
  136. REAL*8 T2
  137. REAL*8 T3
  138. character bmat*1
  139. logical initv
  140. integer ido, ierr, itry, j, ldv, n
  141. REAL*8
  142. & rnorm
  143. c
  144. c %-----------------%
  145. c | Array Arguments |
  146. c %-----------------%
  147. c
  148. integer ipntr(3)
  149. REAL*8
  150. & resid(n), v(ldv,j), workd(2*n)
  151. c
  152. c %------------%
  153. c | Parameters |
  154. c %------------%
  155. c
  156. REAL*8
  157. & one, zero
  158. parameter (one = 1.0D+0, zero = 0.0D+0)
  159. c
  160. c %------------------------%
  161. c | Local Scalars & Arrays |
  162. c %------------------------%
  163. c
  164. logical first, inits, orth
  165. integer idist, iseed(4), iter, msglvl, jj
  166. REAL*8
  167. & rnorm0
  168. save first, iseed, inits, iter, msglvl, orth, rnorm0
  169. c
  170. c %----------------------%
  171. c | External Subroutines |
  172. c %----------------------%
  173. c
  174. c external dlarnv, dvout, dcopy, dgemv, arscnd
  175. c
  176. c %--------------------%
  177. c | External Functions |
  178. c %--------------------%
  179. c
  180. REAL*8
  181. external ddot, dnrm2
  182. c
  183. c %---------------------%
  184. **c | Intrinsic Functions |
  185. **c %---------------------%
  186. **c
  187. ** intrinsic abs, sqrt
  188. **c
  189. **c %-----------------%
  190. **c | Data Statements |
  191. **c %-----------------%
  192. **c
  193. data inits /.true./
  194. **c
  195. **c %-----------------------%
  196. **c | Executable Statements |
  197. c %-----------------------%
  198. T0 = 0.D0
  199. T1 = 0.D0
  200. T2 = 0.D0
  201. T3 = 0.D0
  202. c
  203. c
  204. c %-----------------------------------%
  205. c | Initialize the seed of the LAPACK |
  206. c | random number generator |
  207. c %-----------------------------------%
  208. c
  209. if (inits) then
  210. iseed(1) = 1
  211. iseed(2) = 3
  212. iseed(3) = 5
  213. iseed(4) = 7
  214. inits = .false.
  215. end if
  216. c
  217. if (ido .eq. 0) then
  218. c
  219. c %-------------------------------%
  220. c | Initialize timing statistics |
  221. c | & message level for debugging |
  222. c %-------------------------------%
  223. c
  224. * call arscnd (t0)
  225. msglvl = mgetv0
  226. c
  227. ierr = 0
  228. iter = 0
  229. first = .FALSE.
  230. orth = .FALSE.
  231. c
  232. c %-----------------------------------------------------%
  233. c | Possibly generate a random starting vector in RESID |
  234. c | Use a LAPACK random number generator used by the |
  235. c | matrix generation routines. |
  236. c | idist = 1: uniform (0,1) distribution; |
  237. c | idist = 2: uniform (-1,1) distribution; |
  238. c | idist = 3: normal (0,1) distribution; |
  239. c %-----------------------------------------------------%
  240. c
  241. if (.not.initv) then
  242. idist = 2
  243. call dlarnv (idist, iseed, n, resid)
  244. end if
  245. c
  246. c %----------------------------------------------------------%
  247. c | Force the starting vector into the range of OP to handle |
  248. c | the generalized problem when B is possibly (singular). |
  249. c %----------------------------------------------------------%
  250. c
  251. * call arscnd (t2)
  252. if (bmat .eq. 'G') then
  253. nopx = nopx + 1
  254. ipntr(1) = 1
  255. ipntr(2) = n + 1
  256. call dcopy (n, resid, 1, workd, 1)
  257. ido = -1
  258. go to 9000
  259. end if
  260. end if
  261. c
  262. c %-----------------------------------------%
  263. c | Back from computing OP*(initial-vector) |
  264. c %-----------------------------------------%
  265. c
  266. if (first) go to 20
  267. c
  268. c %-----------------------------------------------%
  269. c | Back from computing B*(orthogonalized-vector) |
  270. c %-----------------------------------------------%
  271. c
  272. if (orth) go to 40
  273. c
  274. if (bmat .eq. 'G') then
  275. * call arscnd (t3)
  276. tmvopx = tmvopx + (t3 - t2)
  277. end if
  278. c
  279. c %------------------------------------------------------%
  280. c | Starting vector is now in the range of OP; r = OP*r; |
  281. c | Compute B-norm of starting vector. |
  282. c %------------------------------------------------------%
  283. c
  284. * call arscnd (t2)
  285. first = .TRUE.
  286. if (bmat .eq. 'G') then
  287. nbx = nbx + 1
  288. call dcopy (n, workd(n+1), 1, resid, 1)
  289. ipntr(1) = n + 1
  290. ipntr(2) = 1
  291. ido = 2
  292. go to 9000
  293. else if (bmat .eq. 'I') then
  294. call dcopy (n, resid, 1, workd, 1)
  295. end if
  296. c
  297. 20 continue
  298. c
  299. if (bmat .eq. 'G') then
  300. * call arscnd (t3)
  301. tmvbx = tmvbx + (t3 - t2)
  302. end if
  303. c
  304. first = .FALSE.
  305. if (bmat .eq. 'G') then
  306. rnorm0 = ddot (n, resid, 1, workd, 1)
  307. rnorm0 = sqrt(abs(rnorm0))
  308. else if (bmat .eq. 'I') then
  309. rnorm0 = dnrm2(n, resid, 1)
  310. end if
  311. rnorm = rnorm0
  312. c
  313. c %---------------------------------------------%
  314. c | Exit if this is the very first Arnoldi step |
  315. c %---------------------------------------------%
  316. c
  317. if (j .eq. 1) go to 50
  318. c
  319. c %----------------------------------------------------------------
  320. c | Otherwise need to B-orthogonalize the starting vector against |
  321. c | the current Arnoldi basis using Gram-Schmidt with iter. ref. |
  322. c | This is the case where an invariant subspace is encountered |
  323. c | in the middle of the Arnoldi factorization. |
  324. c | |
  325. c | s = V^{T}*B*r; r = r - V*s; |
  326. c | |
  327. c | Stopping criteria used for iter. ref. is discussed in |
  328. c | Parlett's book, page 107 and in Gragg & Reichel TOMS paper. |
  329. c %---------------------------------------------------------------%
  330. c
  331. orth = .TRUE.
  332. 30 continue
  333.  
  334. call dgemv ('T', n, j-1, one, v, ldv, workd, 1,
  335. & zero, workd(n+1), 1)
  336. call dgemv ('N', n, j-1, (-one), v, ldv, workd(n+1), 1,
  337. & one, resid, 1)
  338. c
  339. c %----------------------------------------------------------%
  340. c | Compute the B-norm of the orthogonalized starting vector |
  341. c %----------------------------------------------------------%
  342. c
  343. * call arscnd (t2)
  344. if (bmat .eq. 'G') then
  345. nbx = nbx + 1
  346. call dcopy (n, resid, 1, workd(n+1), 1)
  347. ipntr(1) = n + 1
  348. ipntr(2) = 1
  349. ido = 2
  350. go to 9000
  351. else if (bmat .eq. 'I') then
  352. call dcopy (n, resid, 1, workd, 1)
  353. end if
  354. c
  355. 40 continue
  356. c
  357. if (bmat .eq. 'G') then
  358. * call arscnd (t3)
  359. tmvbx = tmvbx + (t3 - t2)
  360. end if
  361. c
  362. if (bmat .eq. 'G') then
  363. rnorm = ddot (n, resid, 1, workd, 1)
  364. rnorm = sqrt(abs(rnorm))
  365. else if (bmat .eq. 'I') then
  366. rnorm = dnrm2(n, resid, 1)
  367. end if
  368. c
  369. c %--------------------------------------%
  370. c | Check for further orthogonalization. |
  371. c %--------------------------------------%
  372. c
  373. c if (msglvl .gt. 2) then
  374. c call dvout (logfil, 1, rnorm0, ndigit,
  375. c & '_getv0: re-orthonalization ; rnorm0 is')
  376. c call dvout (logfil, 1, rnorm, ndigit,
  377. c & '_getv0: re-orthonalization ; rnorm is')
  378. c end if
  379. c
  380. if (rnorm .gt. 0.717*rnorm0) go to 50
  381. c
  382. iter = iter + 1
  383. if (iter .le. 5) then
  384. c
  385. c %-----------------------------------%
  386. c | Perform iterative refinement step |
  387. c %-----------------------------------%
  388. c
  389. rnorm0 = rnorm
  390. go to 30
  391. else
  392. c
  393. c %------------------------------------%
  394. c | Iterative refinement step "failed" |
  395. c %------------------------------------%
  396. c
  397. do 45 jj = 1, n
  398. resid(jj) = zero
  399. 45 continue
  400. rnorm = zero
  401. ierr = -1
  402. end if
  403. c
  404. 50 continue
  405.  
  406. c if (msglvl .gt. 0) then
  407. c call dvout (logfil, 1, rnorm, ndigit,
  408. c & '_getv0: B-norm of initial - restarted starting vector')
  409. c end if
  410. c if (msglvl .gt. 3) then
  411. c call dvout (logfil, n, resid, ndigit,
  412. c & '_getv0: initial - restarted starting vector')
  413. c end if
  414. ido = 99
  415. c
  416. * call arscnd (t1)
  417. tgetv0 = tgetv0 + (t1 - t0)
  418. c
  419. 9000 continue
  420. return
  421. c
  422. c %---------------%
  423. c | End of dgetv0 |
  424. c %---------------%
  425. c
  426. end
  427.  
  428.  
  429.  
  430.  

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