1
/* igraphdnaitr.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 integer c__1 = 1;
40
static logical c_false = FALSE_;
41
static doublereal c_b25 = 1.;
42
static doublereal c_b47 = 0.;
43
static doublereal c_b50 = -1.;
44
static integer c__2 = 2;
46
/* ----------------------------------------------------------------------- */
49
/* \Name: igraphdnaitr */
52
/* Reverse communication interface for applying NP additional steps to */
53
/* a K step nonsymmetric 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 igraphdnaupd. The B-norm of r_{k+p} is also */
64
/* computed and returned. */
67
/* call igraphdnaitr */
68
/* ( IDO, BMAT, N, K, NP, NB, 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 do not need to be */
92
/* recompute in forming OP * Q. */
94
/* BMAT Character*1. (INPUT) */
95
/* BMAT specifies the type of the matrix B that defines the */
96
/* semi-inner product for the operator OP. See igraphdnaupd. */
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 size of V and H. */
106
/* NP Integer. (INPUT) */
107
/* Number of additional Arnoldi steps to take. */
109
/* NB Integer. (INPUT) */
110
/* Blocksize to be used in the recurrence. */
111
/* Only work for NB = 1 right now. The goal is to have a */
112
/* program that implement both the block and non-block method. */
114
/* RESID Double precision array of length N. (INPUT/OUTPUT) */
115
/* On INPUT: RESID contains the residual vector r_{k}. */
116
/* On OUTPUT: RESID contains the residual vector r_{k+p}. */
118
/* RNORM Double precision scalar. (INPUT/OUTPUT) */
119
/* B-norm of the starting residual on input. */
120
/* B-norm of the updated residual r_{k+p} on output. */
122
/* V Double precision N by K+NP array. (INPUT/OUTPUT) */
123
/* On INPUT: V contains the Arnoldi vectors in the first K */
125
/* On OUTPUT: V contains the new NP Arnoldi vectors in the next */
126
/* NP columns. The first K columns are unchanged. */
128
/* LDV Integer. (INPUT) */
129
/* Leading dimension of V exactly as declared in the calling */
132
/* H Double precision (K+NP) by (K+NP) array. (INPUT/OUTPUT) */
133
/* H is used to store the generated upper Hessenberg matrix. */
135
/* LDH Integer. (INPUT) */
136
/* Leading dimension of H exactly as declared in the calling */
139
/* IPNTR Integer array of length 3. (OUTPUT) */
140
/* Pointer to mark the starting locations in the WORK for */
141
/* vectors used by the Arnoldi iteration. */
142
/* ------------------------------------------------------------- */
143
/* IPNTR(1): pointer to the current operand vector X. */
144
/* IPNTR(2): pointer to the current result vector Y. */
145
/* IPNTR(3): pointer to the vector B * X when used in the */
146
/* shift-and-invert mode. X is the current operand. */
147
/* ------------------------------------------------------------- */
149
/* WORKD Double precision work array of length 3*N. (REVERSE COMMUNICATION) */
150
/* Distributed array to be used in the basic Arnoldi iteration */
151
/* for reverse communication. The calling program should not */
152
/* use WORKD as temporary workspace during the iteration !!!!!! */
153
/* On input, WORKD(1:N) = B*RESID and is used to save some */
154
/* computation at the first step. */
156
/* INFO Integer. (OUTPUT) */
157
/* = 0: Normal exit. */
158
/* > 0: Size of the spanning invariant subspace of OP found. */
162
/* ----------------------------------------------------------------------- */
166
/* \Local variables: */
170
/* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
171
/* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
173
/* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
174
/* Restarted Arnoldi Iteration", Rice University Technical Report */
175
/* TR95-13, Department of Computational and Applied Mathematics. */
177
/* \Routines called: */
178
/* dgetv0 ARPACK routine to generate the initial vector. */
179
/* ivout ARPACK utility routine that prints integers. */
180
/* second ARPACK utility routine for timing. */
181
/* dmout ARPACK utility routine that prints matrices */
182
/* dvout ARPACK utility routine that prints vectors. */
183
/* dlabad LAPACK routine that computes machine constants. */
184
/* dlamch LAPACK routine that determines machine constants. */
185
/* dlascl LAPACK routine for careful scaling of a matrix. */
186
/* dlanhs LAPACK routine that computes various norms of a matrix. */
187
/* dgemv Level 2 BLAS routine for matrix vector multiplication. */
188
/* daxpy Level 1 BLAS that computes a vector triad. */
189
/* dscal Level 1 BLAS that scales a vector. */
190
/* dcopy Level 1 BLAS that copies one vector to another . */
191
/* ddot Level 1 BLAS that computes the scalar product of two vectors. */
192
/* dnrm2 Level 1 BLAS that computes the norm of a vector. */
195
/* Danny Sorensen Phuong Vu */
196
/* Richard Lehoucq CRPC / Rice University */
197
/* Dept. of Computational & Houston, Texas */
198
/* Applied Mathematics */
199
/* Rice University */
202
/* \Revision history: */
203
/* xx/xx/92: Version ' 2.4' */
205
/* \SCCS Information: @(#) */
206
/* FILE: naitr.F SID: 2.4 DATE OF SID: 8/27/96 RELEASE: 2 */
209
/* The algorithm implemented is: */
211
/* restart = .false. */
212
/* Given V_{k} = [v_{1}, ..., v_{k}], r_{k}; */
213
/* r_{k} contains the initial residual vector even for k = 0; */
214
/* Also assume that rnorm = || B*r_{k} || and B*r_{k} are already */
215
/* computed by the calling program. */
217
/* betaj = rnorm ; p_{k+1} = B*r_{k} ; */
218
/* For j = k+1, ..., k+np Do */
219
/* 1) if ( betaj < tol ) stop or restart depending on j. */
220
/* ( At present tol is zero ) */
221
/* if ( restart ) generate a new starting vector. */
222
/* 2) v_{j} = r(j-1)/betaj; V_{j} = [V_{j-1}, v_{j}]; */
223
/* p_{j} = p_{j}/betaj */
224
/* 3) r_{j} = OP*v_{j} where OP is defined as in igraphdnaupd */
225
/* For shift-invert mode p_{j} = B*v_{j} is already available. */
226
/* wnorm = || OP*v_{j} || */
227
/* 4) Compute the j-th step residual vector. */
228
/* w_{j} = V_{j}^T * B * OP * v_{j} */
229
/* r_{j} = OP*v_{j} - V_{j} * w_{j} */
230
/* H(:,j) = w_{j}; */
231
/* H(j,j-1) = rnorm */
232
/* rnorm = || r_(j) || */
233
/* If (rnorm > 0.717*wnorm) accept step and go back to 1) */
234
/* 5) Re-orthogonalization step: */
235
/* s = V_{j}'*B*r_{j} */
236
/* r_{j} = r_{j} - V_{j}*s; rnorm1 = || r_{j} || */
237
/* alphaj = alphaj + s_{j}; */
238
/* 6) Iterative refinement step: */
239
/* If (rnorm1 > 0.717*rnorm) then */
241
/* accept step and go back to 1) */
244
/* If this is the first time in step 6), go to 5) */
245
/* Else r_{j} lies in the span of V_{j} numerically. */
246
/* Set r_{j} = 0 and rnorm = 0; go to 1) */
252
/* ----------------------------------------------------------------------- */
254
/* Subroutine */ int igraphdnaitr_(integer *ido, char *bmat, integer *n, integer *k,
255
integer *np, integer *nb, doublereal *resid, doublereal *rnorm,
256
doublereal *v, integer *ldv, doublereal *h__, integer *ldh, integer *
257
ipntr, doublereal *workd, integer *info)
259
/* Initialized data */
261
static logical first = TRUE_;
263
/* System generated locals */
264
integer h_dim1, h_offset, v_dim1, v_offset, i__1, i__2;
265
doublereal d__1, d__2;
267
/* Builtin functions */
268
double sqrt(doublereal);
270
/* Local variables */
271
static integer i__, j;
272
static real t0, t1, t2, t3, t4, t5;
273
static integer jj, ipj, irj, ivj;
274
static doublereal ulp, tst1;
275
extern doublereal igraphddot_(integer *, doublereal *, integer *, doublereal *,
277
static integer ierr, iter;
278
static doublereal unfl, ovfl;
280
extern doublereal igraphdnrm2_(integer *, doublereal *, integer *);
281
static doublereal temp1;
282
static logical orth1, orth2, step3, step4;
283
static doublereal betaj;
284
extern /* Subroutine */ int igraphdscal_(integer *, doublereal *, doublereal *,
285
integer *), igraphdgemv_(char *, integer *, integer *, doublereal *,
286
doublereal *, integer *, doublereal *, integer *, doublereal *,
287
doublereal *, integer *);
288
static integer infol;
289
extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *,
290
doublereal *, integer *), igraphdaxpy_(integer *, doublereal *,
291
doublereal *, integer *, doublereal *, integer *), igraphdmout_(integer
292
*, integer *, integer *, doublereal *, integer *, integer *, char
294
static doublereal xtemp[2];
295
extern /* Subroutine */ int igraphdvout_(integer *, integer *, doublereal *,
297
static doublereal wnorm;
298
extern /* Subroutine */ int igraphivout_(integer *, integer *, integer *,
299
integer *, char *), igraphdgetv0_(integer *, char *, integer *,
300
logical *, integer *, integer *, doublereal *, integer *,
301
doublereal *, doublereal *, integer *, doublereal *, integer *
302
), igraphdlabad_(doublereal *, doublereal *);
303
static doublereal rnorm1;
304
extern doublereal igraphdlamch_(char *);
305
extern /* Subroutine */ int igraphdlascl_(char *, integer *, integer *,
306
doublereal *, doublereal *, integer *, integer *, doublereal *,
307
integer *, integer *);
308
extern doublereal igraphdlanhs_(char *, integer *, doublereal *, integer *,
310
extern /* Subroutine */ int igraphsecond_(real *);
311
static logical rstart;
312
static integer msglvl;
313
static doublereal smlnum;
316
/* %----------------------------------------------------% */
317
/* | Include files for debugging and timing information | */
318
/* %----------------------------------------------------% */
321
/* \SCCS Information: @(#) */
322
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
324
/* %---------------------------------% */
325
/* | See debug.doc for documentation | */
326
/* %---------------------------------% */
328
/* %------------------% */
329
/* | Scalar Arguments | */
330
/* %------------------% */
332
/* %--------------------------------% */
333
/* | See stat.doc for documentation | */
334
/* %--------------------------------% */
336
/* \SCCS Information: @(#) */
337
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
341
/* %-----------------% */
342
/* | Array Arguments | */
343
/* %-----------------% */
351
/* %---------------% */
352
/* | Local Scalars | */
353
/* %---------------% */
356
/* %-----------------------% */
357
/* | Local Array Arguments | */
358
/* %-----------------------% */
361
/* %----------------------% */
362
/* | External Subroutines | */
363
/* %----------------------% */
366
/* %--------------------% */
367
/* | External Functions | */
368
/* %--------------------% */
371
/* %---------------------% */
372
/* | Intrinsic Functions | */
373
/* %---------------------% */
376
/* %-----------------% */
377
/* | Data statements | */
378
/* %-----------------% */
380
/* Parameter adjustments */
384
v_offset = 1 + v_dim1;
387
h_offset = 1 + h_dim1;
393
/* %-----------------------% */
394
/* | Executable Statements | */
395
/* %-----------------------% */
399
/* %-----------------------------------------% */
400
/* | Set machine-dependent constants for the | */
401
/* | the splitting and deflation criterion. | */
402
/* | If norm(H) <= sqrt(OVFL), | */
403
/* | overflow should not occur. | */
404
/* | REFERENCE: LAPACK subroutine dlahqr | */
405
/* %-----------------------------------------% */
407
unfl = igraphdlamch_("safe minimum");
409
igraphdlabad_(&unfl, &ovfl);
410
ulp = igraphdlamch_("precision");
411
smlnum = unfl * (*n / ulp);
417
/* %-------------------------------% */
418
/* | Initialize timing statistics | */
419
/* | & message level for debugging | */
420
/* %-------------------------------% */
423
msglvl = debug_1.mnaitr;
425
/* %------------------------------% */
426
/* | Initial call to this routine | */
427
/* %------------------------------% */
441
/* %-------------------------------------------------% */
442
/* | When in reverse communication mode one of: | */
443
/* | STEP3, STEP4, ORTH1, ORTH2, RSTART | */
444
/* | will be .true. when .... | */
445
/* | STEP3: return from computing OP*v_{j}. | */
446
/* | STEP4: return from computing B-norm of OP*v_{j} | */
447
/* | ORTH1: return from computing B-norm of r_{j+1} | */
448
/* | ORTH2: return from computing B-norm of | */
449
/* | correction to the residual vector. | */
450
/* | RSTART: return from OP computations needed by | */
452
/* %-------------------------------------------------% */
470
/* %-----------------------------% */
471
/* | Else this is the first step | */
472
/* %-----------------------------% */
474
/* %--------------------------------------------------------------% */
476
/* | A R N O L D I I T E R A T I O N L O O P | */
478
/* | Note: B*r_{j-1} is already in WORKD(1:N)=WORKD(IPJ:IPJ+N-1) | */
479
/* %--------------------------------------------------------------% */
483
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_naitr: generat"
484
"ing Arnoldi vector number");
485
igraphdvout_(&debug_1.logfil, &c__1, rnorm, &debug_1.ndigit, "_naitr: B-no"
486
"rm of the current residual is");
489
/* %---------------------------------------------------% */
490
/* | STEP 1: Check if the B norm of j-th residual | */
491
/* | vector is zero. Equivalent to determing whether | */
492
/* | an exact j-step Arnoldi factorization is present. | */
493
/* %---------------------------------------------------% */
500
/* %---------------------------------------------------% */
501
/* | Invariant subspace found, generate a new starting | */
502
/* | vector which is orthogonal to the current Arnoldi | */
503
/* | basis and continue the iteration. | */
504
/* %---------------------------------------------------% */
507
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_naitr: ****** "
508
"RESTART AT STEP ******");
511
/* %---------------------------------------------% */
512
/* | ITRY is the loop variable that controls the | */
513
/* | maximum amount of times that a restart is | */
514
/* | attempted. NRSTRT is used by stat.h | */
515
/* %---------------------------------------------% */
525
/* %--------------------------------------% */
526
/* | If in reverse communication mode and | */
527
/* | RSTART = .true. flow returns here. | */
528
/* %--------------------------------------% */
530
igraphdgetv0_(ido, bmat, &itry, &c_false, n, &j, &v[v_offset], ldv, &resid[1],
531
rnorm, &ipntr[1], &workd[1], &ierr);
541
/* %------------------------------------------------% */
542
/* | Give up after several restart attempts. | */
543
/* | Set INFO to the size of the invariant subspace | */
544
/* | which spans OP and exit. | */
545
/* %------------------------------------------------% */
549
timing_1.tnaitr += t1 - t0;
556
/* %---------------------------------------------------------% */
557
/* | STEP 2: v_{j} = r_{j-1}/rnorm and p_{j} = p_{j}/rnorm | */
558
/* | Note that p_{j} = B*r_{j-1}. In order to avoid overflow | */
559
/* | when reciprocating a small RNORM, test against lower | */
560
/* | machine bound. | */
561
/* %---------------------------------------------------------% */
563
igraphdcopy_(n, &resid[1], &c__1, &v[j * v_dim1 + 1], &c__1);
564
if (*rnorm >= unfl) {
566
igraphdscal_(n, &temp1, &v[j * v_dim1 + 1], &c__1);
567
igraphdscal_(n, &temp1, &workd[ipj], &c__1);
570
/* %-----------------------------------------% */
571
/* | To scale both v_{j} and p_{j} carefully | */
572
/* | use LAPACK routine SLASCL | */
573
/* %-----------------------------------------% */
575
igraphdlascl_("General", &i__, &i__, rnorm, &c_b25, n, &c__1, &v[j * v_dim1
577
igraphdlascl_("General", &i__, &i__, rnorm, &c_b25, n, &c__1, &workd[ipj],
581
/* %------------------------------------------------------% */
582
/* | STEP 3: r_{j} = OP*v_{j}; Note that p_{j} = B*v_{j} | */
583
/* | Note that this is not quite yet r_{j}. See STEP 4 | */
584
/* %------------------------------------------------------% */
589
igraphdcopy_(n, &v[j * v_dim1 + 1], &c__1, &workd[ivj], &c__1);
595
/* %-----------------------------------% */
596
/* | Exit in order to compute OP*v_{j} | */
597
/* %-----------------------------------% */
602
/* %----------------------------------% */
603
/* | Back from reverse communication; | */
604
/* | WORKD(IRJ:IRJ+N-1) := OP*v_{j} | */
605
/* | if step3 = .true. | */
606
/* %----------------------------------% */
609
timing_1.tmvopx += t3 - t2;
612
/* %------------------------------------------% */
613
/* | Put another copy of OP*v_{j} into RESID. | */
614
/* %------------------------------------------% */
616
igraphdcopy_(n, &workd[irj], &c__1, &resid[1], &c__1);
618
/* %---------------------------------------% */
619
/* | STEP 4: Finish extending the Arnoldi | */
620
/* | factorization to length j. | */
621
/* %---------------------------------------% */
624
if (*(unsigned char *)bmat == 'G') {
631
/* %-------------------------------------% */
632
/* | Exit in order to compute B*OP*v_{j} | */
633
/* %-------------------------------------% */
636
} else if (*(unsigned char *)bmat == 'I') {
637
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
641
/* %----------------------------------% */
642
/* | Back from reverse communication; | */
643
/* | WORKD(IPJ:IPJ+N-1) := B*OP*v_{j} | */
644
/* | if step4 = .true. | */
645
/* %----------------------------------% */
647
if (*(unsigned char *)bmat == 'G') {
649
timing_1.tmvbx += t3 - t2;
654
/* %-------------------------------------% */
655
/* | The following is needed for STEP 5. | */
656
/* | Compute the B-norm of OP*v_{j}. | */
657
/* %-------------------------------------% */
659
if (*(unsigned char *)bmat == 'G') {
660
wnorm = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
661
wnorm = sqrt((abs(wnorm)));
662
} else if (*(unsigned char *)bmat == 'I') {
663
wnorm = igraphdnrm2_(n, &resid[1], &c__1);
666
/* %-----------------------------------------% */
667
/* | Compute the j-th residual corresponding | */
668
/* | to the j step factorization. | */
669
/* | Use Classical Gram Schmidt and compute: | */
670
/* | w_{j} <- V_{j}^T * B * OP * v_{j} | */
671
/* | r_{j} <- OP*v_{j} - V_{j} * w_{j} | */
672
/* %-----------------------------------------% */
675
/* %------------------------------------------% */
676
/* | Compute the j Fourier coefficients w_{j} | */
677
/* | WORKD(IPJ:IPJ+N-1) contains B*OP*v_{j}. | */
678
/* %------------------------------------------% */
680
igraphdgemv_("T", n, &j, &c_b25, &v[v_offset], ldv, &workd[ipj], &c__1, &c_b47,
681
&h__[j * h_dim1 + 1], &c__1);
683
/* %--------------------------------------% */
684
/* | Orthogonalize r_{j} against V_{j}. | */
685
/* | RESID contains OP*v_{j}. See STEP 3. | */
686
/* %--------------------------------------% */
688
igraphdgemv_("N", n, &j, &c_b50, &v[v_offset], ldv, &h__[j * h_dim1 + 1], &c__1,
689
&c_b25, &resid[1], &c__1);
692
h__[j + (j - 1) * h_dim1] = betaj;
700
if (*(unsigned char *)bmat == 'G') {
702
igraphdcopy_(n, &resid[1], &c__1, &workd[irj], &c__1);
707
/* %----------------------------------% */
708
/* | Exit in order to compute B*r_{j} | */
709
/* %----------------------------------% */
712
} else if (*(unsigned char *)bmat == 'I') {
713
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
717
/* %---------------------------------------------------% */
718
/* | Back from reverse communication if ORTH1 = .true. | */
719
/* | WORKD(IPJ:IPJ+N-1) := B*r_{j}. | */
720
/* %---------------------------------------------------% */
722
if (*(unsigned char *)bmat == 'G') {
724
timing_1.tmvbx += t3 - t2;
729
/* %------------------------------% */
730
/* | Compute the B-norm of r_{j}. | */
731
/* %------------------------------% */
733
if (*(unsigned char *)bmat == 'G') {
734
*rnorm = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
735
*rnorm = sqrt((abs(*rnorm)));
736
} else if (*(unsigned char *)bmat == 'I') {
737
*rnorm = igraphdnrm2_(n, &resid[1], &c__1);
740
/* %-----------------------------------------------------------% */
741
/* | STEP 5: Re-orthogonalization / Iterative refinement phase | */
742
/* | Maximum NITER_ITREF tries. | */
744
/* | s = V_{j}^T * B * r_{j} | */
745
/* | r_{j} = r_{j} - V_{j}*s | */
746
/* | alphaj = alphaj + s_{j} | */
748
/* | The stopping criteria used for iterative refinement is | */
749
/* | discussed in Parlett's book SEP, page 107 and in Gragg & | */
750
/* | Reichel ACM TOMS paper; Algorithm 686, Dec. 1990. | */
751
/* | Determine if we need to correct the residual. The goal is | */
752
/* | to enforce ||v(:,1:j)^T * r_{j}|| .le. eps * || r_{j} || | */
753
/* | The following test determines whether the sine of the | */
754
/* | angle between OP*x and the computed residual is less | */
755
/* | than or equal to 0.717. | */
756
/* %-----------------------------------------------------------% */
758
if (*rnorm > wnorm * .717f) {
764
/* %---------------------------------------------------% */
765
/* | Enter the Iterative refinement phase. If further | */
766
/* | refinement is necessary, loop back here. The loop | */
767
/* | variable is ITER. Perform a step of Classical | */
768
/* | Gram-Schmidt using all the Arnoldi vectors V_{j} | */
769
/* %---------------------------------------------------% */
776
igraphdvout_(&debug_1.logfil, &c__2, xtemp, &debug_1.ndigit, "_naitr: re-o"
777
"rthonalization; wnorm and rnorm are");
778
igraphdvout_(&debug_1.logfil, &j, &h__[j * h_dim1 + 1], &debug_1.ndigit,
779
"_naitr: j-th column of H");
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_b25, &v[v_offset], ldv, &workd[ipj], &c__1, &c_b47,
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. | */
795
/* %---------------------------------------------% */
797
igraphdgemv_("N", n, &j, &c_b50, &v[v_offset], ldv, &workd[irj], &c__1, &c_b25,
799
igraphdaxpy_(&j, &c_b25, &workd[irj], &c__1, &h__[j * h_dim1 + 1], &c__1);
803
if (*(unsigned char *)bmat == 'G') {
805
igraphdcopy_(n, &resid[1], &c__1, &workd[irj], &c__1);
810
/* %-----------------------------------% */
811
/* | Exit in order to compute B*r_{j}. | */
812
/* | r_{j} is the corrected residual. | */
813
/* %-----------------------------------% */
816
} else if (*(unsigned char *)bmat == 'I') {
817
igraphdcopy_(n, &resid[1], &c__1, &workd[ipj], &c__1);
821
/* %---------------------------------------------------% */
822
/* | Back from reverse communication if ORTH2 = .true. | */
823
/* %---------------------------------------------------% */
825
if (*(unsigned char *)bmat == 'G') {
827
timing_1.tmvbx += t3 - t2;
830
/* %-----------------------------------------------------% */
831
/* | Compute the B-norm of the corrected residual r_{j}. | */
832
/* %-----------------------------------------------------% */
834
if (*(unsigned char *)bmat == 'G') {
835
rnorm1 = igraphddot_(n, &resid[1], &c__1, &workd[ipj], &c__1);
836
rnorm1 = sqrt((abs(rnorm1)));
837
} else if (*(unsigned char *)bmat == 'I') {
838
rnorm1 = igraphdnrm2_(n, &resid[1], &c__1);
841
if (msglvl > 0 && iter > 0) {
842
igraphivout_(&debug_1.logfil, &c__1, &j, &debug_1.ndigit, "_naitr: Iterati"
843
"ve refinement for Arnoldi residual");
847
igraphdvout_(&debug_1.logfil, &c__2, xtemp, &debug_1.ndigit, "_naitr: "
848
"iterative refinement ; rnorm and rnorm1 are");
852
/* %-----------------------------------------% */
853
/* | Determine if we need to perform another | */
854
/* | step of re-orthogonalization. | */
855
/* %-----------------------------------------% */
857
if (rnorm1 > *rnorm * .717f) {
859
/* %---------------------------------------% */
860
/* | No need for further refinement. | */
861
/* | The cosine of the angle between the | */
862
/* | corrected residual vector and the old | */
863
/* | residual vector is greater than 0.717 | */
864
/* | In other words the corrected residual | */
865
/* | and the old residual vector share an | */
866
/* | angle of less than arcCOS(0.717) | */
867
/* %---------------------------------------% */
873
/* %-------------------------------------------% */
874
/* | Another step of iterative refinement step | */
875
/* | is required. NITREF is used by stat.h | */
876
/* %-------------------------------------------% */
885
/* %-------------------------------------------------% */
886
/* | Otherwise RESID is numerically in the span of V | */
887
/* %-------------------------------------------------% */
890
for (jj = 1; jj <= i__1; ++jj) {
897
/* %----------------------------------------------% */
898
/* | Branch here directly if iterative refinement | */
899
/* | wasn't necessary or after at most NITER_REF | */
900
/* | steps of iterative refinement. | */
901
/* %----------------------------------------------% */
909
timing_1.titref += t5 - t4;
911
/* %------------------------------------% */
912
/* | STEP 6: Update j = j+1; Continue | */
913
/* %------------------------------------% */
918
timing_1.tnaitr += t1 - t0;
921
for (i__ = max(1,*k); i__ <= i__1; ++i__) {
923
/* %--------------------------------------------% */
924
/* | Check for splitting and deflation. | */
925
/* | Use a standard test as in the QR algorithm | */
926
/* | REFERENCE: LAPACK subroutine dlahqr | */
927
/* %--------------------------------------------% */
929
tst1 = (d__1 = h__[i__ + i__ * h_dim1], abs(d__1)) + (d__2 = h__[
930
i__ + 1 + (i__ + 1) * h_dim1], abs(d__2));
933
tst1 = igraphdlanhs_("1", &i__2, &h__[h_offset], ldh, &workd[*n + 1]);
937
if ((d__1 = h__[i__ + 1 + i__ * h_dim1], abs(d__1)) <= max(d__2,
939
h__[i__ + 1 + i__ * h_dim1] = 0.;
947
igraphdmout_(&debug_1.logfil, &i__1, &i__2, &h__[h_offset], ldh, &
948
debug_1.ndigit, "_naitr: Final upper Hessenberg matrix H"
955
/* %--------------------------------------------------------% */
956
/* | Loop back to extend the factorization by another step. | */
957
/* %--------------------------------------------------------% */
961
/* %---------------------------------------------------------------% */
963
/* | E N D O F M A I N I T E R A T I O N L O O P | */
965
/* %---------------------------------------------------------------% */
970
/* %---------------% */
971
/* | End of igraphdnaitr | */
972
/* %---------------% */
974
} /* igraphdnaitr_ */
added
removed
src/arpack/dnaitr.c
