fsdavidson.cpp

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

               }
        }
    if (!alpha) return;
    if (lsame(trans, 'n')){
//  Form  y := alpha*A*x + y.
       jx=kx;
       if (incy==1){
       	for (ntyp j=1; j<=n; j++){
               if (x[jx -1]){
               	 temp=alpha*x[jx -1];
                   for (ntyp i=1; i<=m; i++) y[i -1]+=temp*a(i,j);
               }
               jx += incx;
           }
       } else {
           for (ntyp j=1; j<=n; j++){
               if (x[jx -1]){
                   temp=alpha*x[jx -1];
                   iy=ky;
                   for (ntyp i=1; i<=m; i++){
                       y[iy -1]+=temp*a(i,j);
                       iy+=incy;
                   }
               }
               jx+=incx;
           }
       }
    } else {
    	//Form  y := alpha*A'*x + y.
        jy=ky;
        if (incx==1){
           for (ntyp j=1; j<=n; j++){
                temp=0;
                for (ntyp i=1; i<=m; i++) temp+=a(i,j)*x[i -1];
                y[jy -1]+=alpha*temp;
                jy+=incy;
           }
        } else {
               for (ntyp j=1; j<=n; j++){
                   temp=0;
                   ix=kx;
                   for (ntyp i=1; i<=m; i++){
                       temp+=a(i,j)*x[ix -1];
                       ix+=incx;
                   }
                   y[jy -1]+=alpha*temp;
                   jy+=incy;
                }
        }
    }
    return;
}

#undef a //(j,i)
//---------------------------------------------------------------------------
#define a(j,i) a[((((i)-1)*(lda))+((j)-1))]

void fsDavidson::dger(ntyp  m, ntyp n, ftyp alpha, ftyp x[], ntyp incx, ftyp y[],
 		  ntyp incy, ftyp a[], ntyp lda )
/*  Purpose
*  =======
*
*  DGER   performs the rank 1 operation
*
*     A := alpha*x*y' + A,
*
*  where alpha is a scalar, x is an m element vector, y is an n element
*  vector and A is an m by n matrix.
*
*  Parameters
*  ==========
*
*  M      - INTEGER.
*           On entry, M specifies the number of rows of the matrix A.
*           M must be at least zero.
*           Unchanged on exit.
*
*  N      - INTEGER.
*           On entry, N specifies the number of columns of the matrix A.
*           N must be at least zero.
*           Unchanged on exit.
*
*  ALPHA  - DOUBLE PRECISION.
*           On entry, ALPHA specifies the scalar alpha.
*           Unchanged on exit.
*
*  X      - DOUBLE PRECISION array of dimension at least
*           ( 1 + ( m - 1 )*abs( INCX ) ).
*           Before entry, the incremented array X must contain the m
*           element vector x.
*           Unchanged on exit.
*
*  INCX   - INTEGER.
*           On entry, INCX specifies the increment for the elements of
*           X. INCX must not be zero.
*           Unchanged on exit.
*
*  Y      - DOUBLE PRECISION array of dimension at least
*           ( 1 + ( n - 1 )*abs( INCY ) ).
*           Before entry, the incremented array Y must contain the n
*           element vector y.
*           Unchanged on exit.
*
*  INCY   - INTEGER.
*           On entry, INCY specifies the increment for the elements of
*           Y. INCY must not be zero.
*           Unchanged on exit.
*
*  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
*           Before entry, the leading m by n part of the array A must
*           contain the matrix of coefficients. On exit, A is
*           overwritten by the updated matrix.
*
*  LDA    - INTEGER.
*           On entry, LDA specifies the first dimension of A as declared
*           in the calling (sub) program. LDA must be at least
*           max( 1, m ).
*           Unchanged on exit.
*
*
*  Level 2 Blas routine.
*
*  -- Written on 22-October-1986.
*     Jack Dongarra, Argonne National Lab.
*     Jeremy Du Croz, Nag Central Office.
*     Sven Hammarling, Nag Central Office.
*     Richard Hanson, Sandia National Labs.		*/
{
    ntyp info=0, ix, jy, kx;

    //Test the input parameters.
    if (m<0) info=1;
    	else if (n<0) info=2;
        	else if (!incx) info=5;
            	else if (!incy) info=7;
                	else if (lda< MAX(1,m)) info=9;
    if (info){
        xerbla("dger",info);
        return;
    }

    //Quick return if possible.
    if ((!m) || (!n) || (!alpha)) return;

//  Start the operations. In this version the elements of A are
//  accessed sequentially with one pass through A.
    if (incy>0) jy=1;
    	else jy = 1 - ( n - 1 )*incy;
    if (incx==1){
        for (ntyp j=1; j<=n; j++){
            if (y[jy -1]){
                ftyp temp=alpha*y[jy -1];
                for (ntyp i=1; i<=m; i++) a(i,j)+=x[i -1]*temp;
            }
            jy += incy;
        }
    } else {
        if (incx>0) kx=1;
            else kx = 1 - ( m - 1 )*incx;
        for (ntyp j=1; j<=n; j++){
            if (y[jy -1]){
                ftyp temp=alpha*y[jy -1];
                ix=kx;
                for (ntyp i=1; i<=m; i++){
                    a(i,j) += x[ix -1]*temp;
                    ix     += incx;
                }
            }
            jy += incy;
        }
    }
    return;
}

#undef a //(j,i)
//---------------------------------------------------------------------------
void fsDavidson::dinit(ntyp n, ftyp a, ftyp x[], ntyp incx )
/*=======================================================================
*       PURPOSE ... INITIALIZES DOUBLE PRECISION VECTOR TO              
*                   A CONSTANT VALUE 'A'                                
*=======================================================================*/
{
    if  ( incx == 1 )
        for (ntyp i=1; i<=n; i++) x[i -1] = a;
    else {
        ntyp xaddr = 1;
        if  ( incx < 0 ) xaddr = (-n+1)*incx + 1;
        for (ntyp i=1; i<=n; i++) {
            x[xaddr -1] = a;
            xaddr += incx;
        }
    }
	return;
}
//---------------------------------------------------------------------------
void fsDavidson::dlae2(ftyp a, ftyp b, ftyp c, ftyp *rt1, ftyp *rt2 )
/*  -- 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
*  =======
*
*  DLAE2  computes the eigenvalues of a 2-by-2 symmetric matrix
*     [  A   B  ]
*     [  B   C  ].
*  On return, RT1 is the eigenvalue of larger absolute value, and RT2
*  is the eigenvalue of smaller absolute value.
*
*  Arguments
*  =========
*
*  A       (input) DOUBLE PRECISION
*          The (1,1) element of the 2-by-2 matrix.
*
*  B       (input) DOUBLE PRECISION
*          The (1,2) and (2,1) elements 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.
*
*  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.
*
*  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, rt,
         sm = a + c, df = a - c, adf = fabs( df ), tb = b + b, ab = fabs( tb );

    //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 ));
            //Includes case AB=ADF=0
        	else rt = ab*sqrt( 2.0 );
    if (sm<0) {
        *rt1 = 0.5*( sm-rt );

        /*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) {
        *rt1 = 0.5*( sm+rt );

        /*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;
    }
	return;
}
//---------------------------------------------------------------------------
#define ab(j,i) ab[((((i)-1)*(mmax))+((j)-1))]
#define nab(j,i) nab[((((i)-1)*(mmax))+((j)-1))]

void fsDavidson::dlaebz(ntyp ijob, ntyp nitmax, ntyp n, ntyp mmax, ntyp minp,
			ntyp nbmin, ftyp abstol, ftyp reltol, ftyp pivmin, ftyp d[],
            ftyp e[], ftyp e2[], ntyp nval[], ftyp ab[], ftyp c[], ntyp *mout,
            ntyp nab[], ftyp work[], ntyp iwork[], 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
*     June 30, 1999
*
*  Purpose
*  =======
*
*  DLAEBZ contains the iteration loops which compute and use the
*  function N(w), which is the count of eigenvalues of a symmetric
*  tridiagonal matrix T less than or equal to its argument  w.  It
*  performs a choice of two types of loops:
*
*  IJOB=1, followed by
*  IJOB=2: It takes as input a list of intervals and returns a list of
*          sufficiently small intervals whose union contains the same
*          eigenvalues as the union of the original intervals.
*          The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
*          The output interval (AB(j,1),AB(j,2)] will contain
*          eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT.
*
*  IJOB=3: It performs a binary search in each input interval
*          (AB(j,1),AB(j,2)] for a point  w(j)  such that
*          N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
*          the search.  If such a w(j) is found, then on output
*          AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
*          (AB(j,1),AB(j,2)] will be a small interval containing the
*          point where N(w) jumps through NVAL(j), unless that point
*          lies outside the initial interval.
*
*  Note that the intervals are in all cases half-open intervals,
*  i.e., of the form  (a,b] , which includes  b  but not  a .
*
*  To avoid underflow, the matrix should be scaled so that its largest
*  element is no greater than  overflow**(1/2) * underflow**(1/4)
*  in absolute value.  To assure the most accurate computation
*  of small eigenvalues, the matrix should be scaled to be
*  not much smaller than that, either.
*
*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
*  Matrix", Report CS41, Computer Science Dept., Stanford
*  University, July 21, 1966
*
*  Note: the arguments are, in general, *not* checked for unreasonable
*  values.
*
*  Arguments
*  =========
*
*  IJOB    (input) INTEGER
*          Specifies what is to be done:
*          = 1:  Compute NAB for the initial intervals.
*          = 2:  Perform bisection iteration to find eigenvalues of T.
*          = 3:  Perform bisection iteration to invert N(w), i.e.,
*                to find a point which has a specified number of
*                eigenvalues of T to its left.
*          Other values will cause DLAEBZ to return with INFO=-1.
*
*  NITMAX  (input) INTEGER
*          The maximum number of "levels" of bisection to be
*          performed, i.e., an interval of width W will not be made
*          smaller than 2^(-NITMAX) * W.  If not all intervals
*          have converged after NITMAX iterations, then INFO is set
*          to the number of non-converged intervals.
*
*  N       (input) INTEGER
*          The dimension n of the tridiagonal matrix T.  It must be at
*          least 1.
*
*  MMAX    (input) INTEGER
*          The maximum number of intervals.  If more than MMAX intervals
*          are generated, then DLAEBZ will quit with INFO=MMAX+1.
*
*  MINP    (input) INTEGER