1
c-----------------------------------------------------------------------
7
c Intermediate level interface called by dsaupd.
11
c ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD,
12
c ISHIFT, MXITER, V, LDV, H, LDH, RITZ, BOUNDS, Q, LDQ, WORKL,
13
c IPNTR, WORKD, INFO )
17
c IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in dsaupd.
18
c MODE, ISHIFT, MXITER: see the definition of IPARAM in dsaupd.
20
c NP Integer. (INPUT/OUTPUT)
21
c Contains the number of implicit shifts to apply during
22
c each Arnoldi/Lanczos iteration.
23
c If ISHIFT=1, NP is adjusted dynamically at each iteration
24
c to accelerate convergence and prevent stagnation.
25
c This is also roughly equal to the number of matrix-vector
26
c products (involving the operator OP) per Arnoldi iteration.
27
c The logic for adjusting is contained within the current
29
c If ISHIFT=0, NP is the number of shifts the user needs
30
c to provide via reverse comunication. 0 < NP < NCV-NEV.
31
c NP may be less than NCV-NEV since a leading block of the current
32
c upper Tridiagonal matrix has split off and contains "unwanted"
34
c Upon termination of the IRA iteration, NP contains the number
35
c of "converged" wanted Ritz values.
37
c IUPD Integer. (INPUT)
38
c IUPD .EQ. 0: use explicit restart instead implicit update.
39
c IUPD .NE. 0: use implicit update.
41
c V Double precision N by (NEV+NP) array. (INPUT/OUTPUT)
42
c The Lanczos basis vectors.
44
c LDV Integer. (INPUT)
45
c Leading dimension of V exactly as declared in the calling
48
c H Double precision (NEV+NP) by 2 array. (OUTPUT)
49
c H is used to store the generated symmetric tridiagonal matrix
50
c The subdiagonal is stored in the first column of H starting
51
c at H(2,1). The main diagonal is stored in the second column
52
c of H starting at H(1,2). If dsaup2 converges store the
53
c B-norm of the final residual vector in H(1,1).
55
c LDH Integer. (INPUT)
56
c Leading dimension of H exactly as declared in the calling
59
c RITZ Double precision array of length NEV+NP. (OUTPUT)
60
c RITZ(1:NEV) contains the computed Ritz values of OP.
62
c BOUNDS Double precision array of length NEV+NP. (OUTPUT)
63
c BOUNDS(1:NEV) contain the error bounds corresponding to RITZ.
65
c Q Double precision (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.
69
c LDQ Integer. (INPUT)
70
c Leading dimension of Q exactly as declared in the calling
73
c WORKL Double precision array of length at least 3*(NEV+NP). (INPUT/WORKSPACE)
74
c Private (replicated) array on each PE or array allocated on
75
c the front end. It is used in the computation of the
76
c tridiagonal eigenvalue problem, the calculation and
77
c application of the shifts and convergence checking.
78
c If ISHIFT .EQ. O and IDO .EQ. 3, the first NP locations
79
c of WORKL are used in reverse communication to hold the user
82
c IPNTR Integer array of length 3. (OUTPUT)
83
c Pointer to mark the starting locations in the WORKD for
84
c vectors used by the Lanczos iteration.
85
c -------------------------------------------------------------
86
c IPNTR(1): pointer to the current operand vector X.
87
c IPNTR(2): pointer to the current result vector Y.
88
c IPNTR(3): pointer to the vector B * X when used in one of
89
c the spectral transformation modes. X is the current
91
c -------------------------------------------------------------
93
c WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
94
c Distributed array to be used in the basic Lanczos iteration
95
c for reverse communication. The user should not use WORKD
96
c as temporary workspace during the iteration !!!!!!!!!!
97
c See Data Distribution Note in dsaupd.
99
c INFO Integer. (INPUT/OUTPUT)
100
c If INFO .EQ. 0, a randomly initial residual vector is used.
101
c If INFO .NE. 0, RESID contains the initial residual vector,
102
c possibly from a previous run.
103
c Error flag on output.
104
c = 0: Normal return.
105
c = 1: All possible eigenvalues of OP has been found.
106
c NP returns the size of the invariant subspace
107
c spanning the operator OP.
108
c = 2: No shifts could be applied.
109
c = -8: Error return from trid. eigenvalue calculation;
110
c This should never happen.
111
c = -9: Starting vector is zero.
112
c = -9999: Could not build an Lanczos factorization.
113
c Size that was built in returned in NP.
117
c-----------------------------------------------------------------------
122
c 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
123
c a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
125
c 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
126
c Restarted Arnoldi Iteration", Rice University Technical Report
127
c TR95-13, Department of Computational and Applied Mathematics.
128
c 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
130
c 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
131
c Computer Physics Communications, 53 (1989), pp 169-179.
132
c 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
133
c Implement the Spectral Transformation", Math. Comp., 48 (1987),
135
c 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
136
c Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
137
c SIAM J. Matr. Anal. Apps., January (1993).
138
c 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
139
c for Updating the QR decomposition", ACM TOMS, December 1990,
140
c Volume 16 Number 4, pp 369-377.
143
c dgetv0 ARPACK initial vector generation routine.
144
c dsaitr ARPACK Lanczos factorization routine.
145
c dsapps ARPACK application of implicit shifts routine.
146
c dsconv ARPACK convergence of Ritz values routine.
147
c dseigt ARPACK compute Ritz values and error bounds routine.
148
c dsgets ARPACK reorder Ritz values and error bounds routine.
149
c dsortr ARPACK sorting routine.
150
c ivout ARPACK utility routine that prints integers.
151
c second ARPACK utility routine for timing.
152
c dvout ARPACK utility routine that prints vectors.
153
c dlamch LAPACK routine that determines machine constants.
154
c dcopy Level 1 BLAS that copies one vector to another.
155
c ddot Level 1 BLAS that computes the scalar product of two vectors.
156
c dnrm2 Level 1 BLAS that computes the norm of a vector.
157
c dscal Level 1 BLAS that scales a vector.
158
c dswap Level 1 BLAS that swaps two vectors.
161
c Danny Sorensen Phuong Vu
162
c Richard Lehoucq CRPC / Rice University
163
c Dept. of Computational & Houston, Texas
164
c Applied Mathematics
169
c 12/15/93: Version ' 2.4'
170
c xx/xx/95: Version ' 2.4'. (R.B. Lehoucq)
172
c\SCCS Information: @(#)
173
c FILE: saup2.F SID: 2.7 DATE OF SID: 5/19/98 RELEASE: 2
177
c-----------------------------------------------------------------------
180
& ( ido, bmat, n, which, nev, np, tol, resid, mode, iupd,
181
& ishift, mxiter, v, ldv, h, ldh, ritz, bounds,
182
& q, ldq, workl, ipntr, workd, info )
184
c %----------------------------------------------------%
185
c | Include files for debugging and timing information |
186
c %----------------------------------------------------%
191
c %------------------%
192
c | Scalar Arguments |
193
c %------------------%
195
character bmat*1, which*2
196
integer ido, info, ishift, iupd, ldh, ldq, ldv, mxiter,
201
c %-----------------%
202
c | Array Arguments |
203
c %-----------------%
207
& bounds(nev+np), h(ldh,2), q(ldq,nev+np), resid(n),
208
& ritz(nev+np), v(ldv,nev+np), workd(3*n),
217
parameter (one = 1.0D+0, zero = 0.0D+0)
224
logical cnorm, getv0, initv, update, ushift
225
integer ierr, iter, j, kplusp, msglvl, nconv, nevbef, nev0,
226
& np0, nptemp, nevd2, nevm2, kp(3)
229
save cnorm, getv0, initv, update, ushift,
230
& iter, kplusp, msglvl, nconv, nev0, np0,
233
c %----------------------%
234
c | External Subroutines |
235
c %----------------------%
237
external dcopy, dgetv0, dsaitr, dscal, dsconv, dseigt, dsgets,
238
& dsapps, dsortr, dvout, ivout, second, dswap
240
c %--------------------%
241
c | External Functions |
242
c %--------------------%
245
& ddot, dnrm2, dlamch
246
external ddot, dnrm2, dlamch
248
c %---------------------%
249
c | Intrinsic Functions |
250
c %---------------------%
254
c %-----------------------%
255
c | Executable Statements |
256
c %-----------------------%
260
c %-------------------------------%
261
c | Initialize timing statistics |
262
c | & message level for debugging |
263
c %-------------------------------%
268
c %---------------------------------%
269
c | Set machine dependent constant. |
270
c %---------------------------------%
272
eps23 = dlamch('Epsilon-Machine')
273
eps23 = eps23**(2.0D+0/3.0D+0)
275
c %-------------------------------------%
276
c | nev0 and np0 are integer variables |
277
c | hold the initial values of NEV & NP |
278
c %-------------------------------------%
283
c %-------------------------------------%
284
c | kplusp is the bound on the largest |
285
c | Lanczos factorization built. |
286
c | nconv is the current number of |
287
c | "converged" eigenvlues. |
288
c | iter is the counter on the current |
289
c | iteration step. |
290
c %-------------------------------------%
296
c %--------------------------------------------%
297
c | Set flags for computing the first NEV steps |
298
c | of the Lanczos factorization. |
299
c %--------------------------------------------%
306
if (info .ne. 0) then
308
c %--------------------------------------------%
309
c | User provides the initial residual vector. |
310
c %--------------------------------------------%
319
c %---------------------------------------------%
320
c | Get a possibly random starting vector and |
321
c | force it into the range of the operator OP. |
322
c %---------------------------------------------%
327
call dgetv0 (ido, bmat, 1, initv, n, 1, v, ldv, resid, rnorm,
328
& ipntr, workd, info)
330
if (ido .ne. 99) go to 9000
332
if (rnorm .eq. zero) then
334
c %-----------------------------------------%
335
c | The initial vector is zero. Error exit. |
336
c %-----------------------------------------%
345
c %------------------------------------------------------------%
346
c | Back from reverse communication: continue with update step |
347
c %------------------------------------------------------------%
351
c %-------------------------------------------%
352
c | Back from computing user specified shifts |
353
c %-------------------------------------------%
357
c %-------------------------------------%
358
c | Back from computing residual norm |
359
c | at the end of the current iteration |
360
c %-------------------------------------%
364
c %----------------------------------------------------------%
365
c | Compute the first NEV steps of the Lanczos factorization |
366
c %----------------------------------------------------------%
368
call dsaitr (ido, bmat, n, 0, nev0, mode, resid, rnorm, v, ldv,
369
& h, ldh, ipntr, workd, info)
371
c %---------------------------------------------------%
372
c | ido .ne. 99 implies use of reverse communication |
373
c | to compute operations involving OP and possibly B |
374
c %---------------------------------------------------%
376
if (ido .ne. 99) go to 9000
378
if (info .gt. 0) then
380
c %-----------------------------------------------------%
381
c | dsaitr was unable to build an Lanczos factorization |
382
c | of length NEV0. INFO is returned with the size of |
383
c | the factorization built. Exit main loop. |
384
c %-----------------------------------------------------%
392
c %--------------------------------------------------------------%
394
c | M A I N LANCZOS I T E R A T I O N L O O P |
395
c | Each iteration implicitly restarts the Lanczos |
396
c | factorization in place. |
398
c %--------------------------------------------------------------%
404
if (msglvl .gt. 0) then
405
call ivout (logfil, 1, iter, ndigit,
406
& '_saup2: **** Start of major iteration number ****')
408
if (msglvl .gt. 1) then
409
call ivout (logfil, 1, nev, ndigit,
410
& '_saup2: The length of the current Lanczos factorization')
411
call ivout (logfil, 1, np, ndigit,
412
& '_saup2: Extend the Lanczos factorization by')
415
c %------------------------------------------------------------%
416
c | Compute NP additional steps of the Lanczos factorization. |
417
c %------------------------------------------------------------%
423
call dsaitr (ido, bmat, n, nev, np, mode, resid, rnorm, v,
424
& ldv, h, ldh, ipntr, workd, info)
426
c %---------------------------------------------------%
427
c | ido .ne. 99 implies use of reverse communication |
428
c | to compute operations involving OP and possibly B |
429
c %---------------------------------------------------%
431
if (ido .ne. 99) go to 9000
433
if (info .gt. 0) then
435
c %-----------------------------------------------------%
436
c | dsaitr was unable to build an Lanczos factorization |
437
c | of length NEV0+NP0. INFO is returned with the size |
438
c | of the factorization built. Exit main loop. |
439
c %-----------------------------------------------------%
448
if (msglvl .gt. 1) then
449
call dvout (logfil, 1, rnorm, ndigit,
450
& '_saup2: Current B-norm of residual for factorization')
453
c %--------------------------------------------------------%
454
c | Compute the eigenvalues and corresponding error bounds |
455
c | of the current symmetric tridiagonal matrix. |
456
c %--------------------------------------------------------%
458
call dseigt (rnorm, kplusp, h, ldh, ritz, bounds, workl, ierr)
460
if (ierr .ne. 0) then
465
c %----------------------------------------------------%
466
c | Make a copy of eigenvalues and corresponding error |
467
c | bounds obtained from _seigt. |
468
c %----------------------------------------------------%
470
call dcopy(kplusp, ritz, 1, workl(kplusp+1), 1)
471
call dcopy(kplusp, bounds, 1, workl(2*kplusp+1), 1)
473
c %---------------------------------------------------%
474
c | Select the wanted Ritz values and their bounds |
475
c | to be used in the convergence test. |
476
c | The selection is based on the requested number of |
477
c | eigenvalues instead of the current NEV and NP to |
478
c | prevent possible misconvergence. |
479
c | * Wanted Ritz values := RITZ(NP+1:NEV+NP) |
480
c | * Shifts := RITZ(1:NP) := WORKL(1:NP) |
481
c %---------------------------------------------------%
485
call dsgets (ishift, which, nev, np, ritz, bounds, workl)
487
c %-------------------%
488
c | Convergence test. |
489
c %-------------------%
491
call dcopy (nev, bounds(np+1), 1, workl(np+1), 1)
492
call dsconv (nev, ritz(np+1), workl(np+1), tol, nconv)
494
if (msglvl .gt. 2) then
498
call ivout (logfil, 3, kp, ndigit,
499
& '_saup2: NEV, NP, NCONV are')
500
call dvout (logfil, kplusp, ritz, ndigit,
501
& '_saup2: The eigenvalues of H')
502
call dvout (logfil, kplusp, bounds, ndigit,
503
& '_saup2: Ritz estimates of the current NCV Ritz values')
506
c %---------------------------------------------------------%
507
c | Count the number of unwanted Ritz values that have zero |
508
c | Ritz estimates. If any Ritz estimates are equal to zero |
509
c | then a leading block of H of order equal to at least |
510
c | the number of Ritz values with zero Ritz estimates has |
511
c | split off. None of these Ritz values may be removed by |
512
c | shifting. Decrease NP the number of shifts to apply. If |
513
c | no shifts may be applied, then prepare to exit |
514
c %---------------------------------------------------------%
518
if (bounds(j) .eq. zero) then
524
if ( (nconv .ge. nev0) .or.
525
& (iter .gt. mxiter) .or.
528
c %------------------------------------------------%
529
c | Prepare to exit. Put the converged Ritz values |
530
c | and corresponding bounds in RITZ(1:NCONV) and |
531
c | BOUNDS(1:NCONV) respectively. Then sort. Be |
532
c | careful when NCONV > NP since we don't want to |
533
c | swap overlapping locations. |
534
c %------------------------------------------------%
536
if (which .eq. 'BE') then
538
c %-----------------------------------------------------%
539
c | Both ends of the spectrum are requested. |
540
c | Sort the eigenvalues into algebraically decreasing |
541
c | order first then swap low end of the spectrum next |
542
c | to high end in appropriate locations. |
543
c | NOTE: when np < floor(nev/2) be careful not to swap |
544
c | overlapping locations. |
545
c %-----------------------------------------------------%
548
call dsortr (wprime, .true., kplusp, ritz, bounds)
551
if ( nev .gt. 1 ) then
552
call dswap ( min(nevd2,np), ritz(nevm2+1), 1,
553
& ritz( max(kplusp-nevd2+1,kplusp-np+1) ), 1)
554
call dswap ( min(nevd2,np), bounds(nevm2+1), 1,
555
& bounds( max(kplusp-nevd2+1,kplusp-np+1)), 1)
560
c %--------------------------------------------------%
561
c | LM, SM, LA, SA case. |
562
c | Sort the eigenvalues of H into the an order that |
563
c | is opposite to WHICH, and apply the resulting |
564
c | order to BOUNDS. The eigenvalues are sorted so |
565
c | that the wanted part are always within the first |
567
c %--------------------------------------------------%
569
if (which .eq. 'LM') wprime = 'SM'
570
if (which .eq. 'SM') wprime = 'LM'
571
if (which .eq. 'LA') wprime = 'SA'
572
if (which .eq. 'SA') wprime = 'LA'
574
call dsortr (wprime, .true., kplusp, ritz, bounds)
578
c %--------------------------------------------------%
579
c | Scale the Ritz estimate of each Ritz value |
580
c | by 1 / max(eps23,magnitude of the Ritz value). |
581
c %--------------------------------------------------%
584
temp = max( eps23, abs(ritz(j)) )
585
bounds(j) = bounds(j)/temp
588
c %----------------------------------------------------%
589
c | Sort the Ritz values according to the scaled Ritz |
590
c | esitmates. This will push all the converged ones |
591
c | towards the front of ritzr, ritzi, bounds |
592
c | (in the case when NCONV < NEV.) |
593
c %----------------------------------------------------%
596
call dsortr(wprime, .true., nev0, bounds, ritz)
598
c %----------------------------------------------%
599
c | Scale the Ritz estimate back to its original |
601
c %----------------------------------------------%
604
temp = max( eps23, abs(ritz(j)) )
605
bounds(j) = bounds(j)*temp
608
c %--------------------------------------------------%
609
c | Sort the "converged" Ritz values again so that |
610
c | the "threshold" values and their associated Ritz |
611
c | estimates appear at the appropriate position in |
612
c | ritz and bound. |
613
c %--------------------------------------------------%
615
if (which .eq. 'BE') then
617
c %------------------------------------------------%
618
c | Sort the "converged" Ritz values in increasing |
619
c | order. The "threshold" values are in the |
621
c %------------------------------------------------%
624
call dsortr(wprime, .true., nconv, ritz, bounds)
628
c %----------------------------------------------%
629
c | In LM, SM, LA, SA case, sort the "converged" |
630
c | Ritz values according to WHICH so that the |
631
c | "threshold" value appears at the front of |
633
c %----------------------------------------------%
635
call dsortr(which, .true., nconv, ritz, bounds)
639
c %------------------------------------------%
640
c | Use h( 1,1 ) as storage to communicate |
641
c | rnorm to _seupd if needed |
642
c %------------------------------------------%
646
if (msglvl .gt. 1) then
647
call dvout (logfil, kplusp, ritz, ndigit,
648
& '_saup2: Sorted Ritz values.')
649
call dvout (logfil, kplusp, bounds, ndigit,
650
& '_saup2: Sorted ritz estimates.')
653
c %------------------------------------%
654
c | Max iterations have been exceeded. |
655
c %------------------------------------%
657
if (iter .gt. mxiter .and. nconv .lt. nev) info = 1
659
c %---------------------%
660
c | No shifts to apply. |
661
c %---------------------%
663
if (np .eq. 0 .and. nconv .lt. nev0) info = 2
668
else if (nconv .lt. nev .and. ishift .eq. 1) then
670
c %---------------------------------------------------%
671
c | Do not have all the requested eigenvalues yet. |
672
c | To prevent possible stagnation, adjust the number |
673
c | of Ritz values and the shifts. |
674
c %---------------------------------------------------%
677
nev = nev + min (nconv, np/2)
678
if (nev .eq. 1 .and. kplusp .ge. 6) then
680
else if (nev .eq. 1 .and. kplusp .gt. 2) then
685
c %---------------------------------------%
686
c | If the size of NEV was just increased |
687
c | resort the eigenvalues. |
688
c %---------------------------------------%
691
& call dsgets (ishift, which, nev, np, ritz, bounds,
696
if (msglvl .gt. 0) then
697
call ivout (logfil, 1, nconv, ndigit,
698
& '_saup2: no. of "converged" Ritz values at this iter.')
699
if (msglvl .gt. 1) then
702
call ivout (logfil, 2, kp, ndigit,
703
& '_saup2: NEV and NP are')
704
call dvout (logfil, nev, ritz(np+1), ndigit,
705
& '_saup2: "wanted" Ritz values.')
706
call dvout (logfil, nev, bounds(np+1), ndigit,
707
& '_saup2: Ritz estimates of the "wanted" values ')
712
if (ishift .eq. 0) then
714
c %-----------------------------------------------------%
715
c | User specified shifts: reverse communication to |
716
c | compute the shifts. They are returned in the first |
717
c | NP locations of WORKL. |
718
c %-----------------------------------------------------%
727
c %------------------------------------%
728
c | Back from reverse communication; |
729
c | User specified shifts are returned |
730
c | in WORKL(1:*NP) |
731
c %------------------------------------%
736
c %---------------------------------------------------------%
737
c | Move the NP shifts to the first NP locations of RITZ to |
738
c | free up WORKL. This is for the non-exact shift case; |
739
c | in the exact shift case, dsgets already handles this. |
740
c %---------------------------------------------------------%
742
if (ishift .eq. 0) call dcopy (np, workl, 1, ritz, 1)
744
if (msglvl .gt. 2) then
745
call ivout (logfil, 1, np, ndigit,
746
& '_saup2: The number of shifts to apply ')
747
call dvout (logfil, np, workl, ndigit,
748
& '_saup2: shifts selected')
749
if (ishift .eq. 1) then
750
call dvout (logfil, np, bounds, ndigit,
751
& '_saup2: corresponding Ritz estimates')
755
c %---------------------------------------------------------%
756
c | Apply the NP0 implicit shifts by QR bulge chasing. |
757
c | Each shift is applied to the entire tridiagonal matrix. |
758
c | The first 2*N locations of WORKD are used as workspace. |
759
c | After dsapps is done, we have a Lanczos |
760
c | factorization of length NEV. |
761
c %---------------------------------------------------------%
763
call dsapps (n, nev, np, ritz, v, ldv, h, ldh, resid, q, ldq,
766
c %---------------------------------------------%
767
c | Compute the B-norm of the updated residual. |
768
c | Keep B*RESID in WORKD(1:N) to be used in |
769
c | the first step of the next call to dsaitr. |
770
c %---------------------------------------------%
774
if (bmat .eq. 'G') then
776
call dcopy (n, resid, 1, workd(n+1), 1)
781
c %----------------------------------%
782
c | Exit in order to compute B*RESID |
783
c %----------------------------------%
786
else if (bmat .eq. 'I') then
787
call dcopy (n, resid, 1, workd, 1)
792
c %----------------------------------%
793
c | Back from reverse communication; |
794
c | WORKD(1:N) := B*RESID |
795
c %----------------------------------%
797
if (bmat .eq. 'G') then
799
tmvbx = tmvbx + (t3 - t2)
802
if (bmat .eq. 'G') then
803
rnorm = ddot (n, resid, 1, workd, 1)
804
rnorm = sqrt(abs(rnorm))
805
else if (bmat .eq. 'I') then
806
rnorm = dnrm2(n, resid, 1)
811
if (msglvl .gt. 2) then
812
call dvout (logfil, 1, rnorm, ndigit,
813
& '_saup2: B-norm of residual for NEV factorization')
814
call dvout (logfil, nev, h(1,2), ndigit,
815
& '_saup2: main diagonal of compressed H matrix')
816
call dvout (logfil, nev-1, h(2,1), ndigit,
817
& '_saup2: subdiagonal of compressed H matrix')
822
c %---------------------------------------------------------------%
824
c | E N D O F M A I N I T E R A T I O N L O O P |
826
c %---------------------------------------------------------------%
added
removed
scipy/sandbox/arpack/ARPACK/SRC/dsaup2.f
