1 /*
2
3 Copyright (C) 2007 Burkhard Militzer
4 with some modifications by Lucas K. Wagner
5 and extensions to Pfaffian algebra by Michal Bajdich
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License along
18 with this program; if not, write to the Free Software Foundation, Inc.,
19 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
20
21 */
22 #ifndef MATRIXALGEBRA_H_INCLUDED
23 #define MATRIXALGEBRA_H_INCLUDED
24
25 #include "Array.h"
26 #include "Qmc_std.h"
27
28
29 /*!
30 This is basically a replacement for doubles, but always with a logarithmic representation.
31 I've tried to make it so it automatically does the right thing with standard
32 multiplication. Automatic downconversion to double/complex is possible, but I've excluded it to
33 prevent accidental loss of precision. Use the explicit val() function instead.
34 Note also that there may be performance problems with using these automatic conversions, so if performance is critical, one may need to write some specialized code.
35 */
36
37 template <class T> struct log_value {
38 T logval;
39 int sign;
vallog_value40 T val() { return doublevar(sign)*exp(logval); } //doublevar() a bit ugly, but it works ok
log_valuelog_value41 log_value() { logval=0; sign=1; }
log_valuelog_value42 log_value(T t) { }
43 log_value<T> & operator *=(const log_value<T> & right) {
44 this->logval+=right.logval;
45 this->sign*=right.sign;
46 return *this;
47 }
48 };
49
log_value(doublevar t)50 template<> inline log_value<doublevar>::log_value(doublevar t) {
51 doublevar ft=fabs(t);
52 if(ft > 0) logval=log(ft);
53 else logval=-1e201;
54 sign=t<0?-1:1;
55 }
56
log_value(dcomplex t)57 template<> inline log_value<dcomplex>::log_value(dcomplex t) {
58 sign=1;
59 doublevar ft=abs(t);
60 if(ft > 0) {
61 logval=dcomplex(log(abs(t)),arg(t));
62 }
63 else logval=dcomplex(-1e201,0.0);
64 }
65 typedef log_value<doublevar> log_real_value;
66 typedef log_value<dcomplex> log_complex_value;
67
68 //--------Finished with the class definitions
69 //Here are a few helper functions.
70 template <class T> inline log_value<T> operator*(T t, log_value<T> u) {
71 log_value<T> v(t);
72 v*=u;
73 return v;
74 }
75
76 template <class T> inline log_value<T> operator*(log_value<T> u,double t) { return t*u; }
77 template <class T> inline log_value<T> operator*(log_value<T> t, log_value<T> u) {
78 log_value<T> v;
79 v.logval=t.logval+u.logval;
80 v.sign=t.sign*u.sign;
81 return v;
82 }
83
84 inline log_value<dcomplex> operator*(log_value<doublevar> & t, log_value<dcomplex> & u) {
85 log_value<dcomplex> v;
86 v.logval=t.logval+u.logval;
87 v.sign=t.sign*u.sign;
88 return v;
89 }
90
91 //Try to safely sum a series of log_values
sum(const Array1<log_value<T>> & vec)92 template <class T> inline log_value<T> sum(const Array1 <log_value<T> > & vec) {
93 T s=0;
94 int n=vec.GetDim(0);
95 int piv=0;
96 for(int i=0; i< n; i++)
97 s+=T(vec(i).sign*vec(piv).sign)*exp(vec(i).logval-vec(piv).logval);
98 return s*vec(piv);
99 }
100
101 //--------------------------------------------------------------
102
103 doublevar dot(const Array1 <doublevar> & a, const Array1 <doublevar> & b);
104 dcomplex dot(const Array1 <doublevar> & a, const Array1 <dcomplex> & b);
105 dcomplex dot(const Array1 <dcomplex> & a, const Array1 <doublevar> & b);
106 dcomplex dot(const Array1 <dcomplex> & a, const Array1 <dcomplex> & b);
107
108
109 doublevar Determinant(const Array2 <doublevar> & a, const int n);
110 dcomplex Determinant(const Array2 <dcomplex> & a, const int n);
111
112 void MultiplyMatrices(const Array2 <doublevar> & a, const Array2 <doublevar> & b,
113 Array2 <doublevar> & c, int n);
114
115 void InvertMatrix(const Array2 <doublevar> & a, Array2 <doublevar> & a1, const int n);
116 doublevar InverseUpdateRow(Array2 <doublevar> & a1, const Array2 <doublevar> & a,
117 const int lRow, const int n);
118 doublevar InverseUpdateRow(Array2 <doublevar> & a1, const Array1 <doublevar> & newRow,
119 const int lRow, const int n);
120 //doublevar InverseGetNewRatio(const Array2 <doublevar> & a1, const Array1 <doublevar> & newCol,
121 // const int lCol, const int n);
122 //Get the new ratio without updating the inverse..
123
InverseGetNewRatio(const Array2<T> & a1,const Array1<T> & newCol,const int lCol,const int n)124 template <class T> T InverseGetNewRatio(const Array2 <T> & a1,
125 const Array1 <T> & newCol,
126 const int lCol, const int n) {
127 T f=T(0.0);
128 for(int i=0;i<n;++i) {
129 f += a1(lCol,i)*newCol[i];
130 }
131 f =1.0/f;
132
133 return f;
134
135 }
136
137
138 doublevar InverseGetNewRatioRow(const Array2 <doublevar> & a1, const Array1 <doublevar> & newRow,
139 const int lRow, const int n);
140 doublevar InverseUpdateColumn(Array2 <doublevar> & a1, const Array2 <doublevar> & a,
141 const int lCol, const int n);
142 doublevar InverseUpdateColumn(Array2 <doublevar> & a1, const Array1 <doublevar> & newCol,
143 const int lCol, const int n);
144
145 void TransposeMatrix(Array2 <doublevar> & a, const int n);
146 log_real_value TransposeInverseMatrix(const Array2 <doublevar> & a, Array2 <doublevar> & a1, const int n);
147 doublevar TransposeInverseUpdateRow(Array2 <doublevar> & a1, const Array2 <doublevar> & a,
148 const int lCol, const int n);
149 doublevar TransposeInverseUpdateColumn(Array2 <doublevar> & a1, const Array2 <doublevar> & a,
150 const int lRow, const int n);
151 doublevar Pfaffian_nopivot(const Array2 <doublevar> & a);
152
153 doublevar Pfaffian_partialpivot(const Array2 <doublevar> & a);
154
155 doublevar UpdateInversePfaffianMatrix(Array2 <doublevar> & a, Array1 <doublevar> & row,
156 Array1 <doublevar> & column, int n);
157 doublevar GetUpdatedPfaffianValue(Array2 <doublevar> & in,
158 Array1 <doublevar> & row,
159 int e);
160 doublevar PfaffianInverseMatrix(const Array2 <doublevar> & a, Array2 <doublevar> & a1);
161
162
cabs(dcomplex & a)163 inline doublevar cabs(dcomplex & a) {
164 return sqrt(a.real()*a.real()+a.imag()*a.imag());
165 }
166
167 int ludcmp(Array2 <doublevar> & a, const int n, Array1 <int> & indx,
168 doublevar & d);
169 void lubksb(Array2 <doublevar> & a, int n, Array1 <int> & indx,
170 Array1 <doublevar> & b);
171
172 void EigenSystemSolverRealSymmetricMatrix(const Array2 < doublevar > & Ain, Array1 < doublevar> & evals, Array2 < doublevar> & evecs);
173
174 void DGGEV(Array2 <doublevar> & A, Array2 <doublevar> & B,
175 Array1 <dcomplex> & alpha,
176 Array2 <doublevar> & VL, Array2 <doublevar> & VR);
177
178 void GeneralizedEigenSystemSolverRealSymmetricMatrices(const Array2 < doublevar > & Ain, const Array2 < doublevar> & Bin, Array1 < doublevar> & evals, Array2 < doublevar> & evecs);
179 void GeneralizedEigenSystemSolverComplexGeneralMatrices(Array2 < dcomplex > & Ain,
180 Array1 <dcomplex> & W, Array2 <dcomplex> & VL, Array2 <dcomplex> & VR);
181 void GeneralizedEigenSystemSolverRealGeneralMatrices(Array2 < doublevar > & Ain,
182 Array1 <dcomplex> & W, Array2 <doublevar> & VL, Array2 <doublevar> & VR);
183
184
185 void sort_abs_values_descending(const Array1 <doublevar> & vals, Array1 <doublevar> & newvals, Array1 <int> & list);
186 void sort_according_list(const Array1 <doublevar> & vals, Array1 <doublevar> & newvals, const Array1 <int> & list);
187 //------------------------------------------------------------------------
188 //IF Lapack was used
189 #ifdef USE_LAPACK
190 extern "C" {
191 void dsyevr_(char *JOBZp, char *RANGEp, char *UPLOp, int *Np,
192 double *A, int *LDAp, double *VLp, double *VUp,
193 int *ILp, int *IUp, double *ABSTOLp, int *Mp,
194 double *W, double *Z, int *LDZp, int *ISUPPZ,
195 double *WORK, int *LWORKp, int *IWORK, int *LIWORKp,
196 int *INFOp);
197 double dlamch_(char *CMACHp);
198 void dsygv_(int *itype, char *jobz, char *uplo, int *n, double *a, int *lda,
199 double *b, int *ldb, double *w, double *WORK, int *IWORK, int *info);
200 void dgetrf_(int * m, int * n, double * A, int * lda, int * ipiv, int *info);
201 void dgetrs_(char * trans, int *n, int * nrhs, double * A, int * lda,
202 int *ipiv, double * B, int * ldb, int * info);
203
204 void dgeev_(char * JOBVL, char * JOBVR,int * N,double * A,int * lda,
205 double * WR, double * WI,double *VL, int *LDVL,double * VR, int *LDVR,
206 double * WORK, int * LWORK, int * INFO);
207
208 void zgeev_(char * JOBVL, char * JOBVR,int * N,void * A,int * lda,
209 void * W, void *VL, int *LDVL,void * VR, int *LDVR,
210 void * WORK, int * LWORK,void * RWORK, int * INFO);
211
212 void dggev_(char * JOBVL, char * JOBVR, int * N, void * A, int * lda,
213 void * B, int * ldb, void * alphar, void * alphai, void * beta,
214 void* vl, int * ldvl, void * vr, int * ldvr, void * work, int * lwork, int * info);
215 };
216 #endif
217
218 //----------------------------------------------------------------------
219 //complex versions of determinant stuff
220
221 int ludcmp(Array2 <dcomplex > & a, const int n,
222 Array1 <int> & indx, doublevar & d);
223 void lubksb(Array2 < dcomplex > & a, int n,
224 Array1 <int> & indx, Array1 <dcomplex > & b);
225
226 log_value <dcomplex>
227 TransposeInverseMatrix(const Array2 <dcomplex > & a,
228 Array2 <dcomplex> & a1,
229 const int n);
230 dcomplex
231 InverseUpdateColumn(Array2 <dcomplex > & a1,
232 const Array1 <dcomplex > & newCol,
233 const int lCol, const int n);
234 dcomplex InverseUpdateRow(Array2 <dcomplex> & a1, const Array1 <dcomplex> & newRow,
235 const int lRow, const int n);
236 void Jacobi(const Array2 < dcomplex > & Ain, Array1 <doublevar> & evals, Array2 < dcomplex > & evecs);
237
238 #endif // MATRIXALGEBRA_H_INCLUDED
239 //--------------------------------------------------------------------------
240