1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 1998-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_dSparse_h) 27 #define octave_dSparse_h 1 28 29 #include "octave-config.h" 30 31 #include "CColVector.h" 32 #include "CMatrix.h" 33 #include "DET.h" 34 #include "MSparse.h" 35 #include "MatrixType.h" 36 #include "Sparse-op-decls.h" 37 #include "dColVector.h" 38 #include "dMatrix.h" 39 #include "dNDArray.h" 40 41 class PermMatrix; 42 class DiagMatrix; 43 class SparseComplexMatrix; 44 class SparseBoolMatrix; 45 46 class 47 OCTAVE_API 48 SparseMatrix : public MSparse<double> 49 { 50 public: 51 52 // Corresponding dense matrix type for this sparse matrix type. 53 typedef Matrix dense_matrix_type; 54 55 typedef void (*solve_singularity_handler) (double rcond); 56 SparseMatrix(void)57 SparseMatrix (void) : MSparse<double> () { } 58 SparseMatrix(octave_idx_type r,octave_idx_type c)59 SparseMatrix (octave_idx_type r, octave_idx_type c) 60 : MSparse<double> (r, c) { } 61 62 SparseMatrix (const dim_vector& dv, octave_idx_type nz = 0) : 63 MSparse<double> (dv, nz) { } 64 SparseMatrix(octave_idx_type r,octave_idx_type c,double val)65 explicit SparseMatrix (octave_idx_type r, octave_idx_type c, double val) 66 : MSparse<double> (r, c, val) { } 67 SparseMatrix(const SparseMatrix & a)68 SparseMatrix (const SparseMatrix& a) : MSparse<double> (a) { } 69 SparseMatrix(const SparseMatrix & a,const dim_vector & dv)70 SparseMatrix (const SparseMatrix& a, const dim_vector& dv) 71 : MSparse<double> (a, dv) { } 72 SparseMatrix(const MSparse<double> & a)73 SparseMatrix (const MSparse<double>& a) : MSparse<double> (a) { } 74 SparseMatrix(const Sparse<double> & a)75 SparseMatrix (const Sparse<double>& a) : MSparse<double> (a) { } 76 77 explicit SparseMatrix (const SparseBoolMatrix& a); 78 SparseMatrix(const Matrix & a)79 explicit SparseMatrix (const Matrix& a) : MSparse<double> (a) { } 80 SparseMatrix(const NDArray & a)81 explicit SparseMatrix (const NDArray& a) : MSparse<double> (a) { } 82 83 SparseMatrix (const Array<double>& a, const idx_vector& r, 84 const idx_vector& c, octave_idx_type nr = -1, 85 octave_idx_type nc = -1, bool sum_terms = true, 86 octave_idx_type nzm = -1) 87 : MSparse<double> (a, r, c, nr, nc, sum_terms, nzm) { } 88 89 explicit SparseMatrix (const DiagMatrix& a); 90 SparseMatrix(const PermMatrix & a)91 explicit SparseMatrix (const PermMatrix& a) : MSparse<double>(a) { } 92 SparseMatrix(octave_idx_type r,octave_idx_type c,octave_idx_type num_nz)93 SparseMatrix (octave_idx_type r, octave_idx_type c, 94 octave_idx_type num_nz) : MSparse<double> (r, c, num_nz) { } 95 96 SparseMatrix& operator = (const SparseMatrix& a) 97 { 98 MSparse<double>::operator = (a); 99 return *this; 100 } 101 102 bool operator == (const SparseMatrix& a) const; 103 bool operator != (const SparseMatrix& a) const; 104 105 bool issymmetric (void) const; 106 107 SparseMatrix max (int dim = -1) const; 108 SparseMatrix max (Array<octave_idx_type>& index, int dim = -1) const; 109 SparseMatrix min (int dim = -1) const; 110 SparseMatrix min (Array<octave_idx_type>& index, int dim = -1) const; 111 112 // destructive insert/delete/reorder operations 113 114 SparseMatrix& insert (const SparseMatrix& a, octave_idx_type r, 115 octave_idx_type c); 116 117 SparseMatrix& insert (const SparseMatrix& a, 118 const Array<octave_idx_type>& indx); 119 120 SparseMatrix concat (const SparseMatrix& rb, 121 const Array<octave_idx_type>& ra_idx); 122 SparseComplexMatrix concat (const SparseComplexMatrix& rb, 123 const Array<octave_idx_type>& ra_idx); 124 125 friend OCTAVE_API SparseMatrix real (const SparseComplexMatrix& a); 126 friend OCTAVE_API SparseMatrix imag (const SparseComplexMatrix& a); 127 transpose(void)128 SparseMatrix transpose (void) const 129 { 130 return MSparse<double>::transpose (); 131 } hermitian(void)132 SparseMatrix hermitian (void) const { return transpose (); } 133 134 // extract row or column i. 135 136 RowVector row (octave_idx_type i) const; 137 138 ColumnVector column (octave_idx_type i) const; 139 140 private: 141 SparseMatrix dinverse (MatrixType& mattype, octave_idx_type& info, 142 double& rcond, const bool force = false, 143 const bool calccond = true) const; 144 145 SparseMatrix tinverse (MatrixType& mattype, octave_idx_type& info, 146 double& rcond, const bool force = false, 147 const bool calccond = true) const; 148 149 public: 150 SparseMatrix inverse (void) const; 151 SparseMatrix inverse (MatrixType& mattype) const; 152 SparseMatrix inverse (MatrixType& mattype, octave_idx_type& info) const; 153 SparseMatrix inverse (MatrixType& mattype, octave_idx_type& info, 154 double& rcond, bool force = false, 155 bool calc_cond = true) const; 156 157 DET determinant (void) const; 158 DET determinant (octave_idx_type& info) const; 159 DET determinant (octave_idx_type& info, double& rcond, 160 bool calc_cond = true) const; 161 162 private: 163 // Diagonal matrix solvers 164 Matrix dsolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 165 double& rcond, solve_singularity_handler sing_handler, 166 bool calc_cond = false) const; 167 168 ComplexMatrix dsolve (MatrixType& typ, const ComplexMatrix& b, 169 octave_idx_type& info, double& rcond, 170 solve_singularity_handler sing_handler, 171 bool calc_cond = false) const; 172 173 SparseMatrix dsolve (MatrixType& typ, const SparseMatrix& b, 174 octave_idx_type& info, double& rcond, 175 solve_singularity_handler sing_handler, 176 bool calc_cond = false) const; 177 178 SparseComplexMatrix dsolve (MatrixType& typ, const SparseComplexMatrix& b, 179 octave_idx_type& info, double& rcond, 180 solve_singularity_handler sing_handler, 181 bool calc_cond = false) const; 182 183 // Upper triangular matrix solvers 184 Matrix utsolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 185 double& rcond, solve_singularity_handler sing_handler, 186 bool calc_cond = false) const; 187 188 ComplexMatrix utsolve (MatrixType& typ, const ComplexMatrix& b, 189 octave_idx_type& info, double& rcond, 190 solve_singularity_handler sing_handler, 191 bool calc_cond = false) const; 192 193 SparseMatrix utsolve (MatrixType& typ, const SparseMatrix& b, 194 octave_idx_type& info, double& rcond, 195 solve_singularity_handler sing_handler, 196 bool calc_cond = false) const; 197 198 SparseComplexMatrix utsolve (MatrixType& typ, const SparseComplexMatrix& b, 199 octave_idx_type& info, double& rcond, 200 solve_singularity_handler sing_handler, 201 bool calc_cond = false) const; 202 203 // Lower triangular matrix solvers 204 Matrix ltsolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 205 double& rcond, solve_singularity_handler sing_handler, 206 bool calc_cond = false) const; 207 208 ComplexMatrix ltsolve (MatrixType& typ, const ComplexMatrix& b, 209 octave_idx_type& info, double& rcond, 210 solve_singularity_handler sing_handler, 211 bool calc_cond = false) const; 212 213 SparseMatrix ltsolve (MatrixType& typ, const SparseMatrix& b, 214 octave_idx_type& info, double& rcond, 215 solve_singularity_handler sing_handler, 216 bool calc_cond = false) const; 217 218 SparseComplexMatrix ltsolve (MatrixType& typ, const SparseComplexMatrix& b, 219 octave_idx_type& info, double& rcond, 220 solve_singularity_handler sing_handler, 221 bool calc_cond = false) const; 222 223 // Tridiagonal matrix solvers 224 Matrix trisolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 225 double& rcond, solve_singularity_handler sing_handler, 226 bool calc_cond = false) const; 227 228 ComplexMatrix trisolve (MatrixType& typ, const ComplexMatrix& b, 229 octave_idx_type& info, double& rcond, 230 solve_singularity_handler sing_handler, 231 bool calc_cond = false) const; 232 233 SparseMatrix trisolve (MatrixType& typ, const SparseMatrix& b, 234 octave_idx_type& info, double& rcond, 235 solve_singularity_handler sing_handler, 236 bool calc_cond = false) const; 237 238 SparseComplexMatrix trisolve (MatrixType& typ, const SparseComplexMatrix& b, 239 octave_idx_type& info, double& rcond, 240 solve_singularity_handler sing_handler, 241 bool calc_cond = false) const; 242 243 // Banded matrix solvers (umfpack/cholesky) 244 Matrix bsolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 245 double& rcond, solve_singularity_handler sing_handler, 246 bool calc_cond = false) const; 247 248 ComplexMatrix bsolve (MatrixType& typ, const ComplexMatrix& b, 249 octave_idx_type& info, double& rcond, 250 solve_singularity_handler sing_handler, 251 bool calc_cond = false) const; 252 253 SparseMatrix bsolve (MatrixType& typ, const SparseMatrix& b, 254 octave_idx_type& info, double& rcond, 255 solve_singularity_handler sing_handler, 256 bool calc_cond = false) const; 257 258 SparseComplexMatrix bsolve (MatrixType& typ, const SparseComplexMatrix& b, 259 octave_idx_type& info, double& rcond, 260 solve_singularity_handler sing_handler, 261 bool calc_cond = false) const; 262 263 // Full matrix solvers (umfpack/cholesky) 264 void * factorize (octave_idx_type& err, double& rcond, Matrix& Control, 265 Matrix& Info, solve_singularity_handler sing_handler, 266 bool calc_cond = false) const; 267 268 Matrix fsolve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 269 double& rcond, solve_singularity_handler sing_handler, 270 bool calc_cond = false) const; 271 272 ComplexMatrix fsolve (MatrixType& typ, const ComplexMatrix& b, 273 octave_idx_type& info, double& rcond, 274 solve_singularity_handler sing_handler, 275 bool calc_cond = false) const; 276 277 SparseMatrix fsolve (MatrixType& typ, const SparseMatrix& b, 278 octave_idx_type& info, double& rcond, 279 solve_singularity_handler sing_handler, 280 bool calc_cond = false) const; 281 282 SparseComplexMatrix fsolve (MatrixType& typ, const SparseComplexMatrix& b, 283 octave_idx_type& info, double& rcond, 284 solve_singularity_handler sing_handler, 285 bool calc_cond = false) const; 286 287 public: 288 // Generic interface to solver with no probing of type 289 Matrix solve (MatrixType& typ, const Matrix& b) const; 290 Matrix solve (MatrixType& typ, const Matrix& b, octave_idx_type& info) const; 291 Matrix solve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 292 double& rcond) const; 293 Matrix solve (MatrixType& typ, const Matrix& b, octave_idx_type& info, 294 double& rcond, solve_singularity_handler sing_handler, 295 bool singular_fallback = true) const; 296 297 ComplexMatrix solve (MatrixType& typ, const ComplexMatrix& b) const; 298 ComplexMatrix solve (MatrixType& typ, const ComplexMatrix& b, 299 octave_idx_type& info) const; 300 ComplexMatrix solve (MatrixType& typ, const ComplexMatrix& b, 301 octave_idx_type& info, double& rcond) const; 302 ComplexMatrix solve (MatrixType& typ, const ComplexMatrix& b, 303 octave_idx_type& info, double& rcond, 304 solve_singularity_handler sing_handler, 305 bool singular_fallback = true) const; 306 307 SparseMatrix solve (MatrixType& typ, const SparseMatrix& b) const; 308 SparseMatrix solve (MatrixType& typ, const SparseMatrix& b, 309 octave_idx_type& info) const; 310 SparseMatrix solve (MatrixType& typ, const SparseMatrix& b, 311 octave_idx_type& info, double& rcond) const; 312 SparseMatrix solve (MatrixType& typ, const SparseMatrix& b, 313 octave_idx_type& info, double& rcond, 314 solve_singularity_handler sing_handler, 315 bool singular_fallback = true) const; 316 317 SparseComplexMatrix solve (MatrixType& typ, 318 const SparseComplexMatrix& b) const; 319 SparseComplexMatrix solve (MatrixType& typ, const SparseComplexMatrix& b, 320 octave_idx_type& info) const; 321 SparseComplexMatrix solve (MatrixType& typ, const SparseComplexMatrix& b, 322 octave_idx_type& info, double& rcond) const; 323 SparseComplexMatrix solve (MatrixType& typ, const SparseComplexMatrix& b, 324 octave_idx_type& info, double& rcond, 325 solve_singularity_handler sing_handler, 326 bool singular_fallabck = true) const; 327 328 ColumnVector solve (MatrixType& typ, const ColumnVector& b) const; 329 ColumnVector solve (MatrixType& typ, const ColumnVector& b, 330 octave_idx_type& info) const; 331 ColumnVector solve (MatrixType& typ, const ColumnVector& b, 332 octave_idx_type& info, double& rcond) const; 333 ColumnVector solve (MatrixType& typ, const ColumnVector& b, 334 octave_idx_type& info, double& rcond, 335 solve_singularity_handler sing_handler) const; 336 337 ComplexColumnVector solve (MatrixType& typ, 338 const ComplexColumnVector& b) const; 339 ComplexColumnVector solve (MatrixType& typ, const ComplexColumnVector& b, 340 octave_idx_type& info) const; 341 ComplexColumnVector solve (MatrixType& typ, const ComplexColumnVector& b, 342 octave_idx_type& info, double& rcond) const; 343 ComplexColumnVector solve (MatrixType& typ, const ComplexColumnVector& b, 344 octave_idx_type& info, double& rcond, 345 solve_singularity_handler sing_handler) const; 346 347 // Generic interface to solver with probing of type 348 Matrix solve (const Matrix& b) const; 349 Matrix solve (const Matrix& b, octave_idx_type& info) const; 350 Matrix solve (const Matrix& b, octave_idx_type& info, double& rcond) const; 351 Matrix solve (const Matrix& b, octave_idx_type& info, double& rcond, 352 solve_singularity_handler sing_handler) const; 353 354 ComplexMatrix solve (const ComplexMatrix& b) const; 355 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info) const; 356 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, 357 double& rcond) const; 358 ComplexMatrix solve (const ComplexMatrix& b, octave_idx_type& info, 359 double& rcond, 360 solve_singularity_handler sing_handler) const; 361 362 SparseMatrix solve (const SparseMatrix& b) const; 363 SparseMatrix solve (const SparseMatrix& b, octave_idx_type& info) const; 364 SparseMatrix solve (const SparseMatrix& b, octave_idx_type& info, 365 double& rcond) const; 366 SparseMatrix solve (const SparseMatrix& b, octave_idx_type& info, 367 double& rcond, 368 solve_singularity_handler sing_handler) const; 369 370 SparseComplexMatrix solve (const SparseComplexMatrix& b) const; 371 SparseComplexMatrix solve (const SparseComplexMatrix& b, 372 octave_idx_type& info) const; 373 SparseComplexMatrix solve (const SparseComplexMatrix& b, 374 octave_idx_type& info, double& rcond) const; 375 SparseComplexMatrix solve (const SparseComplexMatrix& b, 376 octave_idx_type& info, double& rcond, 377 solve_singularity_handler sing_handler) const; 378 379 ColumnVector solve (const ColumnVector& b) const; 380 ColumnVector solve (const ColumnVector& b, octave_idx_type& info) const; 381 ColumnVector solve (const ColumnVector& b, octave_idx_type& info, 382 double& rcond) const; 383 ColumnVector solve (const ColumnVector& b, octave_idx_type& info, 384 double& rcond, 385 solve_singularity_handler sing_handler) const; 386 387 ComplexColumnVector solve (const ComplexColumnVector& b) const; 388 ComplexColumnVector solve (const ComplexColumnVector& b, 389 octave_idx_type& info) const; 390 ComplexColumnVector solve (const ComplexColumnVector& b, 391 octave_idx_type& info, double& rcond) const; 392 ComplexColumnVector solve (const ComplexColumnVector& b, 393 octave_idx_type& info, double& rcond, 394 solve_singularity_handler sing_handler) const; 395 396 // other operations 397 398 bool any_element_is_negative (bool = false) const; 399 bool any_element_is_nan (void) const; 400 bool any_element_is_inf_or_nan (void) const; 401 bool any_element_not_one_or_zero (void) const; 402 bool all_elements_are_zero (void) const; 403 bool all_elements_are_int_or_inf_or_nan (void) const; 404 bool all_integers (double& max_val, double& min_val) const; 405 bool too_large_for_float (void) const; 406 407 SparseBoolMatrix operator ! (void) const; 408 409 SparseBoolMatrix all (int dim = -1) const; 410 SparseBoolMatrix any (int dim = -1) const; 411 412 SparseMatrix cumprod (int dim = -1) const; 413 SparseMatrix cumsum (int dim = -1) const; 414 SparseMatrix prod (int dim = -1) const; 415 SparseMatrix sum (int dim = -1) const; 416 SparseMatrix sumsq (int dim = -1) const; 417 SparseMatrix abs (void) const; 418 419 SparseMatrix diag (octave_idx_type k = 0) const; 420 421 Matrix matrix_value (void) const; 422 423 SparseMatrix squeeze (void) const; 424 425 SparseMatrix reshape (const dim_vector& new_dims) const; 426 427 SparseMatrix permute (const Array<octave_idx_type>& vec, 428 bool inv = false) const; 429 430 SparseMatrix ipermute (const Array<octave_idx_type>& vec) const; 431 432 // i/o 433 434 friend OCTAVE_API std::ostream& operator << (std::ostream& os, 435 const SparseMatrix& a); 436 friend OCTAVE_API std::istream& operator >> (std::istream& is, 437 SparseMatrix& a); 438 439 }; 440 441 // Publish externally used friend functions. 442 443 extern OCTAVE_API SparseMatrix real (const SparseComplexMatrix& a); 444 extern OCTAVE_API SparseMatrix imag (const SparseComplexMatrix& a); 445 446 // Other operators. 447 448 extern OCTAVE_API SparseMatrix operator * (const SparseMatrix& a, 449 const SparseMatrix& b); 450 extern OCTAVE_API Matrix operator * (const Matrix& a, 451 const SparseMatrix& b); 452 extern OCTAVE_API Matrix mul_trans (const Matrix& a, 453 const SparseMatrix& b); 454 extern OCTAVE_API Matrix operator * (const SparseMatrix& a, 455 const Matrix& b); 456 extern OCTAVE_API Matrix trans_mul (const SparseMatrix& a, 457 const Matrix& b); 458 459 extern OCTAVE_API SparseMatrix operator * (const DiagMatrix&, 460 const SparseMatrix&); 461 extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, 462 const DiagMatrix&); 463 464 extern OCTAVE_API SparseMatrix operator + (const DiagMatrix&, 465 const SparseMatrix&); 466 extern OCTAVE_API SparseMatrix operator + (const SparseMatrix&, 467 const DiagMatrix&); 468 extern OCTAVE_API SparseMatrix operator - (const DiagMatrix&, 469 const SparseMatrix&); 470 extern OCTAVE_API SparseMatrix operator - (const SparseMatrix&, 471 const DiagMatrix&); 472 473 extern OCTAVE_API SparseMatrix operator * (const PermMatrix&, 474 const SparseMatrix&); 475 extern OCTAVE_API SparseMatrix operator * (const SparseMatrix&, 476 const PermMatrix&); 477 478 extern OCTAVE_API SparseMatrix min (double d, const SparseMatrix& m); 479 extern OCTAVE_API SparseMatrix min (const SparseMatrix& m, double d); 480 extern OCTAVE_API SparseMatrix min (const SparseMatrix& a, 481 const SparseMatrix& b); 482 483 extern OCTAVE_API SparseMatrix max (double d, const SparseMatrix& m); 484 extern OCTAVE_API SparseMatrix max (const SparseMatrix& m, double d); 485 extern OCTAVE_API SparseMatrix max (const SparseMatrix& a, 486 const SparseMatrix& b); 487 488 SPARSE_SMS_CMP_OP_DECLS (SparseMatrix, double, OCTAVE_API) 489 SPARSE_SMS_BOOL_OP_DECLS (SparseMatrix, double, OCTAVE_API) 490 491 SPARSE_SSM_CMP_OP_DECLS (double, SparseMatrix, OCTAVE_API) 492 SPARSE_SSM_BOOL_OP_DECLS (double, SparseMatrix, OCTAVE_API) 493 494 SPARSE_SMSM_CMP_OP_DECLS (SparseMatrix, SparseMatrix, OCTAVE_API) 495 SPARSE_SMSM_BOOL_OP_DECLS (SparseMatrix, SparseMatrix, OCTAVE_API) 496 497 SPARSE_FORWARD_DEFS (MSparse, SparseMatrix, Matrix, double) 498 499 #endif 500