Mercurial > octave-nkf

comparison liboctave/CSparse.cc @ 5506:b4cfbb0ec8c4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-10-23 19:09:32 by dbateman]
author dbateman
date Sun, 23 Oct 2005 19:09:33 +0000
parents ed08548b9054
children 8c56849b1509
comparison
equal deleted inserted replaced
5505:17682e3fba2a 5506:b4cfbb0ec8c4
39 #include "boolSparse.h" 39 #include "boolSparse.h"
40 #include "dSparse.h" 40 #include "dSparse.h"
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"
45 #include "SparseCmplxCHOL.h"
44 46
45 // Fortran functions we call. 47 // Fortran functions we call.
46 extern "C" 48 extern "C"
47 { 49 {
48 F77_RET_T 50 F77_RET_T
583 SparseComplexMatrix 585 SparseComplexMatrix
584 SparseComplexMatrix::inverse (void) const 586 SparseComplexMatrix::inverse (void) const
585 { 587 {
586 octave_idx_type info; 588 octave_idx_type info;
587 double rcond; 589 double rcond;
588 return inverse (info, rcond, 0, 0); 590 SparseType mattype (*this);
591 return inverse (mattype, info, rcond, 0, 0);
589 } 592 }
590 593
591 SparseComplexMatrix 594 SparseComplexMatrix
592 SparseComplexMatrix::inverse (octave_idx_type& info) const 595 SparseComplexMatrix::inverse (SparseType& mattype) const
593 { 596 {
597 octave_idx_type info;
594 double rcond; 598 double rcond;
595 return inverse (info, rcond, 0, 0); 599 return inverse (mattype, info, rcond, 0, 0);
596 } 600 }
597 601
598 SparseComplexMatrix 602 SparseComplexMatrix
599 SparseComplexMatrix::inverse (octave_idx_type& info, double& rcond, int force, 603 SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info) const
600 int calc_cond) const 604 {
601 { 605 double rcond;
602 info = -1; 606 return inverse (mattype, info, rcond, 0, 0);
603 (*current_liboctave_error_handler) 607 }
604 ("SparseComplexMatrix::inverse not implemented yet"); 608
605 return SparseComplexMatrix (); 609 SparseComplexMatrix
610 SparseComplexMatrix::dinverse (SparseType &mattyp, octave_idx_type& info,
611 double& rcond, const bool force,
612 const bool calccond) const
613 {
614 SparseComplexMatrix retval;
615
616 octave_idx_type nr = rows ();
617 octave_idx_type nc = cols ();
618 info = 0;
619
620 if (nr == 0 || nc == 0 || nr != nc)
621 (*current_liboctave_error_handler) ("inverse requires square matrix");
622 else
623 {
624 // Print spparms("spumoni") info if requested
625 int typ = mattyp.type ();
626 mattyp.info ();
627
628 if (typ == SparseType::Diagonal ||
629 typ == SparseType::Permuted_Diagonal)
630 {
631 if (typ == SparseType::Permuted_Diagonal)
632 retval = transpose();
633 else
634 retval = *this;
635
636 // Force make_unique to be called
637 Complex *v = retval.data();
638
639 if (calccond)
640 {
641 double dmax = 0., dmin = octave_Inf;
642 for (octave_idx_type i = 0; i < nr; i++)
643 {
644 double tmp = std::abs(v[i]);
645 if (tmp > dmax)
646 dmax = tmp;
647 if (tmp < dmin)
648 dmin = tmp;
649 }
650 rcond = dmin / dmax;
651 }
652
653 for (octave_idx_type i = 0; i < nr; i++)
654 v[i] = 1.0 / v[i];
655 }
656 else
657 (*current_liboctave_error_handler) ("incorrect matrix type");
658 }
659
660 return retval;
661 }
662
663 SparseComplexMatrix
664 SparseComplexMatrix::tinverse (SparseType &mattyp, octave_idx_type& info,
665 double& rcond, const bool force,
666 const bool calccond) const
667 {
668 SparseComplexMatrix retval;
669
670 octave_idx_type nr = rows ();
671 octave_idx_type nc = cols ();
672 info = 0;
673
674 if (nr == 0 || nc == 0 || nr != nc)
675 (*current_liboctave_error_handler) ("inverse requires square matrix");
676 else
677 {
678 // Print spparms("spumoni") info if requested
679 int typ = mattyp.type ();
680 mattyp.info ();
681
682 if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper ||
683 typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
684 {
685 double anorm = 0.;
686 double ainvnorm = 0.;
687
688 if (calccond)
689 {
690 // Calculate the 1-norm of matrix for rcond calculation
691 for (octave_idx_type j = 0; j < nr; j++)
692 {
693 double atmp = 0.;
694 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
695 atmp += std::abs(data(i));
696 if (atmp > anorm)
697 anorm = atmp;
698 }
699 }
700
701 if (typ == SparseType::Upper || typ == SparseType::Lower)
702 {
703 octave_idx_type nz = nnz();
704 octave_idx_type cx = 0;
705 octave_idx_type nz2 = nz;
706 retval = SparseComplexMatrix (nr, nc, nz2);
707
708 for (octave_idx_type i = 0; i < nr; i++)
709 {
710 OCTAVE_QUIT;
711 // place the 1 in the identity position
712 octave_idx_type cx_colstart = cx;
713
714 if (cx == nz2)
715 {
716 nz2 *= 2;
717 retval.change_capacity (nz2);
718 }
719
720 retval.xcidx(i) = cx;
721 retval.xridx(cx) = i;
722 retval.xdata(cx) = 1.0;
723 cx++;
724
725 // iterate accross columns of input matrix
726 for (octave_idx_type j = i+1; j < nr; j++)
727 {
728 Complex v = 0.;
729 // iterate to calculate sum
730 octave_idx_type colXp = retval.xcidx(i);
731 octave_idx_type colUp = cidx(j);
732 octave_idx_type rpX, rpU;
733 do
734 {
735 OCTAVE_QUIT;
736 rpX = retval.xridx(colXp);
737 rpU = ridx(colUp);
738
739 if (rpX < rpU)
740 colXp++;
741 else if (rpX > rpU)
742 colUp++;
743 else
744 {
745 v -= retval.xdata(colXp) * data(colUp);
746 colXp++;
747 colUp++;
748 }
749 } while ((rpX<j) && (rpU<j) &&
750 (colXp<cx) && (colUp<nz));
751
752 // get A(m,m)
753 colUp = cidx(j+1) - 1;
754 Complex pivot = data(colUp);
755 if (pivot == 0.)
756 (*current_liboctave_error_handler)
757 ("division by zero");
758
759 if (v != 0.)
760 {
761 if (cx == nz2)
762 {
763 nz2 *= 2;
764 retval.change_capacity (nz2);
765 }
766
767 retval.xridx(cx) = j;
768 retval.xdata(cx) = v / pivot;
769 cx++;
770 }
771 }
772
773 // get A(m,m)
774 octave_idx_type colUp = cidx(i+1) - 1;
775 Complex pivot = data(colUp);
776 if (pivot == 0.)
777 (*current_liboctave_error_handler) ("division by zero");
778
779 if (pivot != 1.0)
780 for (octave_idx_type j = cx_colstart; j < cx; j++)
781 retval.xdata(j) /= pivot;
782 }
783 retval.xcidx(nr) = cx;
784 retval.maybe_compress ();
785 }
786 else
787 {
788 octave_idx_type nz = nnz();
789 octave_idx_type cx = 0;
790 octave_idx_type nz2 = nz;
791 retval = SparseComplexMatrix (nr, nc, nz2);
792
793 OCTAVE_LOCAL_BUFFER (Complex, work, nr);
794 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
795
796 octave_idx_type *perm = mattyp.triangular_perm();
797 if (typ == SparseType::Permuted_Upper)
798 {
799 for (octave_idx_type i = 0; i < nr; i++)
800 rperm[perm[i]] = i;
801 }
802 else
803 {
804 for (octave_idx_type i = 0; i < nr; i++)
805 rperm[i] = perm[i];
806 for (octave_idx_type i = 0; i < nr; i++)
807 perm[rperm[i]] = i;
808 }
809
810 for (octave_idx_type i = 0; i < nr; i++)
811 {
812 OCTAVE_QUIT;
813 octave_idx_type iidx = rperm[i];
814
815 for (octave_idx_type j = 0; j < nr; j++)
816 work[j] = 0.;
817
818 // place the 1 in the identity position
819 work[iidx] = 1.0;
820
821 // iterate accross columns of input matrix
822 for (octave_idx_type j = iidx+1; j < nr; j++)
823 {
824 Complex v = 0.;
825 octave_idx_type jidx = perm[j];
826 // iterate to calculate sum
827 for (octave_idx_type k = cidx(jidx);
828 k < cidx(jidx+1); k++)
829 {
830 OCTAVE_QUIT;
831 v -= work[ridx(k)] * data(k);
832 }
833
834 // get A(m,m)
835 Complex pivot = data(cidx(jidx+1) - 1);
836 if (pivot == 0.)
837 (*current_liboctave_error_handler)
838 ("division by zero");
839
840 work[j] = v / pivot;
841 }
842
843 // get A(m,m)
844 octave_idx_type colUp = cidx(perm[iidx]+1) - 1;
845 Complex pivot = data(colUp);
846 if (pivot == 0.)
847 (*current_liboctave_error_handler)
848 ("division by zero");
849
850 octave_idx_type new_cx = cx;
851 for (octave_idx_type j = iidx; j < nr; j++)
852 if (work[j] != 0.0)
853 {
854 new_cx++;
855 if (pivot != 1.0)
856 work[j] /= pivot;
857 }
858
859 if (cx < new_cx)
860 {
861 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
862 retval.change_capacity (nz2);
863 }
864
865 retval.xcidx(i) = cx;
866 for (octave_idx_type j = iidx; j < nr; j++)
867 if (work[j] != 0.)
868 {
869 retval.xridx(cx) = j;
870 retval.xdata(cx++) = work[j];
871 }
872 }
873
874 retval.xcidx(nr) = cx;
875 retval.maybe_compress ();
876 }
877
878 if (calccond)
879 {
880 // Calculate the 1-norm of inverse matrix for rcond calculation
881 for (octave_idx_type j = 0; j < nr; j++)
882 {
883 double atmp = 0.;
884 for (octave_idx_type i = retval.cidx(j);
885 i < retval.cidx(j+1); i++)
886 atmp += std::abs(retval.data(i));
887 if (atmp > ainvnorm)
888 ainvnorm = atmp;
889 }
890
891 rcond = 1. / ainvnorm / anorm;
892 }
893 }
894 else
895 (*current_liboctave_error_handler) ("incorrect matrix type");
896 }
897
898 return retval;
899 }
900
901 SparseComplexMatrix
902 SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info,
903 double& rcond, int force, int calc_cond) const
904 {
905 int typ = mattype.type (false);
906 SparseComplexMatrix ret;
907
908 if (typ == SparseType::Unknown)
909 typ = mattype.type (*this);
910
911 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
912 ret = dinverse (mattype, info, rcond, true, calc_cond);
913 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
914 ret = tinverse (mattype, info, rcond, true, calc_cond).transpose();
915 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
916 ret = transpose().tinverse (mattype, info, rcond, true, calc_cond);
917 else if (typ != SparseType::Rectangular)
918 {
919 if (mattype.is_hermitian())
920 {
921 SparseType tmp_typ (SparseType::Upper);
922 SparseComplexCHOL fact (*this, info, false);
923 rcond = fact.rcond();
924 if (info == 0)
925 {
926 double rcond2;
927 SparseMatrix Q = fact.Q();
928 SparseComplexMatrix InvL = fact.L().transpose().
929 tinverse(tmp_typ, info, rcond2, true, false);
930 ret = Q * InvL.hermitian() * InvL * Q.transpose();
931 }
932 else
933 {
934 // Matrix is either singular or not positive definite
935 mattype.mark_as_unsymmetric ();
936 typ = SparseType::Full;
937 }
938 }
939
940 if (!mattype.is_hermitian())
941 {
942 octave_idx_type n = rows();
943 ColumnVector Qinit(n);
944 for (octave_idx_type i = 0; i < n; i++)
945 Qinit(i) = i;
946
947 SparseType tmp_typ (SparseType::Upper);
948 SparseComplexLU fact (*this, Qinit, -1.0, false);
949 rcond = fact.rcond();
950 double rcond2;
951 SparseComplexMatrix InvL = fact.L().transpose().
952 tinverse(tmp_typ, info, rcond2, true, false);
953 SparseComplexMatrix InvU = fact.U().
954 tinverse(tmp_typ, info, rcond2, true, false).transpose();
955 ret = fact.Pc().transpose() * InvU * InvL * fact.Pr();
956 }
957 }
958 else
959 (*current_liboctave_error_handler) ("inverse requires square matrix");
960
961 return ret;
606 } 962 }
607 963
608 ComplexDET 964 ComplexDET
609 SparseComplexMatrix::determinant (void) const 965 SparseComplexMatrix::determinant (void) const
610 { 966 {
4644 volatile int typ = mattype.type (); 5000 volatile int typ = mattype.type ();
4645 mattype.info (); 5001 mattype.info ();
4646 5002
4647 if (typ == SparseType::Hermitian) 5003 if (typ == SparseType::Hermitian)
4648 { 5004 {
4649 // XXX FIXME XXX Write the cholesky solver and only fall 5005 #ifdef HAVE_CHOLMOD
4650 // through if cholesky factorization fails 5006 cholmod_common Common;
4651 5007 cholmod_common *cm = &Common;
5008
5009 // Setup initial parameters
5010 CHOLMOD_NAME(start) (cm);
5011 cm->prefer_zomplex = FALSE;
5012
5013 double spu = Voctave_sparse_controls.get_key ("spumoni");
5014 if (spu == 0.)
5015 {
5016 cm->print = -1;
5017 cm->print_function = NULL;
5018 }
5019 else
5020 {
5021 cm->print = (int)spu + 2;
5022 cm->print_function =&SparseCholPrint;
5023 }
5024
5025 cm->error_handler = &SparseCholError;
5026 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5027 cm->hypotenuse = CHOLMOD_NAME(hypot);
5028
5029 #ifdef HAVE_METIS
5030 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5031 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5032 // which mxMalloc's a huge block of memory (and then immediately
5033 // mxFree's it) before calling METIS
5034 cm->metis_memory = 2.0;
5035
5036 #if defined(METIS_VERSION)
5037 #if (METIS_VERSION >= METIS_VER(4,0,2))
5038 // METIS 4.0.2 uses function pointers for malloc and free
5039 METIS_malloc = cm->malloc_memory;
5040 METIS_free = cm->free_memory;
5041 // Turn off METIS memory guard. It is not needed, because mxMalloc
5042 // will safely terminate the mexFunction and free any workspace
5043 // without killing all of octave.
5044 cm->metis_memory = 0.0;
5045 #endif
5046 #endif
5047 #endif
5048
5049 cm->final_ll = TRUE;
5050
5051 cholmod_sparse Astore;
5052 cholmod_sparse *A = &Astore;
5053 double dummy;
5054 A->nrow = nr;
5055 A->ncol = nc;
5056
5057 A->p = cidx();
5058 A->i = ridx();
5059 A->nzmax = nonzero();
5060 A->packed = TRUE;
5061 A->sorted = TRUE;
5062 A->nz = NULL;
5063 #ifdef IDX_TYPE_LONG
5064 A->itype = CHOLMOD_LONG;
5065 #else
5066 A->itype = CHOLMOD_INT;
5067 #endif
5068 A->dtype = CHOLMOD_DOUBLE;
5069 A->stype = 1;
5070 A->xtype = CHOLMOD_COMPLEX;
5071
5072 if (nr < 1)
5073 A->x = &dummy;
5074 else
5075 A->x = data();
5076
5077 cholmod_dense Bstore;
5078 cholmod_dense *B = &Bstore;
5079 B->nrow = b.rows();
5080 B->ncol = b.cols();
5081 B->d = B->nrow;
5082 B->nzmax = B->nrow * B->ncol;
5083 B->dtype = CHOLMOD_DOUBLE;
5084 B->xtype = CHOLMOD_REAL;
5085 if (nc < 1 || b.cols() < 1)
5086 B->x = &dummy;
5087 else
5088 // We won't alter it, honest :-)
5089 B->x = const_cast<double *>(b.fortran_vec());
5090
5091 cholmod_factor *L;
5092 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5093 L = CHOLMOD_NAME(analyze) (A, cm);
5094 CHOLMOD_NAME(factorize) (A, L, cm);
5095 rcond = CHOLMOD_NAME(rcond)(L, cm);
5096 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5097
5098 if (rcond == 0.0)
5099 {
5100 // Either its indefinite or singular. Try UMFPACK
5101 mattype.mark_as_unsymmetric ();
5102 typ = SparseType::Full;
5103 }
5104 else
5105 {
5106 volatile double rcond_plus_one = rcond + 1.0;
5107
5108 if (rcond_plus_one == 1.0 || xisnan (rcond))
5109 {
5110 err = -2;
5111
5112 if (sing_handler)
5113 sing_handler (rcond);
5114 else
5115 (*current_liboctave_error_handler)
5116 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5117 rcond);
5118
5119 #ifdef HAVE_LSSOLVE
5120 return retval;
5121 #endif
5122 }
5123
5124 cholmod_dense *X;
5125 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5126 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5127 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5128
5129 retval.resize (b.rows (), b.cols());
5130 for (octave_idx_type j = 0; j < b.cols(); j++)
5131 {
5132 octave_idx_type jr = j * b.rows();
5133 for (octave_idx_type i = 0; i < b.rows(); i++)
5134 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i];
5135 }
5136
5137 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5138 CHOLMOD_NAME(free_dense) (&X, cm);
5139 CHOLMOD_NAME(free_factor) (&L, cm);
5140 CHOLMOD_NAME(finish) (cm);
5141 CHOLMOD_NAME(print_common) (" ", cm);
5142 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5143 }
5144 #else
4652 (*current_liboctave_warning_handler) 5145 (*current_liboctave_warning_handler)
4653 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5146 ("CHOLMOD not installed");
4654 5147
4655 mattype.mark_as_unsymmetric (); 5148 mattype.mark_as_unsymmetric ();
4656 typ = SparseType::Full; 5149 typ = SparseType::Full;
5150 #endif
4657 } 5151 }
4658 5152
4659 if (typ == SparseType::Full) 5153 if (typ == SparseType::Full)
4660 { 5154 {
4661 #ifdef HAVE_UMFPACK 5155 #ifdef HAVE_UMFPACK
4770 (*current_liboctave_error_handler) 5264 (*current_liboctave_error_handler)
4771 ("matrix dimension mismatch solution of linear equations"); 5265 ("matrix dimension mismatch solution of linear equations");
4772 else 5266 else
4773 { 5267 {
4774 // Print spparms("spumoni") info if requested 5268 // Print spparms("spumoni") info if requested
4775 int typ = mattype.type (); 5269 volatile int typ = mattype.type ();
4776 mattype.info (); 5270 mattype.info ();
4777 5271
4778 if (typ == SparseType::Hermitian) 5272 if (typ == SparseType::Hermitian)
4779 { 5273 {
4780 // XXX FIXME XXX Write the cholesky solver and only fall 5274 #ifdef HAVE_CHOLMOD
4781 // through if cholesky factorization fails 5275 cholmod_common Common;
4782 5276 cholmod_common *cm = &Common;
5277
5278 // Setup initial parameters
5279 CHOLMOD_NAME(start) (cm);
5280 cm->prefer_zomplex = FALSE;
5281
5282 double spu = Voctave_sparse_controls.get_key ("spumoni");
5283 if (spu == 0.)
5284 {
5285 cm->print = -1;
5286 cm->print_function = NULL;
5287 }
5288 else
5289 {
5290 cm->print = (int)spu + 2;
5291 cm->print_function =&SparseCholPrint;
5292 }
5293
5294 cm->error_handler = &SparseCholError;
5295 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5296 cm->hypotenuse = CHOLMOD_NAME(hypot);
5297
5298 #ifdef HAVE_METIS
5299 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5300 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5301 // which mxMalloc's a huge block of memory (and then immediately
5302 // mxFree's it) before calling METIS
5303 cm->metis_memory = 2.0;
5304
5305 #if defined(METIS_VERSION)
5306 #if (METIS_VERSION >= METIS_VER(4,0,2))
5307 // METIS 4.0.2 uses function pointers for malloc and free
5308 METIS_malloc = cm->malloc_memory;
5309 METIS_free = cm->free_memory;
5310 // Turn off METIS memory guard. It is not needed, because mxMalloc
5311 // will safely terminate the mexFunction and free any workspace
5312 // without killing all of octave.
5313 cm->metis_memory = 0.0;
5314 #endif
5315 #endif
5316 #endif
5317
5318 cm->final_ll = TRUE;
5319
5320 cholmod_sparse Astore;
5321 cholmod_sparse *A = &Astore;
5322 double dummy;
5323 A->nrow = nr;
5324 A->ncol = nc;
5325
5326 A->p = cidx();
5327 A->i = ridx();
5328 A->nzmax = nonzero();
5329 A->packed = TRUE;
5330 A->sorted = TRUE;
5331 A->nz = NULL;
5332 #ifdef IDX_TYPE_LONG
5333 A->itype = CHOLMOD_LONG;
5334 #else
5335 A->itype = CHOLMOD_INT;
5336 #endif
5337 A->dtype = CHOLMOD_DOUBLE;
5338 A->stype = 1;
5339 A->xtype = CHOLMOD_COMPLEX;
5340
5341 if (nr < 1)
5342 A->x = &dummy;
5343 else
5344 A->x = data();
5345
5346 cholmod_sparse Bstore;
5347 cholmod_sparse *B = &Bstore;
5348 B->nrow = b.rows();
5349 B->ncol = b.cols();
5350 B->p = b.cidx();
5351 B->i = b.ridx();
5352 B->nzmax = b.nonzero();
5353 B->packed = TRUE;
5354 B->sorted = TRUE;
5355 B->nz = NULL;
5356 #ifdef IDX_TYPE_LONG
5357 B->itype = CHOLMOD_LONG;
5358 #else
5359 B->itype = CHOLMOD_INT;
5360 #endif
5361 B->dtype = CHOLMOD_DOUBLE;
5362 B->stype = 0;
5363 B->xtype = CHOLMOD_REAL;
5364
5365 if (b.rows() < 1 || b.cols() < 1)
5366 B->x = &dummy;
5367 else
5368 B->x = b.data();
5369
5370 cholmod_factor *L;
5371 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5372 L = CHOLMOD_NAME(analyze) (A, cm);
5373 CHOLMOD_NAME(factorize) (A, L, cm);
5374 rcond = CHOLMOD_NAME(rcond)(L, cm);
5375 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5376
5377 if (rcond == 0.0)
5378 {
5379 // Either its indefinite or singular. Try UMFPACK
5380 mattype.mark_as_unsymmetric ();
5381 typ = SparseType::Full;
5382 }
5383 else
5384 {
5385 volatile double rcond_plus_one = rcond + 1.0;
5386
5387 if (rcond_plus_one == 1.0 || xisnan (rcond))
5388 {
5389 err = -2;
5390
5391 if (sing_handler)
5392 sing_handler (rcond);
5393 else
5394 (*current_liboctave_error_handler)
5395 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5396 rcond);
5397
5398 #ifdef HAVE_LSSOLVE
5399 return retval;
5400 #endif
5401 }
5402
5403 cholmod_sparse *X;
5404 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5405 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
5406 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5407
5408 retval = SparseComplexMatrix
5409 (static_cast<octave_idx_type>(X->nrow),
5410 static_cast<octave_idx_type>(X->ncol),
5411 static_cast<octave_idx_type>(X->nzmax));
5412 for (octave_idx_type j = 0;
5413 j <= static_cast<octave_idx_type>(X->ncol); j++)
5414 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
5415 for (octave_idx_type j = 0;
5416 j < static_cast<octave_idx_type>(X->nzmax); j++)
5417 {
5418 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
5419 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
5420 }
5421
5422 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5423 CHOLMOD_NAME(free_sparse) (&X, cm);
5424 CHOLMOD_NAME(free_factor) (&L, cm);
5425 CHOLMOD_NAME(finish) (cm);
5426 CHOLMOD_NAME(print_common) (" ", cm);
5427 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5428 }
5429 #else
4783 (*current_liboctave_warning_handler) 5430 (*current_liboctave_warning_handler)
4784 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5431 ("CHOLMOD not installed");
4785 5432
4786 mattype.mark_as_unsymmetric (); 5433 mattype.mark_as_unsymmetric ();
4787 typ = SparseType::Full; 5434 typ = SparseType::Full;
5435 #endif
4788 } 5436 }
4789 5437
4790 if (typ == SparseType::Full) 5438 if (typ == SparseType::Full)
4791 { 5439 {
4792 #ifdef HAVE_UMFPACK 5440 #ifdef HAVE_UMFPACK
4930 (*current_liboctave_error_handler) 5578 (*current_liboctave_error_handler)
4931 ("matrix dimension mismatch solution of linear equations"); 5579 ("matrix dimension mismatch solution of linear equations");
4932 else 5580 else
4933 { 5581 {
4934 // Print spparms("spumoni") info if requested 5582 // Print spparms("spumoni") info if requested
4935 int typ = mattype.type (); 5583 volatile int typ = mattype.type ();
4936 mattype.info (); 5584 mattype.info ();
4937 5585
4938 if (typ == SparseType::Hermitian) 5586 if (typ == SparseType::Hermitian)
4939 { 5587 {
4940 // XXX FIXME XXX Write the cholesky solver and only fall 5588 #ifdef HAVE_CHOLMOD
4941 // through if cholesky factorization fails 5589 cholmod_common Common;
4942 5590 cholmod_common *cm = &Common;
5591
5592 // Setup initial parameters
5593 CHOLMOD_NAME(start) (cm);
5594 cm->prefer_zomplex = FALSE;
5595
5596 double spu = Voctave_sparse_controls.get_key ("spumoni");
5597 if (spu == 0.)
5598 {
5599 cm->print = -1;
5600 cm->print_function = NULL;
5601 }
5602 else
5603 {
5604 cm->print = (int)spu + 2;
5605 cm->print_function =&SparseCholPrint;
5606 }
5607
5608 cm->error_handler = &SparseCholError;
5609 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5610 cm->hypotenuse = CHOLMOD_NAME(hypot);
5611
5612 #ifdef HAVE_METIS
5613 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5614 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5615 // which mxMalloc's a huge block of memory (and then immediately
5616 // mxFree's it) before calling METIS
5617 cm->metis_memory = 2.0;
5618
5619 #if defined(METIS_VERSION)
5620 #if (METIS_VERSION >= METIS_VER(4,0,2))
5621 // METIS 4.0.2 uses function pointers for malloc and free
5622 METIS_malloc = cm->malloc_memory;
5623 METIS_free = cm->free_memory;
5624 // Turn off METIS memory guard. It is not needed, because mxMalloc
5625 // will safely terminate the mexFunction and free any workspace
5626 // without killing all of octave.
5627 cm->metis_memory = 0.0;
5628 #endif
5629 #endif
5630 #endif
5631
5632 cm->final_ll = TRUE;
5633
5634 cholmod_sparse Astore;
5635 cholmod_sparse *A = &Astore;
5636 double dummy;
5637 A->nrow = nr;
5638 A->ncol = nc;
5639
5640 A->p = cidx();
5641 A->i = ridx();
5642 A->nzmax = nonzero();
5643 A->packed = TRUE;
5644 A->sorted = TRUE;
5645 A->nz = NULL;
5646 #ifdef IDX_TYPE_LONG
5647 A->itype = CHOLMOD_LONG;
5648 #else
5649 A->itype = CHOLMOD_INT;
5650 #endif
5651 A->dtype = CHOLMOD_DOUBLE;
5652 A->stype = 1;
5653 A->xtype = CHOLMOD_COMPLEX;
5654
5655 if (nr < 1)
5656 A->x = &dummy;
5657 else
5658 A->x = data();
5659
5660 cholmod_dense Bstore;
5661 cholmod_dense *B = &Bstore;
5662 B->nrow = b.rows();
5663 B->ncol = b.cols();
5664 B->d = B->nrow;
5665 B->nzmax = B->nrow * B->ncol;
5666 B->dtype = CHOLMOD_DOUBLE;
5667 B->xtype = CHOLMOD_COMPLEX;
5668 if (nc < 1 || b.cols() < 1)
5669 B->x = &dummy;
5670 else
5671 // We won't alter it, honest :-)
5672 B->x = const_cast<Complex *>(b.fortran_vec());
5673
5674 cholmod_factor *L;
5675 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5676 L = CHOLMOD_NAME(analyze) (A, cm);
5677 CHOLMOD_NAME(factorize) (A, L, cm);
5678 rcond = CHOLMOD_NAME(rcond)(L, cm);
5679 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5680
5681 if (rcond == 0.0)
5682 {
5683 // Either its indefinite or singular. Try UMFPACK
5684 mattype.mark_as_unsymmetric ();
5685 typ = SparseType::Full;
5686 }
5687 else
5688 {
5689 volatile double rcond_plus_one = rcond + 1.0;
5690
5691 if (rcond_plus_one == 1.0 || xisnan (rcond))
5692 {
5693 err = -2;
5694
5695 if (sing_handler)
5696 sing_handler (rcond);
5697 else
5698 (*current_liboctave_error_handler)
5699 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5700 rcond);
5701
5702 #ifdef HAVE_LSSOLVE
5703 return retval;
5704 #endif
5705 }
5706
5707 cholmod_dense *X;
5708 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5709 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5710 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5711
5712 retval.resize (b.rows (), b.cols());
5713 for (octave_idx_type j = 0; j < b.cols(); j++)
5714 {
5715 octave_idx_type jr = j * b.rows();
5716 for (octave_idx_type i = 0; i < b.rows(); i++)
5717 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i];
5718 }
5719
5720 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5721 CHOLMOD_NAME(free_dense) (&X, cm);
5722 CHOLMOD_NAME(free_factor) (&L, cm);
5723 CHOLMOD_NAME(finish) (cm);
5724 CHOLMOD_NAME(print_common) (" ", cm);
5725 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5726 }
5727 #else
4943 (*current_liboctave_warning_handler) 5728 (*current_liboctave_warning_handler)
4944 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5729 ("CHOLMOD not installed");
4945 5730
4946 mattype.mark_as_unsymmetric (); 5731 mattype.mark_as_unsymmetric ();
4947 typ = SparseType::Full; 5732 typ = SparseType::Full;
5733 #endif
4948 } 5734 }
4949 5735
4950 if (typ == SparseType::Full) 5736 if (typ == SparseType::Full)
4951 { 5737 {
4952 #ifdef HAVE_UMFPACK 5738 #ifdef HAVE_UMFPACK
5039 (*current_liboctave_error_handler) 5825 (*current_liboctave_error_handler)
5040 ("matrix dimension mismatch solution of linear equations"); 5826 ("matrix dimension mismatch solution of linear equations");
5041 else 5827 else
5042 { 5828 {
5043 // Print spparms("spumoni") info if requested 5829 // Print spparms("spumoni") info if requested
5044 int typ = mattype.type (); 5830 volatile int typ = mattype.type ();
5045 mattype.info (); 5831 mattype.info ();
5046 5832
5047 if (typ == SparseType::Hermitian) 5833 if (typ == SparseType::Hermitian)
5048 { 5834 {
5049 // XXX FIXME XXX Write the cholesky solver and only fall 5835 #ifdef HAVE_CHOLMOD
5050 // through if cholesky factorization fails 5836 cholmod_common Common;
5051 5837 cholmod_common *cm = &Common;
5838
5839 // Setup initial parameters
5840 CHOLMOD_NAME(start) (cm);
5841 cm->prefer_zomplex = FALSE;
5842
5843 double spu = Voctave_sparse_controls.get_key ("spumoni");
5844 if (spu == 0.)
5845 {
5846 cm->print = -1;
5847 cm->print_function = NULL;
5848 }
5849 else
5850 {
5851 cm->print = (int)spu + 2;
5852 cm->print_function =&SparseCholPrint;
5853 }
5854
5855 cm->error_handler = &SparseCholError;
5856 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5857 cm->hypotenuse = CHOLMOD_NAME(hypot);
5858
5859 #ifdef HAVE_METIS
5860 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5861 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5862 // which mxMalloc's a huge block of memory (and then immediately
5863 // mxFree's it) before calling METIS
5864 cm->metis_memory = 2.0;
5865
5866 #if defined(METIS_VERSION)
5867 #if (METIS_VERSION >= METIS_VER(4,0,2))
5868 // METIS 4.0.2 uses function pointers for malloc and free
5869 METIS_malloc = cm->malloc_memory;
5870 METIS_free = cm->free_memory;
5871 // Turn off METIS memory guard. It is not needed, because mxMalloc
5872 // will safely terminate the mexFunction and free any workspace
5873 // without killing all of octave.
5874 cm->metis_memory = 0.0;
5875 #endif
5876 #endif
5877 #endif
5878
5879 cm->final_ll = TRUE;
5880
5881 cholmod_sparse Astore;
5882 cholmod_sparse *A = &Astore;
5883 double dummy;
5884 A->nrow = nr;
5885 A->ncol = nc;
5886
5887 A->p = cidx();
5888 A->i = ridx();
5889 A->nzmax = nonzero();
5890 A->packed = TRUE;
5891 A->sorted = TRUE;
5892 A->nz = NULL;
5893 #ifdef IDX_TYPE_LONG
5894 A->itype = CHOLMOD_LONG;
5895 #else
5896 A->itype = CHOLMOD_INT;
5897 #endif
5898 A->dtype = CHOLMOD_DOUBLE;
5899 A->stype = 1;
5900 A->xtype = CHOLMOD_COMPLEX;
5901
5902 if (nr < 1)
5903 A->x = &dummy;
5904 else
5905 A->x = data();
5906
5907 cholmod_sparse Bstore;
5908 cholmod_sparse *B = &Bstore;
5909 B->nrow = b.rows();
5910 B->ncol = b.cols();
5911 B->p = b.cidx();
5912 B->i = b.ridx();
5913 B->nzmax = b.nonzero();
5914 B->packed = TRUE;
5915 B->sorted = TRUE;
5916 B->nz = NULL;
5917 #ifdef IDX_TYPE_LONG
5918 B->itype = CHOLMOD_LONG;
5919 #else
5920 B->itype = CHOLMOD_INT;
5921 #endif
5922 B->dtype = CHOLMOD_DOUBLE;
5923 B->stype = 0;
5924 B->xtype = CHOLMOD_COMPLEX;
5925
5926 if (b.rows() < 1 || b.cols() < 1)
5927 B->x = &dummy;
5928 else
5929 B->x = b.data();
5930
5931 cholmod_factor *L;
5932 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5933 L = CHOLMOD_NAME(analyze) (A, cm);
5934 CHOLMOD_NAME(factorize) (A, L, cm);
5935 rcond = CHOLMOD_NAME(rcond)(L, cm);
5936 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5937
5938 if (rcond == 0.0)
5939 {
5940 // Either its indefinite or singular. Try UMFPACK
5941 mattype.mark_as_unsymmetric ();
5942 typ = SparseType::Full;
5943 }
5944 else
5945 {
5946 volatile double rcond_plus_one = rcond + 1.0;
5947
5948 if (rcond_plus_one == 1.0 || xisnan (rcond))
5949 {
5950 err = -2;
5951
5952 if (sing_handler)
5953 sing_handler (rcond);
5954 else
5955 (*current_liboctave_error_handler)
5956 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5957 rcond);
5958
5959 #ifdef HAVE_LSSOLVE
5960 return retval;
5961 #endif
5962 }
5963
5964 cholmod_sparse *X;
5965 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5966 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
5967 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5968
5969 retval = SparseComplexMatrix
5970 (static_cast<octave_idx_type>(X->nrow),
5971 static_cast<octave_idx_type>(X->ncol),
5972 static_cast<octave_idx_type>(X->nzmax));
5973 for (octave_idx_type j = 0;
5974 j <= static_cast<octave_idx_type>(X->ncol); j++)
5975 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
5976 for (octave_idx_type j = 0;
5977 j < static_cast<octave_idx_type>(X->nzmax); j++)
5978 {
5979 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
5980 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
5981 }
5982
5983 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5984 CHOLMOD_NAME(free_sparse) (&X, cm);
5985 CHOLMOD_NAME(free_factor) (&L, cm);
5986 CHOLMOD_NAME(finish) (cm);
5987 CHOLMOD_NAME(print_common) (" ", cm);
5988 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5989 }
5990 #else
5052 (*current_liboctave_warning_handler) 5991 (*current_liboctave_warning_handler)
5053 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5992 ("CHOLMOD not installed");
5054 5993
5055 mattype.mark_as_unsymmetric (); 5994 mattype.mark_as_unsymmetric ();
5056 typ = SparseType::Full; 5995 typ = SparseType::Full;
5996 #endif
5057 } 5997 }
5058 5998
5059 if (typ == SparseType::Full) 5999 if (typ == SparseType::Full)
5060 { 6000 {
5061 #ifdef HAVE_UMFPACK 6001 #ifdef HAVE_UMFPACK