← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/arpack/ARPACK/SRC/snaup2.f

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

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