Commit 18b31a95 by Francois Gygi

rel1_33_1


git-svn-id: http://qboxcode.org/svn/qb/trunk@497 cba15fb0-1239-40c8-b417-11db7ca47a34
parent aea8b9f5
////////////////////////////////////////////////////////////////////////////////
//
// BasisMapping.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: BasisMapping.h,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#ifndef BASISMAPPING_H
#define BASISMAPPING_H
#include <complex>
#include <vector>
class Basis;
class Context;
class BasisMapping
{
private:
const Context& ctxt_;
const Basis& basis_;
int nprocs_, myproc_;
int np0_, np1_, np2_, np012_, np012loc_;
int nvec_;
std::vector<int> np2_loc_; // np2_loc_[iproc], iproc=0, nprocs_-1
std::vector<int> np2_first_; // np2_first_[iproc], iproc=0, nprocs_-1
std::vector<int> scounts, sdispl, rcounts, rdispl;
std::vector<std::complex<double> > sbuf, rbuf;
std::vector<int> ip_, im_;
std::vector<int> ipack_, iunpack_;
public:
BasisMapping (const Basis &basis);
int np0(void) const { return np0_; }
int np1(void) const { return np1_; }
int np2(void) const { return np2_; }
int np2loc(void) const { return np2_loc_[myproc_]; }
int np012(void) const { return np012_; }
int np012loc(void) const { return np012loc_; }
int nvec(void) const { return nvec_; }
int zvec_size(void) const { return nvec_ * np2_; }
const Context& context(void) const { return ctxt_; }
// map a function c(G) to zvec_
void vector_to_zvec(const std::complex<double> *c,
std::complex<double> *zvec);
// map zvec_ to a function c(G)
void zvec_to_vector(const std::complex<double> *zvec,
std::complex<double> *c);
void transpose_fwd(const std::complex<double> *zvec,
std::complex<double> *ct);
void transpose_bwd(const std::complex<double> *ct,
std::complex<double> *zvec);
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// ComputeMLWFCmd.C:
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ComputeMLWFCmd.C,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#include "ComputeMLWFCmd.h"
#include<iostream>
#include "Context.h"
#include "SlaterDet.h"
using namespace std;
int ComputeMLWFCmd::action(int argc, char **argv)
{
Wavefunction& wf = s->wf;
SlaterDet& sd = *(wf.sd(0,0));
mlwft = new MLWFTransform(sd);
mlwft->compute_transform();
mlwft->apply_transform(sd);
}
////////////////////////////////////////////////////////////////////////////////
//
// ComputeMLWFCmd.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: ComputeMLWFCmd.h,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#ifndef COMPUTEMLWFCMD_H
#define COMPUTEMLWFCMD_H
#include <iostream>
#include <stdlib.h>
#include <string>
#include "UserInterface.h"
#include "Sample.h"
#include "MLWFTransform.h"
class ComputeMLWFCmd : public Cmd
{
private:
MLWFTransform* mlwft;
public:
Sample *s;
ComputeMLWFCmd(Sample *sample) : s(sample) {};
char *name(void) const { return "compute_mlwf"; }
char *help_msg(void) const
{
return
"\n compute_mlwf\n\n"
" syntax: compute_mlwf\n\n"
" The compute_mlwf command computes maximally localized \n"
" Wannier functions.\n\n";
}
int action(int argc, char **argv);
ComputeMLWFCmd();
~ComputeMLWFCmd();
};
#endif
////////////////////////////////////////////////////////////////////////////////
//
// MLWFTransform.h
//
////////////////////////////////////////////////////////////////////////////////
// $Id: MLWFTransform.h,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#ifndef MLWFTRANSFORM_H
#define MLWFTRANSFORM_H
#include <vector>
#include <complex>
class SlaterDet;
class UnitCell;
class DoubleMatrix;
#include "D3vector.h"
#include "BasisMapping.h"
class MLWFTransform
{
private:
const SlaterDet& sd_;
const UnitCell& cell_;
const Context& ctxt_;
BasisMapping bm_;
std::vector<DoubleMatrix*> a_; // cosine and sine matrices
DoubleMatrix* u_; // orthogonal transformation
std::vector<std::vector<double> > adiag_; // diagonal elements
public:
void compute_transform(void);
void compute_sincos(const int n, const std::complex<double>* f,
std::complex<double>* fc, std::complex<double>* fs);
void apply_transform(SlaterDet& sd);
double spread2(int i, int j);
double spread2(int i);
double spread2(void);
double spread(int i);
double spread(void);
D3vector center(int i);
D3vector center(void);
MLWFTransform(const SlaterDet& sd);
~MLWFTransform(void);
};
#endif
This diff is collapsed. Click to expand it.
#include <vector>
int jade(int maxsweep, double tol, std::vector<DoubleMatrix*> a,
DoubleMatrix& u, std::vector<std::vector<double> > &adiag);
......@@ -23,6 +23,25 @@ Bug in scalapack pdgetri.f (inverse of square matrix)
Added bugfix from J.Langou into local modified copy of pdgetri.f
Must link to local pdgetri.o before the scalapack lib
--------------------------------------------------------------------------------
Fedora 6 on melodie: rebuilt blacs and scalapack. Use gfortran instead of g77
(g77 is not included in FC6)
--------------------------------------------------------------------------------
Check for potential problem in MDIonicStepper.C: random numbers used in the
Lowe and Andersen thermostats may differ on different nodes.
drand48() is used only in SlaterDet::randomize, where it is initialized
with srand48(ctxt.myproc()).
Could reinitialize the seed in MDIonicsStepper.C using a common value identical
on all processors.
--------------------------------------------------------------------------------
rel1_33_1
MLWFTransform.C: fixed bug in MLWF calc when cell sizes in x and y direction
differ.
--------------------------------------------------------------------------------
rel1_33_0 candidate
Implemented compute_mlwf command
Changed spline.C and sinft.C with new implementations.
Modified Species.C to use new spline and sinft interface.
--------------------------------------------------------------------------------
rel1_32_0
Dt.h, MDIonicStepper.C: modified to allow for negative time step.
Removed Sample.h dependency in WavefunctionStepper and derived classes.
......
......@@ -3,7 +3,7 @@
// qb.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: qb.C,v 1.52 2006-08-22 15:14:18 fgygi Exp $
// $Id: qb.C,v 1.53 2007-08-13 21:26:04 fgygi Exp $
#include <iostream>
#include <string>
......@@ -31,6 +31,7 @@ using namespace std;
#include "AngleCmd.h"
#include "AtomCmd.h"
#include "ComputeMLWFCmd.h"
#include "ConstraintCmd.h"
#include "DistanceCmd.h"
#include "HelpCmd.h"
......@@ -212,6 +213,7 @@ int main(int argc, char **argv, char **envp)
ui.addCmd(new AngleCmd(s));
ui.addCmd(new AtomCmd(s));
ui.addCmd(new ComputeMLWFCmd(s));
ui.addCmd(new ConstraintCmd(s));
ui.addCmd(new DistanceCmd(s));
ui.addCmd(new HelpCmd(s));
......
......@@ -3,10 +3,10 @@
// release.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: release.C,v 1.33 2007-03-17 01:14:00 fgygi Exp $
// $Id: release.C,v 1.34 2007-08-13 21:26:04 fgygi Exp $
#include "release.h"
std::string release(void)
{
return std::string("1.32.0");
return std::string("1.33.1c");
}
////////////////////////////////////////////////////////////////////////////////
//
// testBasisMapping.C
//
////////////////////////////////////////////////////////////////////////////////
// $Id: testBasisMapping.C,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#include "Context.h"
#include "Basis.h"
#include "UnitCell.h"
#include "BasisMapping.h"
#include "Timer.h"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cassert>
using namespace std;
#ifdef USE_MPI
#include <mpi.h>
#endif
int main(int argc, char **argv)
{
#if USE_MPI
MPI_Init(&argc,&argv);
#endif
{
// use: testBasisMapping a0 a1 a2 b0 b1 b2 c0 c1 c2 ecut npr npc
double err;
assert(argc==13);
D3vector a(atof(argv[1]),atof(argv[2]),atof(argv[3]));
D3vector b(atof(argv[4]),atof(argv[5]),atof(argv[6]));
D3vector c(atof(argv[7]),atof(argv[8]),atof(argv[9]));
UnitCell cell(a,b,c);
double ecut = atof(argv[10]);
D3vector kpoint;
int npr = atoi(argv[11]);
int npc = atoi(argv[12]);
Timer tm;
Context ctxt(npr,npc);
Basis basis(ctxt,kpoint);
basis.resize(cell,cell,ecut);
cout << " np0=" << basis.np(0)
<< " np1=" << basis.np(1)
<< " np2=" << basis.np(2) << endl;
cout << " basis.size=" << basis.size() << endl;
BasisMapping bmap(basis);
cout << " zvec_size=" << bmap.zvec_size() << endl;
cout << " np012loc=" << bmap.np012loc() << endl;
vector<complex<double> > zvec(bmap.zvec_size());
vector<complex<double> > ct(bmap.np012loc());
vector<complex<double> > f(basis.localsize());
for ( int i = 0; i < f.size(); i++ )
{
f[i] = exp(-basis.g2(i));
cout << ctxt.mype() << ": "
<< i << " " << basis.idx(3*i) << " "
<< basis.idx(3*i+1) << " " << basis.idx(3*i+2) << " "
<< f[i] << endl;
}
bmap.vector_to_zvec(&f[0],&zvec[0]);
bmap.transpose_fwd(&zvec[0],&ct[0]);
for ( int k = 0; k < bmap.np2loc(); k++ )
for ( int j = 0; j < bmap.np1(); j++ )
for ( int i = 0; i < bmap.np0(); i++ )
{
int index = i + bmap.np0() * ( j + bmap.np1() * k );
cout << ctxt.mype() << ": "
<< i << " " << j << " " << k << " " << ct[index] << endl;
}
// transpose back to zvec
for ( int i = 0; i < zvec.size(); i++ )
zvec[i] = 0.0;
bmap.transpose_bwd(&ct[0],&zvec[0]);
// transpose back to array f2
vector<complex<double> > f2(basis.localsize());
bmap.zvec_to_vector(&zvec[0],&f2[0]);
double sum = 0.0;
for ( int i = 0; i < f.size(); i++ )
{
sum += abs(f[i]-f2[i]);
}
cout << " total error: " << sum << endl;
}
#if USE_MPI
MPI_Finalize();
#endif
}
////////////////////////////////////////////////////////////////////////////////
//
// testjade.C
//
// use: testjade nprow npcol m mb
//
////////////////////////////////////////////////////////////////////////////////
// $Id: testjade.C,v 1.1 2007-08-13 21:26:27 fgygi Exp $
#include <cassert>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <string>
#include <fstream>
#include <valarray>
#include <algorithm>
using namespace std;
#include "Timer.h"
#ifdef USE_MPI
#include <mpi.h>
#endif
#include "Context.h"
#include "Matrix.h"
#include "jade.h"
const bool print = false;
enum MatrixType { FRANK, SCALED_FRANK, LAPLACIAN, SMALL };
const MatrixType mtype = SCALED_FRANK;
double aa(int i, int j) { return 1.0/(i+1)+2.0*i/(j+1); }
double bb(int i, int j) { return i-j-3; }
double frank(int n, int i, int j) { return n - max(i,j); }
double aijf(int n, int k, int i, int j)
{
switch ( mtype )
{
case FRANK :
if ( k == 0 )
return frank(n,i,j);
else
return frank(n,n-i,n-j);
case SCALED_FRANK :
if ( k == 0 )
return frank(n,i,j)/(n*n);
else
return frank(n,n-i,n-j)/(n*n);
case LAPLACIAN :
if ( i == j )
{
return 2.0 + k;
}
else if ( (i-j)*(i-j) == 1 )
{
return 1.0;
}
case SMALL :
{
if ( k == 0 )
{
if ( i == j && i == 0 )
return 1.0;
else if ( i == j && i > 0 )
return 4.0;
else
return 0.0;
}
else
{
if ( i == j && i == 0 )
return 1.0;
else if ( i == j && i > 0 )
return 1.0;
else
return 0.5;
}
}
}
return 0.0;
}
int main(int argc, char **argv)
{
// use: testjade nprow npcol n nb
const int maxsweep = 30;
const double tol = 1.e-8;
const int nmat = 2;
int mype;
int npes;
#ifdef USE_MPI
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD, &npes);
MPI_Comm_rank(MPI_COMM_WORLD, &mype);
#else
npes=1;
mype=0;
#endif
Timer tm;
assert(argc==5);
int nprow=atoi(argv[1]);
int npcol=atoi(argv[2]);
int m_a=atoi(argv[3]);
int mb_a=atoi(argv[4]);
int n_a = m_a;
int nb_a = mb_a;
if(mype == 0)
{
//infile >> nprow >> npcol;
cout << "nprow=" << nprow << ", npcol=" << npcol << endl;
//infile >> m_a >> n_a >> mb_a >> nb_a >> ta;
cout << "nmat=" << nmat << " m_a=" << m_a << ", n_a=" << n_a << endl;
cout << "tol=" << tol << endl;
}
#ifdef USE_MPI
MPI_Bcast(&nprow, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&npcol, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&m_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&n_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&mb_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Bcast(&nb_a, 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif
{
Context ctxt(nprow,npcol);
if ( mype == 0 )
{
cout << " Context " << ctxt.ictxt()
<< ": " << ctxt.nprow() << "x" << ctxt.npcol() << endl;
}
vector<DoubleMatrix*> a(nmat);
vector<vector<double> > adiag(nmat);
for ( int k = 0; k < nmat; k++ )
{
a[k] = new DoubleMatrix(ctxt,m_a,n_a,mb_a,nb_a);
adiag[k].resize(n_a);
}
DoubleMatrix u(ctxt,m_a,n_a,mb_a,nb_a);
if ( mype == 0 )
{
cout << " m_a x n_a / mb_a x nb_a / ta = "
<< a[0]->m() << "x" << a[0]->n() << " / "
<< a[0]->mb() << "x" << a[0]->nb() << endl;
}
for ( int k = 0; k < nmat; k++ )
{
for ( int m = 0; m < a[k]->nblocks(); m++ )
for ( int l = 0; l < a[k]->mblocks(); l++ )
for ( int y = 0; y < a[k]->nbs(m); y++ )
for ( int x = 0; x < a[k]->mbs(l); x++ )
{
int i = a[k]->i(l,x);
int j = a[k]->j(m,y);
//double aij = a.i(l,x) * 10 + a.j(m,y);
double aij = aijf(a[k]->n(),k,i,j);
int iii = x + l*a[k]->mb();
int jjj = y + m*a[k]->nb();
int ival = iii + jjj * a[k]->mloc();
(*a[k])[ival] = aij;
}
//cout << " a[" << k << "]=" << endl;
//cout << (*a[k]);
}
tm.start();
int nsweep = jade(maxsweep,tol,a,u,adiag);
tm.stop();
if ( ctxt.onpe0() )
{
cout << " m=" << m_a << " mb=" << mb_a
<< " ctxt: " << ctxt.nprow() << "x" << ctxt.npcol()
<< " nsweep=" << nsweep << " time: " << tm.real() << endl;
}
for ( int k = 0; k < nmat; k++ )
sort(adiag[k].begin(),adiag[k].end());
if ( nmat == 1 && (mtype == FRANK || mtype == SCALED_FRANK) )
{
vector<double> e_exact(adiag[0].size());
for ( int i = 0; i < e_exact.size(); i++ )
{
e_exact[i] = 1.0 / ( 2.0 * ( 1.0 - cos( ((2*i+1)*M_PI)/(2*n_a+1)) ) );
if ( mtype == SCALED_FRANK )
{
int n = a[0]->n();
e_exact[i] /= (n*n);
}
}
sort(e_exact.begin(),e_exact.end());
if ( mype == 0 )
{
for ( int k = 0; k < nmat; k++ )
{
double asum = 0.0;
for ( int i = 0; i < e_exact.size(); i++ )
{
//cout << ctxt.mype() << ": eig[" << k << "][" << i << "]= "
// << adiag[k][i]
// << " " << e_exact[i] << endl;
asum += fabs(adiag[k][i]-e_exact[i]);
}
cout << "a[" << k << "] sum of abs eigenvalue errors: "
<< asum << endl;
}
}
}
if ( print )
{
// a[k] contains AU.
// compute the product C = U^T A U
DoubleMatrix c(*a[0]);
for ( int k = 0; k < nmat; k++ )
{
c.gemm('t','n',1.0,u,(*a[k]),0.0);
if ( mype == 0 )
{
cout << " a[" << k << "]=" << endl;
cout << c;
}
}
cout << " u=" << endl;
cout << u;
}
}
#ifdef USE_MPI
MPI_Finalize();
#endif
}
......@@ -3,11 +3,11 @@
# x8664_gcc.mk
#
#-------------------------------------------------------------------------------
# $Id: x8664_gcc.mk,v 1.1 2006-07-21 17:56:50 fgygi Exp $
# $Id: x8664_gcc.mk,v 1.2 2007-08-13 21:26:04 fgygi Exp $
#
PLT=Linux_x8664
#-------------------------------------------------------------------------------
GCCDIR=/usr/lib/gcc-lib/x86_64-redhat-linux/3.2.3
#GCCDIR=/usr/lib/gcc/x86_64-redhat-linux/4.1.1
#MPIDIR=$(HOME)/software/mpich/mpich-1.2.6
MPIDIR=/opt/mpich-1.2.6
XERCESCDIR=$(HOME)/software/xml/Linux_x8664/xerces-c-src_2_5_0
......@@ -28,12 +28,12 @@
CXXFLAGS= -g -D$(PLT) $(INCLUDE) $(PLTFLAGS) $(DFLAGS)
LIBPATH = -L$(GCCDIR)/lib -L$(FFTWDIR)/.libs -L/usr/X11R6/lib \
-L$(MPIDIR)/lib -L$(BLASDIR) -L/usr/lib \
-L$(XERCESCDIR)/lib -L$(HOME)/lib
-L$(MPIDIR)/lib -L$(BLASDIR) \
-L$(XERCESCDIR)/lib
LIBS = $(PLIBS) $(GCCDIR)/libg2c.a -lfftw \
-llapack -lf77blas -latlas -lm -lmpich \
$(XERCESCDIR)/lib/libxerces-c.a
LIBS = $(PLIBS) -lfftw \
-llapack -lf77blas -latlas -lm -lmpich -lgfortran \
$(XERCESCDIR)/lib/libxerces-c.a
LDFLAGS = $(LIBPATH) $(LIBS)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment