matrixeigenstateproblem.cpp

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

#include <fstream>
#include <ctype.h>
#include "../davidson/matrixeigenstateproblem.h"

#define b(j,i) b[ (((i)*(n))+(j)) ]
#define c(j,i) c[ (((i)*(n))+(j)) ]

MatrixEigenstateProblem::MatrixEigenstateProblem(const GridParam &gp, const PotentialBase &V)
{
    MatrixEigenstateProblem::gp=&gp;
	diag = new ftyp[MatrixEigenstateProblem::gp->GetNmax()];
    switch (gp.GetDim()) {
        case 1  : Diag1D(V);
                  break;
        case 2  : Diag2D(V);
                  break;
        default : Diag3D(V);
                  break;  
    }
}

void MatrixEigenstateProblem::Diag1D(const PotentialBase &V)
{
    ftyp s2=-0.5,dx=gp->GetDx(),ax=s2/(dx*dx);
    ntyp nx=gp->GetNx();
    
    ax*=-2;
    for (ntyp j=1; j<=nx; j++) {
    	ftyp x=gp->GetXmin()+(j-1)*dx;
        diag[j -1]=ax+V(x);
    }
    return;
}

void MatrixEigenstateProblem::Diag2D(const PotentialBase &V)
{
    ftyp s2=-0.5, dx=gp->GetDx(), dy=gp->GetDy(), ax=s2/(dx*dx),
    	 ay=s2/(dy*dy), Axy=-2*(ax+ay);

    ntyp nx=gp->GetNx(),ny=gp->GetNy();
    for (ntyp k=1; k<=ny; k++){
        ftyp y = gp->GetYmin()+(k-1)*dy;
        for (ntyp j=1; j<=nx; j++){
            ftyp x = gp->GetXmin()+(j-1)*dx;
            diag[ny*(j-1)+k -1]=Axy+V(x, y);
        }
    }
    return;
}

void MatrixEigenstateProblem::Diag3D(const PotentialBase &V)
{
    ftyp s2=-0.5, dx=gp->GetDx(), dy=gp->GetDy(), dz=gp->GetDz(),
         ax=s2/(dx*dx), ay=s2/(dy*dy), az=s2/(dz*dz), Axyz=-2*(ax+ay+az);

    ntyp nx=gp->GetNx(), ny=gp->GetNy(), nz=gp->GetNz();
    for (ntyp l=1; l<=nz; l++){
        ftyp z = gp->GetZmin()+(l-1)*dz;
        for (ntyp k=1; k<=ny; k++){
            ftyp y = gp->GetYmin()+(k-1)*dy;
            for (ntyp j=1; j<=nx; j++){
                ftyp x = gp->GetXmin()+(j-1)*dx;
                diag[(k-1)*nz+(j-1)*ny*nz+l -1]=Axyz+V(x,y,z);
            }
        }
    }
    return;
}


MatrixEigenstateProblem::~MatrixEigenstateProblem()
{
	delete[] diag;
}

void MatrixEigenstateProblem::dssbmv(ntyp m, ftyp b[], ftyp c[], ntyp n) const
{
    switch (gp->GetDim()) {
        case 1  :   dssbmv1D(m, b, c, n);
                    break;
        case 2  :   dssbmv2D(m, b, c, n);
                    break;
        case 3  :   dssbmv3D(m, b, c, n);
                    break;
    }
	return;
}

void MatrixEigenstateProblem::dssbmv1D(ntyp m, ftyp b[], ftyp c[], ntyp n) const
{
    for(ntyp j=0; j<m; j++)
        for(ntyp i=0; i<n; i++)
            c(i,j)=diag[i]*b(i,j);

    ftyp s2=-0.5,
         ax=s2/(gp->GetDx()*gp->GetDx());
    ntyp temp1=n-1,
         temp2=temp1-1;

    for(ntyp k=0; k<m; k++){
        c(0,k)+=ax*b(1,k);
        c(temp1,k)+=ax*b(temp2,k);
        for(ntyp i=1; i<n-1; i++)
            c(i,k)+=ax*(b(i-1,k)+b(i+1,k));
    }
	return;
}

void MatrixEigenstateProblem::dssbmv2D(ntyp m, ftyp b[], ftyp c[], ntyp n) const
{
    for(ntyp j=0; j<m; j++)
        for(ntyp i=0; i<n; i++)
            c(i,j)=diag[i]*b(i,j);

    ftyp s2=-0.5,
         ax=s2/(gp->GetDx()*gp->GetDx()),
         ay=s2/(gp->GetDy()*gp->GetDy());
    ntyp ny=gp->GetNy(),
         temp1=n-1,
         temp2=temp1-1,
         temp3=n-1-ny;
    for(ntyp k=0; k<m; k++){
        c(0,k)+=ay*b(1,k)+ax*b(ny,k);
        c(temp1,k)+=ay*b(temp2,k)+ax*b(temp3,k);
        for(ntyp i=1; i<ny; i++)
            c(i,k)+=ax*b(i+ny,k)+ay*(b(i-1,k)+b(i+1,k));
        for(ntyp i=ny; i<n-1; i++){
            c(i,k)+=ay*(b(i-1,k)+b(i+1,k));
            for(ntyp j=i-ny; j<n; j++)
                if ((j==i-ny) || (j==i+ny)){
                    c(i,k)+=ax*b(j,k);
                    if (j==i+ny) break;
                    j=i+ny-1;
                }
        }
    }
	return;
}

void MatrixEigenstateProblem::dssbmv3D(ntyp m, ftyp b[], ftyp c[], ntyp n) const
{
    for(ntyp j=0; j<m; j++)
        for(ntyp i=0; i<n; i++)
            c(i,j)=diag[i]*b(i,j);

    ftyp s2=-0.5,
         ax=s2/(gp->GetDx()*gp->GetDx()),
         ay=s2/(gp->GetDy()*gp->GetDy()),
         az=s2/(gp->GetDz()*gp->GetDz());
    ntyp ny=gp->GetNy(),
         nz=gp->GetNz(),
         temp1=n-1,
         temp2=temp1-1,
         temp3=n-1-nz,
         temp4=n-1-ny*nz;
    for(ntyp k=0; k<m; k++){
        c(0,k)+=az*b(1,k)+ay*b(nz,k);
        c(temp1,k)+=az*b(temp2,k)+ay*b(temp3,k);
        for(ntyp i=1; i<nz; i++)
            c(i,k)+=ay*b(i+nz,k)+az*(b(i-1,k)+b(i+1,k));
        for(ntyp i=nz; i<n-1; i++){
            c(i,k)+=az*(b(i-1,k)+b(i+1,k));
            for(ntyp j=i-nz; j<n; j++)
                if ((j==i-nz) || (j==i+nz)){
                    c(i,k)+=ay*b(j,k);
                    if (j==i+nz) break;
                    j=i+nz-1;
                }
        }
    }

    ntyp NYZ=nz*ny;
    for(ntyp k=0; k<m; k++){
        c(0,k)+=ax*b(NYZ,k);
        c(temp1,k)+=ax*b(temp4,k);
        for(ntyp i=1; i<NYZ; i++)
            c(i,k)+=ax*b(i+NYZ,k);
        for(ntyp i=NYZ; i<n-1; i++)
            for(ntyp j=i-NYZ; j<n; j++)
                if ((j==i-NYZ) || (j==i+NYZ)){
                    c(i,k)+=ax*b(j,k);
                    if (j==i+NYZ) break;
                    j=i+NYZ-1;
                }
    }

	return;
}

ftyp *MatrixEigenstateProblem::Diagonal() const
{
	return diag;
}

ntyp MatrixEigenstateProblem::GetN() const
{
	return gp->GetNmax();
}

#undef b
#undef c