1
/* igraphdnaup2.f -- translated by f2c (version 20050501).
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
14
#include "arpack_internal.h"
17
/* Common Block Declarations */
20
integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps,
21
msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets,
22
mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd;
25
#define debug_1 debug_
28
integer nopx, nbx, nrorth, nitref, nrstrt;
29
real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd,
30
tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2,
31
tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0,
35
#define timing_1 timing_
37
/* Table of constant values */
39
static doublereal c_b3 = .66666666666666663;
40
static integer c__1 = 1;
41
static integer c__0 = 0;
42
static integer c__4 = 4;
43
static logical c_true = TRUE_;
44
static integer c__2 = 2;
48
/* \Name: igraphdnaup2 */
51
/* Intermediate level interface called by igraphdnaupd . */
54
/* call igraphdnaup2 */
55
/* ( IDO, BMAT, N, WHICH, NEV, NP, TOL, RESID, MODE, IUPD, */
56
/* ISHIFT, MXITER, V, LDV, H, LDH, RITZR, RITZI, BOUNDS, */
57
/* Q, LDQ, WORKL, IPNTR, WORKD, INFO ) */
61
/* IDO, BMAT, N, WHICH, NEV, TOL, RESID: same as defined in igraphdnaupd . */
62
/* MODE, ISHIFT, MXITER: see the definition of IPARAM in igraphdnaupd . */
64
/* NP Integer. (INPUT/OUTPUT) */
65
/* Contains the number of implicit shifts to apply during */
66
/* each Arnoldi iteration. */
67
/* If ISHIFT=1, NP is adjusted dynamically at each iteration */
68
/* to accelerate convergence and prevent stagnation. */
69
/* This is also roughly equal to the number of matrix-vector */
70
/* products (involving the operator OP) per Arnoldi iteration. */
71
/* The logic for adjusting is contained within the current */
73
/* If ISHIFT=0, NP is the number of shifts the user needs */
74
/* to provide via reverse comunication. 0 < NP < NCV-NEV. */
75
/* NP may be less than NCV-NEV for two reasons. The first, is */
76
/* to keep complex conjugate pairs of "wanted" Ritz values */
77
/* together. The second, is that a leading block of the current */
78
/* upper Hessenberg matrix has split off and contains "unwanted" */
80
/* Upon termination of the IRA iteration, NP contains the number */
81
/* of "converged" wanted Ritz values. */
83
/* IUPD Integer. (INPUT) */
84
/* IUPD .EQ. 0: use explicit restart instead implicit update. */
85
/* IUPD .NE. 0: use implicit update. */
87
/* V Double precision N by (NEV+NP) array. (INPUT/OUTPUT) */
88
/* The Arnoldi basis vectors are returned in the first NEV */
91
/* LDV Integer. (INPUT) */
92
/* Leading dimension of V exactly as declared in the calling */
95
/* H Double precision (NEV+NP) by (NEV+NP) array. (OUTPUT) */
96
/* H is used to store the generated upper Hessenberg matrix */
98
/* LDH Integer. (INPUT) */
99
/* Leading dimension of H exactly as declared in the calling */
102
/* RITZR, Double precision arrays of length NEV+NP. (OUTPUT) */
103
/* RITZI RITZR(1:NEV) (resp. RITZI(1:NEV)) contains the real (resp. */
104
/* imaginary) part of the computed Ritz values of OP. */
106
/* BOUNDS Double precision array of length NEV+NP. (OUTPUT) */
107
/* BOUNDS(1:NEV) contain the error bounds corresponding to */
108
/* the computed Ritz values. */
110
/* Q Double precision (NEV+NP) by (NEV+NP) array. (WORKSPACE) */
111
/* Private (replicated) work array used to accumulate the */
112
/* rotation in the shift application step. */
114
/* LDQ Integer. (INPUT) */
115
/* Leading dimension of Q exactly as declared in the calling */
118
/* WORKL Double precision work array of length at least */
119
/* (NEV+NP)**2 + 3*(NEV+NP). (INPUT/WORKSPACE) */
120
/* Private (replicated) array on each PE or array allocated on */
121
/* the front end. It is used in shifts calculation, shifts */
122
/* application and convergence checking. */
124
/* On exit, the last 3*(NEV+NP) locations of WORKL contain */
125
/* the Ritz values (real,imaginary) and associated Ritz */
126
/* estimates of the current Hessenberg matrix. They are */
127
/* listed in the same order as returned from dneigh . */
129
/* If ISHIFT .EQ. O and IDO .EQ. 3, the first 2*NP locations */
130
/* of WORKL are used in reverse communication to hold the user */
131
/* supplied shifts. */
133
/* IPNTR Integer array of length 3. (OUTPUT) */
134
/* Pointer to mark the starting locations in the WORKD for */
135
/* vectors used by the Arnoldi iteration. */
136
/* ------------------------------------------------------------- */
137
/* IPNTR(1): pointer to the current operand vector X. */
138
/* IPNTR(2): pointer to the current result vector Y. */
139
/* IPNTR(3): pointer to the vector B * X when used in the */
140
/* shift-and-invert mode. X is the current operand. */
141
/* ------------------------------------------------------------- */
143
/* WORKD Double precision work array of length 3*N. (WORKSPACE) */
144
/* Distributed array to be used in the basic Arnoldi iteration */
145
/* for reverse communication. The user should not use WORKD */
146
/* as temporary workspace during the iteration !!!!!!!!!! */
147
/* See Data Distribution Note in DNAUPD. */
149
/* INFO Integer. (INPUT/OUTPUT) */
150
/* If INFO .EQ. 0, a randomly initial residual vector is used. */
151
/* If INFO .NE. 0, RESID contains the initial residual vector, */
152
/* possibly from a previous run. */
153
/* Error flag on output. */
154
/* = 0: Normal return. */
155
/* = 1: Maximum number of iterations taken. */
156
/* All possible eigenvalues of OP has been found. */
157
/* NP returns the number of converged Ritz values. */
158
/* = 2: No shifts could be applied. */
159
/* = -8: Error return from LAPACK eigenvalue calculation; */
160
/* This should never happen. */
161
/* = -9: Starting vector is zero. */
162
/* = -9999: Could not build an Arnoldi factorization. */
163
/* Size that was built in returned in NP. */
167
/* ----------------------------------------------------------------------- */
171
/* \Local variables: */
175
/* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
176
/* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
178
/* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
179
/* Restarted Arnoldi Iteration", Rice University Technical Report */
180
/* TR95-13, Department of Computational and Applied Mathematics. */
182
/* \Routines called: */
183
/* dgetv0 ARPACK initial vector generation routine. */
184
/* igraphdnaitr ARPACK Arnoldi factorization routine. */
185
/* igraphdnapps ARPACK application of implicit shifts routine. */
186
/* igraphdnconv ARPACK convergence of Ritz values routine. */
187
/* dneigh ARPACK compute Ritz values and error bounds routine. */
188
/* dngets ARPACK reorder Ritz values and error bounds routine. */
189
/* dsortc ARPACK sorting routine. */
190
/* ivout ARPACK utility routine that prints integers. */
191
/* second ARPACK utility routine for timing. */
192
/* dmout ARPACK utility routine that prints matrices */
193
/* dvout ARPACK utility routine that prints vectors. */
194
/* dlamch LAPACK routine that determines machine constants. */
195
/* dlapy2 LAPACK routine to compute sqrt(x**2+y**2) carefully. */
196
/* dcopy Level 1 BLAS that copies one vector to another . */
197
/* ddot Level 1 BLAS that computes the scalar product of two vectors. */
198
/* dnrm2 Level 1 BLAS that computes the norm of a vector. */
199
/* dswap Level 1 BLAS that swaps two vectors. */
202
/* Danny Sorensen Phuong Vu */
203
/* Richard Lehoucq CRPC / Rice University */
204
/* Dept. of Computational & Houston, Texas */
205
/* Applied Mathematics */
206
/* Rice University */
209
/* \SCCS Information: @(#) */
210
/* FILE: naup2.F SID: 2.8 DATE OF SID: 10/17/00 RELEASE: 2 */
217
/* ----------------------------------------------------------------------- */
219
/* Subroutine */ int igraphdnaup2_(integer *ido, char *bmat, integer *n, char *
220
which, integer *nev, integer *np, doublereal *tol, doublereal *resid,
221
integer *mode, integer *iupd, integer *ishift, integer *mxiter,
222
doublereal *v, integer *ldv, doublereal *h__, integer *ldh,
223
doublereal *ritzr, doublereal *ritzi, doublereal *bounds, doublereal *
224
q, integer *ldq, doublereal *workl, integer *ipntr, doublereal *workd,
227
/* System generated locals */
228
integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2;
229
doublereal d__1, d__2;
231
/* Builtin functions */
232
double igraphpow_dd(doublereal *, doublereal *);
233
integer igraphs_cmp(char *, char *, ftnlen, ftnlen);
234
/* Subroutine */ int igraphs_copy(char *, char *, ftnlen, ftnlen);
235
double sqrt(doublereal);
237
/* Local variables */
239
static real t0, t1, t2, t3;
240
static integer kp[4], np0, nev0;
241
extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *,
243
static doublereal eps23;
244
static integer ierr, iter;
245
static doublereal temp;
246
extern doublereal igraphdnrm2_(integer *, doublereal *, integer *);
247
static logical getv0, cnorm;
248
extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *,
249
doublereal *, integer *);
250
static integer nconv;
251
extern /* Subroutine */ int igraphdmout_(integer *, integer *, integer *,
252
doublereal *, integer *, integer *, char *);
253
static logical initv;
254
static doublereal rnorm;
255
extern /* Subroutine */ int igraphdvout_(integer *, integer *, doublereal *,
256
integer *, char *), igraphivout_(integer *, integer *, integer *
257
, integer *, char *), igraphdgetv0_(integer *, char *, integer *
258
, logical *, integer *, integer *, doublereal *, integer *,
259
doublereal *, doublereal *, integer *, doublereal *, integer *);
260
extern doublereal igraphdlapy2_(doublereal *, doublereal *), igraphdlamch_(char *);
261
extern /* Subroutine */ int igraphdneigh_(doublereal *, integer *, doublereal *,
262
integer *, doublereal *, doublereal *, doublereal *, doublereal *
263
, integer *, doublereal *, integer *);
264
static integer nevbef;
265
extern /* Subroutine */ int igraphsecond_(real *);
266
static logical update;
267
static char wprime[2];
268
static logical ushift;
269
static integer kplusp, msglvl, nptemp, numcnv;
270
extern /* Subroutine */ int igraphdnaitr_(integer *, char *, integer *, integer
271
*, integer *, integer *, doublereal *, doublereal *, doublereal *,
272
integer *, doublereal *, integer *, integer *, doublereal *,
273
integer *), igraphdnconv_(integer *, doublereal *, doublereal *,
274
doublereal *, doublereal *, integer *), igraphdngets_(integer *, char *
275
, integer *, integer *, doublereal *, doublereal *, doublereal *,
276
doublereal *, doublereal *), igraphdnapps_(integer *, integer *,
277
integer *, doublereal *, doublereal *, doublereal *, integer *,
278
doublereal *, integer *, doublereal *, doublereal *, integer *,
279
doublereal *, doublereal *), igraphdsortc_(char *, logical *, integer *,
280
doublereal *, doublereal *, doublereal *);
283
/* %----------------------------------------------------% */
284
/* | Include files for debugging and timing information | */
285
/* %----------------------------------------------------% */
288
/* \SCCS Information: @(#) */
289
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
291
/* %---------------------------------% */
292
/* | See debug.doc for documentation | */
293
/* %---------------------------------% */
295
/* %------------------% */
296
/* | Scalar Arguments | */
297
/* %------------------% */
299
/* %--------------------------------% */
300
/* | See stat.doc for documentation | */
301
/* %--------------------------------% */
303
/* \SCCS Information: @(#) */
304
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
308
/* %-----------------% */
309
/* | Array Arguments | */
310
/* %-----------------% */
318
/* %---------------% */
319
/* | Local Scalars | */
320
/* %---------------% */
323
/* %-----------------------% */
324
/* | Local array arguments | */
325
/* %-----------------------% */
328
/* %----------------------% */
329
/* | External Subroutines | */
330
/* %----------------------% */
333
/* %--------------------% */
334
/* | External Functions | */
335
/* %--------------------% */
338
/* %---------------------% */
339
/* | Intrinsic Functions | */
340
/* %---------------------% */
343
/* %-----------------------% */
344
/* | Executable Statements | */
345
/* %-----------------------% */
347
/* Parameter adjustments */
355
v_offset = 1 + v_dim1;
358
h_offset = 1 + h_dim1;
361
q_offset = 1 + q_dim1;
370
msglvl = debug_1.mnaup2;
372
/* %-------------------------------------% */
373
/* | Get the machine dependent constant. | */
374
/* %-------------------------------------% */
376
eps23 = igraphdlamch_("Epsilon-Machine");
377
eps23 = igraphpow_dd(&eps23, &c_b3);
382
/* %-------------------------------------% */
383
/* | kplusp is the bound on the largest | */
384
/* | Lanczos factorization built. | */
385
/* | nconv is the current number of | */
386
/* | "converged" eigenvlues. | */
387
/* | iter is the counter on the current | */
388
/* | iteration step. | */
389
/* %-------------------------------------% */
395
/* %---------------------------------------% */
396
/* | Set flags for computing the first NEV | */
397
/* | steps of the Arnoldi factorization. | */
398
/* %---------------------------------------% */
407
/* %--------------------------------------------% */
408
/* | User provides the initial residual vector. | */
409
/* %--------------------------------------------% */
418
/* %---------------------------------------------% */
419
/* | Get a possibly random starting vector and | */
420
/* | force it into the range of the operator OP. | */
421
/* %---------------------------------------------% */
426
igraphdgetv0_(ido, bmat, &c__1, &initv, n, &c__1, &v[v_offset], ldv, &resid[
427
1], &rnorm, &ipntr[1], &workd[1], info);
435
/* %-----------------------------------------% */
436
/* | The initial vector is zero. Error exit. | */
437
/* %-----------------------------------------% */
446
/* %-----------------------------------% */
447
/* | Back from reverse communication : | */
448
/* | continue with update step | */
449
/* %-----------------------------------% */
455
/* %-------------------------------------------% */
456
/* | Back from computing user specified shifts | */
457
/* %-------------------------------------------% */
463
/* %-------------------------------------% */
464
/* | Back from computing residual norm | */
465
/* | at the end of the current iteration | */
466
/* %-------------------------------------% */
472
/* %----------------------------------------------------------% */
473
/* | Compute the first NEV steps of the Arnoldi factorization | */
474
/* %----------------------------------------------------------% */
476
igraphdnaitr_(ido, bmat, n, &c__0, nev, mode, &resid[1], &rnorm, &v[v_offset],
477
ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], info);
479
/* %---------------------------------------------------% */
480
/* | ido .ne. 99 implies use of reverse communication | */
481
/* | to compute operations involving OP and possibly B | */
482
/* %---------------------------------------------------% */
495
/* %--------------------------------------------------------------% */
497
/* | M A I N ARNOLDI I T E R A T I O N L O O P | */
498
/* | Each iteration implicitly restarts the Arnoldi | */
499
/* | factorization in place. | */
501
/* %--------------------------------------------------------------% */
508
igraphivout_(&debug_1.logfil, &c__1, &iter, &debug_1.ndigit, "_naup2: ****"
509
" Start of major iteration number ****");
512
/* %-----------------------------------------------------------% */
513
/* | Compute NP additional steps of the Arnoldi factorization. | */
514
/* | Adjust NP since NEV might have been updated by last call | */
515
/* | to the shift application routine igraphdnapps . | */
516
/* %-----------------------------------------------------------% */
521
igraphivout_(&debug_1.logfil, &c__1, nev, &debug_1.ndigit, "_naup2: The le"
522
"ngth of the current Arnoldi factorization");
523
igraphivout_(&debug_1.logfil, &c__1, np, &debug_1.ndigit, "_naup2: Extend "
524
"the Arnoldi factorization by");
527
/* %-----------------------------------------------------------% */
528
/* | Compute NP additional steps of the Arnoldi factorization. | */
529
/* %-----------------------------------------------------------% */
535
igraphdnaitr_(ido, bmat, n, nev, np, mode, &resid[1], &rnorm, &v[v_offset], ldv,
536
&h__[h_offset], ldh, &ipntr[1], &workd[1], info);
538
/* %---------------------------------------------------% */
539
/* | ido .ne. 99 implies use of reverse communication | */
540
/* | to compute operations involving OP and possibly B | */
541
/* %---------------------------------------------------% */
556
igraphdvout_(&debug_1.logfil, &c__1, &rnorm, &debug_1.ndigit, "_naup2: Cor"
557
"responding B-norm of the residual");
560
/* %--------------------------------------------------------% */
561
/* | Compute the eigenvalues and corresponding error bounds | */
562
/* | of the current upper Hessenberg matrix. | */
563
/* %--------------------------------------------------------% */
565
igraphdneigh_(&rnorm, &kplusp, &h__[h_offset], ldh, &ritzr[1], &ritzi[1], &
566
bounds[1], &q[q_offset], ldq, &workl[1], &ierr);
573
/* %----------------------------------------------------% */
574
/* | Make a copy of eigenvalues and corresponding error | */
575
/* | bounds obtained from dneigh . | */
576
/* %----------------------------------------------------% */
578
/* Computing 2nd power */
580
igraphdcopy_(&kplusp, &ritzr[1], &c__1, &workl[i__1 * i__1 + 1], &c__1);
581
/* Computing 2nd power */
583
igraphdcopy_(&kplusp, &ritzi[1], &c__1, &workl[i__1 * i__1 + kplusp + 1], &c__1)
585
/* Computing 2nd power */
587
igraphdcopy_(&kplusp, &bounds[1], &c__1, &workl[i__1 * i__1 + (kplusp << 1) + 1]
590
/* %---------------------------------------------------% */
591
/* | Select the wanted Ritz values and their bounds | */
592
/* | to be used in the convergence test. | */
593
/* | The wanted part of the spectrum and corresponding | */
594
/* | error bounds are in the last NEV loc. of RITZR, | */
595
/* | RITZI and BOUNDS respectively. The variables NEV | */
596
/* | and NP may be updated if the NEV-th wanted Ritz | */
597
/* | value has a non zero imaginary part. In this case | */
598
/* | NEV is increased by one and NP decreased by one. | */
599
/* | NOTE: The last two arguments of dngets are no | */
600
/* | longer used as of version 2.1. | */
601
/* %---------------------------------------------------% */
606
igraphdngets_(ishift, which, nev, np, &ritzr[1], &ritzi[1], &bounds[1], &workl[
607
1], &workl[*np + 1]);
608
if (*nev == nev0 + 1) {
612
/* %-------------------% */
613
/* | Convergence test. | */
614
/* %-------------------% */
616
igraphdcopy_(nev, &bounds[*np + 1], &c__1, &workl[(*np << 1) + 1], &c__1);
617
igraphdnconv_(nev, &ritzr[*np + 1], &ritzi[*np + 1], &workl[(*np << 1) + 1],
625
igraphivout_(&debug_1.logfil, &c__4, kp, &debug_1.ndigit, "_naup2: NEV, NP"
626
", NUMCNV, NCONV are");
627
igraphdvout_(&debug_1.logfil, &kplusp, &ritzr[1], &debug_1.ndigit, "_naup2"
628
": Real part of the eigenvalues of H");
629
igraphdvout_(&debug_1.logfil, &kplusp, &ritzi[1], &debug_1.ndigit, "_naup2"
630
": Imaginary part of the eigenvalues of H");
631
igraphdvout_(&debug_1.logfil, &kplusp, &bounds[1], &debug_1.ndigit, "_naup"
632
"2: Ritz estimates of the current NCV Ritz values")
636
/* %---------------------------------------------------------% */
637
/* | Count the number of unwanted Ritz values that have zero | */
638
/* | Ritz estimates. If any Ritz estimates are equal to zero | */
639
/* | then a leading block of H of order equal to at least | */
640
/* | the number of Ritz values with zero Ritz estimates has | */
641
/* | split off. None of these Ritz values may be removed by | */
642
/* | shifting. Decrease NP the number of shifts to apply. If | */
643
/* | no shifts may be applied, then prepare to exit | */
644
/* %---------------------------------------------------------% */
648
for (j = 1; j <= i__1; ++j) {
649
if (bounds[j] == 0.) {
656
if (nconv >= numcnv || iter > *mxiter || *np == 0) {
659
/* Computing 2nd power */
661
igraphdvout_(&debug_1.logfil, &kplusp, &workl[i__1 * i__1 + 1], &
662
debug_1.ndigit, "_naup2: Real part of the eig computed b"
664
/* Computing 2nd power */
666
igraphdvout_(&debug_1.logfil, &kplusp, &workl[i__1 * i__1 + kplusp + 1],
667
&debug_1.ndigit, "_naup2: Imag part of the eig computed"
669
/* Computing 2nd power */
671
igraphdvout_(&debug_1.logfil, &kplusp, &workl[i__1 * i__1 + (kplusp <<
672
1) + 1], &debug_1.ndigit, "_naup2: Ritz eistmates comput"
676
/* %------------------------------------------------% */
677
/* | Prepare to exit. Put the converged Ritz values | */
678
/* | and corresponding bounds in RITZ(1:NCONV) and | */
679
/* | BOUNDS(1:NCONV) respectively. Then sort. Be | */
680
/* | careful when NCONV > NP | */
681
/* %------------------------------------------------% */
683
/* %------------------------------------------% */
684
/* | Use h( 3,1 ) as storage to communicate | */
685
/* | rnorm to _neupd if needed | */
686
/* %------------------------------------------% */
687
h__[h_dim1 + 3] = rnorm;
689
/* %----------------------------------------------% */
690
/* | To be consistent with dngets , we first do a | */
691
/* | pre-processing sort in order to keep complex | */
692
/* | conjugate pairs together. This is similar | */
693
/* | to the pre-processing sort used in dngets | */
694
/* | except that the sort is done in the opposite | */
696
/* %----------------------------------------------% */
698
if (igraphs_cmp(which, "LM", (ftnlen)2, (ftnlen)2) == 0) {
699
igraphs_copy(wprime, "SR", (ftnlen)2, (ftnlen)2);
701
if (igraphs_cmp(which, "SM", (ftnlen)2, (ftnlen)2) == 0) {
702
igraphs_copy(wprime, "LR", (ftnlen)2, (ftnlen)2);
704
if (igraphs_cmp(which, "LR", (ftnlen)2, (ftnlen)2) == 0) {
705
igraphs_copy(wprime, "SM", (ftnlen)2, (ftnlen)2);
707
if (igraphs_cmp(which, "SR", (ftnlen)2, (ftnlen)2) == 0) {
708
igraphs_copy(wprime, "LM", (ftnlen)2, (ftnlen)2);
710
if (igraphs_cmp(which, "LI", (ftnlen)2, (ftnlen)2) == 0) {
711
igraphs_copy(wprime, "SM", (ftnlen)2, (ftnlen)2);
713
if (igraphs_cmp(which, "SI", (ftnlen)2, (ftnlen)2) == 0) {
714
igraphs_copy(wprime, "LM", (ftnlen)2, (ftnlen)2);
717
igraphdsortc_(wprime, &c_true, &kplusp, &ritzr[1], &ritzi[1], &bounds[1]);
719
/* %----------------------------------------------% */
720
/* | Now sort Ritz values so that converged Ritz | */
721
/* | values appear within the first NEV locations | */
722
/* | of ritzr, ritzi and bounds, and the most | */
723
/* | desired one appears at the front. | */
724
/* %----------------------------------------------% */
726
if (igraphs_cmp(which, "LM", (ftnlen)2, (ftnlen)2) == 0) {
727
igraphs_copy(wprime, "SM", (ftnlen)2, (ftnlen)2);
729
if (igraphs_cmp(which, "SM", (ftnlen)2, (ftnlen)2) == 0) {
730
igraphs_copy(wprime, "LM", (ftnlen)2, (ftnlen)2);
732
if (igraphs_cmp(which, "LR", (ftnlen)2, (ftnlen)2) == 0) {
733
igraphs_copy(wprime, "SR", (ftnlen)2, (ftnlen)2);
735
if (igraphs_cmp(which, "SR", (ftnlen)2, (ftnlen)2) == 0) {
736
igraphs_copy(wprime, "LR", (ftnlen)2, (ftnlen)2);
738
if (igraphs_cmp(which, "LI", (ftnlen)2, (ftnlen)2) == 0) {
739
igraphs_copy(wprime, "SI", (ftnlen)2, (ftnlen)2);
741
if (igraphs_cmp(which, "SI", (ftnlen)2, (ftnlen)2) == 0) {
742
igraphs_copy(wprime, "LI", (ftnlen)2, (ftnlen)2);
745
igraphdsortc_(wprime, &c_true, &kplusp, &ritzr[1], &ritzi[1], &bounds[1]);
747
/* %--------------------------------------------------% */
748
/* | Scale the Ritz estimate of each Ritz value | */
749
/* | by 1 / max(eps23,magnitude of the Ritz value). | */
750
/* %--------------------------------------------------% */
753
for (j = 1; j <= i__1; ++j) {
755
d__1 = eps23, d__2 = igraphdlapy2_(&ritzr[j], &ritzi[j]);
756
temp = max(d__1,d__2);
761
/* %----------------------------------------------------% */
762
/* | Sort the Ritz values according to the scaled Ritz | */
763
/* | esitmates. This will push all the converged ones | */
764
/* | towards the front of ritzr, ritzi, bounds | */
765
/* | (in the case when NCONV < NEV.) | */
766
/* %----------------------------------------------------% */
768
igraphs_copy(wprime, "LR", (ftnlen)2, (ftnlen)2);
769
igraphdsortc_(wprime, &c_true, &numcnv, &bounds[1], &ritzr[1], &ritzi[1]);
771
/* %----------------------------------------------% */
772
/* | Scale the Ritz estimate back to its original | */
774
/* %----------------------------------------------% */
777
for (j = 1; j <= i__1; ++j) {
779
d__1 = eps23, d__2 = igraphdlapy2_(&ritzr[j], &ritzi[j]);
780
temp = max(d__1,d__2);
785
/* %------------------------------------------------% */
786
/* | Sort the converged Ritz values again so that | */
787
/* | the "threshold" value appears at the front of | */
788
/* | ritzr, ritzi and bound. | */
789
/* %------------------------------------------------% */
791
igraphdsortc_(which, &c_true, &nconv, &ritzr[1], &ritzi[1], &bounds[1]);
794
igraphdvout_(&debug_1.logfil, &kplusp, &ritzr[1], &debug_1.ndigit,
795
"_naup2: Sorted real part of the eigenvalues");
796
igraphdvout_(&debug_1.logfil, &kplusp, &ritzi[1], &debug_1.ndigit,
797
"_naup2: Sorted imaginary part of the eigenvalues");
798
igraphdvout_(&debug_1.logfil, &kplusp, &bounds[1], &debug_1.ndigit,
799
"_naup2: Sorted ritz estimates.");
802
/* %------------------------------------% */
803
/* | Max iterations have been exceeded. | */
804
/* %------------------------------------% */
806
if (iter > *mxiter && nconv < numcnv) {
810
/* %---------------------% */
811
/* | No shifts to apply. | */
812
/* %---------------------% */
814
if (*np == 0 && nconv < numcnv) {
821
} else if (nconv < numcnv && *ishift == 1) {
823
/* %-------------------------------------------------% */
824
/* | Do not have all the requested eigenvalues yet. | */
825
/* | To prevent possible stagnation, adjust the size | */
827
/* %-------------------------------------------------% */
831
i__1 = nconv, i__2 = *np / 2;
832
*nev += min(i__1,i__2);
833
if (*nev == 1 && kplusp >= 6) {
835
} else if (*nev == 1 && kplusp > 3) {
840
/* %---------------------------------------% */
841
/* | If the size of NEV was just increased | */
842
/* | resort the eigenvalues. | */
843
/* %---------------------------------------% */
846
igraphdngets_(ishift, which, nev, np, &ritzr[1], &ritzi[1], &bounds[1],
847
&workl[1], &workl[*np + 1]);
853
igraphivout_(&debug_1.logfil, &c__1, &nconv, &debug_1.ndigit, "_naup2: no."
854
" of \"converged\" Ritz values at this iter.");
858
igraphivout_(&debug_1.logfil, &c__2, kp, &debug_1.ndigit, "_naup2: NEV"
860
igraphdvout_(&debug_1.logfil, nev, &ritzr[*np + 1], &debug_1.ndigit,
861
"_naup2: \"wanted\" Ritz values -- real part");
862
igraphdvout_(&debug_1.logfil, nev, &ritzi[*np + 1], &debug_1.ndigit,
863
"_naup2: \"wanted\" Ritz values -- imag part")
865
igraphdvout_(&debug_1.logfil, nev, &bounds[*np + 1], &debug_1.ndigit,
866
"_naup2: Ritz estimates of the \"wanted\" values ");
872
/* %-------------------------------------------------------% */
873
/* | User specified shifts: reverse comminucation to | */
874
/* | compute the shifts. They are returned in the first | */
875
/* | 2*NP locations of WORKL. | */
876
/* %-------------------------------------------------------% */
885
/* %------------------------------------% */
886
/* | Back from reverse communication; | */
887
/* | User specified shifts are returned | */
888
/* | in WORKL(1:2*NP) | */
889
/* %------------------------------------% */
895
/* %----------------------------------% */
896
/* | Move the NP shifts from WORKL to | */
897
/* | RITZR, RITZI to free up WORKL | */
898
/* | for non-exact shift case. | */
899
/* %----------------------------------% */
901
igraphdcopy_(np, &workl[1], &c__1, &ritzr[1], &c__1);
902
igraphdcopy_(np, &workl[*np + 1], &c__1, &ritzi[1], &c__1);
906
igraphivout_(&debug_1.logfil, &c__1, np, &debug_1.ndigit, "_naup2: The num"
907
"ber of shifts to apply ");
908
igraphdvout_(&debug_1.logfil, np, &ritzr[1], &debug_1.ndigit, "_naup2: Rea"
909
"l part of the shifts");
910
igraphdvout_(&debug_1.logfil, np, &ritzi[1], &debug_1.ndigit, "_naup2: Ima"
911
"ginary part of the shifts");
913
igraphdvout_(&debug_1.logfil, np, &bounds[1], &debug_1.ndigit, "_naup2"
914
": Ritz estimates of the shifts");
918
/* %---------------------------------------------------------% */
919
/* | Apply the NP implicit shifts by QR bulge chasing. | */
920
/* | Each shift is applied to the whole upper Hessenberg | */
922
/* | The first 2*N locations of WORKD are used as workspace. | */
923
/* %---------------------------------------------------------% */
925
igraphdnapps_(n, nev, np, &ritzr[1], &ritzi[1], &v[v_offset], ldv, &h__[
926
h_offset], ldh, &resid[1], &q[q_offset], ldq, &workl[1], &workd[1]
929
/* %---------------------------------------------% */
930
/* | Compute the B-norm of the updated residual. | */
931
/* | Keep B*RESID in WORKD(1:N) to be used in | */
932
/* | the first step of the next call to igraphdnaitr . | */
933
/* %---------------------------------------------% */
937
if (*(unsigned char *)bmat == 'G') {
939
igraphdcopy_(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
944
/* %----------------------------------% */
945
/* | Exit in order to compute B*RESID | */
946
/* %----------------------------------% */
949
} else if (*(unsigned char *)bmat == 'I') {
950
igraphdcopy_(n, &resid[1], &c__1, &workd[1], &c__1);
955
/* %----------------------------------% */
956
/* | Back from reverse communication; | */
957
/* | WORKD(1:N) := B*RESID | */
958
/* %----------------------------------% */
960
if (*(unsigned char *)bmat == 'G') {
962
timing_1.tmvbx += t3 - t2;
965
if (*(unsigned char *)bmat == 'G') {
966
rnorm = igraphddot_(n, &resid[1], &c__1, &workd[1], &c__1);
967
rnorm = sqrt((abs(rnorm)));
968
} else if (*(unsigned char *)bmat == 'I') {
969
rnorm = igraphdnrm2_(n, &resid[1], &c__1);
974
igraphdvout_(&debug_1.logfil, &c__1, &rnorm, &debug_1.ndigit, "_naup2: B-n"
975
"orm of residual for compressed factorization");
976
igraphdmout_(&debug_1.logfil, nev, nev, &h__[h_offset], ldh, &
977
debug_1.ndigit, "_naup2: Compressed upper Hessenberg matrix H");
982
/* %---------------------------------------------------------------% */
984
/* | E N D O F M A I N I T E R A T I O N L O O P | */
986
/* %---------------------------------------------------------------% */
1001
timing_1.tnaup2 = t1 - t0;
1005
/* %---------------% */
1006
/* | End of igraphdnaup2 | */
1007
/* %---------------% */
1010
} /* igraphdnaup2_ */
added
removed
src/arpack/dnaup2.c
