1 /* arpack/dsapps.f -- translated by f2c (version 20090411).
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
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Common Block Declarations */
19
20 /*Extern struct { */
21 /* integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps, */
22 /* msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, */
23 /* mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd; */
24 /*} debug_; */
25
26 /*#define debug_1 debug_ */
27
28 /*Extern struct { */
29 /* integer nopx, nbx, nrorth, nitref, nrstrt; */
30 /* real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, */
31 /* tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2, */
32 /* tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0, */
33 /* titref, trvec; */
34 /*} timing_; */
35
36 /*#define timing_1 timing_ */
37
38 /* Table of constant values */
39
40 static doublereal c_b4 = 0.;
41 static doublereal c_b5 = 1.;
42 static doublereal c_b14 = -1.;
43 static integer c__1 = 1;
44
45 /* ----------------------------------------------------------------------- */
46 /* \BeginDoc */
47
48 /* \Name: dsapps */
49
50 /* \Description: */
51 /* Given the Arnoldi factorization */
52
53 /* A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T, */
54
55 /* apply NP shifts implicitly resulting in */
56
57 /* A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q */
58
59 /* where Q is an orthogonal matrix of order KEV+NP. Q is the product of */
60 /* rotations resulting from the NP bulge chasing sweeps. The updated Arnoldi */
61 /* factorization becomes: */
62
63 /* A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T. */
64
65 /* \Usage: */
66 /* call dsapps */
67 /* ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ, WORKD ) */
68
69 /* \Arguments */
70 /* N Integer. (INPUT) */
71 /* Problem size, i.e. dimension of matrix A. */
72
73 /* KEV Integer. (INPUT) */
74 /* INPUT: KEV+NP is the size of the input matrix H. */
75 /* OUTPUT: KEV is the size of the updated matrix HNEW. */
76
77 /* NP Integer. (INPUT) */
78 /* Number of implicit shifts to be applied. */
79
80 /* SHIFT Double precision array of length NP. (INPUT) */
81 /* The shifts to be applied. */
82
83 /* V Double precision N by (KEV+NP) array. (INPUT/OUTPUT) */
84 /* INPUT: V contains the current KEV+NP Arnoldi vectors. */
85 /* OUTPUT: VNEW = V(1:n,1:KEV); the updated Arnoldi vectors */
86 /* are in the first KEV columns of V. */
87
88 /* LDV Integer. (INPUT) */
89 /* Leading dimension of V exactly as declared in the calling */
90 /* program. */
91
92 /* H Double precision (KEV+NP) by 2 array. (INPUT/OUTPUT) */
93 /* INPUT: H contains the symmetric tridiagonal matrix of the */
94 /* Arnoldi factorization with the subdiagonal in the 1st column */
95 /* starting at H(2,1) and the main diagonal in the 2nd column. */
96 /* OUTPUT: H contains the updated tridiagonal matrix in the */
97 /* KEV leading submatrix. */
98
99 /* LDH Integer. (INPUT) */
100 /* Leading dimension of H exactly as declared in the calling */
101 /* program. */
102
103 /* RESID Double precision array of length (N). (INPUT/OUTPUT) */
104 /* INPUT: RESID contains the the residual vector r_{k+p}. */
105 /* OUTPUT: RESID is the updated residual vector rnew_{k}. */
106
107 /* Q Double precision KEV+NP by KEV+NP work array. (WORKSPACE) */
108 /* Work array used to accumulate the rotations during the bulge */
109 /* chase sweep. */
110
111 /* LDQ Integer. (INPUT) */
112 /* Leading dimension of Q exactly as declared in the calling */
113 /* program. */
114
115 /* WORKD Double precision work array of length 2*N. (WORKSPACE) */
116 /* Distributed array used in the application of the accumulated */
117 /* orthogonal matrix Q. */
118
119 /* \EndDoc */
120
121 /* ----------------------------------------------------------------------- */
122
123 /* \BeginLib */
124
125 /* \Local variables: */
126 /* xxxxxx real */
127
128 /* \References: */
129 /* 1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
130 /* a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
131 /* pp 357-385. */
132 /* 2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
133 /* Restarted Arnoldi Iteration", Rice University Technical Report */
134 /* TR95-13, Department of Computational and Applied Mathematics. */
135
136 /* \Routines called: */
137 /* second ARPACK utility routine for timing. */
138 /* dlamch LAPACK routine that determines machine constants. */
139 /* dlartg LAPACK Givens rotation construction routine. */
140 /* dlacpy LAPACK matrix copy routine. */
141 /* dlaset LAPACK matrix initialization routine. */
142 /* dgemv Level 2 BLAS routine for matrix vector multiplication. */
143 /* daxpy Level 1 BLAS that computes a vector triad. */
144 /* dcopy Level 1 BLAS that copies one vector to another. */
145 /* dscal Level 1 BLAS that scales a vector. */
146
147 /* \Author */
148 /* Danny Sorensen Phuong Vu */
149 /* Richard Lehoucq CRPC / Rice University */
150 /* Dept. of Computational & Houston, Texas */
151 /* Applied Mathematics */
152 /* Rice University */
153 /* Houston, Texas */
154
155 /* \Revision history: */
156 /* 12/16/93: Version ' 2.1' */
157
158 /* \SCCS Information: @(#) */
159 /* FILE: sapps.F SID: 2.5 DATE OF SID: 4/19/96 RELEASE: 2 */
160
161 /* \Remarks */
162 /* 1. In this version, each shift is applied to all the subblocks of */
163 /* the tridiagonal matrix H and not just to the submatrix that it */
164 /* comes from. This routine assumes that the subdiagonal elements */
165 /* of H that are stored in h(1:kev+np,1) are nonegative upon input */
166 /* and enforce this condition upon output. This version incorporates */
167 /* deflation. See code for documentation. */
168
169 /* \EndLib */
170
171 /* ----------------------------------------------------------------------- */
172
173 /*< >*/
dsapps_(integer * n,integer * kev,integer * np,doublereal * shift,doublereal * v,integer * ldv,doublereal * h__,integer * ldh,doublereal * resid,doublereal * q,integer * ldq,doublereal * workd)174 /* Subroutine */ int dsapps_(integer *n, integer *kev, integer *np,
175 doublereal *shift, doublereal *v, integer *ldv, doublereal *h__,
176 integer *ldh, doublereal *resid, doublereal *q, integer *ldq,
177 doublereal *workd)
178 {
179 /* Initialized data */
180
181 static logical first = TRUE_;
182
183 /* System generated locals */
184 integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
185 i__3, i__4;
186 doublereal d__1, d__2;
187
188 /* Local variables */
189 doublereal c__, f, g;
190 integer i__, j;
191 doublereal r__, s, a1, a2, a3, a4;
192 /* static real t0, t1; */
193 integer jj;
194 doublereal big;
195 integer iend, itop;
196 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
197 integer *), dgemv_(char *, integer *, integer *, doublereal *,
198 doublereal *, integer *, doublereal *, integer *, doublereal *,
199 doublereal *, integer *, ftnlen), dcopy_(integer *, doublereal *,
200 integer *, doublereal *, integer *), daxpy_(integer *, doublereal
201 *, doublereal *, integer *, doublereal *, integer *);
202 extern doublereal dlamch_(char *, ftnlen);
203 extern /* Subroutine */ int second_(real *);
204 static doublereal epsmch;
205 integer istart, kplusp /*, msglvl */;
206 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
207 doublereal *, integer *, doublereal *, integer *, ftnlen),
208 dlartg_(doublereal *, doublereal *, doublereal *, doublereal *,
209 doublereal *), dlaset_(char *, integer *, integer *, doublereal *,
210 doublereal *, doublereal *, integer *, ftnlen);
211
212
213 /* %----------------------------------------------------% */
214 /* | Include files for debugging and timing information | */
215 /* %----------------------------------------------------% */
216
217 /*< include 'debug.h' >*/
218 /*< include 'stat.h' >*/
219
220 /* \SCCS Information: @(#) */
221 /* FILE: debug.h SID: 2.3 DATE OF SID: 11/16/95 RELEASE: 2 */
222
223 /* %---------------------------------% */
224 /* | See debug.doc for documentation | */
225 /* %---------------------------------% */
226 /*< >*/
227 /*< integer kev, ldh, ldq, ldv, n, np >*/
228
229 /* %------------------% */
230 /* | Scalar Arguments | */
231 /* %------------------% */
232
233 /* %--------------------------------% */
234 /* | See stat.doc for documentation | */
235 /* %--------------------------------% */
236
237 /* \SCCS Information: @(#) */
238 /* FILE: stat.h SID: 2.2 DATE OF SID: 11/16/95 RELEASE: 2 */
239
240 /*< save t0, t1, t2, t3, t4, t5 >*/
241
242 /*< integer nopx, nbx, nrorth, nitref, nrstrt >*/
243 /*< >*/
244 /*< >*/
245
246 /* %-----------------% */
247 /* | Array Arguments | */
248 /* %-----------------% */
249
250 /*< >*/
251
252 /* %------------% */
253 /* | Parameters | */
254 /* %------------% */
255
256 /*< >*/
257 /*< parameter (one = 1.0D+0, zero = 0.0D+0) >*/
258
259 /* %---------------% */
260 /* | Local Scalars | */
261 /* %---------------% */
262
263 /*< integer i, iend, istart, itop, j, jj, kplusp, msglvl >*/
264 /*< logical first >*/
265 /*< >*/
266 /*< save epsmch, first >*/
267
268
269 /* %----------------------% */
270 /* | External Subroutines | */
271 /* %----------------------% */
272
273 /*< >*/
274
275 /* %--------------------% */
276 /* | External Functions | */
277 /* %--------------------% */
278
279 /*< >*/
280 /*< external dlamch >*/
281
282 /* %----------------------% */
283 /* | Intrinsics Functions | */
284 /* %----------------------% */
285
286 /*< intrinsic abs >*/
287
288 /* %----------------% */
289 /* | Data statments | */
290 /* %----------------% */
291
292 /*< data first / .true. / >*/
293 /* Parameter adjustments */
294 --workd;
295 --resid;
296 --shift;
297 v_dim1 = *ldv;
298 v_offset = 1 + v_dim1;
299 v -= v_offset;
300 h_dim1 = *ldh;
301 h_offset = 1 + h_dim1;
302 h__ -= h_offset;
303 q_dim1 = *ldq;
304 q_offset = 1 + q_dim1;
305 q -= q_offset;
306
307 /* Function Body */
308
309 /* %-----------------------% */
310 /* | Executable Statements | */
311 /* %-----------------------% */
312
313 /*< if (first) then >*/
314 if (first) {
315 /*< epsmch = dlamch('Epsilon-Machine') >*/
316 epsmch = dlamch_("Epsilon-Machine", (ftnlen)15);
317 /*< first = .false. >*/
318 first = FALSE_;
319 /*< end if >*/
320 }
321 /*< itop = 1 >*/
322 itop = 1;
323
324 /* %-------------------------------% */
325 /* | Initialize timing statistics | */
326 /* | & message level for debugging | */
327 /* %-------------------------------% */
328
329 /*< call second (t0) >*/
330 /* second_(&t0); */
331 /*< msglvl = msapps >*/
332 /* msglvl = debug_1.msapps; */
333
334 /*< kplusp = kev + np >*/
335 kplusp = *kev + *np;
336
337 /* %----------------------------------------------% */
338 /* | Initialize Q to the identity matrix of order | */
339 /* | kplusp used to accumulate the rotations. | */
340 /* %----------------------------------------------% */
341
342 /*< call dlaset ('All', kplusp, kplusp, zero, one, q, ldq) >*/
343 dlaset_("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq, (ftnlen)
344 3);
345
346 /* %----------------------------------------------% */
347 /* | Quick return if there are no shifts to apply | */
348 /* %----------------------------------------------% */
349
350 /*< if (np .eq. 0) go to 9000 >*/
351 if (*np == 0) {
352 goto L9000;
353 }
354
355 /* %----------------------------------------------------------% */
356 /* | Apply the np shifts implicitly. Apply each shift to the | */
357 /* | whole matrix and not just to the submatrix from which it | */
358 /* | comes. | */
359 /* %----------------------------------------------------------% */
360
361 /*< do 90 jj = 1, np >*/
362 i__1 = *np;
363 for (jj = 1; jj <= i__1; ++jj) {
364
365 /*< istart = itop >*/
366 istart = itop;
367
368 /* %----------------------------------------------------------% */
369 /* | Check for splitting and deflation. Currently we consider | */
370 /* | an off-diagonal element h(i+1,1) negligible if | */
371 /* | h(i+1,1) .le. epsmch*( |h(i,2)| + |h(i+1,2)| ) | */
372 /* | for i=1:KEV+NP-1. | */
373 /* | If above condition tests true then we set h(i+1,1) = 0. | */
374 /* | Note that h(1:KEV+NP,1) are assumed to be non negative. | */
375 /* %----------------------------------------------------------% */
376
377 /*< 20 continue >*/
378 L20:
379
380 /* %------------------------------------------------% */
381 /* | The following loop exits early if we encounter | */
382 /* | a negligible off diagonal element. | */
383 /* %------------------------------------------------% */
384
385 /*< do 30 i = istart, kplusp-1 >*/
386 i__2 = kplusp - 1;
387 for (i__ = istart; i__ <= i__2; ++i__) {
388 /*< big = abs(h(i,2)) + abs(h(i+1,2)) >*/
389 big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[
390 i__ + 1 + (h_dim1 << 1)], abs(d__2));
391 /*< if (h(i+1,1) .le. epsmch*big) then >*/
392 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
393 /* if (msglvl .gt. 0) then */
394 /* call ivout (logfil, 1, i, ndigit, */
395 /* & '_sapps: deflation at row/column no.') */
396 /* call ivout (logfil, 1, jj, ndigit, */
397 /* & '_sapps: occurred before shift number.') */
398 /* call dvout (logfil, 1, h(i+1,1), ndigit, */
399 /* & '_sapps: the corresponding off diagonal element') */
400 /* end if */
401 /*< h(i+1,1) = zero >*/
402 h__[i__ + 1 + h_dim1] = 0.;
403 /*< iend = i >*/
404 iend = i__;
405 /*< go to 40 >*/
406 goto L40;
407 /*< end if >*/
408 }
409 /*< 30 continue >*/
410 /* L30: */
411 }
412 /*< iend = kplusp >*/
413 iend = kplusp;
414 /*< 40 continue >*/
415 L40:
416
417 /*< if (istart .lt. iend) then >*/
418 if (istart < iend) {
419
420 /* %--------------------------------------------------------% */
421 /* | Construct the plane rotation G'(istart,istart+1,theta) | */
422 /* | that attempts to drive h(istart+1,1) to zero. | */
423 /* %--------------------------------------------------------% */
424
425 /*< f = h(istart,2) - shift(jj) >*/
426 f = h__[istart + (h_dim1 << 1)] - shift[jj];
427 /*< g = h(istart+1,1) >*/
428 g = h__[istart + 1 + h_dim1];
429 /*< call dlartg (f, g, c, s, r) >*/
430 dlartg_(&f, &g, &c__, &s, &r__);
431
432 /* %-------------------------------------------------------% */
433 /* | Apply rotation to the left and right of H; | */
434 /* | H <- G' * H * G, where G = G(istart,istart+1,theta). | */
435 /* | This will create a "bulge". | */
436 /* %-------------------------------------------------------% */
437
438 /*< a1 = c*h(istart,2) + s*h(istart+1,1) >*/
439 a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
440 h_dim1];
441 /*< a2 = c*h(istart+1,1) + s*h(istart+1,2) >*/
442 a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
443 h_dim1 << 1)];
444 /*< a4 = c*h(istart+1,2) - s*h(istart+1,1) >*/
445 a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
446 h_dim1];
447 /*< a3 = c*h(istart+1,1) - s*h(istart,2) >*/
448 a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
449 1)];
450 /*< h(istart,2) = c*a1 + s*a2 >*/
451 h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
452 /*< h(istart+1,2) = c*a4 - s*a3 >*/
453 h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
454 /*< h(istart+1,1) = c*a3 + s*a4 >*/
455 h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;
456
457 /* %----------------------------------------------------% */
458 /* | Accumulate the rotation in the matrix Q; Q <- Q*G | */
459 /* %----------------------------------------------------% */
460
461 /*< do 60 j = 1, min(istart+jj,kplusp) >*/
462 /* Computing MIN */
463 i__3 = istart + jj;
464 i__2 = min(i__3,kplusp);
465 for (j = 1; j <= i__2; ++j) {
466 /*< a1 = c*q(j,istart) + s*q(j,istart+1) >*/
467 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
468 q_dim1];
469 /*< q(j,istart+1) = - s*q(j,istart) + c*q(j,istart+1) >*/
470 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
471 c__ * q[j + (istart + 1) * q_dim1];
472 /*< q(j,istart) = a1 >*/
473 q[j + istart * q_dim1] = a1;
474 /*< 60 continue >*/
475 /* L60: */
476 }
477
478
479 /* %----------------------------------------------% */
480 /* | The following loop chases the bulge created. | */
481 /* | Note that the previous rotation may also be | */
482 /* | done within the following loop. But it is | */
483 /* | kept separate to make the distinction among | */
484 /* | the bulge chasing sweeps and the first plane | */
485 /* | rotation designed to drive h(istart+1,1) to | */
486 /* | zero. | */
487 /* %----------------------------------------------% */
488
489 /*< do 70 i = istart+1, iend-1 >*/
490 i__2 = iend - 1;
491 for (i__ = istart + 1; i__ <= i__2; ++i__) {
492
493 /* %----------------------------------------------% */
494 /* | Construct the plane rotation G'(i,i+1,theta) | */
495 /* | that zeros the i-th bulge that was created | */
496 /* | by G(i-1,i,theta). g represents the bulge. | */
497 /* %----------------------------------------------% */
498
499 /*< f = h(i,1) >*/
500 f = h__[i__ + h_dim1];
501 /*< g = s*h(i+1,1) >*/
502 g = s * h__[i__ + 1 + h_dim1];
503
504 /* %----------------------------------% */
505 /* | Final update with G(i-1,i,theta) | */
506 /* %----------------------------------% */
507
508 /*< h(i+1,1) = c*h(i+1,1) >*/
509 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
510 /*< call dlartg (f, g, c, s, r) >*/
511 dlartg_(&f, &g, &c__, &s, &r__);
512
513 /* %-------------------------------------------% */
514 /* | The following ensures that h(1:iend-1,1), | */
515 /* | the first iend-2 off diagonal of elements | */
516 /* | H, remain non negative. | */
517 /* %-------------------------------------------% */
518
519 /*< if (r .lt. zero) then >*/
520 if (r__ < 0.) {
521 /*< r = -r >*/
522 r__ = -r__;
523 /*< c = -c >*/
524 c__ = -c__;
525 /*< s = -s >*/
526 s = -s;
527 /*< end if >*/
528 }
529
530 /* %--------------------------------------------% */
531 /* | Apply rotation to the left and right of H; | */
532 /* | H <- G * H * G', where G = G(i,i+1,theta) | */
533 /* %--------------------------------------------% */
534
535 /*< h(i,1) = r >*/
536 h__[i__ + h_dim1] = r__;
537
538 /*< a1 = c*h(i,2) + s*h(i+1,1) >*/
539 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
540 h_dim1];
541 /*< a2 = c*h(i+1,1) + s*h(i+1,2) >*/
542 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
543 << 1)];
544 /*< a3 = c*h(i+1,1) - s*h(i,2) >*/
545 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
546 ];
547 /*< a4 = c*h(i+1,2) - s*h(i+1,1) >*/
548 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
549 h_dim1];
550
551 /*< h(i,2) = c*a1 + s*a2 >*/
552 h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
553 /*< h(i+1,2) = c*a4 - s*a3 >*/
554 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
555 /*< h(i+1,1) = c*a3 + s*a4 >*/
556 h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;
557
558 /* %----------------------------------------------------% */
559 /* | Accumulate the rotation in the matrix Q; Q <- Q*G | */
560 /* %----------------------------------------------------% */
561
562 /*< do 50 j = 1, min( j+jj, kplusp ) >*/
563 /* Computing MIN */
564 i__4 = j + jj;
565 i__3 = min(i__4,kplusp);
566 for (j = 1; j <= i__3; ++j) {
567 /*< a1 = c*q(j,i) + s*q(j,i+1) >*/
568 a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
569 q_dim1];
570 /*< q(j,i+1) = - s*q(j,i) + c*q(j,i+1) >*/
571 q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
572 c__ * q[j + (i__ + 1) * q_dim1];
573 /*< q(j,i) = a1 >*/
574 q[j + i__ * q_dim1] = a1;
575 /*< 50 continue >*/
576 /* L50: */
577 }
578
579 /*< 70 continue >*/
580 /* L70: */
581 }
582
583 /*< end if >*/
584 }
585
586 /* %--------------------------% */
587 /* | Update the block pointer | */
588 /* %--------------------------% */
589
590 /*< istart = iend + 1 >*/
591 istart = iend + 1;
592
593 /* %------------------------------------------% */
594 /* | Make sure that h(iend,1) is non-negative | */
595 /* | If not then set h(iend,1) <-- -h(iend,1) | */
596 /* | and negate the last column of Q. | */
597 /* | We have effectively carried out a | */
598 /* | similarity on transformation H | */
599 /* %------------------------------------------% */
600
601 /*< if (h(iend,1) .lt. zero) then >*/
602 if (h__[iend + h_dim1] < 0.) {
603 /*< h(iend,1) = -h(iend,1) >*/
604 h__[iend + h_dim1] = -h__[iend + h_dim1];
605 /*< call dscal(kplusp, -one, q(1,iend), 1) >*/
606 dscal_(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
607 /*< end if >*/
608 }
609
610 /* %--------------------------------------------------------% */
611 /* | Apply the same shift to the next block if there is any | */
612 /* %--------------------------------------------------------% */
613
614 /*< if (iend .lt. kplusp) go to 20 >*/
615 if (iend < kplusp) {
616 goto L20;
617 }
618
619 /* %-----------------------------------------------------% */
620 /* | Check if we can increase the the start of the block | */
621 /* %-----------------------------------------------------% */
622
623 /*< do 80 i = itop, kplusp-1 >*/
624 i__2 = kplusp - 1;
625 for (i__ = itop; i__ <= i__2; ++i__) {
626 /*< if (h(i+1,1) .gt. zero) go to 90 >*/
627 if (h__[i__ + 1 + h_dim1] > 0.) {
628 goto L90;
629 }
630 /*< itop = itop + 1 >*/
631 ++itop;
632 /*< 80 continue >*/
633 /* L80: */
634 }
635
636 /* %-----------------------------------% */
637 /* | Finished applying the jj-th shift | */
638 /* %-----------------------------------% */
639
640 /*< 90 continue >*/
641 L90:
642 ;
643 }
644
645 /* %------------------------------------------% */
646 /* | All shifts have been applied. Check for | */
647 /* | more possible deflation that might occur | */
648 /* | after the last shift is applied. | */
649 /* %------------------------------------------% */
650
651 /*< do 100 i = itop, kplusp-1 >*/
652 i__1 = kplusp - 1;
653 for (i__ = itop; i__ <= i__1; ++i__) {
654 /*< big = abs(h(i,2)) + abs(h(i+1,2)) >*/
655 big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[i__
656 + 1 + (h_dim1 << 1)], abs(d__2));
657 /*< if (h(i+1,1) .le. epsmch*big) then >*/
658 if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
659 /* if (msglvl .gt. 0) then */
660 /* call ivout (logfil, 1, i, ndigit, */
661 /* & '_sapps: deflation at row/column no.') */
662 /* call dvout (logfil, 1, h(i+1,1), ndigit, */
663 /* & '_sapps: the corresponding off diagonal element') */
664 /* end if */
665 /*< h(i+1,1) = zero >*/
666 h__[i__ + 1 + h_dim1] = 0.;
667 /*< end if >*/
668 }
669 /*< 100 continue >*/
670 /* L100: */
671 }
672
673 /* %-------------------------------------------------% */
674 /* | Compute the (kev+1)-st column of (V*Q) and | */
675 /* | temporarily store the result in WORKD(N+1:2*N). | */
676 /* | This is not necessary if h(kev+1,1) = 0. | */
677 /* %-------------------------------------------------% */
678
679 /*< >*/
680 if (h__[*kev + 1 + h_dim1] > 0.) {
681 dgemv_("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
682 q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1, (ftnlen)1);
683 }
684
685 /* %-------------------------------------------------------% */
686 /* | Compute column 1 to kev of (V*Q) in backward order | */
687 /* | taking advantage that Q is an upper triangular matrix | */
688 /* | with lower bandwidth np. | */
689 /* | Place results in v(:,kplusp-kev:kplusp) temporarily. | */
690 /* %-------------------------------------------------------% */
691
692 /*< do 130 i = 1, kev >*/
693 i__1 = *kev;
694 for (i__ = 1; i__ <= i__1; ++i__) {
695 /*< >*/
696 i__2 = kplusp - i__ + 1;
697 dgemv_("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
698 q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1, (ftnlen)1);
699 /*< call dcopy (n, workd, 1, v(1,kplusp-i+1), 1) >*/
700 dcopy_(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
701 c__1);
702 /*< 130 continue >*/
703 /* L130: */
704 }
705
706 /* %-------------------------------------------------% */
707 /* | Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). | */
708 /* %-------------------------------------------------% */
709
710 /*< call dlacpy ('All', n, kev, v(1,np+1), ldv, v, ldv) >*/
711 dlacpy_("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv,
712 (ftnlen)3);
713
714 /* %--------------------------------------------% */
715 /* | Copy the (kev+1)-st column of (V*Q) in the | */
716 /* | appropriate place if h(kev+1,1) .ne. zero. | */
717 /* %--------------------------------------------% */
718
719 /*< >*/
720 if (h__[*kev + 1 + h_dim1] > 0.) {
721 dcopy_(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
722 }
723
724 /* %-------------------------------------% */
725 /* | Update the residual vector: | */
726 /* | r <- sigmak*r + betak*v(:,kev+1) | */
727 /* | where | */
728 /* | sigmak = (e_{kev+p}'*Q)*e_{kev} | */
729 /* | betak = e_{kev+1}'*H*e_{kev} | */
730 /* %-------------------------------------% */
731
732 /*< call dscal (n, q(kplusp,kev), resid, 1) >*/
733 dscal_(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
734 /*< >*/
735 if (h__[*kev + 1 + h_dim1] > 0.) {
736 daxpy_(n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
737 &resid[1], &c__1);
738 }
739
740 /* if (msglvl .gt. 1) then */
741 /* call dvout (logfil, 1, q(kplusp,kev), ndigit, */
742 /* & '_sapps: sigmak of the updated residual vector') */
743 /* call dvout (logfil, 1, h(kev+1,1), ndigit, */
744 /* & '_sapps: betak of the updated residual vector') */
745 /* call dvout (logfil, kev, h(1,2), ndigit, */
746 /* & '_sapps: updated main diagonal of H for next iteration') */
747 /* if (kev .gt. 1) then */
748 /* call dvout (logfil, kev-1, h(2,1), ndigit, */
749 /* & '_sapps: updated sub diagonal of H for next iteration') */
750 /* end if */
751 /* end if */
752
753 /*< call second (t1) >*/
754 /* second_(&t1); */
755 /*< tsapps = tsapps + (t1 - t0) >*/
756 /* timing_1.tsapps += t1 - t0; */
757
758 /*< 9000 continue >*/
759 L9000:
760 /*< return >*/
761 return 0;
762
763 /* %---------------% */
764 /* | End of dsapps | */
765 /* %---------------% */
766
767 /*< end >*/
768 } /* dsapps_ */
769
770 #ifdef __cplusplus
771 }
772 #endif
773