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