GNU Octave  8.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
dMatrix.h
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1994-2023 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if ! defined (octave_dMatrix_h)
27 #define octave_dMatrix_h 1
28 
29 #include "octave-config.h"
30 
31 #include "DET.h"
32 #include "MArray.h"
33 #include "MDiagArray2.h"
34 #include "MatrixType.h"
35 #include "dNDArray.h"
36 #include "mx-defs.h"
37 #include "mx-op-decl.h"
38 
39 class
41 Matrix : public NDArray
42 {
43 public:
44 
47 
50 
53 
56 
57  typedef double real_elt_type;
59 
60  typedef void (*solve_singularity_handler) (double rcon);
61 
62  Matrix (void) = default;
63 
64  Matrix (const Matrix& a) = default;
65 
66  Matrix& operator = (const Matrix& a) = default;
67 
68  ~Matrix (void) = default;
69 
71  : NDArray (dim_vector (r, c)) { }
72 
74  : NDArray (dim_vector (r, c), val) { }
75 
76  Matrix (const dim_vector& dv) : NDArray (dv.redim (2)) { }
77 
78  Matrix (const dim_vector& dv, double val)
79  : NDArray (dv.redim (2), val) { }
80 
81  template <typename U>
82  Matrix (const MArray<U>& a) : NDArray (a.as_matrix ()) { }
83 
84  template <typename U>
85  Matrix (const Array<U>& a) : NDArray (a.as_matrix ()) { }
86 
87  explicit OCTAVE_API Matrix (const RowVector& rv);
88 
89  explicit OCTAVE_API Matrix (const ColumnVector& cv);
90 
91  explicit OCTAVE_API Matrix (const DiagMatrix& a);
92 
93  explicit OCTAVE_API Matrix (const MDiagArray2<double>& a);
94 
95  explicit OCTAVE_API Matrix (const DiagArray2<double>& a);
96 
97  explicit OCTAVE_API Matrix (const PermMatrix& a);
98 
99  explicit OCTAVE_API Matrix (const boolMatrix& a);
100 
101  explicit Matrix (const charMatrix& a);
102 
103  OCTAVE_API bool operator == (const Matrix& a) const;
104  OCTAVE_API bool operator != (const Matrix& a) const;
105 
106  OCTAVE_API bool issymmetric (void) const;
107 
108  // destructive insert/delete/reorder operations
109 
111  insert (const Matrix& a, octave_idx_type r, octave_idx_type c);
113  insert (const RowVector& a, octave_idx_type r, octave_idx_type c);
115  insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c);
117  insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c);
118 
119  OCTAVE_API Matrix& fill (double val);
121  fill (double val, octave_idx_type r1, octave_idx_type c1,
123 
124  OCTAVE_API Matrix append (const Matrix& a) const;
125  OCTAVE_API Matrix append (const RowVector& a) const;
126  OCTAVE_API Matrix append (const ColumnVector& a) const;
127  OCTAVE_API Matrix append (const DiagMatrix& a) const;
128 
129  OCTAVE_API Matrix stack (const Matrix& a) const;
130  OCTAVE_API Matrix stack (const RowVector& a) const;
131  OCTAVE_API Matrix stack (const ColumnVector& a) const;
132  OCTAVE_API Matrix stack (const DiagMatrix& a) const;
133 
134  friend OCTAVE_API Matrix real (const ComplexMatrix& a);
135  friend OCTAVE_API Matrix imag (const ComplexMatrix& a);
136 
137  friend class ComplexMatrix;
138 
139  Matrix hermitian (void) const { return MArray<double>::transpose (); }
140  Matrix transpose (void) const { return MArray<double>::transpose (); }
141 
142  // resize is the destructive equivalent for this one
143 
146  octave_idx_type r2, octave_idx_type c2) const;
147 
150  octave_idx_type nr, octave_idx_type nc) const;
151 
152  // extract row or column i.
153 
155 
157 
158  void resize (octave_idx_type nr, octave_idx_type nc, double rfv = 0)
159  {
160  MArray<double>::resize (dim_vector (nr, nc), rfv);
161  }
162 
163 private:
164  Matrix tinverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
165  bool force, bool calc_cond) const;
166 
167  Matrix finverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
168  bool force, bool calc_cond) const;
169 
170 public:
171  OCTAVE_API Matrix inverse (void) const;
174  inverse (octave_idx_type& info, double& rcon, bool force = false,
175  bool calc_cond = true) const;
176 
177  OCTAVE_API Matrix inverse (MatrixType& mattype) const;
178  OCTAVE_API Matrix inverse (MatrixType& mattype, octave_idx_type& info) const;
180  inverse (MatrixType& mattype, octave_idx_type& info, double& rcon,
181  bool force = false, bool calc_cond = true) const;
182 
183  OCTAVE_API Matrix pseudo_inverse (double tol = 0.0) const;
184 
185  OCTAVE_API ComplexMatrix fourier (void) const;
186  OCTAVE_API ComplexMatrix ifourier (void) const;
187 
188  OCTAVE_API ComplexMatrix fourier2d (void) const;
189  OCTAVE_API ComplexMatrix ifourier2d (void) const;
190 
191  OCTAVE_API DET determinant (void) const;
194  determinant (octave_idx_type& info, double& rcon,
195  bool calc_cond = true) const;
197  determinant (MatrixType& mattype, octave_idx_type& info,
198  double& rcon, bool calc_cond = true) const;
199 
200  OCTAVE_API double rcond (void) const;
201  OCTAVE_API double rcond (MatrixType& mattype) const;
202 
203 private:
204  // Upper triangular matrix solvers
205  Matrix utsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
206  double& rcon, solve_singularity_handler sing_handler,
207  bool calc_cond = false,
208  blas_trans_type transt = blas_no_trans) const;
209 
210  // Lower triangular matrix solvers
211  Matrix ltsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
212  double& rcon, solve_singularity_handler sing_handler,
213  bool calc_cond = false,
214  blas_trans_type transt = blas_no_trans) const;
215 
216  // Full matrix solvers (lu/cholesky)
217  Matrix fsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
218  double& rcon, solve_singularity_handler sing_handler,
219  bool calc_cond = false) const;
220 
221 public:
222  // Generic interface to solver with no probing of type
223  OCTAVE_API Matrix solve (MatrixType& mattype, const Matrix& b) const;
225  solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info) const;
227  solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
228  double& rcon) const;
230  solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info,
231  double& rcon, solve_singularity_handler sing_handler,
232  bool singular_fallback = true,
233  blas_trans_type transt = blas_no_trans) const;
234 
236  solve (MatrixType& mattype, const ComplexMatrix& b) const;
238  solve (MatrixType& mattype, const ComplexMatrix& b,
239  octave_idx_type& info) const;
241  solve (MatrixType& mattype, const ComplexMatrix& b,
242  octave_idx_type& info, double& rcon) const;
244  solve (MatrixType& mattype, const ComplexMatrix& b,
245  octave_idx_type& info, double& rcon,
246  solve_singularity_handler sing_handler,
247  bool singular_fallback = true,
248  blas_trans_type transt = blas_no_trans) const;
249 
251  solve (MatrixType& mattype, const ColumnVector& b) const;
253  solve (MatrixType& mattype, const ColumnVector& b,
254  octave_idx_type& info) const;
256  solve (MatrixType& mattype, const ColumnVector& b,
257  octave_idx_type& info, double& rcon) const;
259  solve (MatrixType& mattype, const ColumnVector& b, octave_idx_type& info,
260  double& rcon, solve_singularity_handler sing_handler,
261  blas_trans_type transt = blas_no_trans) const;
262 
264  solve (MatrixType& mattype, const ComplexColumnVector& b) const;
266  solve (MatrixType& mattype, const ComplexColumnVector& b,
267  octave_idx_type& info) const;
269  solve (MatrixType& mattype, const ComplexColumnVector& b,
270  octave_idx_type& info, double& rcon) const;
272  solve (MatrixType& mattype, const ComplexColumnVector& b,
273  octave_idx_type& info, double& rcon,
274  solve_singularity_handler sing_handler,
275  blas_trans_type transt = blas_no_trans) const;
276 
277  // Generic interface to solver with probing of type
278  OCTAVE_API Matrix solve (const Matrix& b) const;
279  OCTAVE_API Matrix solve (const Matrix& b, octave_idx_type& info) const;
281  solve (const Matrix& b, octave_idx_type& info, double& rcon) const;
283  solve (const Matrix& b, octave_idx_type& info, double& rcon,
284  solve_singularity_handler sing_handler,
285  blas_trans_type transt = blas_no_trans) const;
286 
287  OCTAVE_API ComplexMatrix solve (const ComplexMatrix& b) const;
289  solve (const ComplexMatrix& b, octave_idx_type& info) const;
291  solve (const ComplexMatrix& b, octave_idx_type& info,
292  double& rcon) const;
294  solve (const ComplexMatrix& b, octave_idx_type& info, double& rcon,
295  solve_singularity_handler sing_handler,
296  blas_trans_type transt = blas_no_trans) const;
297 
298  OCTAVE_API ColumnVector solve (const ColumnVector& b) const;
300  solve (const ColumnVector& b, octave_idx_type& info) const;
302  solve (const ColumnVector& b, octave_idx_type& info, double& rcon) const;
304  solve (const ColumnVector& b, octave_idx_type& info, double& rcon,
305  solve_singularity_handler sing_handler,
306  blas_trans_type transt = blas_no_trans) const;
307 
310  solve (const ComplexColumnVector& b, octave_idx_type& info) const;
312  solve (const ComplexColumnVector& b, octave_idx_type& info,
313  double& rcon) const;
315  solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcon,
316  solve_singularity_handler sing_handler,
317  blas_trans_type transt = blas_no_trans) const;
318 
319  // Singular solvers
320  OCTAVE_API Matrix lssolve (const Matrix& b) const;
321  OCTAVE_API Matrix lssolve (const Matrix& b, octave_idx_type& info) const;
323  lssolve (const Matrix& b, octave_idx_type& info,
324  octave_idx_type& rank) const;
326  lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank,
327  double& rcon) const;
328 
331  lssolve (const ComplexMatrix& b, octave_idx_type& info) const;
333  lssolve (const ComplexMatrix& b, octave_idx_type& info,
334  octave_idx_type& rank) const;
336  lssolve (const ComplexMatrix& b, octave_idx_type& info,
337  octave_idx_type& rank, double& rcon) const;
338 
339  OCTAVE_API ColumnVector lssolve (const ColumnVector& b) const;
341  lssolve (const ColumnVector& b, octave_idx_type& info) const;
343  lssolve (const ColumnVector& b, octave_idx_type& info,
344  octave_idx_type& rank) const;
346  lssolve (const ColumnVector& b, octave_idx_type& info,
347  octave_idx_type& rank, double& rcon) const;
348 
351  lssolve (const ComplexColumnVector& b, octave_idx_type& info) const;
354  octave_idx_type& rank) const;
357  octave_idx_type& rank, double& rcon) const;
358 
361 
362  // unary operations
363 
364  // other operations
365 
366  OCTAVE_API boolMatrix all (int dim = -1) const;
367  OCTAVE_API boolMatrix any (int dim = -1) const;
368 
369  OCTAVE_API Matrix cumprod (int dim = -1) const;
370  OCTAVE_API Matrix cumsum (int dim = -1) const;
371  OCTAVE_API Matrix prod (int dim = -1) const;
372  OCTAVE_API Matrix sum (int dim = -1) const;
373  OCTAVE_API Matrix sumsq (int dim = -1) const;
374  OCTAVE_API Matrix abs (void) const;
375 
376  OCTAVE_API Matrix diag (octave_idx_type k = 0) const;
377 
379 
380  OCTAVE_API ColumnVector row_min (void) const;
381  OCTAVE_API ColumnVector row_max (void) const;
382 
385 
386  OCTAVE_API RowVector column_min (void) const;
387  OCTAVE_API RowVector column_max (void) const;
388 
391 
392  // i/o
393 
394  friend OCTAVE_API std::ostream&
395  operator << (std::ostream& os, const Matrix& a);
396  friend OCTAVE_API std::istream&
397  operator >> (std::istream& is, Matrix& a);
398 };
399 
400 // Publish externally used friend functions.
401 
402 extern OCTAVE_API Matrix real (const ComplexMatrix& a);
403 extern OCTAVE_API Matrix imag (const ComplexMatrix& a);
404 
405 // column vector by row vector -> matrix operations
406 
407 extern OCTAVE_API Matrix operator * (const ColumnVector& a,
408  const RowVector& b);
409 
410 // Other functions.
411 
412 extern OCTAVE_API Matrix Givens (double, double);
413 
414 extern OCTAVE_API Matrix Sylvester (const Matrix&, const Matrix&,
415  const Matrix&);
416 
417 extern OCTAVE_API Matrix xgemm (const Matrix& a, const Matrix& b,
419  blas_trans_type transb = blas_no_trans);
420 
421 extern OCTAVE_API Matrix operator * (const Matrix& a, const Matrix& b);
422 
423 extern OCTAVE_API Matrix min (double d, const Matrix& m);
424 extern OCTAVE_API Matrix min (const Matrix& m, double d);
425 extern OCTAVE_API Matrix min (const Matrix& a, const Matrix& b);
426 
427 extern OCTAVE_API Matrix max (double d, const Matrix& m);
428 extern OCTAVE_API Matrix max (const Matrix& m, double d);
429 extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b);
430 
431 extern OCTAVE_API Matrix linspace (const ColumnVector& x1,
432  const ColumnVector& x2,
434 
437 
440 
443 
445 
446 template <typename T>
447 void read_int (std::istream& is, bool swap_bytes, T& val);
448 
449 #endif
template OCTAVE_API std::ostream & operator<<(std::ostream &, const Array< bool > &)
#define MARRAY_FORWARD_DEFS(B, R, T)
Definition: MArray.h:130
std::istream & operator>>(std::istream &is, SparseBoolMatrix &a)
Definition: boolSparse.cc:279
static void swap_bytes(void *ptr, unsigned int i, unsigned int j)
Definition: byte-swap.h:32
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:130
OCTARRAY_API Array< T, Alloc > index(const octave::idx_vector &i) const
Indexing without resizing.
Definition: Array-base.cc:719
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Definition: Array-base.cc:1032
ComplexMatrix fsolve(MatrixType &mattype, const ComplexMatrix &b, octave_idx_type &info, double &rcon, solve_singularity_handler sing_handler, bool calc_cond=false) const
Definition: CMatrix.cc:1722
OCTAVE_API ComplexMatrix lssolve(const Matrix &b) const
Definition: CMatrix.cc:2220
OCTAVE_API ComplexMatrix extract(octave_idx_type r1, octave_idx_type c1, octave_idx_type r2, octave_idx_type c2) const
Definition: CMatrix.cc:683
OCTAVE_API ComplexColumnVector row_max(void) const
Definition: CMatrix.cc:2961
OCTAVE_API boolMatrix all(int dim=-1) const
Definition: CMatrix.cc:2782
OCTAVE_API ComplexMatrix fourier2d(void) const
Definition: CMatrix.cc:1103
OCTAVE_API ComplexRowVector column_min(void) const
Definition: CMatrix.cc:3036
OCTAVE_API ComplexColumnVector row_min(void) const
Definition: CMatrix.cc:2886
ComplexMatrix finverse(MatrixType &mattype, octave_idx_type &info, double &rcon, bool force, bool calc_cond) const
Definition: CMatrix.cc:836
OCTAVE_API double rcond(void) const
Definition: CMatrix.cc:1350
OCTAVE_API ComplexMatrix diag(octave_idx_type k=0) const
Definition: CMatrix.cc:2830
OCTAVE_API ComplexMatrix pseudo_inverse(double tol=0.0) const
Definition: CMatrix.cc:1000
OCTAVE_API ComplexMatrix inverse(void) const
Definition: CMatrix.cc:737
OCTAVE_API ComplexDET determinant(void) const
Definition: CMatrix.cc:1171
OCTAVE_API ComplexMatrix ifourier(void) const
Definition: CMatrix.cc:1074
OCTAVE_API ComplexMatrix ifourier2d(void) const
Definition: CMatrix.cc:1117
OCTAVE_API ComplexRowVector row(octave_idx_type i) const
Definition: CMatrix.cc:702
OCTAVE_API ComplexMatrix solve(MatrixType &mattype, const Matrix &b) const
Definition: CMatrix.cc:1931
OCTAVE_API ComplexRowVector column_max(void) const
Definition: CMatrix.cc:3111
OCTAVE_API ComplexMatrix sumsq(int dim=-1) const
Definition: CMatrix.cc:2818
void(* solve_singularity_handler)(double rcon)
Definition: CMatrix.h:61
OCTAVE_API ComplexMatrix extract_n(octave_idx_type r1, octave_idx_type c1, octave_idx_type nr, octave_idx_type nc) const
Definition: CMatrix.cc:693
OCTAVE_API ComplexMatrix cumprod(int dim=-1) const
Definition: CMatrix.cc:2794
OCTAVE_API ComplexMatrix fourier(void) const
Definition: CMatrix.cc:1045
OCTAVE_API ComplexMatrix prod(int dim=-1) const
Definition: CMatrix.cc:2806
OCTAVE_API ComplexMatrix cumsum(int dim=-1) const
Definition: CMatrix.cc:2800
ComplexMatrix tinverse(MatrixType &mattype, octave_idx_type &info, double &rcon, bool force, bool calc_cond) const
Definition: CMatrix.cc:777
OCTAVE_API boolMatrix any(int dim=-1) const
Definition: CMatrix.cc:2788
OCTAVE_API ComplexMatrix sum(int dim=-1) const
Definition: CMatrix.cc:2812
OCTAVE_API ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:708
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:63
MArray< T > transpose(void) const
Definition: MArray.h:99
Definition: dMatrix.h:42
Matrix(const Array< U > &a)
Definition: dMatrix.h:85
DiagMatrix real_diag_matrix_type
Definition: dMatrix.h:54
RowVector row_vector_type
Definition: dMatrix.h:46
Matrix hermitian(void) const
Definition: dMatrix.h:139
ComplexDiagMatrix complex_diag_matrix_type
Definition: dMatrix.h:55
Matrix real_matrix_type
Definition: dMatrix.h:51
ComplexMatrix complex_matrix_type
Definition: dMatrix.h:52
RowVector real_row_vector_type
Definition: dMatrix.h:49
Complex complex_elt_type
Definition: dMatrix.h:58
Matrix(octave_idx_type r, octave_idx_type c)
Definition: dMatrix.h:70
Matrix(const Matrix &a)=default
Matrix(void)=default
Matrix transpose(void) const
Definition: dMatrix.h:140
Matrix(octave_idx_type r, octave_idx_type c, double val)
Definition: dMatrix.h:73
Matrix(const dim_vector &dv)
Definition: dMatrix.h:76
~Matrix(void)=default
ColumnVector real_column_vector_type
Definition: dMatrix.h:48
double real_elt_type
Definition: dMatrix.h:57
ColumnVector column_vector_type
Definition: dMatrix.h:45
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:158
Matrix(const MArray< U > &a)
Definition: dMatrix.h:82
Matrix(const dim_vector &dv, double val)
Definition: dMatrix.h:78
Definition: DET.h:39
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
OCTAVE_API Matrix min(double d, const Matrix &m)
Definition: dMatrix.cc:2884
Matrix operator+=(Matrix &x, const double &y)
Definition: dMatrix.h:444
void read_int(std::istream &is, bool swap_bytes, T &val)
OCTAVE_API Matrix real(const ComplexMatrix &a)
Definition: dMatrix.cc:385
OCTAVE_API Matrix Sylvester(const Matrix &, const Matrix &, const Matrix &)
Definition: dMatrix.cc:2704
OCTAVE_API Matrix max(double d, const Matrix &m)
Definition: dMatrix.cc:2948
OCTAVE_API Matrix imag(const ComplexMatrix &a)
Definition: dMatrix.cc:391
OCTAVE_API Matrix linspace(const ColumnVector &x1, const ColumnVector &x2, octave_idx_type n)
Definition: dMatrix.cc:3011
OCTAVE_API Matrix operator*(const ColumnVector &a, const RowVector &b)
Definition: dMatrix.cc:2330
Matrix operator-=(Matrix &x, const double &y)
Definition: dMatrix.h:444
OCTAVE_API Matrix xgemm(const Matrix &a, const Matrix &b, blas_trans_type transa=blas_no_trans, blas_trans_type transb=blas_no_trans)
Definition: dMatrix.cc:2786
OCTAVE_API Matrix Givens(double, double)
Definition: dMatrix.cc:2687
bool operator!=(const dim_vector &a, const dim_vector &b)
Definition: dim-vector.h:536
bool operator==(const dim_vector &a, const dim_vector &b)
Definition: dim-vector.h:520
static M utsolve(const SM &U, const ColumnVector &Q, const M &m)
Definition: eigs-base.cc:259
static M ltsolve(const SM &L, const ColumnVector &Q, const M &m)
Definition: eigs-base.cc:239
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE * d
#define OCTAVE_API
Definition: main.in.cc:55
blas_trans_type
Definition: mx-defs.h:80
@ blas_no_trans
Definition: mx-defs.h:81
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
T octave_idx_type m
Definition: mx-inlines.cc:773
octave_idx_type n
Definition: mx-inlines.cc:753
T * r
Definition: mx-inlines.cc:773
#define SM_BOOL_OP_DECLS(S, M, API)
Definition: mx-op-decl.h:122
#define MS_CMP_OP_DECLS(M, S, API)
Definition: mx-op-decl.h:89
#define SM_CMP_OP_DECLS(S, M, API)
Definition: mx-op-decl.h:114
#define MS_BOOL_OP_DECLS(M, S, API)
Definition: mx-op-decl.h:97
#define MM_CMP_OP_DECLS(M1, M2, API)
Definition: mx-op-decl.h:139
#define MM_BOOL_OP_DECLS(M1, M2, API)
Definition: mx-op-decl.h:147
std::complex< double > Complex
Definition: oct-cmplx.h:33
static T abs(T x)
Definition: pr-output.cc:1678