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