INMOST
A toolkit for distributed mathematical modeling
inmost_dense.h
1 #ifndef INMOST_DENSE_INCLUDED
2 #define INMOST_DENSE_INCLUDED
3 #include "inmost_common.h"
4 #include "inmost_expression.h"
5 #include <iomanip>
6 #include <stdarg.h>
7 
8 namespace INMOST
9 {
10 
12  template<class A> struct SelfPromote;
13  template<class A> struct ComplexType;
14  template<class A, class B> struct Promote;
15  template<> struct ComplexType<INMOST_DATA_INTEGER_TYPE> { typedef INMOST_DATA_INTEGER_TYPE type; };
16  template<> struct ComplexType<INMOST_DATA_REAL_TYPE> { typedef INMOST_DATA_REAL_TYPE type; };
17  template<> struct ComplexType<INMOST_DATA_CPLX_TYPE> { typedef INMOST_DATA_REAL_TYPE type; };
18  template<typename T> inline typename ComplexType<T>::type real_part(T const& val) { return std::real(val); }
19  template<> struct SelfPromote<INMOST_DATA_INTEGER_TYPE> { typedef INMOST_DATA_INTEGER_TYPE type; };
20  template<> struct SelfPromote<INMOST_DATA_REAL_TYPE> { typedef INMOST_DATA_REAL_TYPE type; };
21  template<> struct SelfPromote<INMOST_DATA_CPLX_TYPE> { typedef INMOST_DATA_CPLX_TYPE type; };
22  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_INTEGER_TYPE type;};
23  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
24  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, INMOST_DATA_CPLX_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
25  template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
26  template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
27  template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_CPLX_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
28  template<> struct Promote<INMOST_DATA_CPLX_TYPE, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
29  template<> struct Promote<INMOST_DATA_CPLX_TYPE, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
30  template<> struct Promote<INMOST_DATA_CPLX_TYPE, INMOST_DATA_CPLX_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
31 #if defined(USE_FP64)
32  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, float> {typedef INMOST_DATA_REAL_TYPE type;};
33  template<> struct Promote<float, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
34  template<> struct Promote<INMOST_DATA_REAL_TYPE, float> {typedef INMOST_DATA_REAL_TYPE type;};
35  template<> struct Promote<float, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
36  template<> struct Promote<INMOST_DATA_CPLX_TYPE, float> {typedef INMOST_DATA_CPLX_TYPE type;};
37  template<> struct Promote<float, INMOST_DATA_CPLX_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
38  template<> struct Promote<float, float> { typedef float type; };
39 #else
40  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, double> {typedef INMOST_DATA_REAL_TYPE type;};
41  template<> struct Promote<double, INMOST_DATA_INTEGER_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
42  template<> struct Promote<INMOST_DATA_REAL_TYPE, double> {typedef INMOST_DATA_REAL_TYPE type;};
43  template<> struct Promote<double, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
44  template<> struct Promote<INMOST_DATA_CPLX_TYPE, double> {typedef INMOST_DATA_CPLX_TYPE type;};
45  template<> struct Promote<double, INMOST_DATA_CPLX_TYPE> {typedef INMOST_DATA_CPLX_TYPE type;};
46  template<> struct Promote<double, double> { typedef double type; };
47 #endif
48 #if defined(USE_AUTODIFF)
49  template<> struct SelfPromote<unknown> { typedef variable type; };
50  template<> struct SelfPromote<variable> { typedef variable type; };
51  template<> struct SelfPromote<value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
52  template<> struct SelfPromote<multivar_expression_reference> { typedef variable type; };
53  template<> struct SelfPromote<hessian_multivar_expression_reference> { typedef hessian_variable type; };
54  template<> struct SelfPromote<hessian_variable> { typedef hessian_variable type; };
55  template<> struct ComplexType<unknown> { typedef unknown type; };
56  template<> struct ComplexType<variable> { typedef variable type; };
57  template<> struct ComplexType<value_reference> { typedef value_reference type; };
58  template<> struct ComplexType<multivar_expression_reference> { typedef variable type; };
59  template<> struct ComplexType<hessian_multivar_expression_reference> { typedef hessian_variable type; };
60  template<> struct ComplexType<hessian_variable> { typedef hessian_variable type; };
61  template<> struct ComplexType< std::complex<unknown> > { typedef unknown type; };
62  template<> struct ComplexType< std::complex<variable> > { typedef variable type; };
63  template<> struct ComplexType< std::complex<value_reference> > { typedef value_reference type; };
64  template<> struct ComplexType< std::complex<multivar_expression_reference> > { typedef multivar_expression_reference type; };
65  template<> struct ComplexType< std::complex<hessian_multivar_expression_reference> > { typedef hessian_multivar_expression_reference type; };
66  template<> struct ComplexType< std::complex<hessian_variable> > { typedef hessian_variable type; };
67  template<> inline typename ComplexType<unknown>::type real_part(unknown const& a) { return a; }
68  template<> inline typename ComplexType<variable>::type real_part(variable const& a) { return a; }
69  template<> inline typename ComplexType<value_reference>::type real_part(value_reference const& a) { return a; }
70  template<> inline typename ComplexType<multivar_expression_reference>::type real_part(multivar_expression_reference const& a) { return a; }
71  template<> inline typename ComplexType<hessian_multivar_expression_reference>::type real_part(hessian_multivar_expression_reference const& a) { return a; }
72  template<> inline typename ComplexType<hessian_variable>::type real_part(hessian_variable const& a) { return a; }
73  template<> inline typename ComplexType< std::complex<unknown> >::type real_part(std::complex<unknown> const& a) { return a.real(); }
74  template<> inline typename ComplexType< std::complex<variable> >::type real_part(std::complex<variable> const& a) { return a.real(); }
75  template<> inline typename ComplexType< std::complex<value_reference> >::type real_part(std::complex<value_reference> const& a) { return a.real(); }
76  template<> inline typename ComplexType< std::complex<multivar_expression_reference> >::type real_part(std::complex<multivar_expression_reference> const& a) { return a.real(); }
77  template<> inline typename ComplexType< std::complex<hessian_multivar_expression_reference> >::type real_part(std::complex<hessian_multivar_expression_reference> const& a) { return a.real(); }
78  template<> inline typename ComplexType< std::complex<hessian_variable> >::type real_part(std::complex<hessian_variable> const& a) { return a.real(); }
79 
80  //For INMOST_DATA_INTEGER_TYPE
81  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, unknown> {typedef variable type;};
82  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, variable> {typedef variable type;};
83  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
84  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, multivar_expression_reference> {typedef variable type;};
85  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, hessian_multivar_expression_reference> {typedef hessian_variable type;};
86  template<> struct Promote<INMOST_DATA_INTEGER_TYPE, hessian_variable> {typedef hessian_variable type;};
87  //For INMOST_DATA_REAL_TYPE
88  template<> struct Promote<INMOST_DATA_REAL_TYPE, unknown> {typedef variable type;};
89  template<> struct Promote<INMOST_DATA_REAL_TYPE, variable> {typedef variable type;};
90  template<> struct Promote<INMOST_DATA_REAL_TYPE, value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
91  template<> struct Promote<INMOST_DATA_REAL_TYPE, multivar_expression_reference> {typedef variable type;};
92  template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_multivar_expression_reference> {typedef hessian_variable type;};
93  template<> struct Promote<INMOST_DATA_REAL_TYPE, hessian_variable> {typedef hessian_variable type;};
94  //For INMOST_DATA_CPLX_TYPE
95  template<> struct Promote<INMOST_DATA_CPLX_TYPE, unknown> { typedef std::complex<variable> type; };
96  template<> struct Promote<INMOST_DATA_CPLX_TYPE, variable> { typedef std::complex<variable> type; };
97  template<> struct Promote<INMOST_DATA_CPLX_TYPE, value_reference> { typedef std::complex<INMOST_DATA_REAL_TYPE> type; };
98  template<> struct Promote<INMOST_DATA_CPLX_TYPE, multivar_expression_reference> { typedef std::complex<variable> type; };
99  template<> struct Promote<INMOST_DATA_CPLX_TYPE, hessian_multivar_expression_reference> { typedef std::complex<hessian_variable> type; };
100  template<> struct Promote<INMOST_DATA_CPLX_TYPE, hessian_variable> { typedef std::complex<hessian_variable> type; };
101  //for unknown
102  template<> struct Promote<unknown, INMOST_DATA_INTEGER_TYPE> {typedef variable type;};
103  template<> struct Promote<unknown, INMOST_DATA_REAL_TYPE> {typedef variable type;};
104  template<> struct Promote<unknown, unknown> {typedef variable type;};
105  template<> struct Promote<unknown, variable> {typedef variable type;};
106  template<> struct Promote<unknown, value_reference> { typedef variable type; };
107  template<> struct Promote<unknown, multivar_expression_reference> {typedef variable type;};
108  template<> struct Promote<unknown, hessian_multivar_expression_reference> {typedef hessian_variable type;};
109  template<> struct Promote<unknown, hessian_variable> {typedef hessian_variable type;};
110  //for value_reference
111  template<> struct Promote<value_reference, unknown> { typedef variable type; };
112  template<> struct Promote<value_reference, variable> { typedef variable type; };
113  template<> struct Promote<value_reference, value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
114  template<> struct Promote<value_reference, multivar_expression_reference> { typedef variable type; };
115  template<> struct Promote<value_reference, hessian_multivar_expression_reference> { typedef hessian_variable type; };
116  template<> struct Promote<value_reference, hessian_variable> { typedef hessian_variable type; };
117  //For multivar_expression_reference
118  template<> struct Promote<multivar_expression_reference, INMOST_DATA_INTEGER_TYPE> {typedef variable type;};
119  template<> struct Promote<multivar_expression_reference, INMOST_DATA_REAL_TYPE> {typedef variable type;};
120  template<> struct Promote<multivar_expression_reference, unknown> {typedef variable type;};
121  template<> struct Promote<multivar_expression_reference, variable> {typedef variable type;};
122  template<> struct Promote<multivar_expression_reference, value_reference> { typedef variable type; };
123  template<> struct Promote<multivar_expression_reference, multivar_expression_reference> {typedef variable type;};
124  template<> struct Promote<multivar_expression_reference, hessian_multivar_expression_reference> {typedef hessian_variable type;};
125  template<> struct Promote<multivar_expression_reference, hessian_variable> {typedef hessian_variable type;};
126  //For variable
127  template<> struct Promote<variable, INMOST_DATA_INTEGER_TYPE> {typedef variable type;};
128  template<> struct Promote<variable, INMOST_DATA_REAL_TYPE> {typedef variable type;};
129  template<> struct Promote<variable, unknown> {typedef variable type;};
130  template<> struct Promote<variable, variable> {typedef variable type;};
131  template<> struct Promote<variable, value_reference> { typedef variable type; };
132  template<> struct Promote<variable, multivar_expression_reference> {typedef variable type;};
133  template<> struct Promote<variable, hessian_multivar_expression_reference> {typedef hessian_variable type;};
134  template<> struct Promote<variable, hessian_variable> {typedef hessian_variable type;};
135  //For hessian_multivar_expression_reference
136  template<> struct Promote<hessian_multivar_expression_reference, INMOST_DATA_INTEGER_TYPE> {typedef hessian_variable type;};
137  template<> struct Promote<hessian_multivar_expression_reference, INMOST_DATA_REAL_TYPE> {typedef hessian_variable type;};
138  template<> struct Promote<hessian_multivar_expression_reference, unknown> {typedef hessian_variable type;};
139  template<> struct Promote<hessian_multivar_expression_reference, variable> {typedef hessian_variable type;};
140  template<> struct Promote<hessian_multivar_expression_reference, value_reference> { typedef hessian_variable type; };
141  template<> struct Promote<hessian_multivar_expression_reference, multivar_expression_reference> {typedef hessian_variable type;};
142  template<> struct Promote<hessian_multivar_expression_reference, hessian_multivar_expression_reference> {typedef hessian_variable type;};
143  template<> struct Promote<hessian_multivar_expression_reference, hessian_variable> {typedef hessian_variable type;};
144  //For hessian_variable
145  template<> struct Promote<hessian_variable, INMOST_DATA_INTEGER_TYPE> {typedef hessian_variable type;};
146  template<> struct Promote<hessian_variable, INMOST_DATA_REAL_TYPE> {typedef hessian_variable type;};
147  template<> struct Promote<hessian_variable, unknown> {typedef hessian_variable type;};
148  template<> struct Promote<hessian_variable, variable> {typedef hessian_variable type;};
149  template<> struct Promote<hessian_variable, value_reference> { typedef hessian_variable type; };
150  template<> struct Promote<hessian_variable, multivar_expression_reference> {typedef hessian_variable type;};
151  template<> struct Promote<hessian_variable, hessian_multivar_expression_reference> {typedef hessian_variable type;};
152  template<> struct Promote<hessian_variable, hessian_variable> {typedef hessian_variable type;};
153 #if defined(USE_FP64)
154  template<> struct Promote<unknown, float> {typedef variable type;};
155  template<> struct Promote<float, unknown> { typedef variable type; };
156  template<> struct Promote<value_reference, float> {typedef INMOST_DATA_REAL_TYPE type;};
157  template<> struct Promote<float, value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
158  template<> struct Promote<multivar_expression_reference, float> { typedef variable type; };
159  template<> struct Promote<float, multivar_expression_reference> { typedef variable type; };
160  template<> struct Promote<variable, float> {typedef variable type;};
161  template<> struct Promote<float, variable> { typedef variable type; };
162  template<> struct Promote<hessian_multivar_expression_reference, float> {typedef hessian_variable type;};
163  template<> struct Promote<float, hessian_multivar_expression_reference> { typedef hessian_variable type; };
164  template<> struct Promote<hessian_variable, float> {typedef hessian_variable type;};
165  template<> struct Promote<float, hessian_variable> { typedef hessian_variable type; };
166 #else
167  template<> struct Promote<unknown, double> {typedef variable type;};
168  template<> struct Promote<double, unknown> { typedef variable type; };
169  template<> struct Promote<value_reference, double> {typedef INMOST_DATA_REAL_TYPE type;};
170  template<> struct Promote<double, value_reference> { typedef INMOST_DATA_REAL_TYPE type; };
171  template<> struct Promote<multivar_expression_reference, double> { typedef variable type; };
172  template<> struct Promote<double, multivar_expression_reference> { typedef variable type; };
173  template<> struct Promote<variable, double> {typedef variable type;};
174  template<> struct Promote<double, variable> { typedef variable type; };
175  template<> struct Promote<hessian_multivar_expression_reference, double> { typedef hessian_variable type; };
176  template<> struct Promote<double, hessian_multivar_expression_reference> {typedef hessian_variable type;};
177  template<> struct Promote<hessian_variable, double> {typedef hessian_variable type;};
178  template<> struct Promote<double, hessian_variable> { typedef hessian_variable type; };
179 #endif
180 #endif
181 
182  class AbstractMatrixBase
183  {
184  protected:
185  static thread_private<Sparse::RowMerger> merger;
186  };
187 
188 
189  template<typename Var>
190  class AbstractMatrixReadOnly : public AbstractMatrixBase
191  {
192  private:
194  static Var sign_func(const Var& a, const Var& b) { return (real_part(get_value(b)) >= 0.0 ? fabs(a) : -fabs(a)); }
196  static INMOST_DATA_REAL_TYPE max_func(INMOST_DATA_REAL_TYPE x, INMOST_DATA_REAL_TYPE y) { return x > y ? x : y; }
198  static Var pythag(const Var& a, const Var& b)
199  {
200  Var at = fabs(a), bt = fabs(b), ct;
201  if (real_part(get_value(at)) > real_part(get_value(bt))) { ct = bt / at; return at * sqrt(1.0 + ct * ct); }
202  else if (real_part(get_value(bt)) > 0.0) { ct = at / bt; return bt * sqrt(1.0 + ct * ct); }
203  else return 0.0;
204  //return sqrt(a * conj(a) + b * conj(b));
205  }
206  public:
207  typedef unsigned enumerator;
210  virtual enumerator Rows() const = 0;
213  virtual enumerator Cols() const = 0;
219  virtual Var compute(enumerator i, enumerator j) const = 0;
222  bool CheckNans() const { return Matrix<Var>(*this).CheckNans(); }
225  bool CheckInfs() const { return Matrix<Var>(*this).CheckInfs();}
228  bool CheckNansInfs() const { return Matrix<Var>(*this).CheckNansInfs(); }
230  Var Det() const
231  {
232  assert(Rows() == Cols());
233  const double eps = 1.0e-13;
234  const enumerator n = Rows();
235  Matrix<Var> A(*this);
236  Matrix<Var> row_visited(n,1,-1);
237  Var ret = 1.0, coef;
238  double sign;
239  for (enumerator d = 0; d < n; ++d)
240  {
241  enumerator r = 0;
242  sign = 1;
243  // find unvisited row with non-zero value in column d
244  while(row_visited(r,0) != -1 || fabs(get_value(A(r, d))) < eps)
245  {
246  if(r == n-1) return Var(0);
247  if(row_visited(r,0) == -1) sign *= -1;
248  ++r;
249  }
250  row_visited(r,0) = d;
251  ret *= sign * A(r, d);
252  for (enumerator i = 0; i < n; ++i) if(row_visited(i,0) == -1)
253  {
254  coef = A(i, d) / A(r, d);
255  if(fabs(coef) > eps)
256  for (enumerator j = 0; j < n; ++j)
257  A(i, j) = A(i, j) - coef * A(r, j);
258  }
259  }
260  return ret;
261  }
273  bool SVD(AbstractMatrix<Var>& U, AbstractMatrix<Var>& Sigma, AbstractMatrix<Var>& V, bool order_singular_values = true, bool nonnegative = true) const;
276  bool cSVD(AbstractMatrix<Var>& U, AbstractMatrix<Var>& Sigma, AbstractMatrix<Var>& V) const;
279  //Matrix<Var> Transpose() const;
280  ConstMatrixTranspose<Var> Transpose() const { return ConstMatrixTranspose<Var>(*this); }
283  ConstMatrixConjugateTranspose<Var> ConjugateTranspose() const { return ConstMatrixConjugateTranspose<Var>(*this); }
286  ConstMatrixConjugate<Var> Conjugate() const { return ConstMatrixConjugate<Var>(*this); }
292  template<typename typeB>
293  Matrix<typename Promote<Var, typeB>::type>
294  CrossProduct(const AbstractMatrixReadOnly<typeB>& other) const;
299  template<typename typeB>
300  Matrix<typename Promote<Var, typeB>::type>
301  Transform(const AbstractMatrixReadOnly<typeB>& other) const;
305  template<typename typeB>
306  //Matrix<typename Promote<Var, typeB>::type>
307  MatrixDifference<Var,typeB>
308  operator-(const AbstractMatrixReadOnly<typeB>& other) const;
312  template<typename typeB>
313  //Matrix<typename Promote<Var, typeB>::type>
314  MatrixSum<Var,typeB>
315  operator+(const AbstractMatrixReadOnly<typeB>& other) const;
319  template<typename typeB>
320  //Matrix<typename Promote<Var, typeB>::type>
321  MatrixMul<Var,typeB, typename Promote<Var,typeB>::type>
322  operator*(const AbstractMatrixReadOnly<typeB>& other) const;
324  MatrixUnaryMinus<Var> operator-() const { return MatrixUnaryMinus<Var>(*this); }
328  template<typename typeB>
329  //Matrix<typename Promote<Var, typeB>::type>
330  KroneckerProduct<Var, typeB>
331  Kronecker(const AbstractMatrixReadOnly<typeB>& other) const;
342  Matrix<Var> Invert(int* ierr = NULL) const;
344  Matrix<Var> CholeskyInvert(int* ierr = NULL) const;
362  template<typename typeB>
363  Matrix<typename Promote<Var, typeB>::type>
364  Solve(const AbstractMatrixReadOnly<typeB>& B, int* ierr = NULL) const;
374  template<typename typeB>
375  Matrix<typename Promote<Var, typeB>::type>
376  CholeskySolve(const AbstractMatrixReadOnly<typeB>& B, int* ierr = NULL) const;
379  Var Trace() const
380  {
381  assert(Cols() == Rows());
382  Var ret = 0.0;
383  for (enumerator i = 0; i < Rows(); ++i) ret += compute(i, i);
384  return ret;
385  }
389  void Print(INMOST_DATA_REAL_TYPE threshold = 1.0e-10, std::ostream & sout = std::cout) const
390  {
391  for (enumerator k = 0; k < Rows(); ++k)
392  {
393  for (enumerator l = 0; l < Cols(); ++l)
394  {
395  //if (__isinf__(real_part(get_value((*this)(k, l)))))
396  // sout << std::setw(12) << "inf";
397  //else if (std::isnan(get_value((*this)(k, l))))
398  // sout << std::setw(12) << "nan";
399  //else
400  if (fabs(get_value(compute(k, l))) > threshold)
401  sout << std::setw(12) << get_value(compute(k, l));
402  else
403  sout << std::setw(12) << 0;
404  sout << " ";
405  }
406  sout << std::endl;
407  }
408  }
411  bool isSymmetric(double eps = 1.0e-7) const
412  {
413  if (Rows() != Cols()) return false;
414  for (enumerator k = 0; k < Rows(); ++k)
415  {
416  for (enumerator l = k + 1; l < Rows(); ++l)
417  if (fabs(compute(k, l) - compute(l, k)) > eps)
418  return false;
419  }
420  return true;
421  }
426  template<typename typeB>
427  typename Promote<Var, typeB>::type
428  DotProduct(const AbstractMatrixReadOnly<typeB>& other) const
429  {
430  assert(Cols() == other.Cols());
431  assert(Rows() == other.Rows());
432  typename Promote<Var, typeB>::type ret = 0.0;
433  for (enumerator i = 0; i < Rows(); ++i)
434  for (enumerator j = 0; j < Cols(); ++j)
435  ret += compute(i, j) * other.compute(i, j);
436  return ret;
437  }
442  template<typename typeB>
443  typename Promote<Var, typeB>::type operator ^(const AbstractMatrixReadOnly<typeB>& other) const
444  {
445  return DotProduct(other);
446  }
449  typename SelfPromote<Var>::type FrobeniusNorm() const
450  {
451  typename SelfPromote<Var>::type ret = 0;
452  for (enumerator i = 0; i < Rows(); ++i)
453  for (enumerator j = 0; j < Cols(); ++j)
454  ret += pow(compute(i, j), 2);
455  return sqrt(ret);
456  }
459  Var MaxNorm() const
460  {
461  Var ret = 0;
462  for (enumerator i = 0; i < Rows(); ++i)
463  for (enumerator j = 0; j < Cols(); ++j)
464  ret = std::max<Var>(ret, compute(i, j));
465  return ret;
466  }
474  Matrix<Var> PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, int* ierr = NULL) const;
483  Matrix<Var> Power(INMOST_DATA_REAL_TYPE n, int * ierr = NULL) const;
491  Matrix<Var> Root(INMOST_DATA_ENUM_TYPE iter = 25, INMOST_DATA_REAL_TYPE tol = 1.0e-7, int* ierr = NULL) const;
502  template<typename typeB>
503  Matrix<typename Promote<Var, typeB>::type>
504  PseudoSolve(const AbstractMatrixReadOnly<typeB>& B, INMOST_DATA_REAL_TYPE tol = 0, int* ierr = NULL) const;
514  Matrix<Var> ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const;
519  ConstMatrixRepack<Var> Repack(enumerator rows, enumerator cols) const {return ConstMatrixRepack<Var>(*this, rows, cols);}
529  ConstSubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const;
539  ConstBlockOfMatrix<Var> BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col) const;
543  template<typename typeB>
544  //Matrix<typename Promote<Var, typeB>::type>
545  MatrixMulCoef<Var, typeB, typename Promote<Var, typeB>::type>
546  operator*(const typeB& coef) const;
547 #if defined(USE_AUTODIFF)
551  template<class A>
552  //Matrix<typename Promote<Var, variable>::type>
553  MatrixMulShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>
554  operator*(shell_expression<A> const& coef) const;// { return operator*(variable(coef)); }
555 #endif //USE_AUTODIFF
559  template<typename typeB>
560  //Matrix<typename Promote<Var, typeB>::type>
561  MatrixDivCoef<Var, typeB, typename Promote<Var, typeB>::type>
562  operator/(const typeB & coef) const;
563 #if defined(USE_AUTODIFF)
567  template<class A>
568  //Matrix<typename Promote<Var, variable>::type>
569  MatrixDivShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>
570  operator/(shell_expression<A> const& coef) const;// { return operator/(variable(coef)); }
571 #endif //USE_AUTODIFF
578  template<typename typeB>
579  Matrix<typename Promote<Var, typeB>::type>
580  operator/(const AbstractMatrixReadOnly<typeB>& other) const
581  {
582  Matrix<typename Promote<Var, typeB>::type> ret(other.Cols(), Cols());
583  ret = other.Solve(*this);
584  return ret;
585  }
592  ConstMatrixConcatCols<Var> ConcatCols(const AbstractMatrixReadOnly<Var>& B) const
593  {
594  return ConstMatrixConcatCols<Var>(*this, B);
595  }
602  template<typename VarB>
603  ConstMatrixConcatCols2<Var, VarB, typename Promote<Var, VarB>::type>
604  ConcatCols(const AbstractMatrixReadOnly<VarB>& B) const
605  {
606  return ConstMatrixConcatCols2<Var, VarB, typename Promote<Var, VarB>::type>(*this, B);
607  }
614  ConstMatrixConcatRows<Var> ConcatRows(const AbstractMatrixReadOnly<Var>& B) const
615  {
616  return ConstMatrixConcatRows<Var>(*this, B);
617  }
624  template<typename VarB>
625  ConstMatrixConcatRows2<Var,VarB,typename Promote<Var,VarB>::type>
626  ConcatRows(const AbstractMatrixReadOnly<VarB>& B) const
627  {
628  return ConstMatrixConcatRows2<Var, VarB, typename Promote<Var, VarB>::type>(*this, B);
629  }
631  virtual ~AbstractMatrixReadOnly() {}
633  __INLINE virtual INMOST_DATA_ENUM_TYPE GetMatrixCount() const
634  {
635  INMOST_DATA_ENUM_TYPE cnt = 0;
636  for (INMOST_DATA_ENUM_TYPE i = 0; i < Rows(); ++i)
637  for (INMOST_DATA_ENUM_TYPE j = 0; j < Cols(); ++j)
638  cnt += GetCount(compute(i, j));
639  return cnt;
640  }
641  };
645  template<typename Var>
646  class AbstractMatrix : public AbstractMatrixReadOnly<Var>
647  {
648  public:
649  using AbstractMatrixReadOnly<Var>::compute;
650  using AbstractMatrixReadOnly<Var>::Rows;
651  using AbstractMatrixReadOnly<Var>::Cols;
652  using AbstractMatrixReadOnly<Var>::Transpose;
653  using AbstractMatrixReadOnly<Var>::BlockOf;
654  using AbstractMatrixReadOnly<Var>::Repack;
655  using AbstractMatrixReadOnly<Var>::ConcatRows;
656  using AbstractMatrixReadOnly<Var>::ConcatCols;
657 
658  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator;
660  AbstractMatrix() {}
661  /*
664  AbstractMatrix(const AbstractMatrix & b)
665  {
666  Resize(b.Rows(),b.Cols());
667  for(enumerator i = 0; i < b.Rows(); ++i)
668  for(enumerator j = 0; j < b.Cols(); ++j)
669  (*this)(i,j) = b(i,j);
670  }
673  template<typename typeB>
674  AbstractMatrix(const AbstractMatrix<typeB> & b)
675  {
676  Resize(b.Rows(),b.Cols());
677  for(enumerator i = 0; i < b.Rows(); ++i)
678  for(enumerator j = 0; j < b.Cols(); ++j)
679  assign((*this)(i,j),b(i,j));
680  }
681  */
685  AbstractMatrix & operator =(AbstractMatrix const & other)
686  {
687  Resize(other.Rows(),other.Cols());
688  for(enumerator i = 0; i < other.Rows(); ++i)
689  for(enumerator j = 0; j < other.Cols(); ++j)
690  assign((*this)(i,j),other.get(i,j));
691  return *this;
692  }
696  template<typename typeB>
697  AbstractMatrix & operator =(AbstractMatrixReadOnly<typeB> const & other)
698  {
699  Resize(other.Rows(),other.Cols());
700  for(enumerator i = 0; i < other.Rows(); ++i)
701  for(enumerator j = 0; j < other.Cols(); ++j)
702  assign((*this)(i,j),other.compute(i,j));
703  return *this;
704  }
708  AbstractMatrix & operator =(Var const & b)
709  {
710  for(enumerator i = 0; i < Rows(); ++i)
711  for(enumerator j = 0; j < Cols(); ++j)
712  assign((*this)(i,j),b);
713  return *this;
714  }
719  virtual Var & operator () (enumerator i, enumerator j) = 0;
724  virtual const Var & get(enumerator i, enumerator j) const = 0;
728  virtual void Resize(enumerator rows, enumerator cols) = 0;
730  void Zero()
731  {
732  for(enumerator i = 0; i < Rows(); ++i)
733  for(enumerator j = 0; j < Cols(); ++j)
734  (*this)(i,j) = 0.0;
735  }
737  virtual void Swap(AbstractMatrix<Var> & b);
741  template<typename typeB>
742  AbstractMatrix & operator-=(const AbstractMatrixReadOnly<typeB> & other);
746  template<typename typeB>
747  AbstractMatrix& operator-=(const AbstractMatrix<typeB>& other);
751  template<typename typeB>
752  AbstractMatrix & operator+=(const AbstractMatrixReadOnly<typeB> & other);
756  template<typename typeB>
757  AbstractMatrix& operator+=(const AbstractMatrix<typeB>& other);
761  template<typename typeB>
762  AbstractMatrix & operator*=(const AbstractMatrixReadOnly<typeB> & B);
766  template<typename typeB>
767  AbstractMatrix & operator*=(typeB coef);
771  template<typename typeB>
772  AbstractMatrix & operator/=(typeB coef);
782  SubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
792  BlockOfMatrix<Var> BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col);
795  //Matrix<Var> Transpose() const;
796  MatrixTranspose<Var> Transpose() { return MatrixTranspose<Var>(*this); }
801  MatrixRepack<Var> Repack(enumerator rows, enumerator cols) { return MatrixRepack<Var>(*this, rows, cols); }
808  MatrixConcatCols<Var> ConcatCols(AbstractMatrix<Var>& B) { return MatrixConcatCols<Var>(*this, B); }
815  MatrixConcatRows<Var> ConcatRows(AbstractMatrix<Var>& B) { return MatrixConcatRows<Var>(*this, B); }
818  bool CheckNans() const
819  {
820  for (enumerator i = 0; i < Rows(); ++i)
821  for (enumerator j = 0; j < Cols(); ++j)
822  if (check_nans(compute(i, j))) return true;
823  return false;
824  }
827  bool CheckInfs() const
828  {
829  for (enumerator i = 0; i < Rows(); ++i)
830  for (enumerator j = 0; j < Cols(); ++j)
831  if (check_infs(compute(i, j))) return true;
832  return false;
833  }
836  bool CheckNansInfs() const
837  {
838  for (enumerator i = 0; i < Rows(); ++i)
839  for (enumerator j = 0; j < Cols(); ++j)
840  if (check_nans_infs(compute(i, j))) return true;
841  return false;
842  }
847  template<typename typeB>
848  Matrix<typename Promote<Var, typeB>::type>
849  Transform(const AbstractMatrix<typeB>& other) const;
854  template<typename typeB>
855  Matrix<typename Promote<Var, typeB>::type>
856  Transform(const AbstractMatrixReadOnly<typeB>& other) const;
862  template<typename typeB>
863  Matrix<typename Promote<Var, typeB>::type>
864  CrossProduct(const AbstractMatrix<typeB>& other) const;
870  template<typename typeB>
871  Matrix<typename Promote<Var, typeB>::type>
872  CrossProduct(const AbstractMatrixReadOnly<typeB>& other) const;
877  template<typename typeB>
878  typename Promote<Var, typeB>::type
879  DotProduct(const AbstractMatrix<typeB>& other) const
880  {
881  assert(Cols() == other.Cols());
882  assert(Rows() == other.Rows());
883  typename Promote<Var, typeB>::type ret = 0.0;
884  for (enumerator i = 0; i < Rows(); ++i)
885  for (enumerator j = 0; j < Cols(); ++j)
886  ret += get(i, j) * other.get(i, j);
887  return ret;
888  }
893  template<typename typeB>
894  typename Promote<Var, typeB>::type
895  DotProduct(const AbstractMatrixReadOnly<typeB>& other) const
896  {
897  assert(Cols() == other.Cols());
898  assert(Rows() == other.Rows());
899  typename Promote<Var, typeB>::type ret = 0.0;
900  for (enumerator i = 0; i < Rows(); ++i)
901  for (enumerator j = 0; j < Cols(); ++j)
902  ret += get(i, j) * other.compute(i, j);
903  return ret;
904  }
909  template<typename typeB>
910  typename Promote<Var, typeB>::type operator ^(const AbstractMatrix<typeB>& other) const
911  {
912  return DotProduct(other);
913  }
918  template<typename typeB>
919  typename Promote<Var, typeB>::type operator ^(const AbstractMatrixReadOnly<typeB>& other) const
920  {
921  return DotProduct(other);
922  }
925  typename SelfPromote<Var>::type FrobeniusNorm() const
926  {
927  typename SelfPromote<Var>::type ret = 0;
928  for (enumerator i = 0; i < Rows(); ++i)
929  for (enumerator j = 0; j < Cols(); ++j)
930  ret += pow(get(i, j), 2);
931  return sqrt(ret);
932  }
935  Var MaxNorm() const
936  {
937  Var ret = 0;
938  for (enumerator i = 0; i < Rows(); ++i)
939  for (enumerator j = 0; j < Cols(); ++j)
940  ret = std::max<Var>(ret, get(i, j));
941  return ret;
942  }
945  Var Trace() const
946  {
947  assert(Cols() == Rows());
948  Var ret = 0.0;
949  for (enumerator i = 0; i < Rows(); ++i) ret += get(i, i);
950  return ret;
951  }
955  void Print(INMOST_DATA_REAL_TYPE threshold = 1.0e-10, std::ostream& sout = std::cout) const
956  {
957  for (enumerator k = 0; k < Rows(); ++k)
958  {
959  for (enumerator l = 0; l < Cols(); ++l)
960  {
961  //if (__isinf__(real_part(get_value((*this)(k, l)))))
962  // sout << std::setw(12) << "inf";
963  //else if (std::isnan(get_value((*this)(k, l))))
964  // sout << std::setw(12) << "nan";
965  //else
966  if (fabs(get_value(compute(k, l))) > threshold)
967  sout << std::setw(12) << get_value(get(k, l));
968  else
969  sout << std::setw(12) << 0;
970  sout << " ";
971  }
972  sout << std::endl;
973  }
974  }
977  bool isSymmetric(double eps = 1.0e-7) const
978  {
979  if (Rows() != Cols()) return false;
980  for (enumerator k = 0; k < Rows(); ++k)
981  {
982  for (enumerator l = k + 1; l < Rows(); ++l)
983  if (fabs(get(k, l) - get(l, k)) > eps)
984  return false;
985  }
986  return true;
987  }
997  void MPT(INMOST_DATA_ENUM_TYPE* Perm, INMOST_DATA_REAL_TYPE* SL = NULL, INMOST_DATA_REAL_TYPE* SR = NULL) const;
999  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
1000  {
1001  INMOST_DATA_ENUM_TYPE cnt = 0;
1002  for (INMOST_DATA_ENUM_TYPE i = 0; i < Rows(); ++i)
1003  for (INMOST_DATA_ENUM_TYPE j = 0; j < Cols(); ++j)
1004  cnt += GetCount(get(i, j));
1005  return cnt;
1006  }
1008  virtual ~AbstractMatrix() {};
1009  };
1010 
1011 
1012  template<typename Var, typename storage_type>
1013  class SymmetricMatrix : public AbstractMatrix<Var>
1014  {
1015  public:
1016  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
1017  using AbstractMatrix<Var>::operator();
1018  using AbstractMatrix<Var>::get;
1019  protected:
1020  storage_type space; //< Array of row-wise stored elements.
1021  enumerator n; //< Number of rows.
1022  public:
1024  void Swap(AbstractMatrix<Var> & b)
1025  {
1026 #if defined(_CPPRTTI) || defined(__GXX_RTTI)
1027  SymmetricMatrix<Var,storage_type> * bb = dynamic_cast<SymmetricMatrix<Var,storage_type> *>(&b);
1028  if( bb != NULL )
1029  {
1030  space.swap((*bb).space);
1031  std::swap(n,(*bb).n);
1032  }
1033  else AbstractMatrix<Var>::Swap(b);
1034 #else //_CPPRTTI
1035  AbstractMatrix<Var>::Swap(b);
1036 #endif //_CPPRTTI
1037  }
1039  SymmetricMatrix() : space(), n(0) {}
1044  SymmetricMatrix(const Var * pspace, enumerator pn) : space(pspace,pspace+pn*(pn+1)/2), n(pn) {}
1052  SymmetricMatrix(const storage_type & pspace, enumerator pn) : space(pspace), n(pn) {}
1058  SymmetricMatrix(const storage_type & pspace) : space(pspace), n(0) {}
1063  SymmetricMatrix(enumerator pn) : space(pn*(pn+1)/2), n(pn) {}
1068  static SymmetricMatrix Make(enumerator pn, ...)
1069  {
1070  SymmetricMatrix A(pn);
1071  va_list argptr;
1072  va_start(argptr, pn);
1073  for (enumerator j = 0; j < pn; ++j)
1074  for (enumerator i = j; i < pn; ++i)
1075  {
1076  Var val = va_arg(argptr, Var);
1077  A(i, j) = val;
1078  }
1079  va_end(argptr);
1080  return A;
1081  }
1086  SymmetricMatrix(enumerator pn, const Var & c) : space(pn*(pn+1)/2,c), n(pn) {}
1089  SymmetricMatrix(const SymmetricMatrix & other) : space(other.space), n(other.n)
1090  {
1091  //for(enumerator i = 0; i < n*m; ++i)
1092  // space[i] = other.space[i];
1093  }
1100  template<typename typeB>
1101  SymmetricMatrix(const AbstractMatrixReadOnly<typeB> & other) : space(other.Rows()*(other.Rows()+1)/2), n(other.Rows())
1102  {
1103  assert(other.Rows() == other.Cols());
1104  for(enumerator i = 0; i < n; ++i)
1105  for(enumerator j = i; j < n; ++j)
1106  assign((*this)(i,j),other.compute(i,j));
1107  }
1114  template<typename typeB>
1115  SymmetricMatrix(const AbstractMatrix<typeB>& other) : space(other.Rows()* (other.Rows() + 1) / 2), n(other.Rows())
1116  {
1117  assert(other.Rows() == other.Cols());
1118  for (enumerator i = 0; i < n; ++i)
1119  for (enumerator j = i; j < n; ++j)
1120  assign((*this)(i, j), other.get(i, j));
1121  }
1123  ~SymmetricMatrix() {}
1128  void Resize(enumerator nrows, enumerator ncols)
1129  {
1130  (void)ncols;
1131  assert(nrows == ncols);
1132  if( space.size() != (nrows+1)*nrows/2 )
1133  space.resize((nrows+1)*nrows/2);
1134  n = nrows;
1135  }
1139  SymmetricMatrix & operator =(SymmetricMatrix const & other)
1140  {
1141  if( this != &other )
1142  {
1143  if( n != other.n )
1144  space.resize(other.n*(other.n+1)/2);
1145  for(enumerator i = 0; i < other.n*(other.n+1)/2; ++i)
1146  space[i] = other.space[i];
1147  n = other.n;
1148  }
1149  return *this;
1150  }
1156  template<typename typeB>
1157  SymmetricMatrix& operator =(AbstractMatrix<typeB> const& other)
1158  {
1159  assert(other.Rows() == other.Cols());
1160  if (Rows() != other.Rows())
1161  space.resize((other.Rows() + 1) * other.Rows() / 2);
1162  n = other.Rows();
1163  for (enumerator i = 0; i < other.Rows(); ++i)
1164  for (enumerator j = i; j < other.Cols(); ++j)
1165  assign((*this)(i, j), other.get(i, j));
1166  return *this;
1167  }
1173  template<typename typeB>
1174  SymmetricMatrix & operator =(AbstractMatrixReadOnly<typeB> const & other)
1175  {
1176  assert(other.Rows() == other.Cols());
1177  if (Rows() != other.Rows())
1178  space.resize((other.Rows() + 1) * other.Rows() / 2);
1179  n = other.Rows();
1180  for (enumerator i = 0; i < other.Rows(); ++i)
1181  for (enumerator j = i; j < other.Cols(); ++j)
1182  assign((*this)(i, j), other.compute(i, j));
1183  return *this;
1184  }
1189  __INLINE Var & operator()(enumerator i, enumerator j)
1190  {
1191  assert(i < n);
1192  assert(j < n);
1193  if( i > j ) std::swap(i,j);
1194  return space[j+n*i-i*(i+1)/2];
1195  }
1200  __INLINE const Var& get(enumerator i, enumerator j) const
1201  {
1202  assert(i < n);
1203  assert(j < n);
1204  if (i > j) std::swap(i, j);
1205  return space[j + n * i - i * (i + 1) / 2];
1206  }
1212  __INLINE Var compute(enumerator i, enumerator j) const
1213  {
1214  assert(i < n);
1215  assert(j < n);
1216  if( i > j ) std::swap(i,j);
1217  return space[j+n*i-i*(i+1)/2];
1218  }
1219 
1222  __INLINE Var * data() {return space.data();}
1226  __INLINE const Var * data() const {return space.data();}
1229  __INLINE enumerator Rows() const {return n;}
1232  __INLINE enumerator Cols() const {return n;}
1235  __INLINE enumerator & Rows() {return n;}
1238  __INLINE enumerator & Cols() {return n;}
1257  static SymmetricMatrix<Var> FromTensor(const Var * K, enumerator size, enumerator matsize = 3)
1258  {
1259  SymmetricMatrix<Var> Kc(matsize,matsize);
1260  if( matsize == 1 )
1261  {
1262  assert(size == 1);
1263  Kc(0,0) = K[0];
1264  }
1265  if( matsize == 2 )
1266  {
1267  assert(size == 1 || size == 2 || size == 3 || size == 4);
1268  switch(size)
1269  {
1270  case 1: //scalar
1271  Kc(0,0) = Kc(1,1) = K[0];
1272  break;
1273  case 2: //diagonal
1274  Kc(0,0) = K[0]; // KXX
1275  Kc(1,1) = K[1]; // KYY
1276  break;
1277  case 3: //symmetric
1278  Kc(0,0) = K[0]; // KXX
1279  Kc(0,1) = K[1]; //KXY
1280  Kc(1,1) = K[2]; //KYY
1281  break;
1282  }
1283  }
1284  else if( matsize == 3 )
1285  {
1286  assert(size == 1 || size == 3 || size == 6 || size == 9);
1287  switch(size)
1288  {
1289  case 1: //scalar permeability tensor
1290  Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0];
1291  break;
1292  case 3: //diagonal permeability tensor
1293  Kc(0,0) = K[0]; //KXX
1294  Kc(1,1) = K[1]; //KYY
1295  Kc(2,2) = K[2]; //KZZ
1296  break;
1297  case 6: //symmetric permeability tensor
1298  Kc(0,0) = K[0]; //KXX
1299  Kc(0,1) = K[1]; //KXY
1300  Kc(0,2) = K[2]; //KXZ
1301  Kc(1,1) = K[3]; //KYY
1302  Kc(1,2) = K[4]; //KYZ
1303  Kc(2,2) = K[5]; //KZZ
1304  break;
1305  }
1306  }
1307  else if( matsize == 6 )
1308  {
1309  assert(size == 1 || size == 6 || size == 21 || size == 36);
1310  switch(size)
1311  {
1312  case 1: //scalar elasticity tensor
1313  Kc(0,0) = Kc(1,1) = Kc(2,2) = Kc(3,3) = Kc(4,4) = Kc(5,5) = K[0];
1314  break;
1315  case 6: //diagonal elasticity tensor
1316  Kc(0,0) = K[0]; //KXX
1317  Kc(1,1) = K[1]; //KYY
1318  Kc(2,2) = K[2]; //KZZ
1319  break;
1320  case 21: //symmetric elasticity tensor (note - diagonal first, then off-diagonal rows)
1321  {
1322  Kc(0,0) = K[0]; //c11
1323  Kc(0,1) = K[1]; //c12
1324  Kc(0,2) = K[2]; //c13
1325  Kc(0,3) = K[3]; //c14
1326  Kc(0,4) = K[4]; //c15
1327  Kc(0,5) = K[5]; //c16
1328  Kc(1,1) = K[6]; //c22
1329  Kc(1,2) = K[7]; //c23
1330  Kc(1,3) = K[8]; //c24
1331  Kc(1,4) = K[9]; //c25
1332  Kc(1,5) = K[10]; //c26
1333  Kc(2,2) = K[11]; //c33
1334  Kc(2,3) = K[12]; //c34
1335  Kc(2,4) = K[13]; //c35
1336  Kc(2,5) = K[14]; //c36
1337  Kc(3,3) = K[15]; //c44
1338  Kc(3,4) = K[16]; //c45
1339  Kc(3,5) = K[17]; //c46
1340  Kc(4,4) = K[18]; //c55
1341  Kc(4,5) = K[19]; //c56
1342  Kc(5,5) = K[20]; //c66
1343  break;
1344  }
1345  }
1346  }
1347  return Kc;
1348  }
1353  static SymmetricMatrix<Var> FromDiagonal(const Var * r, enumerator size)
1354  {
1355  SymmetricMatrix<Var> ret(size);
1356  ret.Zero();
1357  for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k];
1358  return ret;
1359  }
1364  static SymmetricMatrix<Var> FromDiagonalInverse(const Var * r, enumerator size)
1365  {
1366  SymmetricMatrix<Var> ret(size);
1367  ret.Zero();
1368  for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
1369  return ret;
1370  }
1376  static SymmetricMatrix<Var> Unit(enumerator pn, const Var & c = 1.0)
1377  {
1378  SymmetricMatrix<Var> ret(pn,0.0);
1379  for(enumerator i = 0; i < pn; ++i) ret(i,i) = c;
1380  return ret;
1381  }
1382 
1392  //::INMOST::SubMatrix<Var,storage_type> SubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
1393 
1403  //::INMOST::SubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
1404 
1405  //::INMOST::ConstSubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const;
1406  };
1407 
1432  template<typename Var, typename storage_type>
1433  class Matrix : public AbstractMatrix<Var>
1434  {
1435  public:
1436  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
1437  using AbstractMatrix<Var>::operator();
1438  using AbstractMatrix<Var>::get;
1439  protected:
1440  storage_type space; //< Array of row-wise stored elements.
1441  enumerator n; //< Number of rows.
1442  enumerator m; //< Number of columns.
1443  public:
1444  Var& operator [](enumerator k) { return space[k]; }
1445  const Var& operator[](enumerator k) const { return space[k]; }
1448  void RemoveRow(enumerator row)
1449  {
1450  assert(row < n);
1451  for(enumerator k = row+1; k < n; ++k)
1452  {
1453  for(enumerator l = 0; l < m; ++l)
1454  (*this)(k-1,l) = (*this)(k,l);
1455  }
1456  space.resize((n-1)*m);
1457  --n;
1458  }
1464  void RemoveRows(enumerator first, enumerator last)
1465  {
1466  assert(first < n);
1467  assert(last <= n);
1468  assert(first <= last);
1469  enumerator shift = last-first;
1470  for(enumerator k = last; k < n; ++k)
1471  {
1472  for(enumerator l = 0; l < m; ++l)
1473  (*this)(k-shift,l) = (*this)(k,l);
1474  }
1475  space.resize((n-shift)*m);
1476  n-=shift;
1477  }
1480  void RemoveColumn(enumerator col)
1481  {
1482  assert(col < m);
1483  Matrix<Var> tmp(n,m-1);
1484  for(enumerator k = 0; k < n; ++k)
1485  {
1486  for(enumerator l = 0; l < col; ++l)
1487  tmp(k,l) = (*this)(k,l);
1488  for(enumerator l = col+1; l < m; ++l)
1489  tmp(k,l-1) = (*this)(k,l);
1490  }
1491  this->Swap(tmp);
1492  }
1498  void RemoveColumns(enumerator first, enumerator last)
1499  {
1500  assert(first < m);
1501  assert(last <= m);
1502  assert(first <= last);
1503  enumerator shift = last-first;
1504  Matrix<Var> tmp(n,m-shift);
1505  for(enumerator k = 0; k < n; ++k)
1506  {
1507  for(enumerator l = 0; l < first; ++l)
1508  tmp(k,l) = (*this)(k,l);
1509  for(enumerator l = last; l < m; ++l)
1510  tmp(k,l-shift) = (*this)(k,l);
1511  }
1512  this->Swap(tmp);
1513  }
1519  void RemoveSubset(enumerator firstrow, enumerator lastrow, enumerator firstcol, enumerator lastcol)
1520  {
1521  enumerator shiftrow = lastrow-firstrow;
1522  enumerator shiftcol = lastcol-firstcol;
1523  Matrix<Var> tmp(n-shiftrow, m-shiftcol);
1524  for(enumerator k = 0; k < firstrow; ++k)
1525  {
1526  for(enumerator l = 0; l < firstcol; ++l)
1527  tmp(k,l) = (*this)(k,l);
1528  for(enumerator l = lastcol; l < m; ++l)
1529  tmp(k,l-shiftcol) = (*this)(k,l);
1530  }
1531  for(enumerator k = lastrow; k < n; ++k)
1532  {
1533  for(enumerator l = 0; l < firstcol; ++l)
1534  tmp(k-shiftrow,l) = (*this)(k,l);
1535  for(enumerator l = lastcol; l < m; ++l)
1536  tmp(k-shiftrow,l-shiftcol) = (*this)(k,l);
1537  }
1538  this->Swap(tmp);
1539  }
1541  void Swap(AbstractMatrix<Var> & b)
1542  {
1543 #if defined(_CPPRTTI) || defined(__GXX_RTTI)
1544  Matrix<Var,storage_type> * bb = dynamic_cast<Matrix<Var,storage_type> *>(&b);
1545  if( bb != NULL )
1546  {
1547  space.swap((*bb).space);
1548  std::swap(n,(*bb).n);
1549  std::swap(m,(*bb).m);
1550  }
1551  else AbstractMatrix<Var>::Swap(b);
1552 #else //_CPPRTTI
1553  AbstractMatrix<Var>::Swap(b);
1554 #endif //_CPPRTTI
1555  }
1557  Matrix() : space(), n(0), m(0) {}
1562  Matrix(const Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {}
1570  Matrix(const storage_type & pspace, enumerator pn, enumerator pm) : space(pspace), n(pn), m(pm) {}
1576  Matrix(const storage_type & pspace) : space(pspace), n(0), m(0) {}
1581  Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {}
1586  Matrix(enumerator pn, enumerator pm, const Var & c) : space(pn*pm,c), n(pn), m(pm) {}
1591  static Matrix Make(enumerator pn, enumerator pm, ...)
1592  {
1593  Matrix A(pn, pm);
1594  va_list argptr;
1595  va_start(argptr, pm);
1596  for (enumerator i = 0; i < pn; ++i)
1597  for (enumerator j = 0; j < pm; ++j)
1598  {
1599  Var val = va_arg(argptr, Var);
1600  A(i, j) = val;
1601  }
1602  va_end(argptr);
1603  return A;
1604  }
1607  Matrix(const Matrix & other) : space(other.space), n(other.n), m(other.m)
1608  {
1609  //for(enumerator i = 0; i < n*m; ++i)
1610  // space[i] = other.space[i];
1611  }
1616  template<typename typeB>
1617  Matrix(const AbstractMatrixReadOnly<typeB> & other) : space(other.Cols()*other.Rows()), n(other.Rows()), m(other.Cols())
1618  {
1619  for(enumerator i = 0; i < n; ++i)
1620  for(enumerator j = 0; j < m; ++j)
1621  assign((*this)(i,j),other.compute(i,j));
1622  }
1627  template<typename typeB>
1628  Matrix(const AbstractMatrix<typeB>& other) : space(other.Cols()* other.Rows()), n(other.Rows()), m(other.Cols())
1629  {
1630  for (enumerator i = 0; i < n; ++i)
1631  for (enumerator j = 0; j < m; ++j)
1632  assign((*this)(i, j), other.get(i, j));
1633  }
1635  ~Matrix() {}
1639  void Resize(enumerator nrows, enumerator mcols)
1640  {
1641  if( space.size() != mcols*nrows )
1642  space.resize(mcols*nrows);
1643  n = nrows;
1644  m = mcols;
1645  }
1649  Matrix & operator =(Matrix const & other)
1650  {
1651  if( this != &other )
1652  {
1653  if( n*m != other.n*other.m )
1654  space.resize(other.n*other.m);
1655  for(enumerator i = 0; i < other.n*other.m; ++i)
1656  space[i] = other.space[i];
1657  n = other.n;
1658  m = other.m;
1659  }
1660  return *this;
1661  }
1665  template<typename typeB>
1666  Matrix & operator =(AbstractMatrixReadOnly<typeB> const & other)
1667  {
1668  if( Cols()*Rows() != other.Cols()*other.Rows() )
1669  space.resize(other.Cols()*other.Rows());
1670  n = other.Rows();
1671  m = other.Cols();
1672  for(enumerator i = 0; i < other.Rows(); ++i)
1673  for(enumerator j = 0; j < other.Cols(); ++j)
1674  assign((*this)(i,j),other.compute(i,j));
1675  return *this;
1676  }
1680  template<typename typeB>
1681  Matrix& operator =(AbstractMatrix<typeB> const& other)
1682  {
1683  if (Cols() * Rows() != other.Cols() * other.Rows())
1684  space.resize(other.Cols() * other.Rows());
1685  n = other.Rows();
1686  m = other.Cols();
1687  for (enumerator i = 0; i < other.Rows(); ++i)
1688  for (enumerator j = 0; j < other.Cols(); ++j)
1689  assign((*this)(i, j), other.get(i, j));
1690  return *this;
1691  }
1696  __INLINE Var & operator()(enumerator i, enumerator j)
1697  {
1698  assert(i < n);
1699  assert(j < m);
1700  assert(i*m+j < n*m); //overflow check?
1701  return space[i*m+j];
1702  }
1707  __INLINE Var compute(enumerator i, enumerator j) const
1708  {
1709  assert(i < n);
1710  assert(j < m);
1711  assert(i*m+j < n*m); //overflow check?
1712  return space[i*m+j];
1713  }
1719  __INLINE const Var & get(enumerator i, enumerator j) const
1720  {
1721  assert(i < n);
1722  assert(j < m);
1723  assert(i * m + j < n* m); //overflow check?
1724  return space[i * m + j];
1725  }
1728  __INLINE Var * data() {return space.data();}
1732  __INLINE const Var * data() const {return space.data();}
1735  __INLINE enumerator Rows() const {return n;}
1738  __INLINE enumerator Cols() const {return m;}
1741  __INLINE enumerator & Rows() {return n;}
1744  __INLINE enumerator & Cols() {return m;}
1754  static Matrix<Var> Permutation(const INMOST_DATA_ENUM_TYPE * Perm, enumerator size)
1755  {
1756  Matrix<Var> Ret(size,size);
1757  Ret.Zero();
1758  for(enumerator k = 0; k < size; ++k)
1759  Ret(k,Perm[k]) = 1;
1760  return Ret;
1761  }
1783  static Matrix<Var> FromTensor(const Var * K, enumerator size, enumerator matsize = 3)
1784  {
1785  Matrix<Var> Kc(matsize,matsize);
1786  if( matsize == 1 )
1787  {
1788  assert(size == 1);
1789  Kc(0,0) = K[0];
1790  }
1791  if( matsize == 2 )
1792  {
1793  assert(size == 1 || size == 2 || size == 3 || size == 4);
1794  switch(size)
1795  {
1796  case 1: //scalar
1797  Kc(0,0) = Kc(1,1) = K[0];
1798  break;
1799  case 2: //diagonal
1800  Kc(0,0) = K[0]; // KXX
1801  Kc(1,1) = K[1]; // KYY
1802  break;
1803  case 3: //symmetric
1804  Kc(0,0) = K[0]; // KXX
1805  Kc(0,1) = Kc(1,0) = K[1]; //KXY
1806  Kc(1,1) = K[2]; //KYY
1807  break;
1808  case 4: //full
1809  Kc(0,0) = K[0]; //KXX
1810  Kc(0,1) = K[1]; //KXY
1811  Kc(1,0) = K[2]; //KYX
1812  Kc(1,1) = K[3]; //KYY
1813  break;
1814  }
1815  }
1816  else if( matsize == 3 )
1817  {
1818  assert(size == 1 || size == 3 || size == 6 || size == 9);
1819  switch(size)
1820  {
1821  case 1: //scalar permeability tensor
1822  Kc(0,0) = Kc(1,1) = Kc(2,2) = K[0];
1823  break;
1824  case 3: //diagonal permeability tensor
1825  Kc(0,0) = K[0]; //KXX
1826  Kc(1,1) = K[1]; //KYY
1827  Kc(2,2) = K[2]; //KZZ
1828  break;
1829  case 6: //symmetric permeability tensor
1830  Kc(0,0) = K[0]; //KXX
1831  Kc(0,1) = Kc(1,0) = K[1]; //KXY
1832  Kc(0,2) = Kc(2,0) = K[2]; //KXZ
1833  Kc(1,1) = K[3]; //KYY
1834  Kc(1,2) = Kc(2,1) = K[4]; //KYZ
1835  Kc(2,2) = K[5]; //KZZ
1836  break;
1837  case 9: //full permeability tensor
1838  Kc(0,0) = K[0]; //KXX
1839  Kc(0,1) = K[1]; //KXY
1840  Kc(0,2) = K[2]; //KXZ
1841  Kc(1,0) = K[3]; //KYX
1842  Kc(1,1) = K[4]; //KYY
1843  Kc(1,2) = K[5]; //KYZ
1844  Kc(2,0) = K[6]; //KZX
1845  Kc(2,1) = K[7]; //KZY
1846  Kc(2,2) = K[8]; //KZZ
1847  break;
1848  }
1849  }
1850  else if( matsize == 6 )
1851  {
1852  assert(size == 1 || size == 6 || size == 21 || size == 36);
1853  switch(size)
1854  {
1855  case 1: //scalar elasticity tensor
1856  Kc(0,0) = Kc(1,1) = Kc(2,2) = Kc(3,3) = Kc(4,4) = Kc(5,5) = K[0];
1857  break;
1858  case 6: //diagonal elasticity tensor
1859  Kc(0,0) = K[0]; //KXX
1860  Kc(1,1) = K[1]; //KYY
1861  Kc(2,2) = K[2]; //KZZ
1862  break;
1863  case 21: //symmetric elasticity tensor (note - diagonal first, then off-diagonal rows)
1864  {
1865  Kc(0,0) = K[0]; //c11
1866  Kc(0,1) = Kc(1,0) = K[1]; //c12
1867  Kc(0,2) = Kc(2,0) = K[2]; //c13
1868  Kc(0,3) = Kc(3,0) = K[3]; //c14
1869  Kc(0,4) = Kc(4,0) = K[4]; //c15
1870  Kc(0,5) = Kc(5,0) = K[5]; //c16
1871  Kc(1,1) = K[6]; //c22
1872  Kc(1,2) = Kc(2,1) = K[7]; //c23
1873  Kc(1,3) = Kc(3,1) = K[8]; //c24
1874  Kc(1,4) = Kc(4,1) = K[9]; //c25
1875  Kc(1,5) = Kc(5,1) = K[10]; //c26
1876  Kc(2,2) = K[11]; //c33
1877  Kc(2,3) = Kc(3,2) = K[12]; //c34
1878  Kc(2,4) = Kc(4,2) = K[13]; //c35
1879  Kc(2,5) = Kc(5,2) = K[14]; //c36
1880  Kc(3,3) = K[15]; //c44
1881  Kc(3,4) = Kc(4,3) = K[16]; //c45
1882  Kc(3,5) = Kc(5,3) = K[17]; //c46
1883  Kc(4,4) = K[18]; //c55
1884  Kc(4,5) = Kc(5,4) = K[19]; //c56
1885  Kc(5,5) = K[20]; //c66
1886  break;
1887  }
1888  case 36: //full elasticity tensor
1889  for(int i = 0; i < 6; ++i)
1890  for(int j = 0; j < 6; ++j)
1891  Kc(i,j) = K[6*i+j];
1892  break;
1893  }
1894  }
1895  return Kc;
1896  }
1901  static Matrix<Var> FromVector(const Var * r, enumerator size)
1902  {
1903  return Matrix<Var>(r,size,1);
1904  }
1909  static Matrix<Var> FromDiagonal(const Var * r, enumerator size)
1910  {
1911  Matrix<Var> ret(size,size);
1912  ret.Zero();
1913  for(enumerator k = 0; k < size; ++k) ret(k,k) = r[k];
1914  return ret;
1915  }
1920  static Matrix<Var> FromDiagonalInverse(const Var * r, enumerator size)
1921  {
1922  Matrix<Var> ret(size,size);
1923  ret.Zero();
1924  for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
1925  return ret;
1926  }
1933  static Matrix<Var> CrossProductMatrix(const Var vec[3])
1934  {
1935  // | 0 -z y |
1936  // | z 0 -x |
1937  // | -y x 0 |
1938  Matrix<Var> ret(3,3);
1939  ret(0,0) = 0.0;
1940  ret(0,1) = -vec[2]; //-z
1941  ret(0,2) = vec[1]; //y
1942  ret(1,0) = vec[2]; //z
1943  ret(1,1) = 0;
1944  ret(1,2) = -vec[0]; //-x
1945  ret(2,0) = -vec[1]; //-y
1946  ret(2,1) = vec[0]; //x
1947  ret(2,2) = 0;
1948  return ret;
1949  }
1955  static MatrixUnit<Var> Unit(enumerator pn, const Var & c = 1.0)
1956  {
1957  return MatrixUnit<Var>(pn, c);
1958  /*
1959  Matrix<Var> ret(pn,pn);
1960  ret.Zero();
1961  for(enumerator i = 0; i < pn; ++i) ret(i,i) = c;
1962  return ret;
1963  */
1964  }
1970  static MatrixRow<Var> Row(enumerator pn, const Var & c = 1.0)
1971  { return MatrixRow<Var>(pn,c); }
1977  static MatrixCol<Var> Col(enumerator pn, const Var & c = 1.0)
1978  { return MatrixCol<Var>(pn, c); }
1990  Matrix<Var> JointDiagonalization(INMOST_DATA_REAL_TYPE threshold = 1.0e-7)
1991  {
1992  enumerator N = Rows();
1993  enumerator M = Cols() / Rows();
1994  Matrix<Var> V((MatrixUnit<Var>(m)));
1995  Matrix<Var> R(2,M);
1996  Matrix<Var> G(2,2);
1997  Matrix & A = *this;
1998  Var ton, toff, theta, c, s, Ap, Aq, Vp, Vq;
1999  bool repeat;
2000  do
2001  {
2002  repeat = false;
2003  for(enumerator p = 0; p < N-1; ++p)
2004  {
2005  for(enumerator q = p+1; q < N; ++q)
2006  {
2007  for(enumerator k = 0; k < M; ++k)
2008  {
2009  R(0,k) = A(p,p + k*N) - A(q,q + k*N);
2010  R(1,k) = A(p,q + k*N) + A(q,p + k*N);
2011  }
2012  G = R*R.Transpose();
2013  Var ton = G(0,0) - G(1,1);
2014  Var toff = G(0,1) + G(1,0);
2015  Var theta = 0.5 * atan2( toff, ton + sqrt(ton*ton + toff*toff) );
2016  Var c = cos(theta);
2017  Var s = sin(theta);
2018  if( fabs(s) > threshold )
2019  {
2020  //std::cout << "p,q: " << p << "," << q << " c,s: " << c << "," << s << std::endl;
2021  repeat = true;
2022  for(enumerator k = 0; k < M; ++k)
2023  {
2024  for(enumerator i = 0; i < N; ++i)
2025  {
2026  Ap = A(i,p + k*N);
2027  Aq = A(i,q + k*N);
2028  A(i,p + k*N) = Ap*c + Aq*s;
2029  A(i,q + k*N) = Aq*c - Ap*s;
2030  }
2031  }
2032  for(enumerator k = 0; k < M; ++k)
2033  {
2034  for(enumerator j = 0; j < N; ++j)
2035  {
2036  Ap = A(p,j + k*N);
2037  Aq = A(q,j + k*N);
2038  A(p,j + k*N) = Ap*c + Aq*s;
2039  A(q,j + k*N) = Aq*c - Ap*s;
2040  }
2041  }
2042  for(enumerator i = 0; i < N; ++i)
2043  {
2044  Vp = V(i,p);
2045  Vq = V(i,q);
2046  V(i,p) = Vp*c + Vq*s;
2047  V(i,q) = Vq*c - Vp*s;
2048  }
2049  }
2050  }
2051  }
2052  //Print();
2053  } while( repeat );
2054  return V;
2055  }
2065  //::INMOST::SubMatrix<Var,storage_type> SubMatrix(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
2066 
2076  //::INMOST::SubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col);
2077 
2078  //::INMOST::ConstSubMatrix<Var> operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const;
2079  };
2081  template<typename Var>
2082  class SubMatrix : public AbstractMatrix<Var>
2083  {
2084  public:
2085  using AbstractMatrix<Var>::operator();
2086  using AbstractMatrix<Var>::operator =;
2087  using AbstractMatrix<Var>::get;
2088  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
2089  private:
2090  AbstractMatrix<Var> * M;
2091  enumerator brow; //< First row in matrix M.
2092  enumerator erow; //< Last row in matrix M.
2093  enumerator bcol; //< First column in matrix M.
2094  enumerator ecol; //< Last column in matrix M.
2095  public:
2098  __INLINE enumerator Rows() const {return erow-brow;}
2101  __INLINE enumerator Cols() const {return ecol-bcol;}
2108  SubMatrix(AbstractMatrix<Var> & rM, enumerator first_row, enumerator last_row, enumerator first_column, enumerator last_column) : M(&rM), brow(first_row), erow(last_row), bcol(first_column), ecol(last_column)
2109  {}
2110  SubMatrix(const SubMatrix & b) : M(b.M), brow(b.brow), erow(b.erow), bcol(b.bcol), ecol(b.ecol) {}
2114  template<typename typeB>
2115  SubMatrix& operator =(AbstractMatrixReadOnly<typeB> const& other)
2116  {
2117  assert(Cols() == other.Cols());
2118  assert(Rows() == other.Rows());
2119  for (enumerator i = 0; i < other.Rows(); ++i)
2120  for (enumerator j = 0; j < other.Cols(); ++j)
2121  assign((*this)(i, j), other.compute(i, j));
2122  return *this;
2123  }
2127  template<typename typeB>
2128  SubMatrix & operator =(AbstractMatrix<typeB> const & other)
2129  {
2130  assert( Cols() == other.Cols() );
2131  assert( Rows() == other.Rows() );
2132  for(enumerator i = 0; i < other.Rows(); ++i)
2133  for(enumerator j = 0; j < other.Cols(); ++j)
2134  assign((*this)(i,j),other.get(i,j));
2135  return *this;
2136  }
2140  template<typename typeB>
2141  SubMatrix & operator =(SubMatrix<typeB> const & other)
2142  {
2143  assert( Cols() == other.Cols() );
2144  assert( Rows() == other.Rows() );
2145  for(enumerator i = 0; i < other.Rows(); ++i)
2146  for(enumerator j = 0; j < other.Cols(); ++j)
2147  assign((*this)(i,j),other.get(i,j));
2148  return *this;
2149  }
2154  __INLINE Var & operator()(enumerator i, enumerator j)
2155  {
2156  assert(i < Rows());
2157  assert(j < Cols());
2158  assert(i*Cols()+j < Rows()*Cols()); //overflow check?
2159  return (*M)(i+brow,j+bcol);
2160  }
2166  __INLINE Var compute(enumerator i, enumerator j) const
2167  {
2168  assert(i < Rows());
2169  assert(j < Cols());
2170  assert(i*Cols()+j < Rows()*Cols()); //overflow check?
2171  return (*M)(i+brow,j+bcol);
2172  }
2178  __INLINE const Var & get(enumerator i, enumerator j) const
2179  {
2180  assert(i < Rows());
2181  assert(j < Cols());
2182  assert(i * Cols() + j < Rows()* Cols()); //overflow check?
2183  return M->get(i + brow, j + bcol);
2184  }
2190  ::INMOST::Matrix<Var> MakeMatrix()
2191  {
2192  ::INMOST::Matrix<Var> ret(Rows(),Cols());
2193  for(enumerator i = 0; i < Rows(); ++i)
2194  for(enumerator j = 0; j < Cols(); ++j)
2195  ret(i,j) = (*this)(i,j);
2196  return ret;
2197  }
2201  void Resize(enumerator rows, enumerator cols)
2202  {
2203  assert(Cols() == cols);
2204  assert(Rows() == rows);
2205  (void)cols; (void)rows;
2206  }
2207  };
2208 
2210  template<typename Var>
2211  class ConstSubMatrix : public AbstractMatrixReadOnly<Var>
2212  {
2213  public:
2214  using AbstractMatrixReadOnly<Var>::compute;
2215  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2216  private:
2217  const AbstractMatrixReadOnly<Var> * M;
2218  enumerator brow; //< First row in matrix M.
2219  enumerator erow; //< Last row in matrix M.
2220  enumerator bcol; //< First column in matrix M.
2221  enumerator ecol; //< Last column in matrix M.
2222  public:
2225  __INLINE enumerator Rows() const {return erow-brow;}
2228  __INLINE enumerator Cols() const {return ecol-bcol;}
2235  ConstSubMatrix(const AbstractMatrixReadOnly<Var> & rM, enumerator first_row, enumerator last_row, enumerator first_column, enumerator last_column) : M(&rM), brow(first_row), erow(last_row), bcol(first_column), ecol(last_column)
2236  {}
2237  ConstSubMatrix(const ConstSubMatrix & b) : M(b.M), brow(b.brow), erow(b.erow), bcol(b.bcol), ecol(b.ecol) {}
2238 
2244  __INLINE Var compute(enumerator i, enumerator j) const
2245  {
2246  assert(i < Rows());
2247  assert(j < Cols());
2248  assert(i*Cols()+j < Rows()*Cols()); //overflow check?
2249  return M->compute(i+brow,j+bcol);
2250  }
2256  ::INMOST::Matrix<Var> MakeMatrix()
2257  {
2258  ::INMOST::Matrix<Var> ret(Rows(),Cols());
2259  for(enumerator i = 0; i < Rows(); ++i)
2260  for(enumerator j = 0; j < Cols(); ++j)
2261  ret(i,j) = (*this)(i,j);
2262  return ret;
2263  }
2264  };
2265 
2267  template<typename Var>
2268  class BlockOfMatrix : public AbstractMatrix<Var>
2269  {
2270  public:
2271  using AbstractMatrix<Var>::operator();
2272  using AbstractMatrix<Var>::operator =;
2273  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
2274  private:
2275  AbstractMatrix<Var> * M;
2276  enumerator nrows; //< Number of rows in larger matrix.
2277  enumerator ncols; //< Number of colums in larger matrix.
2278  enumerator orow; //< Row offset in larger matrix.
2279  enumerator ocol; //< Column offset in larger matrix.
2280  public:
2283  __INLINE enumerator Rows() const {return nrows;}
2286  __INLINE enumerator Cols() const {return ncols;}
2293  BlockOfMatrix(AbstractMatrix<Var> & rM, enumerator num_rows, enumerator num_cols, enumerator offset_row, enumerator offset_col) : M(&rM), nrows(num_rows), ncols(num_cols), orow(offset_row), ocol(offset_col)
2294  {}
2295  BlockOfMatrix(const BlockOfMatrix & b) : M(b.M), nrows(b.nrows), ncols(b.ncols), orow(b.orow), ocol(b.ocol) {}
2300  template<typename typeB>
2301  BlockOfMatrix & operator =(AbstractMatrix<typeB> const & other)
2302  {
2303  assert( Cols() == other.Cols() );
2304  assert( Rows() == other.Rows() );
2305  for(enumerator i = orow; i < orow+M->Rows(); ++i)
2306  for(enumerator j = ocol; j < ocol+M->Cols(); ++j)
2307  assign((*M)(i-orow,j-ocol),other(i,j));
2308  return *this;
2309  }
2313  template<typename typeB>
2314  BlockOfMatrix & operator =(BlockOfMatrix<typeB> const & other)
2315  {
2316  assert( Cols() == other.Cols() );
2317  assert( Rows() == other.Rows() );
2318  for(enumerator i = orow; i < orow+M->Rows(); ++i)
2319  for(enumerator j = ocol; j < ocol+M->Cols(); ++j)
2320  assign((*M)(i-orow,j-ocol),other(i,j));
2321  return *this;
2322  }
2327  __INLINE Var & operator()(enumerator i, enumerator j)
2328  {
2329  assert(i < Rows());
2330  assert(j < Cols());
2331  if (i < orow || i >= orow + M->Rows() || j < ocol || j >= ocol + M->Cols())
2332  {
2333  static Var zero(0.0);
2334  assert(zero == 0.0);
2335  return zero;
2336  }
2337  else
2338  return (*M)(i-orow,j-ocol);
2339  }
2345  __INLINE Var compute(enumerator i, enumerator j) const
2346  {
2347  assert(i < Rows());
2348  assert(j < Cols());
2349  if (i < orow || i >= orow + M->Rows() || j < ocol || j >= ocol + M->Cols())
2350  return Var(0.0);
2351  else
2352  return M->compute(i-orow,j-ocol);
2353  }
2361  __INLINE const Var & get(enumerator i, enumerator j) const
2362  {
2363  assert(i < Rows());
2364  assert(j < Cols());
2365  assert(! (i < orow || i >= orow + M->Rows() || j < ocol || j >= ocol + M->Cols()) );
2366  return M->get(i-orow,j-ocol);
2367  }
2373  ::INMOST::Matrix<Var> MakeMatrix()
2374  {
2375  ::INMOST::Matrix<Var> ret(Rows(),Cols());
2376  for(enumerator i = 0; i < Rows(); ++i)
2377  for(enumerator j = 0; j < Cols(); ++j)
2378  ret(i,j) = (*this)(i,j);
2379  return ret;
2380  }
2384  void Resize(enumerator rows, enumerator cols)
2385  {
2386  assert(Cols() == cols);
2387  assert(Rows() == rows);
2388  (void)cols; (void)rows;
2389  }
2390  };
2391 
2393  template<typename Var>
2394  class ConstBlockOfMatrix : public AbstractMatrixReadOnly<Var>
2395  {
2396  public:
2397  using AbstractMatrixReadOnly<Var>::compute;
2398  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2399  private:
2400  const AbstractMatrixReadOnly<Var> * M;
2401  enumerator nrows; //< Number of rows in larger matrix.
2402  enumerator ncols; //< Number of colums in larger matrix.
2403  enumerator orow; //< Row offset in larger matrix.
2404  enumerator ocol; //< Column offset in larger matrix.
2405  public:
2408  __INLINE enumerator Rows() const {return nrows;}
2411  __INLINE enumerator Cols() const {return ncols;}
2418  ConstBlockOfMatrix(const AbstractMatrixReadOnly<Var> & rM, enumerator num_rows, enumerator num_cols, enumerator offset_row, enumerator offset_col) : M(&rM), nrows(num_rows), ncols(num_cols), orow(offset_row), ocol(offset_col)
2419  {}
2420  ConstBlockOfMatrix(const ConstBlockOfMatrix & b) : M(b.M), nrows(b.nrows), ncols(b.ncols), orow(b.orow), ocol(b.ocol) {}
2426  __INLINE Var compute(enumerator i, enumerator j) const
2427  {
2428  assert(i < Rows());
2429  assert(j < Cols());
2430  if (i < orow || i >= orow + M->Rows() || j < ocol || j >= ocol + M->Cols())
2431  return Var(0.0);
2432  else
2433  return M->compute(i-orow,j-ocol);
2434  }
2440  ::INMOST::Matrix<Var> MakeMatrix()
2441  {
2442  ::INMOST::Matrix<Var> ret(Rows(),Cols());
2443  for(enumerator i = 0; i < Rows(); ++i)
2444  for(enumerator j = 0; j < Cols(); ++j)
2445  ret(i,j) = (*this)(i,j);
2446  return ret;
2447  }
2448 
2449  };
2450 
2451  template<typename Var>
2452  class MatrixUnit : public AbstractMatrixReadOnly< Var >
2453  {
2454  public:
2455  using AbstractMatrixReadOnly<Var>::compute;
2456  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2457  private:
2458  enumerator n;
2459  Var c;
2460  public:
2461  MatrixUnit(enumerator n, const Var & c = Var(1.0)) :n(n), c(c) {}
2462  MatrixUnit(const MatrixUnit& b) : n(b.n), c(b.c) {}
2465  __INLINE enumerator Rows() const { return n; }
2468  __INLINE enumerator Cols() const { return n; }
2469  __INLINE Var compute(enumerator i, enumerator j) const
2470  {
2471  return i == j ? c : Var(0.0);
2472  }
2473  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return GetCount(c); }
2474  };
2475 
2476  template<typename Var>
2477  class MatrixRow : public AbstractMatrixReadOnly< Var >
2478  {
2479  public:
2480  using AbstractMatrixReadOnly<Var>::compute;
2481  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2482  private:
2483  enumerator n;
2484  Var c;
2485  public:
2486  MatrixRow(enumerator n, const Var& c = Var(1.0)) :n(n), c(c) {}
2487  MatrixRow(const MatrixRow& b) : n(b.n), c(b.c) {}
2490  __INLINE enumerator Rows() const { return 1; }
2493  __INLINE enumerator Cols() const { return n; }
2494  __INLINE Var compute(enumerator i, enumerator j) const { return c; }
2495  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return GetCount(c); }
2496  };
2497 
2498  template<typename Var>
2499  class MatrixCol : public AbstractMatrixReadOnly< Var >
2500  {
2501  public:
2502  using AbstractMatrixReadOnly<Var>::compute;
2503  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2504  private:
2505  enumerator n;
2506  Var c;
2507  public:
2508  MatrixCol(enumerator n, const Var& c = Var(1.0)) :n(n), c(c) {}
2509  MatrixCol(const MatrixCol& b) : n(b.n), c(b.c) {}
2512  __INLINE enumerator Rows() const { return n; }
2515  __INLINE enumerator Cols() const { return 1; }
2516  __INLINE Var compute(enumerator i, enumerator j) const { return c; }
2517  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return GetCount(c); }
2518  };
2519 
2520  template<typename Var>
2521  class MatrixDiag : public AbstractMatrixReadOnly< Var >
2522  {
2523  public:
2524  using AbstractMatrixReadOnly<Var>::compute;
2525  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2526  private:
2527  enumerator n;
2528  Var* diag;
2529  public:
2530  MatrixDiag(Var* diag, enumerator n) : diag(diag), n(n) {}
2531  MatrixDiag(const MatrixDiag& b) : diag(b.diag), n(b.n) {}
2534  __INLINE enumerator Rows() const { return n; }
2537  __INLINE enumerator Cols() const { return n; }
2538  __INLINE Var compute(enumerator i, enumerator j) const
2539  {
2540  return i == j ? diag[i] : Var(0.0);
2541  }
2542  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
2543  {
2544  INMOST_DATA_ENUM_TYPE cnt = 0;
2545  for (enumerator k = 0; k < n; ++k)
2546  cnt += GetCount(diag[k]);
2547  return cnt;
2548  }
2549  };
2550 
2551  template<typename VarA, typename VarB>
2552  class MatrixSum : public AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >
2553  {
2554  public:
2555  using AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >::compute;
2556  typedef typename AbstractMatrixReadOnly<typename Promote<VarA,VarB>::type >::enumerator enumerator; //< Integer type for indexes.
2557  private:
2558  //static thread_private< typename Promote<VarA, VarB>::type > tmp;
2559  const AbstractMatrixReadOnly<VarA>* A;
2560  const AbstractMatrixReadOnly<VarB>* B;
2561  public:
2564  __INLINE enumerator Rows() const { return A->Rows(); }
2567  __INLINE enumerator Cols() const { return A->Cols(); }
2568  MatrixSum(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
2569  : A(&rA), B(&rB)
2570  {
2571  assert(A->Rows() == B->Rows());
2572  assert(A->Cols() == B->Cols());
2573  }
2574  MatrixSum(const MatrixSum& b) : A(b.A), B(b.B) {}
2580  __INLINE typename Promote<VarA, VarB>::type compute(enumerator i, enumerator j) const
2581  {
2582  return A->compute(i,j) + B->compute(i,j);
2583  }
2584  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2585  };
2586  //template<typename VarA, typename VarB>
2587  //thread_private< typename Promote<VarA, VarB>::type > MatrixSum<VarA,VarB>::tmp;
2588 
2589  template<typename VarA, typename VarB>
2590  class MatrixDifference : public AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >
2591  {
2592  public:
2593  using AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >::compute;
2594  typedef typename AbstractMatrixReadOnly<typename Promote<VarA, VarB>::type>::enumerator enumerator; //< Integer type for indexes.
2595  private:
2596  //static thread_private< typename Promote<VarA, VarB>::type > tmp;
2597  const AbstractMatrixReadOnly<VarA>* A;
2598  const AbstractMatrixReadOnly<VarB>* B;
2599  public:
2602  __INLINE enumerator Rows() const { return A->Rows(); }
2605  __INLINE enumerator Cols() const { return A->Cols(); }
2606  MatrixDifference(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
2607  : A(&rA), B(&rB)
2608  {
2609  assert(A->Rows() == B->Rows());
2610  assert(A->Cols() == B->Cols());
2611  }
2612  MatrixDifference(const MatrixDifference& b) : A(b.A), B(b.B) {}
2618  __INLINE typename Promote<VarA, VarB>::type compute(enumerator i, enumerator j) const
2619  {
2620  return A->compute(i, j) - B->compute(i, j);
2621  }
2622  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2623  };
2624  //template<typename VarA, typename VarB>
2625  //thread_private< typename Promote<VarA, VarB>::type > MatrixDifference<VarA, VarB>::tmp;
2626 
2627  template<typename Var>
2628  class ConstMatrixTranspose : public AbstractMatrixReadOnly<Var>
2629  {
2630  public:
2631  using AbstractMatrixReadOnly<Var>::compute;
2632  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2633  private:
2634  const AbstractMatrixReadOnly<Var>* A;
2635  public:
2638  __INLINE enumerator Rows() const { return A->Cols(); }
2641  __INLINE enumerator Cols() const { return A->Rows(); }
2642  ConstMatrixTranspose(const AbstractMatrixReadOnly<Var>& rA)
2643  : A(&rA) {}
2644  ConstMatrixTranspose(const ConstMatrixTranspose& b) : A(b.A) {}
2650  __INLINE Var compute(enumerator i, enumerator j) const { return A->compute(j, i); }
2651  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
2652  };
2653 
2654  template<typename Var>
2655  class ConstMatrixConjugateTranspose : public AbstractMatrixReadOnly<Var>
2656  {
2657  public:
2658  using AbstractMatrixReadOnly<Var>::compute;
2659  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2660  private:
2661  const AbstractMatrixReadOnly<Var>* A;
2662  //static thread_private<Var> tmp;
2663  public:
2666  __INLINE enumerator Rows() const { return A->Cols(); }
2669  __INLINE enumerator Cols() const { return A->Rows(); }
2670  ConstMatrixConjugateTranspose(const AbstractMatrixReadOnly<Var>& rA)
2671  : A(&rA) {}
2672  ConstMatrixConjugateTranspose(const ConstMatrixConjugateTranspose& b) : A(b.A) {}
2678  __INLINE Var compute(enumerator i, enumerator j) const { return conj(A->compute(j, i)); }
2679  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
2680  };
2681  //template<typename Var>
2682  //thread_private<Var> ConstMatrixConjugateTranspose<Var>::tmp;
2683 
2684  template<typename Var>
2685  class ConstMatrixConjugate : public AbstractMatrixReadOnly<Var>
2686  {
2687  public:
2688  using AbstractMatrixReadOnly<Var>::compute;
2689  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2690  private:
2691  const AbstractMatrixReadOnly<Var>* A;
2692  //static thread_private<Var> tmp;
2693  public:
2696  __INLINE enumerator Rows() const { return A->Rows(); }
2699  __INLINE enumerator Cols() const { return A->Cols(); }
2700  ConstMatrixConjugate(const AbstractMatrixReadOnly<Var>& rA)
2701  : A(&rA) {}
2702  ConstMatrixConjugate(const ConstMatrixConjugate& b) : A(b.A) {}
2708  __INLINE Var compute(enumerator i, enumerator j) const { return conj(A->compute(i, j)); }
2709  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
2710  };
2711  //template<typename Var>
2712  //thread_private<Var> ConstMatrixConjugate<Var>::tmp;
2713 
2714  template<typename Var>
2715  class MatrixUnaryMinus : public AbstractMatrixReadOnly<Var>
2716  {
2717  public:
2718  using AbstractMatrixReadOnly<Var>::compute;
2719  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2720  private:
2721  //static thread_private<Var> tmp;
2722  const AbstractMatrixReadOnly<Var>* A;
2723  public:
2726  __INLINE enumerator Rows() const { return A->Rows(); }
2729  __INLINE enumerator Cols() const { return A->Cols(); }
2730  MatrixUnaryMinus(const AbstractMatrixReadOnly<Var>& rA)
2731  : A(&rA) {}
2732  MatrixUnaryMinus(const MatrixUnaryMinus& b) : A(b.A) {}
2738  __INLINE Var compute(enumerator i, enumerator j) const { return -A->compute(i, j); }
2739  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
2740  };
2741  //template<typename Var>
2742  //thread_private<Var> MatrixUnaryMinus<Var>::tmp;
2743 
2744  template<typename Var>
2745  class MatrixTranspose : public AbstractMatrix<Var>
2746  {
2747  public:
2748  using AbstractMatrix<Var>::operator ();
2749  using AbstractMatrix<Var>::get;
2750  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
2751  private:
2752  AbstractMatrix<Var>* A;
2753  public:
2756  __INLINE enumerator Rows() const { return A->Cols(); }
2759  __INLINE enumerator Cols() const { return A->Rows(); }
2760  MatrixTranspose(AbstractMatrix<Var>& rA)
2761  : A(&rA) {}
2762  MatrixTranspose(const MatrixTranspose& b) : A(b.A) {}
2768  __INLINE Var compute(enumerator i, enumerator j) const { return A->compute(j, i); }
2774  __INLINE Var& operator()(enumerator i, enumerator j) { return (*A)(j, i); }
2780  __INLINE const Var& get(enumerator i, enumerator j) const { return A->get(j, i); }
2784  void Resize(enumerator rows, enumerator cols)
2785  {
2786  assert(Cols() == cols);
2787  assert(Rows() == rows);
2788  (void)cols; (void)rows;
2789  }
2790  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
2791  };
2792 
2793  template<typename Var>
2794  class ConstMatrixConcatRows : public AbstractMatrixReadOnly<Var>
2795  {
2796  public:
2797  using AbstractMatrixReadOnly<Var>::compute;
2798  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2799  private:
2800  const AbstractMatrixReadOnly<Var>* A;
2801  const AbstractMatrixReadOnly<Var>* B;
2802  public:
2805  __INLINE enumerator Rows() const { return A->Rows() + B->Rows(); }
2808  __INLINE enumerator Cols() const { return A->Cols(); }
2809  ConstMatrixConcatRows(const AbstractMatrixReadOnly<Var>& rA, const AbstractMatrixReadOnly<Var>& rB)
2810  : A(&rA), B(&rB) {
2811  assert(A->Cols() == B->Cols());
2812  }
2813  ConstMatrixConcatRows(const ConstMatrixConcatRows& b) : A(b.A), B(b.B) {}
2819  __INLINE Var compute(enumerator i, enumerator j) const
2820  {
2821  return i < A->Rows() ? A->compute(i, j) : B->compute(i - A->Rows(), j);
2822  }
2823  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2824  };
2825 
2826  template<typename VarA, typename VarB, typename VarR>
2827  class ConstMatrixConcatRows2 : public AbstractMatrixReadOnly<VarR>
2828  {
2829  public:
2830  using AbstractMatrixReadOnly<VarR>::compute;
2831  typedef typename AbstractMatrixReadOnly<VarR>::enumerator enumerator; //< Integer type for indexes.
2832  private:
2833  //static thread_private<VarR> tmp;
2834  const AbstractMatrixReadOnly<VarA>* A;
2835  const AbstractMatrixReadOnly<VarB>* B;
2836  public:
2839  __INLINE enumerator Rows() const { return A->Rows() + B->Rows(); }
2842  __INLINE enumerator Cols() const { return A->Cols(); }
2843  ConstMatrixConcatRows2(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
2844  : A(&rA), B(&rB) {
2845  assert(A->Cols() == B->Cols());
2846  }
2847  ConstMatrixConcatRows2(const ConstMatrixConcatRows2& b) : A(b.A), B(b.B) {}
2853  __INLINE VarR compute(enumerator i, enumerator j) const
2854  {
2855  if (i < A->Rows())
2856  return A->compute(i, j);
2857  else
2858  return B->compute(i - A->Rows(), j);
2859  }
2860  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2861  };
2862  //template<typename VarA, typename VarB, typename VarR>
2863  //thread_private<VarR> ConstMatrixConcatRows2<VarA, VarB, VarR>::tmp;
2864 
2865  template<typename Var>
2866  class MatrixConcatRows : public AbstractMatrix<Var>
2867  {
2868  public:
2869  using AbstractMatrix<Var>::operator();
2870  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
2871  private:
2872  AbstractMatrix<Var>* A;
2873  AbstractMatrix<Var>* B;
2874  public:
2877  __INLINE enumerator Rows() const { return A->Rows() + B->Rows(); }
2880  __INLINE enumerator Cols() const { return A->Cols(); }
2881  MatrixConcatRows(AbstractMatrix<Var>& rA, AbstractMatrix<Var>& rB)
2882  : A(&rA), B(&rB) { assert(A->Cols() == B->Cols()); }
2883  MatrixConcatRows(const MatrixConcatRows& b) : A(b.A), B(b.B) {}
2889  __INLINE Var compute(enumerator i, enumerator j) const
2890  {
2891  return i < A->Rows() ? A->compute(i, j) : B->compute(i - A->Rows(), j);
2892  }
2898  __INLINE Var& operator()(enumerator i, enumerator j)
2899  {
2900  return i < A->Rows() ? (*A)(i, j) : (*B)(i - A->Rows(), j);
2901  }
2907  __INLINE const Var& get(enumerator i, enumerator j) const
2908  {
2909  return i < A->Rows() ? A->get(i, j) : B->get(i - A->Rows(), j);
2910  }
2914  void Resize(enumerator rows, enumerator cols)
2915  {
2916  assert(Cols() == cols);
2917  assert(Rows() == rows);
2918  (void)cols; (void)rows;
2919  }
2920  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2921  };
2922 
2923  template<typename Var>
2924  class ConstMatrixConcatCols : public AbstractMatrixReadOnly<Var>
2925  {
2926  public:
2927  using AbstractMatrixReadOnly<Var>::compute;
2928  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
2929  private:
2930  const AbstractMatrixReadOnly<Var>* A;
2931  const AbstractMatrixReadOnly<Var>* B;
2932  public:
2935  __INLINE enumerator Rows() const { return A->Rows(); }
2938  __INLINE enumerator Cols() const { return A->Cols() + B->Cols(); }
2939  ConstMatrixConcatCols(const AbstractMatrixReadOnly<Var>& rA, const AbstractMatrixReadOnly<Var>& rB)
2940  : A(&rA), B(&rB) {
2941  assert(A->Rows() == B->Rows());
2942  }
2943  ConstMatrixConcatCols(const ConstMatrixConcatCols& b) : A(b.A), B(b.B) {}
2949  __INLINE Var compute(enumerator i, enumerator j) const
2950  {
2951  return j < A->Cols() ? A->compute(i, j) : B->compute(i, j - A->Cols());
2952  }
2953  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2954  };
2955 
2956  template<typename VarA, typename VarB, typename VarR>
2957  class ConstMatrixConcatCols2 : public AbstractMatrixReadOnly<VarR>
2958  {
2959  public:
2960  using AbstractMatrixReadOnly<VarR>::compute;
2961  typedef typename AbstractMatrixReadOnly<VarR>::enumerator enumerator; //< Integer type for indexes.
2962  private:
2963  //static thread_private<VarR> tmp;
2964  const AbstractMatrixReadOnly<VarA>* A;
2965  const AbstractMatrixReadOnly<VarB>* B;
2966  public:
2969  __INLINE enumerator Rows() const { return A->Rows(); }
2972  __INLINE enumerator Cols() const { return A->Cols() + B->Cols(); }
2973  ConstMatrixConcatCols2(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
2974  : A(&rA), B(&rB) {
2975  assert(A->Rows() == B->Rows());
2976  }
2977  ConstMatrixConcatCols2(const ConstMatrixConcatCols2& b) : A(b.A), B(b.B) {}
2983  __INLINE VarR compute(enumerator i, enumerator j) const
2984  {
2985  if (j < A->Cols())
2986  return A->compute(i, j);
2987  else
2988  return B->compute(i, j - A->Cols());
2989  }
2990  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
2991  };
2992  //template<typename VarA, typename VarB, typename VarR>
2993  //thread_private<VarR> ConstMatrixConcatCols2<VarA, VarB, VarR>::tmp;
2994 
2995  template<typename Var>
2996  class MatrixConcatCols : public AbstractMatrix<Var>
2997  {
2998  public:
2999  using AbstractMatrix<Var>::operator();
3000  using AbstractMatrix<Var>::get;
3001  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
3002  private:
3003  AbstractMatrix<Var>* A;
3004  AbstractMatrix<Var>* B;
3005  public:
3008  __INLINE enumerator Rows() const { return A->Rows(); }
3011  __INLINE enumerator Cols() const { return A->Cols() + B->Cols(); }
3012  MatrixConcatCols(AbstractMatrix<Var>& rA, AbstractMatrix<Var>& rB)
3013  : A(&rA), B(&rB) { assert(A->Rows() == B->Rows()); }
3014  MatrixConcatCols(const MatrixConcatCols& b) : A(b.A), B(b.B) {}
3020  __INLINE Var compute(enumerator i, enumerator j) const
3021  {
3022  return j < A->Cols() ? A->compute(i, j) : B->compute(i, j - A->Cols());
3023  }
3029  __INLINE Var& operator()(enumerator i, enumerator j)
3030  {
3031  return j < A->Cols() ? (*A)(i, j) : (*B)(i, j - A->Cols());
3032  }
3038  __INLINE const Var& get(enumerator i, enumerator j) const
3039  {
3040  return j < A->Cols() ? A->get(i, j) : B->get(i, j - A->Cols());
3041  }
3045  void Resize(enumerator rows, enumerator cols)
3046  {
3047  assert(Cols() == cols);
3048  assert(Rows() == rows);
3049  (void)cols; (void)rows;
3050  }
3051  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
3052  };
3053 
3054  template<typename Var>
3055  class ConstMatrixRepack : public AbstractMatrixReadOnly<Var>
3056  {
3057  public:
3058  using AbstractMatrixReadOnly<Var>::compute;
3059  typedef typename AbstractMatrixReadOnly<Var>::enumerator enumerator; //< Integer type for indexes.
3060  private:
3061  const AbstractMatrixReadOnly<Var>* A;
3062  enumerator n, m;
3063  public:
3066  __INLINE enumerator Rows() const { return n; }
3069  __INLINE enumerator Cols() const { return m; }
3070  ConstMatrixRepack(const AbstractMatrixReadOnly<Var>& rA, enumerator n, enumerator m)
3071  : A(&rA), n(n), m(m) { assert(A->Rows() * A->Cols() == n * m); }
3072  ConstMatrixRepack(const ConstMatrixRepack& b) : A(b.A), n(b.n), m(b.m) {}
3078  __INLINE Var compute(enumerator i, enumerator j) const
3079  {
3080  enumerator ind = i * m + j;
3081  return A->compute(ind / A->Cols(), ind % A->Cols());
3082  }
3083  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
3084  };
3085 
3086  template<typename Var>
3087  class MatrixRepack : public AbstractMatrix<Var>
3088  {
3089  public:
3090  using AbstractMatrix<Var>::operator();
3091  using AbstractMatrix<Var>::get;
3092  typedef typename AbstractMatrix<Var>::enumerator enumerator; //< Integer type for indexes.
3093  private:
3094  AbstractMatrix<Var>* A;
3095  enumerator n, m;
3096  public:
3099  __INLINE enumerator Rows() const { return n; }
3102  __INLINE enumerator Cols() const { return m; }
3103  MatrixRepack(AbstractMatrix<Var>& rA, enumerator n, enumerator m)
3104  : A(&rA), n(n), m(m) {
3105  assert(A->Rows() * A->Cols() == n * m);
3106  }
3107  MatrixRepack(const MatrixRepack& b) : A(b.A), n(b.n), m(b.m) {}
3113  __INLINE Var compute(enumerator i, enumerator j) const
3114  {
3115  enumerator p = i * m + j;
3116  return A->compute(p / A->Cols(), p % A->Cols());
3117  //return A->data()[i * m + j];
3118  }
3123  __INLINE Var & operator()(enumerator i, enumerator j)
3124  {
3125  enumerator p = i * m + j;
3126  return (*A)(p / A->Cols(), p % A->Cols());
3127  //return A->data()[i * m + j];
3128  }
3134  __INLINE const Var& get(enumerator i, enumerator j) const
3135  {
3136  enumerator p = i * m + j;
3137  return A->get(p / A->Cols(), p % A->Cols());
3138  //return A->data()[i * m + j];
3139  }
3143  void Resize(enumerator rows, enumerator cols)
3144  {
3145  assert(Cols() == cols);
3146  assert(Rows() == rows);
3147  (void)cols; (void)rows;
3148  }
3149  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount(); }
3150  };
3151 
3152  template<typename VarA, typename VarB, typename VarR>
3153  class MatrixMul : public AbstractMatrix< VarR >
3154  {
3155  public:
3156  using AbstractMatrix< VarR >::operator();
3157  using AbstractMatrix< VarR >::get;
3158  typedef typename AbstractMatrix< VarR >::enumerator enumerator; //< Integer type for indexes.
3159  private:
3160  Matrix<VarR> M;
3161  public:
3165  void Resize(enumerator rows, enumerator cols)
3166  {
3167  assert(Cols() == cols);
3168  assert(Rows() == rows);
3169  (void)cols; (void)rows;
3170  }
3173  __INLINE enumerator Rows() const { return M.Rows(); }
3176  __INLINE enumerator Cols() const { return M.Cols(); }
3177  MatrixMul(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
3178  {
3179  assert(rA.Cols() == rB.Rows());
3180  const AbstractMatrix<VarA>* pA = dynamic_cast<const AbstractMatrix<VarA> *>(&rA);
3181  const AbstractMatrix<VarB>* pB = dynamic_cast<const AbstractMatrix<VarB> *>(&rB);
3182  if (!pA)
3183  {
3184  static thread_private< Matrix<VarA> > tmpA;
3185  *tmpA = rA;
3186  pA = &(*tmpA);
3187  }
3188  if (!pB)
3189  {
3190  static thread_private< Matrix<VarB> > tmpB;
3191  *tmpB = rB;
3192  pB = &(*tmpB);
3193  }
3194  M.Resize(rA.Rows(), rB.Cols());
3195  for (enumerator i = 0; i < pA->Rows(); ++i)
3196  {
3197  for (enumerator k = 0; k < pB->Cols(); ++k)
3198  M(i, k) = pA->get(i, 0) * pB->get(0, k);
3199  for (enumerator j = 1; j < pA->Cols(); ++j)
3200  for (enumerator k = 0; k < pB->Cols(); ++k)
3201  M(i, k) += pA->get(i, j) * pB->get(j, k);
3202  }
3203  }
3204  MatrixMul(const MatrixMul& b) : M(b.M) {}
3210  __INLINE VarR compute(enumerator i, enumerator j) const { return M.get(i, j); }
3215  __INLINE VarR& operator()(enumerator i, enumerator j) { return M(i, j); }
3221  __INLINE const VarR& get(enumerator i, enumerator j) const { return M.get(i, j); }
3222  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return M.GetMatrixCount(); }
3223  };
3224 #if defined(USE_AUTODIFF)
3225  template<>
3226  class MatrixMul<INMOST_DATA_REAL_TYPE, variable, Promote<INMOST_DATA_REAL_TYPE,variable>::type > : public AbstractMatrix< Promote<INMOST_DATA_REAL_TYPE, variable>::type >
3227  {
3228  public:
3229  using AbstractMatrix< Promote<INMOST_DATA_REAL_TYPE, variable>::type >::operator();
3230  typedef typename AbstractMatrix< Promote<INMOST_DATA_REAL_TYPE, variable>::type >::enumerator enumerator; //< Integer type for indexes.
3231  private:
3232  Matrix<Promote<INMOST_DATA_REAL_TYPE, variable>::type> M;
3233  public:
3237  void Resize(enumerator rows, enumerator cols)
3238  {
3239  assert(Cols() == cols);
3240  assert(Rows() == rows);
3241  (void)cols; (void)rows;
3242  }
3245  __INLINE enumerator Rows() const { return M.Rows(); }
3248  __INLINE enumerator Cols() const { return M.Cols(); }
3249  MatrixMul(const AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE>& rA, const AbstractMatrixReadOnly<variable>& rB)
3250  {
3251  assert(rA.Cols() == rB.Rows());
3252  const AbstractMatrix<INMOST_DATA_REAL_TYPE>* pA = dynamic_cast<const AbstractMatrix<INMOST_DATA_REAL_TYPE> *>(&rA);
3253  const AbstractMatrix<variable>* pB = dynamic_cast<const AbstractMatrix<variable> *>(&rB);
3254  if (!pA)
3255  {
3256  static thread_private< Matrix<INMOST_DATA_REAL_TYPE> > tmpA;
3257  *tmpA = rA;
3258  pA = &(*tmpA);
3259  }
3260  if (!pB)
3261  {
3262  static thread_private< Matrix<variable> > tmpB;
3263  *tmpB = rB;
3264  pB = &(*tmpB);
3265  }
3266  M.Resize(rA.Rows(), rB.Cols());
3267  INMOST_DATA_ENUM_TYPE cnt = pB->GetMatrixCount();
3268  if (cnt >= CNT_USE_MERGER)
3269  {
3270  merger->Resize(cnt);
3271  for (enumerator i = 0; i < pA->Rows(); ++i)
3272  {
3273  for (enumerator j = 0; j < pB->Cols(); ++j)
3274  {
3275  INMOST_DATA_REAL_TYPE value = 0.0;
3276  for (enumerator k = 0; k < pA->Cols(); ++k)
3277  {
3278  value += pA->get(i, k) * pB->get(k, j).GetValue();
3279  merger->AddRow(pA->get(i, k), pB->get(k, j).GetRow());
3280  }
3281  M(i, j).SetValue(value);
3282  merger->RetrieveRow(M(i, j).GetRow());
3283  merger->Clear();
3284  }
3285  }
3286  }
3287  else
3288  {
3289  for (enumerator i = 0; i < pA->Rows(); ++i)
3290  {
3291  for (enumerator k = 0; k < pB->Cols(); ++k)
3292  M(i, k) = pA->get(i, 0) * pB->get(0, k);
3293  for (enumerator j = 1; j < pA->Cols(); ++j)
3294  for (enumerator k = 0; k < pB->Cols(); ++k)
3295  M(i, k) += pA->get(i, j) * pB->get(j, k);
3296  }
3297  }
3298  }
3299  MatrixMul(const MatrixMul& b) : M(b.M) {}
3305  __INLINE variable compute(enumerator i, enumerator j) const { return M.get(i, j); }
3306  __INLINE variable& operator()(enumerator i, enumerator j) { return M(i, j); }
3307  __INLINE const variable& get(enumerator i, enumerator j) const { return M.get(i, j); }
3308  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return M.GetMatrixCount(); }
3309  };
3310 
3311  template<>
3312  class MatrixMul<variable, INMOST_DATA_REAL_TYPE, Promote<variable, INMOST_DATA_REAL_TYPE>::type > : public AbstractMatrix< Promote<variable, INMOST_DATA_REAL_TYPE>::type >
3313  {
3314  public:
3315  using AbstractMatrix< Promote<variable, INMOST_DATA_REAL_TYPE>::type >::operator();
3316  typedef typename AbstractMatrix< Promote<variable, INMOST_DATA_REAL_TYPE>::type >::enumerator enumerator; //< Integer type for indexes.
3317  private:
3318  Matrix<Promote<variable, INMOST_DATA_REAL_TYPE>::type> M;
3319  public:
3323  void Resize(enumerator rows, enumerator cols)
3324  {
3325  assert(Cols() == cols);
3326  assert(Rows() == rows);
3327  (void)cols; (void)rows;
3328  }
3331  __INLINE enumerator Rows() const { return M.Rows(); }
3334  __INLINE enumerator Cols() const { return M.Cols(); }
3335  MatrixMul(const AbstractMatrixReadOnly<variable>& rA, const AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE>& rB)
3336  {
3337  assert(rA.Cols() == rB.Rows());
3338  const AbstractMatrix<variable>* pA = dynamic_cast<const AbstractMatrix<variable> *>(&rA);
3339  const AbstractMatrix<INMOST_DATA_REAL_TYPE>* pB = dynamic_cast<const AbstractMatrix<INMOST_DATA_REAL_TYPE> *>(&rB);
3340  if (!pA)
3341  {
3342  static thread_private< Matrix<variable> > tmpA;
3343  *tmpA = rA;
3344  pA = &(*tmpA);
3345  }
3346  if (!pB)
3347  {
3348  static thread_private< Matrix<INMOST_DATA_REAL_TYPE> > tmpB;
3349  *tmpB = rB;
3350  pB = &(*tmpB);
3351  }
3352  M.Resize(rA.Rows(), rB.Cols());
3353  INMOST_DATA_ENUM_TYPE cnt = pA->GetMatrixCount();
3354  if (cnt >= CNT_USE_MERGER)
3355  {
3356  merger->Resize(cnt);
3357  for (enumerator i = 0; i < pA->Rows(); ++i)
3358  {
3359  for (enumerator j = 0; j < pB->Cols(); ++j)
3360  {
3361  INMOST_DATA_REAL_TYPE value = 0.0;
3362  for (enumerator k = 0; k < pA->Cols(); ++k)
3363  {
3364  value += pA->get(i, k).GetValue() * pB->get(k, j);
3365  merger->AddRow(pB->get(k, j), pA->get(i, k).GetRow());
3366  }
3367  M(i, j).SetValue(value);
3368  merger->RetrieveRow(M(i, j).GetRow());
3369  merger->Clear();
3370  }
3371  }
3372  }
3373  else
3374  {
3375  for (enumerator i = 0; i < pA->Rows(); ++i)
3376  {
3377  for (enumerator k = 0; k < pB->Cols(); ++k)
3378  M(i, k) = pA->get(i, 0) * pB->get(0, k);
3379  for (enumerator j = 1; j < pA->Cols(); ++j)
3380  for (enumerator k = 0; k < pB->Cols(); ++k)
3381  M(i, k) += pA->get(i, j) * pB->get(j, k);
3382  }
3383  }
3384  }
3385  MatrixMul(const MatrixMul& b) : M(b.M) {}
3391  __INLINE variable compute(enumerator i, enumerator j) const { return M.get(i, j); }
3392  __INLINE variable & operator()(enumerator i, enumerator j) { return M(i, j); }
3393  __INLINE const variable& get(enumerator i, enumerator j) const { return M.get(i, j); }
3394  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return M.GetMatrixCount(); }
3395  };
3396 
3397  template<>
3398  class MatrixMul<variable, variable, Promote<variable, variable>::type > : public AbstractMatrix< Promote<variable, variable>::type >
3399  {
3400  public:
3401  using AbstractMatrix< Promote<variable, variable>::type >::operator();
3402  typedef typename AbstractMatrix< Promote<variable, variable>::type >::enumerator enumerator; //< Integer type for indexes.
3403  private:
3404  Matrix<Promote<variable, variable>::type> M;
3405  public:
3409  void Resize(enumerator rows, enumerator cols)
3410  {
3411  assert(Cols() == cols);
3412  assert(Rows() == rows);
3413  (void)cols; (void)rows;
3414  }
3417  __INLINE enumerator Rows() const { return M.Rows(); }
3420  __INLINE enumerator Cols() const { return M.Cols(); }
3421  MatrixMul(const AbstractMatrixReadOnly<variable>& rA, const AbstractMatrixReadOnly<variable>& rB)
3422  {
3423  assert(rA.Cols() == rB.Rows());
3424  const AbstractMatrix<variable>* pA = dynamic_cast<const AbstractMatrix<variable>*>(&rA);
3425  const AbstractMatrix<variable>* pB = dynamic_cast<const AbstractMatrix<variable>*>(&rB);
3426  if (!pA)
3427  {
3428  static thread_private< Matrix<variable> > tmpA;
3429  *tmpA = rA;
3430  pA = &(*tmpA);
3431  }
3432  if (!pB)
3433  {
3434  static thread_private< Matrix<variable> > tmpB;
3435  *tmpB = rB;
3436  pB = &(*tmpB);
3437  }
3438  M.Resize(rA.Rows(), rB.Cols());
3439  INMOST_DATA_ENUM_TYPE cnt = 0;
3440  cnt += pA->GetMatrixCount();
3441  cnt += pB->GetMatrixCount();
3442  if (cnt >= CNT_USE_MERGER)
3443  {
3444  merger->Resize(cnt);
3445  for (enumerator i = 0; i < pA->Rows(); ++i)
3446  {
3447  for (enumerator j = 0; j < pB->Cols(); ++j)
3448  {
3449  INMOST_DATA_REAL_TYPE value = 0.0;
3450  for (enumerator k = 0; k < pA->Cols(); ++k)
3451  {
3452  value += pA->get(i, k).GetValue() * pB->get(k, j).GetValue();
3453  merger->AddRow(pA->get(i, k).GetValue(), pB->get(k, j).GetRow());
3454  merger->AddRow(pB->get(k, j).GetValue(), pA->get(i, k).GetRow());
3455  }
3456  M(i, j).SetValue(value);
3457  merger->RetrieveRow(M(i, j).GetRow());
3458  merger->Clear();
3459  }
3460  }
3461  }
3462  else
3463  {
3464  for (enumerator i = 0; i < pA->Rows(); ++i)
3465  {
3466  for (enumerator k = 0; k < pB->Cols(); ++k)
3467  M(i, k) = pA->get(i, 0) * pB->get(0, k);
3468  for (enumerator j = 1; j < pA->Cols(); ++j)
3469  for (enumerator k = 0; k < pB->Cols(); ++k)
3470  M(i, k) += pA->get(i, j) * pB->get(j, k);
3471  }
3472  }
3473  }
3474  MatrixMul(const MatrixMul& b) : M(b.M) {}
3480  __INLINE variable compute(enumerator i, enumerator j) const { return M.get(i, j); }
3481  __INLINE variable & operator()(enumerator i, enumerator j) { return M(i, j); }
3482  __INLINE const variable& get(enumerator i, enumerator j) const { return M.get(i, j); }
3483  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return M.GetMatrixCount(); }
3484  };
3485 #endif //USE_AUTODIFF
3486  //template<typename VarA, typename VarB, typename VarR>
3487  //thread_private< Matrix<VarR> > MatrixMul<VarA, VarB, VarR>::M;
3488 
3489  template<typename VarA, typename VarB, typename VarR>
3490  class MatrixMulCoef : public AbstractMatrixReadOnly< VarR >
3491  {
3492  public:
3493  using AbstractMatrixReadOnly< VarR >::compute;
3494  typedef typename AbstractMatrixReadOnly< VarR >::enumerator enumerator; //< Integer type for indexes.
3495  private:
3496  //static thread_private<VarR> tmp;
3497  const AbstractMatrixReadOnly<VarA>* A;
3498  const VarB* coef;
3499  public:
3502  bool TrivialArguments() const { return false; };
3505  __INLINE enumerator Rows() const { return A->Rows(); }
3508  __INLINE enumerator Cols() const { return A->Cols(); }
3509  MatrixMulCoef(const AbstractMatrixReadOnly<VarA>& rA, const VarB& rcoef)
3510  : A(&rA), coef(&rcoef) {}
3511  MatrixMulCoef(const MatrixMulCoef& b) : A(b.A), coef(b.coef) {}
3517  __INLINE VarR compute(enumerator i, enumerator j) const { return A->compute(i, j) * (*coef); }
3518  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + GetCount(*coef); }
3519  };
3520  //template<typename VarA, typename VarB, typename VarR>
3521  //thread_private<VarR> MatrixMulCoef<VarA, VarB, VarR>::tmp;
3522 
3523  template<typename VarA, typename VarB, typename VarR>
3524  class MatrixDivCoef : public AbstractMatrixReadOnly< VarR >
3525  {
3526  public:
3527  using AbstractMatrixReadOnly< VarR >::compute;
3528  typedef typename AbstractMatrixReadOnly< VarR >::enumerator enumerator; //< Integer type for indexes.
3529  private:
3530  //static thread_private<VarR> tmp;
3531  const AbstractMatrixReadOnly<VarA>* A;
3532  const VarB* coef;
3533  public:
3536  bool TrivialArguments() const { return false; };
3539  __INLINE enumerator Rows() const { return A->Rows(); }
3542  __INLINE enumerator Cols() const { return A->Cols(); }
3543  MatrixDivCoef(const AbstractMatrixReadOnly<VarA>& rA, const VarB& rcoef)
3544  : A(&rA), coef(&rcoef) {}
3545  MatrixDivCoef(const MatrixDivCoef& b)
3546  : A(b.A), coef(b.coef) {}
3552  __INLINE VarR compute(enumerator i, enumerator j) const { return A->compute(i, j) / (*coef); }
3553  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + GetCount(*coef); }
3554  };
3555  //template<typename VarA, typename VarB, typename VarR>
3556  //thread_private<VarR> MatrixDivCoef<VarA, VarB, VarR>::tmp;
3557 #if defined(USE_AUTODIFF)
3558  template<typename VarA, typename VarB, typename VarR>
3559  class MatrixMulShellCoef : public AbstractMatrixReadOnly< VarR >
3560  {
3561  public:
3562  using AbstractMatrixReadOnly< VarR >::compute;
3563  typedef typename AbstractMatrixReadOnly< VarR >::enumerator enumerator; //< Integer type for indexes.
3564  private:
3565  //static thread_private<VarR> tmp;
3566  const AbstractMatrixReadOnly<VarA>* A;
3567  const variable coef;
3568  public:
3571  bool TrivialArguments() const { return false; };
3574  __INLINE enumerator Rows() const { return A->Rows(); }
3577  __INLINE enumerator Cols() const { return A->Cols(); }
3578  MatrixMulShellCoef(const AbstractMatrixReadOnly<VarA>& rA, const VarB& rcoef)
3579  : A(&rA), coef(rcoef) {}
3580  MatrixMulShellCoef(const MatrixMulShellCoef& b) : A(b.A), coef(b.coef) {}
3586  __INLINE VarR compute(enumerator i, enumerator j) const { return A->compute(i, j) * coef; }
3587  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + GetCount(coef); }
3588  };
3589  //template<typename VarA, typename VarB, typename VarR>
3590  //thread_private<VarR> MatrixMulShellCoef<VarA, VarB, VarR>::tmp;
3591 
3592  template<typename VarA, typename VarB, typename VarR>
3593  class MatrixDivShellCoef : public AbstractMatrixReadOnly< VarR >
3594  {
3595  public:
3596  using AbstractMatrixReadOnly< VarR >::compute;
3597  typedef typename AbstractMatrixReadOnly< VarR >::enumerator enumerator; //< Integer type for indexes.
3598  private:
3599  //static thread_private<VarR> tmp;
3600  const AbstractMatrixReadOnly<VarA>* A;
3601  const variable coef;
3602  public:
3605  bool TrivialArguments() const { return false; };
3608  __INLINE enumerator Rows() const { return A->Rows(); }
3611  __INLINE enumerator Cols() const { return A->Cols(); }
3612  MatrixDivShellCoef(const AbstractMatrixReadOnly<VarA>& rA, const VarB& rcoef)
3613  : A(&rA), coef(rcoef) {}
3614  MatrixDivShellCoef(const MatrixDivShellCoef& b)
3615  : A(b.A), coef(b.coef) {}
3621  __INLINE VarR compute(enumerator i, enumerator j) const { return (A->compute(i, j) / coef); }
3622  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + GetCount(coef); }
3623  };
3624  //template<typename VarA, typename VarB, typename VarR>
3625  //thread_private<VarR> MatrixDivShellCoef<VarA, VarB, VarR>::tmp;
3626 #endif //USE_AUTODIFF
3627  template<typename VarA, typename VarB>
3628  class KroneckerProduct : public AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >
3629  {
3630  public:
3631  using AbstractMatrixReadOnly< typename Promote<VarA, VarB>::type >::compute;
3632  typedef typename AbstractMatrixReadOnly<typename Promote<VarA, VarB>::type>::enumerator enumerator; //< Integer type for indexes.
3633  private:
3634  //static thread_private< typename Promote<VarA, VarB>::type > tmp;
3635  const AbstractMatrixReadOnly<VarA>* A;
3636  const AbstractMatrixReadOnly<VarB>* B;
3637  public:
3640  bool TrivialArguments() const { return false; };
3643  __INLINE enumerator Rows() const { return A->Rows() * B->Rows(); }
3646  __INLINE enumerator Cols() const { return A->Cols() * B->Cols(); }
3647  KroneckerProduct(const AbstractMatrixReadOnly<VarA>& rA, const AbstractMatrixReadOnly<VarB>& rB)
3648  : A(&rA), B(&rB) {}
3649  KroneckerProduct(const KroneckerProduct& b) : A(b.A), B(b.B) {}
3655  __INLINE typename Promote<VarA, VarB>::type compute(enumerator i, enumerator j) const
3656  {
3657  return A->compute(i / B->Rows(), j / B->Cols()) * B->compute(i % B->Rows(), j % B->Cols());
3658  }
3659  __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const { return A->GetMatrixCount() + B->GetMatrixCount(); }
3660  };
3661  //template<typename VarA, typename VarB>
3662  //thread_private<typename Promote<VarA, VarB>::type> KroneckerProduct<VarA, VarB>::tmp;
3663 
3664  template<typename Var>
3665  void
3666  AbstractMatrix<Var>::Swap(AbstractMatrix<Var> & b)
3667  {
3668  Matrix<Var> tmp = b;
3669  b = (*this);
3670  (*this) = tmp;
3671  }
3672 
3673  template<typename Var>
3674  template<typename typeB>
3675  Matrix<typename Promote<Var,typeB>::type>
3676  AbstractMatrixReadOnly<Var>::CrossProduct(const AbstractMatrixReadOnly<typeB> & B) const
3677  {
3678  const AbstractMatrixReadOnly<Var> & A = *this;
3679  assert(A.Rows() == 3);
3680  assert(A.Cols() == 1);
3681  assert(B.Rows() == 3);
3682  Matrix<typename Promote<Var,typeB>::type> ret(3,B.Cols()); //check RVO
3683  for(unsigned k = 0; k < B.Cols(); ++k)
3684  {
3685  ret(0,k) = (A.compute(1,0)*B.compute(2,k) - A.compute(2,0)*B.compute(1,k));
3686  ret(1,k) = (A.compute(2,0)*B.compute(0,k) - A.compute(0,0)*B.compute(2,k));
3687  ret(2,k) = (A.compute(0,0)*B.compute(1,k) - A.compute(1,0)*B.compute(0,k));
3688  }
3689  return ret;
3690  }
3691 
3692 
3693  template<typename Var>
3694  template<typename typeB>
3695  Matrix<typename Promote<Var, typeB>::type>
3696  AbstractMatrix<Var>::CrossProduct(const AbstractMatrix<typeB>& B) const
3697  {
3698  const AbstractMatrix<Var>& A = *this;
3699  assert(A.Rows() == 3);
3700  assert(A.Cols() == 1);
3701  assert(B.Rows() == 3);
3702  Matrix<typename Promote<Var, typeB>::type> ret(3, B.Cols()); //check RVO
3703  for (unsigned k = 0; k < B.Cols(); ++k)
3704  {
3705  ret(0, k) = (A.get(1, 0) * B.get(2, k) - A.get(2, 0) * B.get(1, k));
3706  ret(1, k) = (A.get(2, 0) * B.get(0, k) - A.get(0, 0) * B.get(2, k));
3707  ret(2, k) = (A.get(0, 0) * B.get(1, k) - A.get(1, 0) * B.get(0, k));
3708  }
3709  return ret;
3710  }
3711 
3712  template<typename Var>
3713  template<typename typeB>
3714  Matrix<typename Promote<Var, typeB>::type>
3715  AbstractMatrix<Var>::CrossProduct(const AbstractMatrixReadOnly<typeB>& B) const
3716  {
3717  const AbstractMatrix<Var>& A = *this;
3718  assert(A.Rows() == 3);
3719  assert(A.Cols() == 1);
3720  assert(B.Rows() == 3);
3721  Matrix<typename Promote<Var, typeB>::type> ret(3, B.Cols()); //check RVO
3722  for (unsigned k = 0; k < B.Cols(); ++k)
3723  {
3724  ret(0, k) = (A.get(1, 0) * B.compute(2, k) - A.get(2, 0) * B.compute(1, k));
3725  ret(1, k) = (A.get(2, 0) * B.compute(0, k) - A.get(0, 0) * B.compute(2, k));
3726  ret(2, k) = (A.get(0, 0) * B.compute(1, k) - A.get(1, 0) * B.compute(0, k));
3727  }
3728  return ret;
3729  }
3730 
3731 
3732  template<typename Var>
3733  template<typename typeB>
3734  Matrix<typename Promote<Var,typeB>::type>
3735  AbstractMatrixReadOnly<Var>::Transform(const AbstractMatrixReadOnly<typeB> & B) const
3736  {
3737  const AbstractMatrixReadOnly<Var> & A = *this;
3738  assert(A.Rows() == 3);
3739  assert(A.Cols() == 1);
3740  assert(B.Rows() == 3);
3741  assert(B.Cols() == 1);
3742  typedef typename Promote<Var,typeB>::type var_t;
3743  Matrix<var_t> Q(3,3);
3744  var_t x,y,z,w,n,l;
3745 
3746  x = -(B.compute(1,0)*A.compute(2,0) - B.compute(2,0)*A.compute(1,0));
3747  y = -(B.compute(2,0)*A.compute(0,0) - B.compute(0,0)*A.compute(2,0));
3748  z = -(B.compute(0,0)*A.compute(1,0) - B.compute(1,0)*A.compute(0,0));
3749  w = A.FrobeniusNorm()*B.FrobeniusNorm() + A.DotProduct(B);
3750 
3751  l = B.FrobeniusNorm()/A.FrobeniusNorm();
3752  n = 2*l/(x*x+y*y+z*z+w*w);
3753 
3754  Q(0,0) = l - (y*y + z*z)*n;
3755  Q(0,1) = (x*y - z*w)*n;
3756  Q(0,2) = (x*z + y*w)*n;
3757 
3758  Q(1,0) = (x*y + z*w)*n;
3759  Q(1,1) = l - (x*x + z*z)*n;
3760  Q(1,2) = (y*z - x*w)*n;
3761 
3762  Q(2,0) = (x*z - y*w)*n;
3763  Q(2,1) = (y*z + x*w)*n;
3764  Q(2,2) = l - (x*x + y*y)*n;
3765 
3766  return Q;
3767  }
3768 
3769 
3770 
3771  template<typename Var>
3772  template<typename typeB>
3773  Matrix<typename Promote<Var, typeB>::type>
3774  AbstractMatrix<Var>::Transform(const AbstractMatrix<typeB>& B) const
3775  {
3776  const AbstractMatrix<Var>& A = *this;
3777  assert(A.Rows() == 3);
3778  assert(A.Cols() == 1);
3779  assert(B.Rows() == 3);
3780  assert(B.Cols() == 1);
3781  typedef typename Promote<Var, typeB>::type var_t;
3782  Matrix<var_t> Q(3, 3);
3783  var_t x, y, z, w, n, l;
3784 
3785  x = -(B.get(1, 0) * A.get(2, 0) - B.get(2, 0) * A.get(1, 0));
3786  y = -(B.get(2, 0) * A.get(0, 0) - B.get(0, 0) * A.get(2, 0));
3787  z = -(B.get(0, 0) * A.get(1, 0) - B.get(1, 0) * A.get(0, 0));
3788  w = A.FrobeniusNorm() * B.FrobeniusNorm() + A.DotProduct(B);
3789 
3790  l = B.FrobeniusNorm() / A.FrobeniusNorm();
3791  n = 2 * l / (x * x + y * y + z * z + w * w);
3792 
3793  Q(0, 0) = l - (y * y + z * z) * n;
3794  Q(0, 1) = (x * y - z * w) * n;
3795  Q(0, 2) = (x * z + y * w) * n;
3796 
3797  Q(1, 0) = (x * y + z * w) * n;
3798  Q(1, 1) = l - (x * x + z * z) * n;
3799  Q(1, 2) = (y * z - x * w) * n;
3800 
3801  Q(2, 0) = (x * z - y * w) * n;
3802  Q(2, 1) = (y * z + x * w) * n;
3803  Q(2, 2) = l - (x * x + y * y) * n;
3804 
3805  return Q;
3806  }
3807 
3808 
3809  template<typename Var>
3810  template<typename typeB>
3811  Matrix<typename Promote<Var, typeB>::type>
3812  AbstractMatrix<Var>::Transform(const AbstractMatrixReadOnly<typeB>& B) const
3813  {
3814  const AbstractMatrix<Var>& A = *this;
3815  assert(A.Rows() == 3);
3816  assert(A.Cols() == 1);
3817  assert(B.Rows() == 3);
3818  assert(B.Cols() == 1);
3819  typedef typename Promote<Var, typeB>::type var_t;
3820  Matrix<var_t> Q(3, 3);
3821  var_t x, y, z, w, n, l;
3822 
3823  x = -(B.compute(1, 0) * A.get(2, 0) - B.compute(2, 0) * A.get(1, 0));
3824  y = -(B.compute(2, 0) * A.get(0, 0) - B.compute(0, 0) * A.get(2, 0));
3825  z = -(B.compute(0, 0) * A.get(1, 0) - B.compute(1, 0) * A.get(0, 0));
3826  w = A.FrobeniusNorm() * B.FrobeniusNorm() + A.DotProduct(B);
3827 
3828  l = B.FrobeniusNorm() / A.FrobeniusNorm();
3829  n = 2 * l / (x * x + y * y + z * z + w * w);
3830 
3831  Q(0, 0) = l - (y * y + z * z) * n;
3832  Q(0, 1) = (x * y - z * w) * n;
3833  Q(0, 2) = (x * z + y * w) * n;
3834 
3835  Q(1, 0) = (x * y + z * w) * n;
3836  Q(1, 1) = l - (x * x + z * z) * n;
3837  Q(1, 2) = (y * z - x * w) * n;
3838 
3839  Q(2, 0) = (x * z - y * w) * n;
3840  Q(2, 1) = (y * z + x * w) * n;
3841  Q(2, 2) = l - (x * x + y * y) * n;
3842 
3843  return Q;
3844  }
3845 
3846  template<typename Var>
3847  template<typename typeB>
3848  MatrixDifference<Var,typeB>
3849  AbstractMatrixReadOnly<Var>::operator-(const AbstractMatrixReadOnly<typeB> & other) const
3850  {
3851  return MatrixDifference<Var, typeB>(*this, other);
3852  }
3853 
3854  template<typename Var>
3855  template<typename typeB>
3856  AbstractMatrix<Var>&
3857  AbstractMatrix<Var>::operator-=(const AbstractMatrix<typeB>& other)
3858  {
3859  assert(Rows() == other.Rows());
3860  assert(Cols() == other.Cols());
3861  for (enumerator i = 0; i < Rows(); ++i)
3862  for (enumerator j = 0; j < Cols(); ++j)
3863  assign((*this)(i, j), (*this)(i, j) - other.get(i, j));
3864  return *this;
3865  }
3866 
3867  template<typename Var>
3868  template<typename typeB>
3869  AbstractMatrix<Var> &
3870  AbstractMatrix<Var>::operator-=(const AbstractMatrixReadOnly<typeB> & other)
3871  {
3872  assert(Rows() == other.Rows());
3873  assert(Cols() == other.Cols());
3874  for(enumerator i = 0; i < Rows(); ++i)
3875  for(enumerator j = 0; j < Cols(); ++j)
3876  assign((*this)(i,j),(*this)(i,j)-other.compute(i,j));
3877  return *this;
3878  }
3879 
3880  template<typename Var>
3881  template<typename typeB>
3882  MatrixSum<Var,typeB>
3883  AbstractMatrixReadOnly<Var>::operator+(const AbstractMatrixReadOnly<typeB> & other) const
3884  {
3885  return MatrixSum<Var, typeB>(*this, other);
3886  }
3887 
3888  template<typename Var>
3889  template<typename typeB>
3890  AbstractMatrix<Var>&
3891  AbstractMatrix<Var>::operator+=(const AbstractMatrix<typeB>& other)
3892  {
3893  assert(Rows() == other.Rows());
3894  assert(Cols() == other.Cols());
3895  for (enumerator i = 0; i < Rows(); ++i)
3896  for (enumerator j = 0; j < Cols(); ++j)
3897  assign((*this)(i, j), (*this)(i, j) + other.get(i, j));
3898  return *this;
3899  }
3900 
3901  template<typename Var>
3902  template<typename typeB>
3903  AbstractMatrix<Var> &
3904  AbstractMatrix<Var>::operator+=(const AbstractMatrixReadOnly<typeB> & other)
3905  {
3906  assert(Rows() == other.Rows());
3907  assert(Cols() == other.Cols());
3908  for(enumerator i = 0; i < Rows(); ++i)
3909  for(enumerator j = 0; j < Cols(); ++j)
3910  assign((*this)(i,j),(*this)(i,j)+other.compute(i,j));
3911  return *this;
3912  }
3913 
3914  template<typename Var>
3915  template<typename typeB>
3916  MatrixMul<Var, typeB, typename Promote<Var, typeB>::type>
3917  AbstractMatrixReadOnly<Var>::operator*(const AbstractMatrixReadOnly<typeB> & other) const
3918  {
3919  return MatrixMul<Var, typeB, typename Promote<Var, typeB>::type>(*this, other);
3920  }
3921 #if 0
3922 #if defined(USE_AUTODIFF)
3923  template<>
3924  template<>
3925  __INLINE Promote<INMOST_DATA_REAL_TYPE,variable>::type
3926  AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE>::DotProduct<variable>(const AbstractMatrixReadOnly<variable> & other) const
3927  {
3928  assert(Cols() == other.Cols());
3929  assert(Rows() == other.Rows());
3930  Promote<INMOST_DATA_REAL_TYPE,variable>::type ret = 0.0;
3931  /*
3932  INMOST_DATA_ENUM_TYPE cnt = other.GetMatrixCount();
3933  if( cnt >= CNT_USE_MERGER )
3934  {
3935  merger->Resize(cnt);
3936  double value = 0.0;
3937  for(enumerator i = 0; i < Rows(); ++i)
3938  for(enumerator j = 0; j < Cols(); ++j)
3939  {
3940  value += (*this)(i,j)*other(i,j).GetValue();
3941  merger->AddRow((*this)(i,j),other(i,j).GetRow());
3942  }
3943  ret.SetValue(value);
3944  merger->RetrieveRow(ret.GetRow());
3945  merger->Clear();
3946  }
3947  else*/
3948  {
3949  for(enumerator i = 0; i < Rows(); ++i)
3950  for(enumerator j = 0; j < Cols(); ++j)
3951  ret += compute(i, j) * other.compute(i, j);
3952  }
3953  return ret;
3954  }
3955 
3956  template<>
3957  template<>
3958  __INLINE Promote<variable,INMOST_DATA_REAL_TYPE>::type
3959  AbstractMatrixReadOnly<variable>::DotProduct<INMOST_DATA_REAL_TYPE>(const AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE> & other) const
3960  {
3961  assert(Cols() == other.Cols());
3962  assert(Rows() == other.Rows());
3963  Promote<variable,INMOST_DATA_REAL_TYPE>::type ret = 0.0;
3964  /*
3965  INMOST_DATA_ENUM_TYPE cnt = GetMatrixCount();
3966  if( cnt >= CNT_USE_MERGER )
3967  {
3968  merger->Resize(cnt);
3969  double value = 0.0;
3970  for(enumerator i = 0; i < Rows(); ++i)
3971  for(enumerator j = 0; j < Cols(); ++j)
3972  {
3973  value += (*this)(i,j).GetValue()*other(i,j);
3974  merger->AddRow(other(i,j),(*this)(i,j).GetRow());
3975  }
3976  ret.SetValue(value);
3977  merger->RetrieveRow(ret.GetRow());
3978  merger->Clear();
3979  }
3980  else*/
3981  {
3982  for(enumerator i = 0; i < Rows(); ++i)
3983  for(enumerator j = 0; j < Cols(); ++j)
3984  ret += compute(i,j)*other.compute(i,j);
3985  }
3986  return ret;
3987  }
3988 
3989  template<>
3990  template<>
3991  __INLINE Promote<variable,variable>::type
3992  AbstractMatrixReadOnly<variable>::DotProduct<variable>(const AbstractMatrixReadOnly<variable> & other) const
3993  {
3994  assert(Cols() == other.Cols());
3995  assert(Rows() == other.Rows());
3996  Promote<variable,variable>::type ret = 0.0;
3997  /*
3998  INMOST_DATA_ENUM_TYPE cnt = GetMatrixCount() + other.GetMatrixCount();
3999  if( cnt >= CNT_USE_MERGER )
4000  {
4001  merger->Resize(cnt);
4002  double value = 0.0;
4003  for(enumerator i = 0; i < Rows(); ++i)
4004  for(enumerator j = 0; j < Cols(); ++j)
4005  {
4006  value += (*this)(i,j).GetValue()*other(i,j).GetValue();
4007  merger->AddRow(other(i,j).GetValue(),(*this)(i,j).GetRow());
4008  merger->AddRow((*this)(i,j).GetValue(),other(i,j).GetRow());
4009  }
4010  ret.SetValue(value);
4011  merger->RetrieveRow(ret.GetRow());
4012  merger->Clear();
4013  }
4014  else*/
4015  {
4016  for(enumerator i = 0; i < Rows(); ++i)
4017  for(enumerator j = 0; j < Cols(); ++j)
4018  ret += compute(i,j)*other.compute(i,j);
4019  }
4020  return ret;
4021  }
4022 #endif //USE_AUTODIFF
4023 #endif
4024  template<typename Var>
4025  template<typename typeB>
4026  AbstractMatrix<Var> &
4027  AbstractMatrix<Var>::operator*=(const AbstractMatrixReadOnly<typeB> & B)
4028  {
4029  (*this) = (*this)*B;
4030  return *this;
4031  }
4032 
4033  template<typename Var>
4034  template<typename typeB>
4035  MatrixMulCoef<Var, typeB, typename Promote<Var, typeB>::type>
4036  AbstractMatrixReadOnly<Var>::operator*(const typeB& coef) const
4037  {
4038  return MatrixMulCoef<Var, typeB, typename Promote<Var, typeB>::type>(*this, coef);
4039  }
4040 
4041 #if defined(USE_AUTODIFF)
4042  template<typename Var>
4043  template<typename A>
4044  MatrixMulShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>
4045  AbstractMatrixReadOnly<Var>::operator*(shell_expression<A> const & coef) const
4046  {
4047  return MatrixMulShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>(*this, coef);
4048  }
4049 #endif //USE_AUTODIFF
4050 
4051  template<typename Var>
4052  template<typename typeB>
4053  AbstractMatrix<Var> &
4054  AbstractMatrix<Var>::operator*=(typeB coef)
4055  {
4056  for(enumerator i = 0; i < Rows(); ++i)
4057  for(enumerator j = 0; j < Cols(); ++j)
4058  assign((*this)(i,j),(*this)(i,j)*coef);
4059  return *this;
4060  }
4061 #if defined(USE_AUTODIFF)
4062  template<typename Var>
4063  template<typename A>
4064  MatrixDivShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>
4065  AbstractMatrixReadOnly<Var>::operator/(shell_expression<A> const & coef) const
4066  {
4067  return MatrixDivShellCoef<Var, shell_expression<A>, typename Promote<Var, variable>::type>(*this, coef);
4068  }
4069 #endif //USE_AUTODIFF
4070 
4071  template<typename Var>
4072  template<typename typeB>
4073  MatrixDivCoef<Var, typeB, typename Promote<Var, typeB>::type>
4074  AbstractMatrixReadOnly<Var>::operator/(const typeB& coef) const
4075  {
4076  return MatrixDivCoef<Var, typeB, typename Promote<Var, typeB>::type>(*this, coef);
4077  }
4078 
4079  template<typename Var>
4080  template<typename typeB>
4081  AbstractMatrix<Var> &
4082  AbstractMatrix<Var>::operator/=(typeB coef)
4083  {
4084  for(enumerator i = 0; i < Rows(); ++i)
4085  for(enumerator j = 0; j < Cols(); ++j)
4086  assign((*this)(i,j),(*this)(i,j)/coef);
4087  return *this;
4088  }
4089 
4090  template<typename Var>
4091  template<typename typeB>
4092  KroneckerProduct<Var,typeB>
4093  AbstractMatrixReadOnly<Var>::Kronecker(const AbstractMatrixReadOnly<typeB> & other) const
4094  {
4095  return KroneckerProduct<Var, typeB>(*this, other);
4096  }
4097 
4098  template<typename Var>
4099  Matrix<Var>
4100  AbstractMatrixReadOnly<Var>::Invert(int * ierr) const
4101  {
4102  Matrix<Var> ret(Cols(),Rows());
4103  ret = Solve(MatrixUnit<Var>(Rows()),ierr);
4104  return ret;
4105  }
4106 
4107  template<typename Var>
4108  Matrix<Var>
4109  AbstractMatrixReadOnly<Var>::CholeskyInvert(int * ierr) const
4110  {
4111  Matrix<Var> ret(Rows(),Rows());
4112  ret = CholeskySolve(MatrixUnit<Var>(Rows()),ierr);
4113  return ret;
4114  }
4115 
4116 
4117  template<typename Var>
4118  template<typename typeB>
4119  Matrix<typename Promote<Var,typeB>::type>
4120  AbstractMatrixReadOnly<Var>::CholeskySolve(const AbstractMatrixReadOnly<typeB> & B, int * ierr) const
4121  {
4122  const AbstractMatrixReadOnly<Var> & A = *this;
4123  assert(A.Rows() == A.Cols());
4124  assert(A.Rows() == B.Rows());
4125  enumerator n = A.Rows();
4126  enumerator l = B.Cols();
4127  Matrix<typename Promote<Var,typeB>::type> ret(B);
4128  static thread_private< SymmetricMatrix<Var> > L;// (A);
4129  *L = A;
4130 
4131  //SAXPY
4132  /*
4133  for(enumerator i = 0; i < n; ++i)
4134  {
4135  for(enumerator k = 0; k < i; ++k)
4136  for(enumerator j = i; j < n; ++j)
4137  L(i,j) -= L(i,k)*L(j,k);
4138  if( L(i,i) < 0.0 )
4139  {
4140  ret.second = false;
4141  if( print_fail ) std::cout << "Negative diagonal pivot " << get_value(L(i,i)) << std::endl;
4142  return ret;
4143  }
4144 
4145  L(i,i) = sqrt(L(i,i));
4146 
4147  if( fabs(L(i,i)) < 1.0e-24 )
4148  {
4149  ret.second = false;
4150  if( print_fail ) std::cout << "Diagonal pivot is too small " << get_value(L(i,i)) << std::endl;
4151  return ret;
4152  }
4153 
4154  for(enumerator j = i+1; j < n; ++j)
4155  L(i,j) = L(i,j)/L(i,i);
4156  }
4157  */
4158  //Outer product
4159  for(enumerator k = 0; k < n; ++k)
4160  {
4161  if( real_part((*L)(k,k)) < 0.0 )
4162  {
4163  if( ierr )
4164  {
4165  if( *ierr == -1 ) std::cout << "Negative diagonal pivot " << get_value((*L)(k,k)) << " row " << k << std::endl;
4166  *ierr = k+1;
4167  }
4168  else throw MatrixCholeskySolveFail;
4169  return ret;
4170  }
4171 
4172  (*L)(k,k) = sqrt((*L)(k,k));
4173 
4174  if( fabs((*L)(k,k)) < 1.0e-24 )
4175  {
4176  if( ierr )
4177  {
4178  if( *ierr == -1 ) std::cout << "Diagonal pivot is too small " << get_value((*L)(k,k)) << " row " << k << std::endl;
4179  *ierr = k+1;
4180  }
4181  else throw MatrixCholeskySolveFail;
4182  return ret;
4183  }
4184 
4185  for(enumerator i = k+1; i < n; ++i)
4186  (*L)(i,k) /= (*L)(k,k);
4187 
4188  for(enumerator j = k+1; j < n; ++j)
4189  {
4190  for(enumerator i = j; i < n; ++i)
4191  (*L)(i,j) -= (*L)(i,k)*(*L)(j,k);
4192  }
4193  }
4194  // LY=B
4195  Matrix<typename Promote<Var,typeB>::type> & Y = ret;
4196  for(enumerator i = 0; i < n; ++i)
4197  {
4198  for(enumerator k = 0; k < l; ++k)
4199  {
4200  for(enumerator j = 0; j < i; ++j)
4201  Y(i,k) -= Y(j,k)*(*L)(j,i);
4202  Y(i,k) /= (*L)(i,i);
4203  }
4204  }
4205  // L^TX = Y
4206  Matrix<typename Promote<Var,typeB>::type> & X = ret;
4207  for(enumerator it = n; it > 0; --it)
4208  {
4209  enumerator i = it-1;
4210  for(enumerator k = 0; k < l; ++k)
4211  {
4212  for(enumerator jt = n; jt > it; --jt)
4213  {
4214  enumerator j = jt-1;
4215  X(i,k) -= X(j,k)*(*L)(i,j);
4216  }
4217  X(i,k) /= (*L)(i,i);
4218  }
4219  }
4220  if( ierr ) *ierr = 0;
4221  return ret;
4222  }
4223 #if defined(USE_AUTODIFF)
4224  template<>
4225  template<>
4226  __INLINE Matrix<Promote<variable,variable>::type>
4227  AbstractMatrixReadOnly<variable>::CholeskySolve(const AbstractMatrixReadOnly<variable> & B, int * ierr) const
4228  {
4229  const AbstractMatrixReadOnly<variable> & A = *this;
4230  assert(A.Rows() == A.Cols());
4231  assert(A.Rows() == B.Rows());
4232  enumerator n = A.Rows();
4233  enumerator l = B.Cols();
4234  Matrix<Promote<variable,variable>::type> ret(B);
4235  SymmetricMatrix<variable> L(A);
4236  INMOST_DATA_ENUM_TYPE cnt = 0;
4237  cnt += A.GetMatrixCount();
4238  cnt += B.GetMatrixCount();
4239  Sparse::RowMerger* pmerger = NULL;
4240  if (cnt >= CNT_USE_MERGER)
4241  {
4242  merger->Resize(cnt);
4243  pmerger = &(*merger);
4244  }
4245  //Outer product
4246  for(enumerator k = 0; k < n; ++k)
4247  {
4248  if( L(k,k).GetValue() < 0.0 )
4249  {
4250  if( ierr )
4251  {
4252  if( *ierr == -1 ) std::cout << "Negative diagonal pivot " << get_value(L(k,k)) << " row " << k << std::endl;
4253  *ierr = k+1;
4254  }
4255  else throw MatrixCholeskySolveFail;
4256  return ret;
4257  }
4258 
4259  L(k,k) = sqrt(L(k,k));
4260 
4261  if( fabs(L(k,k).GetValue()) < 1.0e-24 )
4262  {
4263  if( ierr )
4264  {
4265  if( *ierr == -1 ) std::cout << "Diagonal pivot is too small " << get_value(L(k,k)) << " row " << k << std::endl;
4266  *ierr = k+1;
4267  }
4268  else throw MatrixCholeskySolveFail;
4269  return ret;
4270  }
4271 
4272  for(enumerator i = k+1; i < n; ++i)
4273  L(i,k) /= L(k,k);
4274 
4275  for(enumerator j = k+1; j < n; ++j)
4276  {
4277  for(enumerator i = j; i < n; ++i)
4278  L(i,j) -= L(i,k)*L(j,k);
4279  }
4280  }
4281  // LY=B
4282  Matrix<Promote<variable,variable>::type> & Y = ret;
4283  for(enumerator i = 0; i < n; ++i)
4284  {
4285  for(enumerator k = 0; k < l; ++k)
4286  {
4287  if( pmerger )
4288  {
4289  INMOST_DATA_REAL_TYPE value = Y(i,k).GetValue();
4290  pmerger->PushRow(1.0,Y(i,k).GetRow());
4291  for(enumerator j = 0; j < i; ++j)
4292  {
4293  value -= Y(j,k).GetValue()*L(j,i).GetValue();
4294  pmerger->AddRow(-Y(j,k).GetValue(),L(j,i).GetRow());
4295  pmerger->AddRow(-L(j,i).GetValue(),Y(j,k).GetRow());
4296  }
4297  Y(i,k).SetValue(value);
4298  pmerger->RetrieveRow(Y(i,k).GetRow());
4299  pmerger->Clear();
4300  }
4301  else
4302  {
4303  for(enumerator j = 0; j < i; ++j)
4304  Y(i,k) -= Y(j,k)*L(j,i);
4305  }
4306  Y(i,k) /= L(i,i);
4307  }
4308  }
4309  // L^TX = Y
4310  Matrix<Promote<variable,variable>::type> & X = ret;
4311  for(enumerator it = n; it > 0; --it)
4312  {
4313  enumerator i = it-1;
4314  for(enumerator k = 0; k < l; ++k)
4315  {
4316  if( pmerger )
4317  {
4318  double value = X(i,k).GetValue();
4319  pmerger->PushRow(1.0,X(i,k).GetRow());
4320  for(enumerator jt = n; jt > it; --jt)
4321  {
4322  enumerator j = jt-1;
4323  value -= X(j,k).GetValue()*L(i,j).GetValue();
4324  pmerger->AddRow(-X(j,k).GetValue(),L(i,j).GetRow());
4325  pmerger->AddRow(-L(i,j).GetValue(),X(j,k).GetRow());
4326  }
4327  X(i,k).SetValue(value);
4328  pmerger->RetrieveRow(X(i,k).GetRow());
4329  pmerger->Clear();
4330  }
4331  else
4332  {
4333  for(enumerator jt = n; jt > it; --jt)
4334  {
4335  enumerator j = jt-1;
4336  X(i,k) -= X(j,k)*L(i,j);
4337  }
4338  }
4339  X(i,k) /= L(i,i);
4340  }
4341  }
4342  if( ierr ) *ierr = 0;
4343  return ret;
4344  }
4345 
4346 
4347  template<>
4348  template<>
4349  __INLINE Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type>
4350  AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE>::CholeskySolve(const AbstractMatrixReadOnly<variable> & B, int * ierr) const
4351  {
4352  const AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE> & A = *this;
4353  assert(A.Rows() == A.Cols());
4354  assert(A.Rows() == B.Rows());
4355  enumerator n = A.Rows();
4356  enumerator l = B.Cols();
4357  Matrix<Promote<INMOST_DATA_REAL_TYPE,variable>::type> ret(B);
4358  SymmetricMatrix<INMOST_DATA_REAL_TYPE> L(A);
4359  INMOST_DATA_ENUM_TYPE cnt = B.GetMatrixCount();
4360  Sparse::RowMerger* pmerger = NULL;
4361  if (cnt >= CNT_USE_MERGER)
4362  {
4363  merger->Resize(cnt);
4364  pmerger = &(*merger);
4365  }
4366  //Outer product
4367  for(enumerator k = 0; k < n; ++k)
4368  {
4369  if( L(k,k) < 0.0 )
4370  {
4371  if( ierr )
4372  {
4373  if( *ierr == -1 ) std::cout << "Negative diagonal pivot " << get_value(L(k,k)) << " row " << k << std::endl;
4374  *ierr = k+1;
4375  }
4376  else throw MatrixCholeskySolveFail;
4377  return ret;
4378  }
4379 
4380  L(k,k) = sqrt(L(k,k));
4381 
4382  if( fabs(L(k,k)) < 1.0e-24 )
4383  {
4384  if( ierr )
4385  {
4386  if( *ierr == -1 ) std::cout << "Diagonal pivot is too small " << get_value(L(k,k)) << " row " << k << std::endl;
4387  *ierr = k+1;
4388  }
4389  else throw MatrixCholeskySolveFail;
4390  return ret;
4391  }
4392 
4393  for(enumerator i = k+1; i < n; ++i)
4394  L(i,k) /= L(k,k);
4395 
4396  for(enumerator j = k+1; j < n; ++j)
4397  {
4398  for(enumerator i = j; i < n; ++i)
4399  L(i,j) -= L(i,k)*L(j,k);
4400  }
4401  }
4402  // LY=B
4403  Matrix<Promote<variable,variable>::type> & Y = ret;
4404  for(enumerator i = 0; i < n; ++i)
4405  {
4406  for(enumerator k = 0; k < l; ++k)
4407  {
4408  if( pmerger )
4409  {
4410  double value = Y(i,k).GetValue();
4411  pmerger->PushRow(1.0,Y(i,k).GetRow());
4412  for(enumerator j = 0; j < i; ++j)
4413  {
4414  value -= Y(j,k).GetValue()*L(j,i);
4415  pmerger->AddRow(-L(j,i),Y(j,k).GetRow());
4416  }
4417  Y(i,k).SetValue(value);
4418  pmerger->RetrieveRow(Y(i,k).GetRow());
4419  pmerger->Clear();
4420  }
4421  else
4422  {
4423  for(enumerator j = 0; j < i; ++j)
4424  Y(i,k) -= Y(j,k)*L(j,i);
4425  }
4426  Y(i,k) /= L(i,i);
4427  }
4428  }
4429  // L^TX = Y
4430  Matrix<Promote<variable,variable>::type> & X = ret;
4431  for(enumerator it = n; it > 0; --it)
4432  {
4433  enumerator i = it-1;
4434  for(enumerator k = 0; k < l; ++k)
4435  {
4436  if( pmerger )
4437  {
4438  double value = X(i,k).GetValue();
4439  pmerger->PushRow(1.0,X(i,k).GetRow());
4440  for(enumerator jt = n; jt > it; --jt)
4441  {
4442  enumerator j = jt-1;
4443  value -= X(j,k).GetValue()*L(i,j);
4444  pmerger->AddRow(-L(i,j),X(j,k).GetRow());
4445  }
4446  X(i,k).SetValue(value);
4447  pmerger->RetrieveRow(X(i,k).GetRow());
4448  pmerger->Clear();
4449  }
4450  else
4451  {
4452  for(enumerator jt = n; jt > it; --jt)
4453  {
4454  enumerator j = jt-1;
4455  X(i,k) -= X(j,k)*L(i,j);
4456  }
4457  }
4458  X(i,k) /= L(i,i);
4459  }
4460  }
4461  if( ierr ) *ierr = 0;
4462  return ret;
4463  }
4464 
4465  template<>
4466  template<>
4467  __INLINE Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type>
4468  AbstractMatrixReadOnly<variable>::CholeskySolve(const AbstractMatrixReadOnly<INMOST_DATA_REAL_TYPE> & B, int * ierr) const
4469  {
4470  const AbstractMatrixReadOnly<variable> & A = *this;
4471  assert(A.Rows() == A.Cols());
4472  assert(A.Rows() == B.Rows());
4473  enumerator n = A.Rows();
4474  enumerator l = B.Cols();
4475  Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type> ret(B);
4476  SymmetricMatrix<variable> L(A);
4477  INMOST_DATA_ENUM_TYPE cnt = A.GetMatrixCount();
4478  Sparse::RowMerger* pmerger = NULL;
4479  if (cnt >= CNT_USE_MERGER)
4480  {
4481  merger->Resize(cnt);
4482  pmerger = &(*merger);
4483  }
4484  //Outer product
4485  for(enumerator k = 0; k < n; ++k)
4486  {
4487  if( L(k,k).GetValue() < 0.0 )
4488  {
4489  if( ierr )
4490  {
4491  if( *ierr == -1 ) std::cout << "Negative diagonal pivot " << get_value(L(k,k)) << " row " << k << std::endl;
4492  *ierr = k+1;
4493  }
4494  else throw MatrixCholeskySolveFail;
4495  return ret;
4496  }
4497 
4498  L(k,k) = sqrt(L(k,k));
4499 
4500  if( fabs(L(k,k).GetValue()) < 1.0e-24 )
4501  {
4502  if( ierr )
4503  {
4504  if( *ierr == -1 ) std::cout << "Diagonal pivot is too small " << get_value(L(k,k)) << " row " << k << std::endl;
4505  *ierr = k+1;
4506  }
4507  else throw MatrixCholeskySolveFail;
4508  return ret;
4509  }
4510 
4511  for(enumerator i = k+1; i < n; ++i)
4512  L(i,k) /= L(k,k);
4513 
4514  for(enumerator j = k+1; j < n; ++j)
4515  {
4516  for(enumerator i = j; i < n; ++i)
4517  L(i,j) -= L(i,k)*L(j,k);
4518  }
4519  }
4520  // LY=B
4521  Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type> & Y = ret;
4522  for(enumerator i = 0; i < n; ++i)
4523  {
4524  for(enumerator k = 0; k < l; ++k)
4525  {
4526  if( pmerger )
4527  {
4528  double value = Y(i,k).GetValue();
4529  pmerger->PushRow(1.0,Y(i,k).GetRow());
4530  for(enumerator j = 0; j < i; ++j)
4531  {
4532  value -= Y(j,k).GetValue()*L(j,i).GetValue();
4533  pmerger->AddRow(-Y(j,k).GetValue(),L(j,i).GetRow());
4534  pmerger->AddRow(-L(j,i).GetValue(),Y(j,k).GetRow());
4535  }
4536  Y(i,k).SetValue(value);
4537  pmerger->RetrieveRow(Y(i,k).GetRow());
4538  pmerger->Clear();
4539  }
4540  else
4541  {
4542  for(enumerator j = 0; j < i; ++j)
4543  Y(i,k) -= Y(j,k)*L(j,i);
4544  }
4545  Y(i,k) /= L(i,i);
4546  }
4547  }
4548  // L^TX = Y
4549  Matrix<Promote<variable,INMOST_DATA_REAL_TYPE>::type> & X = ret;
4550  for(enumerator it = n; it > 0; --it)
4551  {
4552  enumerator i = it-1;
4553  for(enumerator k = 0; k < l; ++k)
4554  {
4555  if( pmerger )
4556  {
4557  double value = X(i,k).GetValue();
4558  pmerger->PushRow(1.0,X(i,k).GetRow());
4559  for(enumerator jt = n; jt > it; --jt)
4560  {
4561  enumerator j = jt-1;
4562  value -= X(j,k).GetValue()*L(i,j).GetValue();
4563  pmerger->AddRow(-X(j,k).GetValue(),L(i,j).GetRow());
4564  pmerger->AddRow(-L(i,j).GetValue(),X(j,k).GetRow());
4565  }
4566  X(i,k).SetValue(value);
4567  pmerger->RetrieveRow(X(i,k).GetRow());
4568  pmerger->Clear();
4569  }
4570  else
4571  {
4572  for(enumerator jt = n; jt > it; --jt)
4573  {
4574  enumerator j = jt-1;
4575  X(i,k) -= X(j,k)*L(i,j);
4576  }
4577  }
4578  X(i,k) /= L(i,i);
4579  }
4580  }
4581  if( ierr ) *ierr = 0;
4582  return ret;
4583  }
4584 #endif //USE_AUTODIFF
4585  template<typename Var>
4586  template<typename typeB>
4587  Matrix<typename Promote<Var,typeB>::type>
4588  AbstractMatrixReadOnly<Var>::Solve(const AbstractMatrixReadOnly<typeB> & B, int * ierr) const
4589  {
4590  // for A^T B
4591  assert(Rows() == B.Rows());
4592 
4593  if( Rows() != Cols() )
4594  {
4595  ConstMatrixTranspose<Var> At = this->Transpose(); //m by n matrix
4596  return (At * (*this)).Solve(At * B, ierr);
4597  }
4598  Matrix<typeB> AtB = B; //m by l matrix
4599  Matrix<Var> AtA = (*this); //m by m matrix
4600  enumerator l = B.Cols();
4601  enumerator m = Cols();
4602  Matrix<typename Promote<Var,typeB>::type> ret(m,l);
4603  //ConstMatrixTranspose<Var> At = this->Transpose(); //m by n matrix
4604  //Matrix<typename Promote<Var,typeB>::type> AtB = At*B; //m by l matrix
4605  //Matrix<Var> AtA = At*(*this); //m by m matrix
4606  assert(l == AtB.Cols());
4607  //enumerator l = AtB.Cols();
4608  //enumerator n = Rows();
4609 
4610  std::vector<enumerator> order(m);
4611 
4612  Var temp;
4613  INMOST_DATA_REAL_TYPE max,v;
4614  typename Promote<Var,typeB>::type tempb;
4615  for(enumerator i = 0; i < m; ++i) order[i] = i;
4616  for(enumerator i = 0; i < m; i++)
4617  {
4618  enumerator maxk = i, maxq = i, temp2;
4619  max = fabs(get_value(AtA(maxk,maxq)));
4620  //Find best pivot
4621  //if( max < 1.0e-8 )
4622  {
4623  for(enumerator k = i; k < m; k++) // over rows
4624  {
4625  for(enumerator q = i; q < m; q++) // over columns
4626  {
4627  v = fabs(get_value(AtA(k,q)));
4628  if( v > max )
4629  {
4630  max = v;
4631  maxk = k;
4632  maxq = q;
4633  }
4634  }
4635  }
4636  //Exchange rows
4637  if( maxk != i )
4638  {
4639  for(enumerator q = 0; q < m; q++) // over columns of A
4640  {
4641  //std::swap(AtA(maxk,q),AtA(i,q));
4642  temp = AtA(maxk,q);
4643  AtA(maxk,q) = AtA(i,q);
4644  AtA(i,q) = temp;
4645  }
4646  //exchange rhs
4647  for(enumerator q = 0; q < l; q++) // over columns of B
4648  {
4649  //std::swap(AtB(maxk,q),AtB(i,q));
4650  tempb = AtB(maxk,q);
4651  AtB(maxk,q) = AtB(i,q);
4652  AtB(i,q) = tempb;
4653  }
4654  }
4655  //Exchange columns
4656  if( maxq != i )
4657  {
4658  for(enumerator k = 0; k < m; k++) //over rows
4659  {
4660  //std::swap(AtA(k,maxq),AtA(k,i));
4661  temp = AtA(k,maxq);
4662  AtA(k,maxq) = AtA(k,i);
4663  AtA(k,i) = temp;
4664  }
4665  //remember order in sol
4666  {
4667  //std::swap(order[maxq],order[i]);
4668  temp2 = order[maxq];
4669  order[maxq] = order[i];
4670  order[i] = temp2;
4671  }
4672  }
4673  }
4674  //Check small entry
4675  if( fabs(get_value(AtA(i,i))) < 1.0e-54 )
4676  {
4677  bool ok = true;
4678  for(enumerator k = 0; k < l; k++) // over columns of B
4679  {
4680  if( fabs(get_value(AtB(i,k))/1.0e-54) > 1 )
4681  {
4682  ok = false;
4683  break;
4684  }
4685  }
4686  if( ok ) AtA(i,i) = real_part(AtA(i,i)) < 0.0 ? - 1.0e-12 : 1.0e-12;
4687  else
4688  {
4689  if( ierr )
4690  {
4691  if( *ierr == -1 )
4692  {
4693  std::cout << "Failed to invert matrix diag " << get_value(AtA(i,i)) << std::endl;
4694  std::cout << "rhs:";
4695  for(enumerator k = 0; k < l; k++)
4696  std::cout << " " << get_value(AtB(i,k));
4697  std::cout << std::endl;
4698  }
4699  *ierr = i+1;
4700  }
4701  else throw MatrixSolveFail;
4702  return ret;
4703  }
4704  }
4705  //Divide row and column by diagonal values
4706  for(enumerator k = i+1; k < m; k++)
4707  {
4708  AtA(i,k) /= AtA(i,i);
4709  AtA(k,i) /= AtA(i,i);
4710  }
4711  //Elimination step for matrix
4712  for(enumerator k = i+1; k < m; k++)
4713  for(enumerator q = i+1; q < m; q++)
4714  {
4715  AtA(k,q) -= AtA(k,i) * AtA(i,i) * AtA(i,q);
4716  }
4717  //Elimination step for right hand side
4718  for(enumerator k = 0; k < l; k++)
4719  {
4720  for(enumerator j = i+1; j < m; j++) //iterate over columns of L
4721  {
4722  AtB(j,k) -= AtB(i,k) * AtA(j,i);
4723  }
4724  AtB(i,k) /= AtA(i,i);
4725  }
4726  }
4727  //Back substitution step
4728  for(enumerator k = 0; k < l; k++)
4729  {
4730  for(enumerator i = m; i-- > 0; ) //iterate over rows of U
4731  for(enumerator j = i+1; j < m; j++)
4732  {
4733  AtB(i,k) -= AtB(j,k) * AtA(i,j);
4734  }
4735  for(enumerator i = 0; i < m; i++)
4736  ret(order[i],k) = AtB(i,k);
4737  }
4738  //delete [] order;
4739  if( ierr ) *ierr = 0;
4740  return ret;
4741  }
4742 
4743  template<typename Var>
4744  Matrix<Var>
4745  AbstractMatrixReadOnly<Var>::ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const
4746  {
4747  assert(ibeg < Rows());
4748  assert(iend < Rows());
4749  assert(jbeg < Cols());
4750  assert(jend < Cols());
4751  Matrix<Var> ret(iend-ibeg,jend-jbeg);
4752  for(enumerator i = ibeg; i < iend; ++i)
4753  {
4754  for(enumerator j = jbeg; j < jend; ++j)
4755  ret(i-ibeg,j-jbeg) = (*this)(i,j);
4756  }
4757  return ret;
4758  }
4759 
4760  template<typename Var>
4761  Matrix<Var>
4762  AbstractMatrixReadOnly<Var>::PseudoInvert(INMOST_DATA_REAL_TYPE tol, int * ierr) const
4763  {
4764  Matrix<Var> ret(Cols(),Rows());
4765  Matrix<Var> U,S,V;
4766  bool success = SVD(U,S,V);
4767  if( !success )
4768  {
4769  if( ierr )
4770  {
4771  if( *ierr == -1 ) std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl;
4772  *ierr = 1;
4773  return ret;
4774  }
4775  else throw MatrixPseudoSolveFail;
4776  }
4777  for(INMOST_DATA_ENUM_TYPE k = 0; k < S.Cols(); ++k)
4778  {
4779  if (get_value(fabs(S(k, k))) > tol)
4780  S(k, k) = 1.0 / S(k, k);
4781  else
4782  S(k, k) = 0.0;
4783  }
4784  //ret = V*S*U.Transpose();
4785  ret = V.ConjugateTranspose().Transpose() * S * U.ConjugateTranspose();
4786  if( ierr ) *ierr = 0;
4787  return ret;
4788  }
4789 
4790  template<typename Var>
4791  Matrix<Var>
4792  AbstractMatrixReadOnly<Var>::Root(INMOST_DATA_ENUM_TYPE iter, INMOST_DATA_REAL_TYPE tol, int *ierr) const
4793  {
4794  assert(Rows() == Cols());
4795  Matrix<Var> ret(Cols(),Cols());
4796  Matrix<Var> ret0(Cols(),Cols());
4797  ret.Zero();
4798  ret0.Zero();
4799  int k = 0;
4800  for(k = 0; k < Cols(); ++k) ret(k,k) = ret0(k,k) = 1;
4801  while(k < iter)
4802  {
4803  ret0 = ret;
4804  ret = 0.5*(ret + (*this)*ret.Invert());
4805  if( (ret - ret0).FrobeniusNorm() < tol ) return ret;
4806  }
4807  if( ierr )
4808  {
4809  if( *ierr == -1 ) std::cout << "Failed to find square root of matrix by Babylonian method" << std::endl;
4810  *ierr = 1;
4811  return ret;
4812  }
4813  return ret;
4814  }
4815  template<typename Var>
4816  Matrix<Var>
4817  AbstractMatrixReadOnly<Var>::Power(INMOST_DATA_REAL_TYPE n, int * ierr) const
4818  {
4819  Matrix<Var> ret(Cols(),Rows());
4820  Matrix<Var> U, S, V;
4821  bool success = SVD(U,S,V);
4822  if( !success )
4823  {
4824  if( ierr )
4825  {
4826  if( *ierr == -1 )
4827  std::cout << "Failed to compute singular value decomposition of the matrix" << std::endl;
4828  *ierr = 1;
4829  return ret;
4830  }
4831  else throw MatrixPseudoSolveFail;
4832  }
4833  for(INMOST_DATA_ENUM_TYPE k = 0; k < S.Cols(); ++k)
4834  if( get_value(S(k,k)) ) S(k,k) = pow(S(k,k),n);
4835  if( n >= 0 )
4836  ret = U*S*V.Transpose();
4837  else
4838  ret = V*S*U.Transpose();
4839  if( ierr ) *ierr = 0;
4840  return ret;
4841  }
4842  template<typename Var>
4843  template<typename typeB>
4844  Matrix<typename Promote<Var,typeB>::type>
4845  AbstractMatrixReadOnly<Var>::PseudoSolve(const AbstractMatrixReadOnly<typeB> & B, INMOST_DATA_REAL_TYPE tol, int * ierr) const
4846  {
4847  Matrix<typename Promote<Var,typeB>::type> ret(Cols(),B.Cols());
4848  Matrix<Var> U,S,V;
4849  bool success = SVD(U,S,V);
4850  if( !success )
4851  {
4852  if( ierr )
4853  {
4854  if( *ierr == -1 ) std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl;
4855  *ierr = 1;
4856  }
4857  else throw MatrixPseudoSolveFail;
4858  }
4859  for(int k = 0; k < (int)S.Cols(); ++k)
4860  {
4861  if( S(k,k) > tol )
4862  S(k,k) = 1.0/S(k,k);
4863  else
4864  S(k,k) = 0.0;
4865  }
4866  ret = V*S*U.Transpose()*B;
4867  if( ierr ) *ierr = 0;
4868  return ret;
4869  }
4870 
4871  template<typename Var>
4872  SubMatrix<Var> AbstractMatrix<Var>::operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col)
4873  {
4874  return SubMatrix<Var>(*this,first_row,last_row,first_col,last_col);
4875  }
4876  template<typename Var>
4877  ConstSubMatrix<Var> AbstractMatrixReadOnly<Var>::operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const
4878  {
4879  return ConstSubMatrix<Var>(*this, first_row, last_row, first_col, last_col);
4880  }
4881 
4882  template<typename Var>
4883  BlockOfMatrix<Var> AbstractMatrix<Var>::BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col)
4884  {
4885  return BlockOfMatrix<Var>(*this,nrows,ncols,offset_row,offset_col);
4886  }
4887  template<typename Var>
4888  ConstBlockOfMatrix<Var> AbstractMatrixReadOnly<Var>::BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col) const
4889  {
4890  return ConstBlockOfMatrix<Var>(*this,nrows,ncols,offset_row,offset_col);
4891  }
4892 
4893 
4894  template<class T> struct make_integer;
4895  template<> struct make_integer<float> {typedef int type;};
4896  template<> struct make_integer<double> {typedef long long type;};
4897 
4898  __INLINE static bool compare(INMOST_DATA_REAL_TYPE * a, INMOST_DATA_REAL_TYPE * b)
4899  {
4900  return (*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(a)) <=
4901  (*reinterpret_cast< make_integer<INMOST_DATA_REAL_TYPE>::type * >(b));
4902  }
4903 
4904  class BinaryHeapDense
4905  {
4906  INMOST_DATA_ENUM_TYPE size_max, size;
4907  std::vector<INMOST_DATA_ENUM_TYPE> heap;
4908  std::vector<INMOST_DATA_ENUM_TYPE> index;
4909  INMOST_DATA_REAL_TYPE * keys;
4910 
4911  void swap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_ENUM_TYPE j)
4912  {
4913  INMOST_DATA_ENUM_TYPE t = heap[i];
4914  heap[i] = heap[j];
4915  heap[j] = t;
4916  index[heap[i]] = i;
4917  index[heap[j]] = j;
4918  }
4919  void bubbleUp(INMOST_DATA_ENUM_TYPE k)
4920  {
4921  while(k > 1 && keys[heap[k/2]] > keys[heap[k]])
4922  {
4923  swap(k, k/2);
4924  k = k/2;
4925  }
4926  }
4927  void bubbleDown(INMOST_DATA_ENUM_TYPE k)
4928  {
4929  size_t j;
4930  while(2*k <= size)
4931  {
4932  j = 2*static_cast<size_t>(k);
4933  if(j < size && keys[heap[j]] > keys[heap[j+1]])
4934  j++;
4935  if(keys[heap[k]] <= keys[heap[j]])
4936  break;
4937  swap(k, static_cast<INMOST_DATA_ENUM_TYPE>(j));
4938  k = static_cast<INMOST_DATA_ENUM_TYPE>(j);
4939  }
4940  }
4941  public:
4942  BinaryHeapDense(INMOST_DATA_REAL_TYPE * pkeys, INMOST_DATA_ENUM_TYPE len)
4943  {
4944  size_max = len;
4945  keys = pkeys;
4946  size = 0;
4947  heap.resize(static_cast<size_t>(size_max)+1);
4948  index.resize(static_cast<size_t>(size_max)+1);
4949  for(INMOST_DATA_ENUM_TYPE i = 0; i <= size_max; i++)
4950  index[i] = ENUMUNDEF;
4951  }
4952  void PushHeap(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
4953  {
4954  size++;
4955  index[i] = size;
4956  heap[size] = i;
4957  keys[i] = key;
4958  bubbleUp(size);
4959  }
4960  INMOST_DATA_ENUM_TYPE PopHeap()
4961  {
4962  if( size == 0 ) return ENUMUNDEF;
4963  INMOST_DATA_ENUM_TYPE min = heap[1];
4964  swap(1, size--);
4965  bubbleDown(1);
4966  index[min] = ENUMUNDEF;
4967  heap[static_cast<size_t>(size)+1] = ENUMUNDEF;
4968  return min;
4969  }
4970  void DecreaseKey(INMOST_DATA_ENUM_TYPE i, INMOST_DATA_REAL_TYPE key)
4971  {
4972  keys[i] = key;
4973  bubbleUp(index[i]);
4974  }
4975  void Clear()
4976  {
4977  while( size ) keys[PopHeap()] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
4978  }
4979  bool Contains(INMOST_DATA_ENUM_TYPE i)
4980  {
4981  return index[i] != ENUMUNDEF;
4982  }
4983  };
4984 
4985  template<typename Var>
4986  void AbstractMatrix<Var>::MPT(INMOST_DATA_ENUM_TYPE * Perm, INMOST_DATA_REAL_TYPE * SL, INMOST_DATA_REAL_TYPE * SR) const
4987  {
4988  const INMOST_DATA_ENUM_TYPE EOL = ENUMUNDEF-1;
4989  int n = Rows();
4990  int m = Cols();
4991  INMOST_DATA_REAL_TYPE u, l;
4992  std::vector<INMOST_DATA_REAL_TYPE> Cmax(m,0.0);
4993  std::vector<INMOST_DATA_REAL_TYPE> U(m,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
4994  std::vector<INMOST_DATA_REAL_TYPE> V(n,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
4995  std::fill(Perm,Perm+m,ENUMUNDEF);
4996  std::vector<INMOST_DATA_ENUM_TYPE> IPerm(std::max(n,m),ENUMUNDEF);
4997  std::vector<INMOST_DATA_ENUM_TYPE> ColumnList(m,ENUMUNDEF);
4998  std::vector<INMOST_DATA_ENUM_TYPE> ColumnPosition(n,ENUMUNDEF);
4999  std::vector<INMOST_DATA_ENUM_TYPE> Parent(n,ENUMUNDEF);
5000  std::vector<INMOST_DATA_REAL_TYPE> Dist(m,std::numeric_limits<INMOST_DATA_REAL_TYPE>::max());
5001  const AbstractMatrix<Var> & A = *this;
5002  Matrix<INMOST_DATA_REAL_TYPE> C(n,m);
5003  INMOST_DATA_ENUM_TYPE Li, Ui;
5004  INMOST_DATA_ENUM_TYPE ColumnBegin, PathEnd, Trace, IPermPrev;
5005  INMOST_DATA_REAL_TYPE ShortestPath, AugmentPath;
5006  BinaryHeapDense Heap(&Dist[0],m);
5007  //Initial LOG transformation to dual problem and initial extreme match
5008  for(int k = 0; k < n; ++k)
5009  {
5010  for(int i = 0; i < m; ++i)
5011  {
5012  C(k,i) = fabs(get_value(A(k,i)));
5013  if( Cmax[i] < C(k,i) ) Cmax[i] = C(k,i);
5014  }
5015  }
5016  for(int k = 0; k < n; ++k)
5017  {
5018  for (int i = 0; i < m; ++i)
5019  {
5020  if( Cmax[i] == 0 || C(k,i) == 0 )
5021  C(k,i) = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
5022  else
5023  {
5024  C(k,i) = log(Cmax[i])-log(C(k,i));
5025  if( C(k,i) < U[i] ) U[i] = C(k,i);
5026  }
5027  }
5028  }
5029  for(int k = 0; k < n; ++k)
5030  {
5031  for (int i = 0; i < m; ++i)
5032  {
5033  u = C(k,i) - U[i];
5034  if( u < V[k] ) V[k] = u;
5035  }
5036  }
5037  // Update cost and match
5038  for(int k = 0; k < n; ++k)
5039  {
5040  for (int i = 0; i < m; ++i)
5041  {
5042  u = fabs(C(k,i) - V[k] - U[i]);
5043  if( u < 1.0e-30 && Perm[i] == ENUMUNDEF && IPerm[k] == ENUMUNDEF )
5044  {
5045  Perm[i] = k;
5046  IPerm[k] = i;
5047  ColumnPosition[k] = i;
5048  }
5049  }
5050  }
5051  //1-step augmentation
5052  for(int k = 0; k < n; ++k)
5053  {
5054  if( IPerm[k] == ENUMUNDEF ) //unmatched row
5055  {
5056  for (int i = 0; i < m && IPerm[k] == ENUMUNDEF; ++i)
5057  {
5058  u = fabs(C(k,i) - V[k] - U[i]);
5059  if( u <= 1.0e-30 )
5060  {
5061  Li = Perm[i];
5062  assert(Li != ENUMUNDEF);
5063  // Search other row in C for 0
5064  for (int Lit = 0; Lit < m; ++Lit)
5065  {
5066  u = fabs(C(Li,Lit) - V[Li] - U[Lit]);
5067  if( u <= 1.0e-30 && Perm[Lit] == ENUMUNDEF )
5068  {
5069  Perm[i] = k;
5070  IPerm[k] = i;
5071  ColumnPosition[k] = i;
5072  Perm[Lit] = Li;
5073  IPerm[Li] = Lit;
5074  ColumnPosition[Li] = Lit;
5075  break;
5076  }
5077  }
5078  }
5079  }
5080  }
5081  }
5082  // Weighted bipartite matching
5083  for(int k = 0; k < n; ++k)
5084  {
5085  if( IPerm[k] != ENUMUNDEF )
5086  continue;
5087  Li = k;
5088  ColumnBegin = EOL;
5089  Parent[Li] = ENUMUNDEF;
5090  PathEnd = ENUMUNDEF;
5091  Trace = k;
5092  ShortestPath = 0;
5093  AugmentPath = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
5094  while(true)
5095  {
5096  for (int Lit = 0; Lit < m; ++Lit)
5097  {
5098  if( ColumnList[Lit] != ENUMUNDEF ) continue;
5099  l = ShortestPath + C(Li,Lit) - V[Li] - U[Lit];
5100  if( l < 0.0 && fabs(l) < 1.0e-12 ) l = 0;
5101  if( l < AugmentPath )
5102  {
5103  if( Perm[Lit] == ENUMUNDEF )
5104  {
5105  PathEnd = Lit;
5106  Trace = Li;
5107  AugmentPath = l;
5108  }
5109  else if( l < Dist[Lit] )
5110  {
5111  Parent[Perm[Lit]] = Li;
5112  if( Heap.Contains(Lit) )
5113  Heap.DecreaseKey(Lit,l);
5114  else
5115  Heap.PushHeap(Lit,l);
5116  }
5117  }
5118  }
5119 
5120  INMOST_DATA_ENUM_TYPE pop_heap_pos = Heap.PopHeap();
5121  if( pop_heap_pos == ENUMUNDEF ) break;
5122 
5123  Ui = pop_heap_pos;
5124  ShortestPath = Dist[Ui];
5125 
5126  if( AugmentPath <= ShortestPath )
5127  {
5128  Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
5129  break;
5130  }
5131 
5132 
5133  ColumnList[Ui] = ColumnBegin;
5134  ColumnBegin = Ui;
5135 
5136  Li = Perm[Ui];
5137 
5138  }
5139  if( PathEnd != ENUMUNDEF )
5140  {
5141  Ui = ColumnBegin;
5142  while(Ui != EOL)
5143  {
5144  U[Ui] += Dist[Ui] - AugmentPath;
5145  if( Perm[Ui] != ENUMUNDEF ) V[Perm[Ui]] = C(Perm[Ui],ColumnPosition[Perm[Ui]]) - U[Ui];
5146  Dist[Ui] = std::numeric_limits<INMOST_DATA_REAL_TYPE>::max();
5147  Li = ColumnList[Ui];
5148  ColumnList[Ui] = ENUMUNDEF;
5149  Ui = Li;
5150  }
5151 
5152  Ui = PathEnd;
5153  while(Trace != ENUMUNDEF)
5154  {
5155  IPermPrev = IPerm[Trace];
5156  Perm[Ui] = Trace;
5157  IPerm[Trace] = Ui;
5158 
5159  ColumnPosition[Trace] = Ui;
5160  V[Trace] = C(Trace,ColumnPosition[Trace]) - U[Ui];
5161 
5162  Ui = IPermPrev;
5163  Trace = Parent[Trace];
5164 
5165  }
5166  Heap.Clear();
5167  }
5168  }
5169 
5170  if( SL || SR )
5171  {
5172  for (int k = 0; k < std::min(n,m); ++k)
5173  {
5174  l = (V[k] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() ? 1 : exp(V[k]));
5175  u = (U[k] == std::numeric_limits<INMOST_DATA_REAL_TYPE>::max() || Cmax[k] == 0 ? 1 : exp(U[k])/Cmax[k]);
5176 
5177  if( l*get_value(A(k,Perm[k]))*u < 0.0 ) l *= -1;
5178 
5179  if( SL ) SL[Perm[k]] = u;
5180  if( SR ) SR[k] = l;
5181  }
5182  if( SR ) for(int k = std::min(n,m); k < n; ++k) SR[k] = 1;
5183  if( SL ) for(int k = std::min(n,m); k < m; ++k) SL[k] = 1;
5184  }
5185  //check that there are no gaps in Perm
5186  std::fill(IPerm.begin(),IPerm.end(),ENUMUNDEF);
5187  std::vector<INMOST_DATA_ENUM_TYPE> gaps;
5188  for(int k = 0; k < m; ++k)
5189  {
5190  if( Perm[k] != ENUMUNDEF )
5191  IPerm[Perm[k]] = 0;
5192  }
5193  for(int k = 0; k < m; ++k)
5194  {
5195  if( IPerm[k] == ENUMUNDEF )
5196  gaps.push_back(k);
5197  }
5198  for(int k = 0; k < m; ++k)
5199  {
5200  if( Perm[k] == ENUMUNDEF )
5201  {
5202  Perm[k] = gaps.back();
5203  gaps.pop_back();
5204  }
5205  }
5206  }
5207 
5208  template<typename Var>
5209  bool AbstractMatrixReadOnly<Var>::cSVD(AbstractMatrix<Var>& U, AbstractMatrix<Var>& Sigma, AbstractMatrix<Var>& V) const
5210  {
5211  int m = Rows();
5212  int n = Cols();
5213  //data eta / 1.1920929e-07 /
5214  //data tol / 1.5e-31 /
5215  if (m > n)
5216  {
5217  bool success = ConjugateTranspose().cSVD(V, Sigma, U);
5218  if (success)
5219  return true;
5220  else return false;
5221  }
5222  Var cs, eta, f, g, h, q, r, sn, w, x, y, z;
5223  int i, j, k, l, p = 0;
5224  std::vector<Var> t, b, c, s;
5225  Matrix<Var> A = *this;
5226  // Householder reduction.
5227  c.resize(n);
5228  t.resize(n);
5229  s.resize(n);
5230  b.resize(n);
5231  c[0] = 0.0;
5232  k = 0;
5233 
5234  while (true)
5235  {
5236  // Elimination of A(I, K), I = K + 1, ..., M.
5237  z = 0.0;
5238  for (i = k; i < m; ++i)
5239  z += A(i, k) * conj(A(i, k));
5240  b[k] = 0.0;
5241 
5242  if (fabs(get_value(z)))
5243  {
5244  z = sqrt(z);
5245  b[k] = z;
5246  w = fabs(A(k, k));
5247  q = 1.0;
5248  if (fabs(get_value(w))) q = A(k, k) / w;
5249  A(k, k) = q * (z + w);
5250 
5251  if (k != n + p - 1)
5252  {
5253  for (j = k + 1; j < n + p; ++j)
5254  {
5255  q = 0.0;
5256  for (i = k; i < m; ++i)
5257  q += conj(A(i, k)) * A(i, j);
5258  q = q / (z * (z + w));
5259  for (i = k; i < m; ++i)
5260  A(i, j) -= q * A(i, k);
5261  }
5262  // Phase transformation.
5263  q = -conj(A(k, k)) / fabs(A(k, k));
5264  for (j = k + 1; j < n + p; ++j)
5265  A(k, j) *= q;
5266  }
5267  }
5268  // Elimination of A(K, J), J = K + 2, ..., N
5269  if (k == n - 1) break;
5270  z = 0.0;
5271  for (j = k + 1; j < n; ++j)
5272  z += conj(A(k, j)) * A(k, j);
5273  c[k + 1] = 0.0;
5274 
5275  if (fabs(get_value(z)))
5276  {
5277  z = sqrt(z);
5278  c[k + 1] = z;
5279  w = fabs(A(k, k + 1));
5280  q = 1.0;
5281  if (fabs(get_value(w))) q = A(k, k + 1) / w;
5282  A(k, k + 1) = q * (z + w);
5283  for (i = k + 1; i < m; ++i)
5284  {
5285  q = 0.0;
5286  for (j = k + 1; j < n; ++j)
5287  q += conj(A(k, j)) * A(i, j);
5288  q = q / (z * (z + w));
5289 
5290  for (j = k + 1; j < n; ++j)
5291  A(i, j) -= q * A(k, j);
5292  }
5293  // Phase transformation.
5294  q = -conj(A(k, k + 1)) / fabs(A(k, k + 1));
5295  for (i = k + 1; i < m; ++i)
5296  A(i, k + 1) = A(i, k + 1) * q;
5297  }
5298  k++;
5299  }
5300  for (k = 0; k < n; k++)
5301  {
5302  s[k] = b[k];
5303  t[k] = c[k];
5304  }
5305  // Initialization of U and V.
5306  U.Resize(m,m);
5307  U.Zero();
5308  for(j = 0; j < m; ++j) U(j,j) = 1.0;
5309  V.Resize(n,n);
5310  V.Zero();
5311  for(j = 0; j < n; ++j) V(j,j) = 1.0;
5312  // QR diagonalization.
5313  for (k = n - 1; k >= 0; k--)
5314  {
5315  // Test for split.
5316  while (true)
5317  {
5318  bool skip = false;
5319  for (l = k; l >= 0; l--)
5320  {
5321  if (!fabs(get_value(t[l])))
5322  {
5323  skip = true;
5324  break;
5325  }
5326  if (!fabs(get_value(s[l - 1]))) break;
5327  }
5328  // Cancellation of E(L).
5329  if (!skip)
5330  {
5331  cs = 0.0;
5332  sn = 1.0;
5333  for (i = l; i <= k; ++i)
5334  {
5335  f = sn * t[i];
5336  t[i] = cs * t[i];
5337  if (!fabs(get_value(f))) break;
5338  h = s[i];
5339  w = sqrt(f * f + h * h);
5340  s[i] = w;
5341  cs = h / w;
5342  sn = -f / w;
5343  for (j = 0; j < n; ++j)
5344  {
5345  x = real_part(U(j, l - 1));
5346  y = real_part(U(j, i));
5347  U(j, l - 1) = x * cs + y * sn;
5348  U(j, i) = y * cs - x * sn;
5349  }
5350 
5351  if (p == 0) continue;
5352 
5353  for (j = n; j < n + p; ++j)
5354  {
5355  q = A(l - 1, j);
5356  r = A(i, j);
5357  A(l - 1, j) = q * cs + r * sn;
5358  A(i, j) = r * cs - q * sn;
5359  }
5360  }
5361  }
5362  // Test for convergence.
5363  w = s[k];
5364  if (l == k) break;
5365  // Origin shift.
5366  x = s[l];
5367  y = s[k - 1];
5368  g = t[k - 1];
5369  h = t[k];
5370  f = ((y - w) * (y + w) + (g - h) * (g + h)) / (2.0 * h * y);
5371  g = sqrt(f * f + 1.0);
5372  if (real_part(get_value(f)) < 0.0) g = -g;
5373  f = ((x - w) * (x + w) + (y / (f + g) - h) * h) / x;
5374  // QR step.
5375  cs = 1.0;
5376  sn = 1.0;
5377  for (i = l+1; i <= k; ++i)
5378  {
5379  g = t[i];
5380  y = s[i];
5381  h = sn * g;
5382  g = cs * g;
5383  w = sqrt(h * h + f * f);
5384  t[i - 1] = w;
5385  cs = f / w;
5386  sn = h / w;
5387  f = x * cs + g * sn;
5388  g = g * cs - x * sn;
5389  h = y * sn;
5390  y = y * cs;
5391  for (j = 0; j < n; ++j)
5392  {
5393  x = real_part(V(j, i - 1));
5394  w = real_part(V(j, i));
5395  V(j, i - 1) = x * cs + w * sn;
5396  V(j, i) = w * cs - x * sn;
5397  }
5398  w = sqrt(h * h + f * f);
5399  s[i - 1] = w;
5400  cs = f / w;
5401  sn = h / w;
5402  f = cs * g + sn * y;
5403  x = cs * y - sn * g;
5404  for (j = 0; j < n; ++j)
5405  {
5406  y = real_part(U(j, i - 1));
5407  w = real_part(U(j, i));
5408  U(j, i - 1) = y * cs + w * sn;
5409  U(j, i) = w * cs - y * sn;
5410  }
5411  if (p == 0) break;
5412  for (j = n; j < n + p; ++j)
5413  {
5414  q = A(i - 1, j);
5415  r = A(i, j);
5416  A(i - 1, j) = q * cs + r * sn;
5417  A(i, j) = r * cs - q * sn;
5418  }
5419  }
5420  t[l] = 0.0;
5421  t[k] = f;
5422  s[k] = x;
5423  }
5424  // Convergence.
5425  if (real_part(get_value(w)) >= 0.0) continue;
5426  s[k] = -w;
5427  for (j = 0; j < n; ++j) V(j, k) = -V(j, k);
5428  }
5429  // Sort the singular values.
5430  for (k = 0; k < n; ++k)
5431  {
5432  g = -1.0;
5433  j = k;
5434  for (i = k; i < n; ++i)
5435  {
5436  if (real_part(get_value(g)) < real_part(get_value(s[i])))
5437  {
5438  g = s[i];
5439  j = i;
5440  }
5441  }
5442  if (j == k) continue;
5443 
5444  s[j] = s[k];
5445  s[k] = g;
5446  // Interchange V(1:N, J) and V(1:N, K).
5447  for (i = 0; i < n; ++i)
5448  {
5449  q = V(i, j);
5450  V(i, j) = V(i, k);
5451  V(i, k) = q;
5452  }
5453  // Interchange U(1:N, J) and U(1:N, K).
5454  for (i = 0; i < n; ++i)
5455  {
5456  q = U(i, j);
5457  U(i, j) = U(i, k);
5458  U(i, k) = q;
5459  }
5460  // Interchange A(J, N1:NP) and A(K, N1:NP).
5461  if (p == 0) continue;
5462  for (i = n; i < n + p; ++i)
5463  {
5464  q = A(j, i);
5465  A(j, i) = A(k, i);
5466  A(k, i) = q;
5467  }
5468  }
5469  // Back transformation.
5470  for (k = n - 1; k >= 0; k--)
5471  {
5472  if (b[k] == 0.0) continue;
5473  q = -A(k, k) / fabs(A(k, k));
5474  for (j = 0; j < m; ++j)
5475  U(k, j) *= q;
5476  for (j = 0; j < m; ++j)
5477  {
5478  q = 0.0;
5479  for (i = k; i < m; ++i)
5480  q += conj(A(i, k)) * U(i, j);
5481  q = q / (fabs(A(k, k)) * b[k]);
5482  for (i = k; i < m; ++i)
5483  U(i, j) -= q * A(i, k);
5484  }
5485  }
5486  if (n > 1)
5487  {
5488  for (k = n - 2; k >= 0; k--)
5489  {
5490  if (c[k+1] == 0.0) continue;
5491  q = -conj(A(k, k+1)) / fabs(A(k, k+1));
5492  for (j = 0; j < n; ++j)
5493  V(k+1, j) *= q;
5494 
5495  for (j = 0; j < n; ++j)
5496  {
5497  q = 0.0;
5498  for (i = k+1; i < n; ++i)
5499  q += A(k, i) * V(i, j);
5500  q = q / (fabs(A(k, k+1)) * c[k+1]);
5501  for (i = k+1; i < n; ++i)
5502  V(i, j) -= q * conj(A(k, i));
5503  }
5504  }
5505  }
5506  Sigma.Resize(m, n);
5507  Sigma.Zero();
5508  for (i = 0; i < n; ++i)
5509  Sigma(i, i) = s[i];
5510  return true;
5511  }
5512 
5513  template<typename Var>
5514  bool AbstractMatrixReadOnly<Var>::SVD(AbstractMatrix<Var>& U, AbstractMatrix<Var>& Sigma, AbstractMatrix<Var>& V, bool order_singular_values, bool nonnegative) const
5515  {
5516  int flag, i, its, j, jj, k, l, nm;
5517  int n = Rows();
5518  int m = Cols();
5519  if (n >= m)
5520  {
5521  U = (*this);
5522  Sigma.Resize(m, m);
5523  Sigma.Zero();
5524  V.Resize(m, m);
5525  }
5526  else // n < m
5527  {
5528  U.Resize(n, n);
5529  Sigma.Resize(n, n);
5530  Sigma.Zero();
5531  V.Resize(m, n);
5532  }
5533  if (n < m)
5534  {
5535  bool success = Transpose().SVD(V, Sigma, U);
5536  if (success)
5537  {
5538  //U.Swap(V);
5539  //U = U.Transpose();
5540  //V = V.Transpose();
5541  return true;
5542  }
5543  else return false;
5544  } //m <= n
5545  Var c, f, h, s, x, y, z;
5546  Var g = 0.0, scale = 0.0;
5547  INMOST_DATA_REAL_TYPE anorm = 0.0;
5548  std::vector<Var> rv1(m);
5549  std::swap(n, m); //this how original algorithm takes it
5550  // Householder reduction to bidiagonal form
5551  for (i = 0; i < n; i++)
5552  {
5553  // left-hand reduction
5554  l = i + 1;
5555  rv1[i] = scale * g;
5556  g = s = scale = 0.0;
5557  if (i < m)
5558  {
5559  for (k = i; k < m; k++) scale += fabs(U(k, i));
5560  if (fabs(get_value(scale)))
5561  {
5562  for (k = i; k < m; k++)
5563  {
5564  U(k, i) /= scale;
5565  //s += U(k, i) * U(k, i); //original
5566  s += U(k, i) * conj(U(k, i)); //for complex
5567  }
5568  f = U(i, i);
5569  g = -sign_func(sqrt(s), f);
5570  h = f * g - s;
5571  U(i, i) = f - g;
5572  if (i != n - 1 && fabs(get_value(h)))
5573  {
5574  for (j = l; j < n; j++)
5575  {
5576  //for (s = 0.0, k = i; k < m; k++) s += U(k, i) * U(k, j);
5577  for (s = 0.0, k = i; k < m; k++) s += conj(U(k, i)) * U(k, j);
5578  f = s / h;
5579  for (k = i; k < m; k++) U(k, j) += f * U(k, i);
5580  }
5581  }
5582  for (k = i; k < m; k++) U(k, i) *= scale;
5583  }
5584  }
5585  Sigma(i, i) = scale * g;
5586  // right-hand reduction
5587  g = s = scale = 0.0;
5588  if (i < m && i != n - 1)
5589  {
5590  for (k = l; k < n; k++) scale += fabs(U(i, k));
5591  if (fabs(get_value(scale)))
5592  {
5593  for (k = l; k < n; k++)
5594  {
5595  U(i, k) = U(i, k) / scale;
5596  //s += U(i, k) * U(i, k); //original
5597  s += U(i, k) * conj(U(i, k));
5598  }
5599  f = U(i, l);
5600  g = -sign_func(sqrt(s), f);
5601  h = f * g - s;
5602  U(i, l) = f - g;
5603  for (k = l; k < n; k++) rv1[k] = U(i, k) / h;
5604  if (i != m - 1)
5605  {
5606  for (j = l; j < m; j++)
5607  {
5608  for (s = 0.0, k = l; k < n; k++) s += conj(U(i, k)) * U(j, k);
5609  for (k = l; k < n; k++) U(j, k) += s * rv1[k];
5610  }
5611  }
5612  for (k = l; k < n; k++) U(i, k) *= scale;
5613  }
5614  }
5615  anorm = max_func(anorm, fabs(get_value(Sigma(i, i))) + fabs(get_value(rv1[i])));
5616  }
5617 
5618  // accumulate the right-hand transformation
5619  for (i = n - 1; i >= 0; i--)
5620  {
5621  if (i < (n - 1))
5622  {
5623  if (fabs(get_value(g)))
5624  {
5625  for (j = l; j < n; j++) V(j, i) = ((U(i, j) / U(i, l)) / g);
5626  // double division to avoid underflow
5627  for (j = l; j < n; j++)
5628  {
5629  for (s = 0.0, k = l; k < n; k++) s += U(i, k) * V(k, j);
5630  for (k = l; k < n; k++) V(k, j) += s * V(k, i);
5631  }
5632  }
5633  for (j = l; j < n; j++) V(i, j) = V(j, i) = 0.0;
5634  }
5635  V(i, i) = 1.0;
5636  g = rv1[i];
5637  l = i;
5638  }
5639 
5640  // accumulate the left-hand transformation
5641  for (i = n - 1; i >= 0; i--)
5642  {
5643  l = i + 1;
5644  g = Sigma(i, i);
5645  if (i < (n - 1))
5646  for (j = l; j < n; j++)
5647  U(i, j) = 0.0;
5648  if (fabs(get_value(g)))
5649  {
5650  g = 1.0 / g;
5651  if (i != n - 1)
5652  {
5653  for (j = l; j < n; j++)
5654  {
5655  for (s = 0.0, k = l; k < m; k++) s += conj(U(k, i)) * U(k, j);
5656  f = (s / U(i, i)) * g;
5657  for (k = i; k < m; k++) U(k, j) += f * U(k, i);
5658  }
5659  }
5660  for (j = i; j < m; j++) U(j, i) *= g;
5661  }
5662  else for (j = i; j < m; j++) U(j, i) = 0.0;
5663  U(i, i) += 1;
5664  }
5665 
5666  // diagonalize the bidiagonal form
5667  for (k = n - 1; k >= 0; k--)
5668  {// loop over singular values
5669  for (its = 0; its < 30; its++)
5670  {// loop over allowed iterations
5671  flag = 1;
5672  for (l = k; l >= 0; l--)
5673  {// test for splitting
5674  nm = l - 1;
5675  if (fabs(get_value(rv1[l])) + anorm == anorm)
5676  {
5677  flag = 0;
5678  break;
5679  }
5680  if (l == 0) break;
5681  if (fabs(get_value(Sigma(nm, nm))) + anorm == anorm)
5682  break;
5683  }
5684  if (flag && l != 0)
5685  {
5686  c = 0.0;
5687  s = 1.0;
5688  for (i = l; i <= k; i++)
5689  {
5690  f = s * rv1[i];
5691  if (fabs(get_value(f)) + anorm != anorm)
5692  {
5693  g = Sigma(i, i);
5694  h = pythag(f, g);
5695  Sigma(i, i) = h;
5696  h = 1.0 / h;
5697  c = g * h;
5698  s = (-f * h);
5699  for (j = 0; j < m; j++)
5700  {
5701  y = U(j, nm);
5702  z = U(j, i);
5703  U(j, nm) = (y * c + z * s);
5704  U(j, i) = (z * c - y * s);
5705  }
5706  }
5707  }
5708  }
5709  z = Sigma(k, k);
5710  if (l == k)
5711  {// convergence
5712  if (real_part(get_value(z)) < 0.0 && nonnegative)
5713  {// make singular value nonnegative
5714  Sigma(k, k) = -z;
5715  for (j = 0; j < n; j++) V(j, k) = -V(j, k);
5716  }
5717  break;
5718  }
5719  if (its >= 30)
5720  {
5721  std::cout << "No convergence after " << its << " iterations" << std::endl;
5722  std::swap(n, m);
5723  return false;
5724  }
5725  // shift from bottom 2 x 2 minor
5726  x = Sigma(l, l);
5727  nm = k - 1;
5728  y = Sigma(nm, nm);
5729  g = rv1[nm];
5730  h = rv1[k];
5731  f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
5732  g = pythag(f, 1.0);
5733  f = ((x - z) * (x + z) + h * ((y / (f + sign_func(g, f))) - h)) / x;
5734  // next QR transformation
5735  c = s = 1.0;
5736  for (j = l; j <= nm; j++)
5737  {
5738  i = j + 1;
5739  g = rv1[i];
5740  y = Sigma(i, i);
5741  h = s * g;
5742  g = c * g;
5743  z = pythag(f, h);
5744  rv1[j] = z;
5745  c = f / z;
5746  s = h / z;
5747  f = x * c + g * s;
5748  g = g * c - x * s;
5749  h = y * s;
5750  y = y * c;
5751  for (jj = 0; jj < n; jj++)
5752  {
5753  x = V(jj, j);
5754  z = V(jj, i);
5755  V(jj, j) = (x * c + z * s);
5756  V(jj, i) = (z * c - x * s);
5757  }
5758  z = pythag(f, h);
5759  Sigma(j, j) = z;
5760  if (fabs(get_value(z)))
5761  {
5762  z = 1.0 / z;
5763  c = f * z;
5764  s = h * z;
5765  }
5766  f = (c * g) + (s * y);
5767  x = (c * y) - (s * g);
5768  for (jj = 0; jj < m; jj++)
5769  {
5770  y = U(jj, j);
5771  z = U(jj, i);
5772  U(jj, j) = (y * c + z * s);
5773  U(jj, i) = (z * c - y * s);
5774  }
5775  }
5776  rv1[l] = 0.0;
5777  rv1[k] = f;
5778  Sigma(k, k) = x;
5779  }
5780  }
5781  //CHECK THIS!
5782  if (order_singular_values)
5783  {
5784  Var temp;
5785  for (i = 0; i < n; i++)
5786  {
5787  k = i;
5788  for (j = i + 1; j < n; ++j)
5789  if (real_part(get_value(Sigma(k, k))) < real_part(get_value(Sigma(j, j)))) k = j;
5790  if (real_part(get_value(Sigma(k, k))) > real_part(get_value(Sigma(i, i))))
5791  {
5792  temp = Sigma(k, k);
5793  Sigma(k, k) = Sigma(i, i);
5794  Sigma(i, i) = temp;
5795  // U is m by n
5796  for (int j = 0; j < m; ++j)
5797  {
5798  temp = U(j, k);
5799  U(j, k) = U(j, i);
5800  U(j, i) = temp;
5801  }
5802  // V is n by n
5803  for (int j = 0; j < n; ++j)
5804  {
5805  temp = V(j, k);
5806  V(j, k) = V(j, i);
5807  V(j, i) = temp;
5808  }
5809  }
5810  }
5811  }
5812 
5813  std::swap(n, m);
5814  return true;
5815  }
5816 
5818  typedef Matrix<INMOST_DATA_INTEGER_TYPE> iMatrix;
5820  typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix;
5822  typedef Matrix<INMOST_DATA_CPLX_TYPE> cMatrix;
5824  typedef Matrix<INMOST_DATA_INTEGER_TYPE,shell<INMOST_DATA_INTEGER_TYPE> > iaMatrix;
5826  typedef Matrix<INMOST_DATA_REAL_TYPE,shell<INMOST_DATA_REAL_TYPE> > raMatrix;
5828  typedef Matrix<INMOST_DATA_CPLX_TYPE, shell<INMOST_DATA_CPLX_TYPE> > caMatrix;
5830  typedef SymmetricMatrix<INMOST_DATA_INTEGER_TYPE> iSymmetricMatrix;
5832  typedef SymmetricMatrix<INMOST_DATA_REAL_TYPE> rSymmetricMatrix;
5834  typedef SymmetricMatrix<INMOST_DATA_CPLX_TYPE> cSymmetricMatrix;
5836  typedef SymmetricMatrix<INMOST_DATA_INTEGER_TYPE,shell<INMOST_DATA_INTEGER_TYPE> > iaSymmetricMatrix;
5838  typedef SymmetricMatrix<INMOST_DATA_REAL_TYPE,shell<INMOST_DATA_REAL_TYPE> > raSymmetricMatrix;
5840  typedef SymmetricMatrix<INMOST_DATA_CPLX_TYPE, shell<INMOST_DATA_CPLX_TYPE> > caSymmetricMatrix;
5842  __INLINE iaMatrix iaMatrixMake(INMOST_DATA_INTEGER_TYPE * p, iaMatrix::enumerator n, iaMatrix::enumerator m) {return iaMatrix(shell<INMOST_DATA_INTEGER_TYPE>(p,n*m),n,m);}
5843  __INLINE raMatrix raMatrixMake(INMOST_DATA_REAL_TYPE * p, raMatrix::enumerator n, raMatrix::enumerator m) {return raMatrix(shell<INMOST_DATA_REAL_TYPE>(p,n*m),n,m);}
5844  __INLINE caMatrix caMatrixMake(INMOST_DATA_CPLX_TYPE* p, raMatrix::enumerator n, caMatrix::enumerator m) { return caMatrix(shell<INMOST_DATA_CPLX_TYPE>(p, n * m), n, m); }
5845  __INLINE iaSymmetricMatrix iaSymmetricMatrixMake(INMOST_DATA_INTEGER_TYPE * p, iaSymmetricMatrix::enumerator n) {return iaSymmetricMatrix(shell<INMOST_DATA_INTEGER_TYPE>(p,n*(n+1)/2),n);}
5846  __INLINE raSymmetricMatrix raSymmetricMatrixMake(INMOST_DATA_REAL_TYPE * p, raSymmetricMatrix::enumerator n) {return raSymmetricMatrix(shell<INMOST_DATA_REAL_TYPE>(p,n*(n+1)/2),n);}
5847  __INLINE caSymmetricMatrix caSymmetricMatrixMake(INMOST_DATA_CPLX_TYPE* p, caSymmetricMatrix::enumerator n) { return caSymmetricMatrix(shell<INMOST_DATA_CPLX_TYPE>(p, n * (n + 1) / 2), n); }
5848 #if defined(USE_AUTODIFF)
5850  typedef Matrix<unknown> uMatrix;
5852  typedef Matrix<variable> vMatrix;
5853  //< shortcut for matrix of variables with first and second order derivatives.
5854  typedef Matrix<hessian_variable> hMatrix;
5856  typedef Matrix<unknown,shell<unknown> > uaMatrix;
5858  typedef Matrix<variable,shell<variable> > vaMatrix;
5860  typedef Matrix<hessian_variable,shell<hessian_variable> > haMatrix;
5862  typedef SymmetricMatrix<unknown,shell<unknown> > uaSymmetricMatrix;
5864  typedef SymmetricMatrix<variable,shell<variable> > vaSymmetricMatrix;
5866  typedef SymmetricMatrix<hessian_variable,shell<hessian_variable> > haSymmetricMatrix;
5867 
5868 
5869 
5870  __INLINE uaMatrix uaMatrixMake(unknown * p, uaMatrix::enumerator n, uaMatrix::enumerator m) {return uaMatrix(shell<unknown>(p,n*m),n,m);}
5871  __INLINE vaMatrix vaMatrixMake(variable * p, vaMatrix::enumerator n, vaMatrix::enumerator m) {return vaMatrix(shell<variable>(p,n*m),n,m);}
5872  __INLINE haMatrix vaMatrixMake(hessian_variable * p, haMatrix::enumerator n, haMatrix::enumerator m) {return haMatrix(shell<hessian_variable>(p,n*m),n,m);}
5873  __INLINE uaSymmetricMatrix uaSymmetricMatrixMake(unknown * p, uaSymmetricMatrix::enumerator n) {return uaSymmetricMatrix(shell<unknown>(p,n*(n+1)/2),n);}
5874  __INLINE vaSymmetricMatrix vaSymmetricMatrixMake(variable * p, vaSymmetricMatrix::enumerator n) {return vaSymmetricMatrix(shell<variable>(p,n*(n+1)/2),n);}
5875  __INLINE haSymmetricMatrix vaSymmetricMatrixMake(hessian_variable * p, haSymmetricMatrix::enumerator n) {return haSymmetricMatrix(shell<hessian_variable>(p,n*(n+1)/2),n);}
5876 #endif
5877 }
5882 template<typename typeB>
5883 //INMOST::Matrix<typename INMOST::Promote<INMOST_DATA_REAL_TYPE,typeB>::type>
5884 INMOST::MatrixMulCoef<typeB, INMOST_DATA_REAL_TYPE, typename INMOST::Promote<INMOST_DATA_REAL_TYPE, typeB>::type>
5885 operator *(const INMOST_DATA_REAL_TYPE& coef, const INMOST::AbstractMatrixReadOnly<typeB>& other)
5886 //{return other*coef;}
5887 {return INMOST::MatrixMulCoef<typeB, INMOST_DATA_REAL_TYPE, typename INMOST::Promote<INMOST_DATA_REAL_TYPE, typeB>::type>(other, coef);}
5888 #if defined(USE_AUTODIFF)
5894 template<typename typeB>
5895 //INMOST::Matrix<typename INMOST::Promote<INMOST::unknown,typeB>::type>
5896 INMOST::MatrixMulCoef<typeB, INMOST::unknown, typename INMOST::Promote<INMOST::unknown, typeB>::type>
5897 operator *(const INMOST::unknown& coef, const INMOST::AbstractMatrixReadOnly<typeB>& other)
5898 //{return other*coef;}
5899 {return INMOST::MatrixMulCoef<typeB, INMOST::unknown, typename INMOST::Promote<INMOST::unknown, typeB>::type>(other, coef);}
5905 template<typename typeB>
5906 //INMOST::Matrix<typename INMOST::Promote<INMOST::variable,typeB>::type>
5907 INMOST::MatrixMulCoef<typeB, INMOST::variable, typename INMOST::Promote<INMOST::variable, typeB>::type>
5908 operator *(const INMOST::variable& coef, const INMOST::AbstractMatrixReadOnly<typeB>& other)
5909 //{return other*coef;}
5910 {return INMOST::MatrixMulCoef<typeB, INMOST::variable, typename INMOST::Promote<INMOST::variable, typeB>::type>(other, coef);}
5917 template<typename typeB>
5918 //INMOST::Matrix<typename INMOST::Promote<INMOST::hessian_variable,typeB>::type>
5919 INMOST::MatrixMulCoef<typeB, INMOST::hessian_variable, typename INMOST::Promote<INMOST::hessian_variable, typeB>::type>
5920 operator *(const INMOST::hessian_variable& coef, const INMOST::AbstractMatrixReadOnly<typeB>& other)
5921 //{return other*coef;}
5922 {return INMOST::MatrixMulCoef<typeB, INMOST::hessian_variable, typename INMOST::Promote<INMOST::hessian_variable, typeB>::type>(other, coef);}
5927 template<class A, typename typeB>
5928 //INMOST::Matrix<typename INMOST::Promote<INMOST::variable,typeB>::type>
5929 INMOST::MatrixMulShellCoef<typeB, INMOST::variable, typename INMOST::Promote<INMOST::variable, typeB>::type>
5930 operator *(INMOST::shell_expression<A> const& coef, const INMOST::AbstractMatrixReadOnly<typeB>& other)
5931 //{return other*INMOST::variable(coef);}
5932 {return INMOST::MatrixMulShellCoef<typeB, INMOST::variable, typename INMOST::Promote<INMOST::variable, typeB>::type>(other, INMOST::variable(coef));}
5933 #endif
5934 
5935 template<typename T>
5936 __INLINE bool check_nans(const INMOST::AbstractMatrixReadOnly<T> & A) {return A.CheckNans();}
5937 template<typename T>
5938 __INLINE bool check_infs(const INMOST::AbstractMatrixReadOnly<T> & A) {return A.CheckInfs();}
5939 template<typename T>
5940 __INLINE bool check_nans_infs(const INMOST::AbstractMatrixReadOnly<T> & A) {return A.CheckNansInfs();}
5941 
5942 #endif //INMOST_DENSE_INCLUDED
INMOST::AbstractMatrixReadOnly::FrobeniusNorm
SelfPromote< Var >::type FrobeniusNorm() const
Computes frobenius norm of the matrix.
Definition: inmost_dense.h:449
INMOST::AbstractMatrixReadOnly::CholeskyInvert
Matrix< Var > CholeskyInvert(int *ierr=NULL) const
Inverts symmetric positive-definite matrix using Cholesky decomposition.
Definition: inmost_dense.h:4109
INMOST::AbstractMatrixReadOnly::operator/
Matrix< typename Promote< Var, typeB >::type > operator/(const AbstractMatrixReadOnly< typeB > &other) const
Performs B^{-1}*A, multiplication by inverse matrix from left.
Definition: inmost_dense.h:580
INMOST::AbstractMatrixReadOnly::ConcatCols
ConstMatrixConcatCols2< Var, VarB, typename Promote< Var, VarB >::type > ConcatCols(const AbstractMatrixReadOnly< VarB > &B) const
Concatenate B matrix as columns of current matrix.
Definition: inmost_dense.h:604
INMOST::AbstractMatrixReadOnly::DotProduct
Promote< Var, typeB >::type DotProduct(const AbstractMatrixReadOnly< typeB > &other) const
Computes dot product by summing up multiplication of entries with the same indices in the current and...
Definition: inmost_dense.h:428
INMOST::AbstractMatrixReadOnly::GetMatrixCount
virtual __INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:633
INMOST::AbstractMatrixReadOnly::Power
Matrix< Var > Power(INMOST_DATA_REAL_TYPE n, int *ierr=NULL) const
Calcuate A^n, where n is some real value.
Definition: inmost_dense.h:4817
INMOST::AbstractMatrixReadOnly::operator-
MatrixUnaryMinus< Var > operator-() const
Unary minus. Change sign of each element of the matrix.
Definition: inmost_dense.h:324
INMOST::AbstractMatrixReadOnly::CheckInfs
bool CheckInfs() const
Check all matrix entries for infinity.
Definition: inmost_dense.h:225
INMOST::AbstractMatrixReadOnly::operator*
MatrixMulShellCoef< Var, shell_expression< A >, typename Promote< Var, variable >::type > operator*(shell_expression< A > const &coef) const
Multiply the matrix by a coefficient.
Definition: inmost_dense.h:4045
INMOST::AbstractMatrixReadOnly::operator-
MatrixDifference< Var, typeB > operator-(const AbstractMatrixReadOnly< typeB > &other) const
Subtract a matrix.
Definition: inmost_dense.h:3849
INMOST::AbstractMatrixReadOnly::Cols
virtual enumerator Cols() const =0
Obtain number of columns.
INMOST::AbstractMatrixReadOnly::MaxNorm
Var MaxNorm() const
Computes maximum absolute value of the matrix.
Definition: inmost_dense.h:459
INMOST::AbstractMatrixReadOnly::compute
virtual Var compute(enumerator i, enumerator j) const =0
Compute element of the matrix by row and column indices without right to change the element.
INMOST::AbstractMatrixReadOnly::Kronecker
KroneckerProduct< Var, typeB > Kronecker(const AbstractMatrixReadOnly< typeB > &other) const
Kronecker product, latex symbol \otimes.
Definition: inmost_dense.h:4093
INMOST::AbstractMatrixReadOnly::ConjugateTranspose
ConstMatrixConjugateTranspose< Var > ConjugateTranspose() const
Transpose and conjugate current matrix.
Definition: inmost_dense.h:283
INMOST::AbstractMatrixReadOnly::Invert
Matrix< Var > Invert(int *ierr=NULL) const
Inverts matrix using Crout-LU decomposition with full pivoting for maximum element.
Definition: inmost_dense.h:4100
INMOST::AbstractMatrixReadOnly::ConcatCols
ConstMatrixConcatCols< Var > ConcatCols(const AbstractMatrixReadOnly< Var > &B) const
Concatenate B matrix as columns of current matrix.
Definition: inmost_dense.h:592
INMOST::AbstractMatrixReadOnly::operator/
MatrixDivCoef< Var, typeB, typename Promote< Var, typeB >::type > operator/(const typeB &coef) const
Divide the matrix by a coefficient of a different type.
Definition: inmost_dense.h:4074
INMOST::AbstractMatrixReadOnly::PseudoInvert
Matrix< Var > PseudoInvert(INMOST_DATA_REAL_TYPE tol=0, int *ierr=NULL) const
Calculates Moore-Penrose pseudo-inverse of the matrix.
Definition: inmost_dense.h:4762
INMOST::AbstractMatrixReadOnly::operator()
ConstSubMatrix< Var > operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col) const
Extract submatrix of a matrix for in-place manipulation of elements.
Definition: inmost_dense.h:4877
INMOST::AbstractMatrixReadOnly::Trace
Var Trace() const
Calculate sum of the diagonal elements of the matrix.
Definition: inmost_dense.h:379
INMOST::AbstractMatrixReadOnly::CholeskySolve
Matrix< typename Promote< Var, typeB >::type > CholeskySolve(const AbstractMatrixReadOnly< typeB > &B, int *ierr=NULL) const
Finds X in A*X=B, where A is a square symmetric positive definite matrix.
Definition: inmost_dense.h:4120
INMOST::AbstractMatrixReadOnly::SVD
bool SVD(AbstractMatrix< Var > &U, AbstractMatrix< Var > &Sigma, AbstractMatrix< Var > &V, bool order_singular_values=true, bool nonnegative=true) const
Singular value decomposition.
Definition: inmost_dense.h:5514
INMOST::AbstractMatrixReadOnly::operator*
MatrixMulCoef< Var, typeB, typename Promote< Var, typeB >::type > operator*(const typeB &coef) const
Multiply the matrix by a coefficient.
Definition: inmost_dense.h:4036
INMOST::AbstractMatrixReadOnly::Print
void Print(INMOST_DATA_REAL_TYPE threshold=1.0e-10, std::ostream &sout=std::cout) const
Output matrix to screen.
Definition: inmost_dense.h:389
INMOST::AbstractMatrixReadOnly::isSymmetric
bool isSymmetric(double eps=1.0e-7) const
Check if the matrix is symmetric.
Definition: inmost_dense.h:411
INMOST::AbstractMatrixReadOnly::ConcatRows
ConstMatrixConcatRows< Var > ConcatRows(const AbstractMatrixReadOnly< Var > &B) const
Concatenate B matrix as rows of current matrix.
Definition: inmost_dense.h:614
INMOST::AbstractMatrixReadOnly::operator^
Promote< Var, typeB >::type operator^(const AbstractMatrixReadOnly< typeB > &other) const
Computes dot product by summing up multiplication of entries with the same indices in the current and...
Definition: inmost_dense.h:443
INMOST::AbstractMatrixReadOnly::CrossProduct
Matrix< typename Promote< Var, typeB >::type > CrossProduct(const AbstractMatrixReadOnly< typeB > &other) const
Cross-product operation for a vector.
Definition: inmost_dense.h:3676
INMOST::AbstractMatrixReadOnly::cSVD
bool cSVD(AbstractMatrix< Var > &U, AbstractMatrix< Var > &Sigma, AbstractMatrix< Var > &V) const
Singular value decomposition.
Definition: inmost_dense.h:5209
INMOST::AbstractMatrixReadOnly::Transpose
ConstMatrixTranspose< Var > Transpose() const
Transpose current matrix.
Definition: inmost_dense.h:280
INMOST::AbstractMatrixReadOnly::Conjugate
ConstMatrixConjugate< Var > Conjugate() const
Conjugate current matrix.
Definition: inmost_dense.h:286
INMOST::AbstractMatrixReadOnly::Rows
virtual enumerator Rows() const =0
Obtain number of rows.
INMOST::AbstractMatrixReadOnly::ExtractSubMatrix
Matrix< Var > ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const
Extract submatrix of a matrix.
Definition: inmost_dense.h:4745
INMOST::AbstractMatrixReadOnly::PseudoSolve
Matrix< typename Promote< Var, typeB >::type > PseudoSolve(const AbstractMatrixReadOnly< typeB > &B, INMOST_DATA_REAL_TYPE tol=0, int *ierr=NULL) const
Solves the system of equations of the form A*X=B, with A and B matrices.
Definition: inmost_dense.h:4845
INMOST::AbstractMatrixReadOnly::ConcatRows
ConstMatrixConcatRows2< Var, VarB, typename Promote< Var, VarB >::type > ConcatRows(const AbstractMatrixReadOnly< VarB > &B) const
Concatenate B matrix as rows of current matrix.
Definition: inmost_dense.h:626
INMOST::AbstractMatrixReadOnly::operator*
MatrixMul< Var, typeB, typename Promote< Var, typeB >::type > operator*(const AbstractMatrixReadOnly< typeB > &other) const
Multiply the matrix by another matrix.
Definition: inmost_dense.h:3917
INMOST::AbstractMatrixReadOnly::BlockOf
ConstBlockOfMatrix< Var > BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col) const
Define matrix as a part of a matrix of larger size with in-place manipulation of elements.
Definition: inmost_dense.h:4888
INMOST::AbstractMatrixReadOnly::~AbstractMatrixReadOnly
virtual ~AbstractMatrixReadOnly()
Destructor.
Definition: inmost_dense.h:631
INMOST::AbstractMatrixReadOnly::CheckNans
bool CheckNans() const
Check all matrix entries for not a number.
Definition: inmost_dense.h:222
INMOST::AbstractMatrixReadOnly::operator+
MatrixSum< Var, typeB > operator+(const AbstractMatrixReadOnly< typeB > &other) const
Add a matrix.
Definition: inmost_dense.h:3883
INMOST::AbstractMatrixReadOnly::Repack
ConstMatrixRepack< Var > Repack(enumerator rows, enumerator cols) const
Change representation of the matrix into matrix of another size.
Definition: inmost_dense.h:519
INMOST::AbstractMatrixReadOnly::CheckNansInfs
bool CheckNansInfs() const
Check all matrix entries for not a number and infinity.
Definition: inmost_dense.h:228
INMOST::AbstractMatrixReadOnly::Transform
Matrix< typename Promote< Var, typeB >::type > Transform(const AbstractMatrixReadOnly< typeB > &other) const
Transformation matrix from current vector to provided vector using shortest arc rotation.
Definition: inmost_dense.h:3735
INMOST::AbstractMatrixReadOnly::Solve
Matrix< typename Promote< Var, typeB >::type > Solve(const AbstractMatrixReadOnly< typeB > &B, int *ierr=NULL) const
Finds X in A*X=B, where A and B are general matrices.
Definition: inmost_dense.h:4588
INMOST::AbstractMatrixReadOnly::operator/
MatrixDivShellCoef< Var, shell_expression< A >, typename Promote< Var, variable >::type > operator/(shell_expression< A > const &coef) const
Divide the matrix by a coefficient of a different type.
Definition: inmost_dense.h:4065
INMOST::AbstractMatrixReadOnly::Det
Var Det() const
Matrix determinant.
Definition: inmost_dense.h:230
INMOST::AbstractMatrixReadOnly::Root
Matrix< Var > Root(INMOST_DATA_ENUM_TYPE iter=25, INMOST_DATA_REAL_TYPE tol=1.0e-7, int *ierr=NULL) const
Calculate square root of A matrix by Babylonian method.
Definition: inmost_dense.h:4792
INMOST::AbstractMatrix
Abstract class for a matrix, used to abstract away all the data storage and access and provide common...
Definition: inmost_dense.h:647
INMOST::AbstractMatrix::CrossProduct
Matrix< typename Promote< Var, typeB >::type > CrossProduct(const AbstractMatrixReadOnly< typeB > &other) const
Cross-product operation for a vector.
Definition: inmost_dense.h:3715
INMOST::AbstractMatrix::Print
void Print(INMOST_DATA_REAL_TYPE threshold=1.0e-10, std::ostream &sout=std::cout) const
Output matrix to screen.
Definition: inmost_dense.h:955
INMOST::AbstractMatrix::DotProduct
Promote< Var, typeB >::type DotProduct(const AbstractMatrixReadOnly< typeB > &other) const
Computes dot product by summing up multiplication of entries with the same indices in the current and...
Definition: inmost_dense.h:895
INMOST::AbstractMatrix::Trace
Var Trace() const
Calculate sum of the diagonal elements of the matrix.
Definition: inmost_dense.h:945
INMOST::AbstractMatrix::CheckNansInfs
bool CheckNansInfs() const
Check all matrix entries for not a number and infinity.
Definition: inmost_dense.h:836
INMOST::AbstractMatrix::BlockOf
BlockOfMatrix< Var > BlockOf(enumerator nrows, enumerator ncols, enumerator offset_row, enumerator offset_col)
Define matrix as a part of a matrix of larger size with in-place manipulation of elements.
Definition: inmost_dense.h:4883
INMOST::AbstractMatrix::ConcatCols
MatrixConcatCols< Var > ConcatCols(AbstractMatrix< Var > &B)
Concatenate B matrix as columns of current matrix.
Definition: inmost_dense.h:808
INMOST::AbstractMatrix::Zero
void Zero()
Set all the elements of the matrix to zero.
Definition: inmost_dense.h:730
INMOST::AbstractMatrix::operator+=
AbstractMatrix & operator+=(const AbstractMatrixReadOnly< typeB > &other)
Add a matrix and store result in the current.
INMOST::AbstractMatrix::Swap
virtual void Swap(AbstractMatrix< Var > &b)
Exchange contents of two matrices.
Definition: inmost_dense.h:3666
INMOST::AbstractMatrix::operator()
virtual Var & operator()(enumerator i, enumerator j)=0
Access element of the matrix by row and column indices.
INMOST::AbstractMatrix::FrobeniusNorm
SelfPromote< Var >::type FrobeniusNorm() const
Computes frobenious norm of the matrix.
Definition: inmost_dense.h:925
INMOST::AbstractMatrix::get
virtual const Var & get(enumerator i, enumerator j) const =0
Access element of the matrix by row and column indices.
INMOST::AbstractMatrix::operator/=
AbstractMatrix & operator/=(typeB coef)
Divide the matrix by the coefficient of the same type and store the result.
INMOST::AbstractMatrix::Transpose
MatrixTranspose< Var > Transpose()
Transpose the current matrix with access to elements.
Definition: inmost_dense.h:796
INMOST::AbstractMatrix::ConcatRows
MatrixConcatRows< Var > ConcatRows(AbstractMatrix< Var > &B)
Concatenate B matrix as rows of current matrix.
Definition: inmost_dense.h:815
INMOST::AbstractMatrix::Transform
Matrix< typename Promote< Var, typeB >::type > Transform(const AbstractMatrix< typeB > &other) const
Transformation matrix from current vector to provided vector using shortest arc rotation.
Definition: inmost_dense.h:3774
INMOST::AbstractMatrix::~AbstractMatrix
virtual ~AbstractMatrix()
Destructor.
Definition: inmost_dense.h:1008
INMOST::AbstractMatrix::operator()
SubMatrix< Var > operator()(enumerator first_row, enumerator last_row, enumerator first_col, enumerator last_col)
Extract submatrix of a matrix for in-place manipulation of elements.
Definition: inmost_dense.h:4872
INMOST::AbstractMatrix::Repack
MatrixRepack< Var > Repack(enumerator rows, enumerator cols)
Change representation of the matrix into matrix of another size.
Definition: inmost_dense.h:801
INMOST::AbstractMatrix::CheckNans
bool CheckNans() const
Check all matrix entries for not a number.
Definition: inmost_dense.h:818
INMOST::AbstractMatrix::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:999
INMOST::AbstractMatrix::Transform
Matrix< typename Promote< Var, typeB >::type > Transform(const AbstractMatrixReadOnly< typeB > &other) const
Transformation matrix from current vector to provided vector using shortest arc rotation.
Definition: inmost_dense.h:3812
INMOST::AbstractMatrix::operator+=
AbstractMatrix & operator+=(const AbstractMatrix< typeB > &other)
Add a matrix and store result in the current.
INMOST::AbstractMatrix::operator*=
AbstractMatrix & operator*=(const AbstractMatrixReadOnly< typeB > &B)
Multiply matrix with another matrix in-place.
INMOST::AbstractMatrix::DotProduct
Promote< Var, typeB >::type DotProduct(const AbstractMatrix< typeB > &other) const
Computes dot product by summing up multiplication of entries with the same indices in the current and...
Definition: inmost_dense.h:879
INMOST::AbstractMatrix::operator*=
AbstractMatrix & operator*=(typeB coef)
Multiply the matrix by the coefficient of the same type and store the result.
INMOST::AbstractMatrix::operator-=
AbstractMatrix & operator-=(const AbstractMatrixReadOnly< typeB > &other)
Subtract a matrix and store result in the current.
INMOST::AbstractMatrix::CrossProduct
Matrix< typename Promote< Var, typeB >::type > CrossProduct(const AbstractMatrix< typeB > &other) const
Cross-product operation for a vector.
Definition: inmost_dense.h:3696
INMOST::AbstractMatrix::AbstractMatrix
AbstractMatrix()
Construct empty matrix.
Definition: inmost_dense.h:660
INMOST::AbstractMatrix::isSymmetric
bool isSymmetric(double eps=1.0e-7) const
Check if the matrix is symmetric.
Definition: inmost_dense.h:977
INMOST::AbstractMatrix::operator-=
AbstractMatrix & operator-=(const AbstractMatrix< typeB > &other)
Subtract a matrix and store result in the current.
INMOST::AbstractMatrix::operator^
Promote< Var, typeB >::type operator^(const AbstractMatrix< typeB > &other) const
Computes dot product by summing up multiplication of entries with the same indices in the current and...
Definition: inmost_dense.h:910
INMOST::AbstractMatrix::MPT
void MPT(INMOST_DATA_ENUM_TYPE *Perm, INMOST_DATA_REAL_TYPE *SL=NULL, INMOST_DATA_REAL_TYPE *SR=NULL) const
Maximum product transversal.
Definition: inmost_dense.h:4986
INMOST::AbstractMatrix::MaxNorm
Var MaxNorm() const
Computes maximum absolute value of the matrix.
Definition: inmost_dense.h:935
INMOST::AbstractMatrix::CheckInfs
bool CheckInfs() const
Check all matrix entries for infinity.
Definition: inmost_dense.h:827
INMOST::AbstractMatrix::Resize
virtual void Resize(enumerator rows, enumerator cols)=0
Resize the matrix into different size.
INMOST::AbstractMatrix::operator=
AbstractMatrix & operator=(AbstractMatrix const &other)
Assign matrix of the same type.
Definition: inmost_dense.h:685
INMOST::BlockOfMatrix
This class allows to address a matrix as a block of an empty matrix of larger size.
Definition: inmost_dense.h:2269
INMOST::BlockOfMatrix::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:2327
INMOST::BlockOfMatrix::operator=
BlockOfMatrix & operator=(AbstractMatrix< typeB > const &other)
Assign matrix of another type to the block of matrix.
Definition: inmost_dense.h:2301
INMOST::BlockOfMatrix::Rows
__INLINE enumerator Rows() const
Number of rows in submatrix.
Definition: inmost_dense.h:2283
INMOST::BlockOfMatrix::MakeMatrix
::INMOST::Matrix< Var > MakeMatrix()
Convert block of matrix into matrix.
Definition: inmost_dense.h:2373
INMOST::BlockOfMatrix::Cols
__INLINE enumerator Cols() const
Number of columns in submatrix.
Definition: inmost_dense.h:2286
INMOST::BlockOfMatrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2345
INMOST::BlockOfMatrix::BlockOfMatrix
BlockOfMatrix(AbstractMatrix< Var > &rM, enumerator num_rows, enumerator num_cols, enumerator offset_row, enumerator offset_col)
Create submatrix for a matrix.
Definition: inmost_dense.h:2293
INMOST::BlockOfMatrix::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:2384
INMOST::BlockOfMatrix::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2361
INMOST::ConstBlockOfMatrix
This class allows to address a matrix as a block of an empty matrix of larger size.
Definition: inmost_dense.h:2395
INMOST::ConstBlockOfMatrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2426
INMOST::ConstBlockOfMatrix::MakeMatrix
::INMOST::Matrix< Var > MakeMatrix()
Convert block of matrix into matrix.
Definition: inmost_dense.h:2440
INMOST::ConstBlockOfMatrix::Cols
__INLINE enumerator Cols() const
Number of columns in submatrix.
Definition: inmost_dense.h:2411
INMOST::ConstBlockOfMatrix::ConstBlockOfMatrix
ConstBlockOfMatrix(const AbstractMatrixReadOnly< Var > &rM, enumerator num_rows, enumerator num_cols, enumerator offset_row, enumerator offset_col)
Create submatrix for a matrix.
Definition: inmost_dense.h:2418
INMOST::ConstBlockOfMatrix::Rows
__INLINE enumerator Rows() const
Number of rows in submatrix.
Definition: inmost_dense.h:2408
INMOST::ConstMatrixConcatCols2::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2983
INMOST::ConstMatrixConcatCols2::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2969
INMOST::ConstMatrixConcatCols2::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2972
INMOST::ConstMatrixConcatCols2::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2990
INMOST::ConstMatrixConcatCols::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2935
INMOST::ConstMatrixConcatCols::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2953
INMOST::ConstMatrixConcatCols::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2949
INMOST::ConstMatrixConcatCols::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2938
INMOST::ConstMatrixConcatRows2::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2853
INMOST::ConstMatrixConcatRows2::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2842
INMOST::ConstMatrixConcatRows2::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2860
INMOST::ConstMatrixConcatRows2::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2839
INMOST::ConstMatrixConcatRows::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2819
INMOST::ConstMatrixConcatRows::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2808
INMOST::ConstMatrixConcatRows::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2805
INMOST::ConstMatrixConcatRows::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2823
INMOST::ConstMatrixConjugateTranspose::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2669
INMOST::ConstMatrixConjugateTranspose::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2678
INMOST::ConstMatrixConjugateTranspose::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2679
INMOST::ConstMatrixConjugateTranspose::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2666
INMOST::ConstMatrixConjugate::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2708
INMOST::ConstMatrixConjugate::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2699
INMOST::ConstMatrixConjugate::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2709
INMOST::ConstMatrixConjugate::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2696
INMOST::ConstMatrixRepack::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3078
INMOST::ConstMatrixRepack::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3069
INMOST::ConstMatrixRepack::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3066
INMOST::ConstMatrixRepack::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3083
INMOST::ConstMatrixTranspose::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2650
INMOST::ConstMatrixTranspose::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2641
INMOST::ConstMatrixTranspose::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2651
INMOST::ConstMatrixTranspose::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2638
INMOST::ConstSubMatrix
This class allows for in-place operations on submatrix of the matrix elements.
Definition: inmost_dense.h:2212
INMOST::ConstSubMatrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2244
INMOST::ConstSubMatrix::Rows
__INLINE enumerator Rows() const
Number of rows in submatrix.
Definition: inmost_dense.h:2225
INMOST::ConstSubMatrix::Cols
__INLINE enumerator Cols() const
Number of columns in submatrix.
Definition: inmost_dense.h:2228
INMOST::ConstSubMatrix::MakeMatrix
::INMOST::Matrix< Var > MakeMatrix()
Convert submatrix into matrix.
Definition: inmost_dense.h:2256
INMOST::ConstSubMatrix::ConstSubMatrix
ConstSubMatrix(const AbstractMatrixReadOnly< Var > &rM, enumerator first_row, enumerator last_row, enumerator first_column, enumerator last_column)
Create submatrix for a matrix.
Definition: inmost_dense.h:2235
INMOST::KroneckerProduct::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3646
INMOST::KroneckerProduct::TrivialArguments
bool TrivialArguments() const
Check that this matrix does not require calculations for element access.
Definition: inmost_dense.h:3640
INMOST::KroneckerProduct::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3643
INMOST::KroneckerProduct::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3659
INMOST::KroneckerProduct::compute
__INLINE Promote< VarA, VarB >::type compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3655
INMOST::MatrixCol::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2512
INMOST::MatrixCol::compute
__INLINE Var compute(enumerator i, enumerator j) const
Compute element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2516
INMOST::MatrixCol::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2515
INMOST::MatrixCol::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2517
INMOST::MatrixConcatCols::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3008
INMOST::MatrixConcatCols::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3038
INMOST::MatrixConcatCols::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3051
INMOST::MatrixConcatCols::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3029
INMOST::MatrixConcatCols::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3020
INMOST::MatrixConcatCols::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3011
INMOST::MatrixConcatCols::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3045
INMOST::MatrixConcatRows::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2880
INMOST::MatrixConcatRows::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2920
INMOST::MatrixConcatRows::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2877
INMOST::MatrixConcatRows::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2889
INMOST::MatrixConcatRows::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2907
INMOST::MatrixConcatRows::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2898
INMOST::MatrixConcatRows::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:2914
INMOST::MatrixDiag::compute
__INLINE Var compute(enumerator i, enumerator j) const
Compute element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2538
INMOST::MatrixDiag::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2534
INMOST::MatrixDiag::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2537
INMOST::MatrixDiag::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2542
INMOST::MatrixDifference::compute
__INLINE Promote< VarA, VarB >::type compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2618
INMOST::MatrixDifference::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2622
INMOST::MatrixDifference::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2605
INMOST::MatrixDifference::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2602
INMOST::MatrixDivCoef::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3552
INMOST::MatrixDivCoef::TrivialArguments
bool TrivialArguments() const
Check that this matrix does not require calculations for element access.
Definition: inmost_dense.h:3536
INMOST::MatrixDivCoef::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3542
INMOST::MatrixDivCoef::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3553
INMOST::MatrixDivCoef::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3539
INMOST::MatrixDivShellCoef::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3621
INMOST::MatrixDivShellCoef::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3622
INMOST::MatrixDivShellCoef::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3608
INMOST::MatrixDivShellCoef::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3611
INMOST::MatrixDivShellCoef::TrivialArguments
bool TrivialArguments() const
Check that this matrix does not require calculations for element access.
Definition: inmost_dense.h:3605
INMOST::MatrixMulCoef::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3505
INMOST::MatrixMulCoef::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3518
INMOST::MatrixMulCoef::TrivialArguments
bool TrivialArguments() const
Check that this matrix does not require calculations for element access.
Definition: inmost_dense.h:3502
INMOST::MatrixMulCoef::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3508
INMOST::MatrixMulCoef::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3517
INMOST::MatrixMulShellCoef::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3577
INMOST::MatrixMulShellCoef::TrivialArguments
bool TrivialArguments() const
Check that this matrix does not require calculations for element access.
Definition: inmost_dense.h:3571
INMOST::MatrixMulShellCoef::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3574
INMOST::MatrixMulShellCoef::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3587
INMOST::MatrixMulShellCoef::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3586
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3308
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::compute
__INLINE variable compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3305
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3245
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::operator()
__INLINE variable & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3306
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3237
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::get
__INLINE const variable & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3307
INMOST::MatrixMul< INMOST_DATA_REAL_TYPE, variable, Promote< INMOST_DATA_REAL_TYPE, variable >::type >::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3248
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3331
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3334
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3323
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::get
__INLINE const variable & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3393
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3394
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::operator()
__INLINE variable & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3392
INMOST::MatrixMul< variable, INMOST_DATA_REAL_TYPE, Promote< variable, INMOST_DATA_REAL_TYPE >::type >::compute
__INLINE variable compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3391
INMOST::MatrixMul< variable, variable, Promote< variable, variable >::type >::get
__INLINE const variable & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3482
INMOST::MatrixMul< variable, variable, Promote< variable, variable >::type >::operator()
__INLINE variable & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3481
INMOST::MatrixMul< variable, variable, Promote< variable, variable >::type >::compute
__INLINE variable compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3480
INMOST::MatrixMul< variable, variable, Promote< variable, variable >::type >::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3483
INMOST::MatrixMul< variable, variable, Promote< variable, variable >::type >::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3409
INMOST::MatrixMul::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3222
INMOST::MatrixMul::operator()
__INLINE VarR & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3215
INMOST::MatrixMul::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3176
INMOST::MatrixMul::get
__INLINE const VarR & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3221
INMOST::MatrixMul::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3173
INMOST::MatrixMul::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3165
INMOST::MatrixMul::compute
__INLINE VarR compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3210
INMOST::MatrixRepack::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3113
INMOST::MatrixRepack::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:3099
INMOST::MatrixRepack::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:3143
INMOST::MatrixRepack::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:3149
INMOST::MatrixRepack::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:3102
INMOST::MatrixRepack::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:3134
INMOST::MatrixRepack::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:3123
INMOST::MatrixRow::compute
__INLINE Var compute(enumerator i, enumerator j) const
Compute element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2494
INMOST::MatrixRow::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2495
INMOST::MatrixRow::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2493
INMOST::MatrixRow::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2490
INMOST::MatrixSum::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2584
INMOST::MatrixSum::compute
__INLINE Promote< VarA, VarB >::type compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2580
INMOST::MatrixSum::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2567
INMOST::MatrixSum::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2564
INMOST::MatrixTranspose::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2768
INMOST::MatrixTranspose::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2774
INMOST::MatrixTranspose::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:2784
INMOST::MatrixTranspose::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2780
INMOST::MatrixTranspose::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2790
INMOST::MatrixTranspose::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2756
INMOST::MatrixTranspose::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2759
INMOST::MatrixUnaryMinus::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2726
INMOST::MatrixUnaryMinus::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2738
INMOST::MatrixUnaryMinus::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2729
INMOST::MatrixUnaryMinus::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2739
INMOST::MatrixUnit::Rows
__INLINE enumerator Rows() const
Number of rows.
Definition: inmost_dense.h:2465
INMOST::MatrixUnit::compute
__INLINE Var compute(enumerator i, enumerator j) const
Compute element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2469
INMOST::MatrixUnit::GetMatrixCount
__INLINE INMOST_DATA_ENUM_TYPE GetMatrixCount() const
Retrieve number of indices of derivatives.
Definition: inmost_dense.h:2473
INMOST::MatrixUnit::Cols
__INLINE enumerator Cols() const
Number of columns.
Definition: inmost_dense.h:2468
INMOST::Matrix
Class for linear algebra operations on dense matrices.
Definition: inmost_dense.h:1434
INMOST::Matrix::Matrix
Matrix(const AbstractMatrix< typeB > &other)
Construct matrix from matrix of different type.
Definition: inmost_dense.h:1628
INMOST::Matrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:1707
INMOST::Matrix::Cols
__INLINE enumerator & Cols()
Obtain number of rows.
Definition: inmost_dense.h:1744
INMOST::Matrix::~Matrix
~Matrix()
Delete matrix.
Definition: inmost_dense.h:1635
INMOST::Matrix::RemoveRow
void RemoveRow(enumerator row)
Erase single row.
Definition: inmost_dense.h:1448
INMOST::Matrix::Permutation
static Matrix< Var > Permutation(const INMOST_DATA_ENUM_TYPE *Perm, enumerator size)
Construct row permutation matrix from array of new positions for rows.
Definition: inmost_dense.h:1754
INMOST::Matrix::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:1719
INMOST::Matrix::RemoveRows
void RemoveRows(enumerator first, enumerator last)
Erase multiple rows.
Definition: inmost_dense.h:1464
INMOST::Matrix::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:1696
INMOST::Matrix::Resize
void Resize(enumerator nrows, enumerator mcols)
Resize the matrix into different size.
Definition: inmost_dense.h:1639
INMOST::Matrix::Matrix
Matrix(enumerator pn, enumerator pm, const Var &c)
Construct a matrix with provided sizes and fills with value.
Definition: inmost_dense.h:1586
INMOST::Matrix::FromDiagonal
static Matrix< Var > FromDiagonal(const Var *r, enumerator size)
Create diagonal matrix from array.
Definition: inmost_dense.h:1909
INMOST::Matrix::Matrix
Matrix(const storage_type &pspace)
Construct the matrix with the provided storage and unknown size.
Definition: inmost_dense.h:1576
INMOST::Matrix::CrossProductMatrix
static Matrix< Var > CrossProductMatrix(const Var vec[3])
Cross-product matrix.
Definition: inmost_dense.h:1933
INMOST::Matrix::JointDiagonalization
Matrix< Var > JointDiagonalization(INMOST_DATA_REAL_TYPE threshold=1.0e-7)
Joint diagonalization algorithm by Cardoso.
Definition: inmost_dense.h:1990
INMOST::Matrix::Matrix
Matrix(const Matrix &other)
Copy matrix.
Definition: inmost_dense.h:1607
INMOST::Matrix::Cols
__INLINE enumerator Cols() const
Obtain number of columns.
Definition: inmost_dense.h:1738
INMOST::Matrix::FromVector
static Matrix< Var > FromVector(const Var *r, enumerator size)
Create column-vector in matrix form from array.
Definition: inmost_dense.h:1901
INMOST::Matrix::FromTensor
static Matrix< Var > FromTensor(const Var *K, enumerator size, enumerator matsize=3)
Convert values in array into square matrix.
Definition: inmost_dense.h:1783
INMOST::Matrix::RemoveColumns
void RemoveColumns(enumerator first, enumerator last)
Erase multiple columns.
Definition: inmost_dense.h:1498
INMOST::Matrix::Matrix
Matrix(enumerator pn, enumerator pm)
Construct a matrix with provided sizes.
Definition: inmost_dense.h:1581
INMOST::Matrix::Matrix
Matrix(const AbstractMatrixReadOnly< typeB > &other)
Construct matrix from matrix of different type.
Definition: inmost_dense.h:1617
INMOST::Matrix::Rows
__INLINE enumerator Rows() const
Obtain number of rows.
Definition: inmost_dense.h:1735
INMOST::Matrix::data
__INLINE const Var * data() const
Return raw pointer to matrix data without right of change, stored in row-wise format.
Definition: inmost_dense.h:1732
INMOST::Matrix::data
__INLINE Var * data()
Return raw pointer to matrix data, stored in row-wise format.
Definition: inmost_dense.h:1728
INMOST::Matrix::Row
static MatrixRow< Var > Row(enumerator pn, const Var &c=1.0)
Matix with 1 row, Create a matrix of size 1 by pn and fills it with c.
Definition: inmost_dense.h:1970
INMOST::Matrix::Matrix
Matrix(const storage_type &pspace, enumerator pn, enumerator pm)
Construct the matrix with the provided storage with known size.
Definition: inmost_dense.h:1570
INMOST::Matrix::operator=
Matrix & operator=(Matrix const &other)
Assign matrix of the same type.
Definition: inmost_dense.h:1649
INMOST::Matrix::Unit
static MatrixUnit< Var > Unit(enumerator pn, const Var &c=1.0)
Unit matrix.
Definition: inmost_dense.h:1955
INMOST::Matrix::Matrix
Matrix()
Construct empty matrix.
Definition: inmost_dense.h:1557
INMOST::Matrix::Matrix
Matrix(const Var *pspace, enumerator pn, enumerator pm)
Construct the matrix from provided array and sizes.
Definition: inmost_dense.h:1562
INMOST::Matrix::Swap
void Swap(AbstractMatrix< Var > &b)
Exchange contents of two matrices.
Definition: inmost_dense.h:1541
INMOST::Matrix::RemoveColumn
void RemoveColumn(enumerator col)
Erase single column.
Definition: inmost_dense.h:1480
INMOST::Matrix::FromDiagonalInverse
static Matrix< Var > FromDiagonalInverse(const Var *r, enumerator size)
Create diagonal matrix from array of values that have to be inversed.
Definition: inmost_dense.h:1920
INMOST::Matrix::Make
static Matrix Make(enumerator pn, enumerator pm,...)
Construct a matrix with provided elements.
Definition: inmost_dense.h:1591
INMOST::Matrix::RemoveSubset
void RemoveSubset(enumerator firstrow, enumerator lastrow, enumerator firstcol, enumerator lastcol)
Erase part of the matrix.
Definition: inmost_dense.h:1519
INMOST::Matrix::Col
static MatrixCol< Var > Col(enumerator pn, const Var &c=1.0)
Matix with 1 column, Create a matrix of size pn by 1 and fills it with c.
Definition: inmost_dense.h:1977
INMOST::Matrix::Rows
__INLINE enumerator & Rows()
Obtain number of rows.
Definition: inmost_dense.h:1741
INMOST::Sparse::RowMerger
This class may be used to sum multiple sparse rows.
Definition: inmost_sparse.h:638
INMOST::Sparse::RowMerger::AddRow
void AddRow(INMOST_DATA_REAL_TYPE coef, const Row &r)
Add a row with a coefficient into non-empty linked list.
INMOST::Sparse::RowMerger::RetrieveRow
void RetrieveRow(Row &r) const
Place entries from linked list into row.
INMOST::Sparse::RowMerger::Resize
void Resize(INMOST_DATA_ENUM_TYPE size)
Resize linked list for new interval.
INMOST::Sparse::RowMerger::Clear
void Clear()
Clear linked list.
INMOST::Sparse::RowMerger::PushRow
void PushRow(INMOST_DATA_REAL_TYPE coef, const Row &r)
Add a row with a coefficient into empty linked list.
INMOST::SubMatrix
This class allows for in-place operations on submatrix of the matrix elements.
Definition: inmost_dense.h:2083
INMOST::SubMatrix::operator=
SubMatrix & operator=(AbstractMatrixReadOnly< typeB > const &other)
Assign matrix of another type to submatrix.
Definition: inmost_dense.h:2115
INMOST::SubMatrix::Cols
__INLINE enumerator Cols() const
Number of columns in submatrix.
Definition: inmost_dense.h:2101
INMOST::SubMatrix::SubMatrix
SubMatrix(AbstractMatrix< Var > &rM, enumerator first_row, enumerator last_row, enumerator first_column, enumerator last_column)
Create submatrix for a matrix.
Definition: inmost_dense.h:2108
INMOST::SubMatrix::MakeMatrix
::INMOST::Matrix< Var > MakeMatrix()
Convert submatrix into matrix.
Definition: inmost_dense.h:2190
INMOST::SubMatrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2166
INMOST::SubMatrix::Resize
void Resize(enumerator rows, enumerator cols)
This is a stub function to fulfill abstract inheritance.
Definition: inmost_dense.h:2201
INMOST::SubMatrix::Rows
__INLINE enumerator Rows() const
Number of rows in submatrix.
Definition: inmost_dense.h:2098
INMOST::SubMatrix::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:2154
INMOST::SubMatrix::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:2178
INMOST::SymmetricMatrix::Cols
__INLINE enumerator Cols() const
Obtain number of columns.
Definition: inmost_dense.h:1232
INMOST::SymmetricMatrix::FromTensor
static SymmetricMatrix< Var > FromTensor(const Var *K, enumerator size, enumerator matsize=3)
Convert values in array into square matrix.
Definition: inmost_dense.h:1257
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const SymmetricMatrix &other)
Copy matrix.
Definition: inmost_dense.h:1089
INMOST::SymmetricMatrix::Rows
__INLINE enumerator Rows() const
Obtain number of rows.
Definition: inmost_dense.h:1229
INMOST::SymmetricMatrix::data
__INLINE Var * data()
Return raw pointer to matrix data, stored in row-wise format.
Definition: inmost_dense.h:1222
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const AbstractMatrixReadOnly< typeB > &other)
Construct matrix from matrix of different type.
Definition: inmost_dense.h:1101
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const storage_type &pspace, enumerator pn)
Construct the matrix with the provided storage with known size.
Definition: inmost_dense.h:1052
INMOST::SymmetricMatrix::Make
static SymmetricMatrix Make(enumerator pn,...)
Construct a matrix with provided elements.
Definition: inmost_dense.h:1068
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const Var *pspace, enumerator pn)
Construct the matrix from provided array and sizes.
Definition: inmost_dense.h:1044
INMOST::SymmetricMatrix::Rows
__INLINE enumerator & Rows()
Obtain number of rows.
Definition: inmost_dense.h:1235
INMOST::SymmetricMatrix::get
__INLINE const Var & get(enumerator i, enumerator j) const
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:1200
INMOST::SymmetricMatrix::Cols
__INLINE enumerator & Cols()
Obtain number of rows.
Definition: inmost_dense.h:1238
INMOST::SymmetricMatrix::Resize
void Resize(enumerator nrows, enumerator ncols)
Resize the matrix into different size.
Definition: inmost_dense.h:1128
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(enumerator pn, const Var &c)
Construct a matrix with provided sizes and fills with value.
Definition: inmost_dense.h:1086
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const storage_type &pspace)
Construct the matrix with the provided storage and unknown size.
Definition: inmost_dense.h:1058
INMOST::SymmetricMatrix::data
__INLINE const Var * data() const
Return raw pointer to matrix data without right of change, stored in row-wise format.
Definition: inmost_dense.h:1226
INMOST::SymmetricMatrix::operator()
__INLINE Var & operator()(enumerator i, enumerator j)
Access element of the matrix by row and column indices.
Definition: inmost_dense.h:1189
INMOST::SymmetricMatrix::FromDiagonal
static SymmetricMatrix< Var > FromDiagonal(const Var *r, enumerator size)
Create diagonal matrix from array.
Definition: inmost_dense.h:1353
INMOST::SymmetricMatrix::FromDiagonalInverse
static SymmetricMatrix< Var > FromDiagonalInverse(const Var *r, enumerator size)
Create diagonal matrix from array of values that have to be inversed.
Definition: inmost_dense.h:1364
INMOST::SymmetricMatrix::~SymmetricMatrix
~SymmetricMatrix()
Delete matrix.
Definition: inmost_dense.h:1123
INMOST::SymmetricMatrix::compute
__INLINE Var compute(enumerator i, enumerator j) const
Access element of the matrix by row and column indices without right to change the element.
Definition: inmost_dense.h:1212
INMOST::SymmetricMatrix::Unit
static SymmetricMatrix< Var > Unit(enumerator pn, const Var &c=1.0)
Unit matrix.
Definition: inmost_dense.h:1376
INMOST::SymmetricMatrix::Swap
void Swap(AbstractMatrix< Var > &b)
Exchange contents of two matrices.
Definition: inmost_dense.h:1024
INMOST::SymmetricMatrix::operator=
SymmetricMatrix & operator=(SymmetricMatrix const &other)
Assign matrix of the same type.
Definition: inmost_dense.h:1139
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix()
Construct empty matrix.
Definition: inmost_dense.h:1039
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(enumerator pn)
Construct a matrix with provided sizes.
Definition: inmost_dense.h:1063
INMOST::SymmetricMatrix::SymmetricMatrix
SymmetricMatrix(const AbstractMatrix< typeB > &other)
Construct matrix from matrix of different type.
Definition: inmost_dense.h:1115
INMOST::hessian_multivar_expression_reference
Definition: inmost_expression.h:956
INMOST::hessian_multivar_expression
A class that represents a variable with multiple first order and second order variations.
Definition: inmost_expression.h:408
INMOST::multivar_expression_reference
Definition: inmost_expression.h:642
INMOST::multivar_expression
A class that represents a variable with multiple first order variations.
Definition: inmost_expression.h:131
INMOST::shell_expression
Definition: inmost_expression.h:52
INMOST::thread_private
Definition: container.hpp:1569
INMOST::value_reference
Definition: inmost_expression.h:871
INMOST::var_expression
Definition: inmost_expression.h:86
INMOST::SelfPromote
Structure that selects desired class, depending on the operation.
Definition: inmost_dense.h:12
INMOST::make_integer
Definition: inmost_dense.h:4894