6
c This subroutine returns the converged approximations to eigenvalues
7
c of A*z = lambda*B*z and (optionally):
9
c (1) The corresponding approximate eigenvectors;
11
c (2) An orthonormal basis for the associated approximate
16
c There is negligible additional cost to obtain eigenvectors. An orthonormal
17
c basis is always computed. There is an additional storage cost of n*nev
18
c if both are requested (in this case a separate array Z must be supplied).
20
c The approximate eigenvalues and eigenvectors of A*z = lambda*B*z
21
c are derived from approximate eigenvalues and eigenvectors of
22
c of the linear operator OP prescribed by the MODE selection in the
23
c call to ZNAUPD . ZNAUPD must be called before this routine is called.
24
c These approximate eigenvalues and vectors are commonly called Ritz
25
c values and Ritz vectors respectively. They are referred to as such
26
c in the comments that follow. The computed orthonormal basis for the
27
c invariant subspace corresponding to these Ritz values is referred to as a
30
c The definition of OP as well as other terms and the relation of computed
31
c Ritz values and vectors of OP with respect to the given problem
32
c A*z = lambda*B*z may be found in the header of ZNAUPD . For a brief
33
c description, see definitions of IPARAM(7), MODE and WHICH in the
34
c documentation of ZNAUPD .
38
c ( RVEC, HOWMNY, SELECT, D, Z, LDZ, SIGMA, WORKEV, BMAT,
39
c N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR, WORKD,
40
c WORKL, LWORKL, RWORK, INFO )
43
c RVEC LOGICAL (INPUT)
44
c Specifies whether a basis for the invariant subspace corresponding
45
c to the converged Ritz value approximations for the eigenproblem
46
c A*z = lambda*B*z is computed.
48
c RVEC = .FALSE. Compute Ritz values only.
50
c RVEC = .TRUE. Compute Ritz vectors or Schur vectors.
53
c HOWMNY Character*1 (INPUT)
54
c Specifies the form of the basis for the invariant subspace
55
c corresponding to the converged Ritz values that is to be computed.
57
c = 'A': Compute NEV Ritz vectors;
58
c = 'P': Compute NEV Schur vectors;
59
c = 'S': compute some of the Ritz vectors, specified
60
c by the logical array SELECT.
62
c SELECT Logical array of dimension NCV. (INPUT)
63
c If HOWMNY = 'S', SELECT specifies the Ritz vectors to be
64
c computed. To select the Ritz vector corresponding to a
65
c Ritz value D(j), SELECT(j) must be set to .TRUE..
66
c If HOWMNY = 'A' or 'P', SELECT need not be initialized
67
c but it is used as internal workspace.
69
c D Complex*16 array of dimension NEV+1. (OUTPUT)
70
c On exit, D contains the Ritz approximations
71
c to the eigenvalues lambda for A*z = lambda*B*z.
73
c Z Complex*16 N by NEV array (OUTPUT)
74
c On exit, if RVEC = .TRUE. and HOWMNY = 'A', then the columns of
75
c Z represents approximate eigenvectors (Ritz vectors) corresponding
76
c to the NCONV=IPARAM(5) Ritz values for eigensystem
79
c If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED.
81
c NOTE: If if RVEC = .TRUE. and a Schur basis is not required,
82
c the array Z may be set equal to first NEV+1 columns of the Arnoldi
83
c basis array V computed by ZNAUPD . In this case the Arnoldi basis
84
c will be destroyed and overwritten with the eigenvector basis.
86
c LDZ Integer. (INPUT)
87
c The leading dimension of the array Z. If Ritz vectors are
88
c desired, then LDZ .ge. max( 1, N ) is required.
89
c In any case, LDZ .ge. 1 is required.
91
c SIGMA Complex*16 (INPUT)
92
c If IPARAM(7) = 3 then SIGMA represents the shift.
93
c Not referenced if IPARAM(7) = 1 or 2.
95
c WORKEV Complex*16 work array of dimension 2*NCV. (WORKSPACE)
97
c **** The remaining arguments MUST be the same as for the ****
98
c **** call to ZNAUPD that was just completed. ****
100
c NOTE: The remaining arguments
102
c BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
103
c WORKD, WORKL, LWORKL, RWORK, INFO
105
c must be passed directly to ZNEUPD following the last call
106
c to ZNAUPD . These arguments MUST NOT BE MODIFIED between
107
c the the last call to ZNAUPD and the call to ZNEUPD .
109
c Three of these parameters (V, WORKL and INFO) are also output parameters:
111
c V Complex*16 N by NCV array. (INPUT/OUTPUT)
113
c Upon INPUT: the NCV columns of V contain the Arnoldi basis
114
c vectors for OP as constructed by ZNAUPD .
116
c Upon OUTPUT: If RVEC = .TRUE. the first NCONV=IPARAM(5) columns
117
c contain approximate Schur vectors that span the
118
c desired invariant subspace.
120
c NOTE: If the array Z has been set equal to first NEV+1 columns
121
c of the array V and RVEC=.TRUE. and HOWMNY= 'A', then the
122
c Arnoldi basis held by V has been overwritten by the desired
123
c Ritz vectors. If a separate array Z has been passed then
124
c the first NCONV=IPARAM(5) columns of V will contain approximate
125
c Schur vectors that span the desired invariant subspace.
127
c WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE)
128
c WORKL(1:ncv*ncv+2*ncv) contains information obtained in
129
c znaupd . They are not changed by zneupd .
130
c WORKL(ncv*ncv+2*ncv+1:3*ncv*ncv+4*ncv) holds the
131
c untransformed Ritz values, the untransformed error estimates of
132
c the Ritz values, the upper triangular matrix for H, and the
133
c associated matrix representation of the invariant subspace for H.
135
c Note: IPNTR(9:13) contains the pointer into WORKL for addresses
136
c of the above information computed by zneupd .
137
c -------------------------------------------------------------
138
c IPNTR(9): pointer to the NCV RITZ values of the
140
c IPNTR(10): Not used
141
c IPNTR(11): pointer to the NCV corresponding error estimates.
142
c IPNTR(12): pointer to the NCV by NCV upper triangular
143
c Schur matrix for H.
144
c IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
145
c of the upper Hessenberg matrix H. Only referenced by
146
c zneupd if RVEC = .TRUE. See Remark 2 below.
147
c -------------------------------------------------------------
149
c INFO Integer. (OUTPUT)
150
c Error flag on output.
153
c = 1: The Schur form computed by LAPACK routine csheqr
154
c could not be reordered by LAPACK routine ztrsen .
155
c Re-enter subroutine zneupd with IPARAM(5)=NCV and
156
c increase the size of the array D to have
157
c dimension at least dimension NCV and allocate at least NCV
158
c columns for Z. NOTE: Not necessary if Z and V share
159
c the same space. Please notify the authors if this error
162
c = -1: N must be positive.
163
c = -2: NEV must be positive.
164
c = -3: NCV-NEV >= 2 and less than or equal to N.
165
c = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
166
c = -6: BMAT must be one of 'I' or 'G'.
167
c = -7: Length of private work WORKL array is not sufficient.
168
c = -8: Error return from LAPACK eigenvalue calculation.
169
c This should never happened.
170
c = -9: Error return from calculation of eigenvectors.
171
c Informational error from LAPACK routine ztrevc .
172
c = -10: IPARAM(7) must be 1,2,3
173
c = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatible.
174
c = -12: HOWMNY = 'S' not yet implemented
175
c = -13: HOWMNY must be one of 'A' or 'P' if RVEC = .true.
176
c = -14: ZNAUPD did not find any eigenvalues to sufficient
178
c = -15: ZNEUPD got a different count of the number of converged
179
c Ritz values than ZNAUPD got. This indicates the user
180
c probably made an error in passing data from ZNAUPD to
181
c ZNEUPD or that the data was modified before entering
187
c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
188
c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
190
c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
191
c Restarted Arnoldi Iteration", Rice University Technical Report
192
c TR95-13, Department of Computational and Applied Mathematics.
193
c 3. B. Nour-Omid, B. N. Parlett, T. Ericsson and P. S. Jensen,
194
c "How to Implement the Spectral Transformation", Math Comp.,
195
c Vol. 48, No. 178, April, 1987 pp. 664-673.
198
c ivout ARPACK utility routine that prints integers.
199
c zmout ARPACK utility routine that prints matrices
200
c zvout ARPACK utility routine that prints vectors.
201
c zgeqr2 LAPACK routine that computes the QR factorization of
203
c zlacpy LAPACK matrix copy routine.
204
c zlahqr LAPACK routine that computes the Schur form of a
205
c upper Hessenberg matrix.
206
c zlaset LAPACK matrix initialization routine.
207
c ztrevc LAPACK routine to compute the eigenvectors of a matrix
208
c in upper triangular form.
209
c ztrsen LAPACK routine that re-orders the Schur form.
210
c zunm2r LAPACK routine that applies an orthogonal matrix in
212
c dlamch LAPACK routine that determines machine constants.
213
c ztrmm Level 3 BLAS matrix times an upper triangular matrix.
214
c zgeru Level 2 BLAS rank one update to a matrix.
215
c zcopy Level 1 BLAS that copies one vector to another .
216
c zscal Level 1 BLAS that scales a vector.
217
c zdscal Level 1 BLAS that scales a complex vector by a real number.
218
c dznrm2 Level 1 BLAS that computes the norm of a complex vector.
222
c 1. Currently only HOWMNY = 'A' and 'P' are implemented.
224
c 2. Schur vectors are an orthogonal representation for the basis of
225
c Ritz vectors. Thus, their numerical properties are often superior.
226
c If RVEC = .true. then the relationship
227
c A * V(:,1:IPARAM(5)) = V(:,1:IPARAM(5)) * T, and
228
c transpose( V(:,1:IPARAM(5)) ) * V(:,1:IPARAM(5)) = I
229
c are approximately satisfied.
230
c Here T is the leading submatrix of order IPARAM(5) of the
231
c upper triangular matrix stored workl(ipntr(12)).
234
c Danny Sorensen Phuong Vu
235
c Richard Lehoucq CRPC / Rice University
236
c Chao Yang Houston, Texas
237
c Dept. of Computational &
238
c Applied Mathematics
242
c\SCCS Information: @(#)
243
c FILE: neupd.F SID: 2.7 DATE OF SID: 09/20/00 RELEASE: 2
247
c-----------------------------------------------------------------------
248
subroutine zneupd (rvec , howmny, select, d ,
249
& z , ldz , sigma , workev,
250
& bmat , n , which , nev ,
251
& tol , resid , ncv , v ,
252
& ldv , iparam, ipntr , workd ,
253
& workl, lworkl, rwork , info )
255
c %----------------------------------------------------%
256
c | Include files for debugging and timing information |
257
c %----------------------------------------------------%
262
c %------------------%
263
c | Scalar Arguments |
264
c %------------------%
266
character bmat, howmny, which*2
268
integer info, ldz, ldv, lworkl, n, ncv, nev
274
c %-----------------%
275
c | Array Arguments |
276
c %-----------------%
278
integer iparam(11), ipntr(14)
283
& d(nev) , resid(n) , v(ldv,ncv),
285
& workd(3*n) , workl(lworkl), workev(2*ncv)
293
parameter (one = (1.0D+0, 0.0D+0) , zero = (0.0D+0, 0.0D+0) )
300
integer bounds, ierr , ih , ihbds, iheig , nconv ,
301
& invsub, iuptri, iwev , j , ldh , ldq ,
302
& mode , msglvl, ritz , wr , k , irz ,
303
& ibd , outncv, iq , np , numcnv, jj ,
308
& conds, sep, rtemp, eps23
311
c %----------------------%
312
c | External Subroutines |
313
c %----------------------%
315
external zcopy , zgeru , zgeqr2 , zlacpy , zmout ,
316
& zunm2r , ztrmm , zvout , ivout,
319
c %--------------------%
320
c | External Functions |
321
c %--------------------%
324
& dznrm2 , dlamch , dlapy2
325
external dznrm2 , dlamch , dlapy2
331
c %-----------------------%
332
c | Executable Statements |
333
c %-----------------------%
335
c %------------------------%
336
c | Set default parameters |
337
c %------------------------%
345
c %---------------------------------%
346
c | Get machine dependent constant. |
347
c %---------------------------------%
349
eps23 = dlamch ('Epsilon-Machine')
350
eps23 = eps23**(2.0D+0 / 3.0D+0 )
352
c %-------------------------------%
354
c | Check for incompatible input |
355
c %-------------------------------%
359
if (nconv .le. 0) then
361
else if (n .le. 0) then
363
else if (nev .le. 0) then
365
else if (ncv .le. nev+1 .or. ncv .gt. n) then
367
else if (which .ne. 'LM' .and.
368
& which .ne. 'SM' .and.
369
& which .ne. 'LR' .and.
370
& which .ne. 'SR' .and.
371
& which .ne. 'LI' .and.
372
& which .ne. 'SI') then
374
else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
376
else if (lworkl .lt. 3*ncv**2 + 4*ncv) then
378
else if ( (howmny .ne. 'A' .and.
379
& howmny .ne. 'P' .and.
380
& howmny .ne. 'S') .and. rvec ) then
382
else if (howmny .eq. 'S' ) then
386
if (mode .eq. 1 .or. mode .eq. 2) then
388
else if (mode .eq. 3 ) then
393
if (mode .eq. 1 .and. bmat .eq. 'G') ierr = -11
399
if (ierr .ne. 0) then
404
c %--------------------------------------------------------%
405
c | Pointer into WORKL for address of H, RITZ, WORKEV, Q |
406
c | etc... and the remaining workspace. |
407
c | Also update pointer to be used on output. |
408
c | Memory is laid out as follows: |
409
c | workl(1:ncv*ncv) := generated Hessenberg matrix |
410
c | workl(ncv*ncv+1:ncv*ncv+ncv) := ritz values |
411
c | workl(ncv*ncv+ncv+1:ncv*ncv+2*ncv) := error bounds |
412
c %--------------------------------------------------------%
414
c %-----------------------------------------------------------%
415
c | The following is used and set by ZNEUPD . |
416
c | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed |
418
c | workl(ncv*ncv+3*ncv+1:ncv*ncv+4*ncv) := The untransformed |
419
c | error bounds of |
420
c | the Ritz values |
421
c | workl(ncv*ncv+4*ncv+1:2*ncv*ncv+4*ncv) := Holds the upper |
422
c | triangular matrix |
424
c | workl(2*ncv*ncv+4*ncv+1: 3*ncv*ncv+4*ncv) := Holds the |
425
c | associated matrix |
426
c | representation of |
428
c | subspace for H. |
429
c | GRAND total of NCV * ( 3 * NCV + 4 ) locations. |
430
c %-----------------------------------------------------------%
441
invsub = iuptri + ldh*ncv
449
c %-----------------------------------------%
450
c | irz points to the Ritz values computed |
451
c | by _neigh before exiting _naup2. |
452
c | ibd points to the Ritz estimates |
453
c | computed by _neigh before exiting |
455
c %-----------------------------------------%
457
irz = ipntr(14) + ncv*ncv
460
c %------------------------------------%
461
c | RNORM is B-norm of the RESID(1:N). |
462
c %------------------------------------%
467
if (msglvl .gt. 2) then
468
call zvout (logfil, ncv, workl(irz), ndigit,
469
& '_neupd: Ritz values passed in from _NAUPD.')
470
call zvout (logfil, ncv, workl(ibd), ndigit,
471
& '_neupd: Ritz estimates passed in from _NAUPD.')
478
c %---------------------------------------------------%
479
c | Use the temporary bounds array to store indices |
480
c | These will be used to mark the select array later |
481
c %---------------------------------------------------%
484
workl(bounds+j-1) = j
488
c %-------------------------------------%
489
c | Select the wanted Ritz values. |
490
c | Sort the Ritz values so that the |
491
c | wanted ones appear at the tailing |
492
c | NEV positions of workl(irr) and |
493
c | workl(iri). Move the corresponding |
494
c | error estimates in workl(ibd) |
496
c %-------------------------------------%
500
call zngets (ishift, which , nev ,
501
& np , workl(irz), workl(bounds))
503
if (msglvl .gt. 2) then
504
call zvout (logfil, ncv, workl(irz), ndigit,
505
& '_neupd: Ritz values after calling _NGETS.')
506
call zvout (logfil, ncv, workl(bounds), ndigit,
507
& '_neupd: Ritz value indices after calling _NGETS.')
510
c %-----------------------------------------------------%
511
c | Record indices of the converged wanted Ritz values |
512
c | Mark the select array for possible reordering |
513
c %-----------------------------------------------------%
518
& dlapy2 ( dble (workl(irz+ncv-j)),
519
& dimag (workl(irz+ncv-j)) ))
520
jj = workl(bounds + ncv - j)
521
if (numcnv .lt. nconv .and.
522
& dlapy2 ( dble (workl(ibd+jj-1)),
523
& dimag (workl(ibd+jj-1)) )
524
& .le. tol*rtemp) then
527
if (jj .gt. nev) reord = .true.
531
c %-----------------------------------------------------------%
532
c | Check the count (numcnv) of converged Ritz values with |
533
c | the number (nconv) reported by dnaupd. If these two |
534
c | are different then there has probably been an error |
535
c | caused by incorrect passing of the dnaupd data. |
536
c %-----------------------------------------------------------%
538
if (msglvl .gt. 2) then
539
call ivout(logfil, 1, numcnv, ndigit,
540
& '_neupd: Number of specified eigenvalues')
541
call ivout(logfil, 1, nconv, ndigit,
542
& '_neupd: Number of "converged" eigenvalues')
545
if (numcnv .ne. nconv) then
550
c %-------------------------------------------------------%
551
c | Call LAPACK routine zlahqr to compute the Schur form |
552
c | of the upper Hessenberg matrix returned by ZNAUPD . |
553
c | Make a copy of the upper Hessenberg matrix. |
554
c | Initialize the Schur vector matrix Q to the identity. |
555
c %-------------------------------------------------------%
557
call zcopy (ldh*ncv, workl(ih), 1, workl(iuptri), 1)
558
call zlaset ('All', ncv, ncv ,
559
& zero , one, workl(invsub),
561
call zlahqr (.true., .true. , ncv ,
562
& 1 , ncv , workl(iuptri),
563
& ldh , workl(iheig) , 1 ,
564
& ncv , workl(invsub), ldq ,
566
call zcopy (ncv , workl(invsub+ncv-1), ldq,
569
if (ierr .ne. 0) then
574
if (msglvl .gt. 1) then
575
call zvout (logfil, ncv, workl(iheig), ndigit,
576
& '_neupd: Eigenvalues of H')
577
call zvout (logfil, ncv, workl(ihbds), ndigit,
578
& '_neupd: Last row of the Schur vector matrix')
579
if (msglvl .gt. 3) then
580
call zmout (logfil , ncv, ncv ,
581
& workl(iuptri), ldh, ndigit,
582
& '_neupd: The upper triangular matrix ')
588
c %-----------------------------------------------%
589
c | Reorder the computed upper triangular matrix. |
590
c %-----------------------------------------------%
592
call ztrsen ('None' , 'V' , select ,
593
& ncv , workl(iuptri), ldh ,
594
& workl(invsub), ldq , workl(iheig),
595
& nconv , conds , sep ,
596
& workev , ncv , ierr)
598
if (ierr .eq. 1) then
603
if (msglvl .gt. 2) then
604
call zvout (logfil, ncv, workl(iheig), ndigit,
605
& '_neupd: Eigenvalues of H--reordered')
606
if (msglvl .gt. 3) then
607
call zmout (logfil , ncv, ncv ,
608
& workl(iuptri), ldq, ndigit,
609
& '_neupd: Triangular matrix after re-ordering')
615
c %---------------------------------------------%
616
c | Copy the last row of the Schur basis matrix |
617
c | to workl(ihbds). This vector will be used |
618
c | to compute the Ritz estimates of converged |
620
c %---------------------------------------------%
622
call zcopy (ncv , workl(invsub+ncv-1), ldq,
625
c %--------------------------------------------%
626
c | Place the computed eigenvalues of H into D |
627
c | if a spectral transformation was not used. |
628
c %--------------------------------------------%
630
if (type .eq. 'REGULR') then
631
call zcopy (nconv, workl(iheig), 1, d, 1)
634
c %----------------------------------------------------------%
635
c | Compute the QR factorization of the matrix representing |
636
c | the wanted invariant subspace located in the first NCONV |
637
c | columns of workl(invsub,ldq). |
638
c %----------------------------------------------------------%
640
call zgeqr2 (ncv , nconv , workl(invsub),
641
& ldq , workev, workev(ncv+1),
644
c %--------------------------------------------------------%
645
c | * Postmultiply V by Q using zunm2r . |
646
c | * Copy the first NCONV columns of VQ into Z. |
647
c | * Postmultiply Z by R. |
648
c | The N by NCONV matrix Z is now a matrix representation |
649
c | of the approximate invariant subspace associated with |
650
c | the Ritz values in workl(iheig). The first NCONV |
651
c | columns of V are now approximate Schur vectors |
652
c | associated with the upper triangular matrix of order |
653
c | NCONV in workl(iuptri). |
654
c %--------------------------------------------------------%
656
call zunm2r ('Right', 'Notranspose', n ,
657
& ncv , nconv , workl(invsub),
659
& ldv , workd(n+1) , ierr)
660
call zlacpy ('All', n, nconv, v, ldv, z, ldz)
664
c %---------------------------------------------------%
665
c | Perform both a column and row scaling if the |
666
c | diagonal element of workl(invsub,ldq) is negative |
667
c | I'm lazy and don't take advantage of the upper |
668
c | triangular form of workl(iuptri,ldq). |
669
c | Note that since Q is orthogonal, R is a diagonal |
670
c | matrix consisting of plus or minus ones. |
671
c %---------------------------------------------------%
673
if ( dble ( workl(invsub+(j-1)*ldq+j-1) ) .lt.
675
call zscal (nconv, -one, workl(iuptri+j-1), ldq)
676
call zscal (nconv, -one, workl(iuptri+(j-1)*ldq), 1)
681
if (howmny .eq. 'A') then
683
c %--------------------------------------------%
684
c | Compute the NCONV wanted eigenvectors of T |
685
c | located in workl(iuptri,ldq). |
686
c %--------------------------------------------%
689
if (j .le. nconv) then
696
call ztrevc ('Right', 'Select' , select ,
697
& ncv , workl(iuptri), ldq ,
698
& vl , 1 , workl(invsub),
699
& ldq , ncv , outncv ,
700
& workev , rwork , ierr)
702
if (ierr .ne. 0) then
707
c %------------------------------------------------%
708
c | Scale the returning eigenvectors so that their |
709
c | Euclidean norms are all one. LAPACK subroutine |
710
c | ztrevc returns each eigenvector normalized so |
711
c | that the element of largest magnitude has |
713
c %------------------------------------------------%
716
rtemp = dznrm2 (ncv, workl(invsub+(j-1)*ldq), 1)
717
rtemp = dble (one) / rtemp
718
call zdscal ( ncv, rtemp,
719
& workl(invsub+(j-1)*ldq), 1 )
721
c %------------------------------------------%
722
c | Ritz estimates can be obtained by taking |
723
c | the inner product of the last row of the |
724
c | Schur basis of H with eigenvectors of T. |
725
c | Note that the eigenvector matrix of T is |
726
c | upper triangular, thus the length of the |
727
c | inner product can be set to j. |
728
c %------------------------------------------%
730
workev(j) = zdotc (j, workl(ihbds), 1,
731
& workl(invsub+(j-1)*ldq), 1)
734
if (msglvl .gt. 2) then
735
call zcopy (nconv, workl(invsub+ncv-1), ldq,
737
call zvout (logfil, nconv, workl(ihbds), ndigit,
738
& '_neupd: Last row of the eigenvector matrix for T')
739
if (msglvl .gt. 3) then
740
call zmout (logfil , ncv, ncv ,
741
& workl(invsub), ldq, ndigit,
742
& '_neupd: The eigenvector matrix for T')
746
c %---------------------------------------%
747
c | Copy Ritz estimates into workl(ihbds) |
748
c %---------------------------------------%
750
call zcopy (nconv, workev, 1, workl(ihbds), 1)
752
c %----------------------------------------------%
753
c | The eigenvector matrix Q of T is triangular. |
755
c %----------------------------------------------%
757
call ztrmm ('Right' , 'Upper' , 'No transpose',
758
& 'Non-unit', n , nconv ,
759
& one , workl(invsub), ldq ,
765
c %--------------------------------------------------%
766
c | An approximate invariant subspace is not needed. |
767
c | Place the Ritz values computed ZNAUPD into D. |
768
c %--------------------------------------------------%
770
call zcopy (nconv, workl(ritz), 1, d, 1)
771
call zcopy (nconv, workl(ritz), 1, workl(iheig), 1)
772
call zcopy (nconv, workl(bounds), 1, workl(ihbds), 1)
776
c %------------------------------------------------%
777
c | Transform the Ritz values and possibly vectors |
778
c | and corresponding error bounds of OP to those |
779
c | of A*x = lambda*B*x. |
780
c %------------------------------------------------%
782
if (type .eq. 'REGULR') then
785
& call zscal (ncv, rnorm, workl(ihbds), 1)
789
c %---------------------------------------%
790
c | A spectral transformation was used. |
791
c | * Determine the Ritz estimates of the |
792
c | Ritz values in the original system. |
793
c %---------------------------------------%
796
& call zscal (ncv, rnorm, workl(ihbds), 1)
799
temp = workl(iheig+k-1)
800
workl(ihbds+k-1) = workl(ihbds+k-1) / temp / temp
805
c %-----------------------------------------------------------%
806
c | * Transform the Ritz values back to the original system. |
807
c | For TYPE = 'SHIFTI' the transformation is |
808
c | lambda = 1/theta + sigma |
810
c | *The Ritz vectors are not affected by the transformation. |
811
c %-----------------------------------------------------------%
813
if (type .eq. 'SHIFTI') then
815
d(k) = one / workl(iheig+k-1) + sigma
819
if (type .ne. 'REGULR' .and. msglvl .gt. 1) then
820
call zvout (logfil, nconv, d, ndigit,
821
& '_neupd: Untransformed Ritz values.')
822
call zvout (logfil, nconv, workl(ihbds), ndigit,
823
& '_neupd: Ritz estimates of the untransformed Ritz values.')
824
else if ( msglvl .gt. 1) then
825
call zvout (logfil, nconv, d, ndigit,
826
& '_neupd: Converged Ritz values.')
827
call zvout (logfil, nconv, workl(ihbds), ndigit,
828
& '_neupd: Associated Ritz estimates.')
831
c %-------------------------------------------------%
832
c | Eigenvector Purification step. Formally perform |
833
c | one of inverse subspace iteration. Only used |
834
c | for MODE = 3. See reference 3. |
835
c %-------------------------------------------------%
837
if (rvec .and. howmny .eq. 'A' .and. type .eq. 'SHIFTI') then
839
c %------------------------------------------------%
840
c | Purify the computed Ritz vectors by adding a |
841
c | little bit of the residual vector: |
843
c | resid(:)*( e s ) / theta |
845
c | where H s = s theta. |
846
c %------------------------------------------------%
849
if (workl(iheig+j-1) .ne. zero) then
850
workev(j) = workl(invsub+(j-1)*ldq+ncv-1) /
855
c %---------------------------------------%
856
c | Perform a rank one update to Z and |
857
c | purify all the Ritz vectors together. |
858
c %---------------------------------------%
860
call zgeru (n, nconv, one, resid, 1, workev, 1, z, ldz)
added
removed
.pc/fix-segfault.patch/ARPACK/SRC/zneupd.f
