1 //////////////////////////////////////////////////////////////////////// 2 // 3 // Copyright (C) 2005-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 "CMatrix.h" 31 #include "CSparse.h" 32 #include "MArray.h" 33 #include "dColVector.h" 34 #include "dMatrix.h" 35 #include "dSparse.h" 36 #include "lo-error.h" 37 #include "oct-locbuf.h" 38 #include "oct-refcount.h" 39 #include "oct-sparse.h" 40 #include "quit.h" 41 #include "sparse-qr.h" 42 43 namespace octave 44 { 45 namespace math 46 { 47 template <typename SPARSE_T> 48 class 49 cxsparse_types 50 { 51 }; 52 53 template <> 54 class 55 cxsparse_types<SparseMatrix> 56 { 57 public: 58 #if defined (HAVE_CXSPARSE) 59 typedef CXSPARSE_DNAME (s) symbolic_type; 60 typedef CXSPARSE_DNAME (n) numeric_type; 61 #else 62 typedef void symbolic_type; 63 typedef void numeric_type; 64 #endif 65 }; 66 67 template <> 68 class 69 cxsparse_types<SparseComplexMatrix> 70 { 71 public: 72 #if defined (HAVE_CXSPARSE) 73 typedef CXSPARSE_ZNAME (s) symbolic_type; 74 typedef CXSPARSE_ZNAME (n) numeric_type; 75 #else 76 typedef void symbolic_type; 77 typedef void numeric_type; 78 #endif 79 }; 80 81 template <typename SPARSE_T> 82 class sparse_qr<SPARSE_T>::sparse_qr_rep 83 { 84 public: 85 86 sparse_qr_rep (const SPARSE_T& a, int order); 87 88 // No copying! 89 90 sparse_qr_rep (const sparse_qr_rep&) = delete; 91 92 sparse_qr_rep& operator = (const sparse_qr_rep&) = delete; 93 94 ~sparse_qr_rep (void); 95 ok(void) const96 bool ok (void) const 97 { 98 #if defined (HAVE_CXSPARSE) 99 return (N && S); 100 #else 101 return false; 102 #endif 103 } 104 105 SPARSE_T V (void) const; 106 107 ColumnVector Pinv (void) const; 108 109 ColumnVector P (void) const; 110 111 SPARSE_T R (bool econ) const; 112 113 typename SPARSE_T::dense_matrix_type 114 C (const typename SPARSE_T::dense_matrix_type& b) const; 115 116 typename SPARSE_T::dense_matrix_type 117 Q (void) const; 118 119 refcount<octave_idx_type> count; 120 121 octave_idx_type nrows; 122 octave_idx_type ncols; 123 124 typename cxsparse_types<SPARSE_T>::symbolic_type *S; 125 typename cxsparse_types<SPARSE_T>::numeric_type *N; 126 127 template <typename RHS_T, typename RET_T> 128 RET_T 129 tall_solve (const RHS_T& b, octave_idx_type& info) const; 130 131 template <typename RHS_T, typename RET_T> 132 RET_T 133 wide_solve (const RHS_T& b, octave_idx_type& info) const; 134 }; 135 136 template <typename SPARSE_T> 137 ColumnVector Pinv(void) const138 sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const 139 { 140 #if defined (HAVE_CXSPARSE) 141 142 ColumnVector ret (N->L->m); 143 144 for (octave_idx_type i = 0; i < N->L->m; i++) 145 ret.xelem (i) = S->pinv[i]; 146 147 return ret; 148 149 #else 150 151 return ColumnVector (); 152 153 #endif 154 } 155 156 template <typename SPARSE_T> 157 ColumnVector P(void) const158 sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const 159 { 160 #if defined (HAVE_CXSPARSE) 161 162 ColumnVector ret (N->L->m); 163 164 for (octave_idx_type i = 0; i < N->L->m; i++) 165 ret.xelem (S->pinv[i]) = i; 166 167 return ret; 168 169 #else 170 171 return ColumnVector (); 172 173 #endif 174 } 175 176 // Specializations. 177 178 // Real-valued matrices. 179 180 template <> sparse_qr_rep(const SparseMatrix & a,int order)181 sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep 182 (const SparseMatrix& a, int order) 183 : count (1), nrows (a.rows ()), ncols (a.columns ()) 184 #if defined (HAVE_CXSPARSE) 185 , S (nullptr), N (nullptr) 186 #endif 187 { 188 #if defined (HAVE_CXSPARSE) 189 190 CXSPARSE_DNAME () A; 191 192 A.nzmax = a.nnz (); 193 A.m = nrows; 194 A.n = ncols; 195 // Cast away const on A, with full knowledge that CSparse won't touch it 196 // Prevents the methods below making a copy of the data. 197 A.p = const_cast<suitesparse_integer *> 198 (to_suitesparse_intptr (a.cidx ())); 199 A.i = const_cast<suitesparse_integer *> 200 (to_suitesparse_intptr (a.ridx ())); 201 A.x = const_cast<double *> (a.data ()); 202 A.nz = -1; 203 204 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 205 S = CXSPARSE_DNAME (_sqr) (order, &A, 1); 206 N = CXSPARSE_DNAME (_qr) (&A, S); 207 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 208 209 if (! N) 210 (*current_liboctave_error_handler) 211 ("sparse_qr: sparse matrix QR factorization filled"); 212 213 #else 214 215 octave_unused_parameter (order); 216 217 (*current_liboctave_error_handler) 218 ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); 219 220 #endif 221 } 222 223 template <> ~sparse_qr_rep(void)224 sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep (void) 225 { 226 #if defined (HAVE_CXSPARSE) 227 CXSPARSE_DNAME (_sfree) (S); 228 CXSPARSE_DNAME (_nfree) (N); 229 #endif 230 } 231 232 template <> 233 SparseMatrix V(void) const234 sparse_qr<SparseMatrix>::sparse_qr_rep::V (void) const 235 { 236 #if defined (HAVE_CXSPARSE) 237 238 // Drop zeros from V and sort 239 // FIXME: Is the double transpose to sort necessary? 240 241 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 242 CXSPARSE_DNAME (_dropzeros) (N->L); 243 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1); 244 CXSPARSE_DNAME (_spfree) (N->L); 245 N->L = CXSPARSE_DNAME (_transpose) (D, 1); 246 CXSPARSE_DNAME (_spfree) (D); 247 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 248 249 octave_idx_type nc = N->L->n; 250 octave_idx_type nz = N->L->nzmax; 251 SparseMatrix ret (N->L->m, nc, nz); 252 253 for (octave_idx_type j = 0; j < nc+1; j++) 254 ret.xcidx (j) = N->L->p[j]; 255 256 for (octave_idx_type j = 0; j < nz; j++) 257 { 258 ret.xridx (j) = N->L->i[j]; 259 ret.xdata (j) = N->L->x[j]; 260 } 261 262 return ret; 263 264 #else 265 266 return SparseMatrix (); 267 268 #endif 269 } 270 271 template <> 272 SparseMatrix R(bool econ) const273 sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const 274 { 275 #if defined (HAVE_CXSPARSE) 276 277 // Drop zeros from R and sort 278 // FIXME: Is the double transpose to sort necessary? 279 280 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 281 CXSPARSE_DNAME (_dropzeros) (N->U); 282 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1); 283 CXSPARSE_DNAME (_spfree) (N->U); 284 N->U = CXSPARSE_DNAME (_transpose) (D, 1); 285 CXSPARSE_DNAME (_spfree) (D); 286 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 287 288 octave_idx_type nc = N->U->n; 289 octave_idx_type nz = N->U->nzmax; 290 291 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); 292 293 for (octave_idx_type j = 0; j < nc+1; j++) 294 ret.xcidx (j) = N->U->p[j]; 295 296 for (octave_idx_type j = 0; j < nz; j++) 297 { 298 ret.xridx (j) = N->U->i[j]; 299 ret.xdata (j) = N->U->x[j]; 300 } 301 302 return ret; 303 304 #else 305 306 octave_unused_parameter (econ); 307 308 return SparseMatrix (); 309 310 #endif 311 } 312 313 template <> 314 Matrix C(const Matrix & b) const315 sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const 316 { 317 #if defined (HAVE_CXSPARSE) 318 319 octave_idx_type b_nr = b.rows (); 320 octave_idx_type b_nc = b.cols (); 321 322 octave_idx_type nc = N->L->n; 323 octave_idx_type nr = nrows; 324 325 const double *bvec = b.fortran_vec (); 326 327 Matrix ret (b_nr, b_nc); 328 double *vec = ret.fortran_vec (); 329 330 if (nr < 0 || nc < 0 || nr != b_nr) 331 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 332 333 if (nr == 0 || nc == 0 || b_nc == 0) 334 ret = Matrix (nc, b_nc, 0.0); 335 else 336 { 337 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 338 339 for (volatile octave_idx_type j = 0, idx = 0; 340 j < b_nc; 341 j++, idx += b_nr) 342 { 343 octave_quit (); 344 345 for (octave_idx_type i = nr; i < S->m2; i++) 346 buf[i] = 0.; 347 348 volatile octave_idx_type nm = (nr < nc ? nr : nc); 349 350 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 351 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); 352 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 353 354 for (volatile octave_idx_type i = 0; i < nm; i++) 355 { 356 octave_quit (); 357 358 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 359 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); 360 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 361 } 362 363 for (octave_idx_type i = 0; i < b_nr; i++) 364 vec[i+idx] = buf[i]; 365 } 366 } 367 368 return ret; 369 370 #else 371 372 octave_unused_parameter (b); 373 374 return Matrix (); 375 376 #endif 377 } 378 379 template <> 380 Matrix Q(void) const381 sparse_qr<SparseMatrix>::sparse_qr_rep::Q (void) const 382 { 383 #if defined (HAVE_CXSPARSE) 384 octave_idx_type nc = N->L->n; 385 octave_idx_type nr = nrows; 386 Matrix ret (nr, nr); 387 double *vec = ret.fortran_vec (); 388 389 if (nr < 0 || nc < 0) 390 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 391 392 if (nr == 0 || nc == 0) 393 ret = Matrix (nc, nr, 0.0); 394 else 395 { 396 OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1); 397 398 for (octave_idx_type i = 0; i < nr; i++) 399 bvec[i] = 0.; 400 401 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 402 403 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) 404 { 405 octave_quit (); 406 407 bvec[j] = 1.0; 408 for (octave_idx_type i = nr; i < S->m2; i++) 409 buf[i] = 0.; 410 411 volatile octave_idx_type nm = (nr < nc ? nr : nc); 412 413 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 414 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); 415 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 416 417 for (volatile octave_idx_type i = 0; i < nm; i++) 418 { 419 octave_quit (); 420 421 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 422 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); 423 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 424 } 425 426 for (octave_idx_type i = 0; i < nr; i++) 427 vec[i+idx] = buf[i]; 428 429 bvec[j] = 0.0; 430 } 431 } 432 433 return ret.transpose (); 434 435 #else 436 437 return Matrix (); 438 439 #endif 440 } 441 442 template <> 443 template <> 444 Matrix tall_solve(const MArray<double> & b,octave_idx_type & info) const445 sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix> 446 (const MArray<double>& b, octave_idx_type& info) const 447 { 448 info = -1; 449 450 #if defined (HAVE_CXSPARSE) 451 452 octave_idx_type nr = nrows; 453 octave_idx_type nc = ncols; 454 455 octave_idx_type b_nc = b.cols (); 456 octave_idx_type b_nr = b.rows (); 457 458 const double *bvec = b.data (); 459 460 Matrix x (nc, b_nc); 461 double *vec = x.fortran_vec (); 462 463 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 464 465 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 466 i++, idx+=nc, bidx+=b_nr) 467 { 468 octave_quit (); 469 470 for (octave_idx_type j = nr; j < S->m2; j++) 471 buf[j] = 0.; 472 473 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 474 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr); 475 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 476 477 for (volatile octave_idx_type j = 0; j < nc; j++) 478 { 479 octave_quit (); 480 481 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 482 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 483 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 484 } 485 486 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 487 CXSPARSE_DNAME (_usolve) (N->U, buf); 488 CXSPARSE_DNAME (_ipvec) (S->q, buf, vec + idx, nc); 489 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 490 } 491 492 info = 0; 493 494 return x; 495 496 #else 497 498 octave_unused_parameter (b); 499 500 return Matrix (); 501 502 #endif 503 } 504 505 template <> 506 template <> 507 Matrix wide_solve(const MArray<double> & b,octave_idx_type & info) const508 sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix> 509 (const MArray<double>& b, octave_idx_type& info) const 510 { 511 info = -1; 512 513 #if defined (HAVE_CXSPARSE) 514 515 // These are swapped because the original matrix was transposed in 516 // sparse_qr<SparseMatrix>::solve. 517 518 octave_idx_type nr = ncols; 519 octave_idx_type nc = nrows; 520 521 octave_idx_type b_nc = b.cols (); 522 octave_idx_type b_nr = b.rows (); 523 524 const double *bvec = b.data (); 525 526 Matrix x (nc, b_nc); 527 double *vec = x.fortran_vec (); 528 529 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 530 531 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 532 533 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 534 i++, idx+=nc, bidx+=b_nr) 535 { 536 octave_quit (); 537 538 for (octave_idx_type j = nr; j < nbuf; j++) 539 buf[j] = 0.; 540 541 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 542 CXSPARSE_DNAME (_pvec) (S->q, bvec + bidx, buf, nr); 543 CXSPARSE_DNAME (_utsolve) (N->U, buf); 544 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 545 546 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 547 { 548 octave_quit (); 549 550 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 551 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 552 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 553 } 554 555 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 556 CXSPARSE_DNAME (_pvec) (S->pinv, buf, vec + idx, nc); 557 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 558 } 559 560 info = 0; 561 562 return x; 563 564 #else 565 566 octave_unused_parameter (b); 567 568 return Matrix (); 569 570 #endif 571 } 572 573 template <> 574 template <> 575 SparseMatrix tall_solve(const SparseMatrix & b,octave_idx_type & info) const576 sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix> 577 (const SparseMatrix& b, octave_idx_type& info) const 578 { 579 info = -1; 580 581 #if defined (HAVE_CXSPARSE) 582 583 octave_idx_type nr = nrows; 584 octave_idx_type nc = ncols; 585 586 octave_idx_type b_nr = b.rows (); 587 octave_idx_type b_nc = b.cols (); 588 589 SparseMatrix x (nc, b_nc, b.nnz ()); 590 x.xcidx (0) = 0; 591 592 volatile octave_idx_type x_nz = b.nnz (); 593 volatile octave_idx_type ii = 0; 594 595 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 596 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 597 598 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 599 { 600 octave_quit (); 601 602 for (octave_idx_type j = 0; j < b_nr; j++) 603 Xx[j] = b.xelem (j,i); 604 605 for (octave_idx_type j = nr; j < S->m2; j++) 606 buf[j] = 0.; 607 608 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 609 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); 610 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 611 612 for (volatile octave_idx_type j = 0; j < nc; j++) 613 { 614 octave_quit (); 615 616 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 617 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 618 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 619 } 620 621 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 622 CXSPARSE_DNAME (_usolve) (N->U, buf); 623 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); 624 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 625 626 for (octave_idx_type j = 0; j < nc; j++) 627 { 628 double tmp = Xx[j]; 629 630 if (tmp != 0.0) 631 { 632 if (ii == x_nz) 633 { 634 // Resize the sparse matrix 635 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 636 sz = (sz > 10 ? sz : 10) + x_nz; 637 x.change_capacity (sz); 638 x_nz = sz; 639 } 640 641 x.xdata (ii) = tmp; 642 x.xridx (ii++) = j; 643 } 644 } 645 646 x.xcidx (i+1) = ii; 647 } 648 649 info = 0; 650 651 return x; 652 653 #else 654 655 octave_unused_parameter (b); 656 657 return SparseMatrix (); 658 659 #endif 660 } 661 662 template <> 663 template <> 664 SparseMatrix wide_solve(const SparseMatrix & b,octave_idx_type & info) const665 sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseMatrix> 666 (const SparseMatrix& b, octave_idx_type& info) const 667 { 668 info = -1; 669 670 #if defined (HAVE_CXSPARSE) 671 672 // These are swapped because the original matrix was transposed in 673 // sparse_qr<SparseMatrix>::solve. 674 675 octave_idx_type nr = ncols; 676 octave_idx_type nc = nrows; 677 678 octave_idx_type b_nr = b.rows (); 679 octave_idx_type b_nc = b.cols (); 680 681 SparseMatrix x (nc, b_nc, b.nnz ()); 682 x.xcidx (0) = 0; 683 684 volatile octave_idx_type x_nz = b.nnz (); 685 volatile octave_idx_type ii = 0; 686 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 687 688 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 689 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 690 691 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 692 { 693 octave_quit (); 694 695 for (octave_idx_type j = 0; j < b_nr; j++) 696 Xx[j] = b.xelem (j,i); 697 698 for (octave_idx_type j = nr; j < nbuf; j++) 699 buf[j] = 0.; 700 701 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 702 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); 703 CXSPARSE_DNAME (_utsolve) (N->U, buf); 704 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 705 706 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 707 { 708 octave_quit (); 709 710 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 711 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 712 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 713 } 714 715 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 716 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); 717 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 718 719 for (octave_idx_type j = 0; j < nc; j++) 720 { 721 double tmp = Xx[j]; 722 723 if (tmp != 0.0) 724 { 725 if (ii == x_nz) 726 { 727 // Resize the sparse matrix 728 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 729 sz = (sz > 10 ? sz : 10) + x_nz; 730 x.change_capacity (sz); 731 x_nz = sz; 732 } 733 734 x.xdata (ii) = tmp; 735 x.xridx (ii++) = j; 736 } 737 } 738 739 x.xcidx (i+1) = ii; 740 } 741 742 info = 0; 743 744 x.maybe_compress (); 745 746 return x; 747 748 #else 749 750 octave_unused_parameter (b); 751 752 return SparseMatrix (); 753 754 #endif 755 } 756 757 template <> 758 template <> 759 ComplexMatrix tall_solve(const MArray<Complex> & b,octave_idx_type & info) const760 sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix> 761 (const MArray<Complex>& b, octave_idx_type& info) const 762 { 763 info = -1; 764 765 #if defined (HAVE_CXSPARSE) 766 767 octave_idx_type nr = nrows; 768 octave_idx_type nc = ncols; 769 770 octave_idx_type b_nc = b.cols (); 771 octave_idx_type b_nr = b.rows (); 772 773 ComplexMatrix x (nc, b_nc); 774 Complex *vec = x.fortran_vec (); 775 776 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 777 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 778 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 779 780 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 781 { 782 octave_quit (); 783 784 for (octave_idx_type j = 0; j < b_nr; j++) 785 { 786 Complex c = b.xelem (j,i); 787 Xx[j] = c.real (); 788 Xz[j] = c.imag (); 789 } 790 791 for (octave_idx_type j = nr; j < S->m2; j++) 792 buf[j] = 0.; 793 794 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 795 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); 796 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 797 798 for (volatile octave_idx_type j = 0; j < nc; j++) 799 { 800 octave_quit (); 801 802 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 803 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 804 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 805 } 806 807 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 808 CXSPARSE_DNAME (_usolve) (N->U, buf); 809 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); 810 811 for (octave_idx_type j = nr; j < S->m2; j++) 812 buf[j] = 0.; 813 814 CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr); 815 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 816 817 for (volatile octave_idx_type j = 0; j < nc; j++) 818 { 819 octave_quit (); 820 821 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 822 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 823 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 824 } 825 826 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 827 CXSPARSE_DNAME (_usolve) (N->U, buf); 828 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc); 829 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 830 831 for (octave_idx_type j = 0; j < nc; j++) 832 vec[j+idx] = Complex (Xx[j], Xz[j]); 833 } 834 835 info = 0; 836 837 return x; 838 839 #else 840 841 octave_unused_parameter (b); 842 843 return ComplexMatrix (); 844 845 #endif 846 } 847 848 template <> 849 template <> 850 ComplexMatrix wide_solve(const MArray<Complex> & b,octave_idx_type & info) const851 sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix> 852 (const MArray<Complex>& b, octave_idx_type& info) const 853 { 854 info = -1; 855 856 #if defined (HAVE_CXSPARSE) 857 858 // These are swapped because the original matrix was transposed in 859 // sparse_qr<SparseMatrix>::solve. 860 861 octave_idx_type nr = ncols; 862 octave_idx_type nc = nrows; 863 864 octave_idx_type b_nc = b.cols (); 865 octave_idx_type b_nr = b.rows (); 866 867 ComplexMatrix x (nc, b_nc); 868 Complex *vec = x.fortran_vec (); 869 870 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 871 872 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 873 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 874 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 875 876 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 877 { 878 octave_quit (); 879 880 for (octave_idx_type j = 0; j < b_nr; j++) 881 { 882 Complex c = b.xelem (j,i); 883 Xx[j] = c.real (); 884 Xz[j] = c.imag (); 885 } 886 887 for (octave_idx_type j = nr; j < nbuf; j++) 888 buf[j] = 0.; 889 890 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 891 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); 892 CXSPARSE_DNAME (_utsolve) (N->U, buf); 893 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 894 895 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 896 { 897 octave_quit (); 898 899 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 900 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 901 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 902 } 903 904 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 905 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); 906 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 907 908 for (octave_idx_type j = nr; j < nbuf; j++) 909 buf[j] = 0.; 910 911 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 912 CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr); 913 CXSPARSE_DNAME (_utsolve) (N->U, buf); 914 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 915 916 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 917 { 918 octave_quit (); 919 920 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 921 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 922 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 923 } 924 925 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 926 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc); 927 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 928 929 for (octave_idx_type j = 0; j < nc; j++) 930 vec[j+idx] = Complex (Xx[j], Xz[j]); 931 } 932 933 info = 0; 934 935 return x; 936 937 #else 938 939 octave_unused_parameter (b); 940 941 return ComplexMatrix (); 942 943 #endif 944 } 945 946 // Complex-valued matrices. 947 948 template <> sparse_qr_rep(const SparseComplexMatrix & a,int order)949 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep 950 (const SparseComplexMatrix& a, int order) 951 : count (1), nrows (a.rows ()), ncols (a.columns ()) 952 #if defined (HAVE_CXSPARSE) 953 , S (nullptr), N (nullptr) 954 #endif 955 { 956 #if defined (HAVE_CXSPARSE) 957 958 CXSPARSE_ZNAME () A; 959 960 A.nzmax = a.nnz (); 961 A.m = nrows; 962 A.n = ncols; 963 // Cast away const on A, with full knowledge that CSparse won't touch it 964 // Prevents the methods below making a copy of the data. 965 A.p = const_cast<suitesparse_integer *> 966 (to_suitesparse_intptr (a.cidx ())); 967 A.i = const_cast<suitesparse_integer *> 968 (to_suitesparse_intptr (a.ridx ())); 969 A.x = const_cast<cs_complex_t *> 970 (reinterpret_cast<const cs_complex_t *> (a.data ())); 971 A.nz = -1; 972 973 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 974 S = CXSPARSE_ZNAME (_sqr) (order, &A, 1); 975 N = CXSPARSE_ZNAME (_qr) (&A, S); 976 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 977 978 if (! N) 979 (*current_liboctave_error_handler) 980 ("sparse_qr: sparse matrix QR factorization filled"); 981 982 #else 983 984 octave_unused_parameter (order); 985 986 (*current_liboctave_error_handler) 987 ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); 988 989 #endif 990 } 991 992 template <> ~sparse_qr_rep(void)993 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep (void) 994 { 995 #if defined (HAVE_CXSPARSE) 996 CXSPARSE_ZNAME (_sfree) (S); 997 CXSPARSE_ZNAME (_nfree) (N); 998 #endif 999 } 1000 1001 template <> 1002 SparseComplexMatrix V(void) const1003 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V (void) const 1004 { 1005 #if defined (HAVE_CXSPARSE) 1006 // Drop zeros from V and sort 1007 // FIXME: Is the double transpose to sort necessary? 1008 1009 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1010 CXSPARSE_ZNAME (_dropzeros) (N->L); 1011 CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1); 1012 CXSPARSE_ZNAME (_spfree) (N->L); 1013 N->L = CXSPARSE_ZNAME (_transpose) (D, 1); 1014 CXSPARSE_ZNAME (_spfree) (D); 1015 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1016 1017 octave_idx_type nc = N->L->n; 1018 octave_idx_type nz = N->L->nzmax; 1019 SparseComplexMatrix ret (N->L->m, nc, nz); 1020 1021 for (octave_idx_type j = 0; j < nc+1; j++) 1022 ret.xcidx (j) = N->L->p[j]; 1023 1024 for (octave_idx_type j = 0; j < nz; j++) 1025 { 1026 ret.xridx (j) = N->L->i[j]; 1027 ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j]; 1028 } 1029 1030 return ret; 1031 1032 #else 1033 1034 return SparseComplexMatrix (); 1035 1036 #endif 1037 } 1038 1039 template <> 1040 SparseComplexMatrix R(bool econ) const1041 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const 1042 { 1043 #if defined (HAVE_CXSPARSE) 1044 // Drop zeros from R and sort 1045 // FIXME: Is the double transpose to sort necessary? 1046 1047 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1048 CXSPARSE_ZNAME (_dropzeros) (N->U); 1049 CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1); 1050 CXSPARSE_ZNAME (_spfree) (N->U); 1051 N->U = CXSPARSE_ZNAME (_transpose) (D, 1); 1052 CXSPARSE_ZNAME (_spfree) (D); 1053 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1054 1055 octave_idx_type nc = N->U->n; 1056 octave_idx_type nz = N->U->nzmax; 1057 1058 SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), 1059 nc, nz); 1060 1061 for (octave_idx_type j = 0; j < nc+1; j++) 1062 ret.xcidx (j) = N->U->p[j]; 1063 1064 for (octave_idx_type j = 0; j < nz; j++) 1065 { 1066 ret.xridx (j) = N->U->i[j]; 1067 ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j]; 1068 } 1069 1070 return ret; 1071 1072 #else 1073 1074 octave_unused_parameter (econ); 1075 1076 return SparseComplexMatrix (); 1077 1078 #endif 1079 } 1080 1081 template <> 1082 ComplexMatrix C(const ComplexMatrix & b) const1083 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const 1084 { 1085 #if defined (HAVE_CXSPARSE) 1086 octave_idx_type b_nr = b.rows (); 1087 octave_idx_type b_nc = b.cols (); 1088 octave_idx_type nc = N->L->n; 1089 octave_idx_type nr = nrows; 1090 const cs_complex_t *bvec 1091 = reinterpret_cast<const cs_complex_t *> (b.fortran_vec ()); 1092 ComplexMatrix ret (b_nr, b_nc); 1093 Complex *vec = ret.fortran_vec (); 1094 1095 if (nr < 0 || nc < 0 || nr != b_nr) 1096 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 1097 1098 if (nr == 0 || nc == 0 || b_nc == 0) 1099 ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); 1100 else 1101 { 1102 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); 1103 1104 for (volatile octave_idx_type j = 0, idx = 0; 1105 j < b_nc; 1106 j++, idx += b_nr) 1107 { 1108 octave_quit (); 1109 1110 volatile octave_idx_type nm = (nr < nc ? nr : nc); 1111 1112 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1113 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx, 1114 reinterpret_cast<cs_complex_t *> (buf), 1115 b_nr); 1116 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1117 1118 for (volatile octave_idx_type i = 0; i < nm; i++) 1119 { 1120 octave_quit (); 1121 1122 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1123 CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], 1124 reinterpret_cast<cs_complex_t *> (buf)); 1125 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1126 } 1127 1128 for (octave_idx_type i = 0; i < b_nr; i++) 1129 vec[i+idx] = buf[i]; 1130 } 1131 } 1132 1133 return ret; 1134 1135 #else 1136 1137 octave_unused_parameter (b); 1138 1139 return ComplexMatrix (); 1140 1141 #endif 1142 } 1143 1144 template <> 1145 ComplexMatrix Q(void) const1146 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const 1147 { 1148 #if defined (HAVE_CXSPARSE) 1149 octave_idx_type nc = N->L->n; 1150 octave_idx_type nr = nrows; 1151 ComplexMatrix ret (nr, nr); 1152 Complex *vec = ret.fortran_vec (); 1153 1154 if (nr < 0 || nc < 0) 1155 (*current_liboctave_error_handler) ("matrix dimension mismatch"); 1156 1157 if (nr == 0 || nc == 0) 1158 ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0)); 1159 else 1160 { 1161 OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr); 1162 1163 for (octave_idx_type i = 0; i < nr; i++) 1164 bvec[i] = cs_complex_t (0.0, 0.0); 1165 1166 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); 1167 1168 for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) 1169 { 1170 octave_quit (); 1171 1172 bvec[j] = cs_complex_t (1.0, 0.0); 1173 1174 volatile octave_idx_type nm = (nr < nc ? nr : nc); 1175 1176 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1177 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec, 1178 reinterpret_cast<cs_complex_t *> (buf), 1179 nr); 1180 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1181 1182 for (volatile octave_idx_type i = 0; i < nm; i++) 1183 { 1184 octave_quit (); 1185 1186 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1187 CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], 1188 reinterpret_cast<cs_complex_t *> (buf)); 1189 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1190 } 1191 1192 for (octave_idx_type i = 0; i < nr; i++) 1193 vec[i+idx] = buf[i]; 1194 1195 bvec[j] = cs_complex_t (0.0, 0.0); 1196 } 1197 } 1198 1199 return ret.hermitian (); 1200 1201 #else 1202 1203 return ComplexMatrix (); 1204 1205 #endif 1206 } 1207 1208 template <> 1209 template <> 1210 SparseComplexMatrix tall_solve(const SparseComplexMatrix & b,octave_idx_type & info) const1211 sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, 1212 SparseComplexMatrix> 1213 (const SparseComplexMatrix& b, octave_idx_type& info) const 1214 { 1215 info = -1; 1216 1217 #if defined (HAVE_CXSPARSE) 1218 1219 octave_idx_type nr = nrows; 1220 octave_idx_type nc = ncols; 1221 1222 octave_idx_type b_nr = b.rows (); 1223 octave_idx_type b_nc = b.cols (); 1224 1225 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 1226 x.xcidx (0) = 0; 1227 1228 volatile octave_idx_type x_nz = b.nnz (); 1229 volatile octave_idx_type ii = 0; 1230 1231 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 1232 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 1233 OCTAVE_LOCAL_BUFFER (double, buf, S->m2); 1234 1235 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1236 { 1237 octave_quit (); 1238 1239 for (octave_idx_type j = 0; j < b_nr; j++) 1240 { 1241 Complex c = b.xelem (j,i); 1242 Xx[j] = c.real (); 1243 Xz[j] = c.imag (); 1244 } 1245 1246 for (octave_idx_type j = nr; j < S->m2; j++) 1247 buf[j] = 0.; 1248 1249 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1250 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); 1251 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1252 1253 for (volatile octave_idx_type j = 0; j < nc; j++) 1254 { 1255 octave_quit (); 1256 1257 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1258 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 1259 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1260 } 1261 1262 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1263 CXSPARSE_DNAME (_usolve) (N->U, buf); 1264 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); 1265 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1266 1267 for (octave_idx_type j = nr; j < S->m2; j++) 1268 buf[j] = 0.; 1269 1270 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1271 CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr); 1272 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1273 1274 for (volatile octave_idx_type j = 0; j < nc; j++) 1275 { 1276 octave_quit (); 1277 1278 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1279 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 1280 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1281 } 1282 1283 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1284 CXSPARSE_DNAME (_usolve) (N->U, buf); 1285 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc); 1286 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1287 1288 for (octave_idx_type j = 0; j < nc; j++) 1289 { 1290 Complex tmp = Complex (Xx[j], Xz[j]); 1291 1292 if (tmp != 0.0) 1293 { 1294 if (ii == x_nz) 1295 { 1296 // Resize the sparse matrix 1297 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 1298 sz = (sz > 10 ? sz : 10) + x_nz; 1299 x.change_capacity (sz); 1300 x_nz = sz; 1301 } 1302 1303 x.xdata (ii) = tmp; 1304 x.xridx (ii++) = j; 1305 } 1306 } 1307 1308 x.xcidx (i+1) = ii; 1309 } 1310 1311 info = 0; 1312 1313 return x; 1314 1315 #else 1316 1317 octave_unused_parameter (b); 1318 1319 return SparseComplexMatrix (); 1320 1321 #endif 1322 } 1323 1324 template <> 1325 template <> 1326 SparseComplexMatrix wide_solve(const SparseComplexMatrix & b,octave_idx_type & info) const1327 sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, 1328 SparseComplexMatrix> 1329 (const SparseComplexMatrix& b, octave_idx_type& info) const 1330 { 1331 info = -1; 1332 1333 #if defined (HAVE_CXSPARSE) 1334 1335 // These are swapped because the original matrix was transposed in 1336 // sparse_qr<SparseMatrix>::solve. 1337 1338 octave_idx_type nr = ncols; 1339 octave_idx_type nc = nrows; 1340 1341 octave_idx_type b_nr = b.rows (); 1342 octave_idx_type b_nc = b.cols (); 1343 1344 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 1345 x.xcidx (0) = 0; 1346 1347 volatile octave_idx_type x_nz = b.nnz (); 1348 volatile octave_idx_type ii = 0; 1349 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 1350 1351 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); 1352 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); 1353 OCTAVE_LOCAL_BUFFER (double, buf, nbuf); 1354 1355 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1356 { 1357 octave_quit (); 1358 1359 for (octave_idx_type j = 0; j < b_nr; j++) 1360 { 1361 Complex c = b.xelem (j,i); 1362 Xx[j] = c.real (); 1363 Xz[j] = c.imag (); 1364 } 1365 1366 for (octave_idx_type j = nr; j < nbuf; j++) 1367 buf[j] = 0.; 1368 1369 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1370 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); 1371 CXSPARSE_DNAME (_utsolve) (N->U, buf); 1372 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1373 1374 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 1375 { 1376 octave_quit (); 1377 1378 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1379 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 1380 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1381 } 1382 1383 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1384 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); 1385 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1386 1387 for (octave_idx_type j = nr; j < nbuf; j++) 1388 buf[j] = 0.; 1389 1390 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1391 CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr); 1392 CXSPARSE_DNAME (_utsolve) (N->U, buf); 1393 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1394 1395 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 1396 { 1397 octave_quit (); 1398 1399 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1400 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); 1401 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1402 } 1403 1404 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1405 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc); 1406 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1407 1408 for (octave_idx_type j = 0; j < nc; j++) 1409 { 1410 Complex tmp = Complex (Xx[j], Xz[j]); 1411 1412 if (tmp != 0.0) 1413 { 1414 if (ii == x_nz) 1415 { 1416 // Resize the sparse matrix 1417 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 1418 sz = (sz > 10 ? sz : 10) + x_nz; 1419 x.change_capacity (sz); 1420 x_nz = sz; 1421 } 1422 1423 x.xdata (ii) = tmp; 1424 x.xridx (ii++) = j; 1425 } 1426 } 1427 1428 x.xcidx (i+1) = ii; 1429 } 1430 1431 info = 0; 1432 1433 x.maybe_compress (); 1434 1435 return x; 1436 1437 #else 1438 1439 octave_unused_parameter (b); 1440 1441 return SparseComplexMatrix (); 1442 1443 #endif 1444 } 1445 1446 template <> 1447 template <> 1448 ComplexMatrix tall_solve(const MArray<double> & b,octave_idx_type & info) const1449 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>, 1450 ComplexMatrix> 1451 (const MArray<double>& b, octave_idx_type& info) const 1452 { 1453 info = -1; 1454 1455 #if defined (HAVE_CXSPARSE) 1456 1457 octave_idx_type nr = nrows; 1458 octave_idx_type nc = ncols; 1459 1460 octave_idx_type b_nc = b.cols (); 1461 octave_idx_type b_nr = b.rows (); 1462 1463 ComplexMatrix x (nc, b_nc); 1464 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); 1465 1466 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); 1467 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); 1468 1469 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1470 { 1471 octave_quit (); 1472 1473 for (octave_idx_type j = 0; j < b_nr; j++) 1474 Xx[j] = b.xelem (j,i); 1475 1476 for (octave_idx_type j = nr; j < S->m2; j++) 1477 buf[j] = cs_complex_t (0.0, 0.0); 1478 1479 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1480 CXSPARSE_ZNAME (_ipvec) (S->pinv, 1481 reinterpret_cast<cs_complex_t *>(Xx), 1482 buf, nr); 1483 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1484 1485 for (volatile octave_idx_type j = 0; j < nc; j++) 1486 { 1487 octave_quit (); 1488 1489 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1490 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); 1491 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1492 } 1493 1494 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1495 CXSPARSE_ZNAME (_usolve) (N->U, buf); 1496 CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc); 1497 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1498 } 1499 1500 info = 0; 1501 1502 return x; 1503 1504 #else 1505 1506 octave_unused_parameter (b); 1507 1508 return ComplexMatrix (); 1509 1510 #endif 1511 } 1512 1513 template <> 1514 template <> 1515 ComplexMatrix wide_solve(const MArray<double> & b,octave_idx_type & info) const1516 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>, 1517 ComplexMatrix> 1518 (const MArray<double>& b, octave_idx_type& info) const 1519 { 1520 info = -1; 1521 1522 #if defined (HAVE_CXSPARSE) 1523 1524 // These are swapped because the original matrix was transposed in 1525 // sparse_qr<SparseComplexMatrix>::solve. 1526 1527 octave_idx_type nr = ncols; 1528 octave_idx_type nc = nrows; 1529 1530 octave_idx_type b_nc = b.cols (); 1531 octave_idx_type b_nr = b.rows (); 1532 1533 ComplexMatrix x (nc, b_nc); 1534 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); 1535 1536 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 1537 1538 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); 1539 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); 1540 OCTAVE_LOCAL_BUFFER (double, B, nr); 1541 1542 for (octave_idx_type i = 0; i < nr; i++) 1543 B[i] = N->B[i]; 1544 1545 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1546 { 1547 octave_quit (); 1548 1549 for (octave_idx_type j = 0; j < b_nr; j++) 1550 Xx[j] = b.xelem (j,i); 1551 1552 for (octave_idx_type j = nr; j < nbuf; j++) 1553 buf[j] = cs_complex_t (0.0, 0.0); 1554 1555 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1556 CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), 1557 buf, nr); 1558 CXSPARSE_ZNAME (_utsolve) (N->U, buf); 1559 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1560 1561 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 1562 { 1563 octave_quit (); 1564 1565 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1566 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); 1567 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1568 } 1569 1570 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1571 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc); 1572 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1573 } 1574 1575 info = 0; 1576 1577 return x; 1578 1579 #else 1580 1581 octave_unused_parameter (b); 1582 1583 return ComplexMatrix (); 1584 1585 #endif 1586 } 1587 1588 template <> 1589 template <> 1590 SparseComplexMatrix tall_solve(const SparseMatrix & b,octave_idx_type & info) const1591 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, 1592 SparseComplexMatrix> 1593 (const SparseMatrix& b, octave_idx_type& info) const 1594 { 1595 info = -1; 1596 1597 #if defined (HAVE_CXSPARSE) 1598 1599 octave_idx_type nr = nrows; 1600 octave_idx_type nc = ncols; 1601 1602 octave_idx_type b_nc = b.cols (); 1603 octave_idx_type b_nr = b.rows (); 1604 1605 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 1606 x.xcidx (0) = 0; 1607 1608 volatile octave_idx_type x_nz = b.nnz (); 1609 volatile octave_idx_type ii = 0; 1610 1611 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 1612 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); 1613 1614 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1615 { 1616 octave_quit (); 1617 1618 for (octave_idx_type j = 0; j < b_nr; j++) 1619 Xx[j] = b.xelem (j,i); 1620 1621 for (octave_idx_type j = nr; j < S->m2; j++) 1622 buf[j] = cs_complex_t (0.0, 0.0); 1623 1624 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1625 CXSPARSE_ZNAME (_ipvec) (S->pinv, 1626 reinterpret_cast<cs_complex_t *>(Xx), 1627 buf, nr); 1628 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1629 1630 for (volatile octave_idx_type j = 0; j < nc; j++) 1631 { 1632 octave_quit (); 1633 1634 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1635 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); 1636 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1637 } 1638 1639 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1640 CXSPARSE_ZNAME (_usolve) (N->U, buf); 1641 CXSPARSE_ZNAME (_ipvec) (S->q, buf, 1642 reinterpret_cast<cs_complex_t *> (Xx), 1643 nc); 1644 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1645 1646 for (octave_idx_type j = 0; j < nc; j++) 1647 { 1648 Complex tmp = Xx[j]; 1649 1650 if (tmp != 0.0) 1651 { 1652 if (ii == x_nz) 1653 { 1654 // Resize the sparse matrix 1655 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 1656 sz = (sz > 10 ? sz : 10) + x_nz; 1657 x.change_capacity (sz); 1658 x_nz = sz; 1659 } 1660 1661 x.xdata (ii) = tmp; 1662 x.xridx (ii++) = j; 1663 } 1664 } 1665 1666 x.xcidx (i+1) = ii; 1667 } 1668 1669 info = 0; 1670 1671 x.maybe_compress (); 1672 1673 return x; 1674 1675 #else 1676 1677 octave_unused_parameter (b); 1678 1679 return SparseComplexMatrix (); 1680 1681 #endif 1682 } 1683 1684 template <> 1685 template <> 1686 SparseComplexMatrix wide_solve(const SparseMatrix & b,octave_idx_type & info) const1687 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, 1688 SparseComplexMatrix> 1689 (const SparseMatrix& b, octave_idx_type& info) const 1690 { 1691 info = -1; 1692 1693 #if defined (HAVE_CXSPARSE) 1694 1695 // These are swapped because the original matrix was transposed in 1696 // sparse_qr<SparseComplexMatrix>::solve. 1697 1698 octave_idx_type nr = ncols; 1699 octave_idx_type nc = nrows; 1700 1701 octave_idx_type b_nc = b.cols (); 1702 octave_idx_type b_nr = b.rows (); 1703 1704 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 1705 x.xcidx (0) = 0; 1706 1707 volatile octave_idx_type x_nz = b.nnz (); 1708 volatile octave_idx_type ii = 0; 1709 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 1710 1711 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 1712 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); 1713 OCTAVE_LOCAL_BUFFER (double, B, nr); 1714 1715 for (octave_idx_type i = 0; i < nr; i++) 1716 B[i] = N->B[i]; 1717 1718 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1719 { 1720 octave_quit (); 1721 1722 for (octave_idx_type j = 0; j < b_nr; j++) 1723 Xx[j] = b.xelem (j,i); 1724 1725 for (octave_idx_type j = nr; j < nbuf; j++) 1726 buf[j] = cs_complex_t (0.0, 0.0); 1727 1728 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1729 CXSPARSE_ZNAME (_pvec) (S->q, 1730 reinterpret_cast<cs_complex_t *>(Xx), 1731 buf, nr); 1732 CXSPARSE_ZNAME (_utsolve) (N->U, buf); 1733 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1734 1735 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 1736 { 1737 octave_quit (); 1738 1739 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1740 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); 1741 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1742 } 1743 1744 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1745 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, 1746 reinterpret_cast<cs_complex_t *> (Xx), 1747 nc); 1748 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1749 1750 for (octave_idx_type j = 0; j < nc; j++) 1751 { 1752 Complex tmp = Xx[j]; 1753 1754 if (tmp != 0.0) 1755 { 1756 if (ii == x_nz) 1757 { 1758 // Resize the sparse matrix 1759 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 1760 sz = (sz > 10 ? sz : 10) + x_nz; 1761 x.change_capacity (sz); 1762 x_nz = sz; 1763 } 1764 1765 x.xdata (ii) = tmp; 1766 x.xridx (ii++) = j; 1767 } 1768 } 1769 1770 x.xcidx (i+1) = ii; 1771 } 1772 1773 info = 0; 1774 1775 x.maybe_compress (); 1776 1777 return x; 1778 1779 #else 1780 1781 octave_unused_parameter (b); 1782 1783 return SparseComplexMatrix (); 1784 1785 #endif 1786 } 1787 1788 template <> 1789 template <> 1790 ComplexMatrix tall_solve(const MArray<Complex> & b,octave_idx_type & info) const1791 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, 1792 ComplexMatrix> 1793 (const MArray<Complex>& b, octave_idx_type& info) const 1794 { 1795 info = -1; 1796 1797 #if defined (HAVE_CXSPARSE) 1798 1799 octave_idx_type nr = nrows; 1800 octave_idx_type nc = ncols; 1801 1802 octave_idx_type b_nc = b.cols (); 1803 octave_idx_type b_nr = b.rows (); 1804 1805 const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *> 1806 (b.fortran_vec ()); 1807 1808 ComplexMatrix x (nc, b_nc); 1809 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> 1810 (x.fortran_vec ()); 1811 1812 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); 1813 1814 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 1815 i++, idx+=nc, bidx+=b_nr) 1816 { 1817 octave_quit (); 1818 1819 for (octave_idx_type j = nr; j < S->m2; j++) 1820 buf[j] = cs_complex_t (0.0, 0.0); 1821 1822 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1823 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr); 1824 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1825 1826 for (volatile octave_idx_type j = 0; j < nc; j++) 1827 { 1828 octave_quit (); 1829 1830 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1831 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); 1832 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1833 } 1834 1835 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1836 CXSPARSE_ZNAME (_usolve) (N->U, buf); 1837 CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc); 1838 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1839 } 1840 1841 info = 0; 1842 1843 return x; 1844 1845 #else 1846 1847 octave_unused_parameter (b); 1848 1849 return ComplexMatrix (); 1850 1851 #endif 1852 } 1853 1854 template <> 1855 template <> 1856 ComplexMatrix wide_solve(const MArray<Complex> & b,octave_idx_type & info) const1857 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, 1858 ComplexMatrix> 1859 (const MArray<Complex>& b, octave_idx_type& info) const 1860 { 1861 info = -1; 1862 1863 #if defined (HAVE_CXSPARSE) 1864 1865 // These are swapped because the original matrix was transposed in 1866 // sparse_qr<SparseComplexMatrix>::solve. 1867 1868 octave_idx_type nr = ncols; 1869 octave_idx_type nc = nrows; 1870 1871 octave_idx_type b_nc = b.cols (); 1872 octave_idx_type b_nr = b.rows (); 1873 1874 const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *> 1875 (b.fortran_vec ()); 1876 1877 ComplexMatrix x (nc, b_nc); 1878 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); 1879 1880 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 1881 1882 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); 1883 OCTAVE_LOCAL_BUFFER (double, B, nr); 1884 1885 for (octave_idx_type i = 0; i < nr; i++) 1886 B[i] = N->B[i]; 1887 1888 for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; 1889 i++, idx+=nc, bidx+=b_nr) 1890 { 1891 octave_quit (); 1892 1893 for (octave_idx_type j = nr; j < nbuf; j++) 1894 buf[j] = cs_complex_t (0.0, 0.0); 1895 1896 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1897 CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr); 1898 CXSPARSE_ZNAME (_utsolve) (N->U, buf); 1899 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1900 1901 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 1902 { 1903 octave_quit (); 1904 1905 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1906 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); 1907 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1908 } 1909 1910 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1911 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc); 1912 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1913 } 1914 1915 info = 0; 1916 1917 return x; 1918 1919 #else 1920 1921 octave_unused_parameter (b); 1922 1923 return ComplexMatrix (); 1924 1925 #endif 1926 } 1927 1928 template <> 1929 template <> 1930 SparseComplexMatrix tall_solve(const SparseComplexMatrix & b,octave_idx_type & info) const1931 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix> 1932 (const SparseComplexMatrix& b, octave_idx_type& info) const 1933 { 1934 info = -1; 1935 1936 #if defined (HAVE_CXSPARSE) 1937 1938 octave_idx_type nr = nrows; 1939 octave_idx_type nc = ncols; 1940 1941 octave_idx_type b_nc = b.cols (); 1942 octave_idx_type b_nr = b.rows (); 1943 1944 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 1945 x.xcidx (0) = 0; 1946 1947 volatile octave_idx_type x_nz = b.nnz (); 1948 volatile octave_idx_type ii = 0; 1949 1950 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 1951 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); 1952 1953 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 1954 { 1955 octave_quit (); 1956 1957 for (octave_idx_type j = 0; j < b_nr; j++) 1958 Xx[j] = b.xelem (j,i); 1959 1960 for (octave_idx_type j = nr; j < S->m2; j++) 1961 buf[j] = cs_complex_t (0.0, 0.0); 1962 1963 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1964 CXSPARSE_ZNAME (_ipvec) (S->pinv, 1965 reinterpret_cast<cs_complex_t *>(Xx), 1966 buf, nr); 1967 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1968 1969 for (volatile octave_idx_type j = 0; j < nc; j++) 1970 { 1971 octave_quit (); 1972 1973 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1974 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); 1975 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1976 } 1977 1978 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1979 CXSPARSE_ZNAME (_usolve) (N->U, buf); 1980 CXSPARSE_ZNAME (_ipvec) (S->q, buf, 1981 reinterpret_cast<cs_complex_t *> (Xx), 1982 nc); 1983 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 1984 1985 for (octave_idx_type j = 0; j < nc; j++) 1986 { 1987 Complex tmp = Xx[j]; 1988 1989 if (tmp != 0.0) 1990 { 1991 if (ii == x_nz) 1992 { 1993 // Resize the sparse matrix 1994 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 1995 sz = (sz > 10 ? sz : 10) + x_nz; 1996 x.change_capacity (sz); 1997 x_nz = sz; 1998 } 1999 2000 x.xdata (ii) = tmp; 2001 x.xridx (ii++) = j; 2002 } 2003 } 2004 2005 x.xcidx (i+1) = ii; 2006 } 2007 2008 info = 0; 2009 2010 x.maybe_compress (); 2011 2012 return x; 2013 2014 #else 2015 2016 octave_unused_parameter (b); 2017 2018 return SparseComplexMatrix (); 2019 2020 #endif 2021 } 2022 2023 template <> 2024 template <> 2025 SparseComplexMatrix wide_solve(const SparseComplexMatrix & b,octave_idx_type & info) const2026 sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix> 2027 (const SparseComplexMatrix& b, octave_idx_type& info) const 2028 { 2029 info = -1; 2030 2031 #if defined (HAVE_CXSPARSE) 2032 2033 // These are swapped because the original matrix was transposed in 2034 // sparse_qr<SparseComplexMatrix>::solve. 2035 2036 octave_idx_type nr = ncols; 2037 octave_idx_type nc = nrows; 2038 2039 octave_idx_type b_nc = b.cols (); 2040 octave_idx_type b_nr = b.rows (); 2041 2042 SparseComplexMatrix x (nc, b_nc, b.nnz ()); 2043 x.xcidx (0) = 0; 2044 2045 volatile octave_idx_type x_nz = b.nnz (); 2046 volatile octave_idx_type ii = 0; 2047 volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); 2048 2049 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); 2050 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); 2051 OCTAVE_LOCAL_BUFFER (double, B, nr); 2052 2053 for (octave_idx_type i = 0; i < nr; i++) 2054 B[i] = N->B[i]; 2055 2056 for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) 2057 { 2058 octave_quit (); 2059 2060 for (octave_idx_type j = 0; j < b_nr; j++) 2061 Xx[j] = b.xelem (j,i); 2062 2063 for (octave_idx_type j = nr; j < nbuf; j++) 2064 buf[j] = cs_complex_t (0.0, 0.0); 2065 2066 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2067 CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), 2068 buf, nr); 2069 CXSPARSE_ZNAME (_utsolve) (N->U, buf); 2070 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2071 2072 for (volatile octave_idx_type j = nr-1; j >= 0; j--) 2073 { 2074 octave_quit (); 2075 2076 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2077 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); 2078 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2079 } 2080 2081 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2082 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, 2083 reinterpret_cast<cs_complex_t *>(Xx), nc); 2084 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 2085 2086 for (octave_idx_type j = 0; j < nc; j++) 2087 { 2088 Complex tmp = Xx[j]; 2089 2090 if (tmp != 0.0) 2091 { 2092 if (ii == x_nz) 2093 { 2094 // Resize the sparse matrix 2095 octave_idx_type sz = x_nz * (b_nc - i) / b_nc; 2096 sz = (sz > 10 ? sz : 10) + x_nz; 2097 x.change_capacity (sz); 2098 x_nz = sz; 2099 } 2100 2101 x.xdata (ii) = tmp; 2102 x.xridx (ii++) = j; 2103 } 2104 } 2105 2106 x.xcidx (i+1) = ii; 2107 } 2108 2109 info = 0; 2110 2111 x.maybe_compress (); 2112 2113 return x; 2114 2115 #else 2116 2117 octave_unused_parameter (b); 2118 2119 return SparseComplexMatrix (); 2120 2121 #endif 2122 } 2123 2124 template <typename SPARSE_T> sparse_qr(void)2125 sparse_qr<SPARSE_T>::sparse_qr (void) 2126 : rep (new sparse_qr_rep (SPARSE_T (), 0)) 2127 { } 2128 2129 template <typename SPARSE_T> sparse_qr(const SPARSE_T & a,int order)2130 sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order) 2131 : rep (new sparse_qr_rep (a, order)) 2132 { } 2133 2134 template <typename SPARSE_T> sparse_qr(const sparse_qr<SPARSE_T> & a)2135 sparse_qr<SPARSE_T>::sparse_qr (const sparse_qr<SPARSE_T>& a) 2136 : rep (a.rep) 2137 { 2138 rep->count++; 2139 } 2140 2141 template <typename SPARSE_T> ~sparse_qr(void)2142 sparse_qr<SPARSE_T>::~sparse_qr (void) 2143 { 2144 if (--rep->count == 0) 2145 delete rep; 2146 } 2147 2148 template <typename SPARSE_T> 2149 sparse_qr<SPARSE_T>& operator =(const sparse_qr<SPARSE_T> & a)2150 sparse_qr<SPARSE_T>::operator = (const sparse_qr<SPARSE_T>& a) 2151 { 2152 if (this != &a) 2153 { 2154 if (--rep->count == 0) 2155 delete rep; 2156 2157 rep = a.rep; 2158 rep->count++; 2159 } 2160 2161 return *this; 2162 } 2163 2164 template <typename SPARSE_T> 2165 bool ok(void) const2166 sparse_qr<SPARSE_T>::ok (void) const 2167 { 2168 return rep->ok (); 2169 } 2170 2171 template <typename SPARSE_T> 2172 SPARSE_T V(void) const2173 sparse_qr<SPARSE_T>::V (void) const 2174 { 2175 return rep->V (); 2176 } 2177 2178 template <typename SPARSE_T> 2179 ColumnVector Pinv(void) const2180 sparse_qr<SPARSE_T>::Pinv (void) const 2181 { 2182 return rep->P (); 2183 } 2184 2185 template <typename SPARSE_T> 2186 ColumnVector P(void) const2187 sparse_qr<SPARSE_T>::P (void) const 2188 { 2189 return rep->P (); 2190 } 2191 2192 template <typename SPARSE_T> 2193 SPARSE_T R(bool econ) const2194 sparse_qr<SPARSE_T>::R (bool econ) const 2195 { 2196 return rep->R (econ); 2197 } 2198 2199 template <typename SPARSE_T> 2200 typename SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type & b) const2201 sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b) const 2202 { 2203 return rep->C (b); 2204 } 2205 2206 template <typename SPARSE_T> 2207 typename SPARSE_T::dense_matrix_type Q(void) const2208 sparse_qr<SPARSE_T>::Q (void) const 2209 { 2210 return rep->Q (); 2211 } 2212 2213 // FIXME: Why is the "order" of the QR calculation as used in the 2214 // CXSparse function sqr 3 for real matrices and 2 for complex? These 2215 // values seem to be required but there was no explanation in David 2216 // Bateman's original code. 2217 2218 template <typename SPARSE_T> 2219 class 2220 cxsparse_defaults 2221 { 2222 public: 2223 enum { order = -1 }; 2224 }; 2225 2226 template <> 2227 class 2228 cxsparse_defaults<SparseMatrix> 2229 { 2230 public: 2231 enum { order = 3 }; 2232 }; 2233 2234 template <> 2235 class 2236 cxsparse_defaults<SparseComplexMatrix> 2237 { 2238 public: 2239 enum { order = 2 }; 2240 }; 2241 2242 template <typename SPARSE_T> 2243 template <typename RHS_T, typename RET_T> 2244 RET_T solve(const SPARSE_T & a,const RHS_T & b,octave_idx_type & info)2245 sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b, 2246 octave_idx_type& info) 2247 { 2248 info = -1; 2249 2250 octave_idx_type nr = a.rows (); 2251 octave_idx_type nc = a.cols (); 2252 2253 octave_idx_type b_nc = b.cols (); 2254 octave_idx_type b_nr = b.rows (); 2255 2256 int order = cxsparse_defaults<SPARSE_T>::order; 2257 2258 if (nr < 0 || nc < 0 || nr != b_nr) 2259 (*current_liboctave_error_handler) 2260 ("matrix dimension mismatch in solution of minimum norm problem"); 2261 2262 if (nr == 0 || nc == 0 || b_nc == 0) 2263 { 2264 info = 0; 2265 2266 return RET_T (nc, b_nc, 0.0); 2267 } 2268 else if (nr >= nc) 2269 { 2270 sparse_qr<SPARSE_T> q (a, order); 2271 2272 return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T (); 2273 } 2274 else 2275 { 2276 sparse_qr<SPARSE_T> q (a.hermitian (), order); 2277 2278 return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T (); 2279 } 2280 } 2281 2282 template <typename SPARSE_T> 2283 template <typename RHS_T, typename RET_T> 2284 RET_T tall_solve(const RHS_T & b,octave_idx_type & info) const2285 sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const 2286 { 2287 return rep->template tall_solve<RHS_T, RET_T> (b, info); 2288 } 2289 2290 template <typename SPARSE_T> 2291 template <typename RHS_T, typename RET_T> 2292 RET_T wide_solve(const RHS_T & b,octave_idx_type & info) const2293 sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const 2294 { 2295 return rep->template wide_solve<RHS_T, RET_T> (b, info); 2296 } 2297 2298 Matrix qrsolve(const SparseMatrix & a,const MArray<double> & b,octave_idx_type & info)2299 qrsolve (const SparseMatrix& a, const MArray<double>& b, 2300 octave_idx_type& info) 2301 { 2302 return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b, 2303 info); 2304 } 2305 2306 SparseMatrix qrsolve(const SparseMatrix & a,const SparseMatrix & b,octave_idx_type & info)2307 qrsolve (const SparseMatrix& a, const SparseMatrix& b, 2308 octave_idx_type& info) 2309 { 2310 return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b, 2311 info); 2312 } 2313 2314 ComplexMatrix qrsolve(const SparseMatrix & a,const MArray<Complex> & b,octave_idx_type & info)2315 qrsolve (const SparseMatrix& a, const MArray<Complex>& b, 2316 octave_idx_type& info) 2317 { 2318 return sparse_qr<SparseMatrix>::solve<MArray<Complex>, 2319 ComplexMatrix> (a, b, info); 2320 } 2321 2322 SparseComplexMatrix qrsolve(const SparseMatrix & a,const SparseComplexMatrix & b,octave_idx_type & info)2323 qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b, 2324 octave_idx_type& info) 2325 { 2326 return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, 2327 SparseComplexMatrix> (a, b, info); 2328 } 2329 2330 ComplexMatrix qrsolve(const SparseComplexMatrix & a,const MArray<double> & b,octave_idx_type & info)2331 qrsolve (const SparseComplexMatrix& a, const MArray<double>& b, 2332 octave_idx_type& info) 2333 { 2334 return sparse_qr<SparseComplexMatrix>::solve<MArray<double>, 2335 ComplexMatrix> (a, b, info); 2336 } 2337 2338 SparseComplexMatrix qrsolve(const SparseComplexMatrix & a,const SparseMatrix & b,octave_idx_type & info)2339 qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b, 2340 octave_idx_type& info) 2341 { 2342 return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix, 2343 SparseComplexMatrix> 2344 (a, b, info); 2345 } 2346 2347 ComplexMatrix qrsolve(const SparseComplexMatrix & a,const MArray<Complex> & b,octave_idx_type & info)2348 qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b, 2349 octave_idx_type& info) 2350 { 2351 return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>, 2352 ComplexMatrix> (a, b, info); 2353 } 2354 2355 SparseComplexMatrix qrsolve(const SparseComplexMatrix & a,const SparseComplexMatrix & b,octave_idx_type & info)2356 qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b, 2357 octave_idx_type& info) 2358 { 2359 return sparse_qr<SparseComplexMatrix>::solve<SparseComplexMatrix, 2360 SparseComplexMatrix> 2361 (a, b, info); 2362 } 2363 2364 // Instantiations we need. 2365 2366 template class sparse_qr<SparseMatrix>; 2367 2368 template class sparse_qr<SparseComplexMatrix>; 2369 } 2370 } 2371