1
/* -- 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 integer c__1 = 1;
41
/* ----------------------------------------------------------------------- */
48
/* Reverse communication interface for the Implicitly Restarted Arnoldi */
49
/* Iteration. For symmetric problems this reduces to a variant of the Lanczos */
50
/* method. This method has been designed to compute approximations to a */
51
/* few eigenpairs of a linear operator OP that is real and symmetric */
52
/* with respect to a real positive semi-definite symmetric matrix B, */
57
/* Another way to express this condition is */
59
/* < x,OPy > = < OPx,y > where < z,w > = z`Bw . */
61
/* In the standard eigenproblem B is the identity matrix. */
62
/* ( A` denotes transpose of A) */
64
/* The computed approximate eigenvalues are called Ritz values and */
65
/* the corresponding approximate eigenvectors are called Ritz vectors. */
67
/* dsaupd is usually called iteratively to solve one of the */
68
/* following problems: */
70
/* Mode 1: A*x = lambda*x, A symmetric */
71
/* ===> OP = A and B = I. */
73
/* Mode 2: A*x = lambda*M*x, A symmetric, M symmetric positive definite */
74
/* ===> OP = inv[M]*A and B = M. */
75
/* ===> (If M can be factored see remark 3 below) */
77
/* Mode 3: K*x = lambda*M*x, K symmetric, M symmetric positive semi-definite */
78
/* ===> OP = (inv[K - sigma*M])*M and B = M. */
79
/* ===> Shift-and-Invert mode */
81
/* Mode 4: K*x = lambda*KG*x, K symmetric positive semi-definite, */
82
/* KG symmetric indefinite */
83
/* ===> OP = (inv[K - sigma*KG])*K and B = K. */
84
/* ===> Buckling mode */
86
/* Mode 5: A*x = lambda*M*x, A symmetric, M symmetric positive semi-definite */
87
/* ===> OP = inv[A - sigma*M]*[A + sigma*M] and B = M. */
88
/* ===> Cayley transformed mode */
90
/* NOTE: The action of w <- inv[A - sigma*M]*v or w <- inv[M]*v */
91
/* should be accomplished either by a direct method */
92
/* using a sparse matrix factorization and solving */
94
/* [A - sigma*M]*w = v or M*w = v, */
96
/* or through an iterative method for solving these */
97
/* systems. If an iterative method is used, the */
98
/* convergence test must be more stringent than */
99
/* the accuracy requirements for the eigenvalue */
100
/* approximations. */
104
/* ( IDO, BMAT, N, WHICH, NEV, TOL, RESID, NCV, V, LDV, IPARAM, */
105
/* IPNTR, WORKD, WORKL, LWORKL, INFO ) */
108
/* IDO Integer. (INPUT/OUTPUT) */
109
/* Reverse communication flag. IDO must be zero on the first */
110
/* call to dsaupd . IDO will be set internally to */
111
/* indicate the type of operation to be performed. Control is */
112
/* then given back to the calling routine which has the */
113
/* responsibility to carry out the requested operation and call */
114
/* dsaupd with the result. The operand is given in */
115
/* WORKD(IPNTR(1)), the result must be put in WORKD(IPNTR(2)). */
116
/* (If Mode = 2 see remark 5 below) */
117
/* ------------------------------------------------------------- */
118
/* IDO = 0: first call to the reverse communication interface */
119
/* IDO = -1: compute Y = OP * X where */
120
/* IPNTR(1) is the pointer into WORKD for X, */
121
/* IPNTR(2) is the pointer into WORKD for Y. */
122
/* This is for the initialization phase to force the */
123
/* starting vector into the range of OP. */
124
/* IDO = 1: compute Y = OP * X where */
125
/* IPNTR(1) is the pointer into WORKD for X, */
126
/* IPNTR(2) is the pointer into WORKD for Y. */
127
/* In mode 3,4 and 5, the vector B * X is already */
128
/* available in WORKD(ipntr(3)). It does not */
129
/* need to be recomputed in forming OP * X. */
130
/* IDO = 2: compute Y = B * X where */
131
/* IPNTR(1) is the pointer into WORKD for X, */
132
/* IPNTR(2) is the pointer into WORKD for Y. */
133
/* IDO = 3: compute the IPARAM(8) shifts where */
134
/* IPNTR(11) is the pointer into WORKL for */
135
/* placing the shifts. See remark 6 below. */
137
/* ------------------------------------------------------------- */
139
/* BMAT Character*1. (INPUT) */
140
/* BMAT specifies the type of the matrix B that defines the */
141
/* semi-inner product for the operator OP. */
142
/* B = 'I' -> standard eigenvalue problem A*x = lambda*x */
143
/* B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x */
145
/* N Integer. (INPUT) */
146
/* Dimension of the eigenproblem. */
148
/* WHICH Character*2. (INPUT) */
149
/* Specify which of the Ritz values of OP to compute. */
151
/* 'LA' - compute the NEV largest (algebraic) eigenvalues. */
152
/* 'SA' - compute the NEV smallest (algebraic) eigenvalues. */
153
/* 'LM' - compute the NEV largest (in magnitude) eigenvalues. */
154
/* 'SM' - compute the NEV smallest (in magnitude) eigenvalues. */
155
/* 'BE' - compute NEV eigenvalues, half from each end of the */
156
/* spectrum. When NEV is odd, compute one more from the */
157
/* high end than from the low end. */
158
/* (see remark 1 below) */
160
/* NEV Integer. (INPUT) */
161
/* Number of eigenvalues of OP to be computed. 0 < NEV < N. */
163
/* TOL Double precision scalar. (INPUT) */
164
/* Stopping criterion: the relative accuracy of the Ritz value */
165
/* is considered acceptable if BOUNDS(I) .LE. TOL*ABS(RITZ(I)). */
166
/* If TOL .LE. 0. is passed a default is set: */
167
/* DEFAULT = DLAMCH ('EPS') (machine precision as computed */
168
/* by the LAPACK auxiliary subroutine DLAMCH ). */
170
/* RESID Double precision array of length N. (INPUT/OUTPUT) */
172
/* If INFO .EQ. 0, a random initial residual vector is used. */
173
/* If INFO .NE. 0, RESID contains the initial residual vector, */
174
/* possibly from a previous run. */
176
/* RESID contains the final residual vector. */
178
/* NCV Integer. (INPUT) */
179
/* Number of columns of the matrix V (less than or equal to N). */
180
/* This will indicate how many Lanczos vectors are generated */
181
/* at each iteration. After the startup phase in which NEV */
182
/* Lanczos vectors are generated, the algorithm generates */
183
/* NCV-NEV Lanczos vectors at each subsequent update iteration. */
184
/* Most of the cost in generating each Lanczos vector is in the */
185
/* matrix-vector product OP*x. (See remark 4 below). */
187
/* V Double precision N by NCV array. (OUTPUT) */
188
/* The NCV columns of V contain the Lanczos basis vectors. */
190
/* LDV Integer. (INPUT) */
191
/* Leading dimension of V exactly as declared in the calling */
194
/* IPARAM Integer array of length 11. (INPUT/OUTPUT) */
195
/* IPARAM(1) = ISHIFT: method for selecting the implicit shifts. */
196
/* The shifts selected at each iteration are used to restart */
197
/* the Arnoldi iteration in an implicit fashion. */
198
/* ------------------------------------------------------------- */
199
/* ISHIFT = 0: the shifts are provided by the user via */
200
/* reverse communication. The NCV eigenvalues of */
201
/* the current tridiagonal matrix T are returned in */
202
/* the part of WORKL array corresponding to RITZ. */
203
/* See remark 6 below. */
204
/* ISHIFT = 1: exact shifts with respect to the reduced */
205
/* tridiagonal matrix T. This is equivalent to */
206
/* restarting the iteration with a starting vector */
207
/* that is a linear combination of Ritz vectors */
208
/* associated with the "wanted" Ritz values. */
209
/* ------------------------------------------------------------- */
211
/* IPARAM(2) = LEVEC */
212
/* No longer referenced. See remark 2 below. */
214
/* IPARAM(3) = MXITER */
215
/* On INPUT: maximum number of Arnoldi update iterations allowed. */
216
/* On OUTPUT: actual number of Arnoldi update iterations taken. */
218
/* IPARAM(4) = NB: blocksize to be used in the recurrence. */
219
/* The code currently works only for NB = 1. */
221
/* IPARAM(5) = NCONV: number of "converged" Ritz values. */
222
/* This represents the number of Ritz values that satisfy */
223
/* the convergence criterion. */
225
/* IPARAM(6) = IUPD */
226
/* No longer referenced. Implicit restarting is ALWAYS used. */
228
/* IPARAM(7) = MODE */
229
/* On INPUT determines what type of eigenproblem is being solved. */
230
/* Must be 1,2,3,4,5; See under \Description of dsaupd for the */
231
/* five modes available. */
234
/* When ido = 3 and the user provides shifts through reverse */
235
/* communication (IPARAM(1)=0), dsaupd returns NP, the number */
236
/* of shifts the user is to provide. 0 < NP <=NCV-NEV. See Remark */
239
/* IPARAM(9) = NUMOP, IPARAM(10) = NUMOPB, IPARAM(11) = NUMREO, */
240
/* OUTPUT: NUMOP = total number of OP*x operations, */
241
/* NUMOPB = total number of B*x operations if BMAT='G', */
242
/* NUMREO = total number of steps of re-orthogonalization. */
244
/* IPNTR Integer array of length 11. (OUTPUT) */
245
/* Pointer to mark the starting locations in the WORKD and WORKL */
246
/* arrays for matrices/vectors used by the Lanczos iteration. */
247
/* ------------------------------------------------------------- */
248
/* IPNTR(1): pointer to the current operand vector X in WORKD. */
249
/* IPNTR(2): pointer to the current result vector Y in WORKD. */
250
/* IPNTR(3): pointer to the vector B * X in WORKD when used in */
251
/* the shift-and-invert mode. */
252
/* IPNTR(4): pointer to the next available location in WORKL */
253
/* that is untouched by the program. */
254
/* IPNTR(5): pointer to the NCV by 2 tridiagonal matrix T in WORKL. */
255
/* IPNTR(6): pointer to the NCV RITZ values array in WORKL. */
256
/* IPNTR(7): pointer to the Ritz estimates in array WORKL associated */
257
/* with the Ritz values located in RITZ in WORKL. */
258
/* IPNTR(11): pointer to the NP shifts in WORKL. See Remark 6 below. */
260
/* Note: IPNTR(8:10) is only referenced by dseupd . See Remark 2. */
261
/* IPNTR(8): pointer to the NCV RITZ values of the original system. */
262
/* IPNTR(9): pointer to the NCV corresponding error bounds. */
263
/* IPNTR(10): pointer to the NCV by NCV matrix of eigenvectors */
264
/* of the tridiagonal matrix T. Only referenced by */
265
/* dseupd if RVEC = .TRUE. See Remarks. */
266
/* ------------------------------------------------------------- */
268
/* WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION) */
269
/* Distributed array to be used in the basic Arnoldi iteration */
270
/* for reverse communication. The user should not use WORKD */
271
/* as temporary workspace during the iteration. Upon termination */
272
/* WORKD(1:N) contains B*RESID(1:N). If the Ritz vectors are desired */
273
/* subroutine dseupd uses this output. */
274
/* See Data Distribution Note below. */
276
/* WORKL Double precision work array of length LWORKL. (OUTPUT/WORKSPACE) */
277
/* Private (replicated) array on each PE or array allocated on */
278
/* the front end. See Data Distribution Note below. */
280
/* LWORKL Integer. (INPUT) */
281
/* LWORKL must be at least NCV**2 + 8*NCV . */
283
/* INFO Integer. (INPUT/OUTPUT) */
284
/* If INFO .EQ. 0, a randomly initial residual vector is used. */
285
/* If INFO .NE. 0, RESID contains the initial residual vector, */
286
/* possibly from a previous run. */
287
/* Error flag on output. */
288
/* = 0: Normal exit. */
289
/* = 1: Maximum number of iterations taken. */
290
/* All possible eigenvalues of OP has been found. IPARAM(5) */
291
/* returns the number of wanted converged Ritz values. */
292
/* = 2: No longer an informational error. Deprecated starting */
293
/* with release 2 of ARPACK. */
294
/* = 3: No shifts could be applied during a cycle of the */
295
/* Implicitly restarted Arnoldi iteration. One possibility */
296
/* is to increase the size of NCV relative to NEV. */
297
/* See remark 4 below. */
298
/* = -1: N must be positive. */
299
/* = -2: NEV must be positive. */
300
/* = -3: NCV must be greater than NEV and less than or equal to N. */
301
/* = -4: The maximum number of Arnoldi update iterations allowed */
302
/* must be greater than zero. */
303
/* = -5: WHICH must be one of 'LM', 'SM', 'LA', 'SA' or 'BE'. */
304
/* = -6: BMAT must be one of 'I' or 'G'. */
305
/* = -7: Length of private work array WORKL is not sufficient. */
306
/* = -8: Error return from trid. eigenvalue calculation; */
307
/* Informatinal error from LAPACK routine dsteqr . */
308
/* = -9: Starting vector is zero. */
309
/* = -10: IPARAM(7) must be 1,2,3,4,5. */
310
/* = -11: IPARAM(7) = 1 and BMAT = 'G' are incompatable. */
311
/* = -12: IPARAM(1) must be equal to 0 or 1. */
312
/* = -13: NEV and WHICH = 'BE' are incompatable. */
313
/* = -9999: Could not build an Arnoldi factorization. */
314
/* IPARAM(5) returns the size of the current Arnoldi */
315
/* factorization. The user is advised to check that */
316
/* enough workspace and array storage has been allocated. */
320
/* 1. The converged Ritz values are always returned in ascending */
321
/* algebraic order. The computed Ritz values are approximate */
322
/* eigenvalues of OP. The selection of WHICH should be made */
323
/* with this in mind when Mode = 3,4,5. After convergence, */
324
/* approximate eigenvalues of the original problem may be obtained */
325
/* with the ARPACK subroutine dseupd . */
327
/* 2. If the Ritz vectors corresponding to the converged Ritz values */
328
/* are needed, the user must call dseupd immediately following completion */
329
/* of dsaupd . This is new starting with version 2.1 of ARPACK. */
331
/* 3. If M can be factored into a Cholesky factorization M = LL` */
332
/* then Mode = 2 should not be selected. Instead one should use */
333
/* Mode = 1 with OP = inv(L)*A*inv(L`). Appropriate triangular */
334
/* linear systems should be solved with L and L` rather */
335
/* than computing inverses. After convergence, an approximate */
336
/* eigenvector z of the original problem is recovered by solving */
337
/* L`z = x where x is a Ritz vector of OP. */
339
/* 4. At present there is no a-priori analysis to guide the selection */
340
/* of NCV relative to NEV. The only formal requrement is that NCV > NEV. */
341
/* However, it is recommended that NCV .ge. 2*NEV. If many problems of */
342
/* the same type are to be solved, one should experiment with increasing */
343
/* NCV while keeping NEV fixed for a given test problem. This will */
344
/* usually decrease the required number of OP*x operations but it */
345
/* also increases the work and storage required to maintain the orthogonal */
346
/* basis vectors. The optimal "cross-over" with respect to CPU time */
347
/* is problem dependent and must be determined empirically. */
349
/* 5. If IPARAM(7) = 2 then in the Reverse commuication interface the user */
350
/* must do the following. When IDO = 1, Y = OP * X is to be computed. */
351
/* When IPARAM(7) = 2 OP = inv(B)*A. After computing A*X the user */
352
/* must overwrite X with A*X. Y is then the solution to the linear set */
353
/* of equations B*Y = A*X. */
355
/* 6. When IPARAM(1) = 0, and IDO = 3, the user needs to provide the */
356
/* NP = IPARAM(8) shifts in locations: */
357
/* 1 WORKL(IPNTR(11)) */
358
/* 2 WORKL(IPNTR(11)+1) */
362
/* NP WORKL(IPNTR(11)+NP-1). */
364
/* The eigenvalues of the current tridiagonal matrix are located in */
365
/* WORKL(IPNTR(6)) through WORKL(IPNTR(6)+NCV-1). They are in the */
366
/* order defined by WHICH. The associated Ritz estimates are located in */
367
/* WORKL(IPNTR(8)), WORKL(IPNTR(8)+1), ... , WORKL(IPNTR(8)+NCV-1). */
369
/* ----------------------------------------------------------------------- */
371
/* \Data Distribution Note: */
373
/* Fortran-D syntax: */
374
/* ================ */
375
/* REAL RESID(N), V(LDV,NCV), WORKD(3*N), WORKL(LWORKL) */
376
/* DECOMPOSE D1(N), D2(N,NCV) */
377
/* ALIGN RESID(I) with D1(I) */
378
/* ALIGN V(I,J) with D2(I,J) */
379
/* ALIGN WORKD(I) with D1(I) range (1:N) */
380
/* ALIGN WORKD(I) with D1(I-N) range (N+1:2*N) */
381
/* ALIGN WORKD(I) with D1(I-2*N) range (2*N+1:3*N) */
382
/* DISTRIBUTE D1(BLOCK), D2(BLOCK,:) */
383
/* REPLICATED WORKL(LWORKL) */
385
/* Cray MPP syntax: */
386
/* =============== */
387
/* REAL RESID(N), V(LDV,NCV), WORKD(N,3), WORKL(LWORKL) */
388
/* SHARED RESID(BLOCK), V(BLOCK,:), WORKD(BLOCK,:) */
389
/* REPLICATED WORKL(LWORKL) */
395
/* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
396
/* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
398
/* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
399
/* Restarted Arnoldi Iteration", Rice University Technical Report */
400
/* TR95-13, Department of Computational and Applied Mathematics. */
401
/* 3. B.N. Parlett, "The Symmetric Eigenvalue Problem". Prentice-Hall, */
403
/* 4. B.N. Parlett, B. Nour-Omid, "Towards a Black Box Lanczos Program", */
404
/* Computer Physics Communications, 53 (1989), pp 169-179. */
405
/* 5. B. Nour-Omid, B.N. Parlett, T. Ericson, P.S. Jensen, "How to */
406
/* Implement the Spectral Transformation", Math. Comp., 48 (1987), */
408
/* 6. R.G. Grimes, J.G. Lewis and H.D. Simon, "A Shifted Block Lanczos */
409
/* Algorithm for Solving Sparse Symmetric Generalized Eigenproblems", */
410
/* SIAM J. Matr. Anal. Apps., January (1993). */
411
/* 7. L. Reichel, W.B. Gragg, "Algorithm 686: FORTRAN Subroutines */
412
/* for Updating the QR decomposition", ACM TOMS, December 1990, */
413
/* Volume 16 Number 4, pp 369-377. */
414
/* 8. R.B. Lehoucq, D.C. Sorensen, "Implementation of Some Spectral */
415
/* Transformations in a k-Step Arnoldi Method". In Preparation. */
417
/* \Routines called: */
418
/* dsaup2 ARPACK routine that implements the Implicitly Restarted */
419
/* Arnoldi Iteration. */
420
/* dstats ARPACK routine that initialize timing and other statistics */
422
/* ivout ARPACK utility routine that prints integers. */
423
/* second ARPACK utility routine for timing. */
424
/* dvout ARPACK utility routine that prints vectors. */
425
/* dlamch LAPACK routine that determines machine constants. */
428
/* Danny Sorensen Phuong Vu */
429
/* Richard Lehoucq CRPC / Rice University */
430
/* Dept. of Computational & Houston, Texas */
431
/* Applied Mathematics */
432
/* Rice University */
435
/* \Revision history: */
436
/* 12/15/93: Version ' 2.4' */
438
/* \SCCS Information: @(#) */
439
/* FILE: saupd.F SID: 2.8 DATE OF SID: 04/10/01 RELEASE: 2 */
446
/* ----------------------------------------------------------------------- */
448
/* Subroutine */ int igraphdsaupd_(integer *ido, char *bmat, integer *n, char
449
*which, integer *nev, doublereal *tol, doublereal *resid, integer *
450
ncv, doublereal *v, integer *ldv, integer *iparam, integer *ipntr,
451
doublereal *workd, doublereal *workl, integer *lworkl, integer *info)
454
/* static char fmt_1000[] = "(//,5x,\002===================================" */
455
/* "=======\002,/5x,\002= Symmetric implicit Arnoldi update code " */
456
/* "=\002,/5x,\002= Version Number:\002,\002 2.4\002,19x,\002 =\002," */
457
/* "/5x,\002= Version Date: \002,\002 07/31/96\002,14x,\002 =\002,/" */
458
/* "5x,\002==========================================\002,/5x,\002= " */
459
/* "Summary of timing statistics =\002,/5x,\002===========" */
460
/* "===============================\002,//)"; */
461
/* static char fmt_1100[] = "(5x,\002Total number update iterations " */
462
/* " = \002,i5,/5x,\002Total number of OP*x operations " */
463
/* " = \002,i5,/5x,\002Total number of B*x operations = " */
464
/* "\002,i5,/5x,\002Total number of reorthogonalization steps = " */
465
/* "\002,i5,/5x,\002Total number of iterative refinement steps = " */
466
/* "\002,i5,/5x,\002Total number of restart steps = " */
467
/* "\002,i5,/5x,\002Total time in user OP*x operation = " */
468
/* "\002,f12.6,/5x,\002Total time in user B*x operation =" */
469
/* " \002,f12.6,/5x,\002Total time in Arnoldi update routine = " */
470
/* "\002,f12.6,/5x,\002Total time in saup2 routine =" */
471
/* " \002,f12.6,/5x,\002Total time in basic Arnoldi iteration loop = " */
472
/* "\002,f12.6,/5x,\002Total time in reorthogonalization phase =" */
473
/* " \002,f12.6,/5x,\002Total time in (re)start vector generation = " */
474
/* "\002,f12.6,/5x,\002Total time in trid eigenvalue subproblem =" */
475
/* " \002,f12.6,/5x,\002Total time in getting the shifts = " */
476
/* "\002,f12.6,/5x,\002Total time in applying the shifts =" */
477
/* " \002,f12.6,/5x,\002Total time in convergence testing = " */
480
/* System generated locals */
481
integer v_dim1, v_offset, i__1, i__2;
483
/* Builtin functions */
484
integer igraphs_cmp(char *, char *, ftnlen, ftnlen)/* , s_wsfe(cilist *), e_wsfe( */
485
/* void), do_fio(integer *, char *, ftnlen) */;
487
/* Local variables */
490
static integer nb, ih, iq, np, iw, ldh, ldq, nev0;
491
extern /* Subroutine */ int igraphdvout_(integer *, integer *, doublereal
492
*, integer *, char *), igraphivout_(integer *, integer *,
493
integer *, integer *, char *);
495
extern /* Subroutine */ int igraphdsaup2_(integer *, char *, integer *,
496
char *, integer *, integer *, doublereal *, doublereal *, integer
497
*, integer *, integer *, integer *, doublereal *, integer *,
498
doublereal *, integer *, doublereal *, doublereal *, doublereal *,
499
integer *, doublereal *, integer *, doublereal *, integer *);
500
static integer ierr, iupd, next;
501
extern doublereal igraphdlamch_(char *);
503
extern /* Subroutine */ int igraphsecond_(real *), igraphdstats_(void);
504
static integer bounds, ishift, msglvl, mxiter;
506
/* Fortran I/O blocks */
507
/* static cilist io___21 = { 0, 6, 0, fmt_1000, 0 }; */
508
/* static cilist io___22 = { 0, 6, 0, fmt_1100, 0 }; */
512
/* %----------------------------------------------------% */
513
/* | Include files for debugging and timing information | */
514
/* %----------------------------------------------------% */
517
/* \SCCS Information: @(#) */
518
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
520
/* %---------------------------------% */
521
/* | See debug.doc for documentation | */
522
/* %---------------------------------% */
524
/* %------------------% */
525
/* | Scalar Arguments | */
526
/* %------------------% */
528
/* %--------------------------------% */
529
/* | See stat.doc for documentation | */
530
/* %--------------------------------% */
532
/* \SCCS Information: @(#) */
533
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
537
/* %-----------------% */
538
/* | Array Arguments | */
539
/* %-----------------% */
547
/* %---------------% */
548
/* | Local Scalars | */
549
/* %---------------% */
552
/* %----------------------% */
553
/* | External Subroutines | */
554
/* %----------------------% */
557
/* %--------------------% */
558
/* | External Functions | */
559
/* %--------------------% */
562
/* %-----------------------% */
563
/* | Executable Statements | */
564
/* %-----------------------% */
566
/* Parameter adjustments */
570
v_offset = 1 + v_dim1;
579
/* %-------------------------------% */
580
/* | Initialize timing statistics | */
581
/* | & message level for debugging | */
582
/* %-------------------------------% */
586
msglvl = debug_1.msaupd;
594
/* %--------------------------------------------% */
595
/* | Revision 2 performs only implicit restart. | */
596
/* %--------------------------------------------% */
601
/* %----------------% */
602
/* | Error checking | */
603
/* %----------------% */
607
} else if (*nev <= 0) {
609
} else if (*ncv <= *nev || *ncv > *n) {
613
/* %----------------------------------------------% */
614
/* | NP is the number of additional steps to | */
615
/* | extend the length NEV Lanczos factorization. | */
616
/* %----------------------------------------------% */
623
if (igraphs_cmp(which, "LM", (ftnlen)2, (ftnlen)2) != 0 && igraphs_cmp(which,
624
"SM", (ftnlen)2, (ftnlen)2) != 0 && igraphs_cmp(which, "LA", (
625
ftnlen)2, (ftnlen)2) != 0 && igraphs_cmp(which, "SA", (ftnlen)2, (
626
ftnlen)2) != 0 && igraphs_cmp(which, "BE", (ftnlen)2, (ftnlen)2) !=
630
if (*(unsigned char *)bmat != 'I' && *(unsigned char *)bmat != 'G') {
634
/* Computing 2nd power */
636
if (*lworkl < i__1 * i__1 + (*ncv << 3)) {
639
if (mode < 1 || mode > 5) {
641
} else if (mode == 1 && *(unsigned char *)bmat == 'G') {
643
} else if (ishift < 0 || ishift > 1) {
645
} else if (*nev == 1 && igraphs_cmp(which, "BE", (ftnlen)2, (ftnlen)2) == 0)
660
/* %------------------------% */
661
/* | Set default parameters | */
662
/* %------------------------% */
668
*tol = igraphdlamch_("EpsMach");
671
/* %----------------------------------------------% */
672
/* | NP is the number of additional steps to | */
673
/* | extend the length NEV Lanczos factorization. | */
674
/* | NEV0 is the local variable designating the | */
675
/* | size of the invariant subspace desired. | */
676
/* %----------------------------------------------% */
681
/* %-----------------------------% */
682
/* | Zero out internal workspace | */
683
/* %-----------------------------% */
685
/* Computing 2nd power */
687
i__1 = i__2 * i__2 + (*ncv << 3);
688
for (j = 1; j <= i__1; ++j) {
693
/* %-------------------------------------------------------% */
694
/* | Pointer into WORKL for address of H, RITZ, BOUNDS, Q | */
695
/* | etc... and the remaining workspace. | */
696
/* | Also update pointer to be used on output. | */
697
/* | Memory is laid out as follows: | */
698
/* | workl(1:2*ncv) := generated tridiagonal matrix | */
699
/* | workl(2*ncv+1:2*ncv+ncv) := ritz values | */
700
/* | workl(3*ncv+1:3*ncv+ncv) := computed error bounds | */
701
/* | workl(4*ncv+1:4*ncv+ncv*ncv) := rotation matrix Q | */
702
/* | workl(4*ncv+ncv*ncv+1:7*ncv+ncv*ncv) := workspace | */
703
/* %-------------------------------------------------------% */
708
ritz = ih + (ldh << 1);
709
bounds = ritz + *ncv;
711
/* Computing 2nd power */
713
iw = iq + i__1 * i__1;
714
next = iw + *ncv * 3;
723
/* %-------------------------------------------------------% */
724
/* | Carry out the Implicitly restarted Lanczos Iteration. | */
725
/* %-------------------------------------------------------% */
727
igraphdsaup2_(ido, bmat, n, which, &nev0, &np, tol, &resid[1], &mode, &
728
iupd, &ishift, &mxiter, &v[v_offset], ldv, &workl[ih], &ldh, &
729
workl[ritz], &workl[bounds], &workl[iq], &ldq, &workl[iw], &ipntr[
730
1], &workd[1], info);
732
/* %--------------------------------------------------% */
733
/* | ido .ne. 99 implies use of reverse communication | */
734
/* | to compute operations involving OP or shifts. | */
735
/* %--------------------------------------------------% */
746
iparam[9] = timing_1.nopx;
747
iparam[10] = timing_1.nbx;
748
iparam[11] = timing_1.nrorth;
750
/* %------------------------------------% */
751
/* | Exit if there was an informational | */
752
/* | error within dsaup2 . | */
753
/* %------------------------------------% */
763
igraphivout_(&debug_1.logfil, &c__1, &mxiter, &debug_1.ndigit, "_sau"
764
"pd: number of update iterations taken");
765
igraphivout_(&debug_1.logfil, &c__1, &np, &debug_1.ndigit, "_saupd: "
766
"number of \"converged\" Ritz values");
767
igraphdvout_(&debug_1.logfil, &np, &workl[ritz], &debug_1.ndigit,
768
"_saupd: final Ritz values");
769
igraphdvout_(&debug_1.logfil, &np, &workl[bounds], &debug_1.ndigit,
770
"_saupd: corresponding error bounds");
774
timing_1.tsaupd = t1 - t0;
778
/* %--------------------------------------------------------% */
779
/* | Version Number & Version Date are defined in version.h | */
780
/* %--------------------------------------------------------% */
782
/* s_wsfe(&io___21); */
784
/* s_wsfe(&io___22); */
785
/* do_fio(&c__1, (char *)&mxiter, (ftnlen)sizeof(integer)); */
786
/* do_fio(&c__1, (char *)&timing_1.nopx, (ftnlen)sizeof(integer)); */
787
/* do_fio(&c__1, (char *)&timing_1.nbx, (ftnlen)sizeof(integer)); */
788
/* do_fio(&c__1, (char *)&timing_1.nrorth, (ftnlen)sizeof(integer)); */
789
/* do_fio(&c__1, (char *)&timing_1.nitref, (ftnlen)sizeof(integer)); */
790
/* do_fio(&c__1, (char *)&timing_1.nrstrt, (ftnlen)sizeof(integer)); */
791
/* do_fio(&c__1, (char *)&timing_1.tmvopx, (ftnlen)sizeof(real)); */
792
/* do_fio(&c__1, (char *)&timing_1.tmvbx, (ftnlen)sizeof(real)); */
793
/* do_fio(&c__1, (char *)&timing_1.tsaupd, (ftnlen)sizeof(real)); */
794
/* do_fio(&c__1, (char *)&timing_1.tsaup2, (ftnlen)sizeof(real)); */
795
/* do_fio(&c__1, (char *)&timing_1.tsaitr, (ftnlen)sizeof(real)); */
796
/* do_fio(&c__1, (char *)&timing_1.titref, (ftnlen)sizeof(real)); */
797
/* do_fio(&c__1, (char *)&timing_1.tgetv0, (ftnlen)sizeof(real)); */
798
/* do_fio(&c__1, (char *)&timing_1.tseigt, (ftnlen)sizeof(real)); */
799
/* do_fio(&c__1, (char *)&timing_1.tsgets, (ftnlen)sizeof(real)); */
800
/* do_fio(&c__1, (char *)&timing_1.tsapps, (ftnlen)sizeof(real)); */
801
/* do_fio(&c__1, (char *)&timing_1.tsconv, (ftnlen)sizeof(real)); */
809
/* %---------------% */
810
/* | End of dsaupd | */
811
/* %---------------% */
813
} /* igraphdsaupd_ */
added
removed
src/arpack/dsaupd.c
