Mercurial > octave

comparison libcruft/arpack/src/zgetv0.f @ 12274:9f5d2ef078e8 release-3-4-x

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