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 doublereal c_b4 = 0.;
40
static doublereal c_b5 = 1.;
41
static integer c__1 = 1;
42
static doublereal c_b20 = -1.;
44
/* ----------------------------------------------------------------------- */
50
/* Given the Arnoldi factorization */
52
/* A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T, */
54
/* apply NP shifts implicitly resulting in */
56
/* A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q */
58
/* where Q is an orthogonal matrix of order KEV+NP. Q is the product of */
59
/* rotations resulting from the NP bulge chasing sweeps. The updated Arnoldi */
60
/* factorization becomes: */
62
/* A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T. */
66
/* ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ, WORKD ) */
69
/* N Integer. (INPUT) */
70
/* Problem size, i.e. dimension of matrix A. */
72
/* KEV Integer. (INPUT) */
73
/* INPUT: KEV+NP is the size of the input matrix H. */
74
/* OUTPUT: KEV is the size of the updated matrix HNEW. */
76
/* NP Integer. (INPUT) */
77
/* Number of implicit shifts to be applied. */
79
/* SHIFT Double precision array of length NP. (INPUT) */
80
/* The shifts to be applied. */
82
/* V Double precision N by (KEV+NP) array. (INPUT/OUTPUT) */
83
/* INPUT: V contains the current KEV+NP Arnoldi vectors. */
84
/* OUTPUT: VNEW = V(1:n,1:KEV); the updated Arnoldi vectors */
85
/* are in the first KEV columns of V. */
87
/* LDV Integer. (INPUT) */
88
/* Leading dimension of V exactly as declared in the calling */
91
/* H Double precision (KEV+NP) by 2 array. (INPUT/OUTPUT) */
92
/* INPUT: H contains the symmetric tridiagonal matrix of the */
93
/* Arnoldi factorization with the subdiagonal in the 1st column */
94
/* starting at H(2,1) and the main diagonal in the 2nd column. */
95
/* OUTPUT: H contains the updated tridiagonal matrix in the */
96
/* KEV leading submatrix. */
98
/* LDH Integer. (INPUT) */
99
/* Leading dimension of H exactly as declared in the calling */
102
/* RESID Double precision array of length (N). (INPUT/OUTPUT) */
103
/* INPUT: RESID contains the the residual vector r_{k+p}. */
104
/* OUTPUT: RESID is the updated residual vector rnew_{k}. */
106
/* Q Double precision KEV+NP by KEV+NP work array. (WORKSPACE) */
107
/* Work array used to accumulate the rotations during the bulge */
110
/* LDQ Integer. (INPUT) */
111
/* Leading dimension of Q exactly as declared in the calling */
114
/* WORKD Double precision work array of length 2*N. (WORKSPACE) */
115
/* Distributed array used in the application of the accumulated */
116
/* orthogonal matrix Q. */
120
/* ----------------------------------------------------------------------- */
124
/* \Local variables: */
128
/* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
129
/* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
131
/* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
132
/* Restarted Arnoldi Iteration", Rice University Technical Report */
133
/* TR95-13, Department of Computational and Applied Mathematics. */
135
/* \Routines called: */
136
/* ivout ARPACK utility routine that prints integers. */
137
/* second ARPACK utility routine for timing. */
138
/* dvout ARPACK utility routine that prints vectors. */
139
/* dlamch LAPACK routine that determines machine constants. */
140
/* dlartg LAPACK Givens rotation construction routine. */
141
/* dlacpy LAPACK matrix copy routine. */
142
/* dlaset LAPACK matrix initialization routine. */
143
/* dgemv Level 2 BLAS routine for matrix vector multiplication. */
144
/* daxpy Level 1 BLAS that computes a vector triad. */
145
/* dcopy Level 1 BLAS that copies one vector to another. */
146
/* dscal Level 1 BLAS that scales a vector. */
149
/* Danny Sorensen Phuong Vu */
150
/* Richard Lehoucq CRPC / Rice University */
151
/* Dept. of Computational & Houston, Texas */
152
/* Applied Mathematics */
153
/* Rice University */
156
/* \Revision history: */
157
/* 12/16/93: Version ' 2.4' */
159
/* \SCCS Information: @(#) */
160
/* FILE: sapps.F SID: 2.6 DATE OF SID: 3/28/97 RELEASE: 2 */
163
/* 1. In this version, each shift is applied to all the subblocks of */
164
/* the tridiagonal matrix H and not just to the submatrix that it */
165
/* comes from. This routine assumes that the subdiagonal elements */
166
/* of H that are stored in h(1:kev+np,1) are nonegative upon input */
167
/* and enforce this condition upon output. This version incorporates */
168
/* deflation. See code for documentation. */
172
/* ----------------------------------------------------------------------- */
174
/* Subroutine */ int igraphdsapps_(integer *n, integer *kev, integer *np,
175
doublereal *shift, doublereal *v, integer *ldv, doublereal *h__,
176
integer *ldh, doublereal *resid, doublereal *q, integer *ldq,
179
/* Initialized data */
181
static logical first = TRUE_;
183
/* System generated locals */
184
integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
186
doublereal d__1, d__2;
188
/* Local variables */
189
static doublereal c__, f, g;
190
static integer i__, j;
191
static doublereal r__, s, a1, a2, a3, a4;
194
static doublereal big;
195
extern /* Subroutine */ int igraphdscal_(integer *, doublereal *,
196
doublereal *, integer *), igraphdgemv_(char *, integer *, integer
197
*, doublereal *, doublereal *, integer *, doublereal *, integer *,
198
doublereal *, doublereal *, integer *), igraphdcopy_(
199
integer *, doublereal *, integer *, doublereal *, integer *),
200
igraphdaxpy_(integer *, doublereal *, doublereal *, integer *,
201
doublereal *, integer *), igraphdvout_(integer *, integer *,
202
doublereal *, integer *, char *), igraphivout_(integer *,
203
integer *, integer *, integer *, char *);
204
static integer iend, itop;
205
extern doublereal igraphdlamch_(char *);
206
extern /* Subroutine */ int igraphsecond_(real *), igraphdlacpy_(char *,
207
integer *, integer *, doublereal *, integer *, doublereal *,
208
integer *), igraphdlartg_(doublereal *, doublereal *,
209
doublereal *, doublereal *, doublereal *), igraphdlaset_(char *,
210
integer *, integer *, doublereal *, doublereal *, doublereal *,
212
static doublereal epsmch;
213
static integer istart, kplusp, msglvl;
216
/* %----------------------------------------------------% */
217
/* | Include files for debugging and timing information | */
218
/* %----------------------------------------------------% */
221
/* \SCCS Information: @(#) */
222
/* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
224
/* %---------------------------------% */
225
/* | See debug.doc for documentation | */
226
/* %---------------------------------% */
228
/* %------------------% */
229
/* | Scalar Arguments | */
230
/* %------------------% */
232
/* %--------------------------------% */
233
/* | See stat.doc for documentation | */
234
/* %--------------------------------% */
236
/* \SCCS Information: @(#) */
237
/* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
241
/* %-----------------% */
242
/* | Array Arguments | */
243
/* %-----------------% */
251
/* %---------------% */
252
/* | Local Scalars | */
253
/* %---------------% */
257
/* %----------------------% */
258
/* | External Subroutines | */
259
/* %----------------------% */
262
/* %--------------------% */
263
/* | External Functions | */
264
/* %--------------------% */
267
/* %----------------------% */
268
/* | Intrinsics Functions | */
269
/* %----------------------% */
272
/* %----------------% */
273
/* | Data statments | */
274
/* %----------------% */
276
/* Parameter adjustments */
281
v_offset = 1 + v_dim1;
284
h_offset = 1 + h_dim1;
287
q_offset = 1 + q_dim1;
292
/* %-----------------------% */
293
/* | Executable Statements | */
294
/* %-----------------------% */
297
epsmch = igraphdlamch_("Epsilon-Machine");
302
/* %-------------------------------% */
303
/* | Initialize timing statistics | */
304
/* | & message level for debugging | */
305
/* %-------------------------------% */
308
msglvl = debug_1.msapps;
312
/* %----------------------------------------------% */
313
/* | Initialize Q to the identity matrix of order | */
314
/* | kplusp used to accumulate the rotations. | */
315
/* %----------------------------------------------% */
317
igraphdlaset_("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
319
/* %----------------------------------------------% */
320
/* | Quick return if there are no shifts to apply | */
321
/* %----------------------------------------------% */
327
/* %----------------------------------------------------------% */
328
/* | Apply the np shifts implicitly. Apply each shift to the | */
329
/* | whole matrix and not just to the submatrix from which it | */
331
/* %----------------------------------------------------------% */
334
for (jj = 1; jj <= i__1; ++jj) {
338
/* %----------------------------------------------------------% */
339
/* | Check for splitting and deflation. Currently we consider | */
340
/* | an off-diagonal element h(i+1,1) negligible if | */
341
/* | h(i+1,1) .le. epsmch*( |h(i,2)| + |h(i+1,2)| ) | */
342
/* | for i=1:KEV+NP-1. | */
343
/* | If above condition tests true then we set h(i+1,1) = 0. | */
344
/* | Note that h(1:KEV+NP,1) are assumed to be non negative. | */
345
/* %----------------------------------------------------------% */
349
/* %------------------------------------------------% */
350
/* | The following loop exits early if we encounter | */
351
/* | a negligible off diagonal element. | */
352
/* %------------------------------------------------% */
355
for (i__ = istart; i__ <= i__2; ++i__) {
356
big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[
357
i__ + 1 + (h_dim1 << 1)], abs(d__2));
358
if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
360
igraphivout_(&debug_1.logfil, &c__1, &i__, &
361
debug_1.ndigit, "_sapps: deflation at row/column"
363
igraphivout_(&debug_1.logfil, &c__1, &jj, &debug_1.ndigit,
364
"_sapps: occured before shift number.");
365
igraphdvout_(&debug_1.logfil, &c__1, &h__[i__ + 1 +
366
h_dim1], &debug_1.ndigit, "_sapps: the correspon"
367
"ding off diagonal element");
369
h__[i__ + 1 + h_dim1] = 0.;
380
/* %--------------------------------------------------------% */
381
/* | Construct the plane rotation G'(istart,istart+1,theta) | */
382
/* | that attempts to drive h(istart+1,1) to zero. | */
383
/* %--------------------------------------------------------% */
385
f = h__[istart + (h_dim1 << 1)] - shift[jj];
386
g = h__[istart + 1 + h_dim1];
387
igraphdlartg_(&f, &g, &c__, &s, &r__);
389
/* %-------------------------------------------------------% */
390
/* | Apply rotation to the left and right of H; | */
391
/* | H <- G' * H * G, where G = G(istart,istart+1,theta). | */
392
/* | This will create a "bulge". | */
393
/* %-------------------------------------------------------% */
395
a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
397
a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
399
a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
401
a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
403
h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
404
h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
405
h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
407
/* %----------------------------------------------------% */
408
/* | Accumulate the rotation in the matrix Q; Q <- Q*G | */
409
/* %----------------------------------------------------% */
413
i__2 = min(i__3,kplusp);
414
for (j = 1; j <= i__2; ++j) {
415
a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
417
q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
418
c__ * q[j + (istart + 1) * q_dim1];
419
q[j + istart * q_dim1] = a1;
424
/* %----------------------------------------------% */
425
/* | The following loop chases the bulge created. | */
426
/* | Note that the previous rotation may also be | */
427
/* | done within the following loop. But it is | */
428
/* | kept separate to make the distinction among | */
429
/* | the bulge chasing sweeps and the first plane | */
430
/* | rotation designed to drive h(istart+1,1) to | */
432
/* %----------------------------------------------% */
435
for (i__ = istart + 1; i__ <= i__2; ++i__) {
437
/* %----------------------------------------------% */
438
/* | Construct the plane rotation G'(i,i+1,theta) | */
439
/* | that zeros the i-th bulge that was created | */
440
/* | by G(i-1,i,theta). g represents the bulge. | */
441
/* %----------------------------------------------% */
443
f = h__[i__ + h_dim1];
444
g = s * h__[i__ + 1 + h_dim1];
446
/* %----------------------------------% */
447
/* | Final update with G(i-1,i,theta) | */
448
/* %----------------------------------% */
450
h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
451
igraphdlartg_(&f, &g, &c__, &s, &r__);
453
/* %-------------------------------------------% */
454
/* | The following ensures that h(1:iend-1,1), | */
455
/* | the first iend-2 off diagonal of elements | */
456
/* | H, remain non negative. | */
457
/* %-------------------------------------------% */
465
/* %--------------------------------------------% */
466
/* | Apply rotation to the left and right of H; | */
467
/* | H <- G * H * G', where G = G(i,i+1,theta) | */
468
/* %--------------------------------------------% */
470
h__[i__ + h_dim1] = r__;
472
a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
474
a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
476
a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
478
a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
481
h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
482
h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
483
h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
485
/* %----------------------------------------------------% */
486
/* | Accumulate the rotation in the matrix Q; Q <- Q*G | */
487
/* %----------------------------------------------------% */
491
i__3 = min(i__4,kplusp);
492
for (j = 1; j <= i__3; ++j) {
493
a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
495
q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
496
c__ * q[j + (i__ + 1) * q_dim1];
497
q[j + i__ * q_dim1] = a1;
506
/* %--------------------------% */
507
/* | Update the block pointer | */
508
/* %--------------------------% */
512
/* %------------------------------------------% */
513
/* | Make sure that h(iend,1) is non-negative | */
514
/* | If not then set h(iend,1) <-- -h(iend,1) | */
515
/* | and negate the last column of Q. | */
516
/* | We have effectively carried out a | */
517
/* | similarity on transformation H | */
518
/* %------------------------------------------% */
520
if (h__[iend + h_dim1] < 0.) {
521
h__[iend + h_dim1] = -h__[iend + h_dim1];
522
igraphdscal_(&kplusp, &c_b20, &q[iend * q_dim1 + 1], &c__1);
525
/* %--------------------------------------------------------% */
526
/* | Apply the same shift to the next block if there is any | */
527
/* %--------------------------------------------------------% */
533
/* %-----------------------------------------------------% */
534
/* | Check if we can increase the the start of the block | */
535
/* %-----------------------------------------------------% */
538
for (i__ = itop; i__ <= i__2; ++i__) {
539
if (h__[i__ + 1 + h_dim1] > 0.) {
546
/* %-----------------------------------% */
547
/* | Finished applying the jj-th shift | */
548
/* %-----------------------------------% */
554
/* %------------------------------------------% */
555
/* | All shifts have been applied. Check for | */
556
/* | more possible deflation that might occur | */
557
/* | after the last shift is applied. | */
558
/* %------------------------------------------% */
561
for (i__ = itop; i__ <= i__1; ++i__) {
562
big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[i__
563
+ 1 + (h_dim1 << 1)], abs(d__2));
564
if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
566
igraphivout_(&debug_1.logfil, &c__1, &i__, &debug_1.ndigit,
567
"_sapps: deflation at row/column no.");
568
igraphdvout_(&debug_1.logfil, &c__1, &h__[i__ + 1 + h_dim1], &
569
debug_1.ndigit, "_sapps: the corresponding off diago"
572
h__[i__ + 1 + h_dim1] = 0.;
577
/* %-------------------------------------------------% */
578
/* | Compute the (kev+1)-st column of (V*Q) and | */
579
/* | temporarily store the result in WORKD(N+1:2*N). | */
580
/* | This is not necessary if h(kev+1,1) = 0. | */
581
/* %-------------------------------------------------% */
583
if (h__[*kev + 1 + h_dim1] > 0.) {
584
igraphdgemv_("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1)
585
* q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1)
589
/* %-------------------------------------------------------% */
590
/* | Compute column 1 to kev of (V*Q) in backward order | */
591
/* | taking advantage that Q is an upper triangular matrix | */
592
/* | with lower bandwidth np. | */
593
/* | Place results in v(:,kplusp-kev:kplusp) temporarily. | */
594
/* %-------------------------------------------------------% */
597
for (i__ = 1; i__ <= i__1; ++i__) {
598
i__2 = kplusp - i__ + 1;
599
igraphdgemv_("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__
600
+ 1) * q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1)
602
igraphdcopy_(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1],
607
/* %-------------------------------------------------% */
608
/* | Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). | */
609
/* %-------------------------------------------------% */
611
igraphdlacpy_("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset]
614
/* %--------------------------------------------% */
615
/* | Copy the (kev+1)-st column of (V*Q) in the | */
616
/* | appropriate place if h(kev+1,1) .ne. zero. | */
617
/* %--------------------------------------------% */
619
if (h__[*kev + 1 + h_dim1] > 0.) {
620
igraphdcopy_(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &
624
/* %-------------------------------------% */
625
/* | Update the residual vector: | */
626
/* | r <- sigmak*r + betak*v(:,kev+1) | */
628
/* | sigmak = (e_{kev+p}'*Q)*e_{kev} | */
629
/* | betak = e_{kev+1}'*H*e_{kev} | */
630
/* %-------------------------------------% */
632
igraphdscal_(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
633
if (h__[*kev + 1 + h_dim1] > 0.) {
634
igraphdaxpy_(n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1],
635
&c__1, &resid[1], &c__1);
639
igraphdvout_(&debug_1.logfil, &c__1, &q[kplusp + *kev * q_dim1], &
640
debug_1.ndigit, "_sapps: sigmak of the updated residual vect"
642
igraphdvout_(&debug_1.logfil, &c__1, &h__[*kev + 1 + h_dim1], &
643
debug_1.ndigit, "_sapps: betak of the updated residual vector"
645
igraphdvout_(&debug_1.logfil, kev, &h__[(h_dim1 << 1) + 1], &
646
debug_1.ndigit, "_sapps: updated main diagonal of H for next"
650
igraphdvout_(&debug_1.logfil, &i__1, &h__[h_dim1 + 2], &
651
debug_1.ndigit, "_sapps: updated sub diagonal of H for n"
657
timing_1.tsapps += t1 - t0;
662
/* %---------------% */
663
/* | End of dsapps | */
664
/* %---------------% */
666
} /* igraphdsapps_ */
added
removed
src/arpack/dsapps.c
