~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/arpack/dsaupd.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

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