Télécharger dgetv0.eso

Retour à la liste

Numérotation des lignes :

  1. C DGETV0 SOURCE BP208322 15/12/17 21:15:04 8750
  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. character bmat*1
  135. logical initv
  136. integer ido, ierr, itry, j, ldv, n
  137. REAL*8
  138. & rnorm
  139. c
  140. c %-----------------%
  141. c | Array Arguments |
  142. c %-----------------%
  143. c
  144. integer ipntr(3)
  145. REAL*8
  146. & resid(n), v(ldv,j), workd(2*n)
  147. c
  148. c %------------%
  149. c | Parameters |
  150. c %------------%
  151. c
  152. REAL*8
  153. & one, zero
  154. parameter (one = 1.0D+0, zero = 0.0D+0)
  155. c
  156. c %------------------------%
  157. c | Local Scalars & Arrays |
  158. c %------------------------%
  159. c
  160. logical first, inits, orth
  161. integer idist, iseed(4), iter, msglvl, jj
  162. REAL*8
  163. & rnorm0
  164. save first, iseed, inits, iter, msglvl, orth, rnorm0
  165. c
  166. c %----------------------%
  167. c | External Subroutines |
  168. c %----------------------%
  169. c
  170. c external dlarnv, dvout, dcopy, dgemv, arscnd
  171. c
  172. c %--------------------%
  173. c | External Functions |
  174. c %--------------------%
  175. c
  176. REAL*8
  177. external ddot, dnrm2
  178. c
  179. c %---------------------%
  180. **c | Intrinsic Functions |
  181. **c %---------------------%
  182. **c
  183. ** intrinsic abs, sqrt
  184. **c
  185. **c %-----------------%
  186. **c | Data Statements |
  187. **c %-----------------%
  188. **c
  189. data inits /.true./
  190. **c
  191. **c %-----------------------%
  192. **c | Executable Statements |
  193. c %-----------------------%
  194. c
  195. c
  196. c %-----------------------------------%
  197. c | Initialize the seed of the LAPACK |
  198. c | random number generator |
  199. c %-----------------------------------%
  200. c
  201. if (inits) then
  202. iseed(1) = 1
  203. iseed(2) = 3
  204. iseed(3) = 5
  205. iseed(4) = 7
  206. inits = .false.
  207. end if
  208. c
  209. if (ido .eq. 0) then
  210. c
  211. c %-------------------------------%
  212. c | Initialize timing statistics |
  213. c | & message level for debugging |
  214. c %-------------------------------%
  215. c
  216. * call arscnd (t0)
  217. msglvl = mgetv0
  218. c
  219. ierr = 0
  220. iter = 0
  221. first = .FALSE.
  222. orth = .FALSE.
  223. c
  224. c %-----------------------------------------------------%
  225. c | Possibly generate a random starting vector in RESID |
  226. c | Use a LAPACK random number generator used by the |
  227. c | matrix generation routines. |
  228. c | idist = 1: uniform (0,1) distribution; |
  229. c | idist = 2: uniform (-1,1) distribution; |
  230. c | idist = 3: normal (0,1) distribution; |
  231. c %-----------------------------------------------------%
  232. c
  233. if (.not.initv) then
  234. idist = 2
  235. call dlarnv (idist, iseed, n, resid)
  236. end if
  237. c
  238. c %----------------------------------------------------------%
  239. c | Force the starting vector into the range of OP to handle |
  240. c | the generalized problem when B is possibly (singular). |
  241. c %----------------------------------------------------------%
  242. c
  243. * call arscnd (t2)
  244. if (bmat .eq. 'G') then
  245. nopx = nopx + 1
  246. ipntr(1) = 1
  247. ipntr(2) = n + 1
  248. call dcopy (n, resid, 1, workd, 1)
  249. ido = -1
  250. go to 9000
  251. end if
  252. end if
  253. c
  254. c %-----------------------------------------%
  255. c | Back from computing OP*(initial-vector) |
  256. c %-----------------------------------------%
  257. c
  258. if (first) go to 20
  259. c
  260. c %-----------------------------------------------%
  261. c | Back from computing B*(orthogonalized-vector) |
  262. c %-----------------------------------------------%
  263. c
  264. if (orth) go to 40
  265. c
  266. if (bmat .eq. 'G') then
  267. * call arscnd (t3)
  268. tmvopx = tmvopx + (t3 - t2)
  269. end if
  270. c
  271. c %------------------------------------------------------%
  272. c | Starting vector is now in the range of OP; r = OP*r; |
  273. c | Compute B-norm of starting vector. |
  274. c %------------------------------------------------------%
  275. c
  276. * call arscnd (t2)
  277. first = .TRUE.
  278. if (bmat .eq. 'G') then
  279. nbx = nbx + 1
  280. call dcopy (n, workd(n+1), 1, resid, 1)
  281. ipntr(1) = n + 1
  282. ipntr(2) = 1
  283. ido = 2
  284. go to 9000
  285. else if (bmat .eq. 'I') then
  286. call dcopy (n, resid, 1, workd, 1)
  287. end if
  288. c
  289. 20 continue
  290. c
  291. if (bmat .eq. 'G') then
  292. * call arscnd (t3)
  293. tmvbx = tmvbx + (t3 - t2)
  294. end if
  295. c
  296. first = .FALSE.
  297. if (bmat .eq. 'G') then
  298. rnorm0 = ddot (n, resid, 1, workd, 1)
  299. rnorm0 = sqrt(abs(rnorm0))
  300. else if (bmat .eq. 'I') then
  301. rnorm0 = dnrm2(n, resid, 1)
  302. end if
  303. rnorm = rnorm0
  304. c
  305. c %---------------------------------------------%
  306. c | Exit if this is the very first Arnoldi step |
  307. c %---------------------------------------------%
  308. c
  309. if (j .eq. 1) go to 50
  310. c
  311. c %----------------------------------------------------------------
  312. c | Otherwise need to B-orthogonalize the starting vector against |
  313. c | the current Arnoldi basis using Gram-Schmidt with iter. ref. |
  314. c | This is the case where an invariant subspace is encountered |
  315. c | in the middle of the Arnoldi factorization. |
  316. c | |
  317. c | s = V^{T}*B*r; r = r - V*s; |
  318. c | |
  319. c | Stopping criteria used for iter. ref. is discussed in |
  320. c | Parlett's book, page 107 and in Gragg & Reichel TOMS paper. |
  321. c %---------------------------------------------------------------%
  322. c
  323. orth = .TRUE.
  324. 30 continue
  325.  
  326. call dgemv ('T', n, j-1, one, v, ldv, workd, 1,
  327. & zero, workd(n+1), 1)
  328. call dgemv ('N', n, j-1, (-one), v, ldv, workd(n+1), 1,
  329. & one, resid, 1)
  330. c
  331. c %----------------------------------------------------------%
  332. c | Compute the B-norm of the orthogonalized starting vector |
  333. c %----------------------------------------------------------%
  334. c
  335. * call arscnd (t2)
  336. if (bmat .eq. 'G') then
  337. nbx = nbx + 1
  338. call dcopy (n, resid, 1, workd(n+1), 1)
  339. ipntr(1) = n + 1
  340. ipntr(2) = 1
  341. ido = 2
  342. go to 9000
  343. else if (bmat .eq. 'I') then
  344. call dcopy (n, resid, 1, workd, 1)
  345. end if
  346. c
  347. 40 continue
  348. c
  349. if (bmat .eq. 'G') then
  350. * call arscnd (t3)
  351. tmvbx = tmvbx + (t3 - t2)
  352. end if
  353. c
  354. if (bmat .eq. 'G') then
  355. rnorm = ddot (n, resid, 1, workd, 1)
  356. rnorm = sqrt(abs(rnorm))
  357. else if (bmat .eq. 'I') then
  358. rnorm = dnrm2(n, resid, 1)
  359. end if
  360. c
  361. c %--------------------------------------%
  362. c | Check for further orthogonalization. |
  363. c %--------------------------------------%
  364. c
  365. c if (msglvl .gt. 2) then
  366. c call dvout (logfil, 1, rnorm0, ndigit,
  367. c & '_getv0: re-orthonalization ; rnorm0 is')
  368. c call dvout (logfil, 1, rnorm, ndigit,
  369. c & '_getv0: re-orthonalization ; rnorm is')
  370. c end if
  371. c
  372. if (rnorm .gt. 0.717*rnorm0) go to 50
  373. c
  374. iter = iter + 1
  375. if (iter .le. 5) then
  376. c
  377. c %-----------------------------------%
  378. c | Perform iterative refinement step |
  379. c %-----------------------------------%
  380. c
  381. rnorm0 = rnorm
  382. go to 30
  383. else
  384. c
  385. c %------------------------------------%
  386. c | Iterative refinement step "failed" |
  387. c %------------------------------------%
  388. c
  389. do 45 jj = 1, n
  390. resid(jj) = zero
  391. 45 continue
  392. rnorm = zero
  393. ierr = -1
  394. end if
  395. c
  396. 50 continue
  397.  
  398. c if (msglvl .gt. 0) then
  399. c call dvout (logfil, 1, rnorm, ndigit,
  400. c & '_getv0: B-norm of initial - restarted starting vector')
  401. c end if
  402. c if (msglvl .gt. 3) then
  403. c call dvout (logfil, n, resid, ndigit,
  404. c & '_getv0: initial - restarted starting vector')
  405. c end if
  406. ido = 99
  407. c
  408. * call arscnd (t1)
  409. tgetv0 = tgetv0 + (t1 - t0)
  410. c
  411. 9000 continue
  412. return
  413. c
  414. c %---------------%
  415. c | End of dgetv0 |
  416. c %---------------%
  417. c
  418. end
  419.  
  420.  
  421.  

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