INMOST Mathematical Modelling Toolkit
inmost_dense.h
Go to the documentation of this file.
1 #ifndef INMOST_DENSE_INCLUDED
2 #define INMOST_DENSE_INCLUDED
3 #include "inmost_common.h"
4 #if defined(USE_AUTODIFF)
5 #include "inmost_expression.h"
6 #endif
7
8 // Matrix with n columns and m rows
9 // __m__
10 // | |
11 // n| |
12 // |_____|
13 //
14 // todo:
15 // 1. expression templates for operations
16 // (???) how to for multiplication?
17 // 2. (ok) template matrix type for AD variables
18
19
20 namespace INMOST
21 {
22
23
24  template<class A, class B> struct Promote;
26 #if defined(USE_AUTODIFF)
27  template<> struct Promote<INMOST_DATA_REAL_TYPE, variable> {typedef variable type;};
28  template<> struct Promote<variable, INMOST_DATA_REAL_TYPE> {typedef variable type;};
29  template<> struct Promote<variable, variable> {typedef variable type;};
30 #else
32 #endif
33
34  template<typename Var>
35  class Matrix
36  {
37  public:
38  typedef unsigned enumerator;
39  protected:
42
43
44  static Var sign_func(const Var & a, const Var & b) {return (b >= 0.0 ? 1.0 : -1.0)*fabs(a);}
45  static Var max_func(const Var & x, const Var & y) { return x > y ? x : y; }
46  static Var pythag(const Var & a, const Var & b)
47  {
48  Var at = fabs(a), bt = fabs(b), ct, result;
49  if (at > bt) { ct = bt / at; result = sqrt(at) * sqrt(at + ct * bt); }
50  else if (bt > 0.0) { ct = at / bt; result = sqrt(bt) * sqrt(bt + ct * at); }
51  else result = 0.0;
52  return result;
53  }
54  public:
55  void Swap(Matrix & b)
56  {
57  space.swap(b.space);
58  std::swap(n,b.n);
59  std::swap(m,b.m);
60  }
61  bool SVD(Matrix & U, Matrix & Sigma, Matrix & VT)
62  {
63  int flag, i, its, j, jj, k, l, nm;
64  Var c, f, h, s, x, y, z;
65  Var anorm = 0.0, g = 0.0, scale = 0.0;
66  if (n < m)
67  {
68  bool success = Transpose().SVD(U,Sigma,VT);
69  if( success )
70  {
71  U.Swap(VT);
72  U = U.Transpose();
73  VT = VT.Transpose();
74  return true;
75  }
76  else return false;
77  } //m <= n
78  array<Var> _rv1(m);
79  shell<Var> rv1(_rv1);
80  U = (*this);
81  Sigma.Resize(m,m);
82  Sigma.Zero();
83  VT.Resize(m,m);
84
85  std::swap(n,m); //this how original algorithm takes it
86  // Householder reduction to bidiagonal form
87  for (i = 0; i < n; i++)
88  {
89  // left-hand reduction
90  l = i + 1;
91  rv1[i] = scale * g;
92  g = s = scale = 0.0;
93  if (i < m)
94  {
95  for (k = i; k < m; k++) scale += fabs(U[k][i]);
96  if (get_value(scale))
97  {
98  for (k = i; k < m; k++)
99  {
100  U[k][i] /= scale;
101  s += U[k][i] * U[k][i];
102  }
103  f = U[i][i];
104  g = -sign_func(sqrt(s), f);
105  h = f * g - s;
106  U[i][i] = f - g;
107  if (i != n - 1)
108  {
109  for (j = l; j < n; j++)
110  {
111  for (s = 0.0, k = i; k < m; k++) s += U[k][i] * U[k][j];
112  f = s / h;
113  for (k = i; k < m; k++) U[k][j] += f * U[k][i];
114  }
115  }
116  for (k = i; k < m; k++) U[k][i] *= scale;
117  }
118  }
119  Sigma[i][i] = scale * g;
120  // right-hand reduction
121  g = s = scale = 0.0;
122  if (i < m && i != n - 1)
123  {
124  for (k = l; k < n; k++) scale += fabs(U[i][k]);
125  if (get_value(scale))
126  {
127  for (k = l; k < n; k++)
128  {
129  U[i][k] = U[i][k]/scale;
130  s += U[i][k] * U[i][k];
131  }
132  f = U[i][l];
133  g = -sign_func(sqrt(s), f);
134  h = f * g - s;
135  U[i][l] = f - g;
136  for (k = l; k < n; k++) rv1[k] = U[i][k] / h;
137  if (i != m - 1)
138  {
139  for (j = l; j < m; j++)
140  {
141  for (s = 0.0, k = l; k < n; k++) s += U[j][k] * U[i][k];
142  for (k = l; k < n; k++) U[j][k] += s * rv1[k];
143  }
144  }
145  for (k = l; k < n; k++) U[i][k] *= scale;
146  }
147  }
148  anorm = max_func(anorm,fabs(Sigma[i][i]) + fabs(rv1[i]));
149  }
150
151  // accumulate the right-hand transformation
152  for (i = n - 1; i >= 0; i--)
153  {
154  if (i < n - 1)
155  {
156  if (get_value(g))
157  {
158  for (j = l; j < n; j++) VT[j][i] = ((U[i][j] / U[i][l]) / g);
159  // double division to avoid underflow
160  for (j = l; j < n; j++)
161  {
162  for (s = 0.0, k = l; k < n; k++) s += U[i][k] * VT[k][j];
163  for (k = l; k < n; k++) VT[k][j] += s * VT[k][i];
164  }
165  }
166  for (j = l; j < n; j++) VT[i][j] = VT[j][i] = 0.0;
167  }
168  VT[i][i] = 1.0;
169  g = rv1[i];
170  l = i;
171  }
172
173  // accumulate the left-hand transformation
174  for (i = n - 1; i >= 0; i--)
175  {
176  l = i + 1;
177  g = Sigma[i][i];
178  if (i < n - 1)
179  for (j = l; j < n; j++)
180  U[i][j] = 0.0;
181  if (get_value(g))
182  {
183  g = 1.0 / g;
184  if (i != n - 1)
185  {
186  for (j = l; j < n; j++)
187  {
188  for (s = 0.0, k = l; k < m; k++) s += (U[k][i] * U[k][j]);
189  f = (s / U[i][i]) * g;
190  for (k = i; k < m; k++) U[k][j] += f * U[k][i];
191  }
192  }
193  for (j = i; j < m; j++) U[j][i] = U[j][i]*g;
194  }
195  else for (j = i; j < m; j++) U[j][i] = 0.0;
196  U[i][i] += 1;
197  }
198
199  // diagonalize the bidiagonal form
200  for (k = n - 1; k >= 0; k--)
201  {// loop over singular values
202  for (its = 0; its < 30; its++)
203  {// loop over allowed iterations
204  flag = 1;
205  for (l = k; l >= 0; l--)
206  {// test for splitting
207  nm = l - 1;
208  if (fabs(rv1[l]) + anorm == anorm)
209  {
210  flag = 0;
211  break;
212  }
213  if (fabs(Sigma[nm][nm]) + anorm == anorm)
214  break;
215  }
216  if (flag)
217  {
218  c = 0.0;
219  s = 1.0;
220  for (i = l; i <= k; i++)
221  {
222  f = s * rv1[i];
223  if (fabs(f) + anorm != anorm)
224  {
225  g = Sigma[i][i];
226  h = pythag(f, g);
227  Sigma[i][i] = h;
228  h = 1.0 / h;
229  c = g * h;
230  s = (- f * h);
231  for (j = 0; j < m; j++)
232  {
233  y = U[j][nm];
234  z = U[j][i];
235  U[j][nm] = (y * c + z * s);
236  U[j][i] = (z * c - y * s);
237  }
238  }
239  }
240  }
241  z = Sigma[k][k];
242  if (l == k)
243  {// convergence
244  if (z < 0.0)
245  {// make singular value nonnegative
246  Sigma[k][k] = -z;
247  for (j = 0; j < n; j++) VT[j][k] = -VT[j][k];
248  }
249  break;
250  }
251  if (its >= 30)
252  {
253  std::cout << "No convergence after " << its << " iterations" << std::endl;
254  std::swap(n,m);
255  return false;
256  }
257  // shift from bottom 2 x 2 minor
258  x = Sigma[l][l];
259  nm = k - 1;
260  y = Sigma[nm][nm];
261  g = rv1[nm];
262  h = rv1[k];
263  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
264  g = pythag(f, 1.0);
265  f = ((x - z) * (x + z) + h * ((y / (f + sign_func(g, f))) - h)) / x;
266  // next QR transformation
267  c = s = 1.0;
268  for (j = l; j <= nm; j++)
269  {
270  i = j + 1;
271  g = rv1[i];
272  y = Sigma[i][i];
273  h = s * g;
274  g = c * g;
275  z = pythag(f, h);
276  rv1[j] = z;
277  c = f / z;
278  s = h / z;
279  f = x * c + g * s;
280  g = g * c - x * s;
281  h = y * s;
282  y = y * c;
283  for (jj = 0; jj < n; jj++)
284  {
285  x = VT[jj][j];
286  z = VT[jj][i];
287  VT[jj][j] = (x * c + z * s);
288  VT[jj][i] = (z * c - x * s);
289  }
290  z = pythag(f, h);
291  Sigma[j][j] = z;
292  if (z)
293  {
294  z = 1.0 / z;
295  c = f * z;
296  s = h * z;
297  }
298  f = (c * g) + (s * y);
299  x = (c * y) - (s * g);
300  for (jj = 0; jj < m; jj++)
301  {
302  y = U[jj][j];
303  z = U[jj][i];
304  U[jj][j] = (y * c + z * s);
305  U[jj][i] = (z * c - y * s);
306  }
307  }
308  rv1[l] = 0.0;
309  rv1[k] = f;
310  Sigma[k][k] = x;
311  }
312  }
313  std::swap(n,m);
314  return true;
315  }
316
317  Matrix(Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {}
318  Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {}
319  Matrix(const Matrix & other) : space(other.n*other.m), n(other.n), m(other.m)
320  {
321  for(enumerator i = 0; i < n*m; ++i)
322  space[i] = other.space[i];
323  }
324  ~Matrix() {}
325  void Resize(enumerator nrows, enumerator mcols)
326  {
327  if( space.size() != mcols*nrows )
328  space.resize(mcols*nrows);
329  n = nrows;
330  m = mcols;
331  }
332  Matrix & operator =(Matrix const & other)
333  {
334  space.resize(other.n*other.m);
335  for(enumerator i = 0; i < other.n*other.m; ++i)
336  space[i] = other.space[i];
337  n = other.n;
338  m = other.m;
339  return *this;
340  }
341  // i is in [0,n] - row index
342  // j is in [0,m] - column index
344  {
345  assert(i >= 0 && i < n);
346  assert(j >= 0 && j < m);
347  assert(i*m+j < n*m); //overflow check?
348  return space[i*m+j];
349  }
351  {
352  assert(i >= 0 && i < n);
353  assert(j >= 0 && j < m);
354  assert(i*m+j < n*m); //overflow check?
355  return space[i*m+j];
356  }
358  {
359  Matrix ret(n,m);
360  for(enumerator k = 0; k < n*m; ++k) ret.space[k] = -space[k];
361  return ret;
362  }
363  template<typename typeB>
365  {
366  assert(Rows() == other.Rows());
367  assert(Cols() == other.Cols());
368  Matrix<typename Promote<Var,typeB>::type> ret(n,m); //check RVO
369  for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]-other.space[k];
370  return ret;
371  }
372  Matrix & operator-=(const Matrix & other)
373  {
374  assert(n == other.n);
375  assert(m == other.m);
376  for(enumerator k = 0; k < n*m; ++k) space[k] -= other.space[k];
377  return *this;
378  }
379  template<typename typeB>
381  {
382  assert(Rows() == other.Rows());
383  assert(Cols() == other.Cols());
384  Matrix<typename Promote<Var,typeB>::type> ret(n,m); //check RVO
385  for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]+other.space[k];
386  return ret;
387  }
388  Matrix & operator+=(const Matrix & other)
389  {
390  assert(n == other.n);
391  assert(m == other.m);
392  for(enumerator k = 0; k < n*m; ++k) space[k] += other.space[k];
393  return *this;
394  }
395  template<typename typeB>
397  {
398  Matrix<typename Promote<Var,typeB>::type> ret(n,m); //check RVO
399  for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]*coef;
400  return ret;
401  }
402  Matrix & operator*=(Var coef)
403  {
404  for(enumerator k = 0; k < n*m; ++k) space[k] *= coef;
405  return *this;
406  }
407  template<typename typeB>
409  {
410  Matrix<typename Promote<Var,typeB>::type> ret(n,m); //check RVO
411  for(enumerator k = 0; k < n*m; ++k) ret.space[k] = space[k]/coef;
412  return ret;
413  }
414  Matrix & operator/=(Var coef)
415  {
416  for(enumerator k = 0; k < n*m; ++k) space[k] /= coef;
417  return *this;
418  }
419  template<typename typeB>
421  {
422  assert(Cols() == other.Rows());
423  Matrix<typename Promote<Var,typeB>::type> ret(Rows(),other.Cols()); //check RVO
424  for(enumerator i = 0; i < Rows(); ++i) //loop rows
425  {
426  for(enumerator j = 0; j < other.Cols(); ++j) //loop columns
427  {
428  typename Promote<Var,typeB>::type tmp = 0.0;
429  for(enumerator k = 0; k < Cols(); ++k)
430  tmp += (*this)(i,k)*other(k,j);
431  ret(i,j) = tmp;
432  }
433  }
434  return ret;
435  }
438  template<typename typeB>
440  {
441  std::pair<Matrix<typeB>,bool> other_inv = other.Invert();
442  assert(other_inv.second);
443  assert(Cols() == other_inv.Rows());
444  Matrix<typename Promote<Var,typeB>::type> ret(n,other.m); //check RVO
445  for(enumerator i = 0; i < Rows(); ++i) //loop rows
446  {
447  for(enumerator j = 0; j < other_inv.Cols(); ++j) //loop columns
448  {
449  typename Promote<Var,typeB>::type tmp = 0.0;
450  for(enumerator k = 0; k < Cols(); ++k)
451  tmp += (*this)(i,k)*other_inv.first(k,j);
452  ret(i,j) = tmp;
453  }
454  }
455  return ret;
456  }
458  {
459  Matrix ret(m,n);
460  for(enumerator i = 0; i < n; ++i)
461  {
462  for(enumerator j = 0; j < m; ++j)
463  {
464  ret(j,i) = (*this)(i,j);
465  }
466  }
467  return ret;
468  }
469  std::pair<Matrix,bool> Invert() const
470  {
471  std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true);
472  Matrix At = Transpose(); //m by n matrix
473  Matrix AtB = At*Unit(n); //m by n matrix
474  Matrix AtA = At*(*this); //m by m matrix
475  enumerator * order = new enumerator [m];
476  for(enumerator i = 0; i < m; ++i) order[i] = i;
477  for(enumerator i = 0; i < m; i++)
478  {
479  enumerator maxk = i, maxq = i, temp2;
480  Var max, temp;
481  max = fabs(AtA(maxk,maxq));
482  //Find best pivot
483  for(enumerator k = i; k < m; k++) // over rows
484  {
485  for(enumerator q = i; q < m; q++) // over columns
486  {
487  if( fabs(AtA(k,q)) > max )
488  {
489  max = fabs(AtA(k,q));
490  maxk = k;
491  maxq = q;
492  }
493  }
494  }
495  //Exchange rows
496  if( maxk != i )
497  {
498  for(enumerator q = 0; q < m; q++) // over columns of A
499  {
500  temp = AtA(maxk,q);
501  AtA(maxk,q) = AtA(i,q);
502  AtA(i,q) = temp;
503  }
504  //exchange rhs
505  for(enumerator q = 0; q < n; q++) // over columns of B
506  {
507  temp = AtB(maxk,q);
508  AtB(maxk,q) = AtB(i,q);
509  AtB(i,q) = temp;
510  }
511  }
512  //Exchange columns
513  if( maxq != i )
514  {
515  for(enumerator k = 0; k < m; k++) //over rows
516  {
517  temp = AtA(k,maxq);
518  AtA(k,maxq) = AtA(k,i);
519  AtA(k,i) = temp;
520  }
521  //remember order in sol
522  {
523  temp2 = order[maxq];
524  order[maxq] = order[i];
525  order[i] = temp2;
526  }
527  }
528
529  if( fabs(AtA(i,i)) < 1.0e-54 )
530  {
531  bool ok = true;
532  for(enumerator k = 0; k < n; k++) // over columns of B
533  {
534  if( fabs(AtB(i,k)/1.0e-54) > 1 )
535  {
536  ok = false;
537  break;
538  }
539  }
540  if( ok ) AtA(i,i) = AtA(i,i) < 0.0 ? - 1.0e-12 : 1.0e-12;
541  else
542  {
543  ret.second = false;
544  return ret;
545  }
546  }
547  for(enumerator k = i+1; k < m; k++)
548  {
549  AtA(i,k) /= AtA(i,i);
550  AtA(k,i) /= AtA(i,i);
551  }
552  for(enumerator k = i+1; k < m; k++)
553  for(enumerator q = i+1; q < m; q++)
554  {
555  AtA(k,q) -= AtA(k,i) * AtA(i,i) * AtA(i,q);
556  }
557  for(enumerator k = 0; k < n; k++)
558  {
559  for(enumerator j = i+1; j < m; j++) //iterate over columns of L
560  {
561  AtB(j,k) -= AtB(i,k) * AtA(j,i);
562  }
563  AtB(i,k) /= AtA(i,i);
564  }
565  }
566
567  for(enumerator k = 0; k < n; k++)
568  {
569  for(enumerator i = m; i-- > 0; ) //iterate over rows of U
570  for(enumerator j = i+1; j < m; j++)
571  {
572  AtB(i,k) -= AtB(j,k) * AtA(i,j);
573  }
574  for(enumerator i = 0; i < m; i++)
575  ret.first(order[i],k) = AtB(i,k);
576  }
577  delete [] order;
578  return ret;
579  }
580  void Zero()
581  {
582  for(enumerator i = 0; i < n*m; ++i) space[i] = 0.0;
583  }
584  Var Trace() const
585  {
586  assert(n == m);
587  Var ret = 0.0;
588  for(enumerator i = 0; i < n; ++i) ret += (*this)(i,i);
589  return ret;
590  }
591  Var * data() {return space.data();}
592  const Var * data() const {return space.data();}
593  enumerator Rows() const {return n;}
594  enumerator Cols() const {return m;}
595  void Print(INMOST_DATA_REAL_TYPE threshold = 1.0e-10) const
596  {
597  for(enumerator k = 0; k < n; ++k)
598  {
599  for(enumerator l = 0; l < m; ++l)
600  {
601  if( fabs(get_value((*this)(k,l))) > threshold )
602  std::cout << (*this)(k,l);
603  else
604  std::cout << 0;
605  std::cout << " ";
606  }
607  std::cout << std::endl;
608  }
609  }
610  bool isSymmetric() const
611  {
612  if( n != m ) return false;
613  for(enumerator k = 0; k < n; ++k)
614  {
615  for(enumerator l = k+1; l < n; ++l)
616  if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 )
617  return false;
618  }
619  return true;
620  }
621  Var DotProduct(const Matrix & other) const
622  {
623  assert(n == other.n);
624  assert(m == other.m);
625  Var ret = 0.0;
626  for(enumerator i = 0; i < n*m; ++i) ret += space[i]*other.space[i];
627  return ret;
628  }
630  {
631  Var ret = 0;
632  for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i];
633  return sqrt(ret);
634  }
635  static Matrix<Var> FromTensor(Var * K, enumerator size)
636  {
637  Matrix<Var> Kc(3,3);
638  switch(size)
639  {
640  case 1: //scalar permeability tensor
641  Kc.Zero();
642  Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0];
643  case 3:
644  Kc.Zero(); //diagonal permeability tensor
645  Kc(0,0) = K[0]; //KXX
646  Kc(1,1) = K[1]; //KYY
647  Kc(2,2) = K[2]; //KZZ
648  break;
649  case 6: //symmetric permeability tensor
650  Kc(0,0) = K[0]; //KXX
651  Kc(0,1) = Kc(1,0) = K[1]; //KXY
652  Kc(0,2) = Kc(2,0) = K[2]; //KXZ
653  Kc(1,1) = K[3]; //KYY
654  Kc(1,2) = Kc(2,1) = K[4]; //KYZ
655  Kc(2,2) = K[5]; //KZZ
656  break;
657  case 9: //full permeability tensor
658  Kc(0,0) = K[0]; //KXX
659  Kc(0,1) = K[1]; //KXY
660  Kc(0,2) = K[2]; //KXZ
661  Kc(1,0) = K[3]; //KYX
662  Kc(1,1) = K[4]; //KYY
663  Kc(1,2) = K[5]; //KYZ
664  Kc(2,0) = K[6]; //KZX
665  Kc(2,1) = K[7]; //KZY
666  Kc(2,2) = K[8]; //KZZ
667  break;
668  }
669  return Kc;
670  }
671  //retrive vector in matrix form from array
672  static Matrix<Var> FromVector(Var * n, enumerator size)
673  {
674  return Matrix(n,size,1);
675  }
676  //create diagonal matrix from array
677  static Matrix<Var> FromDiagonal(Var * r, enumerator size)
678  {
679  Matrix ret(size,size);
680  ret.Zero();
681  for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k];
682  return ret;
683  }
684  //create diagonal matrix from array of inversed values
686  {
687  Matrix ret(size,size);
688  ret.Zero();
689  for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
690  return ret;
691  }
692  //Unit matrix
693  static Matrix Unit(enumerator pn)
694  {
695  Matrix ret(pn,pn);
696  ret.Zero();
697  for(enumerator i = 0; i < pn; ++i) ret(i,i) = 1.0;
698  return ret;
699  }
700  };
701
702
703  typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix
704 #if defined(USE_AUTODIFF)
705  typedef Matrix<variable> vMatrix; //shortcut for matrix with variations
706 #endif
707 }
708
709
710 #endif //INMOST_DENSE_INCLUDED
__INLINE INMOST::unary_custom_variable< INMOST::abs_expression< typename A::Var >, A > fabs(INMOST::shell_dynamic_variable< typename A::Var, A > const &Arg)
Matrix< typename Promote< Var, typeB >::type > operator*(typeB coef) const
Definition: inmost_dense.h:396
__INLINE INMOST_DATA_REAL_TYPE get_value(INMOST::shell_expression< A > const &Arg)
bool isSymmetric() const
Definition: inmost_dense.h:610
enumerator n
Definition: inmost_dense.h:41
Matrix< typename Promote< Var, typeB >::type > operator/(const Matrix< typeB > &other) const
Definition: inmost_dense.h:439
Var & operator()(enumerator i, enumerator j)
Definition: inmost_dense.h:343
static Matrix< Var > FromTensor(Var *K, enumerator size)
Definition: inmost_dense.h:635
Matrix & operator*=(Var coef)
Definition: inmost_dense.h:402
static Var pythag(const Var &a, const Var &b)
Definition: inmost_dense.h:46
Matrix Transpose() const
Definition: inmost_dense.h:457
Var DotProduct(const Matrix &other) const
Definition: inmost_dense.h:621
void Resize(enumerator nrows, enumerator mcols)
Definition: inmost_dense.h:325
#define INMOST_DATA_REAL_TYPE
__INLINE INMOST::sqrt_expression< A > sqrt(INMOST::shell_expression< A > const &Arg)
Matrix< typename Promote< Var, typeB >::type > operator*(const Matrix< typeB > &other) const
Definition: inmost_dense.h:420
Matrix & operator=(Matrix const &other)
Definition: inmost_dense.h:332
Var Trace() const
Definition: inmost_dense.h:584
void swap(array< element > &other)
Definition: container.hpp:339
Matrix(const Matrix &other)
Definition: inmost_dense.h:319
enumerator Rows() const
Definition: inmost_dense.h:593
Matrix & operator+=(const Matrix &other)
Definition: inmost_dense.h:388
void resize(size_type n, element c=element())
Definition: container.hpp:312
Matrix(Var *pspace, enumerator pn, enumerator pm)
Definition: inmost_dense.h:317
const Var * data() const
Definition: inmost_dense.h:592
Matrix< typename Promote< Var, typeB >::type > operator/(typeB coef) const
Definition: inmost_dense.h:408
bool SVD(Matrix &U, Matrix &Sigma, Matrix &VT)
Definition: inmost_dense.h:61
enumerator Cols() const
Definition: inmost_dense.h:594
void Print(INMOST_DATA_REAL_TYPE threshold=1.0e-10) const
Definition: inmost_dense.h:595
static Matrix< Var > FromDiagonalInverse(Var *r, enumerator size)
Definition: inmost_dense.h:685
enumerator m
Definition: inmost_dense.h:41
static Matrix< Var > FromDiagonal(Var *r, enumerator size)
Definition: inmost_dense.h:677
__INLINE size_type size() const
Definition: container.hpp:330
static Matrix Unit(enumerator pn)
Definition: inmost_dense.h:693
static Var max_func(const Var &x, const Var &y)
Definition: inmost_dense.h:45
Matrix< typename Promote< Var, typeB >::type > operator+(const Matrix< typeB > &other) const
Definition: inmost_dense.h:380
void Swap(Matrix &b)
Definition: inmost_dense.h:55
Matrix(enumerator pn, enumerator pm)
Definition: inmost_dense.h:318
Var operator()(enumerator i, enumerator j) const
Definition: inmost_dense.h:350
Matrix< typename Promote< Var, typeB >::type > operator-(const Matrix< typeB > &other) const
Definition: inmost_dense.h:364
Matrix & operator-=(const Matrix &other)
Definition: inmost_dense.h:372
Matrix operator-() const
Definition: inmost_dense.h:357
Matrix< INMOST_DATA_REAL_TYPE > rMatrix
Definition: inmost_dense.h:703
unsigned enumerator
Definition: inmost_dense.h:38
static Matrix< Var > FromVector(Var *n, enumerator size)
Definition: inmost_dense.h:672
Matrix & operator/=(Var coef)
Definition: inmost_dense.h:414
__INLINE element * data()
Definition: container.hpp:159
std::pair< Matrix, bool > Invert() const
Definition: inmost_dense.h:469
Matrix< variable > vMatrix
Definition: inmost_dense.h:705
static Var sign_func(const Var &a, const Var &b)
Definition: inmost_dense.h:44
array< Var > space
Definition: inmost_dense.h:40