fsdavidson.cpp

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

#include "../davidson/fsdavidson.h"
#include "../davidson/defaulttd.h"
#include <math.h>
#include <stdio.h>
#include <iostream>

//---------------------------------------------------------------------------
void fsDavidson::addabs(ntyp n, ntyp lim, bool hiend, ntyp kpass, ntyp nncv,
		    ftyp basis[], ftyp ab[], ftyp s[],
            const MatrixBase *const diag)
/*-----------------------------------------------------------------------
*       Called by: DVDRVR, SETUP                                        
*                                                                       
*       Calculates the new column in the D matrix and the new column    
*       in the S matrix. The new D column is D(new)=AB(new). S has a    
*       new row and column, but being symmetric only the new column is  
*       stored. S(i,kpass+1)=B(i)^T D(kpass+1) for all i.               
*                                                                       
*       subroutines called:                                             
*       MatrixBase::dssbmv, DDOT, DSCAL
*
*   on entry                                                            
*   -------                                                             
*   N           The order of the matrix A                               
*   kpass       The current dimension of the expanding sub-basis        
*   NNCV        Number of new basis vectors.                            
*   Basis       the basis vectors, including the new NNCV ones.
*	diag        The object of MatrixBase class
*
*   on exit                                                             
*   -------                                                             
*   AB          The new matrix D=AB. (with new NNCV columns)            
*   S           The small matrix with NNCV new columns at the last part 
*-----------------------------------------------------------------------*/
{
// The user specified matrix-vector routine is called with the new
// basis vector B(*,kpass+1) and the result is assigned to AB(idstart)
    int idstart=kpass*n+1;

    diag->dssbmv(nncv,&basis[idstart -1],&ab[idstart -1],n);

// If highest pairs are sought, use the negative of the matrix
    if (hiend) dscal(n*nncv,-1.0,&ab[idstart -1],1);

// The new S is calculated by adding the new last columns
// S(new)=B^T D(new).
    int isstart=kpass*(kpass+1)/2;
    for (int iv=1; iv<=nncv; iv++) {
    	int ibstart=1;
        for (int ibv=1; ibv<=kpass+iv; ibv++) {
           	s[isstart + ibv -1]=ddot(n,&basis[ibstart -1],1,&ab[idstart -1],1);
           	ibstart+=n;
        }
        isstart += (kpass+iv);
        idstart += n;
    }
	return;
}
//---------------------------------------------------------------------------
ftyp fsDavidson::dasum(ntyp n,ftyp dx[],ntyp incx)
{

/*   takes the sum of the absolute values.
     jack dongarra, linpack, 3/11/78.
     modified 3/93 to return if incx .le. 0.
     modified 12/3/93, array(1) declarations changed to array(*)  */

    ftyp dtemp=0.0;
    ntyp m,nincx;

    if ((n<=0) || (incx<=0)) return 0;
    if (incx!=1){
    	//code for increment not equal to 1
        nincx = n*incx;
        for (ntyp i=1;i<=nincx;i+=incx)
            dtemp += fabs(dx[i-1]);
        return dtemp;
    }
/*      code for increment equal to 1

        clean-up loop  */
    m=n%6;
    if (m){
        for (ntyp i=1;i<=m;i++) dtemp += fabs(dx[i-1]);
        if (n<6) return dtemp;
    }
    for (ntyp i=m+1;i<=n;i+=6)
        dtemp += fabs(dx[i-1]) + fabs(dx[i]) + fabs(dx[i+1])
            	 + fabs(dx[i+2]) + fabs(dx[i+3]) + fabs(dx[i+4]);

    return dtemp;
}
//---------------------------------------------------------------------------
void fsDavidson::daxpy(ntyp n, ftyp da, ftyp dx[], ntyp incx, ftyp dy[], ntyp incy)
{
/*   constant times a vector plus a vector.
     uses unrolled loops for increments equal to one.
     jack dongarra, linpack, 3/11/78.
     modified 12/3/93, array(1) declarations changed to array(*) */

      if (n<=0) return;
      if (da==0.0) return;
      if (!(incx==1 && incy==1)){
     	 //code for unequal increments or equal increments
         //not equal to 1
         ntyp ix=1, iy=1;
         if (incx<0) ix=(-n+1)*incx+1;
         if (incy<0) iy=(-n+1)*incy+1;
         for (ntyp i=1; i<=n; i++){
             	dy[iy-1] += da*dx[ix-1];
        	ix += incx;
        	iy += incy;
         }
         return;
      }
//    code for both increments equal to 1

//    clean-up loop
      ntyp m=n%4;
      if (m){
         for (ntyp i=1; i<=m; i++) dy[i-1]+=da*dx[i-1];
         if (n<4) return;
      }
      for (ntyp i=m + 1; i<=n; i+=4){
            dy[i-1] += da*dx[i-1];
            dy[i] += da*dx[i];
            dy[i + 1] += da*dx[i + 1];
            dy[i + 2] += da*dx[i + 2];
      }

      return;
}
//---------------------------------------------------------------------------
void fsDavidson::dcopy(ntyp n, ftyp dx[], ntyp incx, ftyp dy[], ntyp incy)
{
/*   copies a vector, x, to a vector, y.
     uses unrolled loops for increments equal to one.
     jack dongarra, linpack, 3/11/78.
     modified 12/3/93, array(1) declarations changed to array(*) */

    if (n<=0) return;
    if (!(incx==1 && incy==1)){
        //code for unequal increments or equal increments
        //not equal to 1
    	ntyp ix=1,iy=1;
        if (incx<0) ix=(-n+1)*incx+1;
        if (incy<0) iy=(-n+1)*incy+1;
        for (ntyp i=1; i<=n; i++){
            dy[iy-1]=dx[ix-1];
            ix+=incx;
            iy+=incy;
        }
        return;
    }
//  code for both increments equal to 1
//
//  clean-up loop
    ntyp m=n%7;
    if (m){
        for (ntyp i=1; i<=m; i++) dy[i-1]=dx[i-1];
        if (n<7) return;
    }
    for (ntyp i=m+1; i<=n; i+=7){
        dy[i-1] = dx[i-1];
        dy[i] = dx[i];
        dy[i + 1] = dx[i + 1];
        dy[i + 2] = dx[i + 2];
        dy[i + 3] = dx[i + 3];
        dy[i + 4] = dx[i + 4];
        dy[i + 5] = dx[i + 5];
    }
    return;
}
//---------------------------------------------------------------------------
ftyp fsDavidson::ddot(ntyp n, ftyp dx[], ntyp incx, ftyp dy[], ntyp incy)
{
/*   forms the dot product of two vectors.
     uses unrolled loops for increments equal to one.
     jack dongarra, linpack, 3/11/78.
     modified 12/3/93, array(1) declarations changed to array(*) */

    ftyp dtemp=0.0;

    if (n<=0) return 0.0;
    if (!(incx==1 && incy==1)){
        //code for unequal increments or equal increments
        //not equal to 1
        ntyp ix=1, iy=1;
        if (incx<0) ix = (-n+1)*incx + 1;
      	if (incy<0) iy = (-n+1)*incy + 1;
        for (ntyp i=1; i<=n; i++){
            dtemp += dx[ix -1]*dy[iy -1];
            ix += incx;
	    iy += incy;
        }
        return dtemp;
    }
//	code for both increments equal to 1
//
//  clean-up loop
    ntyp m=n%5;
    if (m){
            for (ntyp i=1; i<=m; i++) dtemp += dx[i -1]*dy[i -1];
            if (n<5) return dtemp;
    }
    for (ntyp i=m + 1; i<=n; i+=5)
        dtemp += dx[i -1]*dy[i -1] + dx[i]*dy[i] + dx[i + 1]*dy[i + 1]
                 + dx[i + 2]*dy[i + 2] + dx[i + 3]*dy[i + 3];
    return dtemp;
}
//---------------------------------------------------------------------------
#define a(j,i) a[((((i)-1)*(lda))+((j)-1))]

void fsDavidson::dgemv(char trans, ntyp m, ntyp n, ftyp alpha, ftyp a[], ntyp lda,
	   ftyp x[], ntyp incx, ftyp beta, ftyp y[], ntyp incy )
/*  Purpose
*  =======
*
*  DGEMV  performs one of the matrix-vector operations
*
*     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
*
*  where alpha and beta are scalars, x and y are vectors and A is an
*  m by n matrix.
*
*  Parameters
*  ==========
*
*  TRANS  - CHARACTER*1.
*           On entry, TRANS specifies the operation to be performed as
*           follows:
*
*              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
*
*              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
*
*              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
*
*           Unchanged on exit.
*
*  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.
*
*  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.
*           Unchanged on exit.
*
*  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.
*
*  X      - DOUBLE PRECISION array of DIMENSION at least
*           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
*           and at least
*           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
*           Before entry, the incremented array X must contain the
*           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.
*
*  BETA   - DOUBLE PRECISION.
*           On entry, BETA specifies the scalar beta. When BETA is
*           supplied as zero then Y need not be set on input.
*           Unchanged on exit.
*
*  Y      - DOUBLE PRECISION array of DIMENSION at least
*           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
*           and at least
*           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
*           Before entry with BETA non-zero, the incremented array Y
*           must contain the vector y. On exit, Y is overwritten by the
*           updated vector y.
*
*  INCY   - INTEGER.
*           On entry, INCY specifies the increment for the elements of
*           Y. INCY must not be zero.
*           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.    */
{
    ftyp temp;
    ntyp info=0, ix, iy, jx, jy, kx, ky, lenx, leny;

    //Test the input parameters.
    if (!lsame(trans,'n') && !lsame(trans,'t') && !lsame(trans,'c'))
    	info=1;
    else if (m<0) info=2;
    	else if (n<0) info=3;
	     else if (lda<MAX(1,m)) info=6;
            	else if (!incx) info=8;
                     else if (!incy) info=11;
    if (info){
        xerbla("dgemv ", info);
        return;
    }
    //Quick return if possible.
    if ((!m) || (!n) || ((!alpha) && (beta==1.0))) return;

    //Set  LENX  and  LENY, the lengths of the vectors x and y, and set
    //up the start points in  X  and  Y.
    if (lsame(trans,'N')){
        lenx = n;
        leny = m;
    } else {
        lenx = m;
        leny = n;
    }
    if (incx>0) kx=1;
    	else kx=1-(lenx-1)*incx;
    if (incy>0) ky=1;
    	else ky=1-(leny-1)*incy;

/*   Start the operations. In this version the elements of A are
     accessed sequentially with one pass through A.

     First form  y := beta*y.              */
    if (beta!=1.0)
    	if (incy==1){
           if (!beta)
            	for (ntyp i=1; i<=leny; i++) y[i -1]=0;
           else
            	for (ntyp i=1; i<=leny; i++) y[i -1]=beta*y[i -1];
        } else {
               iy = ky;
               if (!beta)
                  for (ntyp i=1; i<=leny; i++){
                      y[iy -1] = 0;
                      iy += incy;
                  }
               else {
                    for (ntyp i=1; i<=leny; i++){
                        y[iy -1] = beta*y[iy -1];
               	        iy += incy;
                    }