davidson_save.cpp

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

#include "../davidson/davidson_save.h"

#ifdef OS_WINDOWS

#include <stdlib.h>
#include <math.h>

#endif

typedef ftyp* pftyp;
typedef pftyp* ppftyp;

Davidson_save::Davidson_save(char katalog[], const ntyp &TypV, const ftyp &a,
							 const ftyp &V0, const ftyp &e0, const ftyp &w,
                             const GridParam &gp, Davidson &dav, ntyp norm)
{
    char pe[100];
    strcpy(pe,katalog);
    strncat(pe,"e.dat",5);

    ofstream e(pe);
    e<<"# ierr = "<<dav.GetIerr()<<" descending order: "<<dav.GetHiend()<<endl;
    if (TypV!=5) e<<"# a = "<<a<<"    "<<"V0 = "<<V0;
    	else e<<"#";
    if (TypV==4) e<<"    e0 = "<<e0<<"     "<<"w = "<<w;
    e<<endl<<"#     eigenvalues               eigval differences                 "
     <<"residuals"<<endl<<"#"<<endl;
    e.setf(ios::scientific);
    for (ntyp i=1; i<=dav.GetNumeMax(); i++){
    	eig_val eigenval;
        dav.EigenValue(eigenval,i);
     	e.width(20); e.fill(' '); e.precision(15);
            e<<eigenval.Eigenvalue<<"           "
            	<<eigenval.Eigval_differences<<"           "
                <<eigenval.Residual<<endl;
    }

    ntyp np=0;
    for ( ntyp j=1; j<=dav.GetNumeMax(); j++ ){
        char pwd[100], pwp[100], pwv[100];
        strcpy(pwd,katalog);
        strcpy(pwp,katalog);
        strcpy(pwv,katalog);

        ++np;

        char nw[10];
        itoa(np, nw, 10);

        strncat(pwd,"vector",6);    strncat(pwp,"vector",6);
        strncat(pwd,nw,strlen(nw)); strncat(pwp,nw,strlen(nw));
        strncat(pwd,".dat",4);      strncat(pwp,".plt",4);
        strncat(pwv,"v.dat",5);

        ofstream wektor(pwp);
        ofstream wektordat(pwd);

        if (dav.norm){
            wektor<<"# norm= "<<dav.norm<<"   ";
            wektordat<<"# norm= "<<dav.norm<<"   ";
            wektor<<"  sum= "<<dav.sum<<"   ";
        } else {
            wektor<<"# norm= "<<0<<"   ";
            wektordat<<"# norm= "<<0<<"   ";
        }
        wektor<<"n= "<<np<<endl;
        wektordat<<"n= "<<np<<endl;
        wektordat<<"#\t  re\t\t\t  im"<<endl;
        switch (gp.GetDim()) {
            case 1  :   SaveCor1D(gp, dav, j, wektor, wektordat, norm);
                        break;
            case 2  :   SaveCor2D(gp, dav, j, wektor, wektordat, norm);
                        break;
            case 3  :   SaveCor3D(gp, dav, j, wektor, wektordat, norm);
                        break;
        }

        wektor.close();
        wektordat.close();
    }
}


void Davidson_save::SaveCor1D(const GridParam &gp, Davidson &dav,
                         ntyp nrw, ofstream &wektor, ofstream &wektordat,
                         ntyp norm)
{
    ntyp nx=gp.GetNx();
    ftyp xmin=gp.GetXmin(),
         dx=gp.GetDx();

    wektor<<"# nx = "<<nx<<"     range = ["<<gp.GetXmax()<<"]"<<endl;
    wektor<<"#\t  x\t\t\t  psi(x)"<<endl;

    ftyp *psi = new ftyp[nx];
    dav.EigenFunction(nrw, gp, psi, norm);

    wektor.setf(ios::scientific);
    wektordat.setf(ios::scientific);
    wektor.width(20); wektor.fill(' '); wektor.precision(15);
    wektordat.width(20); wektordat.fill(' '); wektordat.precision(15);

    for (ntyp jj=0; jj<nx; jj++){
        ftyp x = xmin+jj*dx;
        wektor<<x<<"     "<<psi[jj]/sqrt(dx)<<endl;
        wektordat<<psi[jj]<<"       "<<0<<endl;
    }

    delete[] psi;

    wektordat<<gp.GetXmin()<<endl<<gp.GetXmax();
}

void Davidson_save::SaveCor2D(const GridParam &gp, Davidson &dav,
                         ntyp nrw, ofstream &wektor, ofstream &wektordat,
                         ntyp norm)
{
    ntyp nx=gp.GetNx(),
    	 ny=gp.GetNy();
    ftyp xmin=gp.GetXmin(),
         ymin=gp.GetYmin(),
         dx=gp.GetDx(),
         dy=gp.GetDy(),
         dxdy=sqrt(dx*dy);

    wektor<<"# nx = "<<nx<<"      ny = "<<ny
          <<"     range = ["<<gp.GetXmax()<<","<<gp.GetYmax()<<"]"<<endl;
    wektor<<"#\t   x\t\t\t\ty\t\t\tpsi(x,y)"<<endl;
    ftyp **psi = new pftyp[nx];
    for(int i=0; i<nx; i++) psi[i]=new ftyp[ny];
    dav.EigenFunction(nrw, gp, psi, norm);

    wektor.setf(ios::scientific);
    wektordat.setf(ios::scientific);
    wektor.width(20); wektor.fill(' '); wektor.precision(15);
    wektordat.width(20); wektordat.fill(' '); wektordat.precision(15);

    for (ntyp kk=0; kk<ny; kk++) {
    	ftyp y = ymin+kk*dy;
        for (ntyp jj=0; jj<nx; jj++){
            ftyp x = xmin+jj*dx;
            wektor<<x<<"     "<<y<<"     "<<psi[jj][kk]/dxdy<<endl;
            wektordat<<psi[jj][kk]<<"       "<<0<<endl;
        }
    }
    for(int i=0; i<nx; i++) delete[] psi[i];
    delete[] psi;

    wektordat<<gp.GetXmin()<<endl<<gp.GetXmax()<<endl;
    wektordat<<gp.GetYmin()<<endl<<gp.GetYmax();
}

void Davidson_save::SaveCor3D(const GridParam &gp, Davidson &dav,
                         ntyp nrw, ofstream &wektor, ofstream &wektordat,
                         ntyp norm)
{
    ntyp nx=gp.GetNx(),
         ny=gp.GetNy(),
         nz=gp.GetNz();
    ftyp xmin=gp.GetXmin(),
         ymin=gp.GetYmin(),
         zmin=gp.GetZmin(),
         dx=gp.GetDx(),
         dy=gp.GetDy(),
         dz=gp.GetDz(),
         dxdydz=sqrt(dx*dy*dz);

    wektor<<"# nx = "<<nx<<"      ny = "<<ny<<"      nz = "<<nz
          <<"     range = ["<<gp.GetXmax()<<","<<gp.GetYmax()<<","
          <<gp.GetZmax()<<"]"<<endl;
    wektor<<"#\t   x\t\t\t\ty\t\t\tz\t\t\tpsi(x,y,z)"<<endl;

    ftyp ***psi = new ppftyp[nx];
    for(int i=0; i<nx; i++) {
    	psi[i]=new pftyp[ny];
        for (int j=0; j<ny; j++) psi[i][j]=new ftyp[nz];
    }
    dav.EigenFunction(nrw, gp, psi, norm);

    wektor.setf(ios::scientific);
    wektordat.setf(ios::scientific);
    wektor.width(20); wektor.fill(' '); wektor.precision(15);
    wektordat.width(20); wektordat.fill(' '); wektordat.precision(15);

    for (ntyp l=0; l<nz; l++){
        ftyp z = zmin+l*dz;
        for (ntyp k=0; k<ny; k++){
            ftyp y = ymin+k*dy;
            for (ntyp j=0; j<nx; j++){
                ftyp x = xmin+j*dx;
                wektor<<x<<"     "<<y<<"     "<<z<<"     "<<psi[j][k][l]/dxdydz<<endl;
                wektordat<<psi[j][k][l]<<"       "<<0<<endl;
            }
        }
    }
    for(int i=0; i<nx; i++) {
    	for (int j=0; j<ny; j++) delete[] psi[i][j];
    	delete psi[i];
    }
    delete[] psi;

    wektordat<<gp.GetXmin()<<endl<<gp.GetXmax()<<endl;
    wektordat<<gp.GetYmin()<<endl<<gp.GetYmax()<<endl;
    wektordat<<gp.GetZmin()<<endl<<gp.GetZmax();
}