diff --git a/heisen.h b/heisen.h new file mode 100755 index 0000000..b80fb4f --- /dev/null +++ b/heisen.h @@ -0,0 +1,376 @@ +/* + * heisen.h + * + * Created on: Aug 22, 2011 + * Author: wsttiger + */ + +#ifndef HEISEN_H_ +#define HEISEN_H_ + +#include +#include +#include +#include +#include + +#include "matrix.h" + +using std::cout; +using std::endl; +using std::map; +using std::pair; +using std::string; +using std::bitset; +using std::vector; +using std::complex; +using std::ifstream; + +#define NSIZE 4 + +typedef bitset ketT; + +//***************************************************************************** +void print_state(const ketT& ket) { + printf("|%s>\n",ket.to_string().c_str()); +} +//***************************************************************************** + +//***************************************************************************** +int cntbits(int s1, int s2, ketT state) { + if (s1 == s2) return 0; + + int sum = 0; + if (s1 > s2) { + for (int i = s2;i < s1; i++) { + if (state[i]) sum++; + } + } + else if (s1 < s2) { + for (int i = s1;i < s2; i++) { + if (state[i]) sum++; + } + } + return sum; +} +//***************************************************************************** + +//***************************************************************************** +class States { +public: + States () {} + + States(int nup, int nsites) + : nsites(nsites), nup(nup) { + + printf("states constructor ... begin\n"); + int tnstates = 1 << nsites; + + // create states for up and down spins + for (int i = 0; i < tnstates; i++) { + ketT tstate(i); + if (nup < 0) { + states.push_back(tstate); + } + else if (cntbits(0,nsites, tstate) == nup) { + states.push_back(tstate); + } + } + + for (unsigned int i = 0; i < states.size(); i++) { + umap.insert(pair(states[i].to_ulong(),i)); + } + + nst= states.size(); + printf("states constructor ... end\n"); + } + + int get_nst() const {return nst;} + + ketT get_state(int indx) const { + return states[indx]; + } + + unsigned long get_state_index(const ketT s) const { + return umap.find(s.to_ulong())->second; + } + + int get_nup() {return nup;} + + int get_ndown() {return nsites-nup;} + + void print() { + printf("nsites: %d nup: %d\n", nsites, nup); + for (unsigned int ist = 0; ist < states.size(); ist++) { + printf("%d |%s> %ld\n",ist, + states[ist].to_string().c_str(),states[ist].to_ulong()); + } + } + +private: + vector< ketT > states; + vector< ketT > dstates; + map umap; + + int nsites; + int nup; + int nst; + +}; +//***************************************************************************** + +//***************************************************************************** +struct Bond { + int i, j; + double J; + Bond() {} + Bond(int i, int j, double J) + : i(i), j(j), J(J) {} +}; +typedef vector latticeT; +//***************************************************************************** + +//***************************************************************************** +class Sector { +public: + Sector() {} + + Sector(latticeT lattice, int sites, int nup) + : lattice(lattice) { + states = States(nup, sites); + int nst = states.get_nst(); + } + + virtual vector apply_works_no_debug(const vector& svec) { + vector rvec(svec.size()); + for (auto ist = 0; ist < svec.size(); ist++) { + ketT st = states.get_state(ist); + double val = svec[ist]; + if (std::abs(val) > 1e-15) { + // loop over lattice sites + for (auto is = 0; is < lattice.size(); is++) { + Bond bond = lattice[is]; + int i = bond.i; int j = bond.j; double J = bond.J; + // diagonal + if (st[i] == st[j]) + rvec[ist] += -0.25*J*val; + // off-diagonal + else { + rvec[ist] += 0.25*J*val; + ketT stij = ketT(st); + stij[i] = !st[i]; stij[j] = !st[j]; + unsigned long lstij = states.get_state_index(stij); + rvec[lstij] -= 0.5*J*val; + } + } + } + } + return rvec; + } + + virtual vector apply_doesnt_work(const vector& svec) { + vector rvec(svec.size()); + for (auto ist = 0; ist < svec.size(); ist++) { + ketT st = states.get_state(ist); + double val = svec[ist]; + if (std::abs(val) > 1e-15) { + // loop over lattice sites + for (auto is = 0; is < lattice.size(); is++) { + Bond bond = lattice[is]; + int i = bond.i; int j = bond.j; double J = bond.J; + // off-diagonal + ketT spmop(0); + spmop[i] = true; spmop[j] = true; + int matchbits = cntbits(0,NSIZE,spmop & st); + if (matchbits == 1) { + rvec[ist] += 0.25*J*val; + ketT stij = ketT(st); + stij ^= spmop; + unsigned long lstij = states.get_state_index(stij); + rvec[lstij] -= 0.5*J*val; + } + else if (matchbits == 2) { + rvec[ist] += -0.25*J*val; + } + } + } + } + return rvec; + } + + virtual vector apply_works(const vector& svec) { + vector rvec(svec.size()); + for (auto ist = 0; ist < svec.size(); ist++) { + ketT st = states.get_state(ist); + double val = svec[ist]; + printf("apply(): ist = %d\n",ist); + if (std::abs(val) > 1e-15) { + // loop over lattice sites + for (auto is = 0; is < lattice.size(); is++) { + Bond bond = lattice[is]; + int i = bond.i; int j = bond.j; double J = bond.J; + // off-diagonal + printf("Bond (i,j): %d %d J = %4.2f\n",i,j,J); + if (st[i] != st[j]) { + ketT stij = ketT(st); + stij[i] = !st[i]; stij[j] = !st[j]; + unsigned long lstij = states.get_state_index(stij); + std::cout << lstij << "|" << st << "> -----> | " << stij << ">" << std::endl; + rvec[lstij] -= 0.5*J*val; + printf("rvec[%ld]: %5.3f\n", lstij, rvec[lstij]); + } + // diagonal + printf("J: %10.5f val: %10.5f\n", J, val); + std::cout << "(before) |" << st << "> --> " << rvec[ist] << std::endl; + rvec[ist] += (st[i] == st[j]) ? -0.25*J*val : 0.25*J*val; + printf("Diagonal: \n"); + std::cout << "(after) |" << st << "> --> " << rvec[ist] << std::endl; + printf("\n\n"); + } + } + } + return rvec; + } + + vector make_matrix() + { + int nst= states.get_nst(); + vector rmat(nst*nst,0.0); + + for (int i = 0; i < nst; i++) { + vector v1(nst,0.0); + v1[i] = 1.0; + for (int j = 0; j < nst; j++) { + vector v2(nst,0.0); + v2[j] = 1.0; + + vector rvec = apply_works_no_debug(v1); + double t1 = dotvec(rvec,v2); + rmat[i*nst+j] = t1; + } + } + + // make sure that the matrix is symmetric + double tol = 1e-10; + for (int i = 0; i < nst; i++) { + for (int j = 0; j < i; j++) { + if (fabs(rmat[i*nst+j]-rmat[j*nst+i]) > tol) + { + printf("ERROR: (make_matrix()) elements [%3d,%3d] != [%3d,%3d]\n", i,j,j,i); + exit(EXIT_FAILURE); + } + } + } + + return rmat; + } + + int get_nst() const {return states.get_nst();} +private: + States states; + latticeT lattice; +}; +//***************************************************************************** + +//***************************************************************************** +class HeisenCalculation { +public: + +// HeisenCalculation(const string& sfile) { +// ifstream astream(sfile.c_str()); +// int nup; +// double J_tmp; +// int nx; +// int ny; +// int latticetype; +// std::string lattice_file; +// +// while(!astream.eof()) { +// string s1; +// astream >> s1; +// s1.erase(remove_if(s1.begin(), s1.end(), isspace), s1.end()); +// if (s1.compare("") == 0) {} +// else if (s1.compare("J") == 0) { +// astream >> J_tmp; +// } +// else if (s1.compare("nup") == 0) { +// astream >> nup; +// } +// else if (s1.compare("latticetype") == 0) { +// astream >> latticetype; +// } +// else if (s1.compare("nx") == 0) { +// astream >> nx; +// } +// else if (s1.compare("ny") == 0) { +// astream >> ny; +// } +// else { +// printf("ERROR reading input file: %s\n\n", s1.c_str()); +// exit(EXIT_FAILURE); +// } +// } +// +// if (latticetype == 0) { +// if (ny == 1) +// { +// //create_1d_lattice(); +// } +// else { +// //create_2d_square_lattice(); +// } +// } +// else { +// printf("lattice file not implemented yet!\n\n"); +// assert(false); +// } +// +// states = States(nup,nx*ny); +// +// astream.close(); +// } + + HeisenCalculation(int sites) { + create_ring(sites); + } + + void create_ring(int sites) { +// double J = 1.0; +// for (int i = 0; i < sites-1; i++) { +// lattice.push_back(Bond(i,i+1,J)); +// } +// lattice.push_back(Bond(sites-1,0,J)); +// sectors.push_back(Sector(lattice,sites,-1)); + + double J = 1.0; + for (int i = 0; i < sites-1; i++) { + lattice.push_back(Bond(i,i+1,J)); + } + lattice.push_back(Bond(sites-1,0,J)); + for (int i = 0; i < sites; i++) { + sectors.push_back(Sector(lattice,sites,i)); + } + } + + vector eigenvalues() { + vector e; + for (unsigned int is = 0; is < sectors.size(); is++) { + Sector s = sectors[is]; + int nst = s.get_nst(); + vector mat = s.make_matrix(); + vector es(nst,0.0); + vector ev(nst*nst,0.0); + diag_matrix(mat,nst,es,ev); + std::copy(es.begin(),es.end(),std::back_inserter(e)); + } + std::sort(e.begin(),e.end(),[](const double& a, const double& b) {return a < b;}); + return e; + } + +private: + latticeT lattice; + std::vector sectors; +}; +//***************************************************************************** + +#endif diff --git a/heisen_chain.f90 b/heisen_chain.f90 new file mode 100644 index 0000000..ba9f5d1 --- /dev/null +++ b/heisen_chain.f90 @@ -0,0 +1,177 @@ +!-------------! + module system +!-------------! + + integer :: nn + integer :: nst + + real(8), allocatable :: mat(:,:) + real(8), allocatable :: vec(:,:) + real(8), allocatable :: enr(:) + real(8), allocatable :: spn(:) + real(8), allocatable :: mag(:) + + end module system +!-----------------! + +!================! + program hchain_0 +!================! + use system; implicit none + + open(10,file='read.in',status='old') + read(10,*)nn + close(10) + + nst=2**nn + allocate(mat(0:nst-1,0:nst-1)) + allocate(vec(0:nst-1,0:nst-1)) + allocate(enr(0:nst-1)) + allocate(spn(0:nst-1)) + allocate(mag(0:nst-1)) + + call hamiltonian() + call diagonalize(nst,mat,vec,enr) + + call spinsquared() + call transform(nst,mat,vec,spn) + spn(:)=0.5d0*abs(sqrt(1.d0+4.d0*spn(:))-1.d0) + + call magnetization() + + call writedata() + + deallocate(mat) + deallocate(vec) + deallocate(enr) + deallocate(spn) + deallocate(mag) + + end program hchain_0 +!====================! + +!----------------------! + subroutine writedata() +!----------------------! + use system; implicit none + + integer :: i + + open(10,file='eig.dat',status='unknown') + do i=0,nst-1 + write(10,'(i5,3f18.10)')i,enr(i),spn(i),mag(i) + enddo + close(10) + + end subroutine writedata +!------------------------! + +!------------------------! + subroutine hamiltonian() +!------------------------! + use system; implicit none + + integer :: i,j,a,b + + mat(:,:)=0.d0 + do a=0,nst-1 + do i=0,nn-1 + j=mod(i+1,nn) + if (btest(a,i).eqv.btest(a,j)) then + mat(a,a)=mat(a,a)+0.25d0 + else + mat(a,a)=mat(a,a)-0.25d0 + b=ieor(a,2**i+2**j) + mat(a,b)=mat(a,b)+0.5d0 + endif + enddo + enddo + + + end subroutine hamiltonian +!--------------------------! + +!------------------------! + subroutine spinsquared() +!------------------------! + use system; implicit none + + integer :: i,j,a,b,m + + mat(:,:)=0.d0 + do a=0,nst-1 + m=0 + do i=0,nn-1 + if (btest(a,i)) m=m+1 + enddo + mat(a,a)=dfloat(m-nn/2)**2+0.5d0*dfloat(nn) + do i=0,nn-1 + do j=i+1,nn-1 + if (btest(a,i).neqv.btest(a,j)) then + b=ieor(a,2**i+2**j) + mat(a,b)=mat(a,b)+1.d0 + endif + enddo + enddo + enddo + + end subroutine spinsquared +!--------------------------! + +!--------------------------! + subroutine magnetization() +!--------------------------! + use system; implicit none + + integer :: i,j,a,b + integer, allocatable :: mz(:) + + allocate(mz(0:nst-1)) + do a=0,nst-1 + mz(a)=0 + do i=0,nn-1 + if (btest(a,i)) mz(a)=mz(a)+1 + enddo + enddo + do a=0,nst-1 + mag(a)=0.d0 + do b=0,nst-1 + mag(a)=mag(a)+mz(b)*vec(b,a)**2 + enddo + enddo + mag(:)=(mag(:)-nn/2)*0.5d0 + deallocate(mz) + + end subroutine magnetization +!----------------------------! + +!-------------------------------------! + subroutine diagonalize(n,mat,vec,eig) +!-------------------------------------! + implicit none + + integer :: n,info + real(8) :: mat(n,n),vec(n,n),eig(n),work(n*(3+n/2)) + + vec=mat + call dsyev('V','U',n,vec,n,eig,work,n*(3+n/2),info) + + end subroutine diagonalize +!--------------------------! + +!-----------------------------------! + subroutine transform(n,mat,vec,dia) +!-----------------------------------! + implicit none + + integer :: i,n + real(8) :: mat(n,n),vec(n,n),dia(n) + + mat=matmul(mat,vec) + mat=matmul(transpose(vec),mat) + do i=1,n + dia(i)=mat(i,i) + enddo + + end subroutine transform +!------------------------! diff --git a/hubbard.h b/hubbard.h new file mode 100755 index 0000000..5cfca56 --- /dev/null +++ b/hubbard.h @@ -0,0 +1,1013 @@ +/* + * hubbard.h + * + * Created on: Aug 22, 2011 + * Author: wsttiger + */ + +#ifndef HUBBARD_H_ +#define HUBBARD_H_ + +#include +#include +#include +#include + +#include "matrix.h" +#include "lanczos.h" + +using std::cout; +using std::endl; +using std::map; +using std::pair; +using std::string; +using std::bitset; +using std::vector; +using std::complex; +using std::ifstream; + +#define NSIZE 26 + +//***************************************************************************** +int cntbits(int s1, int s2, bitset state) +{ + if (s1 == s2) return 0; + + int sum = 0; + if (s1 > s2) + { + for (int i = s2;i < s1; i++) + { + if (state[i]) sum++; + } + } + else if (s1 < s2) + { + for (int i = s1;i < s2; i++) + { + if (state[i]) sum++; + } + } + return sum; +} +//***************************************************************************** + +//***************************************************************************** +double total_U_for_state(const bitset& state1, const bitset& state2, + const std::vector U) +{ + bitset state = state1 & state2; + double totalU = 0.0; + for (unsigned int i = 0; i < U.size(); i++) + { + if (state[i]) totalU += U[i]; + } + return totalU; +} +//***************************************************************************** + +//***************************************************************************** +class States +{ +public: + States(int nup, int ndown, int nsites) + { + this->nup = nup; + this->ndown = ndown; + + int tnstates = 1 << nsites; + int usum = 0; + int dsum = 0; + nup2 = nchoosek(nsites,nup); + ndown2 = nchoosek(nsites,ndown); + + // create states for up and down spins + for (int i = 0; i < tnstates; i++) + { + bitset tstate(i); + if (cntbits(0,nsites, tstate) == nup) + { + ustates.push_back(tstate); + usum++; + } + if (cntbits(0,nsites,tstate) == ndown) + { + dstates.push_back(tstate); + dsum++; + } + if (usum > nup2) + { + cout << "[[Warning]] create_states: usum > nup2" << endl; + cout << " usum: " << usum << " nup2: " << nup2 << endl; + exit(EXIT_FAILURE); + } + if (dsum > ndown2) + { + cout << "[[Warning]] create_states: dsum > ndown2" << endl; + cout << " dsum: " << usum << " ndown2: " << nup2 << endl; + exit(EXIT_FAILURE); + } + } + + // create map where the first item is the state of a single spin state + // i.e. 00000000000000000000000011101000 and give the associated index + for (unsigned int i = 0; i < ustates.size(); i++) + { + umap.insert(pair(ustates[i].to_ulong(),i)); + } + for (unsigned int i = 0; i < dstates.size(); i++) + { + dmap.insert(pair(dstates[i].to_ulong(),i)); + } + + nstates = ustates.size()*dstates.size(); + } + + int get_nstates() {return nstates;} + + bitset get_ustate(int indx) + { + return ustates[indx]; + } + + bitset get_dstate(int indx) + { + return dstates[indx]; + } + + unsigned long get_ustate_index(const bitset us) + { + return umap.find(us.to_ulong())->second; + } + + unsigned long get_dstate_index(const bitset ds) + { + return dmap.find(ds.to_ulong())->second; + } + + int get_nup() {return nup;} + + int get_ndown() {return ndown;} + + int get_nup2() {return nup2;} + + int get_ndown2() {return ndown2;} + + void print_states2() + { + printf("\nStates: %d up states %d down states\n",nup,ndown); + for (unsigned int ist = 0; ist < ustates.size()*dstates.size(); ist++) + { + int ust = ist / ndown2; + int dst = ist % ndown2; + printf("%d %d %d [%s , %s]\n",ist,ust,dst, + ustates[ust].to_string().c_str(),dstates[dst].to_string().c_str()); + } + } + + +private: + vector< bitset > ustates; + vector< bitset > dstates; + map umap; + map dmap; + + int nup; + int ndown; + int nup2; + int ndown2; + int nstates; + +}; +//***************************************************************************** + +////***************************************************************************** +//class IModel +//{ +//public: +// virtual vector apply(const vector& v) = 0; +// virtual vector apply(const vector& v, States* s) = 0; +// virtual int get_nstates() const = 0; +//}; +////***************************************************************************** + +//***************************************************************************** +//class HubbardCalculation : public IModel +class HubbardCalculation +{ +private: + int nx; + int ny; +// double U; +// double t; + States* states; + double tol; + bool bdebug; + + // Ground state vector + vector gsv; + + struct neighbor_2d + { + int me; + int n1; + + neighbor_2d(int me_, int n1_) + { + me = me_; + n1 = n1_; + } + }; + + vector lattice; + vector t; + vector U; + +public: + + HubbardCalculation(const string& sfile) + { + ifstream astream(sfile.c_str()); + int nup; + int ndown; + double t_tmp; + double U_tmp; + double nx_tmp; + double ny_tmp; + int latticetype; + std::string lattice_file; + + while(!astream.eof()) + { + string s1; + astream >> s1; + s1.erase(remove_if(s1.begin(), s1.end(), isspace), s1.end()); + if (s1.compare("") == 0) {} + else if (s1.compare("t") == 0) + { + astream >> t_tmp; + } + else if (s1.compare("U") == 0) + { + astream >> U_tmp; + } + else if (s1.compare("nup") == 0) + { + astream >> nup; + } + else if (s1.compare("ndown") == 0) + { + astream >> ndown; + } + else if (s1.compare("latticetype") == 0) + { + astream >> latticetype; + } + else if (s1.compare("nx") == 0) + { + astream >> nx_tmp; + } + else if (s1.compare("ny") == 0) + { + astream >> ny_tmp; + } +// else if (s1.compare("") == 0) +// { +// astream >> ny_tmp; +// } + else + { + printf("ERROR reading input file: %s\n\n", s1.c_str()); + exit(EXIT_FAILURE); + } + } + + if (latticetype == 0) + { + nx = nx_tmp; + ny = ny_tmp; + + if (ny == 1) + { + create_1d_lattice(); + } + else + { + create_2d_square_lattice(); + } + + for (unsigned int il = 0; il < lattice.size(); il++) + { + t.push_back(t_tmp); + } + for (unsigned int isite = 0; isite < nx*ny; isite++) + { + U.push_back(U_tmp); + } + } + else + { + printf("lattice file not implemented yet!\n\n"); + assert(false); + } + + tol = 1e-8; + bdebug = false; + + states = new States(nup,ndown,nx*ny); + //states->print_states2(); + + astream.close(); + } + + HubbardCalculation(int nx, int ny, int nup, int ndown, double U_, double t_ = 1.0) + : nx(nx), ny(ny) + { + if (ny == 1) + { + create_1d_lattice(); + } + else + { + create_2d_square_lattice(); + } + tol = 1e-8; + bdebug = false; + + states = new States(nup,ndown,nx*ny); + + } + + ~HubbardCalculation() + { + delete states; + } + + virtual States* get_states() + { + return states; + } + + virtual int get_nstates() const + { + return states->get_nstates(); + } + +// void print_states() +// { +// printf("\nStates: %d up states %d down states\n",nup,ndown); +// int i = 0; +// for (unsigned int ust = 0; ust < ustates.size(); ust++) +// { +// for (unsigned int dst = 0; dst < dstates.size(); dst++,i++) +// { +// printf("%5d [%s , %s]\n",i,ustates[ust].to_string().c_str(),dstates[dst].to_string().c_str()); +// } +// } +// } + +// void print_states2() +// { +// printf("\nStates: %d up states %d down states\n",nup,ndown); +// for (unsigned int ist = 0; ist < ustates.size()*dstates.size(); ist++) +// { +// int ust = ist / ndown2; +// int dst = ist % ndown2; +// printf("%d %d [%s , %s]\n",ust,dst, +// ustates[ust].to_string().c_str(),dstates[dst].to_string().c_str()); +// } +// } + + void print_state(const vector& vec) + { + printf("\nState:\n"); + for (unsigned int ist = 0; ist < states->get_nstates(); ist++) + { + unsigned int ndown2 = states->get_ndown2(); + unsigned int ust = ist / ndown2; + unsigned int dst = ist % ndown2; + printf("%10.5f [%s , %s]\n",vec[ist], + states->get_ustate(ust).to_string().c_str(), + states->get_dstate(ust).to_string().c_str()); + } + } + + void print_lattice() + { + printf("\nLattice: \n"); + for (unsigned int i = 0; i < lattice.size(); i++) + { + int me = lattice[i].me; + int n1 = lattice[i].n1; + printf("%10d%10d\n", me, n1); + } + } + + void print_matrix_element(int il1, int il2) + { + int nstates = states->get_nstates(); + vector v1(nstates,0.0); + v1[il1] = 1.0; + vector v2(nstates,0.0); + v2[il2] = 1.0; + vector rvec = apply(v2); + double me = dotvec(rvec,v1); + printf("Matrix element for [%5d,%5d] = %7.3f\n",il1,il2,me); + } + + vector make_matrix() + { + return make_matrix(this->states); + } + + vector make_matrix(States* s) + { + int nstates = s->get_nstates(); + vector rmat(nstates*nstates,0.0); + + for (int i = 0; i < nstates; i++) + { + vector v1(nstates,0.0); + v1[i] = 1.0; + for (int j = 0; j < nstates; j++) + { + vector v2(nstates,0.0); + v2[j] = 1.0; + + // DEBUG + //printf("V1: "); + //for (int iv = 0; iv < v1.size(); iv++) printf("%5.3f ",v1[iv]); + //printf("\n\n"); + //printf("V2: "); + //for (int iv = 0; iv < v2.size(); iv++) printf("%5.3f ",v2[iv]); + //printf("\n\n"); + + vector rvec = apply(v1); + double t1 = dotvec(rvec,v2); + rmat[i*nstates+j] = t1; + } + } + + // make sure that the matrix is symmetric + for (int i = 0; i < nstates; i++) + { + for (int j = 0; j < i; j++) + { + printf("%15.8f%15.8f\n",rmat[i*nstates+j],rmat[j*nstates+i]); + if (fabs(rmat[i*nstates+j]-rmat[j*nstates+i]) > tol) + { + printf("ERROR: (make_matrix()) elements [%3d,%3d] != [%3d,%3d]\n", i,j,j,i); + exit(EXIT_FAILURE); + } + } + } + + return rmat; + } + + void print_apply_U(bitset us, bitset ds, double uval) + { + printf("%s %s U = %10.3f\n", us.to_string().c_str(), + ds.to_string().c_str(),uval); + } + + void print_apply_thop(bitset bs1, bitset bs2, int me, int n1, double tval) + { + printf("%s ----> %s [%4d%4d] t = %10.3f\n", bs1.to_string().c_str(), + bs2.to_string().c_str(),me,n1,tval); + } + + void create_2d_square_lattice() + { + // loop over lattice sites (delegate to lattice object?) + for (int si = 0; si < nx; si++) + { + for (int sj = 0; sj < ny; sj++) + { + // get index + int me = si*ny+sj; + + // the (0,-1) case + int si2 = si; + int sj2 = (sj==0) ? (ny-1) : sj-1; + int n1 = si2*ny+sj2; + lattice.push_back(neighbor_2d(me,n1)); + + // the (0,+1) case + si2 = si; + sj2 = (sj==(ny-1)) ? 0 : sj+1; + int n2 = si2*ny+sj2; + lattice.push_back(neighbor_2d(me,n2)); + + // the (-1,0) case + si2 = (si==0) ? (nx-1) : si-1; + sj2 = sj; + int n3 = si2*ny+sj2; + lattice.push_back(neighbor_2d(me,n3)); + + // the (+1,0) case + si2 = (si==(nx-1)) ? 0 : si+1; + sj2 = sj; + int n4 = si2*ny+sj2; + lattice.push_back(neighbor_2d(me,n4)); + } + } + } + + void create_1d_lattice() + { + // loop over lattice sites (delegate to lattice object?) + for (int si = 0; si < nx; si++) + { + // the -1 case + int si2 = (si==0) ? (nx-1) : si-1; + lattice.push_back(neighbor_2d(si,si2)); + + // the +1 case + si2 = (si==(nx-1)) ? 0 : si+1; + lattice.push_back(neighbor_2d(si,si2)); + } + } + + virtual vector apply_w_openmp(const vector& svec1) + { + vector svec2(svec1.size(),0.0); + + // loop of coeffs of state + int nsites = nx*ny; + int ndown2 = states->get_ndown2(); + #pragma omp parallel + { + vector _svec2(svec1.size(),0.0); + #pragma omp for + for (unsigned ist = 0; ist < svec1.size(); ist++) + { + int ust = ist / ndown2; + int dst = ist % ndown2; + bitset us = states->get_ustate(ust); + bitset ds = states->get_dstate(dst); + string ustr = us.to_string(); + string dstr = ds.to_string(); + + double val1 = svec1[ist]; + // double occupancies + if (fabs(val1) > tol) + { + double totalU = total_U_for_state(us, ds, U); + _svec2[ust*ndown2+dst] += totalU*val1; + } + + for (unsigned int li = 0; li < lattice.size(); li++) + { + neighbor_2d nghbr = lattice[li]; + int me = nghbr.me; + int n1 = nghbr.n1; + + // up spin + // destination state (do hopping) + if ((fabs(val1) > tol) && !us[me] && us[n1]) + { + bitset us2 = thop(me,n1,us); + string us2str = us2.to_string(); + // get the index of the destination state + unsigned long du1 = states->get_ustate_index(us2); + // sign + double ut1 = -t[li]*computephase(me,n1,us); + // write destination + _svec2[du1*ndown2+dst] += ut1*val1; + } + + // down spin + // destination state (do hopping) + if ((fabs(val1) > tol) && !ds[me] && ds[n1]) + { + bitset ds2 = thop(me,n1,ds); + string ds2str = ds2.to_string(); + // get the index of the destination state + unsigned long dd1 = states->get_dstate_index(ds2); + // sign + double dt1 = -t[li]*computephase(me,n1,ds); + // write destination + _svec2[ust*ndown2+dd1] += dt1*val1; + } + } + } + #pragma omp critical (update_sum) + { + for (unsigned int i = 0; i < svec2.size(); i++) svec2[i] += _svec2[i]; + } + } + return svec2; + } + + virtual vector apply(const vector& svec1) + { + return apply(svec1, states); + } + + virtual vector apply(const vector& svec1, States* s) + { + // output vector + vector svec2(svec1.size(),0.0); + + // loop of coeffs of state + int ndown2 = s->get_ndown2(); + for (unsigned int ist = 0; ist < svec1.size(); ist++) + { + // get up and down spins portions + int ust = ist / ndown2; + int dst = ist % ndown2; + bitset us = s->get_ustate(ust); + bitset ds = s->get_dstate(dst); + string ustr = us.to_string(); + string dstr = ds.to_string(); + + double val1 = svec1[ist]; + // double occupancies + if (fabs(val1) > tol) + { + double totalU = total_U_for_state(us, ds, U); + svec2[ust*ndown2+dst] += totalU*val1; + } + + for (unsigned int li = 0; li < lattice.size(); li++) + { + neighbor_2d nghbr = lattice[li]; + int me = nghbr.me; + int n1 = nghbr.n1; + + // up spin + // destination state (do hopping) + if ((fabs(val1) > tol) && !us[me] && us[n1]) + { + bitset us2 = thop(me,n1,us); + string us2str = us2.to_string(); + + // get the index of the destination state + unsigned long du1 = s->get_ustate_index(us2); + // sign + double ut1 = -t[li]*computephase(me,n1,us); + // write destination + svec2[du1*ndown2+dst] += ut1*val1; + + // DEBUG + //int ist2 = du1*ndown2+dst; + //if ((ist == 0 && ist2 == 2) || (ist == 2) && (ist2 == 0)) + //{ + // printf("[[[ ist = %d ist2 = %d ]]]\n",ist,ist2); + // printf("*****************************************************\n"); + // cout << "ustr: " << ustr << endl; + // cout << "dstr: " << dstr << endl; + // cout << endl; + + // printf("\nme = %d n1 = %d\n",me,n1); + // cout << "ustr2: " << us2str << endl; + // printf("*****************************************************\n\n"); + //} + + } + + // down spin + // destination state (do hopping) + if ((fabs(val1) > tol) && !ds[me] && ds[n1]) + { + bitset ds2 = thop(me,n1,ds); + string ds2str = ds2.to_string(); + + // get the index of the destination state + unsigned long dd1 = s->get_dstate_index(ds2); + // sign + double dt1 = -t[li]*computephase(me,n1,ds); + // write destination + svec2[ust*ndown2+dd1] += dt1*val1; + + // DEBUG + //int ist2 = dd1*ndown2+dst; + //if ((ist == 0 && ist2 == 2) || (ist == 2) && (ist2 == 0)) + //{ + // printf("[[[ ist = %d ist2 = %d ]]]\n",ist,ist2); + // printf("*****************************************************\n"); + // cout << "ustr: " << ustr << endl; + // cout << "dstr: " << dstr << endl; + // cout << endl; + + // printf("\nme = %d n1 = %d\n",me,n1); + // cout << "dstr2: " << ds2str << endl; + // printf("*****************************************************\n\n"); + //} + + } + } + } + return svec2; + } + + void compute_1p_greens_function_matrix(double omega0, + double omega1, int nomega, double eta) + { + // ground state + int lsize = 100; + Lanczos l(this,lsize); + l.run(); + vector gs = l.lowstate(); + + // create enlarged Hilbert space (up electrons) + int nsites = nx*ny; + States s(states->get_nup()+1,states->get_ndown(), nsites); + + // energy points + vector< complex > gf1p(nomega,0.0); + double domega = (omega1-omega0)/nomega; + + // create Hamiltonian matrix (with enlarged Hilbert space) + vector H = make_matrix(&s); + + // We're going to reexpress ground state in a larger Hilbert space + // where we've created a particle at site isite. + for (int isite = 0; isite < 1; isite++) + { + vector gs2(s.get_nstates(),0.0); + int ndwn = states->get_ndown2(); + for (unsigned int ist = 0; ist < gs.size(); ist++) + { + // if gs[ist] is not big enough to bother + if (fabs(gs[ist]) > tol) + { + int ust = ist / ndwn; + int dst = ist % ndwn; + // need to make copy so we use the copy constructor + bitset us(states->get_ustate(ust)); + bitset ds(states->get_dstate(dst)); + // get string for debugging purposes + string ustr = us.to_string(); + string dstr = ds.to_string(); + + // Now, we'ed like to create a particle at isite + // Is there already an electron at isite with up spin? + if (!us[isite]) + { + // create "up" particle at isite, and translate this + // "new" state to gs2 + us[isite] = true; + unsigned long ust2 = s.get_ustate_index(us); + unsigned int ist2 = ust2*ndwn+dst; + gs2[ist2] = gs[ist]; + } + } + } + + // loop over energy grid + for (int iw = 0; iw < nomega; iw++) + { + // create matrix + int nstates = s.get_nstates(); + vector< complex > zmat(nstates*nstates,0.0); + double omega = omega0 + iw * domega; + for (int iz = 0; iz < nstates; iz++) + { + zmat[iz*nstates+iz] = complex(omega-H[iz*nstates+iz],eta); + for (int jz = 0; jz < iz; jz++) + { + zmat[iz*nstates+jz] = complex(-H[iz*nstates+jz],0.0); + } + } + + // invert matrix + vector< complex > izmat = invert_matrix(zmat,nstates); + // get first element + complex t1 = izmat[0]; + double t2 = 1/(double)nsites; + gf1p[iw] += t1 * t2; + } + } + + for (int iw = 0; iw < nomega; iw++) + { + double omega = omega0 + iw * domega; + char tmpstr[256]; + sprintf(tmpstr, "%15.8f %15.8f %15.8f",omega,real(gf1p[iw]), imag(gf1p[iw])); + cout << tmpstr << std::endl; + } + + } + + void compute_1p_greens_function(double omega0, + double omega1, int nomega, double eta) + { + // ground state + int lsize = 100; + Lanczos l(this,lsize); + l.run(); + vector gs = l.lowstate(); + + // create enlarged Hilbert space (up electrons) + int nsites = nx*ny; + States s(states->get_nup()+1,states->get_ndown(), nsites); + + // energy points + vector< complex > gf1p(nomega,0.0); + double domega = (omega1-omega0)/nomega; + + // We're going to reexpress ground state in a larger Hilbert space + // where we've created a particle at site isite. + for (int isite = 0; isite < 1; isite++) + { + vector gs2(s.get_nstates(),0.0); + int ndwn = states->get_ndown2(); + for (unsigned int ist = 0; ist < gs.size(); ist++) + { + // if gs[ist] is not big enough to bother + if (fabs(gs[ist]) > tol) + { + int ust = ist / ndwn; + int dst = ist % ndwn; + // need to make copy so we use the copy constructor + bitset us(states->get_ustate(ust)); + bitset ds(states->get_dstate(dst)); + // get string for debugging purposes + string ustr = us.to_string(); + string dstr = ds.to_string(); + + // Now, we'ed like to create a particle at isite + // Is there already an electron at isite with up spin? + if (!us[isite]) + { + // create "up" particle at isite, and translate this + // "new" state to gs2 + us[isite] = true; + unsigned long ust2 = s.get_ustate_index(us); + unsigned int ist2 = ust2*ndwn+dst; + gs2[ist2] = gs[ist]; + } + } + } + + // swap enlarged Hilbert space and smaller + States* ststmp = states; + states = &s; + + // we want to construct an approximation for H (i.e. the diagonal + // and off diagonal coefficients) in the Lanczos basis, but we + // want to apply to our N+1 ground state or other + // words c_i^{+} | GS > + Lanczos l2(this,lsize,gs2); + l2.run(); + + // loop over energy grid + for (int iw = 0; iw < nomega; iw++) + { + // create matrix + vector< complex > zmat(lsize*lsize,0.0); + double omega = omega0 + iw * domega; + for (int iz = 0; iz < lsize; iz++) + { + zmat[iz*lsize+iz] = complex(omega-l._a[iz],eta); + if (iz < lsize-1) + { + zmat[(iz+1)*lsize+iz] = complex(-l._b[iz],0.0); + zmat[iz*lsize+iz+1] = complex(-l._b[iz],0.0); + } + } + // invert matrix + vector< complex > izmat = invert_matrix(zmat,lsize); + // get first element + complex t1 = izmat[0]; + double t2 = 1/(double)nsites; + gf1p[iw] += t1 * t2; + } + + // return states + states = ststmp; + } + + std::ofstream fstr("gf1p.txt"); + for (int iw = 0; iw < nomega; iw++) + { + double omega = omega0 + iw * domega; + char tmpstr[256]; + sprintf(tmpstr, "%15.8f %15.8f %15.8f",omega,real(gf1p[iw]), imag(gf1p[iw])); + cout << tmpstr << std::endl; + } + fstr.close(); + + } + + bitset thop(const int& s1, const int& s2, const bitset& state) + { + // tij --> i == s1 and j == s2 + if (!state[s1] && state[s2]) + { + bitset state2(state); + state2.flip(s1); + state2.flip(s2); + return state2; + } + } + + double computephase(int s1, int s2, bitset state) + { + if (s1 == s2) return 0.0; + + int sum = 0; + if (s1 > s2) + { + for (int i = s2+1;i < s1; i++) + { + if (state[i]) sum++; + } + } + else if (s1 < s2) + { + for (int i = s1+1;i < s2; i++) + { + if (state[i]) sum++; + } + } + return ((sum % 2) == 0) ? 1.0 : -1.0; + + } + +// vector test_apply(const vector& v) +// { +// int n = v.size(); +// vector A(n*n,0.0); +// vector v2(n,0.0); +//// for (int j = 0; j < n; j++) +//// { +//// for (int i = 0; i < j+1; i++) +//// { +//// A[i*n+j] = i*n+j; +//// A[j*n+i] = i*n+j; +//// } +//// } +// A = make_matrix(); +// for (int i = 0; i < n; i++) +// { +// for (int j = 0; j < n; j++) +// { +// v2[i] += A[i*n+j] * v[j]; +// } +// } +// +// return v2; +// } + + void lanczos(int n = 100) + { + // create random initial vector + int nstates = states->get_nstates(); + vector v = random_vector(nstates); + +// // save initial vector +// for (int i = 0; i < nstates; i++) +// { +// vinit[i] = v[i]; +// } + + vector vold(nstates,0.0); + vector a(n,0.0); + vector b(n-1,0.0); + vector v2(nstates,0.0); + for (int i = 0; i < n; i++) + { + if (i > 0) + { + v2 = wst_daxpy(1.0,apply(v),-b[i-1],vold); + } + else + { + v2 = apply(v); + } + a[i] = dotvec(v,v2); + if (i < (n-1)) + { + wst_daxpy_inplace(1.0,v2,-a[i],v); + b[i] = norm2(v2); + for (int iv = 0; iv < nstates; iv++) + { + vold[iv] = v[iv]; + v[iv] = v2[iv]/b[i]; + } + } + } + + vector mat(n*n,0.0); + for (int i = 0; i < n-1; i++) + { + mat[i*n+i] = a[i]; + mat[i*n+i+1] = b[i]; + mat[(i+1)*n+i] = b[i]; + } + mat[n*n-1] = a[n-1]; + + vector e(n,0.0); + vector ev(n*n, 0.0); + diag_matrix(mat,n,e,ev); + printf("Lanczos: lowest eigenvalue is %15.8f\n\n", e[0]); + } +}; +//***************************************************************************** + + + +#endif /* HUBBARD_H_ */ diff --git a/hubbard.in b/hubbard.in new file mode 100755 index 0000000..68afbbd --- /dev/null +++ b/hubbard.in @@ -0,0 +1,7 @@ +t 1.0 +U 5.0 +nup 4 +ndown 3 +latticetype 0 +nx 3 +ny 3 diff --git a/lanczos.h b/lanczos.h new file mode 100755 index 0000000..1055d8b --- /dev/null +++ b/lanczos.h @@ -0,0 +1,290 @@ +/* + * lanczos.h + * + * Created on: Aug 22, 2011 + * Author: wsttiger + */ + +#ifndef LANCZOS_H_ +#define LANCZOS_H_ + +#include +#include +#include +#include +#include + +#include "matrix.h" + +using std::cout; +using std::endl; +using std::map; +using std::pair; +using std::string; +using std::bitset; +using std::vector; +using std::complex; +using std::ifstream; + +static double ttt, sss; +void START_TIMER() { +// ttt=wall_time(); sss=cpu_time(); +} + +void END_TIMER(const char* msg) { +// ttt=wall_time()-ttt; sss=cpu_time()-sss; printf("timer: %20.20s %8.2fs %8.2fs\n", msg, sss, ttt); +} + +//***************************************************************************** +struct MatrixLanczos +{ + vector _H; + + unsigned int _nstates; // is the bigger dimension (full size) + unsigned int _nsize; // is the size of the approximate matrix in the Lanczos basis + + vector _vinit; + vector _a; + vector _b; + vector _e; + vector _ev; + + MatrixLanczos(const vector& H, const int& nsize) + { + _H = vector(H); + _nstates = int(floor(sqrt(H.size()))); + printf("MatrixLanczos constructor: H.size() = %5d _nstates = %5d\n\n", int(H.size()), _nstates); + _nsize = nsize; + create_init_random_vector(); + } + + MatrixLanczos(const vector& H, const int& nsize, const vector& vinit) + { + _H = vector(H); + _nstates = int(floor(sqrt(H.size()))); + printf("MatrixLanczos constructor: H.size() = %5d _nstates = %5d\n\n", int(H.size()), _nstates); + _nsize = nsize; + assert(vinit.size() == _nstates); + _vinit = vector(vinit); + } + + void create_init_random_vector() + { + _vinit = random_vector(_nstates); + } + + // run() computes the _a and _b parameters (matrix diagonal + // and off-diagonal) in the Lanczos basis from an initial vector + // Initial vector can be generated or provided + void run() + { + vector v(_nstates,0.0); + for (unsigned int i = 0; i < _nstates; i++) + { + v[i] = _vinit[i]; + } + vector vold(_nstates,0.0); + _a = vector(_nsize,0.0); + _b = vector(_nsize-1,0.0); + vector v2(_nstates,0.0); + for (unsigned int i = 0; i < _nsize; i++) + { + if (i > 0) + { + v2 = wst_daxpy(1.0,matrix_vector_mult(_H,v),-_b[i-1],vold); + } + else + { + v2 = matrix_vector_mult(_H,v); + } + _a[i] = dotvec(v,v2); + if (i < (_nsize-1)) + { + wst_daxpy_inplace(1.0,v2,-_a[i],v); + _b[i] = norm2(v2); + for (unsigned int iv = 0; iv < _nstates; iv++) + { + vold[iv] = v[iv]; + v[iv] = v2[iv]/_b[i]; + } + } + } + + vector mat(_nsize*_nsize,0.0); + for (unsigned int i = 0; i < _nsize-1; i++) + { + mat[i*_nsize+i] = _a[i]; + mat[i*_nsize+i+1] = _b[i]; + mat[(i+1)*_nsize+i] = _b[i]; + } + mat[_nsize*_nsize-1] = _a[_nsize-1]; + + _e = vector(_nsize,0.0); + _ev = vector(_nsize*_nsize, 0.0); + diag_matrix(mat,_nsize,_e,_ev); + printf("Lanczos: lowest eigenvalue is %15.8f\n\n", _e[0]); + } + +}; +//***************************************************************************** + +//***************************************************************************** +template +struct Lanczos +{ + T* _model; + unsigned int _nstates; // is the bigger dimension (full size) + unsigned int _nsize; // is the size of the approximate matrix in the Lanczos basis + + vector _vinit; + vector _a; + vector _b; + vector _e; + vector _ev; + + + + Lanczos(T* model, int nsize) + { + _model = model; + _nstates = _model->get_nstates(); + _nsize = nsize; + create_init_random_vector(); + } + + Lanczos(T* model, int nsize, const vector& vinit) + { + _model = model; + _nstates = _model->get_nstates(); + _nsize = nsize; + assert(vinit.size() == _nstates); + _vinit = vector(_nstates,0.0); + for (unsigned int i = 0; i < _nstates; i++) + { + _vinit[i] = vinit[i]; + } + } + + void create_init_random_vector() + { + _vinit = vector(_nstates,0.0); + for (unsigned int i = 0; i < _nstates; i++) + { + int i1 = rand(); + double t1 = (i1 % 100000)/100000.0; + _vinit[i] = t1; + } + normalize(_vinit); + } + + // run() computes the _a and _b parameters (matrix diagonal + // and off-diagonal) in the Lanczos basis from an initial vector + // Initial vector can be generated or provided + void run() + { + vector v(_nstates,0.0); + for (unsigned int i = 0; i < _nstates; i++) + { + v[i] = _vinit[i]; + } + vector vold(_nstates,0.0); + _a = vector(_nsize,0.0); + _b = vector(_nsize-1,0.0); + vector v2(_nstates,0.0); + for (unsigned int i = 0; i < _nsize; i++) + { + START_TIMER(); + if (i > 0) + { + v2 = wst_daxpy(1.0,_model->apply(v),-_b[i-1],vold); + } + else + { + v2 = _model->apply(v); + } + END_TIMER("lanczos run 1"); + START_TIMER(); + _a[i] = dotvec(v,v2); + END_TIMER("lanczos run 2"); + if (i < (_nsize-1)) + START_TIMER(); + { + wst_daxpy_inplace(1.0,v2,-_a[i],v); + _b[i] = norm2(v2); + for (unsigned int iv = 0; iv < _nstates; iv++) + { + vold[iv] = v[iv]; + v[iv] = v2[iv]/_b[i]; + } + } + END_TIMER("lanczos run 3"); + } + + vector mat(_nsize*_nsize,0.0); + for (unsigned int i = 0; i < _nsize-1; i++) + { + mat[i*_nsize+i] = _a[i]; + mat[i*_nsize+i+1] = _b[i]; + mat[(i+1)*_nsize+i] = _b[i]; + } + mat[_nsize*_nsize-1] = _a[_nsize-1]; + + _e = vector(_nsize,0.0); + _ev = vector(_nsize*_nsize, 0.0); + diag_matrix(mat,_nsize,_e,_ev); + printf("Lanczos: lowest eigenvalue is %15.8f\n\n", _e[0]); + } + + // run(vector) computes the lowest eigenstate using the _a and _b + // parameters that were already generated from a previous run + vector run(const vector& vin) + { + // create return vector + vector rv(_nstates,0.0); + // copy initial vector + vector v(_nstates,0.0); + for (unsigned int i = 0; i < _nstates; i++) + { + v[i] = vin[i]; + } + vector vold(_nstates,0.0); + vector v2(_nstates,0.0); + for (unsigned int i = 0; i < _nsize; i++) + { + // build up return value + //wst_daxpy_inplace(1.0,rv,_ev[i*_nsize],v); + for (unsigned int iv = 0; iv < _nstates; iv++) + { + rv[iv] += _ev[i*_nsize]*v[iv]; + } + if (i > 0) + { + v2 = wst_daxpy(1.0,_model->apply(v),-_b[i-1],vold); + } + else + { + v2 = _model->apply(v); + } + if (i < (_nsize-1)) + { + wst_daxpy_inplace(1.0,v2,-_a[i],v); + for (unsigned int iv = 0; iv < _nstates; iv++) + { + vold[iv] = v[iv]; + v[iv] = v2[iv]/_b[i]; + } + } + } + normalize(rv); + return rv; + } + + vector lowstate() + { + return run(_vinit); + } + +}; +//***************************************************************************** + +#endif diff --git a/lanczos.m b/lanczos.m new file mode 100755 index 0000000..1d6ac4a --- /dev/null +++ b/lanczos.m @@ -0,0 +1,56 @@ +function [B,gsv] = lanczos(A,n) + +m = size(A,1); +v = rand(m,1); +v = v / norm(v); +vold = zeros(m,1); + +% save initial v +vinit = v; + +a = zeros(n,1); +b = zeros(n-1,1); +Vsaved = zeros(m,n); +for i = 1:n + Vsaved(:,i) = v; + v2 = A*v; + if (i > 1) + v2 = v2 - b(i-1)*vold; + end + a(i) = v'*v2; + v2 = v2 - a(i)*v; + if (i < n) + b(i) = norm(v2); + vold = v; + v = v2 / b(i); + end +end + +B = diag(a) + diag(b,1) + diag(b,-1); + +[ev,e] = eig(B); +gsv = Vsaved*ev(:,1); + +v = vinit; +gsv = zeros(m,1); +% a = zeros(n,1); +% b = zeros(n-1,1); +for i = 1:n + % sum for ground state + gsv = gsv + ev(i,1)*v; + v2 = A*v; + if (i > 1) + v2 = v2 - b(i-1)*vold; + end +% a(i) = v'*v2; + v2 = v2 - a(i)*v; + if (i < n) +% b(i) = norm(v2); + vold = v; + v = v2 / b(i); + end +end + +gsv = gsv / norm(gsv); + + diff --git a/main.cc b/main.cc new file mode 100755 index 0000000..4a4678e --- /dev/null +++ b/main.cc @@ -0,0 +1,169 @@ +/* + * main.cc + * + * Created on: Apr 20, 2011 + * Author: wsttiger + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "matrix.h" +#include "hubbard.h" + +using std::cout; +using std::endl; +using std::bitset; +using std::vector; +using std::map; +using std::pair; +using std::string; +using std::complex; + +using namespace Eigen; + +#define NSIZE 26 + + + +//void test_eigen3() +//{ +// MatrixXcf A = MatrixXcf::Random(4,4); +// cout << "Here is a random 4x4 matrix, A:" << endl << A << endl << endl; +// +// ComplexEigenSolver ces; +// ces.compute(A); +// VectorXcf ev = ces.eigenvalues(); +// MatrixXcf evv = ces.eigenvectors(); +// +// printf("WST: the lowest eigenvalue: %15.7f\n\n", real(ev[0])); +// +// cout << "The eigenvalues of A are:" << endl << ces.eigenvalues() << endl; +// cout << "The matrix of eigenvectors, V, is:" << endl << ces.eigenvectors() << endl << endl; +// +// complex lambda = ces.eigenvalues()[0]; +// cout << "Consider the first eigenvalue, lambda = " << lambda << endl; +// VectorXcf v = ces.eigenvectors().col(0); +// cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl; +// cout << "... and A * v = " << endl << A * v << endl << endl; +// +// cout << "Finally, V * D * V^(-1) = " << endl +// << ces.eigenvectors() * ces.eigenvalues().asDiagonal() * ces.eigenvectors().inverse() << endl; +//} + +//void test_diag() +//{ +// int n = 3; +// vector mat(n,0.0); +// +// for (int i = 0; i < n; i++) +// { +// for (int j = 0; j < n; j++) +// { +// mat[i*n+j] = i*n+j; +// } +// } +// +// print_matrix(mat,n,n); +// vector e(n,0.0); +// vector ev(n*n,0.0); +//} + +void hubbard_w_matrix() +{ + // create HubbardCalculation object + HubbardCalculation hc(3,2,3,3,8.0,1.0); + printf("number of states: %d\n\n", hc.get_nstates()); + + int nstates = hc.get_nstates(); + vector hmat = hc.make_matrix(hc.get_states()); + vector e(nstates,0.0); + vector ev(nstates*nstates,0.0); + diag_matrix(hmat,nstates,e,ev); + printf("diag matrix: %15.8f\n\n", e[0]); +} + +void hubbard_w_lanczos() +{ + // create HubbardCalculation object + HubbardCalculation hc(3,2,3,3,8.0,1.0); + printf("number of states: %d\n\n", hc.get_nstates()); + + Lanczos l(&hc,100); + l.run(); + + vector ev1 = l.lowstate(); + //print_vector(ev1); + hc.compute_1p_greens_function(-10.0,10.0,200,0.1); +} + +void hubbard_w_lanczos_from_file() +{ + // create HubbardCalculation object + HubbardCalculation hc("hubbard.in"); + printf("number of states: %d\n\n", hc.get_nstates()); + + //vector mat = hc.make_matrix(); + //print_matrix(mat,hc.get_nstates(), hc.get_nstates()); + //exit(EXIT_FAILURE); + + Lanczos l(&hc,100); + l.run(); + + hc.make_matrix(); + +// vector ev1 = l.lowstate(); +// //print_vector(ev1); +//// hc.compute_1p_greens_function(-32.0,16.0,800,0.1); +// hc.compute_1p_greens_function(-3.0,6.0,8,0.1); +// hc.compute_1p_greens_function_matrix(-3.0,6.0,8,0.1); +} + +void test_hubbard_w_openmp() +{ + // create HubbardCalculation object + HubbardCalculation hc(4,4,6,6,8.0,1.0); + unsigned int nstates = hc.get_nstates(); + printf("number of states: %d\n\n", nstates); + + // create random vector + vector v1 = random_vector(nstates); + // apply hubbard hamiltonian to v1 without openmp + vector rv1 = hc.apply(v1); + + unsigned int ntimes = 200; + for (unsigned int i = 0; i < ntimes; i++) + { + vector rv2 = hc.apply_w_openmp(v1); + string result = (is_equals(rv1,rv2)) ? "PASS!" : "FAIL!"; + cout << "trial: " << i << " " << result << endl; + } +} + +int main(int argc, char** argv) +{ + hubbard_w_lanczos_from_file(); + +// int tid; +// printf("Hello world from threads:\n"); +// #pragma omp parallel private (tid) +// { +// tid = omp_get_thread_num(); +// printf("Hello from thread: <%d>\n", tid); +// } +// printf("I am sequential now.\n"); +// +// test_hubbard_w_openmp(); +// return 0; + + return 0; +} + diff --git a/matrix.h b/matrix.h new file mode 100755 index 0000000..a33bfd7 --- /dev/null +++ b/matrix.h @@ -0,0 +1,368 @@ +/* + * matrix.h + * + * Created on: Aug 22, 2011 + * Author: wsttiger + */ + +#include +#include + +using std::vector; +using std::complex; + +#ifndef MATRIX_H_ +#define MATRIX_H_ + +extern "C" void dsyev_(char *jobz, char *uplo, int *n, double *a, + int *lda, double *w, double *work, int *lwork, + int *info); +extern "C" void dgetri(int* n, double* a, int *lda, int *ipiv, + double* work, int* lwork, int* info); +extern "C" void dgetrf(int* m, int* n, double* a, int* lda, int* ipiv, int* info); +extern "C" void zgetri(int* n, complex* a, int *lda, int *ipiv, + complex* work, int* lwork, int* info); +extern "C" void zgetrf(int* m, int* n, complex* a, int* lda, int* ipiv, int* info); + +//***************************************************************************** +long fact(int n) +{ + long prod = 1; + for (long i = 1; i <= n; i++) + { + prod *= i; + } + return prod; +} +//***************************************************************************** + +//***************************************************************************** +int nchoosek(int n, int k) +{ + long rval = 1; + for (long i = k+1; i <= n; i++) + { + rval *= i; + } + return (rval/fact(n-k)); +} +//***************************************************************************** + +//***************************************************************************** +double dotvec(const vector& v1, const vector& v2) +{ + double rval = 0.0; + for (unsigned int iv = 0; iv < v1.size(); iv++) + { + rval += v1[iv]*v2[iv]; + } + return rval; +} +//***************************************************************************** + +//r = a*x + b*y +//***************************************************************************** +vector wst_daxpy(const double& a, const vector& x, + const double& b, const vector& y) +{ + vector rv(x.size(),0.0); + for (unsigned int i = 0; i < rv.size(); i++) + { + rv[i] = a*x[i]+b*y[i]; + } + return rv; +} +//***************************************************************************** + +// in-place version of daxpy +// x = a*x + b*y +//***************************************************************************** +void wst_daxpy_inplace(const double& a, vector& x, + const double& b, const vector& y) +{ + for (unsigned int i = 0; i < x.size(); i++) + { + x[i] = a*x[i]+b*y[i]; + } +} +//***************************************************************************** + +//***************************************************************************** +double norm2(const vector& v) +{ + double rnorm = 0.0; + for (unsigned int iv = 0; iv < v.size(); iv++) + { + rnorm += v[iv]*v[iv]; + } + return sqrt(rnorm); +} +//***************************************************************************** + +//***************************************************************************** +vector matrix_vector_mult(vector mat, vector v) +{ + unsigned int szmat = mat.size(); + unsigned int sz2 = v.size(); + unsigned int sz1 = szmat / sz1; + vector rv(sz2,0.0); + for (unsigned int i = 0; i < sz2; i++) + { + for (unsigned int j = 0; j < sz1; j++) + { + rv[i] += mat[i*sz1+j]; + } + } + return rv; +} +//***************************************************************************** + +//***************************************************************************** +void normalize(vector& v) +{ + double rnorm = 0.0; + for (unsigned int iv = 0; iv < v.size(); iv++) + { + rnorm += v[iv]*v[iv]; + } + rnorm = sqrt(rnorm); + for (unsigned int iv = 0; iv < v.size(); iv++) + { + v[iv] /= rnorm; + } +} +//***************************************************************************** + +//***************************************************************************** +void print_vector(const vector& v) +{ + for (unsigned int i = 0; i < v.size(); i++) + { + printf(" %10.5f\n",v[i]); + } +} +//***************************************************************************** + +//***************************************************************************** +void print_matrix(const vector& mat, int m, int n) +{ + for (int i = 0; i < m; i++) + { + for (int j = 0; j < n; j++) + { + printf("%8.2f",mat[i*n+j]); + } + printf("\n"); + } +} +//***************************************************************************** + +//***************************************************************************** +void print_matrix(const vector< complex >& mat, int m, int n) +{ + for (int i = 0; i < m; i++) + { + for (int j = 0; j < n; j++) + { + printf(" (%8.4f,%8.4f)",mat[i*n+j].real(),mat[i*n+j].imag()); + } + printf("\n"); + } +} +//***************************************************************************** + +//void diag_matrix2(const vector& mat, int n, +// vector& e, vector & evec) +//{ +// MatrixXcf A = MatrixXcf::Zero(n,n); +// for (int i = 0; i < n; i++) +// { +// for (int j = 0; j < n; j++) +// { +// A(i,j) = complex(mat[i*n+j],0.0); +// } +// } +// +// ComplexEigenSolver ces; +// ces.compute(A); +// VectorXcf et = ces.eigenvalues(); +// MatrixXcf evt = ces.eigenvectors(); +// +// for (int i = 0; i < n; i++) +// { +// e[i] = real(et(i)); +// for (int j = 0; j < n; j++) +// { +// evec[i*n+j] = real(evt(i,j)); +// } +// } +// +//} + +//***************************************************************************** +vector get_col_from_matrix(const vector mat, int col, int m, int n) +{ + vector rv(m,0.0); + for (int i = 0; i < m; i++) + { + rv[i] = mat[i*n+col]; + } + return rv; +} +//***************************************************************************** + +//***************************************************************************** +void diag_matrix(const vector& mat, int n, + vector& e, vector & evec) +{ + char jobz = 'V'; + char uplo = 'U'; + int info; + double* a = new double[n*n]; + int lda = n; + int lwork = 3*n-1; + double *work = new double[lwork]; + double *et = new double[n]; + for (int i = 0; i < n*n; i++) a[i] = mat[i]; + + dsyev_(&jobz, &uplo, &n, a, &lda, et, work, &lwork, &info); + + if (info != 0) + { + printf("[[Error:]] lapack::dsyev failed --- info = %d\n\n", info); + } + else + { + //printf("Eigenvalues:\n"); + for (int i = 0; i < n; i++) + { + e[i] = et[i]; + } + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + evec[j*n+i] = a[i*n+j]; + } + } + +} + delete a; + delete work; + delete et; +} +//***************************************************************************** + +//***************************************************************************** +vector invert_matrix(const vector& mat, int n) +{ + vector rv(n*n,0.0); + int info; + double* a = new double[n*n]; + int lda = n; + int lwork = 3*n-1; + int* ipiv = new int[n]; + for (int i = 0; i < n; i++) ipiv[i] = 0; + double *work = new double[lwork]; + for (int i = 0; i < n*n; i++) a[i] = mat[i]; + + dgetrf(&n, &n, a, &lda, ipiv, &info); + if (info != 0) + { + printf("[[Error:]] lapack::dgetrf failed --- info = %d\n\n", info); + } + dgetri(&n, a, &lda, ipiv, work, &lwork, &info); + + if (info != 0) + { + printf("[[Error:]] lapack::dgetri failed --- info = %d\n\n", info); + } + else + { + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + rv[i*n+j] = a[i*n+j]; + } + } + + } + delete a; + delete work; + delete ipiv; + + return rv; +} +//***************************************************************************** + +//***************************************************************************** +vector< complex > invert_matrix(const vector< complex >& mat, int n) +{ + vector< complex > rv(n*n,0.0); + int info; + complex* a = new complex[n*n]; + int lda = n; + int lwork = 3*n-1; + int* ipiv = new int[n]; + for (int i = 0; i < n; i++) ipiv[i] = 0; + complex* work = new complex[lwork]; + for (int i = 0; i < n*n; i++) a[i] = mat[i]; + + zgetrf(&n, &n, a, &lda, ipiv, &info); + if (info != 0) + { + printf("[[Error:]] lapack::zgetrf failed --- info = %d\n\n", info); + } + zgetri(&n, a, &lda, ipiv, work, &lwork, &info); + + if (info != 0) + { + printf("[[Error:]] lapack::zgetri failed --- info = %d\n\n", info); + } + else + { + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + rv[i*n+j] = a[i*n+j]; + } + } + + } + delete a; + delete work; + delete ipiv; + + return rv; +} +//***************************************************************************** + +//***************************************************************************** +vector random_vector(const int& nsize) +{ + vector v(nsize,0.0); + for (int i = 0; i < nsize; i++) + { + int i1 = rand(); + double t1 = (i1 % 100000)/100000.0; + v[i] = t1; + } + normalize(v); + return v; +} +//***************************************************************************** + +//***************************************************************************** +bool is_equals(const vector& a, const vector& b, double tol = 1e-10) +{ + for (unsigned int i = 0; i < a.size(); i++) + { + if (fabs(a[i]-b[i]) > tol) return false; + } + return true; +} +//***************************************************************************** + +#endif /* MATRIX_H_ */ diff --git a/read.in b/read.in new file mode 100644 index 0000000..b8626c4 --- /dev/null +++ b/read.in @@ -0,0 +1 @@ +4 diff --git a/test_heisen.cc b/test_heisen.cc new file mode 100644 index 0000000..7eb0aeb --- /dev/null +++ b/test_heisen.cc @@ -0,0 +1,22 @@ +#include "heisen.h" + +int main(int argc, char** argv) { +// // for (int i = 0; i <= 4; i++) { +// { int i = -1; +// printf("%d\n", i); +// HeisenCalculation calc(12,i); +// vector mat = calc.make_matrix(); +// int n = calc.get_nstates(); +// vector e(n,0.0); +// vector ev(n*n, 0.0); +// printf("\n"); +// print_matrix(mat,n,n); +// diag_matrix(mat,n,e,ev); +// print_vector(e); +// // printf("Lanczos: lowest eigenvalue is %15.8f\n\n", e[0]); +// } + HeisenCalculation calc(12); + vector e = calc.eigenvalues(); + print_vector(e); + return 0; +}