fsdavidson.cpp

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

*          On entry, the right hand side vector y.
*          On exit, Y is overwritten by the solution vector x.
*
*  TOL     (input/output) DOUBLE PRECISION
*          On entry, with  JOB .lt. 0, TOL should be the minimum
*          perturbation to be made to very small diagonal elements of U.
*          TOL should normally be chosen as about eps*norm(U), where eps
*          is the relative machine precision, but if TOL is supplied as
*          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
*          If  JOB .gt. 0  then TOL is not referenced.
*
*          On exit, TOL is changed as described above, only if TOL is
*          non-positive on entry. Otherwise TOL is unchanged.
*
*  INFO    (output) INTEGER
*          = 0   : successful exit
*          .lt. 0: if INFO = -i, the i-th argument had an illegal value
*          .gt. 0: overflow would occur when computing the INFO(th)
*                  element of the solution vector x. This can only occur
*                  when JOB is supplied as positive and either means
*                  that a diagonal element of U is very small, or that
*                  the elements of the right-hand side vector y are very
*                  large.
*/
{
	ftyp absak, ak, pert, temp;

    *info = 0;
    if (( fabs( job )>2 ) || ( !job )) *info = -1;
    	else if (n<0) *info = -2;
    if (*info){
    		xerbla( "dlagts ", -*info );
            return;
    }
    if (!n) return;

	ftyp eps = dlamch( 'e' ),
    	 sfmin = dlamch( 's' ),
         bignum = 1.0 / sfmin;

	if (job<0) {
    	if (*tol<=0.0) {
        	*tol = fabs( a[0] );
            if (n>1) *tol = MAX( (MAX(*tol, fabs(a[1]))), fabs(b[0]) );
            for (ntyp k=3; k<=n; k++)
            	*tol = MAX( (MAX(*tol, fabs(a[k -1]))), (MAX(fabs(b[k-2]),fabs(d[k-3]))) );
            *tol*=eps;
            if (!(*tol)) *tol = eps;
        }
    }
    if (fabs( job )==1) {
    	for (ntyp k=2; k<=n; k++)
        	if (!in[k-2]) y[k -1] = y[k -1] - c[k-2]*y[k-2];
            	else {
                	temp = y[k-2];
                    y[k-2] = y[k -1];
                    y[k -1] = temp - c[k-2]*y[k -1];
                }
        if (job==1)
        	for (ntyp k=n; k>=1; k--) {
            	if (k<=(n-2)) temp = y[k -1] - b[k -1]*y[k] - d[k -1]*y[k+1];
                	else if (k==(n-1)) temp = y[k -1] - b[k -1]*y[k];
                    	else temp = y[k -1];
            	ak = a[k -1];
            	absak = fabs( ak );
            	if (absak<1.0)
            		if (absak<sfmin) {
                    	if ( !absak || (fabs( temp )*sfmin)>(absak) ) {
                    		*info = k;
                        	return;
                    	} else {
                    		temp *= bignum;
                        	ak *= bignum;
                    	}
                	} else if ( fabs( temp )>(absak*bignum) ) {
                		*info = k;
                    	return;
            		}
            	y[k -1] = temp / ak;
        	}
        else {
        	for (ntyp k=n; k>=1; k--){
            	if (k<=(n-2)) temp = y[k -1] - b[k -1]*y[k] - d[k -1]*y[k+1];
                	else if ( k==(n-1)) temp = y[k -1] - b[k -1]*y[k];
                    	else temp = y[k -1];
                ak = a[k -1];
                pert = sign( tol, &ak );
                dladtsf(&absak, &ak, &pert, &temp, 1.0, 0.0, sfmin, bignum);
                y[k -1] = temp / ak;
            }
        }
    } else {
    	//Come to here if  JOB = 2 or -2
    	if (job==2)
        	for (ntyp k=1; k<=n; k++) {
            	if (k>=3) temp = y[k -1] - b[k-2]*y[k-2] - d[k-3]*y[k-3];
                	else if (k==2) temp = y[k -1] - b[k-2]*y[k-2];
                    	else temp = y[k -1];
                ak = a[k -1];
               	absak = fabs( ak );
                if ( absak<1.0 )
                	if (absak<sfmin) {
                    	if ( (!absak) || (fabs( temp )*sfmin>absak) ) {
                        	*info = k;
                            return;
                        } else {
                        	temp *= bignum;
                            ak *= bignum;
                        }
                    } else if ( fabs( temp )>absak*bignum ) {
                    	*info = k;
                        return;
                    }
                y[k -1] = temp / ak;
        	}
        else
        	for (ntyp k=1; k<=n; k++) {
            	if (k>=3) temp = y[k -1] - b[k-2]*y[k-2] - d[k-3]*y[k-3];
                	else if (k==2) temp = y[k -1] - b[k-2]*y[k-2];
                    	else temp = y[k -1];
                ak = a[k -1];
            	pert = sign( tol, &ak );
                dladtsf(&absak, &ak, &pert, &temp, 1.0, 0.0, sfmin, bignum);
                y[k -1] = temp / ak;

            }
            for (ntyp k=n; k>=2; k--) {
            	if ( !in[k-2] ) y[k-2] -= c[k-2]*y[k -1];
                	else {
                    	temp = y[k-2];
               			y[k-2] = y[k -1];
               			y[k -1] = temp - c[k-2]*y[k -1];
                    }
            }
    }
	return;
}
//---------------------------------------------------------------------------
/* DLAMC3  is intended to force  A  and  B  to be stored prior to doing
*  the addition of  A  and  B ,  for use in situations where optimizers
*  might hold one of these in a register.    */
#define dlamc3(a,b) ((a)+(b))

void dlamc1(ntyp *beta, ntyp *t, bool *rnd, bool *ieee1 )
/*  -- 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
*  =======
*
*  DLAMC1 determines the machine parameters given by BETA, T, RND, and
*  IEEE1.
*
*  Arguments
*  =========
*
*  BETA    (output) INTEGER
*          The base of the machine.
*
*  T       (output) INTEGER
*          The number of ( BETA ) digits in the mantissa.
*
*  RND     (output) LOGICAL
*          Specifies whether proper rounding  ( RND = .TRUE. )  or
*          chopping  ( RND = .FALSE. )  occurs in addition. This may not
*          be a reliable guide to the way in which the machine performs
*          its arithmetic.
*
*  IEEE1   (output) LOGICAL
*          Specifies whether rounding appears to be done in the IEEE
*          'round to nearest' style.
*
*  Further Details
*  ===============
*
*  The routine is based on the routine  ENVRON  by Malcolm and
*  incorporates suggestions by Gentleman and Marovich. See
*
*     Malcolm M. A. (1972) Algorithms to reveal properties of
*        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
*
*     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
*        that reveal properties of floating point arithmetic units.
*        Comms. of the ACM, 17, 276-277.
*/
{
    volatile static bool first=1, lieee1, lrnd;
    volatile static ntyp lbeta, lt;
    volatile ftyp a, b, c, f, one, qtr, savec, t1, t2;

    if (first) {
    	first = 0;
        one = 1;
/*       LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
*        IEEE1, T and RND.
*
*        Throughout this routine  we use the function  DLAMC3  to ensure
*        that relevant values are  stored and not held in registers,  or
*        are not affected by optimizers.
*
*        Compute  a = 2.0**m  with the  smallest positive integer m such
*        that
*
*           fl( a + 1.0 ) = a.	*/
        a = 1;
        c = 1;
        for (;;)
        	if (c==one){
            	a *= 2;
                c = dlamc3( a, one );
                c = dlamc3( c, -a );
            } else break;
/*      Now compute  b = 2.0**m  with the smallest positive integer m
*       such that
*
*       fl( a + b ) .gt. a.                 */
        b = 1;
        c = dlamc3( a, b );
        for (;;)
        	if (c==a) {
            	b *= 2;
                c = dlamc3( a, b );
            } else break;
/*        Now compute the base.  a and c  are neighbouring floating point
*        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
*        their difference is beta. Adding 0.25 to c is to ensure that it
*        is truncated to beta and not ( beta - 1 ).  */
        qtr = one / 4;
        savec = c;
        c = dlamc3( c, -a );
        lbeta = ntyp(c + qtr);
//      Now determine whether rounding or chopping occurs,  by adding a
//      bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
        b = lbeta;
        f = dlamc3( b / 2, -b / 100 );
        c = dlamc3( f, a);
        if (c==a) lrnd = 1;
        	else lrnd = 0;
        f = dlamc3( b / 2, b / 100 );
        c = dlamc3( f, a );
        if( ( lrnd ) && ( c==a ) ) lrnd = 0;
/*       Try and decide whether rounding is done in the  IEEE  'round to
*        nearest' style. B/2 is half a unit in the last place of the two
*        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
*        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
*        A, but adding B/2 to SAVEC should change SAVEC. */
        t1 = dlamc3( b / 2, a );
        t2 = dlamc3( b / 2, savec );
        lieee1 = ( t1==a ) && ( t2>savec ) && lrnd;
/*       Now find  the  mantissa, t.  It should  be the  integer part of
*        log to the base beta of a,  however it is safer to determine  t
*        by powering.  So we find t as the smallest positive integer for
*        which
*
*           fl( beta**t + 1.0 ) = 1.0. */
        lt = 0;
        a = c = 1;
        for (;;)
        	if (c==one) {
            	++lt;
                a *= lbeta;
                c = dlamc3( a, one );
                c = dlamc3( c, -a );
            } else break;
    }
    *beta = lbeta;
    *t = lt;
    *rnd = lrnd;
    *ieee1 = lieee1;

    return;
}

void dlamc4(ntyp *emin, ftyp start, ntyp base )
/*  -- 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
*  =======
*
*  DLAMC4 is a service routine for DLAMC2.
*
*  Arguments
*  =========
*
*  EMIN    (output) EMIN
*          The minimum exponent before (gradual) underflow, computed by
*          setting A = START and dividing by BASE until the previous A
*          can not be recovered.
*
*  START   (input) DOUBLE PRECISION
*          The starting point for determining EMIN.
*
*  BASE    (input) INTEGER
*          The base of the machine.
*/
{
    volatile ftyp a, b1, b2, c1, c2, d1, d2, one, rbase, zero;

    a = start;
    one = 1;
    rbase = one / base;
    zero = 0;
    *emin = 1;
    b1 = dlamc3( a*rbase, zero );
    c1 = a;
    c2 = a;
    d1 = a;
    d2 = a;
    for (;;)
    	if (( c1==a ) && ( c2==a ) && ( d1==a ) && ( d2==a ) ) {
        	--*emin;
         	a = b1;
         	b1 = dlamc3( a / base, zero );
         	c1 = dlamc3( b1*base, zero );
         	d1 = zero;
            for (ntyp i=1; i<=base; i++) d1 = d1 + b1;
            b2 = dlamc3( a*rbase, zero );
         	c2 = dlamc3( b2 / rbase, zero );
         	d2 = zero;
            for (ntyp i=1; i<=base; i++) d2 = d2 + b2;
        } else break;

    return;
}

void dlamc5(ntyp beta, ntyp p, ntyp emin, bool ieee, ntyp *emax, ftyp *rmax )
/*  -- 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
*  =======
*
*  DLAMC5 attempts to compute RMAX, the largest machine floating-point
*  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
*  approximately to a power of 2.  It will fail on machines where this
*  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
*  EMAX = 28718).  It will also fail if the value supplied for EMIN is
*  too large (i.e. too close to zero), probably with overflow.
*
*  Arguments
*  =========
*
*  BETA    (input) INTEGER
*          The base of floating-point arithmetic.
*
*  P       (input) INTEG