← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/arpack/ARPACK/SRC/snaup2.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:
Lib/sandbox/arpack/ARPACK/SRC/snaup2.f
1
 
c\BeginDoc
2
 
c
3
 
c\Name: snaup2
4
 
c
5
 
c\Description: 
6
 
c  Intermediate level interface called by snaupd.
7
 
c
8
 
c\Usage:
9
 
c  call snaup2
10
 
c     ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD,
11
 
c       ISHIFT, MXITER, V, LDV, H, LDH, RITZR, RITZI, BOUNDS, 
12
 
c       Q, LDQ, WORKL, IPNTR, WORKD, INFO )
13
 
c
14
 
c\Arguments
15
 
c
16
 
c  IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in snaupd.
17
 
c  MODE, ISHIFT, MXITER: see the definition of IPARAM in snaupd.
18
 
c
19
 
c  NP      Integer.  (INPUT/OUTPUT)
20
 
c          Contains the number of implicit shifts to apply during 
21
 
c          each Arnoldi iteration.  
22
 
c          If ISHIFT=1, NP is adjusted dynamically at each iteration 
23
 
c          to accelerate convergence and prevent stagnation.
24
 
c          This is also roughly equal to the number of matrix-vector 
25
 
c          products (involving the operator OP) per Arnoldi iteration.
26
 
c          The logic for adjusting is contained within the current
27
 
c          subroutine.
28
 
c          If ISHIFT=0, NP is the number of shifts the user needs
29
 
c          to provide via reverse comunication. 0 < NP < NCV-NEV.
30
 
c          NP may be less than NCV-NEV for two reasons. The first, is
31
 
c          to keep complex conjugate pairs of "wanted" Ritz values 
32
 
c          together. The second, is that a leading block of the current
33
 
c          upper Hessenberg matrix has split off and contains "unwanted"
34
 
c          Ritz values.
35
 
c          Upon termination of the IRA iteration, NP contains the number 
36
 
c          of "converged" wanted Ritz values.
37
 
c
38
 
c  IUPD    Integer.  (INPUT)
39
 
c          IUPD .EQ. 0: use explicit restart instead implicit update.
40
 
c          IUPD .NE. 0: use implicit update.
41
 
c
42
 
c  V       Real N by (NEV+NP) array.  (INPUT/OUTPUT)
43
 
c          The Arnoldi basis vectors are returned in the first NEV 
44
 
c          columns of V.
45
 
c
46
 
c  LDV     Integer.  (INPUT)
47
 
c          Leading dimension of V exactly as declared in the calling 
48
 
c          program.
49
 
c
50
 
c  H       Real (NEV+NP) by (NEV+NP) array.  (OUTPUT)
51
 
c          H is used to store the generated upper Hessenberg matrix
52
 
c
53
 
c  LDH     Integer.  (INPUT)
54
 
c          Leading dimension of H exactly as declared in the calling 
55
 
c          program.
56
 
c
57
 
c  RITZR,  Real arrays of length NEV+NP.  (OUTPUT)
58
 
c  RITZI   RITZR(1:NEV) (resp. RITZI(1:NEV)) contains the real (resp.
59
 
c          imaginary) part of the computed Ritz values of OP.
60
 
c
61
 
c  BOUNDS  Real array of length NEV+NP.  (OUTPUT)
62
 
c          BOUNDS(1:NEV) contain the error bounds corresponding to 
63
 
c          the computed Ritz values.
64
 
c          
65
 
c  Q       Real (NEV+NP) by (NEV+NP) array.  (WORKSPACE)
66
 
c          Private (replicated) work array used to accumulate the
67
 
c          rotation in the shift application step.
68
 
c
69
 
c  LDQ     Integer.  (INPUT)
70
 
c          Leading dimension of Q exactly as declared in the calling
71
 
c          program.
72
 
c
73
 
c  WORKL   Real work array of length at least 
74
 
c          (NEV+NP)**2 + 3*(NEV+NP).  (INPUT/WORKSPACE)
75
 
c          Private (replicated) array on each PE or array allocated on
76
 
c          the front end.  It is used in shifts calculation, shifts
77
 
c          application and convergence checking.
78
 
c
79
 
c          On exit, the last 3*(NEV+NP) locations of WORKL contain
80
 
c          the Ritz values (real,imaginary) and associated Ritz
81
 
c          estimates of the current Hessenberg matrix.  They are
82
 
c          listed in the same order as returned from sneigh.
83
 
c
84
 
c          If ISHIFT .EQ. O and IDO .EQ. 3, the first 2*NP locations
85
 
c          of WORKL are used in reverse communication to hold the user 
86
 
c          supplied shifts.
87
 
c
88
 
c  IPNTR   Integer array of length 3.  (OUTPUT)
89
 
c          Pointer to mark the starting locations in the WORKD for 
90
 
c          vectors used by the Arnoldi iteration.
91
 
c          -------------------------------------------------------------
92
 
c          IPNTR(1): pointer to the current operand vector X.
93
 
c          IPNTR(2): pointer to the current result vector Y.
94
 
c          IPNTR(3): pointer to the vector B * X when used in the 
95
 
c                    shift-and-invert mode.  X is the current operand.
96
 
c          -------------------------------------------------------------
97
 
c          
98
 
c  WORKD   Real work array of length 3*N.  (WORKSPACE)
99
 
c          Distributed array to be used in the basic Arnoldi iteration
100
 
c          for reverse communication.  The user should not use WORKD
101
 
c          as temporary workspace during the iteration !!!!!!!!!!
102
 
c          See Data Distribution Note in DNAUPD.
103
 
c
104
 
c  INFO    Integer.  (INPUT/OUTPUT)
105
 
c          If INFO .EQ. 0, a randomly initial residual vector is used.
106
 
c          If INFO .NE. 0, RESID contains the initial residual vector,
107
 
c                          possibly from a previous run.
108
 
c          Error flag on output.
109
 
c          =     0: Normal return.
110
 
c          =     1: Maximum number of iterations taken.
111
 
c                   All possible eigenvalues of OP has been found.  
112
 
c                   NP returns the number of converged Ritz values.
113
 
c          =     2: No shifts could be applied.
114
 
c          =    -8: Error return from LAPACK eigenvalue calculation;
115
 
c                   This should never happen.
116
 
c          =    -9: Starting vector is zero.
117
 
c          = -9999: Could not build an Arnoldi factorization.
118
 
c                   Size that was built in returned in NP.
119
 
c
120
 
c\EndDoc
121
 
c
122
 
c-----------------------------------------------------------------------
123
 
c
124
 
c\BeginLib
125
 
c
126
 
c\Local variables:
127
 
c     xxxxxx  real
128
 
c
129
 
c\References:
130
 
c  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
131
 
c     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
132
 
c     pp 357-385.
133
 
c  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly 
134
 
c     Restarted Arnoldi Iteration", Rice University Technical Report
135
 
c     TR95-13, Department of Computational and Applied Mathematics.
136
 
c
137
 
c\Routines called:
138
 
c     sgetv0  ARPACK initial vector generation routine. 
139
 
c     snaitr  ARPACK Arnoldi factorization routine.
140
 
c     snapps  ARPACK application of implicit shifts routine.
141
 
c     snconv  ARPACK convergence of Ritz values routine.
142
 
c     sneigh  ARPACK compute Ritz values and error bounds routine.
143
 
c     sngets  ARPACK reorder Ritz values and error bounds routine.
144
 
c     ssortc  ARPACK sorting routine.
145
 
c     ivout   ARPACK utility routine that prints integers.
146
 
c     second  ARPACK utility routine for timing.
147
 
c     smout   ARPACK utility routine that prints matrices
148
 
c     svout   ARPACK utility routine that prints vectors.
149
 
c     slamch  LAPACK routine that determines machine constants.
150
 
c     slapy2  LAPACK routine to compute sqrt(x**2+y**2) carefully.
151
 
c     scopy   Level 1 BLAS that copies one vector to another .
152
 
c     sdot    Level 1 BLAS that computes the scalar product of two vectors. 
153
 
c     snrm2   Level 1 BLAS that computes the norm of a vector.
154
 
c     sswap   Level 1 BLAS that swaps two vectors.
155
 
c
156
 
c\Author
157
 
c     Danny Sorensen               Phuong Vu
158
 
c     Richard Lehoucq              CRPC / Rice University 
159
 
c     Dept. of Computational &     Houston, Texas 
160
 
c     Applied Mathematics
161
 
c     Rice University           
162
 
c     Houston, Texas    
163
 
c 
164
 
c\SCCS Information: @(#) 
165
 
c FILE: naup2.F   SID: 2.8   DATE OF SID: 10/17/00   RELEASE: 2
166
 
c
167
 
c\Remarks
168
 
c     1. None
169
 
c
170
 
c\EndLib
171
 
c
172
 
c-----------------------------------------------------------------------
173
 
c
174
 
      subroutine snaup2
175
 
     &   ( ido, bmat, n, which, nev, np, tol, resid, mode, iupd, 
176
 
     &     ishift, mxiter, v, ldv, h, ldh, ritzr, ritzi, bounds, 
177
 
     &     q, ldq, workl, ipntr, workd, info )
178
 
c
179
 
c     %----------------------------------------------------%
180
 
c     | Include files for debugging and timing information |
181
 
c     %----------------------------------------------------%
182
 
c
183
 
      include   'debug.h'
184
 
      include   'stat.h'
185
 
c
186
 
c     %------------------%
187
 
c     | Scalar Arguments |
188
 
c     %------------------%
189
 
c
190
 
      character  bmat*1, which*2
191
 
      integer    ido, info, ishift, iupd, mode, ldh, ldq, ldv, mxiter,
192
 
     &           n, nev, np
193
 
      Real
194
 
     &           tol
195
 
c
196
 
c     %-----------------%
197
 
c     | Array Arguments |
198
 
c     %-----------------%
199
 
c
200
 
      integer    ipntr(13)
201
 
      Real
202
 
     &           bounds(nev+np), h(ldh,nev+np), q(ldq,nev+np), resid(n),
203
 
     &           ritzi(nev+np), ritzr(nev+np), v(ldv,nev+np), 
204
 
     &           workd(3*n), workl( (nev+np)*(nev+np+3) )
205
 
c
206
 
c     %------------%
207
 
c     | Parameters |
208
 
c     %------------%
209
 
c
210
 
      Real
211
 
     &           one, zero
212
 
      parameter (one = 1.0E+0, zero = 0.0E+0)
213
 
c
214
 
c     %---------------%
215
 
c     | Local Scalars |
216
 
c     %---------------%
217
 
c
218
 
      character  wprime*2
219
 
      logical    cnorm , getv0, initv, update, ushift
220
 
      integer    ierr  , iter , j    , kplusp, msglvl, nconv, 
221
 
     &           nevbef, nev0 , np0  , nptemp, numcnv
222
 
      Real
223
 
     &           rnorm , temp , eps23
224
 
      save       cnorm , getv0, initv, update, ushift,
225
 
     &           rnorm , iter , eps23, kplusp, msglvl, nconv , 
226
 
     &           nevbef, nev0 , np0  , numcnv
227
 
c
228
 
c     %-----------------------%
229
 
c     | Local array arguments |
230
 
c     %-----------------------%
231
 
c
232
 
      integer    kp(4)
233
 
c
234
 
c     %----------------------%
235
 
c     | External Subroutines |
236
 
c     %----------------------%
237
 
c
238
 
      external   scopy , sgetv0, snaitr, snconv, sneigh, 
239
 
     &           sngets, snapps, svout , ivout , second
240
 
c
241
 
c     %--------------------%
242
 
c     | External Functions |
243
 
c     %--------------------%
244
 
c
245
 
      Real
246
 
     &           sdot, snrm2, slapy2, slamch
247
 
      external   sdot, snrm2, slapy2, slamch
248
 
c
249
 
c     %---------------------%
250
 
c     | Intrinsic Functions |
251
 
c     %---------------------%
252
 
c
253
 
      intrinsic    min, max, abs, sqrt
254
 
c
255
 
c     %-----------------------%
256
 
c     | Executable Statements |
257
 
c     %-----------------------%
258
 
c
259
 
      if (ido .eq. 0) then
260
 
c 
261
 
         call second (t0)
262
 
c 
263
 
         msglvl = mnaup2
264
 
c 
265
 
c        %-------------------------------------%
266
 
c        | Get the machine dependent constant. |
267
 
c        %-------------------------------------%
268
 
c
269
 
         eps23 = slamch('Epsilon-Machine')
270
 
         eps23 = eps23**(2.0E+0 / 3.0E+0)
271
 
c
272
 
         nev0   = nev
273
 
         np0    = np
274
 
c
275
 
c        %-------------------------------------%
276
 
c        | kplusp is the bound on the largest  |
277
 
c        |        Lanczos factorization built. |
278
 
c        | nconv is the current number of      |
279
 
c        |        "converged" eigenvlues.      |
280
 
c        | iter is the counter on the current  |
281
 
c        |      iteration step.                |
282
 
c        %-------------------------------------%
283
 
c
284
 
         kplusp = nev + np
285
 
         nconv  = 0
286
 
         iter   = 0
287
 
c 
288
 
c        %---------------------------------------%
289
 
c        | Set flags for computing the first NEV |
290
 
c        | steps of the Arnoldi factorization.   |
291
 
c        %---------------------------------------%
292
 
c
293
 
         getv0    = .true.
294
 
         update   = .false.
295
 
         ushift   = .false.
296
 
         cnorm    = .false.
297
 
c
298
 
         if (info .ne. 0) then
299
 
c
300
 
c           %--------------------------------------------%
301
 
c           | User provides the initial residual vector. |
302
 
c           %--------------------------------------------%
303
 
c
304
 
            initv = .true.
305
 
            info  = 0
306
 
         else
307
 
            initv = .false.
308
 
         end if
309
 
      end if
310
 
c 
311
 
c     %---------------------------------------------%
312
 
c     | Get a possibly random starting vector and   |
313
 
c     | force it into the range of the operator OP. |
314
 
c     %---------------------------------------------%
315
 
c
316
 
   10 continue
317
 
c
318
 
      if (getv0) then
319
 
         call sgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
320
 
     &                ipntr, workd, info)
321
 
c
322
 
         if (ido .ne. 99) go to 9000
323
 
c
324
 
         if (rnorm .eq. zero) then
325
 
c
326
 
c           %-----------------------------------------%
327
 
c           | The initial vector is zero. Error exit. | 
328
 
c           %-----------------------------------------%
329
 
c
330
 
            info = -9
331
 
            go to 1100
332
 
         end if
333
 
         getv0 = .false.
334
 
         ido  = 0
335
 
      end if
336
 
c 
337
 
c     %-----------------------------------%
338
 
c     | Back from reverse communication : |
339
 
c     | continue with update step         |
340
 
c     %-----------------------------------%
341
 
c
342
 
      if (update) go to 20
343
 
c
344
 
c     %-------------------------------------------%
345
 
c     | Back from computing user specified shifts |
346
 
c     %-------------------------------------------%
347
 
c
348
 
      if (ushift) go to 50
349
 
c
350
 
c     %-------------------------------------%
351
 
c     | Back from computing residual norm   |
352
 
c     | at the end of the current iteration |
353
 
c     %-------------------------------------%
354
 
c
355
 
      if (cnorm)  go to 100
356
 
c 
357
 
c     %----------------------------------------------------------%
358
 
c     | Compute the first NEV steps of the Arnoldi factorization |
359
 
c     %----------------------------------------------------------%
360
 
c
361
 
      call snaitr (ido, bmat, n, 0, nev, mode, resid, rnorm, v, ldv, 
362
 
     &             h, ldh, ipntr, workd, info)
363
 
c 
364
 
c     %---------------------------------------------------%
365
 
c     | ido .ne. 99 implies use of reverse communication  |
366
 
c     | to compute operations involving OP and possibly B |
367
 
c     %---------------------------------------------------%
368
 
c
369
 
      if (ido .ne. 99) go to 9000
370
 
c
371
 
      if (info .gt. 0) then
372
 
         np   = info
373
 
         mxiter = iter
374
 
         info = -9999
375
 
         go to 1200
376
 
      end if
377
 
c 
378
 
c     %--------------------------------------------------------------%
379
 
c     |                                                              |
380
 
c     |           M A I N  ARNOLDI  I T E R A T I O N  L O O P       |
381
 
c     |           Each iteration implicitly restarts the Arnoldi     |
382
 
c     |           factorization in place.                            |
383
 
c     |                                                              |
384
 
c     %--------------------------------------------------------------%
385
 
c 
386
 
 1000 continue
387
 
c
388
 
         iter = iter + 1
389
 
c
390
 
         if (msglvl .gt. 0) then
391
 
            call ivout (logfil, 1, iter, ndigit, 
392
 
     &           '_naup2: **** Start of major iteration number ****')
393
 
         end if
394
 
c 
395
 
c        %-----------------------------------------------------------%
396
 
c        | Compute NP additional steps of the Arnoldi factorization. |
397
 
c        | Adjust NP since NEV might have been updated by last call  |
398
 
c        | to the shift application routine snapps.                  |
399
 
c        %-----------------------------------------------------------%
400
 
c
401
 
         np  = kplusp - nev
402
 
c
403
 
         if (msglvl .gt. 1) then
404
 
            call ivout (logfil, 1, nev, ndigit, 
405
 
     &     '_naup2: The length of the current Arnoldi factorization')
406
 
            call ivout (logfil, 1, np, ndigit, 
407
 
     &           '_naup2: Extend the Arnoldi factorization by')
408
 
         end if
409
 
c
410
 
c        %-----------------------------------------------------------%
411
 
c        | Compute NP additional steps of the Arnoldi factorization. |
412
 
c        %-----------------------------------------------------------%
413
 
c
414
 
         ido = 0
415
 
   20    continue
416
 
         update = .true.
417
 
c
418
 
         call snaitr (ido  , bmat, n  , nev, np , mode , resid, 
419
 
     &                rnorm, v   , ldv, h  , ldh, ipntr, workd,
420
 
     &                info)
421
 
c 
422
 
c        %---------------------------------------------------%
423
 
c        | ido .ne. 99 implies use of reverse communication  |
424
 
c        | to compute operations involving OP and possibly B |
425
 
c        %---------------------------------------------------%
426
 
c
427
 
         if (ido .ne. 99) go to 9000
428
 
c
429
 
         if (info .gt. 0) then
430
 
            np = info
431
 
            mxiter = iter
432
 
            info = -9999
433
 
            go to 1200
434
 
         end if
435
 
         update = .false.
436
 
c
437
 
         if (msglvl .gt. 1) then
438
 
            call svout (logfil, 1, rnorm, ndigit, 
439
 
     &           '_naup2: Corresponding B-norm of the residual')
440
 
         end if
441
 
c 
442
 
c        %--------------------------------------------------------%
443
 
c        | Compute the eigenvalues and corresponding error bounds |
444
 
c        | of the current upper Hessenberg matrix.                |
445
 
c        %--------------------------------------------------------%
446
 
c
447
 
         call sneigh (rnorm, kplusp, h, ldh, ritzr, ritzi, bounds,
448
 
     &                q, ldq, workl, ierr)
449
 
c
450
 
         if (ierr .ne. 0) then
451
 
            info = -8
452
 
            go to 1200
453
 
         end if
454
 
c
455
 
c        %----------------------------------------------------%
456
 
c        | Make a copy of eigenvalues and corresponding error |
457
 
c        | bounds obtained from sneigh.                       |
458
 
c        %----------------------------------------------------%
459
 
c
460
 
         call scopy(kplusp, ritzr, 1, workl(kplusp**2+1), 1)
461
 
         call scopy(kplusp, ritzi, 1, workl(kplusp**2+kplusp+1), 1)
462
 
         call scopy(kplusp, bounds, 1, workl(kplusp**2+2*kplusp+1), 1)
463
 
c
464
 
c        %---------------------------------------------------%
465
 
c        | Select the wanted Ritz values and their bounds    |
466
 
c        | to be used in the convergence test.               |
467
 
c        | The wanted part of the spectrum and corresponding |
468
 
c        | error bounds are in the last NEV loc. of RITZR,   |
469
 
c        | RITZI and BOUNDS respectively. The variables NEV  |
470
 
c        | and NP may be updated if the NEV-th wanted Ritz   |
471
 
c        | value has a non zero imaginary part. In this case |
472
 
c        | NEV is increased by one and NP decreased by one.  |
473
 
c        | NOTE: The last two arguments of sngets are no     |
474
 
c        | longer used as of version 2.1.                    |
475
 
c        %---------------------------------------------------%
476
 
c
477
 
         nev = nev0
478
 
         np = np0
479
 
         numcnv = nev
480
 
         call sngets (ishift, which, nev, np, ritzr, ritzi, 
481
 
     &                bounds, workl, workl(np+1))
482
 
         if (nev .eq. nev0+1) numcnv = nev0+1
483
 
c 
484
 
c        %-------------------%
485
 
c        | Convergence test. | 
486
 
c        %-------------------%
487
 
c
488
 
         call scopy (nev, bounds(np+1), 1, workl(2*np+1), 1)
489
 
         call snconv (nev, ritzr(np+1), ritzi(np+1), workl(2*np+1), 
490
 
     &        tol, nconv)
491
 
c 
492
 
         if (msglvl .gt. 2) then
493
 
            kp(1) = nev
494
 
            kp(2) = np
495
 
            kp(3) = numcnv
496
 
            kp(4) = nconv
497
 
            call ivout (logfil, 4, kp, ndigit, 
498
 
     &                  '_naup2: NEV, NP, NUMCNV, NCONV are')
499
 
            call svout (logfil, kplusp, ritzr, ndigit,
500
 
     &           '_naup2: Real part of the eigenvalues of H')
501
 
            call svout (logfil, kplusp, ritzi, ndigit,
502
 
     &           '_naup2: Imaginary part of the eigenvalues of H')
503
 
            call svout (logfil, kplusp, bounds, ndigit, 
504
 
     &          '_naup2: Ritz estimates of the current NCV Ritz values')
505
 
         end if
506
 
c
507
 
c        %---------------------------------------------------------%
508
 
c        | Count the number of unwanted Ritz values that have zero |
509
 
c        | Ritz estimates. If any Ritz estimates are equal to zero |
510
 
c        | then a leading block of H of order equal to at least    |
511
 
c        | the number of Ritz values with zero Ritz estimates has  |
512
 
c        | split off. None of these Ritz values may be removed by  |
513
 
c        | shifting. Decrease NP the number of shifts to apply. If |
514
 
c        | no shifts may be applied, then prepare to exit          |
515
 
c        %---------------------------------------------------------%
516
 
c
517
 
         nptemp = np
518
 
         do 30 j=1, nptemp
519
 
            if (bounds(j) .eq. zero) then
520
 
               np = np - 1
521
 
               nev = nev + 1
522
 
            end if
523
 
 30      continue
524
 
c     
525
 
         if ( (nconv .ge. numcnv) .or. 
526
 
     &        (iter .gt. mxiter) .or.
527
 
     &        (np .eq. 0) ) then
528
 
c
529
 
            if (msglvl .gt. 4) then
530
 
               call svout(logfil, kplusp, workl(kplusp**2+1), ndigit,
531
 
     &             '_naup2: Real part of the eig computed by _neigh:')
532
 
               call svout(logfil, kplusp, workl(kplusp**2+kplusp+1),
533
 
     &                     ndigit,
534
 
     &             '_naup2: Imag part of the eig computed by _neigh:')
535
 
               call svout(logfil, kplusp, workl(kplusp**2+kplusp*2+1),
536
 
     &                     ndigit,
537
 
     &             '_naup2: Ritz eistmates computed by _neigh:')
538
 
            end if
539
 
c     
540
 
c           %------------------------------------------------%
541
 
c           | Prepare to exit. Put the converged Ritz values |
542
 
c           | and corresponding bounds in RITZ(1:NCONV) and  |
543
 
c           | BOUNDS(1:NCONV) respectively. Then sort. Be    |
544
 
c           | careful when NCONV > NP                        |
545
 
c           %------------------------------------------------%
546
 
c
547
 
c           %------------------------------------------%
548
 
c           |  Use h( 3,1 ) as storage to communicate  |
549
 
c           |  rnorm to _neupd if needed               |
550
 
c           %------------------------------------------%
551
 
 
552
 
            h(3,1) = rnorm
553
 
c
554
 
c           %----------------------------------------------%
555
 
c           | To be consistent with sngets, we first do a  |
556
 
c           | pre-processing sort in order to keep complex |
557
 
c           | conjugate pairs together.  This is similar   |
558
 
c           | to the pre-processing sort used in sngets    |
559
 
c           | except that the sort is done in the opposite |
560
 
c           | order.                                       |
561
 
c           %----------------------------------------------%
562
 
c
563
 
            if (which .eq. 'LM') wprime = 'SR'
564
 
            if (which .eq. 'SM') wprime = 'LR'
565
 
            if (which .eq. 'LR') wprime = 'SM'
566
 
            if (which .eq. 'SR') wprime = 'LM'
567
 
            if (which .eq. 'LI') wprime = 'SM'
568
 
            if (which .eq. 'SI') wprime = 'LM'
569
 
c
570
 
            call ssortc (wprime, .true., kplusp, ritzr, ritzi, bounds)
571
 
c
572
 
c           %----------------------------------------------%
573
 
c           | Now sort Ritz values so that converged Ritz  |
574
 
c           | values appear within the first NEV locations |
575
 
c           | of ritzr, ritzi and bounds, and the most     |
576
 
c           | desired one appears at the front.            |
577
 
c           %----------------------------------------------%
578
 
c
579
 
            if (which .eq. 'LM') wprime = 'SM'
580
 
            if (which .eq. 'SM') wprime = 'LM'
581
 
            if (which .eq. 'LR') wprime = 'SR'
582
 
            if (which .eq. 'SR') wprime = 'LR'
583
 
            if (which .eq. 'LI') wprime = 'SI'
584
 
            if (which .eq. 'SI') wprime = 'LI'
585
 
c
586
 
            call ssortc(wprime, .true., kplusp, ritzr, ritzi, bounds)
587
 
c
588
 
c           %--------------------------------------------------%
589
 
c           | Scale the Ritz estimate of each Ritz value       |
590
 
c           | by 1 / max(eps23,magnitude of the Ritz value).   |
591
 
c           %--------------------------------------------------%
592
 
c
593
 
            do 35 j = 1, numcnv
594
 
                temp = max(eps23,slapy2(ritzr(j),
595
 
     &                                   ritzi(j)))
596
 
                bounds(j) = bounds(j)/temp
597
 
 35         continue
598
 
c
599
 
c           %----------------------------------------------------%
600
 
c           | Sort the Ritz values according to the scaled Ritz  |
601
 
c           | esitmates.  This will push all the converged ones  |
602
 
c           | towards the front of ritzr, ritzi, bounds          |
603
 
c           | (in the case when NCONV < NEV.)                    |
604
 
c           %----------------------------------------------------%
605
 
c
606
 
            wprime = 'LR'
607
 
            call ssortc(wprime, .true., numcnv, bounds, ritzr, ritzi)
608
 
c
609
 
c           %----------------------------------------------%
610
 
c           | Scale the Ritz estimate back to its original |
611
 
c           | value.                                       |
612
 
c           %----------------------------------------------%
613
 
c
614
 
            do 40 j = 1, numcnv
615
 
                temp = max(eps23, slapy2(ritzr(j),
616
 
     &                                   ritzi(j)))
617
 
                bounds(j) = bounds(j)*temp
618
 
 40         continue
619
 
c
620
 
c           %------------------------------------------------%
621
 
c           | Sort the converged Ritz values again so that   |
622
 
c           | the "threshold" value appears at the front of  |
623
 
c           | ritzr, ritzi and bound.                        |
624
 
c           %------------------------------------------------%
625
 
c
626
 
            call ssortc(which, .true., nconv, ritzr, ritzi, bounds)
627
 
c
628
 
            if (msglvl .gt. 1) then
629
 
               call svout (logfil, kplusp, ritzr, ndigit,
630
 
     &            '_naup2: Sorted real part of the eigenvalues')
631
 
               call svout (logfil, kplusp, ritzi, ndigit,
632
 
     &            '_naup2: Sorted imaginary part of the eigenvalues')
633
 
               call svout (logfil, kplusp, bounds, ndigit,
634
 
     &            '_naup2: Sorted ritz estimates.')
635
 
            end if
636
 
c
637
 
c           %------------------------------------%
638
 
c           | Max iterations have been exceeded. | 
639
 
c           %------------------------------------%
640
 
c
641
 
            if (iter .gt. mxiter .and. nconv .lt. numcnv) info = 1
642
 
c
643
 
c           %---------------------%
644
 
c           | No shifts to apply. | 
645
 
c           %---------------------%
646
 
c
647
 
            if (np .eq. 0 .and. nconv .lt. numcnv) info = 2
648
 
c
649
 
            np = nconv
650
 
            go to 1100
651
 
c
652
 
         else if ( (nconv .lt. numcnv) .and. (ishift .eq. 1) ) then
653
 
c     
654
 
c           %-------------------------------------------------%
655
 
c           | Do not have all the requested eigenvalues yet.  |
656
 
c           | To prevent possible stagnation, adjust the size |
657
 
c           | of NEV.                                         |
658
 
c           %-------------------------------------------------%
659
 
c
660
 
            nevbef = nev
661
 
            nev = nev + min(nconv, np/2)
662
 
            if (nev .eq. 1 .and. kplusp .ge. 6) then
663
 
               nev = kplusp / 2
664
 
            else if (nev .eq. 1 .and. kplusp .gt. 3) then
665
 
               nev = 2
666
 
            end if
667
 
            np = kplusp - nev
668
 
c     
669
 
c           %---------------------------------------%
670
 
c           | If the size of NEV was just increased |
671
 
c           | resort the eigenvalues.               |
672
 
c           %---------------------------------------%
673
 
c     
674
 
            if (nevbef .lt. nev) 
675
 
     &         call sngets (ishift, which, nev, np, ritzr, ritzi, 
676
 
     &              bounds, workl, workl(np+1))
677
 
c
678
 
         end if              
679
 
c     
680
 
         if (msglvl .gt. 0) then
681
 
            call ivout (logfil, 1, nconv, ndigit, 
682
 
     &           '_naup2: no. of "converged" Ritz values at this iter.')
683
 
            if (msglvl .gt. 1) then
684
 
               kp(1) = nev
685
 
               kp(2) = np
686
 
               call ivout (logfil, 2, kp, ndigit, 
687
 
     &              '_naup2: NEV and NP are')
688
 
               call svout (logfil, nev, ritzr(np+1), ndigit,
689
 
     &              '_naup2: "wanted" Ritz values -- real part')
690
 
               call svout (logfil, nev, ritzi(np+1), ndigit,
691
 
     &              '_naup2: "wanted" Ritz values -- imag part')
692
 
               call svout (logfil, nev, bounds(np+1), ndigit,
693
 
     &              '_naup2: Ritz estimates of the "wanted" values ')
694
 
            end if
695
 
         end if
696
 
c
697
 
         if (ishift .eq. 0) then
698
 
c
699
 
c           %-------------------------------------------------------%
700
 
c           | User specified shifts: reverse comminucation to       |
701
 
c           | compute the shifts. They are returned in the first    |
702
 
c           | 2*NP locations of WORKL.                              |
703
 
c           %-------------------------------------------------------%
704
 
c
705
 
            ushift = .true.
706
 
            ido = 3
707
 
            go to 9000
708
 
         end if
709
 
c 
710
 
   50    continue
711
 
c
712
 
c        %------------------------------------%
713
 
c        | Back from reverse communication;   |
714
 
c        | User specified shifts are returned |
715
 
c        | in WORKL(1:2*NP)                   |
716
 
c        %------------------------------------%
717
 
c
718
 
         ushift = .false.
719
 
c
720
 
         if ( ishift .eq. 0 ) then
721
 
c 
722
 
c            %----------------------------------%
723
 
c            | Move the NP shifts from WORKL to |
724
 
c            | RITZR, RITZI to free up WORKL    |
725
 
c            | for non-exact shift case.        |
726
 
c            %----------------------------------%
727
 
c
728
 
             call scopy (np, workl,       1, ritzr, 1)
729
 
             call scopy (np, workl(np+1), 1, ritzi, 1)
730
 
         end if
731
 
c
732
 
         if (msglvl .gt. 2) then 
733
 
            call ivout (logfil, 1, np, ndigit, 
734
 
     &                  '_naup2: The number of shifts to apply ')
735
 
            call svout (logfil, np, ritzr, ndigit,
736
 
     &                  '_naup2: Real part of the shifts')
737
 
            call svout (logfil, np, ritzi, ndigit,
738
 
     &                  '_naup2: Imaginary part of the shifts')
739
 
            if ( ishift .eq. 1 ) 
740
 
     &          call svout (logfil, np, bounds, ndigit,
741
 
     &                  '_naup2: Ritz estimates of the shifts')
742
 
         end if
743
 
c
744
 
c        %---------------------------------------------------------%
745
 
c        | Apply the NP implicit shifts by QR bulge chasing.       |
746
 
c        | Each shift is applied to the whole upper Hessenberg     |
747
 
c        | matrix H.                                               |
748
 
c        | The first 2*N locations of WORKD are used as workspace. |
749
 
c        %---------------------------------------------------------%
750
 
c
751
 
         call snapps (n, nev, np, ritzr, ritzi, v, ldv, 
752
 
     &                h, ldh, resid, q, ldq, workl, workd)
753
 
c
754
 
c        %---------------------------------------------%
755
 
c        | Compute the B-norm of the updated residual. |
756
 
c        | Keep B*RESID in WORKD(1:N) to be used in    |
757
 
c        | the first step of the next call to snaitr.  |
758
 
c        %---------------------------------------------%
759
 
c
760
 
         cnorm = .true.
761
 
         call second (t2)
762
 
         if (bmat .eq. 'G') then
763
 
            nbx = nbx + 1
764
 
            call scopy (n, resid, 1, workd(n+1), 1)
765
 
            ipntr(1) = n + 1
766
 
            ipntr(2) = 1
767
 
            ido = 2
768
 
c 
769
 
c           %----------------------------------%
770
 
c           | Exit in order to compute B*RESID |
771
 
c           %----------------------------------%
772
 
c 
773
 
            go to 9000
774
 
         else if (bmat .eq. 'I') then
775
 
            call scopy (n, resid, 1, workd, 1)
776
 
         end if
777
 
c 
778
 
  100    continue
779
 
c 
780
 
c        %----------------------------------%
781
 
c        | Back from reverse communication; |
782
 
c        | WORKD(1:N) := B*RESID            |
783
 
c        %----------------------------------%
784
 
c
785
 
         if (bmat .eq. 'G') then
786
 
            call second (t3)
787
 
            tmvbx = tmvbx + (t3 - t2)
788
 
         end if
789
 
c 
790
 
         if (bmat .eq. 'G') then         
791
 
            rnorm = sdot (n, resid, 1, workd, 1)
792
 
            rnorm = sqrt(abs(rnorm))
793
 
         else if (bmat .eq. 'I') then
794
 
            rnorm = snrm2(n, resid, 1)
795
 
         end if
796
 
         cnorm = .false.
797
 
c
798
 
         if (msglvl .gt. 2) then
799
 
            call svout (logfil, 1, rnorm, ndigit, 
800
 
     &      '_naup2: B-norm of residual for compressed factorization')
801
 
            call smout (logfil, nev, nev, h, ldh, ndigit,
802
 
     &        '_naup2: Compressed upper Hessenberg matrix H')
803
 
         end if
804
 
c 
805
 
      go to 1000
806
 
c
807
 
c     %---------------------------------------------------------------%
808
 
c     |                                                               |
809
 
c     |  E N D     O F     M A I N     I T E R A T I O N     L O O P  |
810
 
c     |                                                               |
811
 
c     %---------------------------------------------------------------%
812
 
c
813
 
 1100 continue
814
 
c
815
 
      mxiter = iter
816
 
      nev = numcnv
817
 
c     
818
 
 1200 continue
819
 
      ido = 99
820
 
c
821
 
c     %------------%
822
 
c     | Error Exit |
823
 
c     %------------%
824
 
c
825
 
      call second (t1)
826
 
      tnaup2 = t1 - t0
827
 
c     
828
 
 9000 continue
829
 
c
830
 
c     %---------------%
831
 
c     | End of snaup2 |
832
 
c     %---------------%
833
 
c
834
 
      return
835
 
      end