1
/* -- translated by f2c (version 20100827).
2
You must link the resulting object file with libf2c:
3
on Microsoft Windows system, link with libf2c.lib;
4
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5
or, if you install libf2c.a in a standard place, with -lf2c -lm
6
-- in that order, at the end of the command line, as in
8
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
10
http://www.netlib.org/f2c/libf2c.zip
15
/* Table of constant values */
17
static doublereal c_b3 = .66666666666666663;
18
static integer c__1 = 1;
19
static integer c__0 = 0;
20
static integer c__3 = 3;
21
static logical c_true = TRUE_;
22
static integer c__2 = 2;
24
/* -----------------------------------------------------------------------
30
Intermediate level interface called by dsaupd.
34
( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD,
35
ISHIFT, MXITER, V, LDV, H, LDH, RITZ, BOUNDS, Q, LDQ, WORKL,
40
IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in dsaupd.
41
MODE, ISHIFT, MXITER: see the definition of IPARAM in dsaupd.
43
NP Integer. (INPUT/OUTPUT)
44
Contains the number of implicit shifts to apply during
45
each Arnoldi/Lanczos iteration.
46
If ISHIFT=1, NP is adjusted dynamically at each iteration
47
to accelerate convergence and prevent stagnation.
48
This is also roughly equal to the number of matrix-vector
49
products (involving the operator OP) per Arnoldi iteration.
50
The logic for adjusting is contained within the current
52
If ISHIFT=0, NP is the number of shifts the user needs
53
to provide via reverse comunication. 0 < NP < NCV-NEV.
54
NP may be less than NCV-NEV since a leading block of the current
55
upper Tridiagonal matrix has split off and contains "unwanted"
57
Upon termination of the IRA iteration, NP contains the number
58
of "converged" wanted Ritz values.
61
IUPD .EQ. 0: use explicit restart instead implicit update.
62
IUPD .NE. 0: use implicit update.
64
V Double precision N by (NEV+NP) array. (INPUT/OUTPUT)
65
The Lanczos basis vectors.
68
Leading dimension of V exactly as declared in the calling
71
H Double precision (NEV+NP) by 2 array. (OUTPUT)
72
H is used to store the generated symmetric tridiagonal matrix
73
The subdiagonal is stored in the first column of H starting
74
at H(2,1). The main diagonal is stored in the second column
75
of H starting at H(1,2). If dsaup2 converges store the
76
B-norm of the final residual vector in H(1,1).
79
Leading dimension of H exactly as declared in the calling
82
RITZ Double precision array of length NEV+NP. (OUTPUT)
83
RITZ(1:NEV) contains the computed Ritz values of OP.
85
BOUNDS Double precision array of length NEV+NP. (OUTPUT)
86
BOUNDS(1:NEV) contain the error bounds corresponding to RITZ.
88
Q Double precision (NEV+NP) by (NEV+NP) array. (WORKSPACE)
89
Private (replicated) work array used to accumulate the
90
rotation in the shift application step.
93
Leading dimension of Q exactly as declared in the calling
96
WORKL Double precision array of length at least 3*(NEV+NP). (INPUT/WORKSPACE)
97
Private (replicated) array on each PE or array allocated on
98
the front end. It is used in the computation of the
99
tridiagonal eigenvalue problem, the calculation and
100
application of the shifts and convergence checking.
101
If ISHIFT .EQ. O and IDO .EQ. 3, the first NP locations
102
of WORKL are used in reverse communication to hold the user
105
IPNTR Integer array of length 3. (OUTPUT)
106
Pointer to mark the starting locations in the WORKD for
107
vectors used by the Lanczos iteration.
108
-------------------------------------------------------------
109
IPNTR(1): pointer to the current operand vector X.
110
IPNTR(2): pointer to the current result vector Y.
111
IPNTR(3): pointer to the vector B * X when used in one of
112
the spectral transformation modes. X is the current
114
-------------------------------------------------------------
116
WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION)
117
Distributed array to be used in the basic Lanczos iteration
118
for reverse communication. The user should not use WORKD
119
as temporary workspace during the iteration !!!!!!!!!!
120
See Data Distribution Note in dsaupd.
122
INFO Integer. (INPUT/OUTPUT)
123
If INFO .EQ. 0, a randomly initial residual vector is used.
124
If INFO .NE. 0, RESID contains the initial residual vector,
125
possibly from a previous run.
126
Error flag on output.
128
= 1: All possible eigenvalues of OP has been found.
129
NP returns the size of the invariant subspace
130
spanning the operator OP.
131
= 2: No shifts could be applied.
132
= -8: Error return from trid. eigenvalue calculation;
133
This should never happen.
134
= -9: Starting vector is zero.
135
= -9999: Could not build an Lanczos factorization.
136
Size that was built in returned in NP.
140
-----------------------------------------------------------------------
145
1. D.C. Sorensen, "Implicit Application of Polynomial Filters in
146
a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992),
148
2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly
149
Restarted Arnoldi Iteration", Rice University Technical Report
150
TR95-13, Department of Computational and Applied Mathematics.
151
3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall,
153
4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program",
154
Computer Physics Communications, 53 (1989), pp 169-179.
155
5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to
156
Implement the Spectral Transformation", Math. Comp., 48 (1987),
158
6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos
159
Algorithm for Solving Sparse Symmetric Generalized Eigenproblems",
160
SIAM J. Matr. Anal. Apps., January (1993).
161
7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines
162
for Updating the QR decomposition", ACM TOMS, December 1990,
163
Volume 16 Number 4, pp 369-377.
166
dgetv0 ARPACK initial vector generation routine.
167
dsaitr ARPACK Lanczos factorization routine.
168
dsapps ARPACK application of implicit shifts routine.
169
dsconv ARPACK convergence of Ritz values routine.
170
dseigt ARPACK compute Ritz values and error bounds routine.
171
dsgets ARPACK reorder Ritz values and error bounds routine.
172
dsortr ARPACK sorting routine.
173
ivout ARPACK utility routine that prints integers.
174
second ARPACK utility routine for timing.
175
dvout ARPACK utility routine that prints vectors.
176
dlamch LAPACK routine that determines machine constants.
177
dcopy Level 1 BLAS that copies one vector to another.
178
ddot Level 1 BLAS that computes the scalar product of two vectors.
179
dnrm2 Level 1 BLAS that computes the norm of a vector.
180
dscal Level 1 BLAS that scales a vector.
181
dswap Level 1 BLAS that swaps two vectors.
184
Danny Sorensen Phuong Vu
185
Richard Lehoucq CRPC / Rice University
186
Dept. of Computational & Houston, Texas
192
12/15/93: Version ' 2.4'
193
xx/xx/95: Version ' 2.4'. (R.B. Lehoucq)
195
\SCCS Information: @(#)
196
FILE: saup2.F SID: 2.6 DATE OF SID: 8/16/96 RELEASE: 2
200
-----------------------------------------------------------------------
202
Subroutine */ int igraphdsaup2_(integer *ido, char *bmat, integer *n, char *
203
which, integer *nev, integer *np, doublereal *tol, doublereal *resid,
204
integer *mode, integer *iupd, integer *ishift, integer *mxiter,
205
doublereal *v, integer *ldv, doublereal *h__, integer *ldh,
206
doublereal *ritz, doublereal *bounds, doublereal *q, integer *ldq,
207
doublereal *workl, integer *ipntr, doublereal *workd, integer *info)
209
/* System generated locals */
210
integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
212
doublereal d__1, d__2, d__3;
214
/* Builtin functions */
215
double pow_dd(doublereal *, doublereal *);
216
integer s_cmp(char *, char *, ftnlen, ftnlen);
217
/* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
218
double sqrt(doublereal);
220
/* Local variables */
224
IGRAPH_F77_SAVE integer np0;
226
IGRAPH_F77_SAVE integer nev0;
227
extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *,
229
IGRAPH_F77_SAVE doublereal eps23;
231
IGRAPH_F77_SAVE integer iter;
234
extern doublereal igraphdnrm2_(integer *, doublereal *, integer *);
235
IGRAPH_F77_SAVE logical getv0;
237
IGRAPH_F77_SAVE logical cnorm;
238
extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *,
239
doublereal *, integer *), igraphdswap_(integer *, doublereal *, integer
240
*, doublereal *, integer *);
241
IGRAPH_F77_SAVE integer nconv;
242
IGRAPH_F77_SAVE logical initv;
243
IGRAPH_F77_SAVE doublereal rnorm;
245
extern /* Subroutine */ int igraphdvout_(integer *, integer *, doublereal *,
246
integer *, char *, ftnlen), igraphivout_(integer *, integer *, integer *
247
, integer *, char *, ftnlen), igraphdgetv0_(integer *, char *, integer *
248
, logical *, integer *, integer *, doublereal *, integer *,
249
doublereal *, doublereal *, integer *, doublereal *, integer *);
252
extern doublereal igraphdlamch_(char *);
254
extern /* Subroutine */ int igraphsecond_(real *);
255
integer logfil=0, ndigit;
256
extern /* Subroutine */ int igraphdseigt_(doublereal *, integer *, doublereal *,
257
integer *, doublereal *, doublereal *, doublereal *, integer *);
258
IGRAPH_F77_SAVE logical update;
259
extern /* Subroutine */ int igraphdsaitr_(integer *, char *, integer *, integer
260
*, integer *, integer *, doublereal *, doublereal *, doublereal *,
261
integer *, doublereal *, integer *, integer *, doublereal *,
262
integer *), igraphdsgets_(integer *, char *, integer *, integer
263
*, doublereal *, doublereal *, doublereal *), igraphdsapps_(
264
integer *, integer *, integer *, doublereal *, doublereal *,
265
integer *, doublereal *, integer *, doublereal *, doublereal *,
266
integer *, doublereal *), igraphdsconv_(integer *, doublereal *,
267
doublereal *, doublereal *, integer *);
268
IGRAPH_F77_SAVE logical ushift;
270
IGRAPH_F77_SAVE integer msglvl;
272
extern /* Subroutine */ int igraphdsortr_(char *, logical *, integer *,
273
doublereal *, doublereal *);
274
IGRAPH_F77_SAVE integer kplusp;
277
/* %----------------------------------------------------%
278
| Include files for debugging and timing information |
279
%----------------------------------------------------%
302
%----------------------%
303
| External Subroutines |
304
%----------------------%
307
%--------------------%
308
| External Functions |
309
%--------------------%
312
%---------------------%
313
| Intrinsic Functions |
314
%---------------------%
317
%-----------------------%
318
| Executable Statements |
319
%-----------------------%
321
Parameter adjustments */
328
v_offset = 1 + v_dim1;
331
h_offset = 1 + h_dim1;
334
q_offset = 1 + q_dim1;
341
/* %-------------------------------%
342
| Initialize timing statistics |
343
| & message level for debugging |
344
%-------------------------------% */
349
/* %---------------------------------%
350
| Set machine dependent constant. |
351
%---------------------------------% */
353
eps23 = igraphdlamch_("Epsilon-Machine");
354
eps23 = pow_dd(&eps23, &c_b3);
356
/* %-------------------------------------%
357
| nev0 and np0 are integer variables |
358
| hold the initial values of NEV & NP |
359
%-------------------------------------% */
364
/* %-------------------------------------%
365
| kplusp is the bound on the largest |
366
| Lanczos factorization built. |
367
| nconv is the current number of |
368
| "converged" eigenvlues. |
369
| iter is the counter on the current |
371
%-------------------------------------% */
377
/* %--------------------------------------------%
378
| Set flags for computing the first NEV steps |
379
| of the Lanczos factorization. |
380
%--------------------------------------------% */
389
/* %--------------------------------------------%
390
| User provides the initial residual vector. |
391
%--------------------------------------------% */
400
/* %---------------------------------------------%
401
| Get a possibly random starting vector and |
402
| force it into the range of the operator OP. |
403
%---------------------------------------------%
408
igraphdgetv0_(ido, bmat, &c__1, &initv, n, &c__1, &v[v_offset], ldv, &resid[
409
1], &rnorm, &ipntr[1], &workd[1], info);
417
/* %-----------------------------------------%
418
| The initial vector is zero. Error exit. |
419
%-----------------------------------------% */
428
/* %------------------------------------------------------------%
429
| Back from reverse communication: continue with update step |
430
%------------------------------------------------------------% */
436
/* %-------------------------------------------%
437
| Back from computing user specified shifts |
438
%-------------------------------------------% */
444
/* %-------------------------------------%
445
| Back from computing residual norm |
446
| at the end of the current iteration |
447
%-------------------------------------% */
453
/* %----------------------------------------------------------%
454
| Compute the first NEV steps of the Lanczos factorization |
455
%----------------------------------------------------------% */
457
igraphdsaitr_(ido, bmat, n, &c__0, &nev0, mode, &resid[1], &rnorm, &v[v_offset],
458
ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], info);
460
/* %---------------------------------------------------%
461
| ido .ne. 99 implies use of reverse communication |
462
| to compute operations involving OP and possibly B |
463
%---------------------------------------------------% */
471
/* %-----------------------------------------------------%
472
| dsaitr was unable to build an Lanczos factorization |
473
| of length NEV0. INFO is returned with the size of |
474
| the factorization built. Exit main loop. |
475
%-----------------------------------------------------% */
483
/* %--------------------------------------------------------------%
485
| M A I N LANCZOS I T E R A T I O N L O O P |
486
| Each iteration implicitly restarts the Lanczos |
487
| factorization in place. |
489
%--------------------------------------------------------------% */
496
igraphivout_(&logfil, &c__1, &iter, &ndigit, "_saup2: **** Start of major "
497
"iteration number ****", (ftnlen)49);
500
igraphivout_(&logfil, &c__1, nev, &ndigit, "_saup2: The length of the curr"
501
"ent Lanczos factorization", (ftnlen)55);
502
igraphivout_(&logfil, &c__1, np, &ndigit, "_saup2: Extend the Lanczos fact"
503
"orization by", (ftnlen)43);
506
/* %------------------------------------------------------------%
507
| Compute NP additional steps of the Lanczos factorization. |
508
%------------------------------------------------------------% */
514
igraphdsaitr_(ido, bmat, n, nev, np, mode, &resid[1], &rnorm, &v[v_offset], ldv,
515
&h__[h_offset], ldh, &ipntr[1], &workd[1], info);
517
/* %---------------------------------------------------%
518
| ido .ne. 99 implies use of reverse communication |
519
| to compute operations involving OP and possibly B |
520
%---------------------------------------------------% */
528
/* %-----------------------------------------------------%
529
| dsaitr was unable to build an Lanczos factorization |
530
| of length NEV0+NP0. INFO is returned with the size |
531
| of the factorization built. Exit main loop. |
532
%-----------------------------------------------------% */
542
igraphdvout_(&logfil, &c__1, &rnorm, &ndigit, "_saup2: Current B-norm of r"
543
"esidual for factorization", (ftnlen)52);
546
/* %--------------------------------------------------------%
547
| Compute the eigenvalues and corresponding error bounds |
548
| of the current symmetric tridiagonal matrix. |
549
%--------------------------------------------------------% */
551
igraphdseigt_(&rnorm, &kplusp, &h__[h_offset], ldh, &ritz[1], &bounds[1], &
559
/* %----------------------------------------------------%
560
| Make a copy of eigenvalues and corresponding error |
561
| bounds obtained from _seigt. |
562
%----------------------------------------------------% */
564
igraphdcopy_(&kplusp, &ritz[1], &c__1, &workl[kplusp + 1], &c__1);
565
igraphdcopy_(&kplusp, &bounds[1], &c__1, &workl[(kplusp << 1) + 1], &c__1);
567
/* %---------------------------------------------------%
568
| Select the wanted Ritz values and their bounds |
569
| to be used in the convergence test. |
570
| The selection is based on the requested number of |
571
| eigenvalues instead of the current NEV and NP to |
572
| prevent possible misconvergence. |
573
| * Wanted Ritz values := RITZ(NP+1:NEV+NP) |
574
| * Shifts := RITZ(1:NP) := WORKL(1:NP) |
575
%---------------------------------------------------% */
579
igraphdsgets_(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
581
/* %-------------------%
582
| Convergence test. |
583
%-------------------% */
585
igraphdcopy_(nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
586
igraphdsconv_(nev, &ritz[*np + 1], &workl[*np + 1], tol, &nconv);
592
igraphivout_(&logfil, &c__3, kp, &ndigit, "_saup2: NEV, NP, NCONV are", (
594
igraphdvout_(&logfil, &kplusp, &ritz[1], &ndigit, "_saup2: The eigenvalues"
595
" of H", (ftnlen)28);
596
igraphdvout_(&logfil, &kplusp, &bounds[1], &ndigit, "_saup2: Ritz estimate"
597
"s of the current NCV Ritz values", (ftnlen)53);
600
/* %---------------------------------------------------------%
601
| Count the number of unwanted Ritz values that have zero |
602
| Ritz estimates. If any Ritz estimates are equal to zero |
603
| then a leading block of H of order equal to at least |
604
| the number of Ritz values with zero Ritz estimates has |
605
| split off. None of these Ritz values may be removed by |
606
| shifting. Decrease NP the number of shifts to apply. If |
607
| no shifts may be applied, then prepare to exit |
608
%---------------------------------------------------------% */
612
for (j = 1; j <= i__1; ++j) {
613
if (bounds[j] == 0.) {
620
if (nconv >= nev0 || iter > *mxiter || *np == 0) {
622
/* %------------------------------------------------%
623
| Prepare to exit. Put the converged Ritz values |
624
| and corresponding bounds in RITZ(1:NCONV) and |
625
| BOUNDS(1:NCONV) respectively. Then sort. Be |
626
| careful when NCONV > NP since we don't want to |
627
| swap overlapping locations. |
628
%------------------------------------------------% */
630
if (s_cmp(which, "BE", (ftnlen)2, (ftnlen)2) == 0) {
632
/* %-----------------------------------------------------%
633
| Both ends of the spectrum are requested. |
634
| Sort the eigenvalues into algebraically decreasing |
635
| order first then swap low end of the spectrum next |
636
| to high end in appropriate locations. |
637
| NOTE: when np < floor(nev/2) be careful not to swap |
638
| overlapping locations. |
639
%-----------------------------------------------------% */
641
s_copy(wprime, "SA", (ftnlen)2, (ftnlen)2);
642
igraphdsortr_(wprime, &c_true, &kplusp, &ritz[1], &bounds[1])
645
nevm2 = *nev - nevd2;
647
i__1 = min(nevd2,*np);
649
i__2 = kplusp - nevd2 + 1, i__3 = kplusp - *np + 1;
650
igraphdswap_(&i__1, &ritz[nevm2 + 1], &c__1, &ritz[max(i__2,i__3)],
652
i__1 = min(nevd2,*np);
654
i__2 = kplusp - nevd2 + 1, i__3 = kplusp - *np;
655
igraphdswap_(&i__1, &bounds[nevm2 + 1], &c__1, &bounds[max(i__2,
661
/* %--------------------------------------------------%
662
| LM, SM, LA, SA case. |
663
| Sort the eigenvalues of H into the an order that |
664
| is opposite to WHICH, and apply the resulting |
665
| order to BOUNDS. The eigenvalues are sorted so |
666
| that the wanted part are always within the first |
668
%--------------------------------------------------% */
670
if (s_cmp(which, "LM", (ftnlen)2, (ftnlen)2) == 0) {
671
s_copy(wprime, "SM", (ftnlen)2, (ftnlen)2);
673
if (s_cmp(which, "SM", (ftnlen)2, (ftnlen)2) == 0) {
674
s_copy(wprime, "LM", (ftnlen)2, (ftnlen)2);
676
if (s_cmp(which, "LA", (ftnlen)2, (ftnlen)2) == 0) {
677
s_copy(wprime, "SA", (ftnlen)2, (ftnlen)2);
679
if (s_cmp(which, "SA", (ftnlen)2, (ftnlen)2) == 0) {
680
s_copy(wprime, "LA", (ftnlen)2, (ftnlen)2);
683
igraphdsortr_(wprime, &c_true, &kplusp, &ritz[1], &bounds[1])
688
/* %--------------------------------------------------%
689
| Scale the Ritz estimate of each Ritz value |
690
| by 1 / max(eps23,magnitude of the Ritz value). |
691
%--------------------------------------------------% */
694
for (j = 1; j <= i__1; ++j) {
696
d__2 = eps23, d__3 = (d__1 = ritz[j], abs(d__1));
697
temp = max(d__2,d__3);
702
/* %----------------------------------------------------%
703
| Sort the Ritz values according to the scaled Ritz |
704
| esitmates. This will push all the converged ones |
705
| towards the front of ritzr, ritzi, bounds |
706
| (in the case when NCONV < NEV.) |
707
%----------------------------------------------------% */
709
s_copy(wprime, "LA", (ftnlen)2, (ftnlen)2);
710
igraphdsortr_(wprime, &c_true, &nev0, &bounds[1], &ritz[1]);
712
/* %----------------------------------------------%
713
| Scale the Ritz estimate back to its original |
715
%----------------------------------------------% */
718
for (j = 1; j <= i__1; ++j) {
720
d__2 = eps23, d__3 = (d__1 = ritz[j], abs(d__1));
721
temp = max(d__2,d__3);
726
/* %--------------------------------------------------%
727
| Sort the "converged" Ritz values again so that |
728
| the "threshold" values and their associated Ritz |
729
| estimates appear at the appropriate position in |
731
%--------------------------------------------------% */
733
if (s_cmp(which, "BE", (ftnlen)2, (ftnlen)2) == 0) {
735
/* %------------------------------------------------%
736
| Sort the "converged" Ritz values in increasing |
737
| order. The "threshold" values are in the |
739
%------------------------------------------------% */
741
s_copy(wprime, "LA", (ftnlen)2, (ftnlen)2);
742
igraphdsortr_(wprime, &c_true, &nconv, &ritz[1], &bounds[1]);
746
/* %----------------------------------------------%
747
| In LM, SM, LA, SA case, sort the "converged" |
748
| Ritz values according to WHICH so that the |
749
| "threshold" value appears at the front of |
751
%----------------------------------------------% */
752
igraphdsortr_(which, &c_true, &nconv, &ritz[1], &bounds[1]);
756
/* %------------------------------------------%
757
| Use h( 1,1 ) as storage to communicate |
758
| rnorm to _seupd if needed |
759
%------------------------------------------% */
761
h__[h_dim1 + 1] = rnorm;
764
igraphdvout_(&logfil, &kplusp, &ritz[1], &ndigit, "_saup2: Sorted Ritz"
765
" values.", (ftnlen)27);
766
igraphdvout_(&logfil, &kplusp, &bounds[1], &ndigit, "_saup2: Sorted ri"
767
"tz estimates.", (ftnlen)30);
770
/* %------------------------------------%
771
| Max iterations have been exceeded. |
772
%------------------------------------% */
774
if (iter > *mxiter && nconv < *nev) {
778
/* %---------------------%
779
| No shifts to apply. |
780
%---------------------% */
782
if (*np == 0 && nconv < nev0) {
789
} else if (nconv < *nev && *ishift == 1) {
791
/* %---------------------------------------------------%
792
| Do not have all the requested eigenvalues yet. |
793
| To prevent possible stagnation, adjust the number |
794
| of Ritz values and the shifts. |
795
%---------------------------------------------------% */
799
i__1 = nconv, i__2 = *np / 2;
800
*nev += min(i__1,i__2);
801
if (*nev == 1 && kplusp >= 6) {
803
} else if (*nev == 1 && kplusp > 2) {
808
/* %---------------------------------------%
809
| If the size of NEV was just increased |
810
| resort the eigenvalues. |
811
%---------------------------------------% */
814
igraphdsgets_(ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
820
igraphivout_(&logfil, &c__1, &nconv, &ndigit, "_saup2: no. of \"converge"
821
"d\" Ritz values at this iter.", (ftnlen)52);
825
igraphivout_(&logfil, &c__2, kp, &ndigit, "_saup2: NEV and NP are", (
827
igraphdvout_(&logfil, nev, &ritz[*np + 1], &ndigit, "_saup2: \"wante"
828
"d\" Ritz values.", (ftnlen)29);
829
igraphdvout_(&logfil, nev, &bounds[*np + 1], &ndigit, "_saup2: Ritz es"
830
"timates of the \"wanted\" values ", (ftnlen)46);
836
/* %-----------------------------------------------------%
837
| User specified shifts: reverse communication to |
838
| compute the shifts. They are returned in the first |
839
| NP locations of WORKL. |
840
%-----------------------------------------------------% */
849
/* %------------------------------------%
850
| Back from reverse communication; |
851
| User specified shifts are returned |
853
%------------------------------------% */
858
/* %---------------------------------------------------------%
859
| Move the NP shifts to the first NP locations of RITZ to |
860
| free up WORKL. This is for the non-exact shift case; |
861
| in the exact shift case, dsgets already handles this. |
862
%---------------------------------------------------------% */
865
igraphdcopy_(np, &workl[1], &c__1, &ritz[1], &c__1);
869
igraphivout_(&logfil, &c__1, np, &ndigit, "_saup2: The number of shifts to"
870
" apply ", (ftnlen)38);
871
igraphdvout_(&logfil, np, &workl[1], &ndigit, "_saup2: shifts selected", (
874
igraphdvout_(&logfil, np, &bounds[1], &ndigit, "_saup2: corresponding "
875
"Ritz estimates", (ftnlen)36);
879
/* %---------------------------------------------------------%
880
| Apply the NP0 implicit shifts by QR bulge chasing. |
881
| Each shift is applied to the entire tridiagonal matrix. |
882
| The first 2*N locations of WORKD are used as workspace. |
883
| After dsapps is done, we have a Lanczos |
884
| factorization of length NEV. |
885
%---------------------------------------------------------% */
887
igraphdsapps_(n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
888
resid[1], &q[q_offset], ldq, &workd[1]);
890
/* %---------------------------------------------%
891
| Compute the B-norm of the updated residual. |
892
| Keep B*RESID in WORKD(1:N) to be used in |
893
| the first step of the next call to dsaitr. |
894
%---------------------------------------------% */
898
if (*(unsigned char *)bmat == 'G') {
900
igraphdcopy_(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
905
/* %----------------------------------%
906
| Exit in order to compute B*RESID |
907
%----------------------------------% */
910
} else if (*(unsigned char *)bmat == 'I') {
911
igraphdcopy_(n, &resid[1], &c__1, &workd[1], &c__1);
916
/* %----------------------------------%
917
| Back from reverse communication; |
918
| WORKD(1:N) := B*RESID |
919
%----------------------------------% */
921
if (*(unsigned char *)bmat == 'G') {
926
if (*(unsigned char *)bmat == 'G') {
927
rnorm = igraphddot_(n, &resid[1], &c__1, &workd[1], &c__1);
928
rnorm = sqrt((abs(rnorm)));
929
} else if (*(unsigned char *)bmat == 'I') {
930
rnorm = igraphdnrm2_(n, &resid[1], &c__1);
936
igraphdvout_(&logfil, &c__1, &rnorm, &ndigit, "_saup2: B-norm of residual "
937
"for NEV factorization", (ftnlen)48);
938
igraphdvout_(&logfil, nev, &h__[(h_dim1 << 1) + 1], &ndigit, "_saup2: main"
939
" diagonal of compressed H matrix", (ftnlen)44);
941
igraphdvout_(&logfil, &i__1, &h__[h_dim1 + 2], &ndigit, "_saup2: subdiagon"
942
"al of compressed H matrix", (ftnlen)42);
947
/* %---------------------------------------------------------------%
949
| E N D O F M A I N I T E R A T I O N L O O P |
951
%---------------------------------------------------------------% */
975
} /* igraphdsaup2_ */
added
removed
src/lapack/dsaup2.c
