/* arpack/dsapps.f -- translated by f2c (version 20090411).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#ifdef __cplusplus
extern "C" {
#endif
#include "v3p_netlib.h"

/* Common Block Declarations */

/*Extern struct { */
/*  integer logfil, ndigit, mgetv0, msaupd, msaup2, msaitr, mseigt, msapps, */
/*          msgets, mseupd, mnaupd, mnaup2, mnaitr, mneigh, mnapps, mngets, */
/*          mneupd, mcaupd, mcaup2, mcaitr, mceigh, mcapps, mcgets, mceupd; */
/*} debug_; */

/*#define debug_1 debug_ */

/*Extern struct { */
/*  integer nopx, nbx, nrorth, nitref, nrstrt; */
/*  real tsaupd, tsaup2, tsaitr, tseigt, tsgets, tsapps, tsconv, tnaupd, */
/*          tnaup2, tnaitr, tneigh, tngets, tnapps, tnconv, tcaupd, tcaup2, */
/*          tcaitr, tceigh, tcgets, tcapps, tcconv, tmvopx, tmvbx, tgetv0, */
/*          titref, trvec; */
/*} timing_; */

/*#define timing_1 timing_ */

/* Table of constant values */

static doublereal c_b4 = 0.;
static doublereal c_b5 = 1.;
static doublereal c_b14 = -1.;
static integer c__1 = 1;

/* ----------------------------------------------------------------------- */
/* \BeginDoc */

/* \Name: dsapps */

/* \Description: */
/*  Given the Arnoldi factorization */

/*     A*V_{k} - V_{k}*H_{k} = r_{k+p}*e_{k+p}^T, */

/*  apply NP shifts implicitly resulting in */

/*     A*(V_{k}*Q) - (V_{k}*Q)*(Q^T* H_{k}*Q) = r_{k+p}*e_{k+p}^T * Q */

/*  where Q is an orthogonal matrix of order KEV+NP. Q is the product of */
/*  rotations resulting from the NP bulge chasing sweeps.  The updated Arnoldi */
/*  factorization becomes: */

/*     A*VNEW_{k} - VNEW_{k}*HNEW_{k} = rnew_{k}*e_{k}^T. */

/* \Usage: */
/*  call dsapps */
/*     ( N, KEV, NP, SHIFT, V, LDV, H, LDH, RESID, Q, LDQ, WORKD ) */

/* \Arguments */
/*  N       Integer.  (INPUT) */
/*          Problem size, i.e. dimension of matrix A. */

/*  KEV     Integer.  (INPUT) */
/*          INPUT: KEV+NP is the size of the input matrix H. */
/*          OUTPUT: KEV is the size of the updated matrix HNEW. */

/*  NP      Integer.  (INPUT) */
/*          Number of implicit shifts to be applied. */

/*  SHIFT   Double precision array of length NP.  (INPUT) */
/*          The shifts to be applied. */

/*  V       Double precision N by (KEV+NP) array.  (INPUT/OUTPUT) */
/*          INPUT: V contains the current KEV+NP Arnoldi vectors. */
/*          OUTPUT: VNEW = V(1:n,1:KEV); the updated Arnoldi vectors */
/*          are in the first KEV columns of V. */

/*  LDV     Integer.  (INPUT) */
/*          Leading dimension of V exactly as declared in the calling */
/*          program. */

/*  H       Double precision (KEV+NP) by 2 array.  (INPUT/OUTPUT) */
/*          INPUT: H contains the symmetric tridiagonal matrix of the */
/*          Arnoldi factorization with the subdiagonal in the 1st column */
/*          starting at H(2,1) and the main diagonal in the 2nd column. */
/*          OUTPUT: H contains the updated tridiagonal matrix in the */
/*          KEV leading submatrix. */

/*  LDH     Integer.  (INPUT) */
/*          Leading dimension of H exactly as declared in the calling */
/*          program. */

/*  RESID   Double precision array of length (N).  (INPUT/OUTPUT) */
/*          INPUT: RESID contains the the residual vector r_{k+p}. */
/*          OUTPUT: RESID is the updated residual vector rnew_{k}. */

/*  Q       Double precision KEV+NP by KEV+NP work array.  (WORKSPACE) */
/*          Work array used to accumulate the rotations during the bulge */
/*          chase sweep. */

/*  LDQ     Integer.  (INPUT) */
/*          Leading dimension of Q exactly as declared in the calling */
/*          program. */

/*  WORKD   Double precision work array of length 2*N.  (WORKSPACE) */
/*          Distributed array used in the application of the accumulated */
/*          orthogonal matrix Q. */

/* \EndDoc */

/* ----------------------------------------------------------------------- */

/* \BeginLib */

/* \Local variables: */
/*     xxxxxx  real */

/* \References: */
/*  1. D.C. Sorensen, "Implicit Application of Polynomial Filters in */
/*     a k-Step Arnoldi Method", SIAM J. Matr. Anal. Apps., 13 (1992), */
/*     pp 357-385. */
/*  2. R.B. Lehoucq, "Analysis and Implementation of an Implicitly */
/*     Restarted Arnoldi Iteration", Rice University Technical Report */
/*     TR95-13, Department of Computational and Applied Mathematics. */

/* \Routines called: */
/*     second  ARPACK utility routine for timing. */
/*     dlamch  LAPACK routine that determines machine constants. */
/*     dlartg  LAPACK Givens rotation construction routine. */
/*     dlacpy  LAPACK matrix copy routine. */
/*     dlaset  LAPACK matrix initialization routine. */
/*     dgemv   Level 2 BLAS routine for matrix vector multiplication. */
/*     daxpy   Level 1 BLAS that computes a vector triad. */
/*     dcopy   Level 1 BLAS that copies one vector to another. */
/*     dscal   Level 1 BLAS that scales a vector. */

/* \Author */
/*     Danny Sorensen               Phuong Vu */
/*     Richard Lehoucq              CRPC / Rice University */
/*     Dept. of Computational &     Houston, Texas */
/*     Applied Mathematics */
/*     Rice University */
/*     Houston, Texas */

/* \Revision history: */
/*     12/16/93: Version ' 2.1' */

/* \SCCS Information: @(#) */
/* FILE: sapps.F   SID: 2.5   DATE OF SID: 4/19/96   RELEASE: 2 */

/* \Remarks */
/*  1. In this version, each shift is applied to all the subblocks of */
/*     the tridiagonal matrix H and not just to the submatrix that it */
/*     comes from. This routine assumes that the subdiagonal elements */
/*     of H that are stored in h(1:kev+np,1) are nonegative upon input */
/*     and enforce this condition upon output. This version incorporates */
/*     deflation. See code for documentation. */

/* \EndLib */

/* ----------------------------------------------------------------------- */

/*<        >*/
/* Subroutine */ int 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)
{
    /* Initialized data */

    static logical first = TRUE_;

    /* System generated locals */
    integer h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
            i__3, i__4;
    doublereal d__1, d__2;

    /* Local variables */
    doublereal c__, f, g;
    integer i__, j;
    doublereal r__, s, a1, a2, a3, a4;
/*  static real t0, t1; */
    integer jj;
    doublereal big;
    integer iend, itop;
    extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
            integer *), dgemv_(char *, integer *, integer *, doublereal *,
            doublereal *, integer *, doublereal *, integer *, doublereal *,
            doublereal *, integer *, ftnlen), dcopy_(integer *, doublereal *,
            integer *, doublereal *, integer *), daxpy_(integer *, doublereal
            *, doublereal *, integer *, doublereal *, integer *);
    extern doublereal dlamch_(char *, ftnlen);
    extern /* Subroutine */ int second_(real *);
    static doublereal epsmch;
    integer istart, kplusp /*, msglvl */;
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
            doublereal *, integer *, doublereal *, integer *, ftnlen),
            dlartg_(doublereal *, doublereal *, doublereal *, doublereal *,
            doublereal *), dlaset_(char *, integer *, integer *, doublereal *,
             doublereal *, doublereal *, integer *, ftnlen);


/*     %----------------------------------------------------% */
/*     | Include files for debugging and timing information | */
/*     %----------------------------------------------------% */

/*<       include   'debug.h' >*/
/*<       include   'stat.h' >*/

/* \SCCS Information: @(#) */
/* FILE: debug.h   SID: 2.3   DATE OF SID: 11/16/95   RELEASE: 2 */

/*     %---------------------------------% */
/*     | See debug.doc for documentation | */
/*     %---------------------------------% */
/*<        >*/
/*<       integer    kev, ldh, ldq, ldv, n, np >*/

/*     %------------------% */
/*     | Scalar Arguments | */
/*     %------------------% */

/*     %--------------------------------% */
/*     | See stat.doc for documentation | */
/*     %--------------------------------% */

/* \SCCS Information: @(#) */
/* FILE: stat.h   SID: 2.2   DATE OF SID: 11/16/95   RELEASE: 2 */

/*<       save       t0, t1, t2, t3, t4, t5 >*/

/*<       integer    nopx, nbx, nrorth, nitref, nrstrt >*/
/*<        >*/
/*<        >*/

/*     %-----------------% */
/*     | Array Arguments | */
/*     %-----------------% */

/*<        >*/

/*     %------------% */
/*     | Parameters | */
/*     %------------% */

/*<        >*/
/*<       parameter (one = 1.0D+0, zero = 0.0D+0) >*/

/*     %---------------% */
/*     | Local Scalars | */
/*     %---------------% */

/*<       integer    i, iend, istart, itop, j, jj, kplusp, msglvl >*/
/*<       logical    first >*/
/*<        >*/
/*<       save       epsmch, first >*/


/*     %----------------------% */
/*     | External Subroutines | */
/*     %----------------------% */

/*<        >*/

/*     %--------------------% */
/*     | External Functions | */
/*     %--------------------% */

/*<        >*/
/*<       external   dlamch >*/

/*     %----------------------% */
/*     | Intrinsics Functions | */
/*     %----------------------% */

/*<       intrinsic  abs >*/

/*     %----------------% */
/*     | Data statments | */
/*     %----------------% */

/*<       data       first / .true. / >*/
    /* Parameter adjustments */
    --workd;
    --resid;
    --shift;
    v_dim1 = *ldv;
    v_offset = 1 + v_dim1;
    v -= v_offset;
    h_dim1 = *ldh;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1;
    q -= q_offset;

    /* Function Body */

/*     %-----------------------% */
/*     | Executable Statements | */
/*     %-----------------------% */

/*<       if (first) then >*/
    if (first) {
/*<          epsmch = dlamch('Epsilon-Machine') >*/
        epsmch = dlamch_("Epsilon-Machine", (ftnlen)15);
/*<          first = .false. >*/
        first = FALSE_;
/*<       end if >*/
    }
/*<       itop = 1 >*/
    itop = 1;

/*     %-------------------------------% */
/*     | Initialize timing statistics  | */
/*     | & message level for debugging | */
/*     %-------------------------------% */

/*<       call second (t0) >*/
/*  second_(&t0); */
/*<       msglvl = msapps >*/
/*  msglvl = debug_1.msapps; */

/*<       kplusp = kev + np  >*/
    kplusp = *kev + *np;

/*     %----------------------------------------------% */
/*     | Initialize Q to the identity matrix of order | */
/*     | kplusp used to accumulate the rotations.     | */
/*     %----------------------------------------------% */

/*<       call dlaset ('All', kplusp, kplusp, zero, one, q, ldq) >*/
    dlaset_("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq, (ftnlen)
            3);

/*     %----------------------------------------------% */
/*     | Quick return if there are no shifts to apply | */
/*     %----------------------------------------------% */

/*<       if (np .eq. 0) go to 9000 >*/
    if (*np == 0) {
        goto L9000;
    }

/*     %----------------------------------------------------------% */
/*     | Apply the np shifts implicitly. Apply each shift to the  | */
/*     | whole matrix and not just to the submatrix from which it | */
/*     | comes.                                                   | */
/*     %----------------------------------------------------------% */

/*<       do 90 jj = 1, np >*/
    i__1 = *np;
    for (jj = 1; jj <= i__1; ++jj) {

/*<          istart = itop >*/
        istart = itop;

/*        %----------------------------------------------------------% */
/*        | Check for splitting and deflation. Currently we consider | */
/*        | an off-diagonal element h(i+1,1) negligible if           | */
/*        |         h(i+1,1) .le. epsmch*( |h(i,2)| + |h(i+1,2)| )   | */
/*        | for i=1:KEV+NP-1.                                        | */
/*        | If above condition tests true then we set h(i+1,1) = 0.  | */
/*        | Note that h(1:KEV+NP,1) are assumed to be non negative.  | */
/*        %----------------------------------------------------------% */

/*<    20    continue >*/
L20:

/*        %------------------------------------------------% */
/*        | The following loop exits early if we encounter | */
/*        | a negligible off diagonal element.             | */
/*        %------------------------------------------------% */

/*<          do 30 i = istart, kplusp-1 >*/
        i__2 = kplusp - 1;
        for (i__ = istart; i__ <= i__2; ++i__) {
/*<             big   = abs(h(i,2)) + abs(h(i+1,2)) >*/
            big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[
                    i__ + 1 + (h_dim1 << 1)], abs(d__2));
/*<             if (h(i+1,1) .le. epsmch*big) then >*/
            if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
/*               if (msglvl .gt. 0) then */
/*                  call ivout (logfil, 1, i, ndigit, */
/*     &                 '_sapps: deflation at row/column no.') */
/*                  call ivout (logfil, 1, jj, ndigit, */
/*     &                 '_sapps: occurred before shift number.') */
/*                  call dvout (logfil, 1, h(i+1,1), ndigit, */
/*     &                 '_sapps: the corresponding off diagonal element') */
/*               end if */
/*<                h(i+1,1) = zero >*/
                h__[i__ + 1 + h_dim1] = 0.;
/*<                iend = i >*/
                iend = i__;
/*<                go to 40 >*/
                goto L40;
/*<             end if >*/
            }
/*<    30    continue >*/
/* L30: */
        }
/*<          iend = kplusp >*/
        iend = kplusp;
/*<    40    continue >*/
L40:

/*<          if (istart .lt. iend) then >*/
        if (istart < iend) {

/*           %--------------------------------------------------------% */
/*           | Construct the plane rotation G'(istart,istart+1,theta) | */
/*           | that attempts to drive h(istart+1,1) to zero.          | */
/*           %--------------------------------------------------------% */

/*<              f = h(istart,2) - shift(jj) >*/
            f = h__[istart + (h_dim1 << 1)] - shift[jj];
/*<              g = h(istart+1,1) >*/
            g = h__[istart + 1 + h_dim1];
/*<              call dlartg (f, g, c, s, r) >*/
            dlartg_(&f, &g, &c__, &s, &r__);

/*            %-------------------------------------------------------% */
/*            | Apply rotation to the left and right of H;            | */
/*            | H <- G' * H * G,  where G = G(istart,istart+1,theta). | */
/*            | This will create a "bulge".                           | */
/*            %-------------------------------------------------------% */

/*<              a1 = c*h(istart,2)   + s*h(istart+1,1) >*/
            a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
                    h_dim1];
/*<              a2 = c*h(istart+1,1) + s*h(istart+1,2) >*/
            a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
                    h_dim1 << 1)];
/*<              a4 = c*h(istart+1,2) - s*h(istart+1,1) >*/
            a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
                    h_dim1];
/*<              a3 = c*h(istart+1,1) - s*h(istart,2)  >*/
            a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
                    1)];
/*<              h(istart,2)   = c*a1 + s*a2 >*/
            h__[istart + (h_dim1 << 1)] = c__ * a1 + s * a2;
/*<              h(istart+1,2) = c*a4 - s*a3 >*/
            h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
/*<              h(istart+1,1) = c*a3 + s*a4 >*/
            h__[istart + 1 + h_dim1] = c__ * a3 + s * a4;

/*            %----------------------------------------------------% */
/*            | Accumulate the rotation in the matrix Q;  Q <- Q*G | */
/*            %----------------------------------------------------% */

/*<              do 60 j = 1, min(istart+jj,kplusp) >*/
/* Computing MIN */
            i__3 = istart + jj;
            i__2 = min(i__3,kplusp);
            for (j = 1; j <= i__2; ++j) {
/*<                 a1            =   c*q(j,istart) + s*q(j,istart+1) >*/
                a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
                        q_dim1];
/*<                 q(j,istart+1) = - s*q(j,istart) + c*q(j,istart+1) >*/
                q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
                        c__ * q[j + (istart + 1) * q_dim1];
/*<                 q(j,istart)   = a1 >*/
                q[j + istart * q_dim1] = a1;
/*<    60        continue >*/
/* L60: */
            }


/*            %----------------------------------------------% */
/*            | The following loop chases the bulge created. | */
/*            | Note that the previous rotation may also be  | */
/*            | done within the following loop. But it is    | */
/*            | kept separate to make the distinction among  | */
/*            | the bulge chasing sweeps and the first plane | */
/*            | rotation designed to drive h(istart+1,1) to  | */
/*            | zero.                                        | */
/*            %----------------------------------------------% */

/*<              do 70 i = istart+1, iend-1 >*/
            i__2 = iend - 1;
            for (i__ = istart + 1; i__ <= i__2; ++i__) {

/*               %----------------------------------------------% */
/*               | Construct the plane rotation G'(i,i+1,theta) | */
/*               | that zeros the i-th bulge that was created   | */
/*               | by G(i-1,i,theta). g represents the bulge.   | */
/*               %----------------------------------------------% */

/*<                 f = h(i,1) >*/
                f = h__[i__ + h_dim1];
/*<                 g = s*h(i+1,1) >*/
                g = s * h__[i__ + 1 + h_dim1];

/*               %----------------------------------% */
/*               | Final update with G(i-1,i,theta) | */
/*               %----------------------------------% */

/*<                 h(i+1,1) = c*h(i+1,1) >*/
                h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
/*<                 call dlartg (f, g, c, s, r) >*/
                dlartg_(&f, &g, &c__, &s, &r__);

/*               %-------------------------------------------% */
/*               | The following ensures that h(1:iend-1,1), | */
/*               | the first iend-2 off diagonal of elements | */
/*               | H, remain non negative.                   | */
/*               %-------------------------------------------% */

/*<                 if (r .lt. zero) then >*/
                if (r__ < 0.) {
/*<                    r = -r >*/
                    r__ = -r__;
/*<                    c = -c >*/
                    c__ = -c__;
/*<                    s = -s >*/
                    s = -s;
/*<                 end if >*/
                }

/*               %--------------------------------------------% */
/*               | Apply rotation to the left and right of H; | */
/*               | H <- G * H * G',  where G = G(i,i+1,theta) | */
/*               %--------------------------------------------% */

/*<                 h(i,1) = r >*/
                h__[i__ + h_dim1] = r__;

/*<                 a1 = c*h(i,2)   + s*h(i+1,1) >*/
                a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
                        h_dim1];
/*<                 a2 = c*h(i+1,1) + s*h(i+1,2) >*/
                a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
                        << 1)];
/*<                 a3 = c*h(i+1,1) - s*h(i,2) >*/
                a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
                        ];
/*<                 a4 = c*h(i+1,2) - s*h(i+1,1) >*/
                a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
                        h_dim1];

/*<                 h(i,2)   = c*a1 + s*a2 >*/
                h__[i__ + (h_dim1 << 1)] = c__ * a1 + s * a2;
/*<                 h(i+1,2) = c*a4 - s*a3 >*/
                h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
/*<                 h(i+1,1) = c*a3 + s*a4 >*/
                h__[i__ + 1 + h_dim1] = c__ * a3 + s * a4;

/*               %----------------------------------------------------% */
/*               | Accumulate the rotation in the matrix Q;  Q <- Q*G | */
/*               %----------------------------------------------------% */

/*<                 do 50 j = 1, min( j+jj, kplusp ) >*/
/* Computing MIN */
                i__4 = j + jj;
                i__3 = min(i__4,kplusp);
                for (j = 1; j <= i__3; ++j) {
/*<                    a1       =   c*q(j,i) + s*q(j,i+1) >*/
                    a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
                            q_dim1];
/*<                    q(j,i+1) = - s*q(j,i) + c*q(j,i+1) >*/
                    q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
                            c__ * q[j + (i__ + 1) * q_dim1];
/*<                    q(j,i)   = a1 >*/
                    q[j + i__ * q_dim1] = a1;
/*<    50           continue >*/
/* L50: */
                }

/*<    70        continue >*/
/* L70: */
            }

/*<          end if >*/
        }

/*        %--------------------------% */
/*        | Update the block pointer | */
/*        %--------------------------% */

/*<          istart = iend + 1 >*/
        istart = iend + 1;

/*        %------------------------------------------% */
/*        | Make sure that h(iend,1) is non-negative | */
/*        | If not then set h(iend,1) <-- -h(iend,1) | */
/*        | and negate the last column of Q.         | */
/*        | We have effectively carried out a        | */
/*        | similarity on transformation H           | */
/*        %------------------------------------------% */

/*<          if (h(iend,1) .lt. zero) then >*/
        if (h__[iend + h_dim1] < 0.) {
/*<              h(iend,1) = -h(iend,1) >*/
            h__[iend + h_dim1] = -h__[iend + h_dim1];
/*<              call dscal(kplusp, -one, q(1,iend), 1) >*/
            dscal_(&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
/*<          end if >*/
        }

/*        %--------------------------------------------------------% */
/*        | Apply the same shift to the next block if there is any | */
/*        %--------------------------------------------------------% */

/*<          if (iend .lt. kplusp) go to 20 >*/
        if (iend < kplusp) {
            goto L20;
        }

/*        %-----------------------------------------------------% */
/*        | Check if we can increase the the start of the block | */
/*        %-----------------------------------------------------% */

/*<          do 80 i = itop, kplusp-1 >*/
        i__2 = kplusp - 1;
        for (i__ = itop; i__ <= i__2; ++i__) {
/*<             if (h(i+1,1) .gt. zero) go to 90 >*/
            if (h__[i__ + 1 + h_dim1] > 0.) {
                goto L90;
            }
/*<             itop  = itop + 1 >*/
            ++itop;
/*<    80    continue >*/
/* L80: */
        }

/*        %-----------------------------------% */
/*        | Finished applying the jj-th shift | */
/*        %-----------------------------------% */

/*<    90 continue >*/
L90:
        ;
    }

/*     %------------------------------------------% */
/*     | All shifts have been applied. Check for  | */
/*     | more possible deflation that might occur | */
/*     | after the last shift is applied.         | */
/*     %------------------------------------------% */

/*<       do 100 i = itop, kplusp-1 >*/
    i__1 = kplusp - 1;
    for (i__ = itop; i__ <= i__1; ++i__) {
/*<          big   = abs(h(i,2)) + abs(h(i+1,2)) >*/
        big = (d__1 = h__[i__ + (h_dim1 << 1)], abs(d__1)) + (d__2 = h__[i__
                + 1 + (h_dim1 << 1)], abs(d__2));
/*<          if (h(i+1,1) .le. epsmch*big) then >*/
        if (h__[i__ + 1 + h_dim1] <= epsmch * big) {
/*            if (msglvl .gt. 0) then */
/*               call ivout (logfil, 1, i, ndigit, */
/*     &              '_sapps: deflation at row/column no.') */
/*               call dvout (logfil, 1, h(i+1,1), ndigit, */
/*     &              '_sapps: the corresponding off diagonal element') */
/*            end if */
/*<             h(i+1,1) = zero >*/
            h__[i__ + 1 + h_dim1] = 0.;
/*<          end if >*/
        }
/*<  100  continue >*/
/* L100: */
    }

/*     %-------------------------------------------------% */
/*     | Compute the (kev+1)-st column of (V*Q) and      | */
/*     | temporarily store the result in WORKD(N+1:2*N). | */
/*     | This is not necessary if h(kev+1,1) = 0.         | */
/*     %-------------------------------------------------% */

/*<        >*/
    if (h__[*kev + 1 + h_dim1] > 0.) {
        dgemv_("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
                q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1, (ftnlen)1);
    }

/*     %-------------------------------------------------------% */
/*     | Compute column 1 to kev of (V*Q) in backward order    | */
/*     | taking advantage that Q is an upper triangular matrix | */
/*     | with lower bandwidth np.                              | */
/*     | Place results in v(:,kplusp-kev:kplusp) temporarily.  | */
/*     %-------------------------------------------------------% */

/*<       do 130 i = 1, kev >*/
    i__1 = *kev;
    for (i__ = 1; i__ <= i__1; ++i__) {
/*<        >*/
        i__2 = kplusp - i__ + 1;
        dgemv_("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
                q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1, (ftnlen)1);
/*<          call dcopy (n, workd, 1, v(1,kplusp-i+1), 1) >*/
        dcopy_(n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
                c__1);
/*<   130 continue >*/
/* L130: */
    }

/*     %-------------------------------------------------% */
/*     |  Move v(:,kplusp-kev+1:kplusp) into v(:,1:kev). | */
/*     %-------------------------------------------------% */

/*<       call dlacpy ('All', n, kev, v(1,np+1), ldv, v, ldv) >*/
    dlacpy_("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv,
             (ftnlen)3);

/*     %--------------------------------------------% */
/*     | Copy the (kev+1)-st column of (V*Q) in the | */
/*     | appropriate place if h(kev+1,1) .ne. zero. | */
/*     %--------------------------------------------% */

/*<        >*/
    if (h__[*kev + 1 + h_dim1] > 0.) {
        dcopy_(n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
    }

/*     %-------------------------------------% */
/*     | Update the residual vector:         | */
/*     |    r <- sigmak*r + betak*v(:,kev+1) | */
/*     | where                               | */
/*     |    sigmak = (e_{kev+p}'*Q)*e_{kev}  | */
/*     |    betak = e_{kev+1}'*H*e_{kev}     | */
/*     %-------------------------------------% */

/*<       call dscal (n, q(kplusp,kev), resid, 1) >*/
    dscal_(n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
/*<        >*/
    if (h__[*kev + 1 + h_dim1] > 0.) {
        daxpy_(n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
                 &resid[1], &c__1);
    }

/*      if (msglvl .gt. 1) then */
/*         call dvout (logfil, 1, q(kplusp,kev), ndigit, */
/*     &      '_sapps: sigmak of the updated residual vector') */
/*         call dvout (logfil, 1, h(kev+1,1), ndigit, */
/*     &      '_sapps: betak of the updated residual vector') */
/*         call dvout (logfil, kev, h(1,2), ndigit, */
/*     &      '_sapps: updated main diagonal of H for next iteration') */
/*         if (kev .gt. 1) then */
/*         call dvout (logfil, kev-1, h(2,1), ndigit, */
/*     &      '_sapps: updated sub diagonal of H for next iteration') */
/*         end if */
/*      end if */

/*<       call second (t1) >*/
/*  second_(&t1); */
/*<       tsapps = tsapps + (t1 - t0) >*/
/*  timing_1.tsapps += t1 - t0; */

/*<  9000 continue  >*/
L9000:
/*<       return >*/
    return 0;

/*     %---------------% */
/*     | End of dsapps | */
/*     %---------------% */

/*<       end >*/
} /* dsapps_ */

#ifdef __cplusplus
        }
#endif