~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/sandbox/arpack/ARPACK/SRC/snaupd.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
* 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:
1
 
c\BeginDoc
2
 
c
3
 
c\Name: snaupd
4
 
c
5
 
c\Description: 
6
 
c  Reverse communication interface for the Implicitly Restarted Arnoldi
7
 
c  iteration. This subroutine computes approximations to a few eigenpairs 
8
 
c  of a linear operator "OP" with respect to a semi-inner product defined by 
9
 
c  a symmetric positive semi-definite real matrix B. B may be the identity 
10
 
c  matrix. NOTE: If the linear operator "OP" is real and symmetric 
11
 
c  with respect to the real positive semi-definite symmetric matrix B, 
12
 
c  i.e. B*OP = (OP`)*B, then subroutine ssaupd should be used instead.
13
 
c
14
 
c  The computed approximate eigenvalues are called Ritz values and
15
 
c  the corresponding approximate eigenvectors are called Ritz vectors.
16
 
c
17
 
c  snaupd is usually called iteratively to solve one of the 
18
 
c  following problems:
19
 
c
20
 
c  Mode 1:  A*x = lambda*x.
21
 
c           ===> OP = A  and  B = I.
22
 
c
23
 
c  Mode 2:  A*x = lambda*M*x, M symmetric positive definite
24
 
c           ===> OP = inv[M]*A  and  B = M.
25
 
c           ===> (If M can be factored see remark 3 below)
26
 
c
27
 
c  Mode 3:  A*x = lambda*M*x, M symmetric semi-definite
28
 
c           ===> OP = Real_Part{ inv[A - sigma*M]*M }  and  B = M. 
29
 
c           ===> shift-and-invert mode (in real arithmetic)
30
 
c           If OP*x = amu*x, then 
31
 
c           amu = 1/2 * [ 1/(lambda-sigma) + 1/(lambda-conjg(sigma)) ].
32
 
c           Note: If sigma is real, i.e. imaginary part of sigma is zero;
33
 
c                 Real_Part{ inv[A - sigma*M]*M } == inv[A - sigma*M]*M 
34
 
c                 amu == 1/(lambda-sigma). 
35
 
c  
36
 
c  Mode 4:  A*x = lambda*M*x, M symmetric semi-definite
37
 
c           ===> OP = Imaginary_Part{ inv[A - sigma*M]*M }  and  B = M. 
38
 
c           ===> shift-and-invert mode (in real arithmetic)
39
 
c           If OP*x = amu*x, then 
40
 
c           amu = 1/2i * [ 1/(lambda-sigma) - 1/(lambda-conjg(sigma)) ].
41
 
c
42
 
c  Both mode 3 and 4 give the same enhancement to eigenvalues close to
43
 
c  the (complex) shift sigma.  However, as lambda goes to infinity,
44
 
c  the operator OP in mode 4 dampens the eigenvalues more strongly than
45
 
c  does OP defined in mode 3.
46
 
c
47
 
c  NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v
48
 
c        should be accomplished either by a direct method
49
 
c        using a sparse matrix factorization and solving
50
 
c
51
 
c           [A - sigma*M]*w = v  or M*w = v,
52
 
c
53
 
c        or through an iterative method for solving these
54
 
c        systems.  If an iterative method is used, the
55
 
c        convergence test must be more stringent than
56
 
c        the accuracy requirements for the eigenvalue
57
 
c        approximations.
58
 
c
59
 
c\Usage:
60
 
c  call snaupd
61
 
c     ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM,
62
 
c       IPNTR, WORKD, WORKL, LWORKL, INFO )
63
 
c
64
 
c\Arguments
65
 
c  IDO     Integer.  (INPUT/OUTPUT)
66
 
c          Reverse communication flag.  IDO must be zero on the first 
67
 
c          call to snaupd.  IDO will be set internally to
68
 
c          indicate the type of operation to be performed.  Control is
69
 
c          then given back to the calling routine which has the
70
 
c          responsibility to carry out the requested operation and call
71
 
c          snaupd with the result.  The operand is given in
72
 
c          WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)).
73
 
c          -------------------------------------------------------------
74
 
c          IDO =  0: first call to the reverse communication interface
75
 
c          IDO = -1: compute  Y = OP * X  where
76
 
c                    IPNTR(1) is the pointer into WORKD for X,
77
 
c                    IPNTR(2) is the pointer into WORKD for Y.
78
 
c                    This is for the initialization phase to force the
79
 
c                    starting vector into the range of OP.
80
 
c          IDO =  1: compute  Y = OP * X  where
81
 
c                    IPNTR(1) is the pointer into WORKD for X,
82
 
c                    IPNTR(2) is the pointer into WORKD for Y.
83
 
c                    In mode 3 and 4, the vector B * X is already
84
 
c                    available in WORKD(ipntr(3)).  It does not
85
 
c                    need to be recomputed in forming OP * X.
86
 
c          IDO =  2: compute  Y = B * X  where
87
 
c                    IPNTR(1) is the pointer into WORKD for X,
88
 
c                    IPNTR(2) is the pointer into WORKD for Y.
89
 
c          IDO =  3: compute the IPARAM(8) real and imaginary parts 
90
 
c                    of the shifts where INPTR(14) is the pointer
91
 
c                    into WORKL for placing the shifts. See Remark
92
 
c                    5 below.
93
 
c          IDO = 99: done
94
 
c          -------------------------------------------------------------
95
 
c             
96
 
c  BMAT    Character*1.  (INPUT)
97
 
c          BMAT specifies the type of the matrix B that defines the
98
 
c          semi-inner product for the operator OP.
99
 
c          BMAT = 'I' -> standard eigenvalue problem A*x = lambda*x
100
 
c          BMAT = 'G' -> generalized eigenvalue problem A*x = lambda*B*x
101
 
c
102
 
c  N       Integer.  (INPUT)
103
 
c          Dimension of the eigenproblem.
104
 
c
105
 
c  WHICH   Character*2.  (INPUT)
106
 
c          'LM' -> want the NEV eigenvalues of largest magnitude.
107
 
c          'SM' -> want the NEV eigenvalues of smallest magnitude.
108
 
c          'LR' -> want the NEV eigenvalues of largest real part.
109
 
c          'SR' -> want the NEV eigenvalues of smallest real part.
110
 
c          'LI' -> want the NEV eigenvalues of largest imaginary part.
111
 
c          'SI' -> want the NEV eigenvalues of smallest imaginary part.
112
 
c
113
 
c  NEV     Integer.  (INPUT/OUTPUT)
114
 
c          Number of eigenvalues of OP to be computed. 0 < NEV < N-1.
115
 
c
116
 
c  TOL     Real scalar.  (INPUT)
117
 
c          Stopping criterion: the relative accuracy of the Ritz value 
118
 
c          is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I))
119
 
c          where ABS(RITZ(I)) is the magnitude when RITZ(I) is complex.
120
 
c          DEFAULT = SLAMCH('EPS')  (machine precision as computed
121
 
c                    by the LAPACK auxiliary subroutine SLAMCH).
122
 
c
123
 
c  RESID   Real array of length N.  (INPUT/OUTPUT)
124
 
c          On INPUT: 
125
 
c          If INFO .EQ. 0, a random initial residual vector is used.
126
 
c          If INFO .NE. 0, RESID contains the initial residual vector,
127
 
c                          possibly from a previous run.
128
 
c          On OUTPUT:
129
 
c          RESID contains the final residual vector.
130
 
c
131
 
c  NCV     Integer.  (INPUT)
132
 
c          Number of columns of the matrix V. NCV must satisfy the two
133
 
c          inequalities 2 <= NCV-NEV and NCV <= N.
134
 
c          This will indicate how many Arnoldi vectors are generated 
135
 
c          at each iteration.  After the startup phase in which NEV 
136
 
c          Arnoldi vectors are generated, the algorithm generates 
137
 
c          approximately NCV-NEV Arnoldi vectors at each subsequent update 
138
 
c          iteration. Most of the cost in generating each Arnoldi vector is 
139
 
c          in the matrix-vector operation OP*x. 
140
 
c          NOTE: 2 <= NCV-NEV in order that complex conjugate pairs of Ritz 
141
 
c          values are kept together. (See remark 4 below)
142
 
c
143
 
c  V       Real array N by NCV.  (OUTPUT)
144
 
c          Contains the final set of Arnoldi basis vectors. 
145
 
c
146
 
c  LDV     Integer.  (INPUT)
147
 
c          Leading dimension of V exactly as declared in the calling program.
148
 
c
149
 
c  IPARAM  Integer array of length 11.  (INPUT/OUTPUT)
150
 
c          IPARAM(1) = ISHIFT: method for selecting the implicit shifts.
151
 
c          The shifts selected at each iteration are used to restart
152
 
c          the Arnoldi iteration in an implicit fashion.
153
 
c          -------------------------------------------------------------
154
 
c          ISHIFT = 0: the shifts are provided by the user via
155
 
c                      reverse communication.  The real and imaginary
156
 
c                      parts of the NCV eigenvalues of the Hessenberg
157
 
c                      matrix H are returned in the part of the WORKL 
158
 
c                      array corresponding to RITZR and RITZI. See remark 
159
 
c                      5 below.
160
 
c          ISHIFT = 1: exact shifts with respect to the current
161
 
c                      Hessenberg matrix H.  This is equivalent to 
162
 
c                      restarting the iteration with a starting vector
163
 
c                      that is a linear combination of approximate Schur
164
 
c                      vectors associated with the "wanted" Ritz values.
165
 
c          -------------------------------------------------------------
166
 
c
167
 
c          IPARAM(2) = No longer referenced.
168
 
c
169
 
c          IPARAM(3) = MXITER
170
 
c          On INPUT:  maximum number of Arnoldi update iterations allowed. 
171
 
c          On OUTPUT: actual number of Arnoldi update iterations taken. 
172
 
c
173
 
c          IPARAM(4) = NB: blocksize to be used in the recurrence.
174
 
c          The code currently works only for NB = 1.
175
 
c
176
 
c          IPARAM(5) = NCONV: number of "converged" Ritz values.
177
 
c          This represents the number of Ritz values that satisfy
178
 
c          the convergence criterion.
179
 
c
180
 
c          IPARAM(6) = IUPD
181
 
c          No longer referenced. Implicit restarting is ALWAYS used.  
182
 
c
183
 
c          IPARAM(7) = MODE
184
 
c          On INPUT determines what type of eigenproblem is being solved.
185
 
c          Must be 1,2,3,4; See under \Description of snaupd for the 
186
 
c          four modes available.
187
 
c
188
 
c          IPARAM(8) = NP
189
 
c          When ido = 3 and the user provides shifts through reverse
190
 
c          communication (IPARAM(1)=0), snaupd returns NP, the number
191
 
c          of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark
192
 
c          5 below.
193
 
c
194
 
c          IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO,
195
 
c          OUTPUT: NUMOP  = total number of OP*x operations,
196
 
c                  NUMOPB = total number of B*x operations if BMAT='G',
197
 
c                  NUMREO = total number of steps of re-orthogonalization.        
198
 
c
199
 
c  IPNTR   Integer array of length 14.  (OUTPUT)
200
 
c          Pointer to mark the starting locations in the WORKD and WORKL
201
 
c          arrays for matrices/vectors used by the Arnoldi iteration.
202
 
c          -------------------------------------------------------------
203
 
c          IPNTR(1): pointer to the current operand vector X in WORKD.
204
 
c          IPNTR(2): pointer to the current result vector Y in WORKD.
205
 
c          IPNTR(3): pointer to the vector B * X in WORKD when used in 
206
 
c                    the shift-and-invert mode.
207
 
c          IPNTR(4): pointer to the next available location in WORKL
208
 
c                    that is untouched by the program.
209
 
c          IPNTR(5): pointer to the NCV by NCV upper Hessenberg matrix
210
 
c                    H in WORKL.
211
 
c          IPNTR(6): pointer to the real part of the ritz value array 
212
 
c                    RITZR in WORKL.
213
 
c          IPNTR(7): pointer to the imaginary part of the ritz value array
214
 
c                    RITZI in WORKL.
215
 
c          IPNTR(8): pointer to the Ritz estimates in array WORKL associated
216
 
c                    with the Ritz values located in RITZR and RITZI in WORKL.
217
 
c
218
 
c          IPNTR(14): pointer to the NP shifts in WORKL. See Remark 5 below.
219
 
c
220
 
c          Note: IPNTR(9:13) is only referenced by sneupd. See Remark 2 below.
221
 
c
222
 
c          IPNTR(9):  pointer to the real part of the NCV RITZ values of the 
223
 
c                     original system.
224
 
c          IPNTR(10): pointer to the imaginary part of the NCV RITZ values of 
225
 
c                     the original system.
226
 
c          IPNTR(11): pointer to the NCV corresponding error bounds.
227
 
c          IPNTR(12): pointer to the NCV by NCV upper quasi-triangular
228
 
c                     Schur matrix for H.
229
 
c          IPNTR(13): pointer to the NCV by NCV matrix of eigenvectors
230
 
c                     of the upper Hessenberg matrix H. Only referenced by
231
 
c                     sneupd if RVEC = .TRUE. See Remark 2 below.
232
 
c          -------------------------------------------------------------
233
 
c          
234
 
c  WORKD   Real work array of length 3*N.  (REVERSE COMMUNICATION)
235
 
c          Distributed array to be used in the basic Arnoldi iteration
236
 
c          for reverse communication.  The user should not use WORKD 
237
 
c          as temporary workspace during the iteration. Upon termination
238
 
c          WORKD(1:N) contains B*RESID(1:N). If an invariant subspace
239
 
c          associated with the converged Ritz values is desired, see remark
240
 
c          2 below, subroutine sneupd uses this output.
241
 
c          See Data Distribution Note below.  
242
 
c
243
 
c  WORKL   Real work array of length LWORKL.  (OUTPUT/WORKSPACE)
244
 
c          Private (replicated) array on each PE or array allocated on
245
 
c          the front end.  See Data Distribution Note below.
246
 
c
247
 
c  LWORKL  Integer.  (INPUT)
248
 
c          LWORKL must be at least 3*NCV**2 + 6*NCV.
249
 
c
250
 
c  INFO    Integer.  (INPUT/OUTPUT)
251
 
c          If INFO .EQ. 0, a randomly initial residual vector is used.
252
 
c          If INFO .NE. 0, RESID contains the initial residual vector,
253
 
c                          possibly from a previous run.
254
 
c          Error flag on output.
255
 
c          =  0: Normal exit.
256
 
c          =  1: Maximum number of iterations taken.
257
 
c                All possible eigenvalues of OP has been found. IPARAM(5)  
258
 
c                returns the number of wanted converged Ritz values.
259
 
c          =  2: No longer an informational error. Deprecated starting
260
 
c                with release 2 of ARPACK.
261
 
c          =  3: No shifts could be applied during a cycle of the 
262
 
c                Implicitly restarted Arnoldi iteration. One possibility 
263
 
c                is to increase the size of NCV relative to NEV. 
264
 
c                See remark 4 below.
265
 
c          = -1: N must be positive.
266
 
c          = -2: NEV must be positive.
267
 
c          = -3: NCV-NEV >= 2 and less than or equal to N.
268
 
c          = -4: The maximum number of Arnoldi update iteration 
269
 
c                must be greater than zero.
270
 
c          = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'
271
 
c          = -6: BMAT must be one of 'I' or 'G'.
272
 
c          = -7: Length of private work array is not sufficient.
273
 
c          = -8: Error return from LAPACK eigenvalue calculation;
274
 
c          = -9: Starting vector is zero.
275
 
c          = -10: IPARAM(7) must be 1,2,3,4.
276
 
c          = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable.
277
 
c          = -12: IPARAM(1) must be equal to 0 or 1.
278
 
c          = -9999: Could not build an Arnoldi factorization.
279
 
c                   IPARAM(5) returns the size of the current Arnoldi
280
 
c                   factorization.
281
 
c
282
 
c\Remarks
283
 
c  1. The computed Ritz values are approximate eigenvalues of OP. The
284
 
c     selection of WHICH should be made with this in mind when
285
 
c     Mode = 3 and 4.  After convergence, approximate eigenvalues of the
286
 
c     original problem may be obtained with the ARPACK subroutine sneupd.
287
 
c
288
 
c  2. If a basis for the invariant subspace corresponding to the converged Ritz 
289
 
c     values is needed, the user must call sneupd immediately following 
290
 
c     completion of snaupd. This is new starting with release 2 of ARPACK.
291
 
c
292
 
c  3. If M can be factored into a Cholesky factorization M = LL`
293
 
c     then Mode = 2 should not be selected.  Instead one should use
294
 
c     Mode = 1 with  OP = inv(L)*A*inv(L`).  Appropriate triangular 
295
 
c     linear systems should be solved with L and L` rather
296
 
c     than computing inverses.  After convergence, an approximate
297
 
c     eigenvector z of the original problem is recovered by solving
298
 
c     L`z = x  where x is a Ritz vector of OP.
299
 
c
300
 
c  4. At present there is no a-priori analysis to guide the selection
301
 
c     of NCV relative to NEV.  The only formal requrement is that NCV > NEV + 2.
302
 
c     However, it is recommended that NCV .ge. 2*NEV+1.  If many problems of
303
 
c     the same type are to be solved, one should experiment with increasing
304
 
c     NCV while keeping NEV fixed for a given test problem.  This will 
305
 
c     usually decrease the required number of OP*x operations but it
306
 
c     also increases the work and storage required to maintain the orthogonal
307
 
c     basis vectors.  The optimal "cross-over" with respect to CPU time
308
 
c     is problem dependent and must be determined empirically. 
309
 
c     See Chapter 8 of Reference 2 for further information.
310
 
c
311
 
c  5. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the 
312
 
c     NP = IPARAM(8) real and imaginary parts of the shifts in locations 
313
 
c         real part                  imaginary part
314
 
c         -----------------------    --------------
315
 
c     1   WORKL(IPNTR(14))           WORKL(IPNTR(14)+NP)
316
 
c     2   WORKL(IPNTR(14)+1)         WORKL(IPNTR(14)+NP+1)
317
 
c                        .                          .
318
 
c                        .                          .
319
 
c                        .                          .
320
 
c     NP  WORKL(IPNTR(14)+NP-1)      WORKL(IPNTR(14)+2*NP-1).
321
 
c
322
 
c     Only complex conjugate pairs of shifts may be applied and the pairs 
323
 
c     must be placed in consecutive locations. The real part of the 
324
 
c     eigenvalues of the current upper Hessenberg matrix are located in 
325
 
c     WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1) and the imaginary part 
326
 
c     in WORKL(IPNTR(7)) through WORKL(IPNTR(7)+NCV-1). They are ordered
327
 
c     according to the order defined by WHICH. The complex conjugate
328
 
c     pairs are kept together and the associated Ritz estimates are located in
329
 
c     WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1).
330
 
c
331
 
c-----------------------------------------------------------------------
332
 
c
333
 
c\Data Distribution Note: 
334
 
c
335
 
c  Fortran-D syntax:
336
 
c  ================
337
 
c  Real resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
338
 
c  decompose  d1(n), d2(n,ncv)
339
 
c  align      resid(i) with d1(i)
340
 
c  align      v(i,j)   with d2(i,j)
341
 
c  align      workd(i) with d1(i)     range (1:n)
342
 
c  align      workd(i) with d1(i-n)   range (n+1:2*n)
343
 
c  align      workd(i) with d1(i-2*n) range (2*n+1:3*n)
344
 
c  distribute d1(block), d2(block,:)
345
 
c  replicated workl(lworkl)
346
 
c
347
 
c  Cray MPP syntax:
348
 
c  ===============
349
 
c  Real  resid(n), v(ldv,ncv), workd(n,3), workl(lworkl)
350
 
c  shared     resid(block), v(block,:), workd(block,:)
351
 
c  replicated workl(lworkl)
352
 
c  
353
 
c  CM2/CM5 syntax:
354
 
c  ==============
355
 
c  
356
 
c-----------------------------------------------------------------------
357
 
c
358
 
c     include   'ex-nonsym.doc'
359
 
c
360
 
c-----------------------------------------------------------------------
361
 
c
362
 
c\BeginLib
363
 
c
364
 
c\Local variables:
365
 
c     xxxxxx  real
366
 
c
367
 
c\References:
368
 
c  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
369
 
c     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
370
 
c     pp 357-385.
371
 
c  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly 
372
 
c     Restarted Arnoldi Iteration", Rice University Technical Report
373
 
c     TR95-13, Department of Computational and Applied Mathematics.
374
 
c  3. B.N. Parlett & Y. Saad, "Complex Shift and Invert Strategies for
375
 
c     Real Matrices", Linear Algebra and its Applications, vol 88/89,
376
 
c     pp 575-595, (1987).
377
 
c
378
 
c\Routines called:
379
 
c     snaup2  ARPACK routine that implements the Implicitly Restarted
380
 
c             Arnoldi Iteration.
381
 
c     ivout   ARPACK utility routine that prints integers.
382
 
c     second  ARPACK utility routine for timing.
383
 
c     svout   ARPACK utility routine that prints vectors.
384
 
c     slamch  LAPACK routine that determines machine constants.
385
 
c
386
 
c\Author
387
 
c     Danny Sorensen               Phuong Vu
388
 
c     Richard Lehoucq              CRPC / Rice University
389
 
c     Dept. of Computational &     Houston, Texas
390
 
c     Applied Mathematics
391
 
c     Rice University           
392
 
c     Houston, Texas            
393
 
394
 
c\Revision history:
395
 
c     12/16/93: Version '1.1'
396
 
c
397
 
c\SCCS Information: @(#) 
398
 
c FILE: naupd.F   SID: 2.10   DATE OF SID: 08/23/02   RELEASE: 2
399
 
c
400
 
c\Remarks
401
 
c
402
 
c\EndLib
403
 
c
404
 
c-----------------------------------------------------------------------
405
 
c
406
 
      subroutine snaupd
407
 
     &   ( ido, bmat, n, which, nev, tol, resid, ncv, v, ldv, iparam, 
408
 
     &     ipntr, workd, workl, lworkl, info )
409
 
c
410
 
c     %----------------------------------------------------%
411
 
c     | Include files for debugging and timing information |
412
 
c     %----------------------------------------------------%
413
 
c
414
 
      include   'debug.h'
415
 
      include   'stat.h'
416
 
c
417
 
c     %------------------%
418
 
c     | Scalar Arguments |
419
 
c     %------------------%
420
 
c
421
 
      character  bmat*1, which*2
422
 
      integer    ido, info, ldv, lworkl, n, ncv, nev
423
 
      Real
424
 
     &           tol
425
 
c
426
 
c     %-----------------%
427
 
c     | Array Arguments |
428
 
c     %-----------------%
429
 
c
430
 
      integer    iparam(11), ipntr(14)
431
 
      Real
432
 
     &           resid(n), v(ldv,ncv), workd(3*n), workl(lworkl)
433
 
c
434
 
c     %------------%
435
 
c     | Parameters |
436
 
c     %------------%
437
 
c
438
 
      Real
439
 
     &           one, zero
440
 
      parameter (one = 1.0E+0, zero = 0.0E+0)
441
 
c
442
 
c     %---------------%
443
 
c     | Local Scalars |
444
 
c     %---------------%
445
 
c
446
 
      integer    bounds, ierr, ih, iq, ishift, iupd, iw, 
447
 
     &           ldh, ldq, levec, mode, msglvl, mxiter, nb,
448
 
     &           nev0, next, np, ritzi, ritzr, j
449
 
      save       bounds, ih, iq, ishift, iupd, iw, ldh, ldq,
450
 
     &           levec, mode, msglvl, mxiter, nb, nev0, next,
451
 
     &           np, ritzi, ritzr
452
 
c
453
 
c     %----------------------%
454
 
c     | External Subroutines |
455
 
c     %----------------------%
456
 
c
457
 
      external   snaup2, svout, ivout, second, sstatn
458
 
c
459
 
c     %--------------------%
460
 
c     | External Functions |
461
 
c     %--------------------%
462
 
c
463
 
      Real
464
 
     &           slamch
465
 
      external   slamch
466
 
c
467
 
c     %-----------------------%
468
 
c     | Executable Statements |
469
 
c     %-----------------------%
470
 
471
 
      if (ido .eq. 0) then
472
 
473
 
c        %-------------------------------%
474
 
c        | Initialize timing statistics  |
475
 
c        | & message level for debugging |
476
 
c        %-------------------------------%
477
 
c
478
 
         call sstatn
479
 
         call second (t0)
480
 
         msglvl = mnaupd
481
 
c
482
 
c        %----------------%
483
 
c        | Error checking |
484
 
c        %----------------%
485
 
c
486
 
         ierr   = 0
487
 
         ishift = iparam(1)
488
 
c         levec  = iparam(2)
489
 
         mxiter = iparam(3)
490
 
c         nb     = iparam(4)
491
 
         nb     = 1
492
 
c
493
 
c        %--------------------------------------------%
494
 
c        | Revision 2 performs only implicit restart. |
495
 
c        %--------------------------------------------%
496
 
c
497
 
         iupd   = 1
498
 
         mode   = iparam(7)
499
 
c
500
 
         if (n .le. 0) then
501
 
            ierr = -1
502
 
         else if (nev .le. 0) then
503
 
            ierr = -2
504
 
         else if (ncv .le. nev+1 .or.  ncv .gt. n) then
505
 
            ierr = -3
506
 
         else if (mxiter .le.          0) then
507
 
            ierr = 4
508
 
         else if (which .ne. 'LM' .and.
509
 
     &       which .ne. 'SM' .and.
510
 
     &       which .ne. 'LR' .and.
511
 
     &       which .ne. 'SR' .and.
512
 
     &       which .ne. 'LI' .and.
513
 
     &       which .ne. 'SI') then
514
 
            ierr = -5
515
 
         else if (bmat .ne. 'I' .and. bmat .ne. 'G') then
516
 
            ierr = -6
517
 
         else if (lworkl .lt. 3*ncv**2 + 6*ncv) then
518
 
            ierr = -7
519
 
         else if (mode .lt. 1 .or. mode .gt. 4) then
520
 
            ierr = -10
521
 
         else if (mode .eq. 1 .and. bmat .eq. 'G') then
522
 
            ierr = -11
523
 
         else if (ishift .lt. 0 .or. ishift .gt. 1) then
524
 
            ierr = -12
525
 
         end if
526
 
527
 
c        %------------%
528
 
c        | Error Exit |
529
 
c        %------------%
530
 
c
531
 
         if (ierr .ne. 0) then
532
 
            info = ierr
533
 
            ido  = 99
534
 
            go to 9000
535
 
         end if
536
 
537
 
c        %------------------------%
538
 
c        | Set default parameters |
539
 
c        %------------------------%
540
 
c
541
 
         if (nb .le. 0)                         nb = 1
542
 
         if (tol .le. zero)                     tol = slamch('EpsMach')
543
 
c
544
 
c        %----------------------------------------------%
545
 
c        | NP is the number of additional steps to      |
546
 
c        | extend the length NEV Lanczos factorization. |
547
 
c        | NEV0 is the local variable designating the   |
548
 
c        | size of the invariant subspace desired.      |
549
 
c        %----------------------------------------------%
550
 
c
551
 
         np     = ncv - nev
552
 
         nev0   = nev 
553
 
554
 
c        %-----------------------------%
555
 
c        | Zero out internal workspace |
556
 
c        %-----------------------------%
557
 
c
558
 
         do 10 j = 1, 3*ncv**2 + 6*ncv
559
 
            workl(j) = zero
560
 
  10     continue
561
 
562
 
c        %-------------------------------------------------------------%
563
 
c        | Pointer into WORKL for address of H, RITZ, BOUNDS, Q        |
564
 
c        | etc... and the remaining workspace.                         |
565
 
c        | Also update pointer to be used on output.                   |
566
 
c        | Memory is laid out as follows:                              |
567
 
c        | workl(1:ncv*ncv) := generated Hessenberg matrix             |
568
 
c        | workl(ncv*ncv+1:ncv*ncv+2*ncv) := real and imaginary        |
569
 
c        |                                   parts of ritz values      |
570
 
c        | workl(ncv*ncv+2*ncv+1:ncv*ncv+3*ncv) := error bounds        |
571
 
c        | workl(ncv*ncv+3*ncv+1:2*ncv*ncv+3*ncv) := rotation matrix Q |
572
 
c        | workl(2*ncv*ncv+3*ncv+1:3*ncv*ncv+6*ncv) := workspace       |
573
 
c        | The final workspace is needed by subroutine sneigh called   |
574
 
c        | by snaup2. Subroutine sneigh calls LAPACK routines for      |
575
 
c        | calculating eigenvalues and the last row of the eigenvector |
576
 
c        | matrix.                                                     |
577
 
c        %-------------------------------------------------------------%
578
 
c
579
 
         ldh    = ncv
580
 
         ldq    = ncv
581
 
         ih     = 1
582
 
         ritzr  = ih     + ldh*ncv
583
 
         ritzi  = ritzr  + ncv
584
 
         bounds = ritzi  + ncv
585
 
         iq     = bounds + ncv
586
 
         iw     = iq     + ldq*ncv
587
 
         next   = iw     + ncv**2 + 3*ncv
588
 
c
589
 
         ipntr(4) = next
590
 
         ipntr(5) = ih
591
 
         ipntr(6) = ritzr
592
 
         ipntr(7) = ritzi
593
 
         ipntr(8) = bounds
594
 
         ipntr(14) = iw 
595
 
c
596
 
      end if
597
 
c
598
 
c     %-------------------------------------------------------%
599
 
c     | Carry out the Implicitly restarted Arnoldi Iteration. |
600
 
c     %-------------------------------------------------------%
601
 
c
602
 
      call snaup2 
603
 
     &   ( ido, bmat, n, which, nev0, np, tol, resid, mode, iupd,
604
 
     &     ishift, mxiter, v, ldv, workl(ih), ldh, workl(ritzr), 
605
 
     &     workl(ritzi), workl(bounds), workl(iq), ldq, workl(iw), 
606
 
     &     ipntr, workd, info )
607
 
608
 
c     %--------------------------------------------------%
609
 
c     | ido .ne. 99 implies use of reverse communication |
610
 
c     | to compute operations involving OP or shifts.    |
611
 
c     %--------------------------------------------------%
612
 
c
613
 
      if (ido .eq. 3) iparam(8) = np
614
 
      if (ido .ne. 99) go to 9000
615
 
616
 
      iparam(3) = mxiter
617
 
      iparam(5) = np
618
 
      iparam(9) = nopx
619
 
      iparam(10) = nbx
620
 
      iparam(11) = nrorth
621
 
c
622
 
c     %------------------------------------%
623
 
c     | Exit if there was an informational |
624
 
c     | error within snaup2.               |
625
 
c     %------------------------------------%
626
 
c
627
 
      if (info .lt. 0) go to 9000
628
 
      if (info .eq. 2) info = 3
629
 
c
630
 
      if (msglvl .gt. 0) then
631
 
         call ivout (logfil, 1, mxiter, ndigit,
632
 
     &               '_naupd: Number of update iterations taken')
633
 
         call ivout (logfil, 1, np, ndigit,
634
 
     &               '_naupd: Number of wanted "converged" Ritz values')
635
 
         call svout (logfil, np, workl(ritzr), ndigit, 
636
 
     &               '_naupd: Real part of the final Ritz values')
637
 
         call svout (logfil, np, workl(ritzi), ndigit, 
638
 
     &               '_naupd: Imaginary part of the final Ritz values')
639
 
         call svout (logfil, np, workl(bounds), ndigit, 
640
 
     &               '_naupd: Associated Ritz estimates')
641
 
      end if
642
 
c
643
 
      call second (t1)
644
 
      tnaupd = t1 - t0
645
 
c
646
 
      if (msglvl .gt. 0) then
647
 
c
648
 
c        %--------------------------------------------------------%
649
 
c        | Version Number & Version Date are defined in version.h |
650
 
c        %--------------------------------------------------------%
651
 
c
652
 
         write (6,1000)
653
 
         write (6,1100) mxiter, nopx, nbx, nrorth, nitref, nrstrt,
654
 
     &                  tmvopx, tmvbx, tnaupd, tnaup2, tnaitr, titref,
655
 
     &                  tgetv0, tneigh, tngets, tnapps, tnconv, trvec
656
 
 1000    format (//,
657
 
     &      5x, '=============================================',/
658
 
     &      5x, '= Nonsymmetric implicit Arnoldi update code =',/
659
 
     &      5x, '= Version Number: ', ' 2.4', 21x, ' =',/
660
 
     &      5x, '= Version Date:   ', ' 07/31/96', 16x,   ' =',/
661
 
     &      5x, '=============================================',/
662
 
     &      5x, '= Summary of timing statistics              =',/
663
 
     &      5x, '=============================================',//)
664
 
 1100    format (
665
 
     &      5x, 'Total number update iterations             = ', i5,/
666
 
     &      5x, 'Total number of OP*x operations            = ', i5,/
667
 
     &      5x, 'Total number of B*x operations             = ', i5,/
668
 
     &      5x, 'Total number of reorthogonalization steps  = ', i5,/
669
 
     &      5x, 'Total number of iterative refinement steps = ', i5,/
670
 
     &      5x, 'Total number of restart steps              = ', i5,/
671
 
     &      5x, 'Total time in user OP*x operation          = ', f12.6,/
672
 
     &      5x, 'Total time in user B*x operation           = ', f12.6,/
673
 
     &      5x, 'Total time in Arnoldi update routine       = ', f12.6,/
674
 
     &      5x, 'Total time in naup2 routine                = ', f12.6,/
675
 
     &      5x, 'Total time in basic Arnoldi iteration loop = ', f12.6,/
676
 
     &      5x, 'Total time in reorthogonalization phase    = ', f12.6,/
677
 
     &      5x, 'Total time in (re)start vector generation  = ', f12.6,/
678
 
     &      5x, 'Total time in Hessenberg eig. subproblem   = ', f12.6,/
679
 
     &      5x, 'Total time in getting the shifts           = ', f12.6,/
680
 
     &      5x, 'Total time in applying the shifts          = ', f12.6,/
681
 
     &      5x, 'Total time in convergence testing          = ', f12.6,/
682
 
     &      5x, 'Total time in computing final Ritz vectors = ', f12.6/)
683
 
      end if
684
 
c
685
 
 9000 continue
686
 
c
687
 
      return
688
 
c
689
 
c     %---------------%
690
 
c     | End of snaupd |
691
 
c     %---------------%
692
 
c
693
 
      end