1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1994-2021 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 40 OCTAVE_API 41 Matrix : public NDArray 42 { 43 public: 44 45 typedef ColumnVector column_vector_type; 46 typedef RowVector row_vector_type; 47 48 typedef ColumnVector real_column_vector_type; 49 typedef RowVector real_row_vector_type; 50 51 typedef Matrix real_matrix_type; 52 typedef ComplexMatrix complex_matrix_type; 53 54 typedef DiagMatrix real_diag_matrix_type; 55 typedef ComplexDiagMatrix complex_diag_matrix_type; 56 57 typedef double real_elt_type; 58 typedef Complex complex_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 Matrix(octave_idx_type r,octave_idx_type c)70 Matrix (octave_idx_type r, octave_idx_type c) 71 : NDArray (dim_vector (r, c)) { } 72 Matrix(octave_idx_type r,octave_idx_type c,double val)73 Matrix (octave_idx_type r, octave_idx_type c, double val) 74 : NDArray (dim_vector (r, c), val) { } 75 Matrix(const dim_vector & dv)76 Matrix (const dim_vector& dv) : NDArray (dv.redim (2)) { } 77 Matrix(const dim_vector & dv,double val)78 Matrix (const dim_vector& dv, double val) 79 : NDArray (dv.redim (2), val) { } 80 81 template <typename U> Matrix(const MArray<U> & a)82 Matrix (const MArray<U>& a) : NDArray (a.as_matrix ()) { } 83 84 template <typename U> Matrix(const Array<U> & a)85 Matrix (const Array<U>& a) : NDArray (a.as_matrix ()) { } 86 87 explicit Matrix (const RowVector& rv); 88 89 explicit Matrix (const ColumnVector& cv); 90 91 explicit Matrix (const DiagMatrix& a); 92 93 explicit Matrix (const MDiagArray2<double>& a); 94 95 explicit Matrix (const DiagArray2<double>& a); 96 97 explicit Matrix (const PermMatrix& a); 98 99 explicit Matrix (const boolMatrix& a); 100 101 explicit Matrix (const charMatrix& a); 102 103 bool operator == (const Matrix& a) const; 104 bool operator != (const Matrix& a) const; 105 106 bool issymmetric (void) const; 107 108 // destructive insert/delete/reorder operations 109 110 Matrix& insert (const Matrix& a, octave_idx_type r, octave_idx_type c); 111 Matrix& insert (const RowVector& a, octave_idx_type r, octave_idx_type c); 112 Matrix& insert (const ColumnVector& a, octave_idx_type r, octave_idx_type c); 113 Matrix& insert (const DiagMatrix& a, octave_idx_type r, octave_idx_type c); 114 115 Matrix& fill (double val); 116 Matrix& fill (double val, octave_idx_type r1, octave_idx_type c1, 117 octave_idx_type r2, octave_idx_type c2); 118 119 Matrix append (const Matrix& a) const; 120 Matrix append (const RowVector& a) const; 121 Matrix append (const ColumnVector& a) const; 122 Matrix append (const DiagMatrix& a) const; 123 124 Matrix stack (const Matrix& a) const; 125 Matrix stack (const RowVector& a) const; 126 Matrix stack (const ColumnVector& a) const; 127 Matrix stack (const DiagMatrix& a) const; 128 129 friend OCTAVE_API Matrix real (const ComplexMatrix& a); 130 friend OCTAVE_API Matrix imag (const ComplexMatrix& a); 131 132 friend class ComplexMatrix; 133 hermitian(void)134 Matrix hermitian (void) const { return MArray<double>::transpose (); } transpose(void)135 Matrix transpose (void) const { return MArray<double>::transpose (); } 136 137 // resize is the destructive equivalent for this one 138 139 Matrix extract (octave_idx_type r1, octave_idx_type c1, 140 octave_idx_type r2, octave_idx_type c2) const; 141 142 Matrix extract_n (octave_idx_type r1, octave_idx_type c1, 143 octave_idx_type nr, octave_idx_type nc) const; 144 145 // extract row or column i. 146 147 RowVector row (octave_idx_type i) const; 148 149 ColumnVector column (octave_idx_type i) const; 150 151 void resize (octave_idx_type nr, octave_idx_type nc, double rfv = 0) 152 { 153 MArray<double>::resize (dim_vector (nr, nc), rfv); 154 } 155 156 private: 157 Matrix tinverse (MatrixType& mattype, octave_idx_type& info, double& rcon, 158 bool force, bool calc_cond) const; 159 160 Matrix finverse (MatrixType& mattype, octave_idx_type& info, double& rcon, 161 bool force, bool calc_cond) const; 162 163 public: 164 Matrix inverse (void) const; 165 Matrix inverse (octave_idx_type& info) const; 166 Matrix inverse (octave_idx_type& info, double& rcon, bool force = false, 167 bool calc_cond = true) const; 168 169 Matrix inverse (MatrixType& mattype) const; 170 Matrix inverse (MatrixType& mattype, octave_idx_type& info) const; 171 Matrix inverse (MatrixType& mattype, octave_idx_type& info, double& rcon, 172 bool force = false, bool calc_cond = true) const; 173 174 Matrix pseudo_inverse (double tol = 0.0) const; 175 176 ComplexMatrix fourier (void) const; 177 ComplexMatrix ifourier (void) const; 178 179 ComplexMatrix fourier2d (void) const; 180 ComplexMatrix ifourier2d (void) const; 181 182 DET determinant (void) const; 183 DET determinant (octave_idx_type& info) const; 184 DET determinant (octave_idx_type& info, double& rcon, 185 bool calc_cond = true) const; 186 DET determinant (MatrixType& mattype, octave_idx_type& info, 187 double& rcon, bool calc_cond = true) const; 188 189 double rcond (void) const; 190 double rcond (MatrixType& mattype) const; 191 192 private: 193 // Upper triangular matrix solvers 194 Matrix utsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info, 195 double& rcon, solve_singularity_handler sing_handler, 196 bool calc_cond = false, 197 blas_trans_type transt = blas_no_trans) const; 198 199 // Lower triangular matrix solvers 200 Matrix ltsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info, 201 double& rcon, solve_singularity_handler sing_handler, 202 bool calc_cond = false, 203 blas_trans_type transt = blas_no_trans) const; 204 205 // Full matrix solvers (lu/cholesky) 206 Matrix fsolve (MatrixType& mattype, const Matrix& b, octave_idx_type& info, 207 double& rcon, solve_singularity_handler sing_handler, 208 bool calc_cond = false) const; 209 210 public: 211 // Generic interface to solver with no probing of type 212 Matrix solve (MatrixType& mattype, const Matrix& b) const; 213 Matrix solve (MatrixType& mattype, const Matrix& b, 214 octave_idx_type& info) const; 215 Matrix solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info, 216 double& rcon) const; 217 Matrix solve (MatrixType& mattype, const Matrix& b, octave_idx_type& info, 218 double& rcon, solve_singularity_handler sing_handler, 219 bool singular_fallback = true, 220 blas_trans_type transt = blas_no_trans) const; 221 222 ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b) const; 223 ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b, 224 octave_idx_type& info) const; 225 ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b, 226 octave_idx_type& info, double& rcon) const; 227 ComplexMatrix solve (MatrixType& mattype, const ComplexMatrix& b, 228 octave_idx_type& info, double& rcon, 229 solve_singularity_handler sing_handler, 230 bool singular_fallback = true, 231 blas_trans_type transt = blas_no_trans) const; 232 233 ColumnVector solve (MatrixType& mattype, const ColumnVector& b) const; 234 ColumnVector solve (MatrixType& mattype, const ColumnVector& b, 235 octave_idx_type& info) const; 236 ColumnVector solve (MatrixType& mattype, const ColumnVector& b, 237 octave_idx_type& info, double& rcon) const; 238 ColumnVector solve (MatrixType& mattype, const ColumnVector& b, 239 octave_idx_type& info, double& rcon, 240 solve_singularity_handler sing_handler, 241 blas_trans_type transt = blas_no_trans) const; 242 243 ComplexColumnVector solve (MatrixType& mattype, 244 const ComplexColumnVector& b) const; 245 ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b, 246 octave_idx_type& info) const; 247 ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b, 248 octave_idx_type& info, double& rcon) const; 249 ComplexColumnVector solve (MatrixType& mattype, const ComplexColumnVector& b, 250 octave_idx_type& info, double& rcon, 251 solve_singularity_handler sing_handler, 252 blas_trans_type transt = blas_no_trans) const; 253 254 // Generic interface to solver with probing of type 255 Matrix solve (const Matrix& b) const; 256 Matrix solve (const Matrix& b, octave_idx_type& info) const; 257 Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon) const; 258 Matrix solve (const Matrix& b, octave_idx_type& info, double& rcon, 259 solve_singularity_handler sing_handler, 260 blas_trans_type transt = blas_no_trans) const; 261 262 ComplexMatrix solve (const ComplexMatrix& b) const; 263 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const; 264 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, 265 double& rcon) const; 266 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, 267 double& rcon, 268 solve_singularity_handler sing_handler, 269 blas_trans_type transt = blas_no_trans) const; 270 271 ColumnVector solve (const ColumnVector& b) const; 272 ColumnVector solve (const ColumnVector& b, octave_idx_type& info) const; 273 ColumnVector solve (const ColumnVector& b, octave_idx_type& info, 274 double& rcon) const; 275 ColumnVector solve (const ColumnVector& b, octave_idx_type& info, 276 double& rcon, 277 solve_singularity_handler sing_handler, 278 blas_trans_type transt = blas_no_trans) const; 279 280 ComplexColumnVector solve (const ComplexColumnVector& b) const; 281 ComplexColumnVector solve (const ComplexColumnVector& b, 282 octave_idx_type& info) const; 283 ComplexColumnVector solve (const ComplexColumnVector& b, 284 octave_idx_type& info, double& rcon) const; 285 ComplexColumnVector solve (const ComplexColumnVector& b, 286 octave_idx_type& info, double& rcon, 287 solve_singularity_handler sing_handler, 288 blas_trans_type transt = blas_no_trans) const; 289 290 // Singular solvers 291 Matrix lssolve (const Matrix& b) const; 292 Matrix lssolve (const Matrix& b, octave_idx_type& info) const; 293 Matrix lssolve (const Matrix& b, octave_idx_type& info, 294 octave_idx_type& rank) const; 295 Matrix lssolve (const Matrix& b, octave_idx_type& info, 296 octave_idx_type& rank, double& rcon) const; 297 298 ComplexMatrix lssolve (const ComplexMatrix& b) const; 299 ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info) const; 300 ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info, 301 octave_idx_type& rank) const; 302 ComplexMatrix lssolve (const ComplexMatrix& b, octave_idx_type& info, 303 octave_idx_type& rank, double& rcon) const; 304 305 ColumnVector lssolve (const ColumnVector& b) const; 306 ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info) const; 307 ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info, 308 octave_idx_type& rank) const; 309 ColumnVector lssolve (const ColumnVector& b, octave_idx_type& info, 310 octave_idx_type& rank, double& rcon) const; 311 312 ComplexColumnVector lssolve (const ComplexColumnVector& b) const; 313 ComplexColumnVector lssolve (const ComplexColumnVector& b, 314 octave_idx_type& info) const; 315 ComplexColumnVector lssolve (const ComplexColumnVector& b, 316 octave_idx_type& info, 317 octave_idx_type& rank) const; 318 ComplexColumnVector lssolve (const ComplexColumnVector& b, 319 octave_idx_type& info, 320 octave_idx_type& rank, double& rcon) const; 321 322 Matrix& operator += (const DiagMatrix& a); 323 Matrix& operator -= (const DiagMatrix& a); 324 325 // unary operations 326 327 // other operations 328 329 boolMatrix all (int dim = -1) const; 330 boolMatrix any (int dim = -1) const; 331 332 Matrix cumprod (int dim = -1) const; 333 Matrix cumsum (int dim = -1) const; 334 Matrix prod (int dim = -1) const; 335 Matrix sum (int dim = -1) const; 336 Matrix sumsq (int dim = -1) const; 337 Matrix abs (void) const; 338 339 Matrix diag (octave_idx_type k = 0) const; 340 341 DiagMatrix diag (octave_idx_type m, octave_idx_type n) const; 342 343 ColumnVector row_min (void) const; 344 ColumnVector row_max (void) const; 345 346 ColumnVector row_min (Array<octave_idx_type>& index) const; 347 ColumnVector row_max (Array<octave_idx_type>& index) const; 348 349 RowVector column_min (void) const; 350 RowVector column_max (void) const; 351 352 RowVector column_min (Array<octave_idx_type>& index) const; 353 RowVector column_max (Array<octave_idx_type>& index) const; 354 355 // i/o 356 357 friend OCTAVE_API std::ostream& operator << (std::ostream& os, 358 const Matrix& a); 359 friend OCTAVE_API std::istream& operator >> (std::istream& is, Matrix& a); 360 }; 361 362 // Publish externally used friend functions. 363 364 extern OCTAVE_API Matrix real (const ComplexMatrix& a); 365 extern OCTAVE_API Matrix imag (const ComplexMatrix& a); 366 367 // column vector by row vector -> matrix operations 368 369 extern OCTAVE_API Matrix operator * (const ColumnVector& a, 370 const RowVector& b); 371 372 // Other functions. 373 374 extern OCTAVE_API Matrix Givens (double, double); 375 376 extern OCTAVE_API Matrix Sylvester (const Matrix&, const Matrix&, 377 const Matrix&); 378 379 extern OCTAVE_API Matrix xgemm (const Matrix& a, const Matrix& b, 380 blas_trans_type transa = blas_no_trans, 381 blas_trans_type transb = blas_no_trans); 382 383 extern OCTAVE_API Matrix operator * (const Matrix& a, const Matrix& b); 384 385 extern OCTAVE_API Matrix min (double d, const Matrix& m); 386 extern OCTAVE_API Matrix min (const Matrix& m, double d); 387 extern OCTAVE_API Matrix min (const Matrix& a, const Matrix& b); 388 389 extern OCTAVE_API Matrix max (double d, const Matrix& m); 390 extern OCTAVE_API Matrix max (const Matrix& m, double d); 391 extern OCTAVE_API Matrix max (const Matrix& a, const Matrix& b); 392 393 extern OCTAVE_API Matrix linspace (const ColumnVector& x1, 394 const ColumnVector& x2, 395 octave_idx_type n); 396 397 MS_CMP_OP_DECLS (Matrix, double, OCTAVE_API) 398 MS_BOOL_OP_DECLS (Matrix, double, OCTAVE_API) 399 400 SM_CMP_OP_DECLS (double, Matrix, OCTAVE_API) 401 SM_BOOL_OP_DECLS (double, Matrix, OCTAVE_API) 402 403 MM_CMP_OP_DECLS (Matrix, Matrix, OCTAVE_API) 404 MM_BOOL_OP_DECLS (Matrix, Matrix, OCTAVE_API) 405 406 MARRAY_FORWARD_DEFS (MArray, Matrix, double) 407 408 template <typename T> 409 void read_int (std::istream& is, bool swap_bytes, T& val); 410 411 #endif 412