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;
40
static doublereal c_b24 = 1.;
41
static doublereal c_b26 = 0.;
42
static doublereal c_b29 = -1.;
44
/* ----------------------------------------------------------------------- */
50
/* Generate a random initial residual vector for the Arnoldi process. */
51
/* Force the residual vector to be in the range of the operator OP. */
55
/* ( IDO, BMAT, ITRY, INITV, N, J, V, LDV, RESID, RNORM, */
56
/* IPNTR, WORKD, IERR ) */
59
/* IDO Integer. (INPUT/OUTPUT) */
60
/* Reverse communication flag. IDO must be zero on the first */
62
/* ------------------------------------------------------------- */
63
/* IDO = 0: first call to the reverse communication interface */
64
/* IDO = -1: compute Y = OP * X where */
65
/* IPNTR(1) is the pointer into WORKD for X, */
66
/* IPNTR(2) is the pointer into WORKD for Y. */
67
/* This is for the initialization phase to force the */
68
/* starting vector into the range of OP. */
69
/* IDO = 2: compute Y = B * X where */
70
/* IPNTR(1) is the pointer into WORKD for X, */
71
/* IPNTR(2) is the pointer into WORKD for Y. */
73
/* ------------------------------------------------------------- */
75
/* BMAT Character*1. (INPUT) */
76
/* BMAT specifies the type of the matrix B in the (generalized) */
77
/* eigenvalue problem A*x = lambda*B*x. */
78
/* B = 'I' -> standard eigenvalue problem A*x = lambda*x */
79
/* B = 'G' -> generalized eigenvalue problem A*x = lambda*B*x */
81
/* ITRY Integer. (INPUT) */
82
/* ITRY counts the number of times that dgetv0 is called. */
83
/* It should be set to 1 on the initial call to dgetv0. */
85
/* INITV Logical variable. (INPUT) */
86
/* .TRUE. => the initial residual vector is given in RESID. */
87
/* .FALSE. => generate a random initial residual vector. */
89
/* N Integer. (INPUT) */
90
/* Dimension of the problem. */
92
/* J Integer. (INPUT) */
93
/* Index of the residual vector to be generated, with respect to */
94
/* the Arnoldi process. J > 1 in case of a "restart". */
96
/* V Double precision N by J array. (INPUT) */
97
/* The first J-1 columns of V contain the current Arnoldi basis */
98
/* if this is a "restart". */
100
/* LDV Integer. (INPUT) */
101
/* Leading dimension of V exactly as declared in the calling */
104
/* RESID Double precision array of length N. (INPUT/OUTPUT) */
105
/* Initial residual vector to be generated. If RESID is */
106
/* provided, force RESID into the range of the operator OP. */
108
/* RNORM Double precision scalar. (OUTPUT) */
109
/* B-norm of the generated residual. */
111
/* IPNTR Integer array of length 3. (OUTPUT) */
113
/* WORKD Double precision work array of length 2*N. (REVERSE COMMUNICATION). */
114
/* On exit, WORK(1:N) = B*RESID to be used in SSAITR. */
116
/* IERR Integer. (OUTPUT) */
117
/* = 0: Normal exit. */
118
/* = -1: Cannot generate a nontrivial restarted residual vector */
119
/* in the range of the operator OP. */
123
/* ----------------------------------------------------------------------- */
127
/* \Local variables: */
131
/* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
132
/* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
134
/* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
135
/* Restarted Arnoldi Iteration", Rice University Technical Report */
136
/* TR95-13, Department of Computational and Applied Mathematics. */
138
/* \Routines called: */
139
/* second ARPACK utility routine for timing. */
140
/* dvout ARPACK utility routine for vector output. */
141
/* dlarnv LAPACK routine for generating a random vector. */
142
/* dgemv Level 2 BLAS routine for matrix vector multiplication. */
143
/* dcopy Level 1 BLAS that copies one vector to another. */
144
/* ddot Level 1 BLAS that computes the scalar product of two vectors. */
145
/* dnrm2 Level 1 BLAS that computes the norm of a vector. */
148
/* Danny Sorensen Phuong Vu */
149
/* Richard Lehoucq CRPC / Rice University */
150
/* Dept. of Computational & Houston, Texas */
151
/* Applied Mathematics */
152
/* Rice University */
155
/* \SCCS Information: @(#) */
156
/* FILE: getv0.F SID: 2.7 DATE OF SID: 04/07/99 RELEASE: 2 */
160
/* ----------------------------------------------------------------------- */
162
/* Subroutine */ int igraphdgetv0_(integer *ido, char *bmat, integer *itry,
163
logical *initv, integer *n, integer *j, doublereal *v, integer *ldv,
164
doublereal *resid, doublereal *rnorm, integer *ipntr, doublereal *
165
workd, integer *ierr)
167
/* Initialized data */
169
static logical inits = TRUE_;
171
/* System generated locals */
172
integer v_dim1, v_offset, i__1;
174
/* Builtin functions */
175
double sqrt(doublereal);
177
/* Local variables */
178
static real t0, t1, t2, t3;
180
extern doublereal igraphddot_(integer *, doublereal *, integer *,
181
doublereal *, integer *), igraphdnrm2_(integer *, doublereal *,
183
extern /* Subroutine */ int igraphdgemv_(char *, integer *, integer *,
184
doublereal *, doublereal *, integer *, doublereal *, integer *,
185
doublereal *, doublereal *, integer *), igraphdcopy_(
186
integer *, doublereal *, integer *, doublereal *, integer *),
187
igraphdvout_(integer *, integer *, doublereal *, integer *, char *
191
extern /* Subroutine */ int igraphsecond_(real *), igraphdlarnv_(integer *
192
, integer *, integer *, doublereal *);
193
static integer iseed[4], idist;
194
static logical first;
195
static doublereal rnorm0;
196
static integer msglvl;
199
/* %----------------------------------------------------% */
200
/* | Include files for debugging and timing information | */
201
/* %----------------------------------------------------% */
204
/* \SCCS Information: @(#) */
205
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
207
/* %---------------------------------% */
208
/* | See debug.doc for documentation | */
209
/* %---------------------------------% */
211
/* %------------------% */
212
/* | Scalar Arguments | */
213
/* %------------------% */
215
/* %--------------------------------% */
216
/* | See stat.doc for documentation | */
217
/* %--------------------------------% */
219
/* \SCCS Information: @(#) */
220
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
224
/* %-----------------% */
225
/* | Array Arguments | */
226
/* %-----------------% */
234
/* %------------------------% */
235
/* | Local Scalars & Arrays | */
236
/* %------------------------% */
239
/* %----------------------% */
240
/* | External Subroutines | */
241
/* %----------------------% */
244
/* %--------------------% */
245
/* | External Functions | */
246
/* %--------------------% */
249
/* %---------------------% */
250
/* | Intrinsic Functions | */
251
/* %---------------------% */
254
/* %-----------------% */
255
/* | Data Statements | */
256
/* %-----------------% */
258
/* Parameter adjustments */
262
v_offset = 1 + v_dim1;
268
/* %-----------------------% */
269
/* | Executable Statements | */
270
/* %-----------------------% */
273
/* %-----------------------------------% */
274
/* | Initialize the seed of the LAPACK | */
275
/* | random number generator | */
276
/* %-----------------------------------% */
288
/* %-------------------------------% */
289
/* | Initialize timing statistics | */
290
/* | & message level for debugging | */
291
/* %-------------------------------% */
294
msglvl = debug_1.mgetv0;
301
/* %-----------------------------------------------------% */
302
/* | Possibly generate a random starting vector in RESID | */
303
/* | Use a LAPACK random number generator used by the | */
304
/* | matrix generation routines. | */
305
/* | idist = 1: uniform (0,1) distribution; | */
306
/* | idist = 2: uniform (-1,1) distribution; | */
307
/* | idist = 3: normal (0,1) distribution; | */
308
/* %-----------------------------------------------------% */
312
igraphdlarnv_(&idist, iseed, n, &resid[1]);
315
/* %----------------------------------------------------------% */
316
/* | Force the starting vector into the range of OP to handle | */
317
/* | the generalized problem when B is possibly (singular). | */
318
/* %----------------------------------------------------------% */
321
if (*(unsigned char *)bmat == 'G') {
325
igraphdcopy_(n, &resid[1], &c__1, &workd[1], &c__1);
331
/* %-----------------------------------------% */
332
/* | Back from computing OP*(initial-vector) | */
333
/* %-----------------------------------------% */
339
/* %-----------------------------------------------% */
340
/* | Back from computing B*(orthogonalized-vector) | */
341
/* %-----------------------------------------------% */
347
if (*(unsigned char *)bmat == 'G') {
349
timing_1.tmvopx += t3 - t2;
352
/* %------------------------------------------------------% */
353
/* | Starting vector is now in the range of OP; r = OP*r; | */
354
/* | Compute B-norm of starting vector. | */
355
/* %------------------------------------------------------% */
359
if (*(unsigned char *)bmat == 'G') {
361
igraphdcopy_(n, &workd[*n + 1], &c__1, &resid[1], &c__1);
366
} else if (*(unsigned char *)bmat == 'I') {
367
igraphdcopy_(n, &resid[1], &c__1, &workd[1], &c__1);
372
if (*(unsigned char *)bmat == 'G') {
374
timing_1.tmvbx += t3 - t2;
378
if (*(unsigned char *)bmat == 'G') {
379
rnorm0 = igraphddot_(n, &resid[1], &c__1, &workd[1], &c__1);
380
rnorm0 = sqrt((abs(rnorm0)));
381
} else if (*(unsigned char *)bmat == 'I') {
382
rnorm0 = igraphdnrm2_(n, &resid[1], &c__1);
386
/* %---------------------------------------------% */
387
/* | Exit if this is the very first Arnoldi step | */
388
/* %---------------------------------------------% */
394
/* %---------------------------------------------------------------- */
395
/* | Otherwise need to B-orthogonalize the starting vector against | */
396
/* | the current Arnoldi basis using Gram-Schmidt with iter. ref. | */
397
/* | This is the case where an invariant subspace is encountered | */
398
/* | in the middle of the Arnoldi factorization. | */
400
/* | s = V^{T}*B*r; r = r - V*s; | */
402
/* | Stopping criteria used for iter. ref. is discussed in | */
403
/* | Parlett's book, page 107 and in Gragg & Reichel TOMS paper. | */
404
/* %---------------------------------------------------------------% */
410
igraphdgemv_("T", n, &i__1, &c_b24, &v[v_offset], ldv, &workd[1], &c__1, &
411
c_b26, &workd[*n + 1], &c__1);
413
igraphdgemv_("N", n, &i__1, &c_b29, &v[v_offset], ldv, &workd[*n + 1], &
414
c__1, &c_b24, &resid[1], &c__1);
416
/* %----------------------------------------------------------% */
417
/* | Compute the B-norm of the orthogonalized starting vector | */
418
/* %----------------------------------------------------------% */
421
if (*(unsigned char *)bmat == 'G') {
423
igraphdcopy_(n, &resid[1], &c__1, &workd[*n + 1], &c__1);
428
} else if (*(unsigned char *)bmat == 'I') {
429
igraphdcopy_(n, &resid[1], &c__1, &workd[1], &c__1);
434
if (*(unsigned char *)bmat == 'G') {
436
timing_1.tmvbx += t3 - t2;
439
if (*(unsigned char *)bmat == 'G') {
440
*rnorm = igraphddot_(n, &resid[1], &c__1, &workd[1], &c__1);
441
*rnorm = sqrt((abs(*rnorm)));
442
} else if (*(unsigned char *)bmat == 'I') {
443
*rnorm = igraphdnrm2_(n, &resid[1], &c__1);
446
/* %--------------------------------------% */
447
/* | Check for further orthogonalization. | */
448
/* %--------------------------------------% */
451
igraphdvout_(&debug_1.logfil, &c__1, &rnorm0, &debug_1.ndigit, "_get"
452
"v0: re-orthonalization ; rnorm0 is");
453
igraphdvout_(&debug_1.logfil, &c__1, rnorm, &debug_1.ndigit, "_getv0"
454
": re-orthonalization ; rnorm is");
457
if (*rnorm > rnorm0 * .717f) {
464
/* %-----------------------------------% */
465
/* | Perform iterative refinement step | */
466
/* %-----------------------------------% */
472
/* %------------------------------------% */
473
/* | Iterative refinement step "failed" | */
474
/* %------------------------------------% */
477
for (jj = 1; jj <= i__1; ++jj) {
488
igraphdvout_(&debug_1.logfil, &c__1, rnorm, &debug_1.ndigit, "_getv0"
489
": B-norm of initial / restarted starting vector");
492
igraphdvout_(&debug_1.logfil, n, &resid[1], &debug_1.ndigit, "_getv0"
493
": initial / restarted starting vector");
498
timing_1.tgetv0 += t1 - t0;
503
/* %---------------% */
504
/* | End of dgetv0 | */
505
/* %---------------% */
507
} /* igraphdgetv0_ */
added
removed
src/arpack/dgetv0.c
