comp_with_fs.cpp

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

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include "../davidson/comp_with_fs.h"

CompareWithFS::CompareWithFS()
{
	Nmax=100;
    IBAND=10;
    NZERmax = (IBAND+1)*(Nmax-IBAND/2);

	N=Nmax;

    A = new ftyp[NZERmax];
    IndCol = new ntyp[Nmax];
    IndRow = new ntyp[NZERmax];

    DIAG = new ftyp[Nmax];

    Init();
}

void CompareWithFS::Init()
{
	ntyp icur=1;
	for(ntyp i=1; i<=N; i++){
    	A[icur-1]=i;
        IndRow[icur-1]=i;
        IndCol[i-1]=icur;
        ++icur;
        for(ntyp j=i+1; j<=MIN(N,i+IBAND); j++){
        	A[icur-1]=0.001;
            IndRow[icur-1]=j;
            IndCol[i-1]=icur;
            ++icur;
        }
    }
    lupper=false;

    DIAG[0]=A[0];
    for(ntyp i=2; i<=N; i++)
    	DIAG[i-1]=A[IndCol[i-1 -1]+1 -1];

    return;
}

CompareWithFS::~CompareWithFS()
{
	delete[] A;
    delete[] IndCol;
    delete[] IndRow;
    delete[] DIAG;
}

ntyp CompareWithFS::GetN() const
{
	return N;
}


void CompareWithFS::dssbmv(ntyp m, ftyp b[], ftyp c[], ntyp n) const
{

     ntyp ioffs,istart,icur, numelem;
     ftyp diag;

    ftyp *TEMPB = new ftyp[Nmax];
    ftyp *TEMPC = new ftyp[Nmax];

     ioffs=1;
     if (lupper) ioffs=0;
     istart=1;
     dinit(n*m,0,c,1);
     for (int icol=1; icol<=n; icol++)
     	{
        	numelem = IndCol[icol -1]-istart+1;
            icur=1;
            for (int iv=1; iv<=m; iv++)
            	{
                    gather(numelem,TEMPB,&b[icur -1],&IndRow[istart -1]);
                    gather(numelem,TEMPC,&c[icur -1],&IndRow[istart -1]);
                    diag=c[icur-1+icol -1]+ddot(numelem,&A[istart -1],1,TEMPB,1);
                    daxpy(numelem-1,b[icur-1+icol -1], &A[istart+ioffs -1],1, &TEMPC[ioffs+1 -1],1);
                    scatter(numelem,&c[icur -1],&IndRow[istart -1],TEMPC);
                    c[icur-1+icol -1]=diag;
              		icur=icur+n;
                }
            istart=IndCol[icol -1]+1;
        }

    delete[] TEMPB;
    delete[] TEMPC;

	return;
}


ftyp *CompareWithFS::Diagonal() const
{
	return DIAG;
}

void CompareWithFS::dinit(ntyp n, ftyp a, ftyp x[], ntyp incx) const
/*=======================================================================
*       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 CompareWithFS::gather(ntyp n, ftyp a[], ftyp b[], ntyp index[]) const
/*      This subroutine collects array elements accessed via an     
*       integer pointer to contiguous storage. 	*/
{
    for (ntyp i=1; i<=n; i++)
    	a[i -1] = b[index[i -1] -1];

    return;
}

void CompareWithFS::scatter(ntyp n, ftyp a[], ntyp index[], ftyp b[]) const
/* This subroutine disperses array elements accessed via an integer
*  pointer from contiguous storage to the appropriate location. */
{
    for (ntyp i=1; i<=n; i++)
    	a[index[i -1] -1] = b[i -1];

    return;
}

ftyp CompareWithFS::ddot(ntyp n, ftyp dx[], ntyp incx, ftyp dy[], ntyp incy) const
{
/*   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;
}

void CompareWithFS::daxpy(ntyp n, ftyp da, ftyp dx[], ntyp incx, ftyp dy[], ntyp incy) const
{
/*   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 CompareWithFS::print(Davidson *dav) const
{
    eig_val eigen;

    cout<<"\n\nIERR =  "<<dav->GetIerr()<<"\tMatrix accesses =  "<<dav->loop
    	<<"\tMatrix-Vector products =  "<<dav->GetMatrixVectProd()<<endl<<endl;
    cout<<"  Descending Order: "<<dav->GetHiend()<<endl<<endl;
    cout<<"\n\tEigenvalues\t\t"<<"Eigval Differences\t\t"<<"Residuals"<<endl<<endl;
    cout.setf(ios::scientific);
    cout.precision(15);
    for(ntyp i=1; i<=10; i++){
    	dav->EigenValue(eigen,i);
        cout<<"    "<<eigen.Eigenvalue<<"     "<<eigen.Eigval_differences<<"     "<<eigen.Residual<<endl;
    }

    GridParam grid(-10,10,100);
    ftyp *Psi_value = new ftyp[100];
    cout<<"\n First six components of the eigenvectors"<<endl<<endl;
    for(ntyp i=1; i<=10; i++){
    	dav->EigenFunction(i,grid,Psi_value,0);
        for(ntyp j=0; j<6; j++){
            printf("  %10.3e",Psi_value[j]);
        }
        cout<<endl;
    }
    delete[] Psi_value;
    return;
}