fsdavidson.cpp

An object-oriented C++ implementation of Davidson method for finding a few selected extreme eigenpai

/*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*
*  Purpose
*  =======
*
*  DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
*     [  A   B  ]
*     [  B   C  ].
*  On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
*  eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
*  eigenvector for RT1, giving the decomposition
*
*     [ CS1  SN1 ] [  A   B  ] [ CS1 -SN1 ]  =  [ RT1  0  ]
*     [-SN1  CS1 ] [  B   C  ] [ SN1  CS1 ]     [  0  RT2 ].
*
*  Arguments
*  =========
*
*  A       (input) DOUBLE PRECISION
*          The (1,1) element of the 2-by-2 matrix.
*
*  B       (input) DOUBLE PRECISION
*          The (1,2) element and the conjugate of the (2,1) element of
*          the 2-by-2 matrix.
*
*  C       (input) DOUBLE PRECISION
*          The (2,2) element of the 2-by-2 matrix.
*
*  RT1     (output) DOUBLE PRECISION
*          The eigenvalue of larger absolute value.
*
*  RT2     (output) DOUBLE PRECISION
*          The eigenvalue of smaller absolute value.
*
*  CS1     (output) DOUBLE PRECISION
*  SN1     (output) DOUBLE PRECISION
*          The vector (CS1, SN1) is a unit right eigenvector for RT1.
*
*  Further Details
*  ===============
*
*  RT1 is accurate to a few ulps barring over/underflow.
*
*  RT2 may be inaccurate if there is massive cancellation in the
*  determinant A*C-B*B; higher precision or correctly rounded or
*  correctly truncated arithmetic would be needed to compute RT2
*  accurately in all cases.
*
*  CS1 and SN1 are accurate to a few ulps barring over/underflow.
*
*  Overflow is possible only if RT1 is within a factor of 5 of overflow.
*  Underflow is harmless if the input data is 0 or exceeds
*     underflow_threshold / macheps.
*/
{
    ftyp acmn=a, acmx=c, acs, cs, ct, rt, tn,
         sm = a + c, df = a - c, adf = fabs( df ), tb = b + b, ab = fabs( tb );
    ntyp sgn1, sgn2;

    //Compute the eigenvalues
    if (fabs(a)>fabs(c)) {
        acmx = a;
        acmn = c;
    }
    if (adf>ab) rt = adf*sqrt( 1.0+( ab / adf )*( ab / adf ));
    	else if (adf<ab) rt = ab*sqrt( 1.0+( adf / ab )*( adf / ab ));
        	else
            	//Includes case AB=ADF=0
            	rt = ab*sqrt( 2.0 );
    if (sm<0.0) {
        *rt1 = 0.5*( sm-rt );
        sgn1 = -1;
        /*Order of execution important.
          To get fully accurate smaller eigenvalue,
          next line needs to be executed in higher precision.  */
        *rt2 = ( acmx / *rt1 )*acmn - ( b / *rt1 )*b;
    } else if (sm>0.0) {
        *rt1 = 0.5*( sm+rt );
        sgn1 = 1;
        /*Order of execution important.
          To get fully accurate smaller eigenvalue,
          next line needs to be executed in higher precision.  */
        *rt2 = ( acmx / *rt1 )*acmn - ( b / *rt1 )*b;
    } else {
    	//Includes case RT1 = RT2 = 0
        *rt1 = 0.5*rt;
        *rt2 = -0.5*rt;
        sgn1 = 1;
    }
    //Compute the eigenvector
	if (df>=0.0) {
        cs = df + rt;
        sgn2 = 1;
    } else {
        cs = df - rt;
        sgn2 = -1;
    }
    acs = fabs( cs );
    if (acs>ab) {
        ct = -tb / cs;
        *sn1 = 1.0 / sqrt( 1.0+ct*ct );
        *cs1 = ct**sn1;
    } else if (!ab) {
        *cs1 = 1.0;
        *sn1 = 0.0;
    } else {
        tn = -cs / tb;
        *cs1 = 1.0 / sqrt( 1.0+tn*tn );
        *sn1 = tn**cs1;
    }
    if (sgn1==sgn2) {
        tn = *cs1;
        *cs1 = -*sn1;
        *sn1 = tn;
    }
    return;
}
//---------------------------------------------------------------------------
void fsDavidson::dlagtf(ntyp n, ftyp a[], ftyp lambda, ftyp b[], ftyp c[],
			ftyp tol, ftyp d[], ntyp in[], ntyp *info)
/*  -- LAPACK routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*
*  Purpose
*  =======
*
*  DLAGTF factorizes the matrix (T - lambda*I), where T is an n by n
*  tridiagonal matrix and lambda is a scalar, as
*
*     T - lambda*I = PLU,
*
*  where P is a permutation matrix, L is a unit lower tridiagonal matrix
*  with at most one non-zero sub-diagonal elements per column and U is
*  an upper triangular matrix with at most two non-zero super-diagonal
*  elements per column.
*
*  The factorization is obtained by Gaussian elimination with partial
*  pivoting and implicit row scaling.
*
*  The parameter LAMBDA is included in the routine so that DLAGTF may
*  be used, in conjunction with DLAGTS, to obtain eigenvectors of T by
*  inverse iteration.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix T.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (N)
*          On entry, A must contain the diagonal elements of T.
*
*          On exit, A is overwritten by the n diagonal elements of the
*          upper triangular matrix U of the factorization of T.
*
*  LAMBDA  (input) DOUBLE PRECISION
*          On entry, the scalar lambda.
*
*  B       (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, B must contain the (n-1) super-diagonal elements of
*          T.
*
*          On exit, B is overwritten by the (n-1) super-diagonal
*          elements of the matrix U of the factorization of T.
*
*  C       (input/output) DOUBLE PRECISION array, dimension (N-1)
*          On entry, C must contain the (n-1) sub-diagonal elements of
*          T.
*
*          On exit, C is overwritten by the (n-1) sub-diagonal elements
*          of the matrix L of the factorization of T.
*
*  TOL     (input) DOUBLE PRECISION
*          On entry, a relative tolerance used to indicate whether or
*          not the matrix (T - lambda*I) is nearly singular. TOL should
*          normally be chose as approximately the largest relative error
*          in the elements of T. For example, if the elements of T are
*          correct to about 4 significant figures, then TOL should be
*          set to about 5*10**(-4). If TOL is supplied as less than eps,
*          where eps is the relative machine precision, then the value
*          eps is used in place of TOL.
*
*  D       (output) DOUBLE PRECISION array, dimension (N-2)
*          On exit, D is overwritten by the (n-2) second super-diagonal
*          elements of the matrix U of the factorization of T.
*
*  IN      (output) INTEGER array, dimension (N)
*          On exit, IN contains details of the permutation matrix P. If
*          an interchange occurred at the kth step of the elimination,
*          then IN(k) = 1, otherwise IN(k) = 0. The element IN(n)
*          returns the smallest positive integer j such that
*
*             abs( u(j,j) ).le. norm( (T - lambda*I)(j) )*TOL,
*
*          where norm( A(j) ) denotes the sum of the absolute values of
*          the jth row of the matrix A. If no such j exists then IN(n)
*          is returned as zero. If IN(n) is returned as positive, then a
*          diagonal element of U is small, indicating that
*          (T - lambda*I) is singular or nearly singular,
*
*  INFO    (output) INTEGER
*          = 0   : successful exit
*          .lt. 0: if INFO = -k, the kth argument had an illegal value
*/
{
	ftyp piv1, piv2;

    *info=0;
    if (n<0) {
        *info=-1;
        xerbla("dlagtf ",-*info);
        return;
    }
    if (!n) return;
    a[0] -= lambda;
    in[ n  -1] = 0;
    if (n==1)
    	if (!a[0]) {
            in[0]=1;
            return;
        }
    ftyp eps = dlamch( 'e' ),
         tl = MAX( tol, eps ),
         scale1 = fabs( a[0] ) + fabs( b[0] );
    for (ntyp k=1; k<=n-1; k++) {
        a[k] -= lambda;
        ftyp scale2 = fabs( c[ k  -1] ) + fabs( a[k] );
        if (k<(n-1)) scale2 += fabs( b[k] );
        if (!a[k -1]) piv1 = 0.0;
            else piv1 = fabs( a[ k  -1] ) / scale1;
        if (!c[k -1]) {
            in[k -1] = 0;
            piv2 = 0.0;
            scale1 = scale2;
            if (k<(n-1)) d[k -1]=0.0;
        } else {
            piv2 = fabs( c[k -1] ) / scale2;
            if (piv2<=piv1) {
                in[ k  -1] = 0;
                scale1 = scale2;
                c[ k  -1] /= a[ k  -1];
                a[k] -= c[ k  -1]*b[ k  -1];
                if (k<(n-1)) d[k -1]=0.0;
            } else {
                in[k -1]=1;
                ftyp mult = a[ k  -1] / c[ k  -1],
                     temp = a[k];
                a[ k  -1] = c[ k  -1];
                a[k] = b[ k  -1] - mult*temp;
                if (k<(n-1)) {
                    d[ k  -1] = b[k];
                    b[k] = -mult*d[ k  -1];
                }
                b[ k  -1] = temp;
                c[ k  -1] = mult;
            }
        }
        if ((MAX(piv1,piv2)<=tl) && (!in[n -1])) in[n -1]=k;
    }
	if ((fabs(a[n -1])<=scale1*tl) && (!in[n -1])) in[n -1]=n;
    return;
}
//---------------------------------------------------------------------------
void fsDavidson::dladtsf(ftyp *absak, ftyp *ak, ftyp *pert, ftyp *temp,
			 ftyp one, ftyp zero, ftyp sfmin, ftyp bignum)
{
	*absak = fabs( *ak );
    if (*absak<one)
    	if (*absak<sfmin){
        	if ( (*absak==zero) || (fabs( *temp )*sfmin>*absak) ) {
            	*ak += *pert;
                *pert *= 2;
                dladtsf(absak, ak, pert, temp, one, zero, sfmin, bignum);
            } else {
            	*temp *= bignum;
                *ak *= bignum;
            }
        } else if ( fabs( *temp )>(*absak*bignum) ) {
        	*ak += *pert;
            *pert *= 2;
            dladtsf(absak, ak, pert, temp, one, zero, sfmin, bignum);
        }
     return;
}


void fsDavidson::dlagts(ntyp job, ntyp n, ftyp a[], ftyp b[], ftyp c[], ftyp d[],
			ntyp in[], ftyp y[], ftyp *tol, ntyp *info )
/*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     October 31, 1992
*
*
*  Purpose
*  =======
*
*  DLAGTS may be used to solve one of the systems of equations
*
*     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
*
*  where T is an n by n tridiagonal matrix, for x, following the
*  factorization of (T - lambda*I) as
*
*     (T - lambda*I) = P*L*U ,
*
*  by routine DLAGTF. The choice of equation to be solved is
*  controlled by the argument JOB, and in each case there is an option
*  to perturb zero or very small diagonal elements of U, this option
*  being intended for use in applications such as inverse iteration.
*
*  Arguments
*  =========
*
*  JOB     (input) INTEGER
*          Specifies the job to be performed by DLAGTS as follows:
*          =  1: The equations  (T - lambda*I)x = y  are to be solved,
*                but diagonal elements of U are not to be perturbed.
*          = -1: The equations  (T - lambda*I)x = y  are to be solved
*                and, if overflow would otherwise occur, the diagonal
*                elements of U are to be perturbed. See argument TOL
*                below.
*          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
*                but diagonal elements of U are not to be perturbed.
*          = -2: The equations  (T - lambda*I)'x = y  are to be solved
*                and, if overflow would otherwise occur, the diagonal
*                elements of U are to be perturbed. See argument TOL
*                below.
*
*  N       (input) INTEGER
*          The order of the matrix T.
*
*  A       (input) DOUBLE PRECISION array, dimension (N)
*          On entry, A must contain the diagonal elements of U as
*          returned from DLAGTF.
*
*  B       (input) DOUBLE PRECISION array, dimension (N-1)
*          On entry, B must contain the first super-diagonal elements of
*          U as returned from DLAGTF.
*
*  C       (input) DOUBLE PRECISION array, dimension (N-1)
*          On entry, C must contain the sub-diagonal elements of L as
*          returned from DLAGTF.
*
*  D       (input) DOUBLE PRECISION array, dimension (N-2)
*          On entry, D must contain the second super-diagonal elements
*          of U as returned from DLAGTF.
*
*  IN      (input) INTEGER array, dimension (N)
*          On entry, IN must contain details of the matrix P as returned
*          from DLAGTF.
*
*  Y       (input/output) DOUBLE PRECISION array, dimension (N)