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 logical c_false = FALSE_;
41
static doublereal c_b24 = 1.;
42
static doublereal c_b49 = 0.;
43
static doublereal c_b57 = -1.;
44
static integer c__2 = 2;
46
/* ----------------------------------------------------------------------- */
52
/* Reverse communication interface for applying NP additional steps to */
53
/* a K step symmetric Arnoldi factorization. */
55
/* Input: OP*V_{k} - V_{k}*H = r_{k}*e_{k}^T */
57
/* with (V_{k}^T)*B*V_{k} = I, (V_{k}^T)*B*r_{k} = 0. */
59
/* Output: OP*V_{k+p} - V_{k+p}*H = r_{k+p}*e_{k+p}^T */
61
/* with (V_{k+p}^T)*B*V_{k+p} = I, (V_{k+p}^T)*B*r_{k+p} = 0. */
63
/* where OP and B are as in dsaupd. The B-norm of r_{k+p} is also */
64
/* computed and returned. */
68
/* ( IDO, BMAT, N, K, NP, MODE, RESID, RNORM, V, LDV, H, LDH, */
69
/* IPNTR, WORKD, INFO ) */
72
/* IDO Integer. (INPUT/OUTPUT) */
73
/* Reverse communication flag. */
74
/* ------------------------------------------------------------- */
75
/* IDO = 0: first call to the reverse communication interface */
76
/* IDO = -1: compute Y = OP * X where */
77
/* IPNTR(1) is the pointer into WORK for X, */
78
/* IPNTR(2) is the pointer into WORK for Y. */
79
/* This is for the restart phase to force the new */
80
/* starting vector into the range of OP. */
81
/* IDO = 1: compute Y = OP * X where */
82
/* IPNTR(1) is the pointer into WORK for X, */
83
/* IPNTR(2) is the pointer into WORK for Y, */
84
/* IPNTR(3) is the pointer into WORK for B * X. */
85
/* IDO = 2: compute Y = B * X where */
86
/* IPNTR(1) is the pointer into WORK for X, */
87
/* IPNTR(2) is the pointer into WORK for Y. */
89
/* ------------------------------------------------------------- */
90
/* When the routine is used in the "shift-and-invert" mode, the */
91
/* vector B * Q is already available and does not need to be */
92
/* recomputed in forming OP * Q. */
94
/* BMAT Character*1. (INPUT) */
95
/* BMAT specifies the type of matrix B that defines the */
96
/* semi-inner product for the operator OP. See dsaupd. */
97
/* B = 'I' -> standard eigenvalue problem A*x = lambda*x */
98
/* B = 'G' -> generalized eigenvalue problem A*x = lambda*M*x */
100
/* N Integer. (INPUT) */
101
/* Dimension of the eigenproblem. */
103
/* K Integer. (INPUT) */
104
/* Current order of H and the number of columns of V. */
106
/* NP Integer. (INPUT) */
107
/* Number of additional Arnoldi steps to take. */
109
/* MODE Integer. (INPUT) */
110
/* Signifies which form for "OP". If MODE=2 then */
111
/* a reduction in the number of B matrix vector multiplies */
112
/* is possible since the B-norm of OP*x is equivalent to */
113
/* the inv(B)-norm of A*x. */
115
/* RESID Double precision array of length N. (INPUT/OUTPUT) */
116
/* On INPUT: RESID contains the residual vector r_{k}. */
117
/* On OUTPUT: RESID contains the residual vector r_{k+p}. */
119
/* RNORM Double precision scalar. (INPUT/OUTPUT) */
120
/* On INPUT the B-norm of r_{k}. */
121
/* On OUTPUT the B-norm of the updated residual r_{k+p}. */
123
/* V Double precision N by K+NP array. (INPUT/OUTPUT) */
124
/* On INPUT: V contains the Arnoldi vectors in the first K */
126
/* On OUTPUT: V contains the new NP Arnoldi vectors in the next */
127
/* NP columns. The first K columns are unchanged. */
129
/* LDV Integer. (INPUT) */
130
/* Leading dimension of V exactly as declared in the calling */
133
/* H Double precision (K+NP) by 2 array. (INPUT/OUTPUT) */
134
/* H is used to store the generated symmetric tridiagonal matrix */
135
/* with the subdiagonal in the first column starting at H(2,1) */
136
/* and the main diagonal in the second column. */
138
/* LDH Integer. (INPUT) */
139
/* Leading dimension of H exactly as declared in the calling */
142
/* IPNTR Integer array of length 3. (OUTPUT) */
143
/* Pointer to mark the starting locations in the WORK for */
144
/* vectors used by the Arnoldi iteration. */
145
/* ------------------------------------------------------------- */
146
/* IPNTR(1): pointer to the current operand vector X. */
147
/* IPNTR(2): pointer to the current result vector Y. */
148
/* IPNTR(3): pointer to the vector B * X when used in the */
149
/* shift-and-invert mode. X is the current operand. */
150
/* ------------------------------------------------------------- */
152
/* WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION) */
153
/* Distributed array to be used in the basic Arnoldi iteration */
154
/* for reverse communication. The calling program should not */
155
/* use WORKD as temporary workspace during the iteration !!!!!! */
156
/* On INPUT, WORKD(1:N) = B*RESID where RESID is associated */
157
/* with the K step Arnoldi factorization. Used to save some */
158
/* computation at the first step. */
159
/* On OUTPUT, WORKD(1:N) = B*RESID where RESID is associated */
160
/* with the K+NP step Arnoldi factorization. */
162
/* INFO Integer. (OUTPUT) */
163
/* = 0: Normal exit. */
164
/* > 0: Size of an invariant subspace of OP is found that is */
165
/* less than K + NP. */
169
/* ----------------------------------------------------------------------- */
173
/* \Local variables: */
176
/* \Routines called: */
177
/* dgetv0 ARPACK routine to generate the initial vector. */
178
/* ivout ARPACK utility routine that prints integers. */
179
/* dmout ARPACK utility routine that prints matrices. */
180
/* dvout ARPACK utility routine that prints vectors. */
181
/* dlamch LAPACK routine that determines machine constants. */
182
/* dlascl LAPACK routine for careful scaling of a matrix. */
183
/* dgemv Level 2 BLAS routine for matrix vector multiplication. */
184
/* daxpy Level 1 BLAS that computes a vector triad. */
185
/* dscal Level 1 BLAS that scales a vector. */
186
/* dcopy Level 1 BLAS that copies one vector to another . */
187
/* ddot Level 1 BLAS that computes the scalar product of two vectors. */
188
/* dnrm2 Level 1 BLAS that computes the norm of a vector. */
191
/* Danny Sorensen Phuong Vu */
192
/* Richard Lehoucq CRPC / Rice University */
193
/* Dept. of Computational & Houston, Texas */
194
/* Applied Mathematics */
195
/* Rice University */
198
/* \Revision history: */
199
/* xx/xx/93: Version ' 2.4' */
201
/* \SCCS Information: @(#) */
202
/* FILE: saitr.F SID: 2.6 DATE OF SID: 8/28/96 RELEASE: 2 */
205
/* The algorithm implemented is: */
207
/* restart = .false. */
208
/* Given V_{k} = [v_{1}, ..., v_{k}], r_{k}; */
209
/* r_{k} contains the initial residual vector even for k = 0; */
210
/* Also assume that rnorm = || B*r_{k} || and B*r_{k} are already */
211
/* computed by the calling program. */
213
/* betaj = rnorm ; p_{k+1} = B*r_{k} ; */
214
/* For j = k+1, ..., k+np Do */
215
/* 1) if ( betaj < tol ) stop or restart depending on j. */
216
/* if ( restart ) generate a new starting vector. */
217
/* 2) v_{j} = r(j-1)/betaj; V_{j} = [V_{j-1}, v_{j}]; */
218
/* p_{j} = p_{j}/betaj */
219
/* 3) r_{j} = OP*v_{j} where OP is defined as in dsaupd */
220
/* For shift-invert mode p_{j} = B*v_{j} is already available. */
221
/* wnorm = || OP*v_{j} || */
222
/* 4) Compute the j-th step residual vector. */
223
/* w_{j} = V_{j}^T * B * OP * v_{j} */
224
/* r_{j} = OP*v_{j} - V_{j} * w_{j} */
225
/* alphaj <- j-th component of w_{j} */
226
/* rnorm = || r_{j} || */
227
/* betaj+1 = rnorm */
228
/* If (rnorm > 0.717*wnorm) accept step and go back to 1) */
229
/* 5) Re-orthogonalization step: */
230
/* s = V_{j}'*B*r_{j} */
231
/* r_{j} = r_{j} - V_{j}*s; rnorm1 = || r_{j} || */
232
/* alphaj = alphaj + s_{j}; */
233
/* 6) Iterative refinement step: */
234
/* If (rnorm1 > 0.717*rnorm) then */
236
/* accept step and go back to 1) */
239
/* If this is the first time in step 6), go to 5) */
240
/* Else r_{j} lies in the span of V_{j} numerically. */
241
/* Set r_{j} = 0 and rnorm = 0; go to 1) */
247
/* ----------------------------------------------------------------------- */
249
/* Subroutine */ int igraphdsaitr_(integer *ido, char *bmat, integer *n,
250
integer *k, integer *np, integer *mode, doublereal *resid, doublereal
251
*rnorm, doublereal *v, integer *ldv, doublereal *h__, integer *ldh,
252
integer *ipntr, doublereal *workd, integer *info)
254
/* Initialized data */
256
static logical first = TRUE_;
258
/* System generated locals */
259
integer h_dim1, h_offset, v_dim1, v_offset, i__1;
261
/* Builtin functions */
262
double sqrt(doublereal);
264
/* Local variables */
265
static integer i__, j;
266
static real t0, t1, t2, t3, t4, t5;
268
extern doublereal igraphddot_(integer *, doublereal *, integer *,
269
doublereal *, integer *), igraphdnrm2_(integer *, doublereal *,
271
static integer ipj, irj, ivj;
272
extern /* Subroutine */ int igraphdscal_(integer *, doublereal *,
273
doublereal *, integer *), igraphdgemv_(char *, integer *, integer
274
*, doublereal *, doublereal *, integer *, doublereal *, integer *,
275
doublereal *, doublereal *, integer *), igraphdcopy_(
276
integer *, doublereal *, integer *, doublereal *, integer *),
277
igraphdvout_(integer *, integer *, doublereal *, integer *, char *
278
), igraphivout_(integer *, integer *, integer *, integer *
279
, char *), igraphdgetv0_(integer *, char *, integer *,
280
logical *, integer *, integer *, doublereal *, integer *,
281
doublereal *, doublereal *, integer *, doublereal *, integer *);
282
static integer ierr, iter;
283
extern doublereal igraphdlamch_(char *);
285
extern /* Subroutine */ int igraphdlascl_(char *, integer *, integer *,
286
doublereal *, doublereal *, integer *, integer *, doublereal *,
287
integer *, integer *), igraphsecond_(real *);
288
static doublereal temp1;
289
static logical orth1, orth2, step3, step4;
290
static integer infol;
291
static doublereal xtemp[2], wnorm, rnorm1, safmin;
292
static logical rstart;
293
static integer msglvl;
296
/* %----------------------------------------------------% */
297
/* | Include files for debugging and timing information | */
298
/* %----------------------------------------------------% */
301
/* \SCCS Information: @(#) */
302
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
304
/* %---------------------------------% */
305
/* | See debug.doc for documentation | */
306
/* %---------------------------------% */
308
/* %------------------% */
309
/* | Scalar Arguments | */
310
/* %------------------% */
312
/* %--------------------------------% */
313
/* | See stat.doc for documentation | */
314
/* %--------------------------------% */
316
/* \SCCS Information: @(#) */
317
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
321
/* %-----------------% */
322
/* | Array Arguments | */
323
/* %-----------------% */
331
/* %---------------% */
332
/* | Local Scalars | */
333
/* %---------------% */
336
/* %-----------------------% */
337
/* | Local Array Arguments | */
338
/* %-----------------------% */
341
/* %----------------------% */
342
/* | External Subroutines | */
343
/* %----------------------% */
346
/* %--------------------% */
347
/* | External Functions | */
348
/* %--------------------% */
351
/* %-----------------% */
352
/* | Data statements | */
353
/* %-----------------% */
355
/* Parameter adjustments */
359
v_offset = 1 + v_dim1;
362
h_offset = 1 + h_dim1;
368
/* %-----------------------% */
369
/* | Executable Statements | */
370
/* %-----------------------% */
375
/* %--------------------------------% */
376
/* | safmin = safe minimum is such | */
377
/* | that 1/sfmin does not overflow | */
378
/* %--------------------------------% */
380
safmin = igraphdlamch_("safmin");
385
/* %-------------------------------% */
386
/* | Initialize timing statistics | */
387
/* | & message level for debugging | */
388
/* %-------------------------------% */
391
msglvl = debug_1.msaitr;
393
/* %------------------------------% */
394
/* | Initial call to this routine | */
395
/* %------------------------------% */
404
/* %--------------------------------% */
405
/* | Pointer to the current step of | */
406
/* | the factorization to build | */
407
/* %--------------------------------% */
411
/* %------------------------------------------% */
412
/* | Pointers used for reverse communication | */
413
/* | when using WORKD. | */
414
/* %------------------------------------------% */
421
/* %-------------------------------------------------% */
422
/* | When in reverse communication mode one of: | */
423
/* | STEP3, STEP4, ORTH1, ORTH2, RSTART | */
424
/* | will be .true. | */
425
/* | STEP3: return from computing OP*v_{j}. | */
426
/* | STEP4: return from computing B-norm of OP*v_{j} | */
427
/* | ORTH1: return from computing B-norm of r_{j+1} | */
428
/* | ORTH2: return from computing B-norm of | */
429
/* | correction to the residual vector. | */
430
/* | RSTART: return from OP computations needed by | */
432
/* %-------------------------------------------------% */
450
/* %------------------------------% */
451
/* | Else this is the first step. | */
452
/* %------------------------------% */
454
/* %--------------------------------------------------------------% */
456
/* | A R N O L D I I T E R A T I O N L O O P | */
458
/* | Note: B*r_{j-1} is already in WORKD(1:N)=WORKD(IPJ:IPJ+N-1) | */
459
/* %--------------------------------------------------------------% */
464
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_saitr: g"
465
"enerating Arnoldi vector no.");
466
igraphdvout_(&debug_1.logfil, &c__1, rnorm, &debug_1.ndigit, "_saitr"
467
": B-norm of the current residual =");
470
/* %---------------------------------------------------------% */
471
/* | Check for exact zero. Equivalent to determing whether a | */
472
/* | j-step Arnoldi factorization is present. | */
473
/* %---------------------------------------------------------% */
479
/* %---------------------------------------------------% */
480
/* | Invariant subspace found, generate a new starting | */
481
/* | vector which is orthogonal to the current Arnoldi | */
482
/* | basis and continue the iteration. | */
483
/* %---------------------------------------------------% */
486
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_saitr: *"
487
"***** restart at step ******");
490
/* %---------------------------------------------% */
491
/* | ITRY is the loop variable that controls the | */
492
/* | maximum amount of times that a restart is | */
493
/* | attempted. NRSTRT is used by stat.h | */
494
/* %---------------------------------------------% */
503
/* %--------------------------------------% */
504
/* | If in reverse communication mode and | */
505
/* | RSTART = .true. flow returns here. | */
506
/* %--------------------------------------% */
508
igraphdgetv0_(ido, bmat, &itry, &c_false, n, &j, &v[v_offset], ldv, &
509
resid[1], rnorm, &ipntr[1], &workd[1], &ierr);
519
/* %------------------------------------------------% */
520
/* | Give up after several restart attempts. | */
521
/* | Set INFO to the size of the invariant subspace | */
522
/* | which spans OP and exit. | */
523
/* %------------------------------------------------% */
527
timing_1.tsaitr += t1 - t0;
534
/* %---------------------------------------------------------% */
535
/* | STEP 2: v_{j} = r_{j-1}/rnorm and p_{j} = p_{j}/rnorm | */
536
/* | Note that p_{j} = B*r_{j-1}. In order to avoid overflow | */
537
/* | when reciprocating a small RNORM, test against lower | */
538
/* | machine bound. | */
539
/* %---------------------------------------------------------% */
541
igraphdcopy_(n, &resid[1], &c__1, &v[j * v_dim1 + 1], &c__1);
542
if (*rnorm >= safmin) {
544
igraphdscal_(n, &temp1, &v[j * v_dim1 + 1], &c__1);
545
igraphdscal_(n, &temp1, &workd[ipj], &c__1);
548
/* %-----------------------------------------% */
549
/* | To scale both v_{j} and p_{j} carefully | */
550
/* | use LAPACK routine SLASCL | */
551
/* %-----------------------------------------% */
553
igraphdlascl_("General", &i__, &i__, rnorm, &c_b24, n, &c__1, &v[j *
554
v_dim1 + 1], n, &infol);
555
igraphdlascl_("General", &i__, &i__, rnorm, &c_b24, n, &c__1, &workd[
559
/* %------------------------------------------------------% */
560
/* | STEP 3: r_{j} = OP*v_{j}; Note that p_{j} = B*v_{j} | */
561
/* | Note that this is not quite yet r_{j}. See STEP 4 | */
562
/* %------------------------------------------------------% */
567
igraphdcopy_(n, &v[j * v_dim1 + 1], &c__1, &workd[ivj], &c__1);
573
/* %-----------------------------------% */
574
/* | Exit in order to compute OP*v_{j} | */
575
/* %-----------------------------------% */
580
/* %-----------------------------------% */
581
/* | Back from reverse communication; | */
582
/* | WORKD(IRJ:IRJ+N-1) := OP*v_{j}. | */
583
/* %-----------------------------------% */
586
timing_1.tmvopx += t3 - t2;
590
/* %------------------------------------------% */
591
/* | Put another copy of OP*v_{j} into RESID. | */
592
/* %------------------------------------------% */
594
igraphdcopy_(n, &workd[irj], &c__1, &resid[1], &c__1);
596
/* %-------------------------------------------% */
597
/* | STEP 4: Finish extending the symmetric | */
598
/* | Arnoldi to length j. If MODE = 2 | */
599
/* | then B*OP = B*inv(B)*A = A and | */
600
/* | we don't need to compute B*OP. | */
601
/* | NOTE: If MODE = 2 WORKD(IVJ:IVJ+N-1) is | */
602
/* | assumed to have A*v_{j}. | */
603
/* %-------------------------------------------% */
609
if (*(unsigned char *)bmat == 'G') {
616
/* %-------------------------------------% */
617
/* | Exit in order to compute B*OP*v_{j} | */
618
/* %-------------------------------------% */
621
} else if (*(unsigned char *)bmat == 'I') {
622
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
626
/* %-----------------------------------% */
627
/* | Back from reverse communication; | */
628
/* | WORKD(IPJ:IPJ+N-1) := B*OP*v_{j}. | */
629
/* %-----------------------------------% */
631
if (*(unsigned char *)bmat == 'G') {
633
timing_1.tmvbx += t3 - t2;
638
/* %-------------------------------------% */
639
/* | The following is needed for STEP 5. | */
640
/* | Compute the B-norm of OP*v_{j}. | */
641
/* %-------------------------------------% */
646
/* %----------------------------------% */
647
/* | Note that the B-norm of OP*v_{j} | */
648
/* | is the inv(B)-norm of A*v_{j}. | */
649
/* %----------------------------------% */
651
wnorm = igraphddot_(n, &resid[1], &c__1, &workd[ivj], &c__1);
652
wnorm = sqrt((abs(wnorm)));
653
} else if (*(unsigned char *)bmat == 'G') {
654
wnorm = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
655
wnorm = sqrt((abs(wnorm)));
656
} else if (*(unsigned char *)bmat == 'I') {
657
wnorm = igraphdnrm2_(n, &resid[1], &c__1);
660
/* %-----------------------------------------% */
661
/* | Compute the j-th residual corresponding | */
662
/* | to the j step factorization. | */
663
/* | Use Classical Gram Schmidt and compute: | */
664
/* | w_{j} <- V_{j}^T * B * OP * v_{j} | */
665
/* | r_{j} <- OP*v_{j} - V_{j} * w_{j} | */
666
/* %-----------------------------------------% */
669
/* %------------------------------------------% */
670
/* | Compute the j Fourier coefficients w_{j} | */
671
/* | WORKD(IPJ:IPJ+N-1) contains B*OP*v_{j}. | */
672
/* %------------------------------------------% */
675
igraphdgemv_("T", n, &j, &c_b24, &v[v_offset], ldv, &workd[ipj], &
676
c__1, &c_b49, &workd[irj], &c__1);
677
} else if (*mode == 2) {
678
igraphdgemv_("T", n, &j, &c_b24, &v[v_offset], ldv, &workd[ivj], &
679
c__1, &c_b49, &workd[irj], &c__1);
682
/* %--------------------------------------% */
683
/* | Orthgonalize r_{j} against V_{j}. | */
684
/* | RESID contains OP*v_{j}. See STEP 3. | */
685
/* %--------------------------------------% */
687
igraphdgemv_("N", n, &j, &c_b57, &v[v_offset], ldv, &workd[irj], &c__1, &
688
c_b24, &resid[1], &c__1);
690
/* %--------------------------------------% */
691
/* | Extend H to have j rows and columns. | */
692
/* %--------------------------------------% */
694
h__[j + (h_dim1 << 1)] = workd[irj + j - 1];
695
if (j == 1 || rstart) {
696
h__[j + h_dim1] = 0.;
698
h__[j + h_dim1] = *rnorm;
706
if (*(unsigned char *)bmat == 'G') {
708
igraphdcopy_(n, &resid[1], &c__1, &workd[irj], &c__1);
713
/* %----------------------------------% */
714
/* | Exit in order to compute B*r_{j} | */
715
/* %----------------------------------% */
718
} else if (*(unsigned char *)bmat == 'I') {
719
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
723
/* %---------------------------------------------------% */
724
/* | Back from reverse communication if ORTH1 = .true. | */
725
/* | WORKD(IPJ:IPJ+N-1) := B*r_{j}. | */
726
/* %---------------------------------------------------% */
728
if (*(unsigned char *)bmat == 'G') {
730
timing_1.tmvbx += t3 - t2;
735
/* %------------------------------% */
736
/* | Compute the B-norm of r_{j}. | */
737
/* %------------------------------% */
739
if (*(unsigned char *)bmat == 'G') {
740
*rnorm = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
741
*rnorm = sqrt((abs(*rnorm)));
742
} else if (*(unsigned char *)bmat == 'I') {
743
*rnorm = igraphdnrm2_(n, &resid[1], &c__1);
746
/* %-----------------------------------------------------------% */
747
/* | STEP 5: Re-orthogonalization / Iterative refinement phase | */
748
/* | Maximum NITER_ITREF tries. | */
750
/* | s = V_{j}^T * B * r_{j} | */
751
/* | r_{j} = r_{j} - V_{j}*s | */
752
/* | alphaj = alphaj + s_{j} | */
754
/* | The stopping criteria used for iterative refinement is | */
755
/* | discussed in Parlett's book SEP, page 107 and in Gragg & | */
756
/* | Reichel ACM TOMS paper; Algorithm 686, Dec. 1990. | */
757
/* | Determine if we need to correct the residual. The goal is | */
758
/* | to enforce ||v(:,1:j)^T * r_{j}|| .le. eps * || r_{j} || | */
759
/* %-----------------------------------------------------------% */
761
if (*rnorm > wnorm * .717f) {
766
/* %---------------------------------------------------% */
767
/* | Enter the Iterative refinement phase. If further | */
768
/* | refinement is necessary, loop back here. The loop | */
769
/* | variable is ITER. Perform a step of Classical | */
770
/* | Gram-Schmidt using all the Arnoldi vectors V_{j} | */
771
/* %---------------------------------------------------% */
778
igraphdvout_(&debug_1.logfil, &c__2, xtemp, &debug_1.ndigit, "_saitr"
779
": re-orthonalization ; wnorm and rnorm are");
782
/* %----------------------------------------------------% */
783
/* | Compute V_{j}^T * B * r_{j}. | */
784
/* | WORKD(IRJ:IRJ+J-1) = v(:,1:J)'*WORKD(IPJ:IPJ+N-1). | */
785
/* %----------------------------------------------------% */
787
igraphdgemv_("T", n, &j, &c_b24, &v[v_offset], ldv, &workd[ipj], &c__1, &
788
c_b49, &workd[irj], &c__1);
790
/* %----------------------------------------------% */
791
/* | Compute the correction to the residual: | */
792
/* | r_{j} = r_{j} - V_{j} * WORKD(IRJ:IRJ+J-1). | */
793
/* | The correction to H is v(:,1:J)*H(1:J,1:J) + | */
794
/* | v(:,1:J)*WORKD(IRJ:IRJ+J-1)*e'_j, but only | */
795
/* | H(j,j) is updated. | */
796
/* %----------------------------------------------% */
798
igraphdgemv_("N", n, &j, &c_b57, &v[v_offset], ldv, &workd[irj], &c__1, &
799
c_b24, &resid[1], &c__1);
801
if (j == 1 || rstart) {
802
h__[j + h_dim1] = 0.;
804
h__[j + (h_dim1 << 1)] += workd[irj + j - 1];
808
if (*(unsigned char *)bmat == 'G') {
810
igraphdcopy_(n, &resid[1], &c__1, &workd[irj], &c__1);
815
/* %-----------------------------------% */
816
/* | Exit in order to compute B*r_{j}. | */
817
/* | r_{j} is the corrected residual. | */
818
/* %-----------------------------------% */
821
} else if (*(unsigned char *)bmat == 'I') {
822
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
826
/* %---------------------------------------------------% */
827
/* | Back from reverse communication if ORTH2 = .true. | */
828
/* %---------------------------------------------------% */
830
if (*(unsigned char *)bmat == 'G') {
832
timing_1.tmvbx += t3 - t2;
835
/* %-----------------------------------------------------% */
836
/* | Compute the B-norm of the corrected residual r_{j}. | */
837
/* %-----------------------------------------------------% */
839
if (*(unsigned char *)bmat == 'G') {
840
rnorm1 = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
841
rnorm1 = sqrt((abs(rnorm1)));
842
} else if (*(unsigned char *)bmat == 'I') {
843
rnorm1 = igraphdnrm2_(n, &resid[1], &c__1);
846
if (msglvl > 0 && iter > 0) {
847
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_saitr: I"
848
"terative refinement for Arnoldi residual");
852
igraphdvout_(&debug_1.logfil, &c__2, xtemp, &debug_1.ndigit,
853
"_saitr: iterative refinement ; rnorm and rnorm1 are");
857
/* %-----------------------------------------% */
858
/* | Determine if we need to perform another | */
859
/* | step of re-orthogonalization. | */
860
/* %-----------------------------------------% */
862
if (rnorm1 > *rnorm * .717f) {
864
/* %--------------------------------% */
865
/* | No need for further refinement | */
866
/* %--------------------------------% */
872
/* %-------------------------------------------% */
873
/* | Another step of iterative refinement step | */
874
/* | is required. NITREF is used by stat.h | */
875
/* %-------------------------------------------% */
884
/* %-------------------------------------------------% */
885
/* | Otherwise RESID is numerically in the span of V | */
886
/* %-------------------------------------------------% */
889
for (jj = 1; jj <= i__1; ++jj) {
896
/* %----------------------------------------------% */
897
/* | Branch here directly if iterative refinement | */
898
/* | wasn't necessary or after at most NITER_REF | */
899
/* | steps of iterative refinement. | */
900
/* %----------------------------------------------% */
908
timing_1.titref += t5 - t4;
910
/* %----------------------------------------------------------% */
911
/* | Make sure the last off-diagonal element is non negative | */
912
/* | If not perform a similarity transformation on H(1:j,1:j) | */
913
/* | and scale v(:,j) by -1. | */
914
/* %----------------------------------------------------------% */
916
if (h__[j + h_dim1] < 0.) {
917
h__[j + h_dim1] = -h__[j + h_dim1];
919
igraphdscal_(n, &c_b57, &v[(j + 1) * v_dim1 + 1], &c__1);
921
igraphdscal_(n, &c_b57, &resid[1], &c__1);
925
/* %------------------------------------% */
926
/* | STEP 6: Update j = j+1; Continue | */
927
/* %------------------------------------% */
932
timing_1.tsaitr += t1 - t0;
937
igraphdvout_(&debug_1.logfil, &i__1, &h__[(h_dim1 << 1) + 1], &
938
debug_1.ndigit, "_saitr: main diagonal of matrix H of st"
942
igraphdvout_(&debug_1.logfil, &i__1, &h__[h_dim1 + 2], &
943
debug_1.ndigit, "_saitr: sub diagonal of matrix H of"
951
/* %--------------------------------------------------------% */
952
/* | Loop back to extend the factorization by another step. | */
953
/* %--------------------------------------------------------% */
957
/* %---------------------------------------------------------------% */
959
/* | E N D O F M A I N I T E R A T I O N L O O P | */
961
/* %---------------------------------------------------------------% */
966
/* %---------------% */
967
/* | End of dsaitr | */
968
/* %---------------% */
970
} /* igraphdsaitr_ */
added
removed
src/arpack/dsaitr.c
