Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMatrix4x4.h
Go to the documentation of this file.
1 #include <align_16.h>
2 
3 #ifndef USE_SSE2
4 #include "DMatrixDSym.h"
5 // Matrix class without SIMD instructions
6 
7 class DMatrix4x4{
8  public:
10  for (unsigned int i=0;i<4;i++){
11  for (unsigned int j=0;j<4;j++){
12  mA[i][j]=0.;
13  }
14  }
15  }
16  DMatrix4x4(double c11, double c12, double c13, double c14,
17  double c21, double c22, double c23, double c24,
18  double c31, double c32, double c33, double c34,
19  double c41, double c42, double c43, double c44
20  ){
21  mA[0][0]=c11;
22  mA[0][1]=c12;
23  mA[0][2]=c13;
24  mA[0][3]=c14;
25  mA[1][0]=c21;
26  mA[1][1]=c22;
27  mA[1][2]=c23;
28  mA[1][3]=c24;
29  mA[2][0]=c31;
30  mA[2][1]=c32;
31  mA[2][2]=c33;
32  mA[2][3]=c34;
33  mA[3][0]=c41;
34  mA[3][1]=c42;
35  mA[3][2]=c43;
36  mA[3][3]=c44;
37 
38  }
39  DMatrix4x4(const DMatrix2x2 &m1,const DMatrix2x2 &m2,
40  const DMatrix2x2 &m3,const DMatrix2x2 &m4){
41  mA[0][0]=m1(0,0);
42  mA[0][1]=m1(0,1);
43  mA[1][0]=m1(1,0);
44  mA[1][1]=m1(1,1);
45  mA[0][2]=m2(0,0);
46  mA[0][3]=m2(0,1);
47  mA[1][2]=m2(1,0);
48  mA[1][3]=m2(1,1);
49  mA[2][0]=m3(0,0);
50  mA[2][1]=m3(0,1);
51  mA[3][0]=m3(1,0);
52  mA[3][1]=m3(1,1);
53  mA[2][2]=m4(0,0);
54  mA[2][3]=m4(0,1);
55  mA[3][2]=m4(1,0);
56  mA[3][3]=m4(1,1);
57  }
58 
60 
61  double &operator() (int row, int col){
62  return mA[row][col];
63  }
64  double operator() (int row, int col) const{
65  return mA[row][col];
66  }
67  // Assignment
69  for (unsigned int i=0;i<4;i++){
70  for (unsigned int j=0;j<4;j++){
71  mA[i][j]=m1(i,j);
72  }
73  }
74  return *this;
75  }
76 
77 
78 
79 
80  // unary minus
82  return DMatrix4x4(-mA[0][0],-mA[0][1],-mA[0][2],-mA[0][3],
83  -mA[1][0],-mA[1][1],-mA[1][2],-mA[1][3],
84  -mA[2][0],-mA[2][1],-mA[2][2],-mA[2][3],
85  -mA[3][0],-mA[3][1],-mA[3][2],-mA[3][3]
86  );
87  }
88  // Matrix Addition
90  return DMatrix4x4(mA[0][0]+m2(0,0),mA[0][1]+m2(0,1),mA[0][2]+m2(0,2),mA[0][3]+m2(0,3),
91  mA[1][0]+m2(1,0),mA[1][1]+m2(1,1),mA[1][2]+m2(1,2),mA[1][3]+m2(1,3),
92  mA[2][0]+m2(2,0),mA[2][1]+m2(2,1),mA[2][2]+m2(2,2),mA[2][3]+m2(2,3),
93  mA[3][0]+m2(3,0),mA[3][1]+m2(3,1),mA[3][2]+m2(3,2),mA[3][3]+m2(3,3)
94  );
95  }
96 
97 
98  // Matrix subtraction
100  return DMatrix4x4(mA[0][0]-m2(0,0),mA[0][1]-m2(0,1),mA[0][2]-m2(0,2),mA[0][3]-m2(0,3),
101  mA[1][0]-m2(1,0),mA[1][1]-m2(1,1),mA[1][2]-m2(1,2),mA[1][3]-m2(1,3),
102  mA[2][0]-m2(2,0),mA[2][1]-m2(2,1),mA[2][2]-m2(2,2),mA[2][3]-m2(2,3),
103  mA[3][0]-m2(3,0),mA[3][1]-m2(3,1),mA[3][2]-m2(3,2),mA[3][3]-m2(3,3)
104  );
105  }
106 
107  // Matrix multiplication: (4x4) x (4x4)
109  return DMatrix4x4(mA[0][0]*m2(0,0)+mA[0][1]*m2(1,0)+mA[0][2]*m2(2,0)+mA[0][3]*m2(3,0),
110  mA[0][0]*m2(0,1)+mA[0][1]*m2(1,1)+mA[0][2]*m2(2,1)+mA[0][3]*m2(3,1),
111  mA[0][0]*m2(0,2)+mA[0][1]*m2(1,2)+mA[0][2]*m2(2,2)+mA[0][3]*m2(3,2),
112  mA[0][0]*m2(0,3)+mA[0][1]*m2(1,3)+mA[0][2]*m2(2,3)+mA[0][3]*m2(3,3),
113 
114  mA[1][0]*m2(0,0)+mA[1][1]*m2(1,0)+mA[1][2]*m2(2,0)+mA[1][3]*m2(3,0),
115  mA[1][0]*m2(0,1)+mA[1][1]*m2(1,1)+mA[1][2]*m2(2,1)+mA[1][3]*m2(3,1),
116  mA[1][0]*m2(0,2)+mA[1][1]*m2(1,2)+mA[1][2]*m2(2,2)+mA[1][3]*m2(3,2),
117  mA[1][0]*m2(0,3)+mA[1][1]*m2(1,3)+mA[1][2]*m2(2,3)+mA[1][3]*m2(3,3),
118 
119  mA[2][0]*m2(0,0)+mA[2][1]*m2(1,0)+mA[2][2]*m2(2,0)+mA[2][3]*m2(3,0),
120  mA[2][0]*m2(0,1)+mA[2][1]*m2(1,1)+mA[2][2]*m2(2,1)+mA[2][3]*m2(3,1),
121  mA[2][0]*m2(0,2)+mA[2][1]*m2(1,2)+mA[2][2]*m2(2,2)+mA[2][3]*m2(3,2),
122  mA[2][0]*m2(0,3)+mA[2][1]*m2(1,3)+mA[2][2]*m2(2,3)+mA[2][3]*m2(3,3),
123 
124  mA[3][0]*m2(0,0)+mA[3][1]*m2(1,0)+mA[3][2]*m2(2,0)+mA[3][3]*m2(3,0),
125  mA[3][0]*m2(0,1)+mA[3][1]*m2(1,1)+mA[3][2]*m2(2,1)+mA[3][3]*m2(3,1),
126  mA[3][0]*m2(0,2)+mA[3][1]*m2(1,2)+mA[3][2]*m2(2,2)+mA[3][3]*m2(3,2),
127  mA[3][0]*m2(0,3)+mA[3][1]*m2(1,3)+mA[3][2]*m2(2,3)+mA[3][3]*m2(3,3)
128  );
129  }
130 
131  // Matrix multiplication: (4x4) x (4x2)
133  return DMatrix4x2(mA[0][0]*m2(0,0)+mA[0][1]*m2(1,0)+mA[0][2]*m2(2,0)+mA[0][3]*m2(3,0),
134  mA[0][0]*m2(0,1)+mA[0][1]*m2(1,1)+mA[0][2]*m2(2,1)+mA[0][3]*m2(3,1),
135 
136  mA[1][0]*m2(0,0)+mA[1][1]*m2(1,0)+mA[1][2]*m2(2,0)+mA[1][3]*m2(3,0),
137  mA[1][0]*m2(0,1)+mA[1][1]*m2(1,1)+mA[1][2]*m2(2,1)+mA[1][3]*m2(3,1),
138 
139  mA[2][0]*m2(0,0)+mA[2][1]*m2(1,0)+mA[2][2]*m2(2,0)+mA[2][3]*m2(3,0),
140  mA[2][0]*m2(0,1)+mA[2][1]*m2(1,1)+mA[2][2]*m2(2,1)+mA[2][3]*m2(3,1),
141 
142  mA[3][0]*m2(0,0)+mA[3][1]*m2(1,0)+mA[3][2]*m2(2,0)+mA[3][3]*m2(3,0),
143  mA[3][0]*m2(0,1)+mA[3][1]*m2(1,1)+mA[3][2]*m2(2,1)+mA[3][3]*m2(3,1)
144 
145  );
146  }
147 
148 
149  // Matrix multiplication: (4x4) x (4x1)
151  return DMatrix4x1(mA[0][0]*m2(0)+mA[0][1]*m2(1)+mA[0][2]*m2(2)+mA[0][3]*m2(3),
152  mA[1][0]*m2(0)+mA[1][1]*m2(1)+mA[1][2]*m2(2)+mA[1][3]*m2(3),
153  mA[2][0]*m2(0)+mA[2][1]*m2(1)+mA[2][2]*m2(2)+mA[2][3]*m2(3),
154  mA[3][0]*m2(0)+mA[3][1]*m2(1)+mA[3][2]*m2(2)+mA[3][3]*m2(3)
155  );
156  }
157 
158 
159 
160 
162  DMatrix2x2 F(mA[0][0],mA[0][1],mA[1][0],mA[1][1]);
163  DMatrix2x2 Finv=F.Invert();
164  DMatrix2x2 G(mA[0][2],mA[0][3],mA[1][2],mA[1][3]);
165  DMatrix2x2 H(mA[2][0],mA[2][1],mA[3][0],mA[3][1]);
166  DMatrix2x2 J(mA[2][2],mA[2][3],mA[3][2],mA[3][3]);
167  DMatrix2x2 Jinv=J.Invert();
168  DMatrix2x2 FF=(F-G*Jinv*H).Invert();
169  DMatrix2x2 JJ=(J-H*Finv*G).Invert();
170  return DMatrix4x4(FF,-FF*G*Jinv,-JJ*H*Finv,JJ);
171  }
172 
173  // Find the transpose of this matrix
176  for (unsigned int i=0;i<4;i++){
177  for (unsigned int j=0;j<4;j++){
178  temp(i,j)=mA[j][i];
179  }
180  }
181  return temp;
182  }
183 
184  DMatrixDSym GetSub(unsigned int lowerBound, unsigned int upperBound){
185  if (upperBound >= lowerBound) return DMatrixDSym();
186  DMatrixDSym subMatrix(upperBound - lowerBound);
187  for (unsigned int i=lowerBound; i <= upperBound; i++){
188  for (unsigned int j=lowerBound; j <= upperBound; j++){
189  subMatrix(i,j) = mA[i][j];
190  }
191  }
192  return subMatrix;
193  }
194 
195  bool IsPosDef(){
196  if(mA[0][0] > 0. &&
197  GetSub(0,1).Determinant() > 0. && GetSub(0,2).Determinant() > 0. &&
198  GetSub(0,3).Determinant() > 0.)
199  return true;
200  else return false;
201  }
202 
203  void Print(){
204  cout << "DMatrix4x4:" <<endl;
205  cout << " | 0 | 1 | 2 | 3 |" <<endl;
206  cout << "------------------------------------------------------" <<endl;
207 
208  for (unsigned int i=0;i<4;i++){
209  cout << " " << i <<" |";
210  for (unsigned int j=0;j<4;j++){
211  cout <<setw(11)<<setprecision(4)<< mA[i][j] <<" ";
212  }
213  cout << endl;
214  }
215  }
216 
217 
218  private:
219  double mA[4][4];
220 
221 };
222 
223 #else
224 
225 class DMatrix4x4{
226  public:
227  DMatrix4x4()
228  : mA( ALIGNED_16_BLOCK_PTR(union dvec, 4, mA) )
229  {
230  for (unsigned int i=0;i<4;i++){
231  mA[i].v[0]=_mm_setzero_pd();
232  mA[i].v[1]=_mm_setzero_pd();
233  }
234  }
235  DMatrix4x4(double c11, double c12, double c13, double c14,
236  double c21, double c22, double c23, double c24,
237  double c31, double c32, double c33, double c34,
238  double c41, double c42, double c43, double c44)
239  : mA( ALIGNED_16_BLOCK_PTR(union dvec, 4, mA) )
240  {
241  mA[0].v[0]=_mm_setr_pd(c11,c21);
242  mA[0].v[1]=_mm_setr_pd(c31,c41);
243  mA[1].v[0]=_mm_setr_pd(c12,c22);
244  mA[1].v[1]=_mm_setr_pd(c32,c42);
245  mA[2].v[0]=_mm_setr_pd(c13,c23);
246  mA[2].v[1]=_mm_setr_pd(c33,c43);
247  mA[3].v[0]=_mm_setr_pd(c14,c24);
248  mA[3].v[1]=_mm_setr_pd(c34,c44);
249  }
250 
251  DMatrix4x4(__m128d m11, __m128d m12, __m128d m13, __m128d m14,
252  __m128d m21, __m128d m22, __m128d m23, __m128d m24)
253  : mA( ALIGNED_16_BLOCK_PTR(union dvec, 4, mA) )
254  {
255  mA[0].v[0]=m11;
256  mA[1].v[0]=m12;
257  mA[2].v[0]=m13;
258  mA[3].v[0]=m14;
259  mA[0].v[1]=m21;
260  mA[1].v[1]=m22;
261  mA[2].v[1]=m23;
262  mA[3].v[1]=m24;
263  }
264  DMatrix4x4(const DMatrix2x2 &m1,const DMatrix2x2 &m2,
265  const DMatrix2x2 &m3,const DMatrix2x2 &m4)
266  : mA( ALIGNED_16_BLOCK_PTR(union dvec, 4, mA) )
267  {
268  mA[0].v[0]=m1.GetV(0);
269  mA[1].v[0]=m1.GetV(1);
270  mA[0].v[1]=m3.GetV(0);
271  mA[1].v[1]=m3.GetV(1);
272  mA[2].v[0]=m2.GetV(0);
273  mA[3].v[0]=m2.GetV(1);
274  mA[2].v[1]=m4.GetV(0);
275  mA[3].v[1]=m4.GetV(1);
276  }
277  DMatrix4x4(const DMatrix4x4& dm)
278  : mA( ALIGNED_16_BLOCK_PTR(union dvec, 4, mA) )
279  {
280  mA[0].v[0]=dm.mA[0].v[0];
281  mA[1].v[0]=dm.mA[1].v[0];
282  mA[2].v[0]=dm.mA[2].v[0];
283  mA[3].v[0]=dm.mA[3].v[0];
284  mA[0].v[1]=dm.mA[0].v[1];
285  mA[1].v[1]=dm.mA[1].v[1];
286  mA[2].v[1]=dm.mA[2].v[1];
287  mA[3].v[1]=dm.mA[3].v[1];
288  }
289  ~DMatrix4x4(){};
290 
291  __m128d GetV(int pair, int col) const{
292  return mA[col].v[pair];
293  }
294  void SetV(int pair, int col, __m128d v){
295  mA[col].v[pair]=v;
296  }
297 
298  double &operator() (int row, int col){
299  return mA[col].d[row];
300  }
301  double operator() (int row, int col) const{
302  return mA[col].d[row];
303  }
304 
306  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 1, p)
307  __m128d &zero=p[0];
308  zero=_mm_setzero_pd();
309  return DMatrix4x4(_mm_sub_pd(zero,GetV(0,0)),
310  _mm_sub_pd(zero,GetV(0,1)),
311  _mm_sub_pd(zero,GetV(0,2)),
312  _mm_sub_pd(zero,GetV(0,3)),
313  _mm_sub_pd(zero,GetV(1,0)),
314  _mm_sub_pd(zero,GetV(1,1)),
315  _mm_sub_pd(zero,GetV(1,2)),
316  _mm_sub_pd(zero,GetV(1,3)));
317  }
318 
319  DMatrix4x4 operator-(const DMatrix4x4 &m2){
320  return DMatrix4x4(_mm_sub_pd(GetV(0,0),m2.GetV(0,0)),
321  _mm_sub_pd(GetV(0,1),m2.GetV(0,1)),
322  _mm_sub_pd(GetV(0,2),m2.GetV(0,2)),
323  _mm_sub_pd(GetV(0,3),m2.GetV(0,3)),
324  _mm_sub_pd(GetV(1,0),m2.GetV(1,0)),
325  _mm_sub_pd(GetV(1,1),m2.GetV(1,1)),
326  _mm_sub_pd(GetV(1,2),m2.GetV(1,2)),
327  _mm_sub_pd(GetV(1,3),m2.GetV(1,3)));
328  }
329 
330  DMatrix4x4 operator+(const DMatrix4x4 &m2){
331  return DMatrix4x4(_mm_add_pd(GetV(0,0),m2.GetV(0,0)),
332  _mm_add_pd(GetV(0,1),m2.GetV(0,1)),
333  _mm_add_pd(GetV(0,2),m2.GetV(0,2)),
334  _mm_add_pd(GetV(0,3),m2.GetV(0,3)),
335  _mm_add_pd(GetV(1,0),m2.GetV(1,0)),
336  _mm_add_pd(GetV(1,1),m2.GetV(1,1)),
337  _mm_add_pd(GetV(1,2),m2.GetV(1,2)),
338  _mm_add_pd(GetV(1,3),m2.GetV(1,3)));
339  }
340 
341 
342  DMatrix4x2 operator*(const DMatrix4x2 &m2){
343  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 8, p);
344  __m128d &m11=p[0];
345  __m128d &m12=p[1];
346  __m128d &m21=p[2];
347  __m128d &m22=p[3];
348  __m128d &m31=p[4];
349  __m128d &m32=p[5];
350  __m128d &m41=p[6];
351  __m128d &m42=p[7];
352  m11=_mm_set1_pd(m2(0,0));
353  m12=_mm_set1_pd(m2(0,1));
354  m21=_mm_set1_pd(m2(1,0));
355  m22=_mm_set1_pd(m2(1,1));
356  m31=_mm_set1_pd(m2(2,0));
357  m32=_mm_set1_pd(m2(2,1));
358  m41=_mm_set1_pd(m2(3,0));
359  m42=_mm_set1_pd(m2(3,1));
360  return
361  DMatrix4x2(_mm_add_pd(_mm_mul_pd(GetV(0,0),m11),
362  _mm_add_pd(_mm_mul_pd(GetV(0,1),m21),
363  _mm_add_pd(_mm_mul_pd(GetV(0,2),m31),
364  _mm_mul_pd(GetV(0,3),m41)))),
365  _mm_add_pd(_mm_mul_pd(GetV(0,0),m12),
366  _mm_add_pd(_mm_mul_pd(GetV(0,1),m22),
367  _mm_add_pd(_mm_mul_pd(GetV(0,2),m32),
368  _mm_mul_pd(GetV(0,3),m42)))),
369  _mm_add_pd(_mm_mul_pd(GetV(1,0),m11),
370  _mm_add_pd(_mm_mul_pd(GetV(1,1),m21),
371  _mm_add_pd(_mm_mul_pd(GetV(1,2),m31),
372  _mm_mul_pd(GetV(1,3),m41)))),
373  _mm_add_pd(_mm_mul_pd(GetV(1,0),m12),
374  _mm_add_pd(_mm_mul_pd(GetV(1,1),m22),
375  _mm_add_pd(_mm_mul_pd(GetV(1,2),m32),
376  _mm_mul_pd(GetV(1,3),m42)))));
377 
378  }
379 
380 
381  DMatrix4x4 &operator=(const DMatrix4x4 &m1){
382  for (unsigned int i=0;i<4;i++){
383  mA[i].v[0]=m1.GetV(0,i);
384  mA[i].v[1]=m1.GetV(1,i);
385  }
386  return *this;
387  }
388 
389  DMatrix4x4 Invert(){
390  DMatrix2x2 F(GetV(0,0),GetV(0,1));
391  DMatrix2x2 Finv=F.Invert();
392  DMatrix2x2 G(GetV(0,2),GetV(0,3));
393  DMatrix2x2 H(GetV(1,0),GetV(1,1));
394  DMatrix2x2 J(GetV(1,2),GetV(1,3));
395  DMatrix2x2 Jinv=J.Invert();
396  DMatrix2x2 FF=(F-G*Jinv*H).Invert();
397  DMatrix2x2 JJ=(J-H*Finv*G).Invert();
398  return DMatrix4x4(FF,-FF*G*Jinv,-JJ*H*Finv,JJ);
399  }
400 
401  // Find the transpose of this matrix
403 #define SWAP(i,j,k,m) _mm_setr_pd(mA[(j)].d[(i)],mA[(m)].d[(k)])
404 
405  return DMatrix4x4(SWAP(0,0,0,1),SWAP(1,0,1,1),SWAP(2,0,2,1),SWAP(3,0,3,1),
406  SWAP(0,2,0,3),SWAP(1,2,1,3),SWAP(2,2,2,3),SWAP(3,2,3,3));
407 
408  }
409 
410  // Matrix multiplication: (4x4) x (4x1)
411  DMatrix4x1 operator*(const DMatrix4x1 &m2){
412  ALIGNED_16_BLOCK_WITH_PTR(__m128d, 4, p)
413  __m128d &a1=p[0];
414  __m128d &a2=p[1];
415  __m128d &a3=p[2];
416  __m128d &a4=p[3];
417  a1=_mm_set1_pd(m2(0));
418  a2=_mm_set1_pd(m2(1));
419  a3=_mm_set1_pd(m2(2));
420  a4=_mm_set1_pd(m2(3));
421  return
422  DMatrix4x1(_mm_add_pd(_mm_mul_pd(GetV(0,0),a1),
423  _mm_add_pd(_mm_mul_pd(GetV(0,1),a2),
424  _mm_add_pd(_mm_mul_pd(GetV(0,2),a3),
425  _mm_mul_pd(GetV(0,3),a4)))),
426  _mm_add_pd(_mm_mul_pd(GetV(1,0),a1),
427  _mm_add_pd(_mm_mul_pd(GetV(1,1),a2),
428  _mm_add_pd(_mm_mul_pd(GetV(1,2),a3),
429  _mm_mul_pd(GetV(1,3),a4)))));
430  }
431 
432  // Matrix multiplication: (4x4) x (4x4)
433  DMatrix4x4 operator*(const DMatrix4x4 &m2){
434  struct dvec1{
435  __m128d v[4][4];
436  };
437  ALIGNED_16_BLOCK_WITH_PTR(struct dvec1, 1, p)
438  struct dvec1 &temp=p[0];
439  for (unsigned int i=0;i<4;i++){
440  for (unsigned int j=0;j<4;j++){
441  temp.v[i][j]=_mm_set1_pd(m2(i,j));
442  }
443  }
444  // Preprocessor macro for multiplying two __m128d elements together
445 #define MUL(i,j,k) _mm_mul_pd(GetV((i),(j)),temp.v[(j)][(k)])
446 
447  return DMatrix4x4(_mm_add_pd(MUL(0,0,0),
448  _mm_add_pd(MUL(0,1,0),
449  _mm_add_pd(MUL(0,2,0),
450  MUL(0,3,0)))),
451  _mm_add_pd(MUL(0,0,1),
452  _mm_add_pd(MUL(0,1,1),
453  _mm_add_pd(MUL(0,2,1),
454  MUL(0,3,1)))),
455  _mm_add_pd(MUL(0,0,2),
456  _mm_add_pd(MUL(0,1,2),
457  _mm_add_pd(MUL(0,2,2),
458  MUL(0,3,2)))),
459  _mm_add_pd(MUL(0,0,3),
460  _mm_add_pd(MUL(0,1,3),
461  _mm_add_pd(MUL(0,2,3),
462  MUL(0,3,3)))),
463  _mm_add_pd(MUL(1,0,0),
464  _mm_add_pd(MUL(1,1,0),
465  _mm_add_pd(MUL(1,2,0),
466  MUL(1,3,0)))),
467  _mm_add_pd(MUL(1,0,1),
468  _mm_add_pd(MUL(1,1,1),
469  _mm_add_pd(MUL(1,2,1),
470  MUL(1,3,1)))),
471  _mm_add_pd(MUL(1,0,2),
472  _mm_add_pd(MUL(1,1,2),
473  _mm_add_pd(MUL(1,2,2),
474  MUL(1,3,2)))),
475  _mm_add_pd(MUL(1,0,3),
476  _mm_add_pd(MUL(1,1,3),
477  _mm_add_pd(MUL(1,2,3),
478  MUL(1,3,3)))));
479  }
480 
481 
482  void Print(){
483  cout << "DMatrix4x4:" <<endl;
484  cout << " | 0 | 1 | 2 | 3 |" <<endl;
485  cout << "----------------------------------------------------------" <<endl;
486 
487  for (unsigned int i=0;i<4;i++){
488  for (unsigned int j=0;j<4;j++){
489  cout << mA[j].d[i] <<" ";
490  }
491  cout << endl;
492  }
493  }
494 
495  private:
496  union dvec{
497  __m128d v[2];
498  double d[4];
499  };
500  ALIGNED_16_BLOCK(union dvec, 4, mA)
501 };
502 #endif
#define F(x, y, z)
DMatrix4x4(const DMatrix2x2 &m1, const DMatrix2x2 &m2, const DMatrix2x2 &m3, const DMatrix2x2 &m4)
Definition: DMatrix4x4.h:39
DMatrix2x2 Invert()
Definition: DMatrix2x2.h:64
#define ALIGNED_16_BLOCK_PTR(TYPE, NUM, PTR)
Definition: align_16.h:24
void Print()
Definition: DMatrix4x4.h:203
DMatrix4x4 operator-(const DMatrix4x4 &m2)
Definition: DMatrix4x4.h:99
#define ALIGNED_16_BLOCK_WITH_PTR(TYPE, NUM, PTR)
Definition: align_16.h:27
TMatrixDSym DMatrixDSym
Definition: DMatrixDSym.h:13
DMatrix4x4 Transpose()
Definition: DMatrix4x4.h:174
DMatrix4x4 operator-()
Definition: DMatrix4x4.h:81
DMatrix4x4(double c11, double c12, double c13, double c14, double c21, double c22, double c23, double c24, double c31, double c32, double c33, double c34, double c41, double c42, double c43, double c44)
Definition: DMatrix4x4.h:16
double mA[4][4]
Definition: DMatrix4x4.h:219
#define H(x, y, z)
double & operator()(int row, int col)
Definition: DMatrix4x4.h:61
DMatrix4x4 Invert()
Definition: DMatrix4x4.h:161
DMatrix4x1 operator*(const DMatrix4x1 &m2)
Definition: DMatrix4x4.h:150
DMatrix4x4()
Definition: DMatrix4x4.h:9
DMatrix4x4 operator+(const DMatrix4x4 &m2)
Definition: DMatrix4x4.h:89
bool IsPosDef()
Definition: DMatrix4x4.h:195
DMatrix4x2 operator*(const DMatrix4x2 &m2)
Definition: DMatrix4x4.h:132
DMatrix4x4 & operator=(const DMatrix4x4 &m1)
Definition: DMatrix4x4.h:68
DMatrix4x4 operator*(const DMatrix4x4 &m2)
Definition: DMatrix4x4.h:108
#define ALIGNED_16_BLOCK(TYPE, NUM, PTR)
Definition: align_16.h:19
DMatrixDSym GetSub(unsigned int lowerBound, unsigned int upperBound)
Definition: DMatrix4x4.h:184
for(auto locVertexInfo:dStepVertexInfos)
TH2F * temp
#define G(x, y, z)