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 (HAVE_CONFIG_H) 27 # include "config.h" 28 #endif 29 30 #include <cstddef> 31 32 #include "CSparse.h" 33 #include "MatrixType.h" 34 #include "dRowVector.h" 35 #include "dSparse.h" 36 #include "lo-error.h" 37 #include "oct-cmplx.h" 38 #include "oct-refcount.h" 39 #include "oct-sparse.h" 40 #include "oct-spparms.h" 41 #include "quit.h" 42 #include "sparse-chol.h" 43 #include "sparse-util.h" 44 45 namespace octave 46 { 47 namespace math 48 { 49 template <typename chol_type> 50 class sparse_chol<chol_type>::sparse_chol_rep 51 { 52 public: 53 sparse_chol_rep(void)54 sparse_chol_rep (void) 55 : count (1), is_pd (false), minor_p (0), perms (), cond (0) 56 #if defined (HAVE_CHOLMOD) 57 , Lsparse (nullptr), Common () 58 #endif 59 { } 60 sparse_chol_rep(const chol_type & a,bool natural,bool force)61 sparse_chol_rep (const chol_type& a, bool natural, bool force) 62 : count (1), is_pd (false), minor_p (0), perms (), cond (0) 63 #if defined (HAVE_CHOLMOD) 64 , Lsparse (nullptr), Common () 65 #endif 66 { 67 init (a, natural, force); 68 } 69 sparse_chol_rep(const chol_type & a,octave_idx_type & info,bool natural,bool force)70 sparse_chol_rep (const chol_type& a, octave_idx_type& info, 71 bool natural, bool force) 72 : count (1), is_pd (false), minor_p (0), perms (), cond (0) 73 #if defined (HAVE_CHOLMOD) 74 , Lsparse (nullptr), Common () 75 #endif 76 { 77 info = init (a, natural, force); 78 } 79 80 // No copying! 81 82 sparse_chol_rep (const sparse_chol_rep&) = delete; 83 84 sparse_chol_rep& operator = (const sparse_chol_rep&) = delete; 85 ~sparse_chol_rep(void)86 ~sparse_chol_rep (void) 87 { 88 #if defined (HAVE_CHOLMOD) 89 if (Lsparse) 90 CHOLMOD_NAME (free_sparse) (&Lsparse, &Common); 91 92 CHOLMOD_NAME(finish) (&Common); 93 #endif 94 } 95 96 #if defined (HAVE_CHOLMOD) L(void) const97 cholmod_sparse * L (void) const 98 { 99 return Lsparse; 100 } 101 #endif 102 P(void) const103 octave_idx_type P (void) const 104 { 105 #if defined (HAVE_CHOLMOD) 106 return (minor_p == static_cast<octave_idx_type> (Lsparse->ncol) ? 107 0 : minor_p + 1); 108 #else 109 return 0; 110 #endif 111 } 112 perm(void) const113 RowVector perm (void) const { return perms + 1; } 114 115 SparseMatrix Q (void) const; 116 is_positive_definite(void) const117 bool is_positive_definite (void) const { return is_pd; } 118 rcond(void) const119 double rcond (void) const { return cond; } 120 121 refcount<octave_idx_type> count; 122 123 private: 124 125 bool is_pd; 126 127 octave_idx_type minor_p; 128 129 RowVector perms; 130 131 double cond; 132 133 #if defined (HAVE_CHOLMOD) 134 cholmod_sparse *Lsparse; 135 136 cholmod_common Common; 137 138 void drop_zeros (const cholmod_sparse *S); 139 #endif 140 141 octave_idx_type init (const chol_type& a, bool natural, bool force); 142 }; 143 144 #if defined (HAVE_CHOLMOD) 145 146 // Can't use CHOLMOD_NAME(drop)(0.0, S, cm) because it doesn't treat 147 // complex matrices. 148 149 template <typename chol_type> 150 void drop_zeros(const cholmod_sparse * S)151 sparse_chol<chol_type>::sparse_chol_rep::drop_zeros (const cholmod_sparse *S) 152 { 153 if (! S) 154 return; 155 156 octave_idx_type *Sp = static_cast<octave_idx_type *>(S->p); 157 octave_idx_type *Si = static_cast<octave_idx_type *>(S->i); 158 chol_elt *Sx = static_cast<chol_elt *>(S->x); 159 160 octave_idx_type pdest = 0; 161 octave_idx_type ncol = S->ncol; 162 163 for (octave_idx_type k = 0; k < ncol; k++) 164 { 165 octave_idx_type p = Sp[k]; 166 octave_idx_type pend = Sp[k+1]; 167 Sp[k] = pdest; 168 169 for (; p < pend; p++) 170 { 171 chol_elt sik = Sx[p]; 172 173 if (CHOLMOD_IS_NONZERO (sik)) 174 { 175 if (p != pdest) 176 { 177 Si[pdest] = Si[p]; 178 Sx[pdest] = sik; 179 } 180 181 pdest++; 182 } 183 } 184 } 185 186 Sp[ncol] = pdest; 187 } 188 189 // Must provide a specialization for this function. 190 template <typename T> 191 int 192 get_xtype (void); 193 194 template <> 195 inline int get_xtype(void)196 get_xtype<double> (void) 197 { 198 return CHOLMOD_REAL; 199 } 200 201 template <> 202 inline int get_xtype(void)203 get_xtype<Complex> (void) 204 { 205 return CHOLMOD_COMPLEX; 206 } 207 208 #endif 209 210 template <typename chol_type> 211 octave_idx_type init(const chol_type & a,bool natural,bool force)212 sparse_chol<chol_type>::sparse_chol_rep::init (const chol_type& a, 213 bool natural, bool force) 214 { 215 volatile octave_idx_type info = 0; 216 217 #if defined (HAVE_CHOLMOD) 218 219 octave_idx_type a_nr = a.rows (); 220 octave_idx_type a_nc = a.cols (); 221 222 if (a_nr != a_nc) 223 (*current_liboctave_error_handler) 224 ("sparse_chol requires square matrix"); 225 226 cholmod_common *cm = &Common; 227 228 // Setup initial parameters 229 230 CHOLMOD_NAME(start) (cm); 231 cm->prefer_zomplex = false; 232 233 double spu = octave_sparse_params::get_key ("spumoni"); 234 235 if (spu == 0.) 236 { 237 cm->print = -1; 238 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, nullptr); 239 } 240 else 241 { 242 cm->print = static_cast<int> (spu) + 2; 243 SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 244 &SparseCholPrint); 245 } 246 247 cm->error_handler = &SparseCholError; 248 249 SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, 250 divcomplex); 251 252 SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); 253 254 cm->final_asis = false; 255 cm->final_super = false; 256 cm->final_ll = true; 257 cm->final_pack = true; 258 cm->final_monotonic = true; 259 cm->final_resymbol = false; 260 261 cholmod_sparse A; 262 cholmod_sparse *ac = &A; 263 double dummy; 264 265 ac->nrow = a_nr; 266 ac->ncol = a_nc; 267 268 ac->p = a.cidx (); 269 ac->i = a.ridx (); 270 ac->nzmax = a.nnz (); 271 ac->packed = true; 272 ac->sorted = true; 273 ac->nz = nullptr; 274 #if defined (OCTAVE_ENABLE_64) 275 ac->itype = CHOLMOD_LONG; 276 #else 277 ac->itype = CHOLMOD_INT; 278 #endif 279 ac->dtype = CHOLMOD_DOUBLE; 280 ac->stype = 1; 281 ac->xtype = get_xtype<chol_elt> (); 282 283 if (a_nr < 1) 284 ac->x = &dummy; 285 else 286 ac->x = a.data (); 287 288 // use natural ordering if no q output parameter 289 if (natural) 290 { 291 cm->nmethods = 1; 292 cm->method[0].ordering = CHOLMOD_NATURAL; 293 cm->postorder = false; 294 } 295 296 cholmod_factor *Lfactor; 297 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 298 Lfactor = CHOLMOD_NAME(analyze) (ac, cm); 299 CHOLMOD_NAME(factorize) (ac, Lfactor, cm); 300 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 301 302 is_pd = cm->status == CHOLMOD_OK; 303 info = (is_pd ? 0 : cm->status); 304 305 if (is_pd || force) 306 { 307 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 308 cond = CHOLMOD_NAME(rcond) (Lfactor, cm); 309 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 310 311 minor_p = Lfactor->minor; 312 313 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 314 Lsparse = CHOLMOD_NAME(factor_to_sparse) (Lfactor, cm); 315 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 316 317 if (minor_p > 0 && minor_p < a_nr) 318 { 319 std::size_t n1 = a_nr + 1; 320 Lsparse->p = CHOLMOD_NAME(realloc) (minor_p+1, 321 sizeof(octave_idx_type), 322 Lsparse->p, &n1, cm); 323 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 324 CHOLMOD_NAME(reallocate_sparse) 325 (static_cast<octave_idx_type *>(Lsparse->p)[minor_p], 326 Lsparse, cm); 327 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 328 329 Lsparse->ncol = minor_p; 330 } 331 332 drop_zeros (Lsparse); 333 334 if (! natural) 335 { 336 perms.resize (a_nr); 337 for (octave_idx_type i = 0; i < a_nr; i++) 338 perms(i) = static_cast<octave_idx_type *>(Lfactor->Perm)[i]; 339 } 340 } 341 342 // NAME used to prefix statistics report from print_common 343 static char blank_name[] = " "; 344 345 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 346 CHOLMOD_NAME(print_common) (blank_name, cm); 347 CHOLMOD_NAME(free_factor) (&Lfactor, cm); 348 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 349 350 return info; 351 352 #else 353 354 octave_unused_parameter (a); 355 octave_unused_parameter (natural); 356 octave_unused_parameter (force); 357 358 (*current_liboctave_error_handler) 359 ("support for CHOLMOD was unavailable or disabled when liboctave was built"); 360 361 return info; 362 363 #endif 364 } 365 366 template <typename chol_type> 367 SparseMatrix Q(void) const368 sparse_chol<chol_type>::sparse_chol_rep::Q (void) const 369 { 370 #if defined (HAVE_CHOLMOD) 371 372 octave_idx_type n = Lsparse->nrow; 373 SparseMatrix p (n, n, n); 374 375 for (octave_idx_type i = 0; i < n; i++) 376 { 377 p.xcidx (i) = i; 378 p.xridx (i) = static_cast<octave_idx_type> (perms (i)); 379 p.xdata (i) = 1; 380 } 381 382 p.xcidx (n) = n; 383 384 return p; 385 386 #else 387 388 return SparseMatrix (); 389 390 #endif 391 } 392 393 template <typename chol_type> sparse_chol(void)394 sparse_chol<chol_type>::sparse_chol (void) 395 : rep (new typename sparse_chol<chol_type>::sparse_chol_rep ()) 396 { } 397 398 template <typename chol_type> sparse_chol(const chol_type & a,bool natural,bool force)399 sparse_chol<chol_type>::sparse_chol (const chol_type& a, bool natural, 400 bool force) 401 : rep (new typename 402 sparse_chol<chol_type>::sparse_chol_rep (a, natural, force)) 403 { } 404 405 template <typename chol_type> sparse_chol(const chol_type & a,octave_idx_type & info,bool natural,bool force)406 sparse_chol<chol_type>::sparse_chol (const chol_type& a, 407 octave_idx_type& info, 408 bool natural, bool force) 409 : rep (new typename 410 sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, force)) 411 { } 412 413 template <typename chol_type> sparse_chol(const chol_type & a,octave_idx_type & info,bool natural)414 sparse_chol<chol_type>::sparse_chol (const chol_type& a, 415 octave_idx_type& info, 416 bool natural) 417 : rep (new typename 418 sparse_chol<chol_type>::sparse_chol_rep (a, info, natural, false)) 419 { } 420 421 template <typename chol_type> sparse_chol(const chol_type & a,octave_idx_type & info)422 sparse_chol<chol_type>::sparse_chol (const chol_type& a, 423 octave_idx_type& info) 424 : rep (new typename 425 sparse_chol<chol_type>::sparse_chol_rep (a, info, false, false)) 426 { } 427 428 template <typename chol_type> sparse_chol(const sparse_chol<chol_type> & a)429 sparse_chol<chol_type>::sparse_chol (const sparse_chol<chol_type>& a) 430 : rep (a.rep) 431 { 432 rep->count++; 433 } 434 435 template <typename chol_type> ~sparse_chol(void)436 sparse_chol<chol_type>::~sparse_chol (void) 437 { 438 if (--rep->count == 0) 439 delete rep; 440 } 441 442 template <typename chol_type> 443 sparse_chol<chol_type>& operator =(const sparse_chol & a)444 sparse_chol<chol_type>::operator = (const sparse_chol& a) 445 { 446 if (this != &a) 447 { 448 if (--rep->count == 0) 449 delete rep; 450 451 rep = a.rep; 452 rep->count++; 453 } 454 455 return *this; 456 } 457 458 template <typename chol_type> 459 chol_type L(void) const460 sparse_chol<chol_type>::L (void) const 461 { 462 #if defined (HAVE_CHOLMOD) 463 464 cholmod_sparse *m = rep->L (); 465 466 octave_idx_type nc = m->ncol; 467 octave_idx_type nnz = m->nzmax; 468 469 chol_type ret (m->nrow, nc, nnz); 470 471 for (octave_idx_type j = 0; j < nc+1; j++) 472 ret.xcidx (j) = static_cast<octave_idx_type *>(m->p)[j]; 473 474 for (octave_idx_type i = 0; i < nnz; i++) 475 { 476 ret.xridx (i) = static_cast<octave_idx_type *>(m->i)[i]; 477 ret.xdata (i) = static_cast<chol_elt *>(m->x)[i]; 478 } 479 480 return ret; 481 482 #else 483 484 return chol_type (); 485 486 #endif 487 } 488 489 template <typename chol_type> 490 octave_idx_type P(void) const491 sparse_chol<chol_type>::P (void) const 492 { 493 return rep->P (); 494 } 495 496 template <typename chol_type> 497 RowVector perm(void) const498 sparse_chol<chol_type>::perm (void) const 499 { 500 return rep->perm (); 501 } 502 503 template <typename chol_type> 504 SparseMatrix Q(void) const505 sparse_chol<chol_type>::Q (void) const 506 { 507 return rep->Q (); 508 } 509 510 template <typename chol_type> 511 bool is_positive_definite(void) const512 sparse_chol<chol_type>::is_positive_definite (void) const 513 { 514 return rep->is_positive_definite (); 515 } 516 517 template <typename chol_type> 518 double rcond(void) const519 sparse_chol<chol_type>::rcond (void) const 520 { 521 return rep->rcond (); 522 } 523 524 template <typename chol_type> 525 chol_type inverse(void) const526 sparse_chol<chol_type>::inverse (void) const 527 { 528 chol_type retval; 529 530 #if defined (HAVE_CHOLMOD) 531 532 cholmod_sparse *m = rep->L (); 533 octave_idx_type n = m->ncol; 534 RowVector perms = rep->perm (); 535 double rcond2; 536 octave_idx_type info; 537 MatrixType mattype (MatrixType::Upper); 538 chol_type linv = L ().hermitian ().inverse (mattype, info, rcond2, 1, 0); 539 540 if (perms.numel () == n) 541 { 542 SparseMatrix Qc = Q (); 543 544 retval = Qc * linv * linv.hermitian () * Qc.transpose (); 545 } 546 else 547 retval = linv * linv.hermitian (); 548 549 #endif 550 551 return retval; 552 } 553 554 template <typename chol_type> 555 chol_type chol2inv(const chol_type & r)556 chol2inv (const chol_type& r) 557 { 558 octave_idx_type r_nr = r.rows (); 559 octave_idx_type r_nc = r.cols (); 560 chol_type retval; 561 562 if (r_nr != r_nc) 563 (*current_liboctave_error_handler) ("U must be a square matrix"); 564 565 MatrixType mattype (r); 566 int typ = mattype.type (false); 567 double rcond; 568 octave_idx_type info; 569 chol_type rtra, multip; 570 571 if (typ == MatrixType::Upper) 572 { 573 rtra = r.transpose (); 574 multip = (rtra*r); 575 } 576 else if (typ == MatrixType::Lower) 577 { 578 rtra = r.transpose (); 579 multip = (r*rtra); 580 } 581 else 582 (*current_liboctave_error_handler) ("U must be a triangular matrix"); 583 584 MatrixType mattypenew (multip); 585 retval = multip.inverse (mattypenew, info, rcond, true, false); 586 return retval; 587 } 588 589 // SparseComplexMatrix specialization (the value for the NATURAL 590 // parameter in the sparse_chol<T>::sparse_chol_rep constructor is 591 // different from the default). 592 593 template <> sparse_chol(const SparseComplexMatrix & a,octave_idx_type & info)594 sparse_chol<SparseComplexMatrix>::sparse_chol (const SparseComplexMatrix& a, 595 octave_idx_type& info) 596 : rep (new sparse_chol<SparseComplexMatrix>::sparse_chol_rep (a, info, 597 true, 598 false)) 599 { } 600 601 // Instantiations we need. 602 603 template class sparse_chol<SparseMatrix>; 604 605 template class sparse_chol<SparseComplexMatrix>; 606 607 template SparseMatrix 608 chol2inv<SparseMatrix> (const SparseMatrix& r); 609 610 template SparseComplexMatrix 611 chol2inv<SparseComplexMatrix> (const SparseComplexMatrix& r); 612 } 613 } 614