1 //////////////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (c) 2014 The Regents of the University of California 4 // 5 // This file is part of Qbox 6 // 7 // Qbox is distributed under the terms of the GNU General Public License 8 // as published by the Free Software Foundation, either version 2 of 9 // the License, or (at your option) any later version. 10 // See the file COPYING in the root directory of this distribution 11 // or <http://www.gnu.org/licenses/>. 12 // 13 //////////////////////////////////////////////////////////////////////////////// 14 // 15 // D3tensor.h 16 // 17 // double 3x3 tensor 18 // 19 //////////////////////////////////////////////////////////////////////////////// 20 21 #ifndef D3TENSOR_H 22 #define D3TENSOR_H 23 #include <iostream> 24 #include <iomanip> 25 #include <cmath> 26 #include <cassert> 27 #include "blas.h" 28 #include "D3vector.h" 29 30 using namespace std; 31 32 class D3tensor 33 { 34 public: 35 36 double a_[9]; 37 a(void)38 double* a(void) { return &a_[0]; } a(void)39 const double* a(void) const { return &a_[0]; } 40 D3tensor(void)41 explicit D3tensor(void) { clear(); } 42 D3tensor(double xx,double yy,double zz)43 explicit D3tensor(double xx, double yy, double zz) 44 { a_[0]=xx; a_[4]=yy; a_[8]=zz; } 45 D3tensor(double xx,double yy,double zz,double xy,double yz,double xz,char & uplo)46 explicit D3tensor(double xx, double yy, double zz, 47 double xy, double yz, double xz, char& uplo) 48 { 49 a_[0] = xx; 50 a_[4] = yy; 51 a_[8] = zz; 52 53 if ( uplo == 'l' ) 54 { 55 a_[1] = xy; 56 a_[2] = xz; 57 a_[5] = yz; 58 } 59 else if ( uplo == 'u' ) 60 { 61 a_[3] = xy; 62 a_[6] = xz; 63 a_[7] = yz; 64 } 65 else if ( uplo == 's' ) 66 { 67 a_[1] = xy; 68 a_[2] = xz; 69 a_[5] = yz; 70 a_[3] = xy; 71 a_[6] = xz; 72 a_[7] = yz; 73 } 74 else 75 assert(false); 76 } 77 D3tensor(const D3vector & diag,const D3vector & offdiag)78 explicit D3tensor(const D3vector& diag, const D3vector& offdiag) 79 { 80 a_[0] = diag[0]; 81 a_[4] = diag[1]; 82 a_[8] = diag[2]; 83 a_[1] = offdiag[0]; 84 a_[5] = offdiag[1]; 85 a_[2] = offdiag[2]; 86 a_[3] = offdiag[0]; 87 a_[7] = offdiag[1]; 88 a_[6] = offdiag[2]; 89 } 90 D3tensor(const double * a)91 explicit D3tensor(const double* a) 92 { 93 for ( int i = 0; i < 9; i++ ) a_[i] = a[i]; 94 } 95 96 double& operator[](int i) 97 { 98 assert(i>=0 && i < 9); 99 return a_[i]; 100 } 101 102 double operator[](int i) const 103 { 104 assert(i>=0 && i < 9); 105 return a_[i]; 106 } 107 setdiag(int i,double b)108 void setdiag(int i, double b) 109 { 110 assert(i>=0 && i<3); 111 a_[i*4] = b; 112 } 113 setdiag(const D3vector & b)114 void setdiag(const D3vector& b) 115 { 116 for ( int i = 0; i < 3; i++ ) 117 a_[i*4] = b[i]; 118 } 119 setoffdiag(int i,double b)120 void setoffdiag(int i, double b) 121 { 122 assert(i>=0 && i<3); 123 if ( i == 0 ) 124 { 125 a_[1] = b; 126 a_[3] = b; 127 } 128 else if ( i == 2 ) 129 { 130 a_[2] = b; 131 a_[6] = b; 132 } 133 else 134 { 135 a_[5] = b; 136 a_[7] = b; 137 } 138 } 139 setoffdiag(const D3vector & b)140 void setoffdiag(const D3vector& b) 141 { 142 a_[1] = b[0]; 143 a_[3] = b[0]; 144 a_[5] = b[1]; 145 a_[7] = b[1]; 146 a_[2] = b[2]; 147 a_[6] = b[2]; 148 } 149 150 bool operator==(const D3tensor &rhs) const 151 { 152 bool eq = true; 153 for ( int i = 0; i < 9; i++ ) 154 { 155 if ( rhs[i] != a_[i] ) 156 { 157 eq &= false; 158 } 159 } 160 return eq; 161 } 162 163 bool operator!=(const D3tensor &rhs) const 164 { 165 bool neq = false; 166 for ( int i = 0; i < 9; i++ ) 167 { 168 if ( rhs[i] != a_[i] ) 169 { 170 neq |= true; 171 } 172 } 173 return neq; 174 } 175 176 D3tensor& operator+=(const D3tensor& rhs) 177 { 178 for ( int i = 0; i < 9; i++ ) 179 a_[i] += rhs[i]; 180 return *this; 181 } 182 183 D3tensor& operator-=(const D3tensor& rhs) 184 { 185 for ( int i = 0; i < 9; i++ ) 186 a_[i] -= rhs[i]; 187 return *this; 188 } 189 190 D3tensor& operator*=(const double& rhs) 191 { 192 for ( int i = 0; i < 9; i++ ) 193 a_[i] *= rhs; 194 return *this; 195 } 196 197 D3tensor& operator/=(const double& rhs) 198 { 199 for ( int i = 0; i < 9; i++ ) 200 a_[i] /= rhs; 201 return *this; 202 } 203 204 friend const D3tensor operator+(const D3tensor& lhs, const D3tensor& rhs) 205 { 206 return D3tensor(lhs) += rhs; 207 } 208 209 friend const D3tensor operator-(const D3tensor& a, const D3tensor& b) 210 { 211 return D3tensor(a) -= b; 212 } 213 214 friend D3tensor operator*(double a, const D3tensor& b) 215 { 216 return D3tensor(b) *= a; 217 } 218 219 friend D3tensor operator*(const D3tensor& a, double b) 220 { 221 return D3tensor(a) *= b; 222 } 223 224 friend D3tensor operator/(const D3tensor& a, double b) 225 { 226 return D3tensor(a) /= b; 227 } 228 229 friend D3tensor operator*(D3tensor& a, D3tensor& b) 230 { 231 D3tensor c; 232 int ithree = 3; 233 double one = 1.0, zero = 0.0; 234 char t = 'n'; 235 dgemm ( &t, &t, &ithree, &ithree, &ithree, &one, &a[0], &ithree, 236 &b[0], &ithree, &zero, &c[0], &ithree ); 237 return c; 238 } 239 240 friend D3tensor operator-(const D3tensor& a) // unary minus 241 { 242 return D3tensor()-a; 243 } 244 norm2(const D3tensor & a)245 double norm2(const D3tensor& a) const 246 { 247 return a[0]*a[0] + a[1]*a[1] + a[2]*a[2] + 248 a[3]*a[3] + a[4]*a[4] + a[5]*a[5] + 249 a[6]*a[6] + a[7]*a[7] + a[8]*a[8]; 250 } 251 norm(const D3tensor & a)252 double norm(const D3tensor& a) const 253 { 254 return sqrt(norm2(a)); 255 } 256 trace(void)257 double trace(void) const 258 { 259 return a_[0]+a_[4]+a_[8]; 260 } 261 traceless(void)262 void traceless(void) 263 { 264 double b = trace() / 3; 265 a_[0] -= b; 266 a_[4] -= b; 267 a_[8] -= b; 268 } 269 clear(void)270 void clear(void) 271 { 272 for ( int i = 0; i < 9; i++ ) 273 a_[i] = 0.0; 274 } 275 identity(void)276 void identity(void) 277 { 278 clear(); 279 a_[0] = 1.0; 280 a_[4] = 1.0; 281 a_[8] = 1.0; 282 } 283 284 friend std::ostream& operator<<(std::ostream& os, const D3tensor& rhs) 285 { 286 const double * const v = rhs.a(); 287 os.setf(ios::fixed,ios::floatfield); 288 os.setf(ios::right,ios::adjustfield); 289 os.precision(8); 290 os << setw(14) << v[0] << " " << setw(14) << v[3] << " " << setw(14) << v[6] 291 << "\n" 292 << setw(14) << v[1] << " " << setw(14) << v[4] << " " << setw(14) << v[7] 293 << "\n" 294 << setw(14) << v[2] << " " << setw(14) << v[5] << " " << setw(14) << v[8] 295 << "\n"; 296 return os; 297 } 298 syev(char uplo,D3vector & eigval,D3tensor & eigvec)299 void syev(char uplo, D3vector& eigval, D3tensor& eigvec) 300 { 301 double w[3]; 302 303 int info; 304 char jobz = 'V'; 305 int lwork=-1; 306 double tmplwork; 307 int n = 3; 308 309 dsyev(&jobz, &uplo, &n, eigvec.a(), &n, &w[0], &tmplwork, &lwork, &info); 310 311 lwork = (int) tmplwork + 1; 312 double* work = new double[lwork]; 313 314 eigvec = *this; 315 dsyev(&jobz, &uplo, &n, eigvec.a(), &n, &w[0], work, &lwork, &info); 316 delete[] work; 317 318 eigval = D3vector(&w[0]); 319 } 320 }; 321 #endif 322