Mercurial > hg > octave-thorsten

comparison libcruft/arpack/src/cneupd.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 9ed4018d538c
comparison
equal deleted inserted replaced
12273:83133b5bf392 12274:9f5d2ef078e8
1 c\BeginDoc
2 c
3 c\Name: cneupd
4 c
5 c\Description:
6 c This subroutine returns the converged approximations to eigenvalues
7 c of A*z = lambda*B*z and (optionally):
8 c
9 c (1) The corresponding approximate eigenvectors;
10 c
11 c (2) An orthonormal basis for the associated approximate
12 c invariant subspace;
13 c
14 c (3) Both.
15 c
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).
19 c
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 CNAUPD. CNAUPD 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
28 c Schur basis.
29 c
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 CNAUPD. For a brief
33 c description, see definitions of IPARAM(7), MODE and WHICH in the
34 c documentation of CNAUPD.
35 c
36 c\Usage:
37 c call cneupd
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 )
41 c
42 c\Arguments:
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.
47 c
48 c RVEC = .FALSE. Compute Ritz values only.
49 c
50 c RVEC = .TRUE. Compute Ritz vectors or Schur vectors.
51 c See Remarks below.
52 c
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.
56 c
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.
61 c
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.
68 c
69 c D Complex 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.
72 c
73 c Z Complex 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
77 c A*z = lambda*B*z.
78 c
79 c If RVEC = .FALSE. or HOWMNY = 'P', then Z is NOT REFERENCED.
80 c
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 CNAUPD. In this case the Arnoldi basis
84 c will be destroyed and overwritten with the eigenvector basis.
85 c
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.
90 c
91 c SIGMA Complex (INPUT)
92 c If IPARAM(7) = 3 then SIGMA represents the shift.
93 c Not referenced if IPARAM(7) = 1 or 2.
94 c
95 c WORKEV Complex work array of dimension 2*NCV. (WORKSPACE)
96 c
97 c **** The remaining arguments MUST be the same as for the ****
98 c **** call to CNAUPD that was just completed. ****
99 c
100 c NOTE: The remaining arguments
101 c
102 c BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, IPNTR,
103 c WORKD, WORKL, LWORKL, RWORK, INFO
104 c
105 c must be passed directly to CNEUPD following the last call
106 c to CNAUPD. These arguments MUST NOT BE MODIFIED between
107 c the the last call to CNAUPD and the call to CNEUPD.
108 c
109 c Three of these parameters (V, WORKL and INFO) are also output parameters:
110 c
111 c V Complex N by NCV array. (INPUT/OUTPUT)
112 c
113 c Upon INPUT: the NCV columns of V contain the Arnoldi basis
114 c vectors for OP as constructed by CNAUPD .
115 c
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.
119 c
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.
126 c
127 c WORKL Real work array of length LWORKL. (OUTPUT/WORKSPACE)
128 c WORKL(1:ncv*ncv+2*ncv) contains information obtained in
129 c cnaupd. They are not changed by cneupd.
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.
134 c
135 c Note: IPNTR(9:13) contains the pointer into WORKL for addresses
136 c of the above information computed by cneupd.
137 c -------------------------------------------------------------
138 c IPNTR(9): pointer to the NCV RITZ values of the
139 c original system.
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 cneupd if RVEC = .TRUE. See Remark 2 below.
147 c -------------------------------------------------------------
148 c
149 c INFO Integer. (OUTPUT)
150 c Error flag on output.
151 c = 0: Normal exit.
152 c
153 c = 1: The Schur form computed by LAPACK routine csheqr
154 c could not be reordered by LAPACK routine ctrsen.
155 c Re-enter subroutine cneupd 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
160 c occurs.
161 c
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 ctrevc.
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: CNAUPD did not find any eigenvalues to sufficient
177 c accuracy.
178 c = -15: CNEUPD got a different count of the number of converged
179 c Ritz values than CNAUPD got. This indicates the user
180 c probably made an error in passing data from CNAUPD to
181 c CNEUPD or that the data was modified before entering
182 c CNEUPD
183 c
184 c\BeginLib
185 c
186 c\References:
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),
189 c pp 357-385.
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.
196 c
197 c\Routines called:
198 c ivout ARPACK utility routine that prints integers.
199 c cmout ARPACK utility routine that prints matrices
200 c cvout ARPACK utility routine that prints vectors.
201 c cgeqr2 LAPACK routine that computes the QR factorization of
202 c a matrix.
203 c clacpy LAPACK matrix copy routine.
204 c clahqr LAPACK routine that computes the Schur form of a
205 c upper Hessenberg matrix.
206 c claset LAPACK matrix initialization routine.
207 c ctrevc LAPACK routine to compute the eigenvectors of a matrix
208 c in upper triangular form.
209 c ctrsen LAPACK routine that re-orders the Schur form.
210 c cunm2r LAPACK routine that applies an orthogonal matrix in
211 c factored form.
212 c slamch LAPACK routine that determines machine constants.
213 c ctrmm Level 3 BLAS matrix times an upper triangular matrix.
214 c cgeru Level 2 BLAS rank one update to a matrix.
215 c ccopy Level 1 BLAS that copies one vector to another .
216 c cscal Level 1 BLAS that scales a vector.
217 c csscal Level 1 BLAS that scales a complex vector by a real number.
218 c scnrm2 Level 1 BLAS that computes the norm of a complex vector.
219 c
220 c\Remarks
221 c
222 c 1. Currently only HOWMNY = 'A' and 'P' are implemented.
223 c
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)).
232 c
233 c\Authors
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
239 c Rice University
240 c Houston, Texas
241 c
242 c\SCCS Information: @(#)
243 c FILE: neupd.F SID: 2.7 DATE OF SID: 09/20/00 RELEASE: 2
244 c
245 c\EndLib
246 c
247 c-----------------------------------------------------------------------
248 subroutine cneupd(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 )
254 c
255 c %----------------------------------------------------%
256 c | Include files for debugging and timing information |
257 c %----------------------------------------------------%
258 c
259 include 'debug.h'
260 include 'stat.h'
261 c
262 c %------------------%
263 c | Scalar Arguments |
264 c %------------------%
265 c
266 character bmat, howmny, which*2
267 logical rvec
268 integer info, ldz, ldv, lworkl, n, ncv, nev
269 Complex
270 & sigma
271 Real
272 & tol
273 c
274 c %-----------------%
275 c | Array Arguments |
276 c %-----------------%
277 c
278 integer iparam(11), ipntr(14)
279 logical select(ncv)
280 Real
281 & rwork(ncv)
282 Complex
283 & d(nev) , resid(n) , v(ldv,ncv),
284 & z(ldz, nev),
285 & workd(3*n) , workl(lworkl), workev(2*ncv)
286 c
287 c %------------%
288 c | Parameters |
289 c %------------%
290 c
291 Complex
292 & one, zero
293 parameter (one = (1.0E+0, 0.0E+0) , zero = (0.0E+0, 0.0E+0) )
294 c
295 c %---------------%
296 c | Local Scalars |
297 c %---------------%
298 c
299 character type*6
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 ,
304 & ishift
305 Complex
306 & rnorm, temp, vl(1)
307 Real
308 & conds, sep, rtemp, eps23
309 logical reord
310 c
311 c %----------------------%
312 c | External Subroutines |
313 c %----------------------%
314 c
315 external ccopy , cgeru, cgeqr2, clacpy, cmout,
316 & cunm2r, ctrmm, cvout, ivout,
317 & clahqr
318 c
319 c %--------------------%
320 c | External Functions |
321 c %--------------------%
322 c
323 Real
324 & scnrm2, slamch, slapy2
325 external scnrm2, slamch, slapy2
326 c
327 Complex
328 & cdotc
329 external cdotc
330 c
331 c %-----------------------%
332 c | Executable Statements |
333 c %-----------------------%
334 c
335 c %------------------------%
336 c | Set default parameters |
337 c %------------------------%
338 c
339 msglvl = mceupd
340 mode = iparam(7)
341 nconv = iparam(5)
342 info = 0
343 c
344 c
345 c %---------------------------------%
346 c | Get machine dependent constant. |
347 c %---------------------------------%
348 c
349 eps23 = slamch('Epsilon-Machine')
350 eps23 = eps23**(2.0E+0 / 3.0E+0 )
351 c
352 c %-------------------------------%
353 c | Quick return |
354 c | Check for incompatible input |
355 c %-------------------------------%
356 c
357 ierr = 0
358 c
359 if (nconv .le. 0) then
360 ierr = -14
361 else if (n .le. 0) then
362 ierr = -1
363 else if (nev .le. 0) then
364 ierr = -2
365 else if (ncv .le. nev+1 .or. ncv .gt. n) then
366 ierr = -3
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
373 ierr = -5
374 else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
375 ierr = -6
376 else if (lworkl .lt. 3*ncv**2 + 4*ncv) then
377 ierr = -7
378 else if ( (howmny .ne. 'A' .and.
379 & howmny .ne. 'P' .and.
380 & howmny .ne. 'S') .and. rvec ) then
381 ierr = -13
382 else if (howmny .eq. 'S' ) then
383 ierr = -12
384 end if
385 c
386 if (mode .eq. 1 .or. mode .eq. 2) then
387 type = 'REGULR'
388 else if (mode .eq. 3 ) then
389 type = 'SHIFTI'
390 else
391 ierr = -10
392 end if
393 if (mode .eq. 1 .and. bmat .eq. 'G') ierr = -11
394 c
395 c %------------%
396 c | Error Exit |
397 c %------------%
398 c
399 if (ierr .ne. 0) then
400 info = ierr
401 go to 9000
402 end if
403 c
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 %--------------------------------------------------------%
413 c
414 c %-----------------------------------------------------------%
415 c | The following is used and set by CNEUPD. |
416 c | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := The untransformed |
417 c | Ritz values. |
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 |
423 c | for H. |
424 c | workl(2*ncv*ncv+4*ncv+1: 3*ncv*ncv+4*ncv) := Holds the |
425 c | associated matrix |
426 c | representation of |
427 c | the invariant |
428 c | subspace for H. |
429 c | GRAND total of NCV * ( 3 * NCV + 4 ) locations. |
430 c %-----------------------------------------------------------%
431 c
432 ih = ipntr(5)
433 ritz = ipntr(6)
434 iq = ipntr(7)
435 bounds = ipntr(8)
436 ldh = ncv
437 ldq = ncv
438 iheig = bounds + ldh
439 ihbds = iheig + ldh
440 iuptri = ihbds + ldh
441 invsub = iuptri + ldh*ncv
442 ipntr(9) = iheig
443 ipntr(11) = ihbds
444 ipntr(12) = iuptri
445 ipntr(13) = invsub
446 wr = 1
447 iwev = wr + ncv
448 c
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 |
454 c | _naup2. |
455 c %-----------------------------------------%
456 c
457 irz = ipntr(14) + ncv*ncv
458 ibd = irz + ncv
459 c
460 c %------------------------------------%
461 c | RNORM is B-norm of the RESID(1:N). |
462 c %------------------------------------%
463 c
464 rnorm = workl(ih+2)
465 workl(ih+2) = zero
466 c
467 if (msglvl .gt. 2) then
468 call cvout(logfil, ncv, workl(irz), ndigit,
469 & '_neupd: Ritz values passed in from _NAUPD.')
470 call cvout(logfil, ncv, workl(ibd), ndigit,
471 & '_neupd: Ritz estimates passed in from _NAUPD.')
472 end if
473 c
474 if (rvec) then
475 c
476 reord = .false.
477 c
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 %---------------------------------------------------%
482 c
483 do 10 j = 1,ncv
484 workl(bounds+j-1) = j
485 select(j) = .false.
486 10 continue
487 c
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) |
495 c | accordingly. |
496 c %-------------------------------------%
497 c
498 np = ncv - nev
499 ishift = 0
500 call cngets(ishift, which , nev ,
501 & np , workl(irz), workl(bounds))
502 c
503 if (msglvl .gt. 2) then
504 call cvout (logfil, ncv, workl(irz), ndigit,
505 & '_neupd: Ritz values after calling _NGETS.')
506 call cvout (logfil, ncv, workl(bounds), ndigit,
507 & '_neupd: Ritz value indices after calling _NGETS.')
508 end if
509 c
510 c %-----------------------------------------------------%
511 c | Record indices of the converged wanted Ritz values |
512 c | Mark the select array for possible reordering |
513 c %-----------------------------------------------------%
514 c
515 numcnv = 0
516 do 11 j = 1,ncv
517 rtemp = max(eps23,
518 & slapy2 ( real (workl(irz+ncv-j)),
519 & aimag(workl(irz+ncv-j)) ))
520 jj = workl(bounds + ncv - j)
521 if (numcnv .lt. nconv .and.
522 & slapy2( real (workl(ibd+jj-1)),
523 & aimag(workl(ibd+jj-1)) )
524 & .le. tol*rtemp) then
525 select(jj) = .true.
526 numcnv = numcnv + 1
527 if (jj .gt. nev) reord = .true.
528 endif
529 11 continue
530 c
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 %-----------------------------------------------------------%
537 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')
543 end if
544 c
545 if (numcnv .ne. nconv) then
546 info = -15
547 go to 9000
548 end if
549 c
550 c %-------------------------------------------------------%
551 c | Call LAPACK routine clahqr to compute the Schur form |
552 c | of the upper Hessenberg matrix returned by CNAUPD. |
553 c | Make a copy of the upper Hessenberg matrix. |
554 c | Initialize the Schur vector matrix Q to the identity. |
555 c %-------------------------------------------------------%
556 c
557 call ccopy(ldh*ncv, workl(ih), 1, workl(iuptri), 1)
558 call claset('All', ncv, ncv ,
559 & zero , one, workl(invsub),
560 & ldq)
561 call clahqr(.true., .true. , ncv ,
562 & 1 , ncv , workl(iuptri),
563 & ldh , workl(iheig) , 1 ,
564 & ncv , workl(invsub), ldq ,
565 & ierr)
566 call ccopy(ncv , workl(invsub+ncv-1), ldq,
567 & workl(ihbds), 1)
568 c
569 if (ierr .ne. 0) then
570 info = -8
571 go to 9000
572 end if
573 c
574 if (msglvl .gt. 1) then
575 call cvout (logfil, ncv, workl(iheig), ndigit,
576 & '_neupd: Eigenvalues of H')
577 call cvout (logfil, ncv, workl(ihbds), ndigit,
578 & '_neupd: Last row of the Schur vector matrix')
579 if (msglvl .gt. 3) then
580 call cmout (logfil , ncv, ncv ,
581 & workl(iuptri), ldh, ndigit,
582 & '_neupd: The upper triangular matrix ')
583 end if
584 end if
585 c
586 if (reord) then
587 c
588 c %-----------------------------------------------%
589 c | Reorder the computed upper triangular matrix. |
590 c %-----------------------------------------------%
591 c
592 call ctrsen('None' , 'V' , select ,
593 & ncv , workl(iuptri), ldh ,
594 & workl(invsub), ldq , workl(iheig),
595 & nconv , conds , sep ,
596 & workev , ncv , ierr)
597 c
598 if (ierr .eq. 1) then
599 info = 1
600 go to 9000
601 end if
602 c
603 if (msglvl .gt. 2) then
604 call cvout (logfil, ncv, workl(iheig), ndigit,
605 & '_neupd: Eigenvalues of H--reordered')
606 if (msglvl .gt. 3) then
607 call cmout(logfil , ncv, ncv ,
608 & workl(iuptri), ldq, ndigit,
609 & '_neupd: Triangular matrix after re-ordering')
610 end if
611 end if
612 c
613 end if
614 c
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 |
619 c | Ritz values. |
620 c %---------------------------------------------%
621 c
622 call ccopy(ncv , workl(invsub+ncv-1), ldq,
623 & workl(ihbds), 1)
624 c
625 c %--------------------------------------------%
626 c | Place the computed eigenvalues of H into D |
627 c | if a spectral transformation was not used. |
628 c %--------------------------------------------%
629 c
630 if (type .eq. 'REGULR') then
631 call ccopy(nconv, workl(iheig), 1, d, 1)
632 end if
633 c
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 %----------------------------------------------------------%
639 c
640 call cgeqr2(ncv , nconv , workl(invsub),
641 & ldq , workev, workev(ncv+1),
642 & ierr)
643 c
644 c %--------------------------------------------------------%
645 c | * Postmultiply V by Q using cunm2r. |
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 %--------------------------------------------------------%
655 c
656 call cunm2r('Right', 'Notranspose', n ,
657 & ncv , nconv , workl(invsub),
658 & ldq , workev , v ,
659 & ldv , workd(n+1) , ierr)
660 call clacpy('All', n, nconv, v, ldv, z, ldz)
661 c
662 do 20 j=1, nconv
663 c
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 %---------------------------------------------------%
672 c
673 if ( real ( workl(invsub+(j-1)*ldq+j-1) ) .lt.
674 & real (zero) ) then
675 call cscal(nconv, -one, workl(iuptri+j-1), ldq)
676 call cscal(nconv, -one, workl(iuptri+(j-1)*ldq), 1)
677 end if
678 c
679 20 continue
680 c
681 if (howmny .eq. 'A') then
682 c
683 c %--------------------------------------------%
684 c | Compute the NCONV wanted eigenvectors of T |
685 c | located in workl(iuptri,ldq). |
686 c %--------------------------------------------%
687 c
688 do 30 j=1, ncv
689 if (j .le. nconv) then
690 select(j) = .true.
691 else
692 select(j) = .false.
693 end if
694 30 continue
695 c
696 call ctrevc('Right', 'Select' , select ,
697 & ncv , workl(iuptri), ldq ,
698 & vl , 1 , workl(invsub),
699 & ldq , ncv , outncv ,
700 & workev , rwork , ierr)
701 c
702 if (ierr .ne. 0) then
703 info = -9
704 go to 9000
705 end if
706 c
707 c %------------------------------------------------%
708 c | Scale the returning eigenvectors so that their |
709 c | Euclidean norms are all one. LAPACK subroutine |
710 c | ctrevc returns each eigenvector normalized so |
711 c | that the element of largest magnitude has |
712 c | magnitude 1. |
713 c %------------------------------------------------%
714 c
715 do 40 j=1, nconv
716 rtemp = scnrm2(ncv, workl(invsub+(j-1)*ldq), 1)
717 rtemp = real (one) / rtemp
718 call csscal ( ncv, rtemp,
719 & workl(invsub+(j-1)*ldq), 1 )
720 c
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 %------------------------------------------%
729 c
730 workev(j) = cdotc(j, workl(ihbds), 1,
731 & workl(invsub+(j-1)*ldq), 1)
732 40 continue
733 c
734 if (msglvl .gt. 2) then
735 call ccopy(nconv, workl(invsub+ncv-1), ldq,
736 & workl(ihbds), 1)
737 call cvout (logfil, nconv, workl(ihbds), ndigit,
738 & '_neupd: Last row of the eigenvector matrix for T')
739 if (msglvl .gt. 3) then
740 call cmout(logfil , ncv, ncv ,
741 & workl(invsub), ldq, ndigit,
742 & '_neupd: The eigenvector matrix for T')
743 end if
744 end if
745 c
746 c %---------------------------------------%
747 c | Copy Ritz estimates into workl(ihbds) |
748 c %---------------------------------------%
749 c
750 call ccopy(nconv, workev, 1, workl(ihbds), 1)
751 c
752 c %----------------------------------------------%
753 c | The eigenvector matrix Q of T is triangular. |
754 c | Form Z*Q. |
755 c %----------------------------------------------%
756 c
757 call ctrmm('Right' , 'Upper' , 'No transpose',
758 & 'Non-unit', n , nconv ,
759 & one , workl(invsub), ldq ,
760 & z , ldz)
761 end if
762 c
763 else
764 c
765 c %--------------------------------------------------%
766 c | An approximate invariant subspace is not needed. |
767 c | Place the Ritz values computed CNAUPD into D. |
768 c %--------------------------------------------------%
769 c
770 call ccopy(nconv, workl(ritz), 1, d, 1)
771 call ccopy(nconv, workl(ritz), 1, workl(iheig), 1)
772 call ccopy(nconv, workl(bounds), 1, workl(ihbds), 1)
773 c
774 end if
775 c
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 %------------------------------------------------%
781 c
782 if (type .eq. 'REGULR') then
783 c
784 if (rvec)
785 & call cscal(ncv, rnorm, workl(ihbds), 1)
786 c
787 else
788 c
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 %---------------------------------------%
794 c
795 if (rvec)
796 & call cscal(ncv, rnorm, workl(ihbds), 1)
797 c
798 do 50 k=1, ncv
799 temp = workl(iheig+k-1)
800 workl(ihbds+k-1) = workl(ihbds+k-1) / temp / temp
801 50 continue
802 c
803 end if
804 c
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 |
809 c | NOTES: |
810 c | *The Ritz vectors are not affected by the transformation. |
811 c %-----------------------------------------------------------%
812 c
813 if (type .eq. 'SHIFTI') then
814 do 60 k=1, nconv
815 d(k) = one / workl(iheig+k-1) + sigma
816 60 continue
817 end if
818 c
819 if (type .ne. 'REGULR' .and. msglvl .gt. 1) then
820 call cvout (logfil, nconv, d, ndigit,
821 & '_neupd: Untransformed Ritz values.')
822 call cvout (logfil, nconv, workl(ihbds), ndigit,
823 & '_neupd: Ritz estimates of the untransformed Ritz values.')
824 else if ( msglvl .gt. 1) then
825 call cvout (logfil, nconv, d, ndigit,
826 & '_neupd: Converged Ritz values.')
827 call cvout (logfil, nconv, workl(ihbds), ndigit,
828 & '_neupd: Associated Ritz estimates.')
829 end if
830 c
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 %-------------------------------------------------%
836 c
837 if (rvec .and. howmny .eq. 'A' .and. type .eq. 'SHIFTI') then
838 c
839 c %------------------------------------------------%
840 c | Purify the computed Ritz vectors by adding a |
841 c | little bit of the residual vector: |
842 c | T |
843 c | resid(:)*( e s ) / theta |
844 c | NCV |
845 c | where H s = s theta. |
846 c %------------------------------------------------%
847 c
848 do 100 j=1, nconv
849 if (workl(iheig+j-1) .ne. zero) then
850 workev(j) = workl(invsub+(j-1)*ldq+ncv-1) /
851 & workl(iheig+j-1)
852 endif
853 100 continue
854
855 c %---------------------------------------%
856 c | Perform a rank one update to Z and |
857 c | purify all the Ritz vectors together. |
858 c %---------------------------------------%
859 c
860 call cgeru (n, nconv, one, resid, 1, workev, 1, z, ldz)
861 c
862 end if
863 c
864 9000 continue
865 c
866 return
867 c
868 c %---------------%
869 c | End of cneupd|
870 c %---------------%
871 c
872 end