6
c Given the Arnoldi factorization
8
c A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T,
10
c apply NP implicit shifts resulting in
12
c A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q
14
c where Q is an orthogonal matrix which is the product of rotations
15
c and reflections resulting from the NP bulge change sweeps.
16
c The updated Arnoldi factorization becomes:
18
c A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T.
22
c ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ,
27
c Problem size, i.e. size of matrix A.
29
c KEV Integer. (INPUT/OUTPUT)
30
c KEV+NP is the size of the input matrix H.
31
c KEV is the size of the updated matrix HNEW.
34
c Number of implicit shifts to be applied.
36
c SHIFT Complex*16 array of length NP. (INPUT)
37
c The shifts to be applied.
39
c V Complex*16 N by (KEV+NP) array. (INPUT/OUTPUT)
40
c On INPUT, V contains the current KEV+NP Arnoldi vectors.
41
c On OUTPUT, V contains the updated KEV Arnoldi vectors
42
c in the first KEV columns of V.
44
c LDV Integer. (INPUT)
45
c Leading dimension of V exactly as declared in the calling
48
c H Complex*16 (KEV+NP) by (KEV+NP) array. (INPUT/OUTPUT)
49
c On INPUT, H contains the current KEV+NP by KEV+NP upper
50
c Hessenberg matrix of the Arnoldi factorization.
51
c On OUTPUT, H contains the updated KEV by KEV upper Hessenberg
52
c matrix in the KEV leading submatrix.
54
c LDH Integer. (INPUT)
55
c Leading dimension of H exactly as declared in the calling
58
c RESID Complex*16 array of length N. (INPUT/OUTPUT)
59
c On INPUT, RESID contains the the residual vector r_{k+p}.
60
c On OUTPUT, RESID is the update residual vector rnew_{k}
61
c in the first KEV locations.
63
c Q Complex*16 KEV+NP by KEV+NP work array. (WORKSPACE)
64
c Work array used to accumulate the rotations and reflections
65
c during the bulge chase sweep.
67
c LDQ Integer. (INPUT)
68
c Leading dimension of Q exactly as declared in the calling
71
c WORKL Complex*16 work array of length (KEV+NP). (WORKSPACE)
72
c Private (replicated) array on each PE or array allocated on
75
c WORKD Complex*16 work array of length 2*N. (WORKSPACE)
76
c Distributed array used in the application of the accumulated
77
c orthogonal matrix Q.
81
c-----------------------------------------------------------------------
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),
94
c ivout ARPACK utility routine that prints integers.
95
c second ARPACK utility routine for timing.
96
c zmout ARPACK utility routine that prints matrices
97
c zvout ARPACK utility routine that prints vectors.
98
c zlacpy LAPACK matrix copy routine.
99
c zlanhs LAPACK routine that computes various norms of a matrix.
100
c zlartg LAPACK Givens rotation construction routine.
101
c zlaset LAPACK matrix initialization routine.
102
c dlabad LAPACK routine for defining the underflow and overflow
104
c dlamch LAPACK routine that determines machine constants.
105
c dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully.
106
c zgemv Level 2 BLAS routine for matrix vector multiplication.
107
c zaxpy Level 1 BLAS that computes a vector triad.
108
c zcopy Level 1 BLAS that copies one vector to another.
109
c zscal Level 1 BLAS that scales a vector.
112
c Danny Sorensen Phuong Vu
113
c Richard Lehoucq CRPC / Rice University
114
c Dept. of Computational & Houston, Texas
115
c Applied Mathematics
119
c\SCCS Information: @(#)
120
c FILE: napps.F SID: 2.3 DATE OF SID: 3/28/97 RELEASE: 2
123
c 1. In this version, each shift is applied to all the sublocks of
124
c the Hessenberg matrix H and not just to the submatrix that it
125
c comes from. Deflation as in LAPACK routine zlahqr (QR algorithm
126
c for upper Hessenberg matrices ) is used.
127
c Upon output, the subdiagonals of H are enforced to be non-negative
132
c-----------------------------------------------------------------------
135
& ( n, kev, np, shift, v, ldv, h, ldh, resid, q, ldq,
138
c %----------------------------------------------------%
139
c | Include files for debugging and timing information |
140
c %----------------------------------------------------%
145
c %------------------%
146
c | Scalar Arguments |
147
c %------------------%
149
integer kev, ldh, ldq, ldv, n, np
151
c %-----------------%
152
c | Array Arguments |
153
c %-----------------%
156
& h(ldh,kev+np), resid(n), shift(np),
157
& v(ldv,kev+np), q(ldq,kev+np), workd(2*n), workl(kev+np)
167
parameter (one = (1.0D+0, 0.0D+0), zero = (0.0D+0, 0.0D+0),
170
c %------------------------%
171
c | Local Scalars & Arrays |
172
c %------------------------%
174
integer i, iend, istart, j, jj, kplusp, msglvl
177
& cdum, f, g, h11, h21, r, s, sigma, t
179
& c, ovfl, smlnum, ulp, unfl, tst1
180
save first, ovfl, smlnum, ulp, unfl
182
c %----------------------%
183
c | External Subroutines |
184
c %----------------------%
186
external zaxpy, zcopy, zgemv, zscal, zlacpy, zlartg,
187
& zvout, zlaset, dlabad, zmout, second, ivout
189
c %--------------------%
190
c | External Functions |
191
c %--------------------%
194
& zlanhs, dlamch, dlapy2
195
external zlanhs, dlamch, dlapy2
197
c %----------------------%
198
c | Intrinsics Functions |
199
c %----------------------%
201
intrinsic abs, dimag, conjg, dcmplx, max, min, dble
203
c %---------------------%
204
c | Statement Functions |
205
c %---------------------%
209
zabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
215
data first / .true. /
217
c %-----------------------%
218
c | Executable Statements |
219
c %-----------------------%
223
c %-----------------------------------------------%
224
c | Set machine-dependent constants for the |
225
c | stopping criterion. If norm(H) <= sqrt(OVFL), |
226
c | overflow should not occur. |
227
c | REFERENCE: LAPACK subroutine zlahqr |
228
c %-----------------------------------------------%
230
unfl = dlamch( 'safe minimum' )
231
ovfl = dble(one / unfl)
232
call dlabad( unfl, ovfl )
233
ulp = dlamch( 'precision' )
234
smlnum = unfl*( n / ulp )
238
c %-------------------------------%
239
c | Initialize timing statistics |
240
c | & message level for debugging |
241
c %-------------------------------%
248
c %--------------------------------------------%
249
c | Initialize Q to the identity to accumulate |
250
c | the rotations and reflections |
251
c %--------------------------------------------%
253
call zlaset ('All', kplusp, kplusp, zero, one, q, ldq)
255
c %----------------------------------------------%
256
c | Quick return if there are no shifts to apply |
257
c %----------------------------------------------%
259
if (np .eq. 0) go to 9000
261
c %----------------------------------------------%
262
c | Chase the bulge with the application of each |
263
c | implicit shift. Each shift is applied to the |
264
c | whole matrix including each block. |
265
c %----------------------------------------------%
270
if (msglvl .gt. 2 ) then
271
call ivout (logfil, 1, jj, ndigit,
272
& '_napps: shift number.')
273
call zvout (logfil, 1, sigma, ndigit,
274
& '_napps: Value of the shift ')
280
do 30 i = istart, kplusp-1
282
c %----------------------------------------%
283
c | Check for splitting and deflation. Use |
284
c | a standard test as in the QR algorithm |
285
c | REFERENCE: LAPACK subroutine zlahqr |
286
c %----------------------------------------%
288
tst1 = zabs1( h( i, i ) ) + zabs1( h( i+1, i+1 ) )
290
& tst1 = zlanhs( '1', kplusp-jj+1, h, ldh, workl )
291
if ( abs(dble(h(i+1,i)))
292
& .le. max(ulp*tst1, smlnum) ) then
293
if (msglvl .gt. 0) then
294
call ivout (logfil, 1, i, ndigit,
295
& '_napps: matrix splitting at row/column no.')
296
call ivout (logfil, 1, jj, ndigit,
297
& '_napps: matrix splitting with shift number.')
298
call zvout (logfil, 1, h(i+1,i), ndigit,
299
& '_napps: off diagonal element.')
309
if (msglvl .gt. 2) then
310
call ivout (logfil, 1, istart, ndigit,
311
& '_napps: Start of current block ')
312
call ivout (logfil, 1, iend, ndigit,
313
& '_napps: End of current block ')
316
c %------------------------------------------------%
317
c | No reason to apply a shift to block of order 1 |
318
c | or if the current block starts after the point |
319
c | of compression since we'll discard this stuff |
320
c %------------------------------------------------%
322
if ( istart .eq. iend .or. istart .gt. kev) go to 100
324
h11 = h(istart,istart)
325
h21 = h(istart+1,istart)
329
do 80 i = istart, iend-1
331
c %------------------------------------------------------%
332
c | Construct the plane rotation G to zero out the bulge |
333
c %------------------------------------------------------%
335
call zlartg (f, g, c, s, r)
336
if (i .gt. istart) then
341
c %---------------------------------------------%
342
c | Apply rotation to the left of H; H <- G'*H |
343
c %---------------------------------------------%
346
t = c*h(i,j) + s*h(i+1,j)
347
h(i+1,j) = -conjg(s)*h(i,j) + c*h(i+1,j)
351
c %---------------------------------------------%
352
c | Apply rotation to the right of H; H <- H*G |
353
c %---------------------------------------------%
355
do 60 j = 1, min(i+2,iend)
356
t = c*h(j,i) + conjg(s)*h(j,i+1)
357
h(j,i+1) = -s*h(j,i) + c*h(j,i+1)
361
c %-----------------------------------------------------%
362
c | Accumulate the rotation in the matrix Q; Q <- Q*G' |
363
c %-----------------------------------------------------%
365
do 70 j = 1, min(i+jj, kplusp)
366
t = c*q(j,i) + conjg(s)*q(j,i+1)
367
q(j,i+1) = - s*q(j,i) + c*q(j,i+1)
371
c %---------------------------%
372
c | Prepare for next rotation |
373
c %---------------------------%
375
if (i .lt. iend-1) then
381
c %-------------------------------%
382
c | Finished applying the shift. |
383
c %-------------------------------%
387
c %---------------------------------------------------------%
388
c | Apply the same shift to the next block if there is any. |
389
c %---------------------------------------------------------%
392
if (iend .lt. kplusp) go to 20
394
c %---------------------------------------------%
395
c | Loop back to the top to get the next shift. |
396
c %---------------------------------------------%
400
c %---------------------------------------------------%
401
c | Perform a similarity transformation that makes |
402
c | sure that the compressed H will have non-negative |
403
c | real subdiagonal elements. |
404
c %---------------------------------------------------%
407
if ( dble( h(j+1,j) ) .lt. rzero .or.
408
& dimag( h(j+1,j) ) .ne. rzero ) then
409
t = h(j+1,j) / dlapy2(dble(h(j+1,j)),dimag(h(j+1,j)))
410
call zscal( kplusp-j+1, conjg(t), h(j+1,j), ldh )
411
call zscal( min(j+2, kplusp), t, h(1,j+1), 1 )
412
call zscal( min(j+np+1,kplusp), t, q(1,j+1), 1 )
413
h(j+1,j) = dcmplx( dble( h(j+1,j) ), rzero )
419
c %--------------------------------------------%
420
c | Final check for splitting and deflation. |
421
c | Use a standard test as in the QR algorithm |
422
c | REFERENCE: LAPACK subroutine zlahqr. |
423
c | Note: Since the subdiagonals of the |
424
c | compressed H are nonnegative real numbers, |
425
c | we take advantage of this. |
426
c %--------------------------------------------%
428
tst1 = zabs1( h( i, i ) ) + zabs1( h( i+1, i+1 ) )
429
if( tst1 .eq. rzero )
430
& tst1 = zlanhs( '1', kev, h, ldh, workl )
431
if( dble( h( i+1,i ) ) .le. max( ulp*tst1, smlnum ) )
435
c %-------------------------------------------------%
436
c | Compute the (kev+1)-st column of (V*Q) and |
437
c | temporarily store the result in WORKD(N+1:2*N). |
438
c | This is needed in the residual update since we |
439
c | cannot GUARANTEE that the corresponding entry |
440
c | of H would be zero as in exact arithmetic. |
441
c %-------------------------------------------------%
443
if ( dble( h(kev+1,kev) ) .gt. rzero )
444
& call zgemv ('N', n, kplusp, one, v, ldv, q(1,kev+1), 1, zero,
447
c %----------------------------------------------------------%
448
c | Compute column 1 to kev of (V*Q) in backward order |
449
c | taking advantage of the upper Hessenberg structure of Q. |
450
c %----------------------------------------------------------%
453
call zgemv ('N', n, kplusp-i+1, one, v, ldv,
454
& q(1,kev-i+1), 1, zero, workd, 1)
455
call zcopy (n, workd, 1, v(1,kplusp-i+1), 1)
458
c %-------------------------------------------------%
459
c | Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). |
460
c %-------------------------------------------------%
462
call zlacpy ('A', n, kev, v(1,kplusp-kev+1), ldv, v, ldv)
464
c %--------------------------------------------------------------%
465
c | Copy the (kev+1)-st column of (V*Q) in the appropriate place |
466
c %--------------------------------------------------------------%
468
if ( dble( h(kev+1,kev) ) .gt. rzero )
469
& call zcopy (n, workd(n+1), 1, v(1,kev+1), 1)
471
c %-------------------------------------%
472
c | Update the residual vector: |
473
c | r <- sigmak*r + betak*v(:,kev+1) |
475
c | sigmak = (e_{kev+p}'*Q)*e_{kev} |
476
c | betak = e_{kev+1}'*H*e_{kev} |
477
c %-------------------------------------%
479
call zscal (n, q(kplusp,kev), resid, 1)
480
if ( dble( h(kev+1,kev) ) .gt. rzero )
481
& call zaxpy (n, h(kev+1,kev), v(1,kev+1), 1, resid, 1)
483
if (msglvl .gt. 1) then
484
call zvout (logfil, 1, q(kplusp,kev), ndigit,
485
& '_napps: sigmak = (e_{kev+p}^T*Q)*e_{kev}')
486
call zvout (logfil, 1, h(kev+1,kev), ndigit,
487
& '_napps: betak = e_{kev+1}^T*H*e_{kev}')
488
call ivout (logfil, 1, kev, ndigit,
489
& '_napps: Order of the final Hessenberg matrix ')
490
if (msglvl .gt. 2) then
491
call zmout (logfil, kev, kev, h, ldh, ndigit,
492
& '_napps: updated Hessenberg matrix H for next iteration')
499
tcapps = tcapps + (t1 - t0)
added
removed
routines/arpack/znapps.f
