cdavidson.cpp

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

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <ctype.h>

#include "../davidson/cdavidson.h"
#include "../davidson/errdavidson.h"

Davidson::Davidson(const MatrixBase *diag, ntyp eigvalmax, ntyp eigvalmin)
{
    Davidson::diag = diag;
    n = diag->GetN();
    tab=false;

    ilow = eigvalmin;
    ihigh = eigvalmax;
    if (!ilow) ilow=ihigh;

    ntyp Nume;
    if ((ilow+ihigh-1)>n) {
        Nume=n-ilow+1;
    } else Nume=ihigh;
    SetDefault(Nume);

    iselec = new ntyp[limmax]; //mona by poprawic na liczbe poszukiwanych wart. wasnych
    iwork  = new ntyp[iiwsz];
    work   = new ftyp[irwsz];
}

Davidson::Davidson(const MatrixBase *diag, ntyp eigval[])
{
    Davidson::diag = diag;
    n = diag->GetN();
    tab=true;

    ilow = 0;
    ntyp eigvallow=eigval[0], eigvalhigh=eigval[0], Nume;
    for (ntyp i=1; eigval[i]>0; i++) {
    	eigvallow=MIN(eigvallow, eigval[i]);
        eigvalhigh=MAX(eigvalhigh, eigval[i]);
    }
    if ((eigvallow+eigvalhigh-1)>n) {
        Nume=n-eigvallow+1;
    } else Nume = eigvalhigh;
    SetDefault(Nume);

    iselec = eigval;   //new ntyp[limmax];
    iwork  = new ntyp[iiwsz];
    work   = new ftyp[irwsz];
}

void Davidson::SetEigSelect(ntyp eigval[])
{
	if (!tab) delete[] iselec;
    iselec =  eigval;

    ilow = 0;
    ntyp eigvallow=eigval[0], eigvalhigh=eigval[0], Nume;
    for (ntyp i=1; eigval[i]>0; i++) {
    	eigvallow=MIN(eigvallow, eigval[i]);
        eigvalhigh=MAX(eigvalhigh, eigval[i]);
    }
    if ((eigvallow+eigvalhigh-1)>n) {
        Nume=n-eigvallow+1;
    } else Nume = eigvalhigh;
    SetDefault(Nume);

    iwork = (ntyp *)realloc(iwork, sizeof(ntyp)*iiwsz);
    work  = (ftyp *)realloc(work, sizeof(ftyp)*irwsz);
    return;
}

Davidson::~Davidson()
{
    if (!tab) delete[] iselec;
    delete[] iwork;
    delete[] work;
}


void Davidson::SetDefault(ntyp Nume)
{
	Init=false;
    ierr=0;
    crite = 1e-15; critc = 1e-12; critr = 1e-8; ortho = 1e-9;
    niv = 0;
    mblock = 1;
    numemax = Nume;
    maxiter = MAX( numemax*40, 200 );
    nmv=0;
    neig = 0;

    limmax=numemax+20;
    iiwsz = 6*limmax + numemax;
    lim = limmax;

    irwsz = 2*n*limmax + limmax*limmax + (numemax+10)*limmax + numemax;

    hiend=false;

    norm=0;
    silence=false;
    loop=0;
    inf=0;
}

void Davidson::SetNumeMax(ntyp NumeMax)
{
    numemax=NumeMax;
    maxiter = MAX( numemax*40, 200 );
    ihigh = numemax;
    limmax=numemax+20;
    iiwsz = 6*limmax + numemax;
    lim = limmax;
    irwsz = 2*n*limmax + limmax*limmax + (numemax+10)*limmax + numemax;

    if (!tab) delete[] iselec;
    delete[] iwork;
    delete[] work;
    if (!tab) iselec = new ntyp[limmax];
    iwork  = new ntyp[iiwsz];
    work   = new ftyp[irwsz];
}

ntyp Davidson::GetNumeMax() const
{
	return numemax;
}

ntyp Davidson::GetLimMax() const
{
	return limmax;
}

void Davidson::SetLimMax(ntyp LimMax)
{
    limmax=LimMax;
    iiwsz = 6*limmax + numemax;
    lim = limmax;
    irwsz = 2*n*limmax + limmax*limmax + (numemax+10)*limmax + numemax;

    if (!tab) delete[] iselec;
    delete[] iwork;
    delete[] work;
    if (!tab) iselec = new ntyp[limmax];
    iwork  = new ntyp[iiwsz];
    work   = new ftyp[irwsz];
}

void Davidson::CheckError(ntyp iselec[])
{
	ntyp err=0;
    if (lim>n) ++err;
    if (lim<=0) err += 2;

    if ((ilow<=0)||(ilow>n)) {
    	//..Look for user choice of eigenpairs in ISELEC
    	if (iselec[0]<=0)
        	//..Nothing is given in ISELEC
        	err += 4;
        else {
        	/* ..Find number of eigenpairs wanted, and their            
               ..min/max indices    */
        	neig=1;
            ilow=iselec[0];
            ihigh=iselec[0];
            for (ntyp i=1; iselec[i]>0; i++) {
            	ilow=MIN(ilow, iselec[i]);
                ihigh=MAX(ihigh, iselec[i]);
                ++neig;
            }
            //..Check if a very large index is asked for
            if (ihigh>n) err += 8;
        }
    } else {
    	/* ..Look for a range between ILOW and IHIGH
           ..Invalid range. IHIGH>N */
    	if (ihigh>n) err += 8;
        neig=ihigh-ilow+1;
        //..Invalid range. IHIGH<ILOW
        if (neig<=0) err += 16;
        if (neig>lim)
        	//..Not enough Basis space. Increase LIM or decrease NEIG
        	err += 32;
        else
            //..Fill in the ISELEC with the required indices
        	for (ntyp i=1; i<=neig; i++) iselec[i-1]=ilow+i-1;
    }
    if (err) throw ErrDavidson(err);

    nume=ihigh;
    //..Identify if few of the highest eigenpairs are wanted.
    if ((ilow+ihigh-1)>n) {
    	hiend=true;
        nume=n-ilow+1;
    	/* ..Change the problem to a minimum eipenpairs one
           ..by picking the corresponding eigenpairs on the
           ..opposite side of the spectrum.    */
        for (ntyp i=1; i<=neig; i++) iselec[i-1]=n-iselec[i-1]+1;
    }

    if (neig>nume) err += 64;  //..duplications in ISELEC
    if ((nume>lim)||((nume==lim)&&(nume!=n))) err += 128; //..Not enough Basis space. Increase LIM or decrease NUME
    if ( (mblock<1)||(mblock>neig) ) err += 256; //..Size of Block out of bounds
    if ((irwsz<(lim*(2*n+lim+(nume+10))+nume)) ||
        (iiwsz<(6*lim+nume))) err += 512; //..Check for enough workspace for Dvdson
    if (err) throw ErrDavidson(err);

	if (niv>lim)
    	//..Check number of initial estimates NIV is lower than LIM.
    	cout<<"\n\nWARNING: Too many initial estimates. \n "
        	  "The routine will pick the appropriate number \n";
    else if (niv<nume && niv>0)
    	//..check if enough initial estimates.
        //..(NIV<1 => program chooses)
    	cout<<"\n\nWARNING: Not enough initial estimates \n"
        	  "The routine will pick the appropriate number \n";
}


void Davidson::Calculate()
{
    try {
        ExecuteDavidson();
    } catch (ErrDavidson e) {
        throw e;
    }
}

void Davidson::ExecuteDavidson()
{
    if (!tab)
    	for (ntyp i=0; i<limmax; i++) iselec[i]=0;
    //for (ntyp i=0; i<iiwsz; i++) iwork[i]=0;
    //for (ntyp i=0; i<irwsz; i++) work[i]=0.0;

    //set of initial vectors
    if (Init) InitVectors();

    try {
        CheckError(iselec);
    } catch (ErrDavidson e){
        throw e;
    }

    //Assigning space for the real work arrays
	ntyp ibasis    =1,
         ieigval   =ibasis  +n*lim,
       	 iab       =ieigval +lim,
     	 is        =iab     +n*lim,
    	 itemps    =is      +lim*(lim+1)/2,
    	 isvec     =itemps  +lim*(lim+1)/2,
    	 iscra1    =isvec   +lim*nume,
    	 ioldval   =iscra1  +8*lim,
    	 iscra2    =1,
    	 iscra3    =iscra2  +5*lim,
    	 iicv      =iscra3  +lim;

    if (hiend) dscal(n,-1.0,diag->Diagonal(),1);       
    ntyp istart=niv;
    setup(n,lim,nume,hiend,&work[iscra1-1], &work[ibasis-1],
          &work[iab-1],&work[is-1],&istart,
     	  diag);
    nmv=istart;
    dvdrvr(n,hiend,lim,mblock, nume,istart,neig,iselec,crite,critc,critr,
           ortho,maxiter,&work[ieigval-1],&work[ibasis-1],&work[iab-1],
           &work[is-1],&work[itemps-1],&work[isvec-1],&work[iscra1-1],
           &iwork[iscra2-1],&iwork[iscra3-1], &iwork[iicv-1],&work[ioldval-1],
     	   &nmv,&ierr, &loop, diag, silence, inf);
    if (hiend) {
    	dscal(n,-1.0,diag->Diagonal(),1);
        dscal(nume,-1.0,&work[ieigval-1],1);
    }
/*  -Copy the eigenvalues after the eigenvectors
*   -Next, copy the difference of eigenvalues between the last two steps  
*   -Next, copy the residuals for the first NUME estimates */
    dcopy(nume,&work[ieigval-1],1,&work[ibasis+n*nume-1],1);
    dcopy(nume,&work[ioldval-1],1,&work[ibasis+(n+1)*nume-1],1);
    dcopy(nume,&work[iscra1-1],1,&work[ibasis+(n+2)*nume-1],1);

    return;
}

ftyp Davidson::suma(ntyp n, ftyp tab[])
{
    ftyp suma=0;
    for(ntyp i=0; i<n; i++)
        suma+=tab[i]*tab[i];
    return suma;
}

ntyp Davidson::normalizacja(ntyp n, ftyp tab[])
{
    ftyp s=suma(n,tab);
    if (!s) return 0;
    ftyp wsp=1/suma(n,tab);

    for(ntyp i=0; i<n; i++)
        tab[i]*=wsp;

    return 1;
}


ftyp Davidson::EigenFunction(ntyp n, const GridParam &gp,
 							 ftyp *Psi_value, ntyp norm)
{
    ntyp j = 1 + (n-1)*gp.GetNmax();

    for (ntyp i=j; i<j+gp.GetNmax(); i++)
    	Psi_value[i-j]=work[i-1];
    if (norm){
    	Davidson::norm=normalizacja(gp.GetNmax(),Psi_value);
        sum=suma(gp.GetNmax(),Psi_value);
    } else {
    	Davidson::norm=0;
    }

    return work[nume*gp.GetNmax()+2*nume+n];
}

ftyp Davidson::EigenFunction(ntyp n, const GridParam &gp,
							 ftyp **Psi_value, ntyp norm)
{
    ftyp *psi = new ftyp[gp.GetNmax()];
    ntyp j = 1 + (n-1)*gp.GetNmax();

    for (ntyp i=j; i<j+gp.GetNmax(); i++)
    	psi[i-j]=work[i-1];
    if (norm){
    	Davidson::norm=normalizacja(gp.GetNmax(),psi);
        sum=suma(gp.GetNmax(),psi);
    } else {
    	Davidson::norm=0;
    }

    for (ntyp kk=0; kk<gp.GetNy(); kk++)
        for (ntyp jj=0; jj<gp.GetNx(); jj++)
        	Psi_value[jj][kk]=psi[gp.GetNy()*jj+kk];

    delete[] psi;
    return work[nume*gp.GetNmax()+2*nume+n];
}

ftyp Davidson::EigenFunction(ntyp n, const GridParam &gp,
 							 ftyp ***Psi_value, ntyp norm)
{
    ftyp *psi = new ftyp[gp.GetNmax()];
    ntyp j = 1 + (n-1)*gp.GetNmax();

    for (ntyp i=j; i<j+gp.GetNmax(); i++)
    	psi[i-j]=work[i-1];
    for (ntyp i=j; i<j+gp.GetNmax(); i++)
    	psi[i-j]=work[i-1];
    
    if (norm){
    	Davidson::norm=normalizacja(gp.GetNmax(),psi);
        sum=suma(gp.GetNmax(),psi);
    } else {
    	Davidson::norm=0;
    }

    for (ntyp ll=0; ll<gp.GetNz(); ll++)
        for (ntyp kk=0; kk<gp.GetNy(); kk++)
            for (ntyp jj=0; jj<gp.GetNx(); jj++)
                Psi_value[jj][kk][ll]=
                	psi[kk*gp.GetNz()+jj*gp.GetNy()*gp.GetNz()+ll];

    delete[] psi;
    return work[nume*gp.GetNmax()+2*nume+n];
}

void Davidson::EigenValue(eig_val &Eigen, ntyp n) const
{   
    Eigen.Eigenvalue=work[nume*Davidson::n+n-1];
    Eigen.Eigval_differences=work[nume*Davidson::n+nume+n-1];
    Eigen.Residual=work[nume*Davidson::n+2*nume+n-1];
}

ntyp Davidson::GetIerr() const
{
	return ierr;
}

bool Davidson::GetHiend() const
{
	return hiend;
}

ntyp Davidson::GetMatrixVectProd() const
{
	return nmv;
}

void Davidson::InitVectors()
{
    char filename[255];
	ifstream files(initfromfiles);
    unsigned int indx=0;
    ntyp noinitvectors=0;

    while (!files.eof() && noinitvectors<numemax) {
    	files.getline(filename,255);
        ifstream  init(filename);

        ntyp nosharp=0;
        char line[255];
        init.getline(line,255);
        while (line[0]=='#'){
        	++nosharp;
            init.getline(line,255);
        }
        init.close();
        init.open(filename,ios::in);
        for (ntyp i=0; i<nosharp; i++)
        	init.getline(line,255);

		for (ntyp i=0; i<n; i++) {
        	ftyp re,im;
            init>>re>>im;
            work[indx]=re;
            ++indx;
        }
        ++noinitvectors;
        init.close();
    }
    niv=noinitvectors;
}

void Davidson::ReadInitVectors(const char initfromfiles[])
{
	strcpy(Davidson::initfromfiles,initfromfiles);
    Init=true;
}

void Davidson::SetNiv(ntyp niv)
{
	Davidson::niv=niv;
}