From b23e336fc8e611e978be9ee947d9835e172b24b1 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Sun, 27 Sep 2015 15:01:53 -0400 Subject: [PATCH] commit with existing files --- heisen.h | 376 +++++++++++++++++ heisen_chain.f90 | 177 ++++++++ hubbard.h | 1013 ++++++++++++++++++++++++++++++++++++++++++++++ hubbard.in | 7 + lanczos.h | 290 +++++++++++++ lanczos.m | 56 +++ main.cc | 169 ++++++++ matrix.h | 368 +++++++++++++++++ read.in | 1 + test_heisen.cc | 22 + 10 files changed, 2479 insertions(+) create mode 100755 heisen.h create mode 100644 heisen_chain.f90 create mode 100755 hubbard.h create mode 100755 hubbard.in create mode 100755 lanczos.h create mode 100755 lanczos.m create mode 100755 main.cc create mode 100755 matrix.h create mode 100644 read.in create mode 100644 test_heisen.cc 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; +}