Mercurial > hg > octave-nkf

comparison liboctave/CSparse.cc @ 5610:9761b7d24e9e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2006-02-09 09:12:02 by dbateman]
author dbateman
date Thu, 09 Feb 2006 09:12:03 +0000
parents 2857357f9d3c
children 512d0d11ae39
comparison
equal deleted inserted replaced
5609:bf96b0f9dbd7 5610:9761b7d24e9e
41 #include "oct-spparms.h" 41 #include "oct-spparms.h"
42 #include "SparseCmplxLU.h" 42 #include "SparseCmplxLU.h"
43 #include "oct-sparse.h" 43 #include "oct-sparse.h"
44 #include "sparse-util.h" 44 #include "sparse-util.h"
45 #include "SparseCmplxCHOL.h" 45 #include "SparseCmplxCHOL.h"
46 #include "SparseCmplxQR.h"
46 47
47 #include "oct-sort.h" 48 #include "oct-sort.h"
48 49
49 // Fortran functions we call. 50 // Fortran functions we call.
50 extern "C" 51 extern "C"
608 return inverse (mattype, info, rcond, 0, 0); 609 return inverse (mattype, info, rcond, 0, 0);
609 } 610 }
610 611
611 SparseComplexMatrix 612 SparseComplexMatrix
612 SparseComplexMatrix::dinverse (SparseType &mattyp, octave_idx_type& info, 613 SparseComplexMatrix::dinverse (SparseType &mattyp, octave_idx_type& info,
613 double& rcond, const bool force, 614 double& rcond, const bool,
614 const bool calccond) const 615 const bool calccond) const
615 { 616 {
616 SparseComplexMatrix retval; 617 SparseComplexMatrix retval;
617 618
618 octave_idx_type nr = rows (); 619 octave_idx_type nr = rows ();
662 return retval; 663 return retval;
663 } 664 }
664 665
665 SparseComplexMatrix 666 SparseComplexMatrix
666 SparseComplexMatrix::tinverse (SparseType &mattyp, octave_idx_type& info, 667 SparseComplexMatrix::tinverse (SparseType &mattyp, octave_idx_type& info,
667 double& rcond, const bool force, 668 double& rcond, const bool,
668 const bool calccond) const 669 const bool calccond) const
669 { 670 {
670 SparseComplexMatrix retval; 671 SparseComplexMatrix retval;
671 672
672 octave_idx_type nr = rows (); 673 octave_idx_type nr = rows ();
900 return retval; 901 return retval;
901 } 902 }
902 903
903 SparseComplexMatrix 904 SparseComplexMatrix
904 SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info, 905 SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info,
905 double& rcond, int force, int calc_cond) const 906 double& rcond, int, int calc_cond) const
906 { 907 {
907 int typ = mattype.type (false); 908 int typ = mattype.type (false);
908 SparseComplexMatrix ret; 909 SparseComplexMatrix ret;
909 910
910 if (typ == SparseType::Unknown) 911 if (typ == SparseType::Unknown)
977 double rcond; 978 double rcond;
978 return determinant (info, rcond, 0); 979 return determinant (info, rcond, 0);
979 } 980 }
980 981
981 ComplexDET 982 ComplexDET
982 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int calc_cond) const 983 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int) const
983 { 984 {
984 ComplexDET retval; 985 ComplexDET retval;
985 #ifdef HAVE_UMFPACK 986 #ifdef HAVE_UMFPACK
986 987
987 octave_idx_type nr = rows (); 988 octave_idx_type nr = rows ();
4934 status = UMFPACK_ZNAME (numeric) (Ap, Ai, 4935 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
4935 X_CAST (const double *, Ax), NULL, 4936 X_CAST (const double *, Ax), NULL,
4936 Symbolic, &Numeric, control, info) ; 4937 Symbolic, &Numeric, control, info) ;
4937 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 4938 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ;
4938 4939
4939 #ifdef HAVE_LSSOLVE
4940 rcond = Info (UMFPACK_RCOND); 4940 rcond = Info (UMFPACK_RCOND);
4941 volatile double rcond_plus_one = rcond + 1.0; 4941 volatile double rcond_plus_one = rcond + 1.0;
4942 4942
4943 if (status == UMFPACK_WARNING_singular_matrix || 4943 if (status == UMFPACK_WARNING_singular_matrix ||
4944 rcond_plus_one == 1.0 || xisnan (rcond)) 4944 rcond_plus_one == 1.0 || xisnan (rcond))
4953 (*current_liboctave_error_handler) 4953 (*current_liboctave_error_handler)
4954 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", 4954 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
4955 rcond); 4955 rcond);
4956 4956
4957 } 4957 }
4958 else 4958 else if (status < 0)
4959 #endif
4960 if (status < 0)
4961 { 4959 {
4962 (*current_liboctave_error_handler) 4960 (*current_liboctave_error_handler)
4963 ("SparseComplexMatrix::solve numeric factorization failed"); 4961 ("SparseComplexMatrix::solve numeric factorization failed");
4964 4962
4965 UMFPACK_ZNAME (report_status) (control, status); 4963 UMFPACK_ZNAME (report_status) (control, status);
5115 else 5113 else
5116 (*current_liboctave_error_handler) 5114 (*current_liboctave_error_handler)
5117 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", 5115 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5118 rcond); 5116 rcond);
5119 5117
5120 #ifdef HAVE_LSSOLVE
5121 return retval; 5118 return retval;
5122 #endif
5123 } 5119 }
5124 5120
5125 cholmod_dense *X; 5121 cholmod_dense *X;
5126 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 5122 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5127 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); 5123 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5214 5210
5215 break; 5211 break;
5216 } 5212 }
5217 } 5213 }
5218 5214
5219 #ifndef HAVE_LSSOLVE
5220 rcond = Info (UMFPACK_RCOND);
5221 volatile double rcond_plus_one = rcond + 1.0;
5222
5223 if (status == UMFPACK_WARNING_singular_matrix ||
5224 rcond_plus_one == 1.0 || xisnan (rcond))
5225 {
5226 err = -2;
5227
5228 if (sing_handler)
5229 sing_handler (rcond);
5230 else
5231 (*current_liboctave_error_handler)
5232 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
5233 rcond);
5234
5235 }
5236 #endif
5237
5238 UMFPACK_ZNAME (report_info) (control, info); 5215 UMFPACK_ZNAME (report_info) (control, info);
5239 5216
5240 UMFPACK_ZNAME (free_numeric) (&Numeric); 5217 UMFPACK_ZNAME (free_numeric) (&Numeric);
5241 } 5218 }
5242 #else 5219 #else
5394 else 5371 else
5395 (*current_liboctave_error_handler) 5372 (*current_liboctave_error_handler)
5396 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", 5373 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5397 rcond); 5374 rcond);
5398 5375
5399 #ifdef HAVE_LSSOLVE
5400 return retval; 5376 return retval;
5401 #endif
5402 } 5377 }
5403 5378
5404 cholmod_sparse *X; 5379 cholmod_sparse *X;
5405 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 5380 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5406 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); 5381 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
5528 retval.xcidx(j+1) = ii; 5503 retval.xcidx(j+1) = ii;
5529 } 5504 }
5530 5505
5531 retval.maybe_compress (); 5506 retval.maybe_compress ();
5532 5507
5533 #ifndef HAVE_LSSOLVE
5534 rcond = Info (UMFPACK_RCOND);
5535 volatile double rcond_plus_one = rcond + 1.0;
5536
5537 if (status == UMFPACK_WARNING_singular_matrix ||
5538 rcond_plus_one == 1.0 || xisnan (rcond))
5539 {
5540 err = -2;
5541
5542 if (sing_handler)
5543 sing_handler (rcond);
5544 else
5545 (*current_liboctave_error_handler)
5546 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
5547 rcond);
5548
5549 }
5550 #endif
5551
5552 UMFPACK_ZNAME (report_info) (control, info); 5508 UMFPACK_ZNAME (report_info) (control, info);
5553 5509
5554 UMFPACK_ZNAME (free_numeric) (&Numeric); 5510 UMFPACK_ZNAME (free_numeric) (&Numeric);
5555 } 5511 }
5556 #else 5512 #else
5698 else 5654 else
5699 (*current_liboctave_error_handler) 5655 (*current_liboctave_error_handler)
5700 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", 5656 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5701 rcond); 5657 rcond);
5702 5658
5703 #ifdef HAVE_LSSOLVE
5704 return retval; 5659 return retval;
5705 #endif
5706 } 5660 }
5707 5661
5708 cholmod_dense *X; 5662 cholmod_dense *X;
5709 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 5663 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5710 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); 5664 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5774 err = -1; 5728 err = -1;
5775 5729
5776 break; 5730 break;
5777 } 5731 }
5778 } 5732 }
5779
5780 #ifndef HAVE_LSSOLVE
5781 rcond = Info (UMFPACK_RCOND);
5782 volatile double rcond_plus_one = rcond + 1.0;
5783
5784 if (status == UMFPACK_WARNING_singular_matrix ||
5785 rcond_plus_one == 1.0 || xisnan (rcond))
5786 {
5787 err = -2;
5788
5789 if (sing_handler)
5790 sing_handler (rcond);
5791 else
5792 (*current_liboctave_error_handler)
5793 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
5794 rcond);
5795
5796 }
5797 #endif
5798 5733
5799 UMFPACK_ZNAME (report_info) (control, info); 5734 UMFPACK_ZNAME (report_info) (control, info);
5800 5735
5801 UMFPACK_ZNAME (free_numeric) (&Numeric); 5736 UMFPACK_ZNAME (free_numeric) (&Numeric);
5802 } 5737 }
5955 else 5890 else
5956 (*current_liboctave_error_handler) 5891 (*current_liboctave_error_handler)
5957 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", 5892 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5958 rcond); 5893 rcond);
5959 5894
5960 #ifdef HAVE_LSSOLVE
5961 return retval; 5895 return retval;
5962 #endif
5963 } 5896 }
5964 5897
5965 cholmod_sparse *X; 5898 cholmod_sparse *X;
5966 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 5899 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5967 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); 5900 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6068 retval.xcidx(j+1) = ii; 6001 retval.xcidx(j+1) = ii;
6069 } 6002 }
6070 6003
6071 retval.maybe_compress (); 6004 retval.maybe_compress ();
6072 6005
6073 #ifndef HAVE_LSSOLVE
6074 rcond = Info (UMFPACK_RCOND); 6006 rcond = Info (UMFPACK_RCOND);
6075 volatile double rcond_plus_one = rcond + 1.0; 6007 volatile double rcond_plus_one = rcond + 1.0;
6076 6008
6077 if (status == UMFPACK_WARNING_singular_matrix || 6009 if (status == UMFPACK_WARNING_singular_matrix ||
6078 rcond_plus_one == 1.0 || xisnan (rcond)) 6010 rcond_plus_one == 1.0 || xisnan (rcond))
6085 (*current_liboctave_error_handler) 6017 (*current_liboctave_error_handler)
6086 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", 6018 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
6087 rcond); 6019 rcond);
6088 6020
6089 } 6021 }
6090 #endif
6091 6022
6092 UMFPACK_ZNAME (report_info) (control, info); 6023 UMFPACK_ZNAME (report_info) (control, info);
6093 6024
6094 UMFPACK_ZNAME (free_numeric) (&Numeric); 6025 UMFPACK_ZNAME (free_numeric) (&Numeric);
6095 } 6026 }
6578 octave_idx_type rank; 6509 octave_idx_type rank;
6579 return lssolve (b, info, rank); 6510 return lssolve (b, info, rank);
6580 } 6511 }
6581 6512
6582 ComplexMatrix 6513 ComplexMatrix
6583 SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const 6514 SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type&) const
6584 { 6515 {
6585 info = -1; 6516 return qrsolve (*this, b, info);
6586 (*current_liboctave_error_handler)
6587 ("SparseComplexMatrix::lssolve not implemented yet");
6588 return ComplexMatrix ();
6589 } 6517 }
6590 6518
6591 SparseComplexMatrix 6519 SparseComplexMatrix
6592 SparseComplexMatrix::lssolve (const SparseMatrix& b) const 6520 SparseComplexMatrix::lssolve (const SparseMatrix& b) const
6593 { 6521 {
6603 return lssolve (b, info, rank); 6531 return lssolve (b, info, rank);
6604 } 6532 }
6605 6533
6606 SparseComplexMatrix 6534 SparseComplexMatrix
6607 SparseComplexMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info, 6535 SparseComplexMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info,
6608 octave_idx_type& rank) const 6536 octave_idx_type&) const
6609 { 6537 {
6610 info = -1; 6538 return qrsolve (*this, b, info);
6611 (*current_liboctave_error_handler)
6612 ("SparseComplexMatrix::lssolve not implemented yet");
6613 return SparseComplexMatrix ();
6614 } 6539 }
6615 6540
6616 ComplexMatrix 6541 ComplexMatrix
6617 SparseComplexMatrix::lssolve (const ComplexMatrix& b) const 6542 SparseComplexMatrix::lssolve (const ComplexMatrix& b) const
6618 { 6543 {
6628 return lssolve (b, info, rank); 6553 return lssolve (b, info, rank);
6629 } 6554 }
6630 6555
6631 ComplexMatrix 6556 ComplexMatrix
6632 SparseComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, 6557 SparseComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info,
6633 octave_idx_type& rank) const 6558 octave_idx_type&) const
6634 { 6559 {
6635 info = -1; 6560 return qrsolve (*this, b, info);
6636 (*current_liboctave_error_handler)
6637 ("SparseComplexMatrix::lssolve not implemented yet");
6638 return ComplexMatrix ();
6639 } 6561 }
6640 6562
6641 SparseComplexMatrix 6563 SparseComplexMatrix
6642 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b) const 6564 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b) const
6643 { 6565 {
6653 return lssolve (b, info, rank); 6575 return lssolve (b, info, rank);
6654 } 6576 }
6655 6577
6656 SparseComplexMatrix 6578 SparseComplexMatrix
6657 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info, 6579 SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info,
6658 octave_idx_type& rank) const 6580 octave_idx_type&) const
6659 { 6581 {
6660 info = -1; 6582 return qrsolve (*this, b, info);
6661 (*current_liboctave_error_handler)
6662 ("SparseComplexMatrix::lssolve not implemented yet");
6663 return SparseComplexMatrix ();
6664 } 6583 }
6665 6584
6666 ComplexColumnVector 6585 ComplexColumnVector
6667 SparseComplexMatrix::lssolve (const ColumnVector& b) const 6586 SparseComplexMatrix::lssolve (const ColumnVector& b) const
6668 { 6587 {
6679 } 6598 }
6680 6599
6681 ComplexColumnVector 6600 ComplexColumnVector
6682 SparseComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const 6601 SparseComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const
6683 { 6602 {
6684 info = -1; 6603 Matrix tmp (b);
6685 (*current_liboctave_error_handler) 6604 return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0));
6686 ("SparseComplexMatrix::lssolve not implemented yet");
6687 return ComplexColumnVector ();
6688 } 6605 }
6689 6606
6690 ComplexColumnVector 6607 ComplexColumnVector
6691 SparseComplexMatrix::lssolve (const ComplexColumnVector& b) const 6608 SparseComplexMatrix::lssolve (const ComplexColumnVector& b) const
6692 { 6609 {
6704 6621
6705 ComplexColumnVector 6622 ComplexColumnVector
6706 SparseComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, 6623 SparseComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info,
6707 octave_idx_type& rank) const 6624 octave_idx_type& rank) const
6708 { 6625 {
6709 info = -1; 6626 ComplexMatrix tmp (b);
6710 (*current_liboctave_error_handler) 6627 return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0));
6711 ("SparseComplexMatrix::lssolve not implemented yet");
6712 return ComplexColumnVector ();
6713 } 6628 }
6714 6629
6715 // unary operations 6630 // unary operations
6716 SparseBoolMatrix 6631 SparseBoolMatrix
6717 SparseComplexMatrix::operator ! (void) const 6632 SparseComplexMatrix::operator ! (void) const