← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/python-scipy/oneiric-proposed

« back to all changes in this revision

Viewing changes to scipy/sandbox/arpack/ARPACK/SRC/zneupd.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sandbox/arpack/ARPACK/SRC/zneupd.f
 
1
c\BeginDoc
 
2
c 
 
3
c\Name: zneupd 
 
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 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 
 
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 ZNAUPD.  For a brief 
 
33
c  description, see definitions of IPARAM(7), MODE and WHICH in the
 
34
c  documentation of ZNAUPD.
 
35
c
 
36
c\Usage:
 
37
c  call zneupd 
 
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*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.
 
72
c
 
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
 
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 ZNAUPD.  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*16  (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*16 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 ZNAUPD 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 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.
 
108
c
 
109
c  Three of these parameters (V, WORKL and INFO) are also output parameters:
 
110
c
 
111
c  V       Complex*16 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 ZNAUPD .
 
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   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.
 
134
c
 
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
 
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                     zneupd 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 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
 
160
c                occurs.
 
161
c
 
162
c          = -1: N must be positive.
 
163
c          = -2: NEV must be positive.
 
164
c          = -3: NCV-NEV >= 1 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
 
177
c                 accuracy.
 
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
 
182
c                 ZNEUPD
 
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     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 
 
202
c             a matrix.
 
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 
 
211
c             factored form.
 
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.
 
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.8   DATE OF SID: 07/21/02   RELEASE: 2
 
244
c
 
245
c\EndLib
 
246
c
 
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  )
 
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*16     
 
270
     &           sigma
 
271
      Double precision 
 
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
      Double precision
 
281
     &           rwork(ncv)
 
282
      Complex*16
 
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*16
 
292
     &           one, zero
 
293
      parameter  (one = (1.0D+0, 0.0D+0), zero = (0.0D+0, 0.0D+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*16
 
306
     &           rnorm, temp, vl(1)
 
307
      Double precision
 
308
     &           conds, sep, rtemp, eps23
 
309
      logical    reord
 
310
c
 
311
c     %----------------------%
 
312
c     | External Subroutines |
 
313
c     %----------------------%
 
314
c
 
315
      external   zcopy , zgeru, zgeqr2, zlacpy, zmout,
 
316
     &           zunm2r, ztrmm, zvout, ivout,
 
317
     &           zlahqr
 
318
c  
 
319
c     %--------------------%
 
320
c     | External Functions |
 
321
c     %--------------------%
 
322
c
 
323
      Double precision
 
324
     &           dznrm2, dlamch, dlapy2
 
325
      external   dznrm2, dlamch, dlapy2
 
326
c
 
327
      Complex*16
 
328
     &           zdotc
 
329
      external   zdotc
 
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 = dlamch('Epsilon-Machine')
 
350
      eps23 = eps23**(2.0D+0 / 3.0D+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 .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 ZNEUPD.                 |
 
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 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.')
 
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 zngets(ishift, which     , nev          ,
 
501
     &                np    , workl(irz), workl(bounds))
 
502
c
 
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.')
 
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
     &                 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
 
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 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        %-------------------------------------------------------%
 
556
c
 
557
         call zcopy(ldh*ncv, workl(ih), 1, workl(iuptri), 1)
 
558
         call zlaset('All', ncv, ncv          , 
 
559
     &                zero , one, workl(invsub),
 
560
     &                ldq)
 
561
         call zlahqr(.true., .true.       , ncv          , 
 
562
     &                1     , ncv          , workl(iuptri),
 
563
     &                ldh   , workl(iheig) , 1            ,
 
564
     &                ncv   , workl(invsub), ldq          ,
 
565
     &                ierr)
 
566
         call zcopy(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 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 ')
 
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 ztrsen('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 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')
 
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 zcopy(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 zcopy(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 zgeqr2(ncv , nconv , workl(invsub),
 
641
     &                ldq , workev, workev(ncv+1),
 
642
     &                ierr)
 
643
c
 
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        %--------------------------------------------------------%
 
655
c
 
656
         call zunm2r('Right', 'Notranspose', n            ,
 
657
     &                ncv    , nconv        , workl(invsub),
 
658
     &                ldq    , workev       , v            ,
 
659
     &                ldv    , workd(n+1)   , ierr)
 
660
         call zlacpy('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 ( dble( workl(invsub+(j-1)*ldq+j-1) ) .lt. 
 
674
     &                  dble(zero) ) then
 
675
               call zscal(nconv, -one, workl(iuptri+j-1), ldq)
 
676
               call zscal(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 ztrevc('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           | ztrevc 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 = 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 )
 
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) = zdotc(j, workl(ihbds), 1,
 
731
     &                        workl(invsub+(j-1)*ldq), 1)
 
732
 40         continue
 
733
c
 
734
            if (msglvl .gt. 2) then
 
735
               call zcopy(nconv, workl(invsub+ncv-1), ldq,
 
736
     &                    workl(ihbds), 1)
 
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')
 
743
               end if
 
744
            end if
 
745
c
 
746
c           %---------------------------------------%
 
747
c           | Copy Ritz estimates into workl(ihbds) |
 
748
c           %---------------------------------------%
 
749
c 
 
750
            call zcopy(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 ztrmm('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 ZNAUPD into D.    |
 
768
c        %--------------------------------------------------%
 
769
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)
 
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 zscal(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 zscal(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 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.')
 
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 zgeru (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 zneupd|
 
870
c     %---------------%
 
871
c
 
872
      end