Mercurial > octave-nkf

comparison liboctave/dSparse.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
40 #include "dSparse.h" 40 #include "dSparse.h"
41 #include "oct-spparms.h" 41 #include "oct-spparms.h"
42 #include "SparsedbleLU.h" 42 #include "SparsedbleLU.h"
43 #include "SparseType.h" 43 #include "SparseType.h"
44 #include "oct-sparse.h" 44 #include "oct-sparse.h"
45 #include "sparse-util.h"
46 #include "SparsedbleCHOL.h"
45 47
46 // Fortran functions we call. 48 // Fortran functions we call.
47 extern "C" 49 extern "C"
48 { 50 {
49 F77_RET_T 51 F77_RET_T
670 SparseMatrix 672 SparseMatrix
671 SparseMatrix::inverse (void) const 673 SparseMatrix::inverse (void) const
672 { 674 {
673 octave_idx_type info; 675 octave_idx_type info;
674 double rcond; 676 double rcond;
675 return inverse (info, rcond, 0, 0); 677 SparseType mattype (*this);
678 return inverse (mattype, info, rcond, 0, 0);
676 } 679 }
677 680
678 SparseMatrix 681 SparseMatrix
679 SparseMatrix::inverse (octave_idx_type& info) const 682 SparseMatrix::inverse (SparseType& mattype) const
680 { 683 {
684 octave_idx_type info;
681 double rcond; 685 double rcond;
682 return inverse (info, rcond, 0, 0); 686 return inverse (mattype, info, rcond, 0, 0);
683 } 687 }
684 688
685 SparseMatrix 689 SparseMatrix
686 SparseMatrix::inverse (octave_idx_type& info, double& rcond, int force, int calc_cond) const 690 SparseMatrix::inverse (SparseType& mattype, octave_idx_type& info) const
687 { 691 {
688 info = -1; 692 double rcond;
689 (*current_liboctave_error_handler) 693 return inverse (mattype, info, rcond, 0, 0);
690 ("SparseMatrix::inverse not implemented yet"); 694 }
691 return SparseMatrix (); 695
696 SparseMatrix
697 SparseMatrix::dinverse (SparseType &mattyp, octave_idx_type& info,
698 double& rcond, const bool force,
699 const bool calccond) const
700 {
701 SparseMatrix retval;
702
703 octave_idx_type nr = rows ();
704 octave_idx_type nc = cols ();
705 info = 0;
706
707 if (nr == 0 || nc == 0 || nr != nc)
708 (*current_liboctave_error_handler) ("inverse requires square matrix");
709 else
710 {
711 // Print spparms("spumoni") info if requested
712 int typ = mattyp.type ();
713 mattyp.info ();
714
715 if (typ == SparseType::Diagonal ||
716 typ == SparseType::Permuted_Diagonal)
717 {
718 if (typ == SparseType::Permuted_Diagonal)
719 retval = transpose();
720 else
721 retval = *this;
722
723 // Force make_unique to be called
724 double *v = retval.data();
725
726 if (calccond)
727 {
728 double dmax = 0., dmin = octave_Inf;
729 for (octave_idx_type i = 0; i < nr; i++)
730 {
731 double tmp = fabs(v[i]);
732 if (tmp > dmax)
733 dmax = tmp;
734 if (tmp < dmin)
735 dmin = tmp;
736 }
737 rcond = dmin / dmax;
738 }
739
740 for (octave_idx_type i = 0; i < nr; i++)
741 v[i] = 1.0 / v[i];
742 }
743 else
744 (*current_liboctave_error_handler) ("incorrect matrix type");
745 }
746
747 return retval;
748 }
749
750 SparseMatrix
751 SparseMatrix::tinverse (SparseType &mattyp, octave_idx_type& info,
752 double& rcond, const bool force,
753 const bool calccond) const
754 {
755 SparseMatrix retval;
756
757 octave_idx_type nr = rows ();
758 octave_idx_type nc = cols ();
759 info = 0;
760
761 if (nr == 0 || nc == 0 || nr != nc)
762 (*current_liboctave_error_handler) ("inverse requires square matrix");
763 else
764 {
765 // Print spparms("spumoni") info if requested
766 int typ = mattyp.type ();
767 mattyp.info ();
768
769 if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper ||
770 typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
771 {
772 double anorm = 0.;
773 double ainvnorm = 0.;
774
775 if (calccond)
776 {
777 // Calculate the 1-norm of matrix for rcond calculation
778 for (octave_idx_type j = 0; j < nr; j++)
779 {
780 double atmp = 0.;
781 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
782 atmp += fabs(data(i));
783 if (atmp > anorm)
784 anorm = atmp;
785 }
786 }
787
788 if (typ == SparseType::Upper || typ == SparseType::Lower)
789 {
790 octave_idx_type nz = nnz();
791 octave_idx_type cx = 0;
792 octave_idx_type nz2 = nz;
793 retval = SparseMatrix (nr, nc, nz2);
794
795 for (octave_idx_type i = 0; i < nr; i++)
796 {
797 OCTAVE_QUIT;
798 // place the 1 in the identity position
799 octave_idx_type cx_colstart = cx;
800
801 if (cx == nz2)
802 {
803 nz2 *= 2;
804 retval.change_capacity (nz2);
805 }
806
807 retval.xcidx(i) = cx;
808 retval.xridx(cx) = i;
809 retval.xdata(cx) = 1.0;
810 cx++;
811
812 // iterate accross columns of input matrix
813 for (octave_idx_type j = i+1; j < nr; j++)
814 {
815 double v = 0.;
816 // iterate to calculate sum
817 octave_idx_type colXp = retval.xcidx(i);
818 octave_idx_type colUp = cidx(j);
819 octave_idx_type rpX, rpU;
820 do
821 {
822 OCTAVE_QUIT;
823 rpX = retval.xridx(colXp);
824 rpU = ridx(colUp);
825
826 if (rpX < rpU)
827 colXp++;
828 else if (rpX > rpU)
829 colUp++;
830 else
831 {
832 v -= retval.xdata(colXp) * data(colUp);
833 colXp++;
834 colUp++;
835 }
836 } while ((rpX<j) && (rpU<j) &&
837 (colXp<cx) && (colUp<nz));
838
839 // get A(m,m)
840 colUp = cidx(j+1) - 1;
841 double pivot = data(colUp);
842 if (pivot == 0.)
843 (*current_liboctave_error_handler)
844 ("division by zero");
845
846 if (v != 0.)
847 {
848 if (cx == nz2)
849 {
850 nz2 *= 2;
851 retval.change_capacity (nz2);
852 }
853
854 retval.xridx(cx) = j;
855 retval.xdata(cx) = v / pivot;
856 cx++;
857 }
858 }
859
860 // get A(m,m)
861 octave_idx_type colUp = cidx(i+1) - 1;
862 double pivot = data(colUp);
863 if (pivot == 0.)
864 (*current_liboctave_error_handler) ("division by zero");
865
866 if (pivot != 1.0)
867 for (octave_idx_type j = cx_colstart; j < cx; j++)
868 retval.xdata(j) /= pivot;
869 }
870 retval.xcidx(nr) = cx;
871 retval.maybe_compress ();
872 }
873 else
874 {
875 octave_idx_type nz = nnz();
876 octave_idx_type cx = 0;
877 octave_idx_type nz2 = nz;
878 retval = SparseMatrix (nr, nc, nz2);
879
880 OCTAVE_LOCAL_BUFFER (double, work, nr);
881 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
882
883 octave_idx_type *perm = mattyp.triangular_perm();
884 if (typ == SparseType::Permuted_Upper)
885 {
886 for (octave_idx_type i = 0; i < nr; i++)
887 rperm[perm[i]] = i;
888 }
889 else
890 {
891 for (octave_idx_type i = 0; i < nr; i++)
892 rperm[i] = perm[i];
893 for (octave_idx_type i = 0; i < nr; i++)
894 perm[rperm[i]] = i;
895 }
896
897 for (octave_idx_type i = 0; i < nr; i++)
898 {
899 OCTAVE_QUIT;
900 octave_idx_type iidx = rperm[i];
901
902 for (octave_idx_type j = 0; j < nr; j++)
903 work[j] = 0.;
904
905 // place the 1 in the identity position
906 work[iidx] = 1.0;
907
908 // iterate accross columns of input matrix
909 for (octave_idx_type j = iidx+1; j < nr; j++)
910 {
911 double v = 0.;
912 octave_idx_type jidx = perm[j];
913 // iterate to calculate sum
914 for (octave_idx_type k = cidx(jidx);
915 k < cidx(jidx+1); k++)
916 {
917 OCTAVE_QUIT;
918 v -= work[ridx(k)] * data(k);
919 }
920
921 // get A(m,m)
922 double pivot = data(cidx(jidx+1) - 1);
923 if (pivot == 0.)
924 (*current_liboctave_error_handler)
925 ("division by zero");
926
927 work[j] = v / pivot;
928 }
929
930 // get A(m,m)
931 octave_idx_type colUp = cidx(perm[iidx]+1) - 1;
932 double pivot = data(colUp);
933 if (pivot == 0.)
934 (*current_liboctave_error_handler)
935 ("division by zero");
936
937 octave_idx_type new_cx = cx;
938 for (octave_idx_type j = iidx; j < nr; j++)
939 if (work[j] != 0.0)
940 {
941 new_cx++;
942 if (pivot != 1.0)
943 work[j] /= pivot;
944 }
945
946 if (cx < new_cx)
947 {
948 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
949 retval.change_capacity (nz2);
950 }
951
952 retval.xcidx(i) = cx;
953 for (octave_idx_type j = iidx; j < nr; j++)
954 if (work[j] != 0.)
955 {
956 retval.xridx(cx) = j;
957 retval.xdata(cx++) = work[j];
958 }
959 }
960
961 retval.xcidx(nr) = cx;
962 retval.maybe_compress ();
963 }
964
965 if (calccond)
966 {
967 // Calculate the 1-norm of inverse matrix for rcond calculation
968 for (octave_idx_type j = 0; j < nr; j++)
969 {
970 double atmp = 0.;
971 for (octave_idx_type i = retval.cidx(j);
972 i < retval.cidx(j+1); i++)
973 atmp += fabs(retval.data(i));
974 if (atmp > ainvnorm)
975 ainvnorm = atmp;
976 }
977
978 rcond = 1. / ainvnorm / anorm;
979 }
980 }
981 else
982 (*current_liboctave_error_handler) ("incorrect matrix type");
983 }
984
985 return retval;
986 }
987
988 SparseMatrix
989 SparseMatrix::inverse (SparseType &mattype, octave_idx_type& info,
990 double& rcond, int force, int calc_cond) const
991 {
992 int typ = mattype.type (false);
993 SparseMatrix ret;
994
995 if (typ == SparseType::Unknown)
996 typ = mattype.type (*this);
997
998 if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal)
999 ret = dinverse (mattype, info, rcond, true, calc_cond);
1000 else if (typ == SparseType::Upper || typ == SparseType::Permuted_Upper)
1001 ret = tinverse (mattype, info, rcond, true, calc_cond).transpose();
1002 else if (typ == SparseType::Lower || typ == SparseType::Permuted_Lower)
1003 ret = transpose().tinverse (mattype, info, rcond, true, calc_cond);
1004 else if (typ != SparseType::Rectangular)
1005 {
1006 if (mattype.is_hermitian())
1007 {
1008 SparseType tmp_typ (SparseType::Upper);
1009 SparseCHOL fact (*this, info, false);
1010 rcond = fact.rcond();
1011 if (info == 0)
1012 {
1013 double rcond2;
1014 SparseMatrix Q = fact.Q();
1015 SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ,
1016 info, rcond2, true, false);
1017 ret = Q * InvL.transpose() * InvL * Q.transpose();
1018 }
1019 else
1020 {
1021 // Matrix is either singular or not positive definite
1022 mattype.mark_as_unsymmetric ();
1023 typ = SparseType::Full;
1024 }
1025 }
1026
1027 if (!mattype.is_hermitian())
1028 {
1029 octave_idx_type n = rows();
1030 ColumnVector Qinit(n);
1031 for (octave_idx_type i = 0; i < n; i++)
1032 Qinit(i) = i;
1033
1034 SparseType tmp_typ (SparseType::Upper);
1035 SparseLU fact (*this, Qinit, -1.0, false);
1036 rcond = fact.rcond();
1037 double rcond2;
1038 SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ,
1039 info, rcond2, true, false);
1040 SparseMatrix InvU = fact.U().tinverse(tmp_typ, info, rcond2,
1041 true, false).transpose();
1042 ret = fact.Pc().transpose() * InvU * InvL * fact.Pr();
1043 }
1044 }
1045 else
1046 (*current_liboctave_error_handler) ("inverse requires square matrix");
1047
1048 return ret;
692 } 1049 }
693 1050
694 DET 1051 DET
695 SparseMatrix::determinant (void) const 1052 SparseMatrix::determinant (void) const
696 { 1053 {
4857 (*current_liboctave_error_handler) 5214 (*current_liboctave_error_handler)
4858 ("matrix dimension mismatch solution of linear equations"); 5215 ("matrix dimension mismatch solution of linear equations");
4859 else 5216 else
4860 { 5217 {
4861 // Print spparms("spumoni") info if requested 5218 // Print spparms("spumoni") info if requested
4862 int typ = mattype.type (); 5219 volatile int typ = mattype.type ();
4863 mattype.info (); 5220 mattype.info ();
4864 5221
4865 if (typ == SparseType::Hermitian) 5222 if (typ == SparseType::Hermitian)
4866 { 5223 {
4867 // XXX FIXME XXX Write the cholesky solver and only fall 5224 #ifdef HAVE_CHOLMOD
4868 // through if cholesky factorization fails 5225 cholmod_common Common;
4869 5226 cholmod_common *cm = &Common;
5227
5228 // Setup initial parameters
5229 CHOLMOD_NAME(start) (cm);
5230 cm->prefer_zomplex = FALSE;
5231
5232 double spu = Voctave_sparse_controls.get_key ("spumoni");
5233 if (spu == 0.)
5234 {
5235 cm->print = -1;
5236 cm->print_function = NULL;
5237 }
5238 else
5239 {
5240 cm->print = (int)spu + 2;
5241 cm->print_function =&SparseCholPrint;
5242 }
5243
5244 cm->error_handler = &SparseCholError;
5245 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5246 cm->hypotenuse = CHOLMOD_NAME(hypot);
5247
5248 #ifdef HAVE_METIS
5249 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5250 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5251 // which mxMalloc's a huge block of memory (and then immediately
5252 // mxFree's it) before calling METIS
5253 cm->metis_memory = 2.0;
5254
5255 #if defined(METIS_VERSION)
5256 #if (METIS_VERSION >= METIS_VER(4,0,2))
5257 // METIS 4.0.2 uses function pointers for malloc and free
5258 METIS_malloc = cm->malloc_memory;
5259 METIS_free = cm->free_memory;
5260 // Turn off METIS memory guard. It is not needed, because mxMalloc
5261 // will safely terminate the mexFunction and free any workspace
5262 // without killing all of octave.
5263 cm->metis_memory = 0.0;
5264 #endif
5265 #endif
5266 #endif
5267
5268 cm->final_ll = TRUE;
5269
5270 cholmod_sparse Astore;
5271 cholmod_sparse *A = &Astore;
5272 double dummy;
5273 A->nrow = nr;
5274 A->ncol = nc;
5275
5276 A->p = cidx();
5277 A->i = ridx();
5278 A->nzmax = nonzero();
5279 A->packed = TRUE;
5280 A->sorted = TRUE;
5281 A->nz = NULL;
5282 #ifdef IDX_TYPE_LONG
5283 A->itype = CHOLMOD_LONG;
5284 #else
5285 A->itype = CHOLMOD_INT;
5286 #endif
5287 A->dtype = CHOLMOD_DOUBLE;
5288 A->stype = 1;
5289 A->xtype = CHOLMOD_REAL;
5290
5291 if (nr < 1)
5292 A->x = &dummy;
5293 else
5294 A->x = data();
5295
5296 cholmod_dense Bstore;
5297 cholmod_dense *B = &Bstore;
5298 B->nrow = b.rows();
5299 B->ncol = b.cols();
5300 B->d = B->nrow;
5301 B->nzmax = B->nrow * B->ncol;
5302 B->dtype = CHOLMOD_DOUBLE;
5303 B->xtype = CHOLMOD_REAL;
5304 if (nc < 1 || b.cols() < 1)
5305 B->x = &dummy;
5306 else
5307 // We won't alter it, honest :-)
5308 B->x = const_cast<double *>(b.fortran_vec());
5309
5310 cholmod_factor *L;
5311 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5312 L = CHOLMOD_NAME(analyze) (A, cm);
5313 CHOLMOD_NAME(factorize) (A, L, cm);
5314 rcond = CHOLMOD_NAME(rcond)(L, cm);
5315 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5316
5317 if (rcond == 0.0)
5318 {
5319 // Either its indefinite or singular. Try UMFPACK
5320 mattype.mark_as_unsymmetric ();
5321 typ = SparseType::Full;
5322 }
5323 else
5324 {
5325 volatile double rcond_plus_one = rcond + 1.0;
5326
5327 if (rcond_plus_one == 1.0 || xisnan (rcond))
5328 {
5329 err = -2;
5330
5331 if (sing_handler)
5332 sing_handler (rcond);
5333 else
5334 (*current_liboctave_error_handler)
5335 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5336 rcond);
5337
5338 #ifdef HAVE_LSSOLVE
5339 return retval;
5340 #endif
5341 }
5342
5343 cholmod_dense *X;
5344 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5345 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5346 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5347
5348 retval.resize (b.rows (), b.cols());
5349 for (octave_idx_type j = 0; j < b.cols(); j++)
5350 {
5351 octave_idx_type jr = j * b.rows();
5352 for (octave_idx_type i = 0; i < b.rows(); i++)
5353 retval.xelem(i,j) = static_cast<double *>(X->x)[jr + i];
5354 }
5355
5356 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5357 CHOLMOD_NAME(free_dense) (&X, cm);
5358 CHOLMOD_NAME(free_factor) (&L, cm);
5359 CHOLMOD_NAME(finish) (cm);
5360 CHOLMOD_NAME(print_common) (" ", cm);
5361 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5362 }
5363 #else
4870 (*current_liboctave_warning_handler) 5364 (*current_liboctave_warning_handler)
4871 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5365 ("CHOLMOD not installed");
4872 5366
4873 mattype.mark_as_unsymmetric (); 5367 mattype.mark_as_unsymmetric ();
4874 typ = SparseType::Full; 5368 typ = SparseType::Full;
5369 #endif
4875 } 5370 }
4876 5371
4877 if (typ == SparseType::Full) 5372 if (typ == SparseType::Full)
4878 { 5373 {
4879 #ifdef HAVE_UMFPACK 5374 #ifdef HAVE_UMFPACK
4961 (*current_liboctave_error_handler) 5456 (*current_liboctave_error_handler)
4962 ("matrix dimension mismatch solution of linear equations"); 5457 ("matrix dimension mismatch solution of linear equations");
4963 else 5458 else
4964 { 5459 {
4965 // Print spparms("spumoni") info if requested 5460 // Print spparms("spumoni") info if requested
4966 int typ = mattype.type (); 5461 volatile int typ = mattype.type ();
4967 mattype.info (); 5462 mattype.info ();
4968 5463
4969 if (typ == SparseType::Hermitian) 5464 if (typ == SparseType::Hermitian)
4970 { 5465 {
4971 // XXX FIXME XXX Write the cholesky solver and only fall 5466 #ifdef HAVE_CHOLMOD
4972 // through if cholesky factorization fails 5467 cholmod_common Common;
4973 5468 cholmod_common *cm = &Common;
5469
5470 // Setup initial parameters
5471 CHOLMOD_NAME(start) (cm);
5472 cm->prefer_zomplex = FALSE;
5473
5474 double spu = Voctave_sparse_controls.get_key ("spumoni");
5475 if (spu == 0.)
5476 {
5477 cm->print = -1;
5478 cm->print_function = NULL;
5479 }
5480 else
5481 {
5482 cm->print = (int)spu + 2;
5483 cm->print_function =&SparseCholPrint;
5484 }
5485
5486 cm->error_handler = &SparseCholError;
5487 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5488 cm->hypotenuse = CHOLMOD_NAME(hypot);
5489
5490 #ifdef HAVE_METIS
5491 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5492 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5493 // which mxMalloc's a huge block of memory (and then immediately
5494 // mxFree's it) before calling METIS
5495 cm->metis_memory = 2.0;
5496
5497 #if defined(METIS_VERSION)
5498 #if (METIS_VERSION >= METIS_VER(4,0,2))
5499 // METIS 4.0.2 uses function pointers for malloc and free
5500 METIS_malloc = cm->malloc_memory;
5501 METIS_free = cm->free_memory;
5502 // Turn off METIS memory guard. It is not needed, because mxMalloc
5503 // will safely terminate the mexFunction and free any workspace
5504 // without killing all of octave.
5505 cm->metis_memory = 0.0;
5506 #endif
5507 #endif
5508 #endif
5509
5510 cm->final_ll = TRUE;
5511
5512 cholmod_sparse Astore;
5513 cholmod_sparse *A = &Astore;
5514 double dummy;
5515 A->nrow = nr;
5516 A->ncol = nc;
5517
5518 A->p = cidx();
5519 A->i = ridx();
5520 A->nzmax = nonzero();
5521 A->packed = TRUE;
5522 A->sorted = TRUE;
5523 A->nz = NULL;
5524 #ifdef IDX_TYPE_LONG
5525 A->itype = CHOLMOD_LONG;
5526 #else
5527 A->itype = CHOLMOD_INT;
5528 #endif
5529 A->dtype = CHOLMOD_DOUBLE;
5530 A->stype = 1;
5531 A->xtype = CHOLMOD_REAL;
5532
5533 if (nr < 1)
5534 A->x = &dummy;
5535 else
5536 A->x = data();
5537
5538 cholmod_sparse Bstore;
5539 cholmod_sparse *B = &Bstore;
5540 B->nrow = b.rows();
5541 B->ncol = b.cols();
5542 B->p = b.cidx();
5543 B->i = b.ridx();
5544 B->nzmax = b.nonzero();
5545 B->packed = TRUE;
5546 B->sorted = TRUE;
5547 B->nz = NULL;
5548 #ifdef IDX_TYPE_LONG
5549 B->itype = CHOLMOD_LONG;
5550 #else
5551 B->itype = CHOLMOD_INT;
5552 #endif
5553 B->dtype = CHOLMOD_DOUBLE;
5554 B->stype = 0;
5555 B->xtype = CHOLMOD_REAL;
5556
5557 if (b.rows() < 1 || b.cols() < 1)
5558 B->x = &dummy;
5559 else
5560 B->x = b.data();
5561
5562 cholmod_factor *L;
5563 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5564 L = CHOLMOD_NAME(analyze) (A, cm);
5565 CHOLMOD_NAME(factorize) (A, L, cm);
5566 rcond = CHOLMOD_NAME(rcond)(L, cm);
5567 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5568
5569 if (rcond == 0.0)
5570 {
5571 // Either its indefinite or singular. Try UMFPACK
5572 mattype.mark_as_unsymmetric ();
5573 typ = SparseType::Full;
5574 }
5575 else
5576 {
5577 volatile double rcond_plus_one = rcond + 1.0;
5578
5579 if (rcond_plus_one == 1.0 || xisnan (rcond))
5580 {
5581 err = -2;
5582
5583 if (sing_handler)
5584 sing_handler (rcond);
5585 else
5586 (*current_liboctave_error_handler)
5587 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5588 rcond);
5589
5590 #ifdef HAVE_LSSOLVE
5591 return retval;
5592 #endif
5593 }
5594
5595 cholmod_sparse *X;
5596 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5597 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
5598 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5599
5600 retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow),
5601 static_cast<octave_idx_type>(X->ncol),
5602 static_cast<octave_idx_type>(X->nzmax));
5603 for (octave_idx_type j = 0;
5604 j <= static_cast<octave_idx_type>(X->ncol); j++)
5605 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
5606 for (octave_idx_type j = 0;
5607 j < static_cast<octave_idx_type>(X->nzmax); j++)
5608 {
5609 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
5610 retval.xdata(j) = static_cast<double *>(X->x)[j];
5611 }
5612
5613 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5614 CHOLMOD_NAME(free_sparse) (&X, cm);
5615 CHOLMOD_NAME(free_factor) (&L, cm);
5616 CHOLMOD_NAME(finish) (cm);
5617 CHOLMOD_NAME(print_common) (" ", cm);
5618 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5619 }
5620 #else
4974 (*current_liboctave_warning_handler) 5621 (*current_liboctave_warning_handler)
4975 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5622 ("CHOLMOD not installed");
4976 5623
4977 mattype.mark_as_unsymmetric (); 5624 mattype.mark_as_unsymmetric ();
4978 typ = SparseType::Full; 5625 typ = SparseType::Full;
5626 #endif
4979 } 5627 }
4980 5628
4981 if (typ == SparseType::Full) 5629 if (typ == SparseType::Full)
4982 { 5630 {
4983 #ifdef HAVE_UMFPACK 5631 #ifdef HAVE_UMFPACK
5097 (*current_liboctave_error_handler) 5745 (*current_liboctave_error_handler)
5098 ("matrix dimension mismatch solution of linear equations"); 5746 ("matrix dimension mismatch solution of linear equations");
5099 else 5747 else
5100 { 5748 {
5101 // Print spparms("spumoni") info if requested 5749 // Print spparms("spumoni") info if requested
5102 int typ = mattype.type (); 5750 volatile int typ = mattype.type ();
5103 mattype.info (); 5751 mattype.info ();
5104 5752
5105 if (typ == SparseType::Hermitian) 5753 if (typ == SparseType::Hermitian)
5106 { 5754 {
5107 // XXX FIXME XXX Write the cholesky solver and only fall 5755 #ifdef HAVE_CHOLMOD
5108 // through if cholesky factorization fails 5756 cholmod_common Common;
5109 5757 cholmod_common *cm = &Common;
5758
5759 // Setup initial parameters
5760 CHOLMOD_NAME(start) (cm);
5761 cm->prefer_zomplex = FALSE;
5762
5763 double spu = Voctave_sparse_controls.get_key ("spumoni");
5764 if (spu == 0.)
5765 {
5766 cm->print = -1;
5767 cm->print_function = NULL;
5768 }
5769 else
5770 {
5771 cm->print = (int)spu + 2;
5772 cm->print_function =&SparseCholPrint;
5773 }
5774
5775 cm->error_handler = &SparseCholError;
5776 cm->complex_divide = CHOLMOD_NAME(divcomplex);
5777 cm->hypotenuse = CHOLMOD_NAME(hypot);
5778
5779 #ifdef HAVE_METIS
5780 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
5781 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
5782 // which mxMalloc's a huge block of memory (and then immediately
5783 // mxFree's it) before calling METIS
5784 cm->metis_memory = 2.0;
5785
5786 #if defined(METIS_VERSION)
5787 #if (METIS_VERSION >= METIS_VER(4,0,2))
5788 // METIS 4.0.2 uses function pointers for malloc and free
5789 METIS_malloc = cm->malloc_memory;
5790 METIS_free = cm->free_memory;
5791 // Turn off METIS memory guard. It is not needed, because mxMalloc
5792 // will safely terminate the mexFunction and free any workspace
5793 // without killing all of octave.
5794 cm->metis_memory = 0.0;
5795 #endif
5796 #endif
5797 #endif
5798
5799 cm->final_ll = TRUE;
5800
5801 cholmod_sparse Astore;
5802 cholmod_sparse *A = &Astore;
5803 double dummy;
5804 A->nrow = nr;
5805 A->ncol = nc;
5806
5807 A->p = cidx();
5808 A->i = ridx();
5809 A->nzmax = nonzero();
5810 A->packed = TRUE;
5811 A->sorted = TRUE;
5812 A->nz = NULL;
5813 #ifdef IDX_TYPE_LONG
5814 A->itype = CHOLMOD_LONG;
5815 #else
5816 A->itype = CHOLMOD_INT;
5817 #endif
5818 A->dtype = CHOLMOD_DOUBLE;
5819 A->stype = 1;
5820 A->xtype = CHOLMOD_REAL;
5821
5822 if (nr < 1)
5823 A->x = &dummy;
5824 else
5825 A->x = data();
5826
5827 cholmod_dense Bstore;
5828 cholmod_dense *B = &Bstore;
5829 B->nrow = b.rows();
5830 B->ncol = b.cols();
5831 B->d = B->nrow;
5832 B->nzmax = B->nrow * B->ncol;
5833 B->dtype = CHOLMOD_DOUBLE;
5834 B->xtype = CHOLMOD_COMPLEX;
5835 if (nc < 1 || b.cols() < 1)
5836 B->x = &dummy;
5837 else
5838 // We won't alter it, honest :-)
5839 B->x = const_cast<Complex *>(b.fortran_vec());
5840
5841 cholmod_factor *L;
5842 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5843 L = CHOLMOD_NAME(analyze) (A, cm);
5844 CHOLMOD_NAME(factorize) (A, L, cm);
5845 rcond = CHOLMOD_NAME(rcond)(L, cm);
5846 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5847
5848 if (rcond == 0.0)
5849 {
5850 // Either its indefinite or singular. Try UMFPACK
5851 mattype.mark_as_unsymmetric ();
5852 typ = SparseType::Full;
5853 }
5854 else
5855 {
5856 volatile double rcond_plus_one = rcond + 1.0;
5857
5858 if (rcond_plus_one == 1.0 || xisnan (rcond))
5859 {
5860 err = -2;
5861
5862 if (sing_handler)
5863 sing_handler (rcond);
5864 else
5865 (*current_liboctave_error_handler)
5866 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
5867 rcond);
5868
5869 #ifdef HAVE_LSSOLVE
5870 return retval;
5871 #endif
5872 }
5873
5874 cholmod_dense *X;
5875 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5876 X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm);
5877 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5878
5879 retval.resize (b.rows (), b.cols());
5880 for (octave_idx_type j = 0; j < b.cols(); j++)
5881 {
5882 octave_idx_type jr = j * b.rows();
5883 for (octave_idx_type i = 0; i < b.rows(); i++)
5884 retval.xelem(i,j) = static_cast<Complex *>(X->x)[jr + i];
5885 }
5886
5887 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5888 CHOLMOD_NAME(free_dense) (&X, cm);
5889 CHOLMOD_NAME(free_factor) (&L, cm);
5890 CHOLMOD_NAME(finish) (cm);
5891 CHOLMOD_NAME(print_common) (" ", cm);
5892 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
5893 }
5894 #else
5110 (*current_liboctave_warning_handler) 5895 (*current_liboctave_warning_handler)
5111 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 5896 ("CHOLMOD not installed");
5112 5897
5113 mattype.mark_as_unsymmetric (); 5898 mattype.mark_as_unsymmetric ();
5114 typ = SparseType::Full; 5899 typ = SparseType::Full;
5900 #endif
5115 } 5901 }
5116 5902
5117 if (typ == SparseType::Full) 5903 if (typ == SparseType::Full)
5118 { 5904 {
5119 #ifdef HAVE_UMFPACK 5905 #ifdef HAVE_UMFPACK
5221 (*current_liboctave_error_handler) 6007 (*current_liboctave_error_handler)
5222 ("matrix dimension mismatch solution of linear equations"); 6008 ("matrix dimension mismatch solution of linear equations");
5223 else 6009 else
5224 { 6010 {
5225 // Print spparms("spumoni") info if requested 6011 // Print spparms("spumoni") info if requested
5226 int typ = mattype.type (); 6012 volatile int typ = mattype.type ();
5227 mattype.info (); 6013 mattype.info ();
5228 6014
5229 if (typ == SparseType::Hermitian) 6015 if (typ == SparseType::Hermitian)
5230 { 6016 {
5231 // XXX FIXME XXX Write the cholesky solver and only fall 6017 #ifdef HAVE_CHOLMOD
5232 // through if cholesky factorization fails 6018 cholmod_common Common;
5233 6019 cholmod_common *cm = &Common;
6020
6021 // Setup initial parameters
6022 CHOLMOD_NAME(start) (cm);
6023 cm->prefer_zomplex = FALSE;
6024
6025 double spu = Voctave_sparse_controls.get_key ("spumoni");
6026 if (spu == 0.)
6027 {
6028 cm->print = -1;
6029 cm->print_function = NULL;
6030 }
6031 else
6032 {
6033 cm->print = (int)spu + 2;
6034 cm->print_function =&SparseCholPrint;
6035 }
6036
6037 cm->error_handler = &SparseCholError;
6038 cm->complex_divide = CHOLMOD_NAME(divcomplex);
6039 cm->hypotenuse = CHOLMOD_NAME(hypot);
6040
6041 #ifdef HAVE_METIS
6042 // METIS 4.0.1 uses malloc and free, and will terminate MATLAB if
6043 // it runs out of memory. Use CHOLMOD's memory guard for METIS,
6044 // which mxMalloc's a huge block of memory (and then immediately
6045 // mxFree's it) before calling METIS
6046 cm->metis_memory = 2.0;
6047
6048 #if defined(METIS_VERSION)
6049 #if (METIS_VERSION >= METIS_VER(4,0,2))
6050 // METIS 4.0.2 uses function pointers for malloc and free
6051 METIS_malloc = cm->malloc_memory;
6052 METIS_free = cm->free_memory;
6053 // Turn off METIS memory guard. It is not needed, because mxMalloc
6054 // will safely terminate the mexFunction and free any workspace
6055 // without killing all of octave.
6056 cm->metis_memory = 0.0;
6057 #endif
6058 #endif
6059 #endif
6060
6061 cm->final_ll = TRUE;
6062
6063 cholmod_sparse Astore;
6064 cholmod_sparse *A = &Astore;
6065 double dummy;
6066 A->nrow = nr;
6067 A->ncol = nc;
6068
6069 A->p = cidx();
6070 A->i = ridx();
6071 A->nzmax = nonzero();
6072 A->packed = TRUE;
6073 A->sorted = TRUE;
6074 A->nz = NULL;
6075 #ifdef IDX_TYPE_LONG
6076 A->itype = CHOLMOD_LONG;
6077 #else
6078 A->itype = CHOLMOD_INT;
6079 #endif
6080 A->dtype = CHOLMOD_DOUBLE;
6081 A->stype = 1;
6082 A->xtype = CHOLMOD_REAL;
6083
6084 if (nr < 1)
6085 A->x = &dummy;
6086 else
6087 A->x = data();
6088
6089 cholmod_sparse Bstore;
6090 cholmod_sparse *B = &Bstore;
6091 B->nrow = b.rows();
6092 B->ncol = b.cols();
6093 B->p = b.cidx();
6094 B->i = b.ridx();
6095 B->nzmax = b.nonzero();
6096 B->packed = TRUE;
6097 B->sorted = TRUE;
6098 B->nz = NULL;
6099 #ifdef IDX_TYPE_LONG
6100 B->itype = CHOLMOD_LONG;
6101 #else
6102 B->itype = CHOLMOD_INT;
6103 #endif
6104 B->dtype = CHOLMOD_DOUBLE;
6105 B->stype = 0;
6106 B->xtype = CHOLMOD_COMPLEX;
6107
6108 if (b.rows() < 1 || b.cols() < 1)
6109 B->x = &dummy;
6110 else
6111 B->x = b.data();
6112
6113 cholmod_factor *L;
6114 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6115 L = CHOLMOD_NAME(analyze) (A, cm);
6116 CHOLMOD_NAME(factorize) (A, L, cm);
6117 rcond = CHOLMOD_NAME(rcond)(L, cm);
6118 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6119
6120 if (rcond == 0.0)
6121 {
6122 // Either its indefinite or singular. Try UMFPACK
6123 mattype.mark_as_unsymmetric ();
6124 typ = SparseType::Full;
6125 }
6126 else
6127 {
6128 volatile double rcond_plus_one = rcond + 1.0;
6129
6130 if (rcond_plus_one == 1.0 || xisnan (rcond))
6131 {
6132 err = -2;
6133
6134 if (sing_handler)
6135 sing_handler (rcond);
6136 else
6137 (*current_liboctave_error_handler)
6138 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
6139 rcond);
6140
6141 #ifdef HAVE_LSSOLVE
6142 return retval;
6143 #endif
6144 }
6145
6146 cholmod_sparse *X;
6147 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6148 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6149 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6150
6151 retval = SparseComplexMatrix
6152 (static_cast<octave_idx_type>(X->nrow),
6153 static_cast<octave_idx_type>(X->ncol),
6154 static_cast<octave_idx_type>(X->nzmax));
6155 for (octave_idx_type j = 0;
6156 j <= static_cast<octave_idx_type>(X->ncol); j++)
6157 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
6158 for (octave_idx_type j = 0;
6159 j < static_cast<octave_idx_type>(X->nzmax); j++)
6160 {
6161 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
6162 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
6163 }
6164
6165 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6166 CHOLMOD_NAME(free_sparse) (&X, cm);
6167 CHOLMOD_NAME(free_factor) (&L, cm);
6168 CHOLMOD_NAME(finish) (cm);
6169 CHOLMOD_NAME(print_common) (" ", cm);
6170 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6171 }
6172 #else
5234 (*current_liboctave_warning_handler) 6173 (*current_liboctave_warning_handler)
5235 ("SparseMatrix::solve XXX FIXME XXX Cholesky code not done"); 6174 ("CHOLMOD not installed");
5236 6175
5237 mattype.mark_as_unsymmetric (); 6176 mattype.mark_as_unsymmetric ();
5238 typ = SparseType::Full; 6177 typ = SparseType::Full;
6178 #endif
5239 } 6179 }
5240 6180
5241 if (typ == SparseType::Full) 6181 if (typ == SparseType::Full)
5242 { 6182 {
5243 #ifdef HAVE_UMFPACK 6183 #ifdef HAVE_UMFPACK