Mercurial > jwe > octave

comparison liboctave/array/dSparse.cc @ 21136:7cac4e7458f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: clean up code around calls to current_liboctave_error_handler. Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code. * Array-util.cc, Array.cc, CColVector.cc, CDiagMatrix.cc, CMatrix.cc, CNDArray.cc, CRowVector.cc, CSparse.cc, DiagArray2.cc, MArray.cc, PermMatrix.cc, Sparse.cc, Sparse.h, chMatrix.cc, chNDArray.cc, dColVector.cc, dDiagMatrix.cc, dMatrix.cc, dNDArray.cc, dRowVector.cc, dSparse.cc, fCColVector.cc, fCDiagMatrix.cc, fCMatrix.cc, fCNDArray.cc, fCRowVector.cc, fColVector.cc, fDiagMatrix.cc, fMatrix.cc, fNDArray.cc, fRowVector.cc, idx-vector.cc, CmplxAEPBAL.cc, CmplxCHOL.cc, CmplxGEPBAL.cc, CmplxHESS.cc, CmplxLU.cc, CmplxQR.cc, CmplxSCHUR.cc, CmplxSVD.cc, DASPK.cc, EIG.cc, LSODE.cc, Quad.cc, SparseCmplxCHOL.cc, SparseCmplxLU.cc, SparseCmplxQR.cc, SparseQR.cc, SparsedbleCHOL.cc, SparsedbleLU.cc, base-lu.cc, bsxfun-defs.cc, dbleAEPBAL.cc, dbleCHOL.cc, dbleGEPBAL.cc, dbleHESS.cc, dbleLU.cc, dbleQR.cc, dbleSCHUR.cc, dbleSVD.cc, eigs-base.cc, fCmplxAEPBAL.cc, fCmplxCHOL.cc, fCmplxLU.cc, fCmplxQR.cc, fCmplxSCHUR.cc, fEIG.cc, floatAEPBAL.cc, floatCHOL.cc, floatGEPBAL.cc, floatHESS.cc, floatLU.cc, floatQR.cc, floatSCHUR.cc, floatSVD.cc, lo-specfun.cc, oct-fftw.cc, oct-rand.cc, oct-spparms.cc, sparse-base-chol.cc, sparse-dmsolve.cc, file-ops.cc, lo-sysdep.cc, mach-info.cc, oct-env.cc, oct-syscalls.cc, cmd-edit.cc, cmd-hist.cc, data-conv.cc, lo-ieee.cc, lo-regexp.cc, oct-base64.cc, oct-shlib.cc, pathsearch.cc, singleton-cleanup.cc, sparse-util.cc, unwind-prot.cc: Remove statements after call to handler that are no longer reachable. Place input validation first and immediately call handler if necessary. Change if/error_handler/else to if/error_handler and re-indent code.
author Rik <rik@octave.org>
date Sat, 23 Jan 2016 13:52:03 -0800
parents 499b851fbfae
children 623fc7d08cc6
comparison
equal deleted inserted replaced
21135:95da3bc8a281 21136:7cac4e7458f2
749 } 749 }
750 750
751 SparseMatrix 751 SparseMatrix
752 atan2 (const SparseMatrix& x, const SparseMatrix& y) 752 atan2 (const SparseMatrix& x, const SparseMatrix& y)
753 { 753 {
754 if ((x.rows () != y.rows ()) || (x.cols () != y.cols ()))
755 (*current_liboctave_error_handler) ("matrix size mismatch");
756
757 octave_idx_type x_nr = x.rows ();
758 octave_idx_type x_nc = x.cols ();
759
760 octave_idx_type y_nr = y.rows ();
761 octave_idx_type y_nc = y.cols ();
762
763 if (x_nr != y_nr || x_nc != y_nc)
764 err_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc);
765
754 SparseMatrix r; 766 SparseMatrix r;
755 767
756 if ((x.rows () == y.rows ()) && (x.cols () == y.cols ())) 768 r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ()));
769
770 octave_idx_type jx = 0;
771 r.cidx (0) = 0;
772 for (octave_idx_type i = 0 ; i < x_nc ; i++)
757 { 773 {
758 octave_idx_type x_nr = x.rows (); 774 octave_idx_type ja = x.cidx (i);
759 octave_idx_type x_nc = x.cols (); 775 octave_idx_type ja_max = x.cidx (i+1);
760 776 bool ja_lt_max= ja < ja_max;
761 octave_idx_type y_nr = y.rows (); 777
762 octave_idx_type y_nc = y.cols (); 778 octave_idx_type jb = y.cidx (i);
763 779 octave_idx_type jb_max = y.cidx (i+1);
764 if (x_nr != y_nr || x_nc != y_nc) 780 bool jb_lt_max = jb < jb_max;
765 err_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc); 781
766 782 while (ja_lt_max || jb_lt_max)
767 r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ())); 783 {
768 784 octave_quit ();
769 octave_idx_type jx = 0; 785 if ((! jb_lt_max)
770 r.cidx (0) = 0; 786 || (ja_lt_max && (x.ridx (ja) < y.ridx (jb))))
771 for (octave_idx_type i = 0 ; i < x_nc ; i++) 787 {
772 { 788 r.ridx (jx) = x.ridx (ja);
773 octave_idx_type ja = x.cidx (i); 789 r.data (jx) = atan2 (x.data (ja), 0.);
774 octave_idx_type ja_max = x.cidx (i+1); 790 jx++;
775 bool ja_lt_max= ja < ja_max; 791 ja++;
776 792 ja_lt_max= ja < ja_max;
777 octave_idx_type jb = y.cidx (i); 793 }
778 octave_idx_type jb_max = y.cidx (i+1); 794 else if ((! ja_lt_max)
779 bool jb_lt_max = jb < jb_max; 795 || (jb_lt_max && (y.ridx (jb) < x.ridx (ja))))
780 796 {
781 while (ja_lt_max || jb_lt_max) 797 jb++;
782 { 798 jb_lt_max= jb < jb_max;
783 octave_quit (); 799 }
784 if ((! jb_lt_max) 800 else
785 || (ja_lt_max && (x.ridx (ja) < y.ridx (jb)))) 801 {
786 { 802 double tmp = atan2 (x.data (ja), y.data (jb));
803 if (tmp != 0.)
804 {
805 r.data (jx) = tmp;
787 r.ridx (jx) = x.ridx (ja); 806 r.ridx (jx) = x.ridx (ja);
788 r.data (jx) = atan2 (x.data (ja), 0.);
789 jx++; 807 jx++;
790 ja++; 808 }
791 ja_lt_max= ja < ja_max; 809 ja++;
792 } 810 ja_lt_max= ja < ja_max;
793 else if ((! ja_lt_max) 811 jb++;
794 || (jb_lt_max && (y.ridx (jb) < x.ridx (ja)))) 812 jb_lt_max= jb < jb_max;
795 { 813 }
796 jb++; 814 }
797 jb_lt_max= jb < jb_max; 815 r.cidx (i+1) = jx;
798 }
799 else
800 {
801 double tmp = atan2 (x.data (ja), y.data (jb));
802 if (tmp != 0.)
803 {
804 r.data (jx) = tmp;
805 r.ridx (jx) = x.ridx (ja);
806 jx++;
807 }
808 ja++;
809 ja_lt_max= ja < ja_max;
810 jb++;
811 jb_lt_max= jb < jb_max;
812 }
813 }
814 r.cidx (i+1) = jx;
815 }
816
817 r.maybe_compress ();
818 } 816 }
819 else 817
820 (*current_liboctave_error_handler) ("matrix size mismatch"); 818 r.maybe_compress ();
821 819
822 return r; 820 return r;
823 } 821 }
824 822
825 SparseMatrix 823 SparseMatrix
857 octave_idx_type nc = cols (); 855 octave_idx_type nc = cols ();
858 info = 0; 856 info = 0;
859 857
860 if (nr == 0 || nc == 0 || nr != nc) 858 if (nr == 0 || nc == 0 || nr != nc)
861 (*current_liboctave_error_handler) ("inverse requires square matrix"); 859 (*current_liboctave_error_handler) ("inverse requires square matrix");
860
861 // Print spparms("spumoni") info if requested
862 int typ = mattyp.type ();
863 mattyp.info ();
864
865 if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
866 (*current_liboctave_error_handler) ("incorrect matrix type");
867
868 if (typ == MatrixType::Permuted_Diagonal)
869 retval = transpose ();
862 else 870 else
871 retval = *this;
872
873 // Force make_unique to be called
874 double *v = retval.data ();
875
876 if (calccond)
863 { 877 {
864 // Print spparms("spumoni") info if requested 878 double dmax = 0.;
865 int typ = mattyp.type (); 879 double dmin = octave_Inf;
866 mattyp.info (); 880 for (octave_idx_type i = 0; i < nr; i++)
867 881 {
868 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) 882 double tmp = fabs (v[i]);
869 { 883 if (tmp > dmax)
870 if (typ == MatrixType::Permuted_Diagonal) 884 dmax = tmp;
871 retval = transpose (); 885 if (tmp < dmin)
872 else 886 dmin = tmp;
873 retval = *this; 887 }
874 888 rcond = dmin / dmax;
875 // Force make_unique to be called
876 double *v = retval.data ();
877
878 if (calccond)
879 {
880 double dmax = 0.;
881 double dmin = octave_Inf;
882 for (octave_idx_type i = 0; i < nr; i++)
883 {
884 double tmp = fabs (v[i]);
885 if (tmp > dmax)
886 dmax = tmp;
887 if (tmp < dmin)
888 dmin = tmp;
889 }
890 rcond = dmin / dmax;
891 }
892
893 for (octave_idx_type i = 0; i < nr; i++)
894 v[i] = 1.0 / v[i];
895 }
896 else
897 (*current_liboctave_error_handler) ("incorrect matrix type");
898 } 889 }
890
891 for (octave_idx_type i = 0; i < nr; i++)
892 v[i] = 1.0 / v[i];
899 893
900 return retval; 894 return retval;
901 } 895 }
902 896
903 SparseMatrix 897 SparseMatrix
911 octave_idx_type nc = cols (); 905 octave_idx_type nc = cols ();
912 info = 0; 906 info = 0;
913 907
914 if (nr == 0 || nc == 0 || nr != nc) 908 if (nr == 0 || nc == 0 || nr != nc)
915 (*current_liboctave_error_handler) ("inverse requires square matrix"); 909 (*current_liboctave_error_handler) ("inverse requires square matrix");
910
911 // Print spparms("spumoni") info if requested
912 int typ = mattyp.type ();
913 mattyp.info ();
914
915 if (typ != MatrixType::Upper && typ != MatrixType::Permuted_Upper
916 && typ != MatrixType::Lower && typ != MatrixType::Permuted_Lower)
917 (*current_liboctave_error_handler) ("incorrect matrix type");
918
919 double anorm = 0.;
920 double ainvnorm = 0.;
921
922 if (calccond)
923 {
924 // Calculate the 1-norm of matrix for rcond calculation
925 for (octave_idx_type j = 0; j < nr; j++)
926 {
927 double atmp = 0.;
928 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
929 atmp += fabs (data (i));
930 if (atmp > anorm)
931 anorm = atmp;
932 }
933 }
934
935 if (typ == MatrixType::Upper || typ == MatrixType::Lower)
936 {
937 octave_idx_type nz = nnz ();
938 octave_idx_type cx = 0;
939 octave_idx_type nz2 = nz;
940 retval = SparseMatrix (nr, nc, nz2);
941
942 for (octave_idx_type i = 0; i < nr; i++)
943 {
944 octave_quit ();
945 // place the 1 in the identity position
946 octave_idx_type cx_colstart = cx;
947
948 if (cx == nz2)
949 {
950 nz2 *= 2;
951 retval.change_capacity (nz2);
952 }
953
954 retval.xcidx (i) = cx;
955 retval.xridx (cx) = i;
956 retval.xdata (cx) = 1.0;
957 cx++;
958
959 // iterate accross columns of input matrix
960 for (octave_idx_type j = i+1; j < nr; j++)
961 {
962 double v = 0.;
963 // iterate to calculate sum
964 octave_idx_type colXp = retval.xcidx (i);
965 octave_idx_type colUp = cidx (j);
966 octave_idx_type rpX, rpU;
967
968 if (cidx (j) == cidx (j+1))
969 (*current_liboctave_error_handler) ("division by zero");
970
971 do
972 {
973 octave_quit ();
974 rpX = retval.xridx (colXp);
975 rpU = ridx (colUp);
976
977 if (rpX < rpU)
978 colXp++;
979 else if (rpX > rpU)
980 colUp++;
981 else
982 {
983 v -= retval.xdata (colXp) * data (colUp);
984 colXp++;
985 colUp++;
986 }
987 }
988 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
989
990 // get A(m,m)
991 if (typ == MatrixType::Upper)
992 colUp = cidx (j+1) - 1;
993 else
994 colUp = cidx (j);
995 double pivot = data (colUp);
996 if (pivot == 0. || ridx (colUp) != j)
997 (*current_liboctave_error_handler) ("division by zero");
998
999 if (v != 0.)
1000 {
1001 if (cx == nz2)
1002 {
1003 nz2 *= 2;
1004 retval.change_capacity (nz2);
1005 }
1006
1007 retval.xridx (cx) = j;
1008 retval.xdata (cx) = v / pivot;
1009 cx++;
1010 }
1011 }
1012
1013 // get A(m,m)
1014 octave_idx_type colUp;
1015 if (typ == MatrixType::Upper)
1016 colUp = cidx (i+1) - 1;
1017 else
1018 colUp = cidx (i);
1019 double pivot = data (colUp);
1020 if (pivot == 0. || ridx (colUp) != i)
1021 (*current_liboctave_error_handler) ("division by zero");
1022
1023 if (pivot != 1.0)
1024 for (octave_idx_type j = cx_colstart; j < cx; j++)
1025 retval.xdata (j) /= pivot;
1026 }
1027 retval.xcidx (nr) = cx;
1028 retval.maybe_compress ();
1029 }
916 else 1030 else
917 { 1031 {
918 // Print spparms("spumoni") info if requested 1032 octave_idx_type nz = nnz ();
919 int typ = mattyp.type (); 1033 octave_idx_type cx = 0;
920 mattyp.info (); 1034 octave_idx_type nz2 = nz;
921 1035 retval = SparseMatrix (nr, nc, nz2);
922 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper 1036
923 || typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) 1037 OCTAVE_LOCAL_BUFFER (double, work, nr);
924 { 1038 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr);
925 double anorm = 0.; 1039
926 double ainvnorm = 0.; 1040 octave_idx_type *perm = mattyp.triangular_perm ();
927 1041 if (typ == MatrixType::Permuted_Upper)
928 if (calccond) 1042 {
929 { 1043 for (octave_idx_type i = 0; i < nr; i++)
930 // Calculate the 1-norm of matrix for rcond calculation 1044 rperm[perm[i]] = i;
931 for (octave_idx_type j = 0; j < nr; j++) 1045 }
932 { 1046 else
933 double atmp = 0.; 1047 {
934 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 1048 for (octave_idx_type i = 0; i < nr; i++)
935 atmp += fabs (data (i)); 1049 rperm[i] = perm[i];
936 if (atmp > anorm) 1050 for (octave_idx_type i = 0; i < nr; i++)
937 anorm = atmp; 1051 perm[rperm[i]] = i;
938 } 1052 }
939 } 1053
940 1054 for (octave_idx_type i = 0; i < nr; i++)
941 if (typ == MatrixType::Upper || typ == MatrixType::Lower) 1055 {
942 { 1056 octave_quit ();
943 octave_idx_type nz = nnz (); 1057 octave_idx_type iidx = rperm[i];
944 octave_idx_type cx = 0; 1058
945 octave_idx_type nz2 = nz; 1059 for (octave_idx_type j = 0; j < nr; j++)
946 retval = SparseMatrix (nr, nc, nz2); 1060 work[j] = 0.;
947 1061
948 for (octave_idx_type i = 0; i < nr; i++) 1062 // place the 1 in the identity position
1063 work[iidx] = 1.0;
1064
1065 // iterate accross columns of input matrix
1066 for (octave_idx_type j = iidx+1; j < nr; j++)
1067 {
1068 double v = 0.;
1069 octave_idx_type jidx = perm[j];
1070 // iterate to calculate sum
1071 for (octave_idx_type k = cidx (jidx);
1072 k < cidx (jidx+1); k++)
949 { 1073 {
950 octave_quit (); 1074 octave_quit ();
951 // place the 1 in the identity position 1075 v -= work[ridx (k)] * data (k);
952 octave_idx_type cx_colstart = cx; 1076 }
953 1077
954 if (cx == nz2) 1078 // get A(m,m)
955 { 1079 double pivot;
956 nz2 *= 2; 1080 if (typ == MatrixType::Permuted_Upper)
957 retval.change_capacity (nz2); 1081 pivot = data (cidx (jidx+1) - 1);
958 } 1082 else
959 1083 pivot = data (cidx (jidx));
960 retval.xcidx (i) = cx; 1084 if (pivot == 0.)
961 retval.xridx (cx) = i; 1085 (*current_liboctave_error_handler) ("division by zero");
962 retval.xdata (cx) = 1.0; 1086
963 cx++; 1087 work[j] = v / pivot;
964 1088 }
965 // iterate accross columns of input matrix 1089
966 for (octave_idx_type j = i+1; j < nr; j++) 1090 // get A(m,m)
967 { 1091 octave_idx_type colUp;
968 double v = 0.; 1092 if (typ == MatrixType::Permuted_Upper)
969 // iterate to calculate sum 1093 colUp = cidx (perm[iidx]+1) - 1;
970 octave_idx_type colXp = retval.xcidx (i);
971 octave_idx_type colUp = cidx (j);
972 octave_idx_type rpX, rpU;
973
974 if (cidx (j) == cidx (j+1))
975 {
976 (*current_liboctave_error_handler)
977 ("division by zero");
978 goto inverse_singular;
979 }
980
981 do
982 {
983 octave_quit ();
984 rpX = retval.xridx (colXp);
985 rpU = ridx (colUp);
986
987 if (rpX < rpU)
988 colXp++;
989 else if (rpX > rpU)
990 colUp++;
991 else
992 {
993 v -= retval.xdata (colXp) * data (colUp);
994 colXp++;
995 colUp++;
996 }
997 }
998 while (rpX < j && rpU < j && colXp < cx && colUp < nz);
999
1000 // get A(m,m)
1001 if (typ == MatrixType::Upper)
1002 colUp = cidx (j+1) - 1;
1003 else
1004 colUp = cidx (j);
1005 double pivot = data (colUp);
1006 if (pivot == 0. || ridx (colUp) != j)
1007 {
1008 (*current_liboctave_error_handler)
1009 ("division by zero");
1010 goto inverse_singular;
1011 }
1012
1013 if (v != 0.)
1014 {
1015 if (cx == nz2)
1016 {
1017 nz2 *= 2;
1018 retval.change_capacity (nz2);
1019 }
1020
1021 retval.xridx (cx) = j;
1022 retval.xdata (cx) = v / pivot;
1023 cx++;
1024 }
1025 }
1026
1027 // get A(m,m)
1028 octave_idx_type colUp;
1029 if (typ == MatrixType::Upper)
1030 colUp = cidx (i+1) - 1;
1031 else
1032 colUp = cidx (i);
1033 double pivot = data (colUp);
1034 if (pivot == 0. || ridx (colUp) != i)
1035 {
1036 (*current_liboctave_error_handler) ("division by zero");
1037 goto inverse_singular;
1038 }
1039
1040 if (pivot != 1.0)
1041 for (octave_idx_type j = cx_colstart; j < cx; j++)
1042 retval.xdata (j) /= pivot;
1043 }
1044 retval.xcidx (nr) = cx;
1045 retval.maybe_compress ();
1046 }
1047 else 1094 else
1048 { 1095 colUp = cidx (perm[iidx]);
1049 octave_idx_type nz = nnz (); 1096
1050 octave_idx_type cx = 0; 1097 double pivot = data (colUp);
1051 octave_idx_type nz2 = nz; 1098 if (pivot == 0.)
1052 retval = SparseMatrix (nr, nc, nz2); 1099 (*current_liboctave_error_handler) ("division by zero");
1053 1100
1054 OCTAVE_LOCAL_BUFFER (double, work, nr); 1101 octave_idx_type new_cx = cx;
1055 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); 1102 for (octave_idx_type j = iidx; j < nr; j++)
1056 1103 if (work[j] != 0.0)
1057 octave_idx_type *perm = mattyp.triangular_perm (); 1104 {
1058 if (typ == MatrixType::Permuted_Upper) 1105 new_cx++;
1059 { 1106 if (pivot != 1.0)
1060 for (octave_idx_type i = 0; i < nr; i++) 1107 work[j] /= pivot;
1061 rperm[perm[i]] = i; 1108 }
1062 } 1109
1063 else 1110 if (cx < new_cx)
1064 { 1111 {
1065 for (octave_idx_type i = 0; i < nr; i++) 1112 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1066 rperm[i] = perm[i]; 1113 retval.change_capacity (nz2);
1067 for (octave_idx_type i = 0; i < nr; i++) 1114 }
1068 perm[rperm[i]] = i; 1115
1069 } 1116 retval.xcidx (i) = cx;
1070 1117 for (octave_idx_type j = iidx; j < nr; j++)
1071 for (octave_idx_type i = 0; i < nr; i++) 1118 if (work[j] != 0.)
1072 { 1119 {
1073 octave_quit (); 1120 retval.xridx (cx) = j;
1074 octave_idx_type iidx = rperm[i]; 1121 retval.xdata (cx++) = work[j];
1075 1122 }
1076 for (octave_idx_type j = 0; j < nr; j++) 1123 }
1077 work[j] = 0.; 1124
1078 1125 retval.xcidx (nr) = cx;
1079 // place the 1 in the identity position 1126 retval.maybe_compress ();
1080 work[iidx] = 1.0; 1127 }
1081 1128
1082 // iterate accross columns of input matrix 1129 if (calccond)
1083 for (octave_idx_type j = iidx+1; j < nr; j++) 1130 {
1084 { 1131 // Calculate the 1-norm of inverse matrix for rcond calculation
1085 double v = 0.; 1132 for (octave_idx_type j = 0; j < nr; j++)
1086 octave_idx_type jidx = perm[j]; 1133 {
1087 // iterate to calculate sum 1134 double atmp = 0.;
1088 for (octave_idx_type k = cidx (jidx); 1135 for (octave_idx_type i = retval.cidx (j);
1089 k < cidx (jidx+1); k++) 1136 i < retval.cidx (j+1); i++)
1090 { 1137 atmp += fabs (retval.data (i));
1091 octave_quit (); 1138 if (atmp > ainvnorm)
1092 v -= work[ridx (k)] * data (k); 1139 ainvnorm = atmp;
1093 } 1140 }
1094 1141
1095 // get A(m,m) 1142 rcond = 1. / ainvnorm / anorm;
1096 double pivot;
1097 if (typ == MatrixType::Permuted_Upper)
1098 pivot = data (cidx (jidx+1) - 1);
1099 else
1100 pivot = data (cidx (jidx));
1101 if (pivot == 0.)
1102 {
1103 (*current_liboctave_error_handler)
1104 ("division by zero");
1105 goto inverse_singular;
1106 }
1107
1108 work[j] = v / pivot;
1109 }
1110
1111 // get A(m,m)
1112 octave_idx_type colUp;
1113 if (typ == MatrixType::Permuted_Upper)
1114 colUp = cidx (perm[iidx]+1) - 1;
1115 else
1116 colUp = cidx (perm[iidx]);
1117
1118 double pivot = data (colUp);
1119 if (pivot == 0.)
1120 {
1121 (*current_liboctave_error_handler)
1122 ("division by zero");
1123 goto inverse_singular;
1124 }
1125
1126 octave_idx_type new_cx = cx;
1127 for (octave_idx_type j = iidx; j < nr; j++)
1128 if (work[j] != 0.0)
1129 {
1130 new_cx++;
1131 if (pivot != 1.0)
1132 work[j] /= pivot;
1133 }
1134
1135 if (cx < new_cx)
1136 {
1137 nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2);
1138 retval.change_capacity (nz2);
1139 }
1140
1141 retval.xcidx (i) = cx;
1142 for (octave_idx_type j = iidx; j < nr; j++)
1143 if (work[j] != 0.)
1144 {
1145 retval.xridx (cx) = j;
1146 retval.xdata (cx++) = work[j];
1147 }
1148 }
1149
1150 retval.xcidx (nr) = cx;
1151 retval.maybe_compress ();
1152 }
1153
1154 if (calccond)
1155 {
1156 // Calculate the 1-norm of inverse matrix for rcond calculation
1157 for (octave_idx_type j = 0; j < nr; j++)
1158 {
1159 double atmp = 0.;
1160 for (octave_idx_type i = retval.cidx (j);
1161 i < retval.cidx (j+1); i++)
1162 atmp += fabs (retval.data (i));
1163 if (atmp > ainvnorm)
1164 ainvnorm = atmp;
1165 }
1166
1167 rcond = 1. / ainvnorm / anorm;
1168 }
1169 }
1170 else
1171 (*current_liboctave_error_handler) ("incorrect matrix type");
1172 } 1143 }
1173 1144
1174 return retval; 1145 return retval;
1175 1146
1176 inverse_singular: 1147 inverse_singular:
1310 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, 1281 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai,
1311 Ax, 0, &Symbolic, control, info); 1282 Ax, 0, &Symbolic, control, info);
1312 1283
1313 if (status < 0) 1284 if (status < 0)
1314 { 1285 {
1286 UMFPACK_DNAME (report_status) (control, status);
1287 UMFPACK_DNAME (report_info) (control, info);
1288
1289 UMFPACK_DNAME (free_symbolic) (&Symbolic);
1290
1315 (*current_liboctave_error_handler) 1291 (*current_liboctave_error_handler)
1316 ("SparseMatrix::determinant symbolic factorization failed"); 1292 ("SparseMatrix::determinant symbolic factorization failed");
1317
1318 UMFPACK_DNAME (report_status) (control, status);
1319 UMFPACK_DNAME (report_info) (control, info);
1320
1321 UMFPACK_DNAME (free_symbolic) (&Symbolic);
1322 } 1293 }
1323 else 1294 else
1324 { 1295 {
1325 UMFPACK_DNAME (report_symbolic) (Symbolic, control); 1296 UMFPACK_DNAME (report_symbolic) (Symbolic, control);
1326 1297
1331 1302
1332 rcond = Info (UMFPACK_RCOND); 1303 rcond = Info (UMFPACK_RCOND);
1333 1304
1334 if (status < 0) 1305 if (status < 0)
1335 { 1306 {
1307 UMFPACK_DNAME (report_status) (control, status);
1308 UMFPACK_DNAME (report_info) (control, info);
1309
1310 UMFPACK_DNAME (free_numeric) (&Numeric);
1336 (*current_liboctave_error_handler) 1311 (*current_liboctave_error_handler)
1337 ("SparseMatrix::determinant numeric factorization failed"); 1312 ("SparseMatrix::determinant numeric factorization failed");
1338
1339 UMFPACK_DNAME (report_status) (control, status);
1340 UMFPACK_DNAME (report_info) (control, info);
1341
1342 UMFPACK_DNAME (free_numeric) (&Numeric);
1343 } 1313 }
1344 else 1314 else
1345 { 1315 {
1346 UMFPACK_DNAME (report_numeric) (Numeric, control); 1316 UMFPACK_DNAME (report_numeric) (Numeric, control);
1347 1317
1350 status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric, 1320 status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric,
1351 info); 1321 info);
1352 1322
1353 if (status < 0) 1323 if (status < 0)
1354 { 1324 {
1325 UMFPACK_DNAME (report_status) (control, status);
1326 UMFPACK_DNAME (report_info) (control, info);
1327
1355 (*current_liboctave_error_handler) 1328 (*current_liboctave_error_handler)
1356 ("SparseMatrix::determinant error calculating determinant"); 1329 ("SparseMatrix::determinant error calculating determinant");
1357
1358 UMFPACK_DNAME (report_status) (control, status);
1359 UMFPACK_DNAME (report_info) (control, info);
1360 } 1330 }
1361 else 1331 else
1362 retval = DET (c10, e10, 10); 1332 retval = DET (c10, e10, 10);
1363 1333
1364 UMFPACK_DNAME (free_numeric) (&Numeric); 1334 UMFPACK_DNAME (free_numeric) (&Numeric);
1389 err = 0; 1359 err = 0;
1390 1360
1391 if (nr != b.rows ()) 1361 if (nr != b.rows ())
1392 (*current_liboctave_error_handler) 1362 (*current_liboctave_error_handler)
1393 ("matrix dimension mismatch solution of linear equations"); 1363 ("matrix dimension mismatch solution of linear equations");
1394 else if (nr == 0 || nc == 0 || b.cols () == 0) 1364
1365 if (nr == 0 || nc == 0 || b.cols () == 0)
1395 retval = Matrix (nc, b.cols (), 0.0); 1366 retval = Matrix (nc, b.cols (), 0.0);
1396 else 1367 else
1397 { 1368 {
1398 // Print spparms("spumoni") info if requested 1369 // Print spparms("spumoni") info if requested
1399 int typ = mattype.type (); 1370 int typ = mattype.type ();
1400 mattype.info (); 1371 mattype.info ();
1401 1372
1402 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) 1373 if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
1403 { 1374 (*current_liboctave_error_handler) ("incorrect matrix type");
1404 retval.resize (nc, b.cols (), 0.); 1375
1405 if (typ == MatrixType::Diagonal) 1376 retval.resize (nc, b.cols (), 0.);
1406 for (octave_idx_type j = 0; j < b.cols (); j++) 1377 if (typ == MatrixType::Diagonal)
1407 for (octave_idx_type i = 0; i < nm; i++) 1378 for (octave_idx_type j = 0; j < b.cols (); j++)
1408 retval(i,j) = b(i,j) / data (i); 1379 for (octave_idx_type i = 0; i < nm; i++)
1409 else 1380 retval(i,j) = b(i,j) / data (i);
1410 for (octave_idx_type j = 0; j < b.cols (); j++)
1411 for (octave_idx_type k = 0; k < nc; k++)
1412 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1413 retval(k,j) = b(ridx (i),j) / data (i);
1414
1415 if (calc_cond)
1416 {
1417 double dmax = 0.;
1418 double dmin = octave_Inf;
1419 for (octave_idx_type i = 0; i < nm; i++)
1420 {
1421 double tmp = fabs (data (i));
1422 if (tmp > dmax)
1423 dmax = tmp;
1424 if (tmp < dmin)
1425 dmin = tmp;
1426 }
1427 rcond = dmin / dmax;
1428 }
1429 else
1430 rcond = 1.;
1431 }
1432 else 1381 else
1433 (*current_liboctave_error_handler) ("incorrect matrix type"); 1382 for (octave_idx_type j = 0; j < b.cols (); j++)
1383 for (octave_idx_type k = 0; k < nc; k++)
1384 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1385 retval(k,j) = b(ridx (i),j) / data (i);
1386
1387 if (calc_cond)
1388 {
1389 double dmax = 0.;
1390 double dmin = octave_Inf;
1391 for (octave_idx_type i = 0; i < nm; i++)
1392 {
1393 double tmp = fabs (data (i));
1394 if (tmp > dmax)
1395 dmax = tmp;
1396 if (tmp < dmin)
1397 dmin = tmp;
1398 }
1399 rcond = dmin / dmax;
1400 }
1401 else
1402 rcond = 1.;
1434 } 1403 }
1435 1404
1436 return retval; 1405 return retval;
1437 } 1406 }
1438 1407
1449 err = 0; 1418 err = 0;
1450 1419
1451 if (nr != b.rows ()) 1420 if (nr != b.rows ())
1452 (*current_liboctave_error_handler) 1421 (*current_liboctave_error_handler)
1453 ("matrix dimension mismatch solution of linear equations"); 1422 ("matrix dimension mismatch solution of linear equations");
1454 else if (nr == 0 || nc == 0 || b.cols () == 0) 1423
1424 if (nr == 0 || nc == 0 || b.cols () == 0)
1455 retval = SparseMatrix (nc, b.cols ()); 1425 retval = SparseMatrix (nc, b.cols ());
1456 else 1426 else
1457 { 1427 {
1458 // Print spparms("spumoni") info if requested 1428 // Print spparms("spumoni") info if requested
1459 int typ = mattype.type (); 1429 int typ = mattype.type ();
1460 mattype.info (); 1430 mattype.info ();
1461 1431
1462 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) 1432 if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
1463 { 1433 (*current_liboctave_error_handler) ("incorrect matrix type");
1464 octave_idx_type b_nc = b.cols (); 1434
1465 octave_idx_type b_nz = b.nnz (); 1435 octave_idx_type b_nc = b.cols ();
1466 retval = SparseMatrix (nc, b_nc, b_nz); 1436 octave_idx_type b_nz = b.nnz ();
1467 1437 retval = SparseMatrix (nc, b_nc, b_nz);
1468 retval.xcidx (0) = 0; 1438
1469 octave_idx_type ii = 0; 1439 retval.xcidx (0) = 0;
1470 if (typ == MatrixType::Diagonal) 1440 octave_idx_type ii = 0;
1471 for (octave_idx_type j = 0; j < b_nc; j++) 1441 if (typ == MatrixType::Diagonal)
1442 for (octave_idx_type j = 0; j < b_nc; j++)
1443 {
1444 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1472 { 1445 {
1473 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 1446 if (b.ridx (i) >= nm)
1474 { 1447 break;
1475 if (b.ridx (i) >= nm) 1448 retval.xridx (ii) = b.ridx (i);
1476 break; 1449 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1477 retval.xridx (ii) = b.ridx (i);
1478 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1479 }
1480 retval.xcidx (j+1) = ii;
1481 } 1450 }
1482 else 1451 retval.xcidx (j+1) = ii;
1483 for (octave_idx_type j = 0; j < b_nc; j++) 1452 }
1484 {
1485 for (octave_idx_type l = 0; l < nc; l++)
1486 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1487 {
1488 bool found = false;
1489 octave_idx_type k;
1490 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1491 if (ridx (i) == b.ridx (k))
1492 {
1493 found = true;
1494 break;
1495 }
1496 if (found)
1497 {
1498 retval.xridx (ii) = l;
1499 retval.xdata (ii++) = b.data (k) / data (i);
1500 }
1501 }
1502 retval.xcidx (j+1) = ii;
1503 }
1504
1505 if (calc_cond)
1506 {
1507 double dmax = 0.;
1508 double dmin = octave_Inf;
1509 for (octave_idx_type i = 0; i < nm; i++)
1510 {
1511 double tmp = fabs (data (i));
1512 if (tmp > dmax)
1513 dmax = tmp;
1514 if (tmp < dmin)
1515 dmin = tmp;
1516 }
1517 rcond = dmin / dmax;
1518 }
1519 else
1520 rcond = 1.;
1521 }
1522 else 1453 else
1523 (*current_liboctave_error_handler) ("incorrect matrix type"); 1454 for (octave_idx_type j = 0; j < b_nc; j++)
1455 {
1456 for (octave_idx_type l = 0; l < nc; l++)
1457 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1458 {
1459 bool found = false;
1460 octave_idx_type k;
1461 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1462 if (ridx (i) == b.ridx (k))
1463 {
1464 found = true;
1465 break;
1466 }
1467 if (found)
1468 {
1469 retval.xridx (ii) = l;
1470 retval.xdata (ii++) = b.data (k) / data (i);
1471 }
1472 }
1473 retval.xcidx (j+1) = ii;
1474 }
1475
1476 if (calc_cond)
1477 {
1478 double dmax = 0.;
1479 double dmin = octave_Inf;
1480 for (octave_idx_type i = 0; i < nm; i++)
1481 {
1482 double tmp = fabs (data (i));
1483 if (tmp > dmax)
1484 dmax = tmp;
1485 if (tmp < dmin)
1486 dmin = tmp;
1487 }
1488 rcond = dmin / dmax;
1489 }
1490 else
1491 rcond = 1.;
1524 } 1492 }
1525 1493
1526 return retval; 1494 return retval;
1527 } 1495 }
1528 1496
1539 err = 0; 1507 err = 0;
1540 1508
1541 if (nr != b.rows ()) 1509 if (nr != b.rows ())
1542 (*current_liboctave_error_handler) 1510 (*current_liboctave_error_handler)
1543 ("matrix dimension mismatch solution of linear equations"); 1511 ("matrix dimension mismatch solution of linear equations");
1544 else if (nr == 0 || nc == 0 || b.cols () == 0) 1512
1513 if (nr == 0 || nc == 0 || b.cols () == 0)
1545 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 1514 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
1546 else 1515 else
1547 { 1516 {
1548 // Print spparms("spumoni") info if requested 1517 // Print spparms("spumoni") info if requested
1549 int typ = mattype.type (); 1518 int typ = mattype.type ();
1550 mattype.info (); 1519 mattype.info ();
1551 1520
1552 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) 1521 if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
1553 { 1522 (*current_liboctave_error_handler) ("incorrect matrix type");
1554 retval.resize (nc, b.cols (), 0); 1523
1555 if (typ == MatrixType::Diagonal) 1524 retval.resize (nc, b.cols (), 0);
1556 for (octave_idx_type j = 0; j < b.cols (); j++) 1525 if (typ == MatrixType::Diagonal)
1557 for (octave_idx_type i = 0; i < nm; i++) 1526 for (octave_idx_type j = 0; j < b.cols (); j++)
1558 retval(i,j) = b(i,j) / data (i); 1527 for (octave_idx_type i = 0; i < nm; i++)
1559 else 1528 retval(i,j) = b(i,j) / data (i);
1560 for (octave_idx_type j = 0; j < b.cols (); j++)
1561 for (octave_idx_type k = 0; k < nc; k++)
1562 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1563 retval(k,j) = b(ridx (i),j) / data (i);
1564
1565 if (calc_cond)
1566 {
1567 double dmax = 0.;
1568 double dmin = octave_Inf;
1569 for (octave_idx_type i = 0; i < nm; i++)
1570 {
1571 double tmp = fabs (data (i));
1572 if (tmp > dmax)
1573 dmax = tmp;
1574 if (tmp < dmin)
1575 dmin = tmp;
1576 }
1577 rcond = dmin / dmax;
1578 }
1579 else
1580 rcond = 1.;
1581 }
1582 else 1529 else
1583 (*current_liboctave_error_handler) ("incorrect matrix type"); 1530 for (octave_idx_type j = 0; j < b.cols (); j++)
1531 for (octave_idx_type k = 0; k < nc; k++)
1532 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1533 retval(k,j) = b(ridx (i),j) / data (i);
1534
1535 if (calc_cond)
1536 {
1537 double dmax = 0.;
1538 double dmin = octave_Inf;
1539 for (octave_idx_type i = 0; i < nm; i++)
1540 {
1541 double tmp = fabs (data (i));
1542 if (tmp > dmax)
1543 dmax = tmp;
1544 if (tmp < dmin)
1545 dmin = tmp;
1546 }
1547 rcond = dmin / dmax;
1548 }
1549 else
1550 rcond = 1.;
1584 } 1551 }
1585 1552
1586 return retval; 1553 return retval;
1587 } 1554 }
1588 1555
1599 err = 0; 1566 err = 0;
1600 1567
1601 if (nr != b.rows ()) 1568 if (nr != b.rows ())
1602 (*current_liboctave_error_handler) 1569 (*current_liboctave_error_handler)
1603 ("matrix dimension mismatch solution of linear equations"); 1570 ("matrix dimension mismatch solution of linear equations");
1604 else if (nr == 0 || nc == 0 || b.cols () == 0) 1571
1572 if (nr == 0 || nc == 0 || b.cols () == 0)
1605 retval = SparseComplexMatrix (nc, b.cols ()); 1573 retval = SparseComplexMatrix (nc, b.cols ());
1606 else 1574 else
1607 { 1575 {
1608 // Print spparms("spumoni") info if requested 1576 // Print spparms("spumoni") info if requested
1609 int typ = mattype.type (); 1577 int typ = mattype.type ();
1610 mattype.info (); 1578 mattype.info ();
1611 1579
1612 if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) 1580 if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal)
1613 { 1581 (*current_liboctave_error_handler) ("incorrect matrix type");
1614 octave_idx_type b_nc = b.cols (); 1582
1615 octave_idx_type b_nz = b.nnz (); 1583 octave_idx_type b_nc = b.cols ();
1616 retval = SparseComplexMatrix (nc, b_nc, b_nz); 1584 octave_idx_type b_nz = b.nnz ();
1617 1585 retval = SparseComplexMatrix (nc, b_nc, b_nz);
1618 retval.xcidx (0) = 0; 1586
1619 octave_idx_type ii = 0; 1587 retval.xcidx (0) = 0;
1620 if (typ == MatrixType::Diagonal) 1588 octave_idx_type ii = 0;
1621 for (octave_idx_type j = 0; j < b.cols (); j++) 1589 if (typ == MatrixType::Diagonal)
1590 for (octave_idx_type j = 0; j < b.cols (); j++)
1591 {
1592 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1622 { 1593 {
1623 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 1594 if (b.ridx (i) >= nm)
1624 { 1595 break;
1625 if (b.ridx (i) >= nm) 1596 retval.xridx (ii) = b.ridx (i);
1626 break; 1597 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1627 retval.xridx (ii) = b.ridx (i);
1628 retval.xdata (ii++) = b.data (i) / data (b.ridx (i));
1629 }
1630 retval.xcidx (j+1) = ii;
1631 } 1598 }
1632 else 1599 retval.xcidx (j+1) = ii;
1633 for (octave_idx_type j = 0; j < b.cols (); j++) 1600 }
1634 {
1635 for (octave_idx_type l = 0; l < nc; l++)
1636 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1637 {
1638 bool found = false;
1639 octave_idx_type k;
1640 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1641 if (ridx (i) == b.ridx (k))
1642 {
1643 found = true;
1644 break;
1645 }
1646 if (found)
1647 {
1648 retval.xridx (ii) = l;
1649 retval.xdata (ii++) = b.data (k) / data (i);
1650 }
1651 }
1652 retval.xcidx (j+1) = ii;
1653 }
1654
1655 if (calc_cond)
1656 {
1657 double dmax = 0.;
1658 double dmin = octave_Inf;
1659 for (octave_idx_type i = 0; i < nm; i++)
1660 {
1661 double tmp = fabs (data (i));
1662 if (tmp > dmax)
1663 dmax = tmp;
1664 if (tmp < dmin)
1665 dmin = tmp;
1666 }
1667 rcond = dmin / dmax;
1668 }
1669 else
1670 rcond = 1.;
1671 }
1672 else 1601 else
1673 (*current_liboctave_error_handler) ("incorrect matrix type"); 1602 for (octave_idx_type j = 0; j < b.cols (); j++)
1603 {
1604 for (octave_idx_type l = 0; l < nc; l++)
1605 for (octave_idx_type i = cidx (l); i < cidx (l+1); i++)
1606 {
1607 bool found = false;
1608 octave_idx_type k;
1609 for (k = b.cidx (j); k < b.cidx (j+1); k++)
1610 if (ridx (i) == b.ridx (k))
1611 {
1612 found = true;
1613 break;
1614 }
1615 if (found)
1616 {
1617 retval.xridx (ii) = l;
1618 retval.xdata (ii++) = b.data (k) / data (i);
1619 }
1620 }
1621 retval.xcidx (j+1) = ii;
1622 }
1623
1624 if (calc_cond)
1625 {
1626 double dmax = 0.;
1627 double dmin = octave_Inf;
1628 for (octave_idx_type i = 0; i < nm; i++)
1629 {
1630 double tmp = fabs (data (i));
1631 if (tmp > dmax)
1632 dmax = tmp;
1633 if (tmp < dmin)
1634 dmin = tmp;
1635 }
1636 rcond = dmin / dmax;
1637 }
1638 else
1639 rcond = 1.;
1674 } 1640 }
1675 1641
1676 return retval; 1642 return retval;
1677 } 1643 }
1678 1644
1690 err = 0; 1656 err = 0;
1691 1657
1692 if (nr != b.rows ()) 1658 if (nr != b.rows ())
1693 (*current_liboctave_error_handler) 1659 (*current_liboctave_error_handler)
1694 ("matrix dimension mismatch solution of linear equations"); 1660 ("matrix dimension mismatch solution of linear equations");
1695 else if (nr == 0 || nc == 0 || b.cols () == 0) 1661
1662 if (nr == 0 || nc == 0 || b.cols () == 0)
1696 retval = Matrix (nc, b.cols (), 0.0); 1663 retval = Matrix (nc, b.cols (), 0.0);
1697 else 1664 else
1698 { 1665 {
1699 // Print spparms("spumoni") info if requested 1666 // Print spparms("spumoni") info if requested
1700 int typ = mattype.type (); 1667 int typ = mattype.type ();
1701 mattype.info (); 1668 mattype.info ();
1702 1669
1703 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper) 1670 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1704 { 1671 (*current_liboctave_error_handler) ("incorrect matrix type");
1705 double anorm = 0.; 1672
1706 double ainvnorm = 0.; 1673 double anorm = 0.;
1707 octave_idx_type b_nc = b.cols (); 1674 double ainvnorm = 0.;
1708 rcond = 1.; 1675 octave_idx_type b_nc = b.cols ();
1676 rcond = 1.;
1677
1678 if (calc_cond)
1679 {
1680 // Calculate the 1-norm of matrix for rcond calculation
1681 for (octave_idx_type j = 0; j < nc; j++)
1682 {
1683 double atmp = 0.;
1684 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1685 atmp += fabs (data (i));
1686 if (atmp > anorm)
1687 anorm = atmp;
1688 }
1689 }
1690
1691 if (typ == MatrixType::Permuted_Upper)
1692 {
1693 retval.resize (nc, b_nc);
1694 octave_idx_type *perm = mattype.triangular_perm ();
1695 OCTAVE_LOCAL_BUFFER (double, work, nm);
1696
1697 for (octave_idx_type j = 0; j < b_nc; j++)
1698 {
1699 for (octave_idx_type i = 0; i < nr; i++)
1700 work[i] = b(i,j);
1701 for (octave_idx_type i = nr; i < nc; i++)
1702 work[i] = 0.;
1703
1704 for (octave_idx_type k = nc-1; k >= 0; k--)
1705 {
1706 octave_idx_type kidx = perm[k];
1707
1708 if (work[k] != 0.)
1709 {
1710 if (ridx (cidx (kidx+1)-1) != k
1711 || data (cidx (kidx+1)-1) == 0.)
1712 {
1713 err = -2;
1714 goto triangular_error;
1715 }
1716
1717 double tmp = work[k] / data (cidx (kidx+1)-1);
1718 work[k] = tmp;
1719 for (octave_idx_type i = cidx (kidx);
1720 i < cidx (kidx+1)-1; i++)
1721 {
1722 octave_idx_type iidx = ridx (i);
1723 work[iidx] = work[iidx] - tmp * data (i);
1724 }
1725 }
1726 }
1727
1728 for (octave_idx_type i = 0; i < nc; i++)
1729 retval.xelem (perm[i], j) = work[i];
1730 }
1709 1731
1710 if (calc_cond) 1732 if (calc_cond)
1711 { 1733 {
1712 // Calculate the 1-norm of matrix for rcond calculation 1734 // Calculation of 1-norm of inv(*this)
1713 for (octave_idx_type j = 0; j < nc; j++) 1735 for (octave_idx_type i = 0; i < nm; i++)
1714 { 1736 work[i] = 0.;
1715 double atmp = 0.; 1737
1716 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 1738 for (octave_idx_type j = 0; j < nr; j++)
1717 atmp += fabs (data (i)); 1739 {
1718 if (atmp > anorm) 1740 work[j] = 1.;
1719 anorm = atmp; 1741
1720 } 1742 for (octave_idx_type k = j; k >= 0; k--)
1721 } 1743 {
1722 1744 octave_idx_type iidx = perm[k];
1723 if (typ == MatrixType::Permuted_Upper)
1724 {
1725 retval.resize (nc, b_nc);
1726 octave_idx_type *perm = mattype.triangular_perm ();
1727 OCTAVE_LOCAL_BUFFER (double, work, nm);
1728
1729 for (octave_idx_type j = 0; j < b_nc; j++)
1730 {
1731 for (octave_idx_type i = 0; i < nr; i++)
1732 work[i] = b(i,j);
1733 for (octave_idx_type i = nr; i < nc; i++)
1734 work[i] = 0.;
1735
1736 for (octave_idx_type k = nc-1; k >= 0; k--)
1737 {
1738 octave_idx_type kidx = perm[k];
1739 1745
1740 if (work[k] != 0.) 1746 if (work[k] != 0.)
1741 { 1747 {
1742 if (ridx (cidx (kidx+1)-1) != k 1748 double tmp = work[k] / data (cidx (iidx+1)-1);
1743 || data (cidx (kidx+1)-1) == 0.) 1749 work[k] = tmp;
1750 for (octave_idx_type i = cidx (iidx);
1751 i < cidx (iidx+1)-1; i++)
1744 { 1752 {
1745 err = -2; 1753 octave_idx_type idx2 = ridx (i);
1746 goto triangular_error; 1754 work[idx2] = work[idx2] - tmp * data (i);
1747 }
1748
1749 double tmp = work[k] / data (cidx (kidx+1)-1);
1750 work[k] = tmp;
1751 for (octave_idx_type i = cidx (kidx);
1752 i < cidx (kidx+1)-1; i++)
1753 {
1754 octave_idx_type iidx = ridx (i);
1755 work[iidx] = work[iidx] - tmp * data (i);
1756 } 1755 }
1757 } 1756 }
1758 } 1757 }
1759 1758 double atmp = 0;
1760 for (octave_idx_type i = 0; i < nc; i++) 1759 for (octave_idx_type i = 0; i < j+1; i++)
1761 retval.xelem (perm[i], j) = work[i]; 1760 {
1762 } 1761 atmp += fabs (work[i]);
1763 1762 work[i] = 0.;
1764 if (calc_cond) 1763 }
1765 { 1764 if (atmp > ainvnorm)
1766 // Calculation of 1-norm of inv(*this) 1765 ainvnorm = atmp;
1767 for (octave_idx_type i = 0; i < nm; i++) 1766 }
1768 work[i] = 0.; 1767 rcond = 1. / ainvnorm / anorm;
1769 1768 }
1770 for (octave_idx_type j = 0; j < nr; j++) 1769 }
1771 { 1770 else
1772 work[j] = 1.; 1771 {
1773 1772 OCTAVE_LOCAL_BUFFER (double, work, nm);
1774 for (octave_idx_type k = j; k >= 0; k--) 1773 retval.resize (nc, b_nc);
1774
1775 for (octave_idx_type j = 0; j < b_nc; j++)
1776 {
1777 for (octave_idx_type i = 0; i < nr; i++)
1778 work[i] = b(i,j);
1779 for (octave_idx_type i = nr; i < nc; i++)
1780 work[i] = 0.;
1781
1782 for (octave_idx_type k = nc-1; k >= 0; k--)
1783 {
1784 if (work[k] != 0.)
1785 {
1786 if (ridx (cidx (k+1)-1) != k
1787 || data (cidx (k+1)-1) == 0.)
1775 { 1788 {
1776 octave_idx_type iidx = perm[k]; 1789 err = -2;
1777 1790 goto triangular_error;
1778 if (work[k] != 0.)
1779 {
1780 double tmp = work[k] / data (cidx (iidx+1)-1);
1781 work[k] = tmp;
1782 for (octave_idx_type i = cidx (iidx);
1783 i < cidx (iidx+1)-1; i++)
1784 {
1785 octave_idx_type idx2 = ridx (i);
1786 work[idx2] = work[idx2] - tmp * data (i);
1787 }
1788 }
1789 } 1791 }
1790 double atmp = 0; 1792
1791 for (octave_idx_type i = 0; i < j+1; i++) 1793 double tmp = work[k] / data (cidx (k+1)-1);
1794 work[k] = tmp;
1795 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1792 { 1796 {
1793 atmp += fabs (work[i]); 1797 octave_idx_type iidx = ridx (i);
1794 work[i] = 0.; 1798 work[iidx] = work[iidx] - tmp * data (i);
1795 } 1799 }
1796 if (atmp > ainvnorm) 1800 }
1797 ainvnorm = atmp; 1801 }
1798 } 1802
1799 rcond = 1. / ainvnorm / anorm; 1803 for (octave_idx_type i = 0; i < nc; i++)
1800 } 1804 retval.xelem (i, j) = work[i];
1801 } 1805 }
1802 else 1806
1803 { 1807 if (calc_cond)
1804 OCTAVE_LOCAL_BUFFER (double, work, nm); 1808 {
1805 retval.resize (nc, b_nc); 1809 // Calculation of 1-norm of inv(*this)
1806 1810 for (octave_idx_type i = 0; i < nm; i++)
1807 for (octave_idx_type j = 0; j < b_nc; j++) 1811 work[i] = 0.;
1808 { 1812
1809 for (octave_idx_type i = 0; i < nr; i++) 1813 for (octave_idx_type j = 0; j < nr; j++)
1810 work[i] = b(i,j); 1814 {
1811 for (octave_idx_type i = nr; i < nc; i++) 1815 work[j] = 1.;
1812 work[i] = 0.; 1816
1813 1817 for (octave_idx_type k = j; k >= 0; k--)
1814 for (octave_idx_type k = nc-1; k >= 0; k--)
1815 { 1818 {
1816 if (work[k] != 0.) 1819 if (work[k] != 0.)
1817 { 1820 {
1818 if (ridx (cidx (k+1)-1) != k
1819 || data (cidx (k+1)-1) == 0.)
1820 {
1821 err = -2;
1822 goto triangular_error;
1823 }
1824
1825 double tmp = work[k] / data (cidx (k+1)-1); 1821 double tmp = work[k] / data (cidx (k+1)-1);
1826 work[k] = tmp; 1822 work[k] = tmp;
1827 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) 1823 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
1828 { 1824 {
1829 octave_idx_type iidx = ridx (i); 1825 octave_idx_type iidx = ridx (i);
1830 work[iidx] = work[iidx] - tmp * data (i); 1826 work[iidx] = work[iidx] - tmp * data (i);
1831 } 1827 }
1832 } 1828 }
1833 } 1829 }
1834 1830 double atmp = 0;
1835 for (octave_idx_type i = 0; i < nc; i++) 1831 for (octave_idx_type i = 0; i < j+1; i++)
1836 retval.xelem (i, j) = work[i]; 1832 {
1837 } 1833 atmp += fabs (work[i]);
1838 1834 work[i] = 0.;
1839 if (calc_cond) 1835 }
1840 { 1836 if (atmp > ainvnorm)
1841 // Calculation of 1-norm of inv(*this) 1837 ainvnorm = atmp;
1842 for (octave_idx_type i = 0; i < nm; i++) 1838 }
1843 work[i] = 0.; 1839 rcond = 1. / ainvnorm / anorm;
1844 1840 }
1845 for (octave_idx_type j = 0; j < nr; j++) 1841 }
1846 { 1842
1847 work[j] = 1.; 1843 triangular_error:
1848 1844 if (err != 0)
1849 for (octave_idx_type k = j; k >= 0; k--) 1845 {
1850 { 1846 if (sing_handler)
1851 if (work[k] != 0.) 1847 {
1852 { 1848 sing_handler (rcond);
1853 double tmp = work[k] / data (cidx (k+1)-1); 1849 mattype.mark_as_rectangular ();
1854 work[k] = tmp; 1850 }
1855 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) 1851 else
1856 { 1852 warn_singular_matrix (rcond);
1857 octave_idx_type iidx = ridx (i); 1853 }
1858 work[iidx] = work[iidx] - tmp * data (i); 1854
1859 } 1855 volatile double rcond_plus_one = rcond + 1.0;
1860 } 1856
1861 } 1857 if (rcond_plus_one == 1.0 || xisnan (rcond))
1862 double atmp = 0; 1858 {
1863 for (octave_idx_type i = 0; i < j+1; i++) 1859 err = -2;
1864 { 1860
1865 atmp += fabs (work[i]); 1861 if (sing_handler)
1866 work[i] = 0.; 1862 {
1867 } 1863 sing_handler (rcond);
1868 if (atmp > ainvnorm) 1864 mattype.mark_as_rectangular ();
1869 ainvnorm = atmp; 1865 }
1870 } 1866 else
1871 rcond = 1. / ainvnorm / anorm; 1867 warn_singular_matrix (rcond);
1872 } 1868 }
1873 }
1874
1875 triangular_error:
1876 if (err != 0)
1877 {
1878 if (sing_handler)
1879 {
1880 sing_handler (rcond);
1881 mattype.mark_as_rectangular ();
1882 }
1883 else
1884 warn_singular_matrix (rcond);
1885 }
1886
1887 volatile double rcond_plus_one = rcond + 1.0;
1888
1889 if (rcond_plus_one == 1.0 || xisnan (rcond))
1890 {
1891 err = -2;
1892
1893 if (sing_handler)
1894 {
1895 sing_handler (rcond);
1896 mattype.mark_as_rectangular ();
1897 }
1898 else
1899 warn_singular_matrix (rcond);
1900 }
1901 }
1902 else
1903 (*current_liboctave_error_handler) ("incorrect matrix type");
1904 } 1869 }
1905 1870
1906 return retval; 1871 return retval;
1907 } 1872 }
1908 1873
1920 err = 0; 1885 err = 0;
1921 1886
1922 if (nr != b.rows ()) 1887 if (nr != b.rows ())
1923 (*current_liboctave_error_handler) 1888 (*current_liboctave_error_handler)
1924 ("matrix dimension mismatch solution of linear equations"); 1889 ("matrix dimension mismatch solution of linear equations");
1925 else if (nr == 0 || nc == 0 || b.cols () == 0) 1890
1891 if (nr == 0 || nc == 0 || b.cols () == 0)
1926 retval = SparseMatrix (nc, b.cols ()); 1892 retval = SparseMatrix (nc, b.cols ());
1927 else 1893 else
1928 { 1894 {
1929 // Print spparms("spumoni") info if requested 1895 // Print spparms("spumoni") info if requested
1930 int typ = mattype.type (); 1896 int typ = mattype.type ();
1931 mattype.info (); 1897 mattype.info ();
1932 1898
1933 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper) 1899 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
1934 { 1900 (*current_liboctave_error_handler) ("incorrect matrix type");
1935 double anorm = 0.; 1901
1936 double ainvnorm = 0.; 1902 double anorm = 0.;
1937 rcond = 1.; 1903 double ainvnorm = 0.;
1904 rcond = 1.;
1905
1906 if (calc_cond)
1907 {
1908 // Calculate the 1-norm of matrix for rcond calculation
1909 for (octave_idx_type j = 0; j < nc; j++)
1910 {
1911 double atmp = 0.;
1912 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
1913 atmp += fabs (data (i));
1914 if (atmp > anorm)
1915 anorm = atmp;
1916 }
1917 }
1918
1919 octave_idx_type b_nc = b.cols ();
1920 octave_idx_type b_nz = b.nnz ();
1921 retval = SparseMatrix (nc, b_nc, b_nz);
1922 retval.xcidx (0) = 0;
1923 octave_idx_type ii = 0;
1924 octave_idx_type x_nz = b_nz;
1925
1926 if (typ == MatrixType::Permuted_Upper)
1927 {
1928 octave_idx_type *perm = mattype.triangular_perm ();
1929 OCTAVE_LOCAL_BUFFER (double, work, nm);
1930
1931 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1932 for (octave_idx_type i = 0; i < nc; i++)
1933 rperm[perm[i]] = i;
1934
1935 for (octave_idx_type j = 0; j < b_nc; j++)
1936 {
1937 for (octave_idx_type i = 0; i < nm; i++)
1938 work[i] = 0.;
1939 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1940 work[b.ridx (i)] = b.data (i);
1941
1942 for (octave_idx_type k = nc-1; k >= 0; k--)
1943 {
1944 octave_idx_type kidx = perm[k];
1945
1946 if (work[k] != 0.)
1947 {
1948 if (ridx (cidx (kidx+1)-1) != k
1949 || data (cidx (kidx+1)-1) == 0.)
1950 {
1951 err = -2;
1952 goto triangular_error;
1953 }
1954
1955 double tmp = work[k] / data (cidx (kidx+1)-1);
1956 work[k] = tmp;
1957 for (octave_idx_type i = cidx (kidx);
1958 i < cidx (kidx+1)-1; i++)
1959 {
1960 octave_idx_type iidx = ridx (i);
1961 work[iidx] = work[iidx] - tmp * data (i);
1962 }
1963 }
1964 }
1965
1966 // Count nonzeros in work vector and adjust space in
1967 // retval if needed
1968 octave_idx_type new_nnz = 0;
1969 for (octave_idx_type i = 0; i < nc; i++)
1970 if (work[i] != 0.)
1971 new_nnz++;
1972
1973 if (ii + new_nnz > x_nz)
1974 {
1975 // Resize the sparse matrix
1976 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
1977 retval.change_capacity (sz);
1978 x_nz = sz;
1979 }
1980
1981 for (octave_idx_type i = 0; i < nc; i++)
1982 if (work[rperm[i]] != 0.)
1983 {
1984 retval.xridx (ii) = i;
1985 retval.xdata (ii++) = work[rperm[i]];
1986 }
1987 retval.xcidx (j+1) = ii;
1988 }
1989
1990 retval.maybe_compress ();
1938 1991
1939 if (calc_cond) 1992 if (calc_cond)
1940 { 1993 {
1941 // Calculate the 1-norm of matrix for rcond calculation 1994 // Calculation of 1-norm of inv(*this)
1942 for (octave_idx_type j = 0; j < nc; j++) 1995 for (octave_idx_type i = 0; i < nm; i++)
1943 { 1996 work[i] = 0.;
1944 double atmp = 0.; 1997
1945 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 1998 for (octave_idx_type j = 0; j < nr; j++)
1946 atmp += fabs (data (i)); 1999 {
1947 if (atmp > anorm) 2000 work[j] = 1.;
1948 anorm = atmp; 2001
1949 } 2002 for (octave_idx_type k = j; k >= 0; k--)
1950 } 2003 {
1951 2004 octave_idx_type iidx = perm[k];
1952 octave_idx_type b_nc = b.cols ();
1953 octave_idx_type b_nz = b.nnz ();
1954 retval = SparseMatrix (nc, b_nc, b_nz);
1955 retval.xcidx (0) = 0;
1956 octave_idx_type ii = 0;
1957 octave_idx_type x_nz = b_nz;
1958
1959 if (typ == MatrixType::Permuted_Upper)
1960 {
1961 octave_idx_type *perm = mattype.triangular_perm ();
1962 OCTAVE_LOCAL_BUFFER (double, work, nm);
1963
1964 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
1965 for (octave_idx_type i = 0; i < nc; i++)
1966 rperm[perm[i]] = i;
1967
1968 for (octave_idx_type j = 0; j < b_nc; j++)
1969 {
1970 for (octave_idx_type i = 0; i < nm; i++)
1971 work[i] = 0.;
1972 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
1973 work[b.ridx (i)] = b.data (i);
1974
1975 for (octave_idx_type k = nc-1; k >= 0; k--)
1976 {
1977 octave_idx_type kidx = perm[k];
1978 2005
1979 if (work[k] != 0.) 2006 if (work[k] != 0.)
1980 { 2007 {
1981 if (ridx (cidx (kidx+1)-1) != k 2008 double tmp = work[k] / data (cidx (iidx+1)-1);
1982 || data (cidx (kidx+1)-1) == 0.) 2009 work[k] = tmp;
2010 for (octave_idx_type i = cidx (iidx);
2011 i < cidx (iidx+1)-1; i++)
1983 { 2012 {
1984 err = -2; 2013 octave_idx_type idx2 = ridx (i);
1985 goto triangular_error; 2014 work[idx2] = work[idx2] - tmp * data (i);
1986 } 2015 }
1987 2016 }
1988 double tmp = work[k] / data (cidx (kidx+1)-1); 2017 }
2018 double atmp = 0;
2019 for (octave_idx_type i = 0; i < j+1; i++)
2020 {
2021 atmp += fabs (work[i]);
2022 work[i] = 0.;
2023 }
2024 if (atmp > ainvnorm)
2025 ainvnorm = atmp;
2026 }
2027 rcond = 1. / ainvnorm / anorm;
2028 }
2029 }
2030 else
2031 {
2032 OCTAVE_LOCAL_BUFFER (double, work, nm);
2033
2034 for (octave_idx_type j = 0; j < b_nc; j++)
2035 {
2036 for (octave_idx_type i = 0; i < nm; i++)
2037 work[i] = 0.;
2038 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2039 work[b.ridx (i)] = b.data (i);
2040
2041 for (octave_idx_type k = nc-1; k >= 0; k--)
2042 {
2043 if (work[k] != 0.)
2044 {
2045 if (ridx (cidx (k+1)-1) != k
2046 || data (cidx (k+1)-1) == 0.)
2047 {
2048 err = -2;
2049 goto triangular_error;
2050 }
2051
2052 double tmp = work[k] / data (cidx (k+1)-1);
2053 work[k] = tmp;
2054 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2055 {
2056 octave_idx_type iidx = ridx (i);
2057 work[iidx] = work[iidx] - tmp * data (i);
2058 }
2059 }
2060 }
2061
2062 // Count nonzeros in work vector and adjust space in
2063 // retval if needed
2064 octave_idx_type new_nnz = 0;
2065 for (octave_idx_type i = 0; i < nc; i++)
2066 if (work[i] != 0.)
2067 new_nnz++;
2068
2069 if (ii + new_nnz > x_nz)
2070 {
2071 // Resize the sparse matrix
2072 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2073 retval.change_capacity (sz);
2074 x_nz = sz;
2075 }
2076
2077 for (octave_idx_type i = 0; i < nc; i++)
2078 if (work[i] != 0.)
2079 {
2080 retval.xridx (ii) = i;
2081 retval.xdata (ii++) = work[i];
2082 }
2083 retval.xcidx (j+1) = ii;
2084 }
2085
2086 retval.maybe_compress ();
2087
2088 if (calc_cond)
2089 {
2090 // Calculation of 1-norm of inv(*this)
2091 for (octave_idx_type i = 0; i < nm; i++)
2092 work[i] = 0.;
2093
2094 for (octave_idx_type j = 0; j < nr; j++)
2095 {
2096 work[j] = 1.;
2097
2098 for (octave_idx_type k = j; k >= 0; k--)
2099 {
2100 if (work[k] != 0.)
2101 {
2102 double tmp = work[k] / data (cidx (k+1)-1);
1989 work[k] = tmp; 2103 work[k] = tmp;
1990 for (octave_idx_type i = cidx (kidx); 2104 for (octave_idx_type i = cidx (k);
1991 i < cidx (kidx+1)-1; i++) 2105 i < cidx (k+1)-1; i++)
1992 { 2106 {
1993 octave_idx_type iidx = ridx (i); 2107 octave_idx_type iidx = ridx (i);
1994 work[iidx] = work[iidx] - tmp * data (i); 2108 work[iidx] = work[iidx] - tmp * data (i);
1995 } 2109 }
1996 } 2110 }
1997 } 2111 }
1998 2112 double atmp = 0;
1999 // Count nonzeros in work vector and adjust space in 2113 for (octave_idx_type i = 0; i < j+1; i++)
2000 // retval if needed 2114 {
2001 octave_idx_type new_nnz = 0; 2115 atmp += fabs (work[i]);
2002 for (octave_idx_type i = 0; i < nc; i++) 2116 work[i] = 0.;
2003 if (work[i] != 0.) 2117 }
2004 new_nnz++; 2118 if (atmp > ainvnorm)
2005 2119 ainvnorm = atmp;
2006 if (ii + new_nnz > x_nz) 2120 }
2007 { 2121 rcond = 1. / ainvnorm / anorm;
2008 // Resize the sparse matrix 2122 }
2009 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 2123 }
2010 retval.change_capacity (sz); 2124
2011 x_nz = sz; 2125 triangular_error:
2012 } 2126 if (err != 0)
2013 2127 {
2014 for (octave_idx_type i = 0; i < nc; i++) 2128 if (sing_handler)
2015 if (work[rperm[i]] != 0.) 2129 {
2016 { 2130 sing_handler (rcond);
2017 retval.xridx (ii) = i; 2131 mattype.mark_as_rectangular ();
2018 retval.xdata (ii++) = work[rperm[i]];
2019 }
2020 retval.xcidx (j+1) = ii;
2021 }
2022
2023 retval.maybe_compress ();
2024
2025 if (calc_cond)
2026 {
2027 // Calculation of 1-norm of inv(*this)
2028 for (octave_idx_type i = 0; i < nm; i++)
2029 work[i] = 0.;
2030
2031 for (octave_idx_type j = 0; j < nr; j++)
2032 {
2033 work[j] = 1.;
2034
2035 for (octave_idx_type k = j; k >= 0; k--)
2036 {
2037 octave_idx_type iidx = perm[k];
2038
2039 if (work[k] != 0.)
2040 {
2041 double tmp = work[k] / data (cidx (iidx+1)-1);
2042 work[k] = tmp;
2043 for (octave_idx_type i = cidx (iidx);
2044 i < cidx (iidx+1)-1; i++)
2045 {
2046 octave_idx_type idx2 = ridx (i);
2047 work[idx2] = work[idx2] - tmp * data (i);
2048 }
2049 }
2050 }
2051 double atmp = 0;
2052 for (octave_idx_type i = 0; i < j+1; i++)
2053 {
2054 atmp += fabs (work[i]);
2055 work[i] = 0.;
2056 }
2057 if (atmp > ainvnorm)
2058 ainvnorm = atmp;
2059 }
2060 rcond = 1. / ainvnorm / anorm;
2061 }
2062 } 2132 }
2063 else 2133 else
2064 { 2134 warn_singular_matrix (rcond);
2065 OCTAVE_LOCAL_BUFFER (double, work, nm); 2135 }
2066 2136
2067 for (octave_idx_type j = 0; j < b_nc; j++) 2137 volatile double rcond_plus_one = rcond + 1.0;
2068 { 2138
2069 for (octave_idx_type i = 0; i < nm; i++) 2139 if (rcond_plus_one == 1.0 || xisnan (rcond))
2070 work[i] = 0.; 2140 {
2071 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 2141 err = -2;
2072 work[b.ridx (i)] = b.data (i); 2142
2073 2143 if (sing_handler)
2074 for (octave_idx_type k = nc-1; k >= 0; k--) 2144 {
2075 { 2145 sing_handler (rcond);
2076 if (work[k] != 0.) 2146 mattype.mark_as_rectangular ();
2077 { 2147 }
2078 if (ridx (cidx (k+1)-1) != k 2148 else
2079 || data (cidx (k+1)-1) == 0.) 2149 warn_singular_matrix (rcond);
2080 { 2150 }
2081 err = -2;
2082 goto triangular_error;
2083 }
2084
2085 double tmp = work[k] / data (cidx (k+1)-1);
2086 work[k] = tmp;
2087 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2088 {
2089 octave_idx_type iidx = ridx (i);
2090 work[iidx] = work[iidx] - tmp * data (i);
2091 }
2092 }
2093 }
2094
2095 // Count nonzeros in work vector and adjust space in
2096 // retval if needed
2097 octave_idx_type new_nnz = 0;
2098 for (octave_idx_type i = 0; i < nc; i++)
2099 if (work[i] != 0.)
2100 new_nnz++;
2101
2102 if (ii + new_nnz > x_nz)
2103 {
2104 // Resize the sparse matrix
2105 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2106 retval.change_capacity (sz);
2107 x_nz = sz;
2108 }
2109
2110 for (octave_idx_type i = 0; i < nc; i++)
2111 if (work[i] != 0.)
2112 {
2113 retval.xridx (ii) = i;
2114 retval.xdata (ii++) = work[i];
2115 }
2116 retval.xcidx (j+1) = ii;
2117 }
2118
2119 retval.maybe_compress ();
2120
2121 if (calc_cond)
2122 {
2123 // Calculation of 1-norm of inv(*this)
2124 for (octave_idx_type i = 0; i < nm; i++)
2125 work[i] = 0.;
2126
2127 for (octave_idx_type j = 0; j < nr; j++)
2128 {
2129 work[j] = 1.;
2130
2131 for (octave_idx_type k = j; k >= 0; k--)
2132 {
2133 if (work[k] != 0.)
2134 {
2135 double tmp = work[k] / data (cidx (k+1)-1);
2136 work[k] = tmp;
2137 for (octave_idx_type i = cidx (k);
2138 i < cidx (k+1)-1; i++)
2139 {
2140 octave_idx_type iidx = ridx (i);
2141 work[iidx] = work[iidx] - tmp * data (i);
2142 }
2143 }
2144 }
2145 double atmp = 0;
2146 for (octave_idx_type i = 0; i < j+1; i++)
2147 {
2148 atmp += fabs (work[i]);
2149 work[i] = 0.;
2150 }
2151 if (atmp > ainvnorm)
2152 ainvnorm = atmp;
2153 }
2154 rcond = 1. / ainvnorm / anorm;
2155 }
2156 }
2157
2158 triangular_error:
2159 if (err != 0)
2160 {
2161 if (sing_handler)
2162 {
2163 sing_handler (rcond);
2164 mattype.mark_as_rectangular ();
2165 }
2166 else
2167 warn_singular_matrix (rcond);
2168 }
2169
2170 volatile double rcond_plus_one = rcond + 1.0;
2171
2172 if (rcond_plus_one == 1.0 || xisnan (rcond))
2173 {
2174 err = -2;
2175
2176 if (sing_handler)
2177 {
2178 sing_handler (rcond);
2179 mattype.mark_as_rectangular ();
2180 }
2181 else
2182 warn_singular_matrix (rcond);
2183 }
2184 }
2185 else
2186 (*current_liboctave_error_handler) ("incorrect matrix type");
2187 } 2151 }
2188 return retval; 2152 return retval;
2189 } 2153 }
2190 2154
2191 ComplexMatrix 2155 ComplexMatrix
2202 err = 0; 2166 err = 0;
2203 2167
2204 if (nr != b.rows ()) 2168 if (nr != b.rows ())
2205 (*current_liboctave_error_handler) 2169 (*current_liboctave_error_handler)
2206 ("matrix dimension mismatch solution of linear equations"); 2170 ("matrix dimension mismatch solution of linear equations");
2207 else if (nr == 0 || nc == 0 || b.cols () == 0) 2171
2172 if (nr == 0 || nc == 0 || b.cols () == 0)
2208 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 2173 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
2209 else 2174 else
2210 { 2175 {
2211 // Print spparms("spumoni") info if requested 2176 // Print spparms("spumoni") info if requested
2212 int typ = mattype.type (); 2177 int typ = mattype.type ();
2213 mattype.info (); 2178 mattype.info ();
2214 2179
2215 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper) 2180 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2216 { 2181 (*current_liboctave_error_handler) ("incorrect matrix type");
2217 double anorm = 0.; 2182
2218 double ainvnorm = 0.; 2183 double anorm = 0.;
2219 octave_idx_type b_nc = b.cols (); 2184 double ainvnorm = 0.;
2220 rcond = 1.; 2185 octave_idx_type b_nc = b.cols ();
2186 rcond = 1.;
2187
2188 if (calc_cond)
2189 {
2190 // Calculate the 1-norm of matrix for rcond calculation
2191 for (octave_idx_type j = 0; j < nc; j++)
2192 {
2193 double atmp = 0.;
2194 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2195 atmp += fabs (data (i));
2196 if (atmp > anorm)
2197 anorm = atmp;
2198 }
2199 }
2200
2201 if (typ == MatrixType::Permuted_Upper)
2202 {
2203 retval.resize (nc, b_nc);
2204 octave_idx_type *perm = mattype.triangular_perm ();
2205 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2206
2207 for (octave_idx_type j = 0; j < b_nc; j++)
2208 {
2209 for (octave_idx_type i = 0; i < nr; i++)
2210 cwork[i] = b(i,j);
2211 for (octave_idx_type i = nr; i < nc; i++)
2212 cwork[i] = 0.;
2213
2214 for (octave_idx_type k = nc-1; k >= 0; k--)
2215 {
2216 octave_idx_type kidx = perm[k];
2217
2218 if (cwork[k] != 0.)
2219 {
2220 if (ridx (cidx (kidx+1)-1) != k
2221 || data (cidx (kidx+1)-1) == 0.)
2222 {
2223 err = -2;
2224 goto triangular_error;
2225 }
2226
2227 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2228 cwork[k] = tmp;
2229 for (octave_idx_type i = cidx (kidx);
2230 i < cidx (kidx+1)-1; i++)
2231 {
2232 octave_idx_type iidx = ridx (i);
2233 cwork[iidx] = cwork[iidx] - tmp * data (i);
2234 }
2235 }
2236 }
2237
2238 for (octave_idx_type i = 0; i < nc; i++)
2239 retval.xelem (perm[i], j) = cwork[i];
2240 }
2221 2241
2222 if (calc_cond) 2242 if (calc_cond)
2223 { 2243 {
2224 // Calculate the 1-norm of matrix for rcond calculation 2244 // Calculation of 1-norm of inv(*this)
2225 for (octave_idx_type j = 0; j < nc; j++) 2245 OCTAVE_LOCAL_BUFFER (double, work, nm);
2226 { 2246 for (octave_idx_type i = 0; i < nm; i++)
2227 double atmp = 0.; 2247 work[i] = 0.;
2228 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 2248
2229 atmp += fabs (data (i)); 2249 for (octave_idx_type j = 0; j < nr; j++)
2230 if (atmp > anorm) 2250 {
2231 anorm = atmp; 2251 work[j] = 1.;
2232 } 2252
2233 } 2253 for (octave_idx_type k = j; k >= 0; k--)
2234 2254 {
2235 if (typ == MatrixType::Permuted_Upper) 2255 octave_idx_type iidx = perm[k];
2236 { 2256
2237 retval.resize (nc, b_nc); 2257 if (work[k] != 0.)
2238 octave_idx_type *perm = mattype.triangular_perm ();
2239 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2240
2241 for (octave_idx_type j = 0; j < b_nc; j++)
2242 {
2243 for (octave_idx_type i = 0; i < nr; i++)
2244 cwork[i] = b(i,j);
2245 for (octave_idx_type i = nr; i < nc; i++)
2246 cwork[i] = 0.;
2247
2248 for (octave_idx_type k = nc-1; k >= 0; k--)
2249 {
2250 octave_idx_type kidx = perm[k];
2251
2252 if (cwork[k] != 0.)
2253 { 2258 {
2254 if (ridx (cidx (kidx+1)-1) != k 2259 double tmp = work[k] / data (cidx (iidx+1)-1);
2255 || data (cidx (kidx+1)-1) == 0.) 2260 work[k] = tmp;
2261 for (octave_idx_type i = cidx (iidx);
2262 i < cidx (iidx+1)-1; i++)
2256 { 2263 {
2257 err = -2; 2264 octave_idx_type idx2 = ridx (i);
2258 goto triangular_error; 2265 work[idx2] = work[idx2] - tmp * data (i);
2259 } 2266 }
2260 2267 }
2261 Complex tmp = cwork[k] / data (cidx (kidx+1)-1); 2268 }
2262 cwork[k] = tmp; 2269 double atmp = 0;
2263 for (octave_idx_type i = cidx (kidx); 2270 for (octave_idx_type i = 0; i < j+1; i++)
2264 i < cidx (kidx+1)-1; i++) 2271 {
2272 atmp += fabs (work[i]);
2273 work[i] = 0.;
2274 }
2275 if (atmp > ainvnorm)
2276 ainvnorm = atmp;
2277 }
2278 rcond = 1. / ainvnorm / anorm;
2279 }
2280 }
2281 else
2282 {
2283 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2284 retval.resize (nc, b_nc);
2285
2286 for (octave_idx_type j = 0; j < b_nc; j++)
2287 {
2288 for (octave_idx_type i = 0; i < nr; i++)
2289 cwork[i] = b(i,j);
2290 for (octave_idx_type i = nr; i < nc; i++)
2291 cwork[i] = 0.;
2292
2293 for (octave_idx_type k = nc-1; k >= 0; k--)
2294 {
2295 if (cwork[k] != 0.)
2296 {
2297 if (ridx (cidx (k+1)-1) != k
2298 || data (cidx (k+1)-1) == 0.)
2299 {
2300 err = -2;
2301 goto triangular_error;
2302 }
2303
2304 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2305 cwork[k] = tmp;
2306 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2307 {
2308 octave_idx_type iidx = ridx (i);
2309 cwork[iidx] = cwork[iidx] - tmp * data (i);
2310 }
2311 }
2312 }
2313
2314 for (octave_idx_type i = 0; i < nc; i++)
2315 retval.xelem (i, j) = cwork[i];
2316 }
2317
2318 if (calc_cond)
2319 {
2320 // Calculation of 1-norm of inv(*this)
2321 OCTAVE_LOCAL_BUFFER (double, work, nm);
2322 for (octave_idx_type i = 0; i < nm; i++)
2323 work[i] = 0.;
2324
2325 for (octave_idx_type j = 0; j < nr; j++)
2326 {
2327 work[j] = 1.;
2328
2329 for (octave_idx_type k = j; k >= 0; k--)
2330 {
2331 if (work[k] != 0.)
2332 {
2333 double tmp = work[k] / data (cidx (k+1)-1);
2334 work[k] = tmp;
2335 for (octave_idx_type i = cidx (k);
2336 i < cidx (k+1)-1; i++)
2265 { 2337 {
2266 octave_idx_type iidx = ridx (i); 2338 octave_idx_type iidx = ridx (i);
2267 cwork[iidx] = cwork[iidx] - tmp * data (i); 2339 work[iidx] = work[iidx] - tmp * data (i);
2268 } 2340 }
2269 } 2341 }
2270 } 2342 }
2271 2343 double atmp = 0;
2272 for (octave_idx_type i = 0; i < nc; i++) 2344 for (octave_idx_type i = 0; i < j+1; i++)
2273 retval.xelem (perm[i], j) = cwork[i]; 2345 {
2274 } 2346 atmp += fabs (work[i]);
2275 2347 work[i] = 0.;
2276 if (calc_cond) 2348 }
2277 { 2349 if (atmp > ainvnorm)
2278 // Calculation of 1-norm of inv(*this) 2350 ainvnorm = atmp;
2279 OCTAVE_LOCAL_BUFFER (double, work, nm); 2351 }
2280 for (octave_idx_type i = 0; i < nm; i++) 2352 rcond = 1. / ainvnorm / anorm;
2281 work[i] = 0.; 2353 }
2282 2354 }
2283 for (octave_idx_type j = 0; j < nr; j++) 2355
2284 { 2356 triangular_error:
2285 work[j] = 1.; 2357 if (err != 0)
2286 2358 {
2287 for (octave_idx_type k = j; k >= 0; k--) 2359 if (sing_handler)
2288 { 2360 {
2289 octave_idx_type iidx = perm[k]; 2361 sing_handler (rcond);
2290 2362 mattype.mark_as_rectangular ();
2291 if (work[k] != 0.)
2292 {
2293 double tmp = work[k] / data (cidx (iidx+1)-1);
2294 work[k] = tmp;
2295 for (octave_idx_type i = cidx (iidx);
2296 i < cidx (iidx+1)-1; i++)
2297 {
2298 octave_idx_type idx2 = ridx (i);
2299 work[idx2] = work[idx2] - tmp * data (i);
2300 }
2301 }
2302 }
2303 double atmp = 0;
2304 for (octave_idx_type i = 0; i < j+1; i++)
2305 {
2306 atmp += fabs (work[i]);
2307 work[i] = 0.;
2308 }
2309 if (atmp > ainvnorm)
2310 ainvnorm = atmp;
2311 }
2312 rcond = 1. / ainvnorm / anorm;
2313 }
2314 } 2363 }
2315 else 2364 else
2316 { 2365 warn_singular_matrix (rcond);
2317 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 2366 }
2318 retval.resize (nc, b_nc); 2367
2319 2368 volatile double rcond_plus_one = rcond + 1.0;
2320 for (octave_idx_type j = 0; j < b_nc; j++) 2369
2321 { 2370 if (rcond_plus_one == 1.0 || xisnan (rcond))
2322 for (octave_idx_type i = 0; i < nr; i++) 2371 {
2323 cwork[i] = b(i,j); 2372 err = -2;
2324 for (octave_idx_type i = nr; i < nc; i++) 2373
2325 cwork[i] = 0.; 2374 if (sing_handler)
2326 2375 {
2327 for (octave_idx_type k = nc-1; k >= 0; k--) 2376 sing_handler (rcond);
2328 { 2377 mattype.mark_as_rectangular ();
2329 if (cwork[k] != 0.) 2378 }
2330 { 2379 else
2331 if (ridx (cidx (k+1)-1) != k 2380 warn_singular_matrix (rcond);
2332 || data (cidx (k+1)-1) == 0.) 2381 }
2333 {
2334 err = -2;
2335 goto triangular_error;
2336 }
2337
2338 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2339 cwork[k] = tmp;
2340 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2341 {
2342 octave_idx_type iidx = ridx (i);
2343 cwork[iidx] = cwork[iidx] - tmp * data (i);
2344 }
2345 }
2346 }
2347
2348 for (octave_idx_type i = 0; i < nc; i++)
2349 retval.xelem (i, j) = cwork[i];
2350 }
2351
2352 if (calc_cond)
2353 {
2354 // Calculation of 1-norm of inv(*this)
2355 OCTAVE_LOCAL_BUFFER (double, work, nm);
2356 for (octave_idx_type i = 0; i < nm; i++)
2357 work[i] = 0.;
2358
2359 for (octave_idx_type j = 0; j < nr; j++)
2360 {
2361 work[j] = 1.;
2362
2363 for (octave_idx_type k = j; k >= 0; k--)
2364 {
2365 if (work[k] != 0.)
2366 {
2367 double tmp = work[k] / data (cidx (k+1)-1);
2368 work[k] = tmp;
2369 for (octave_idx_type i = cidx (k);
2370 i < cidx (k+1)-1; i++)
2371 {
2372 octave_idx_type iidx = ridx (i);
2373 work[iidx] = work[iidx] - tmp * data (i);
2374 }
2375 }
2376 }
2377 double atmp = 0;
2378 for (octave_idx_type i = 0; i < j+1; i++)
2379 {
2380 atmp += fabs (work[i]);
2381 work[i] = 0.;
2382 }
2383 if (atmp > ainvnorm)
2384 ainvnorm = atmp;
2385 }
2386 rcond = 1. / ainvnorm / anorm;
2387 }
2388 }
2389
2390 triangular_error:
2391 if (err != 0)
2392 {
2393 if (sing_handler)
2394 {
2395 sing_handler (rcond);
2396 mattype.mark_as_rectangular ();
2397 }
2398 else
2399 warn_singular_matrix (rcond);
2400 }
2401
2402 volatile double rcond_plus_one = rcond + 1.0;
2403
2404 if (rcond_plus_one == 1.0 || xisnan (rcond))
2405 {
2406 err = -2;
2407
2408 if (sing_handler)
2409 {
2410 sing_handler (rcond);
2411 mattype.mark_as_rectangular ();
2412 }
2413 else
2414 warn_singular_matrix (rcond);
2415 }
2416 }
2417 else
2418 (*current_liboctave_error_handler) ("incorrect matrix type");
2419 } 2382 }
2420 2383
2421 return retval; 2384 return retval;
2422 } 2385 }
2423 2386
2435 err = 0; 2398 err = 0;
2436 2399
2437 if (nr != b.rows ()) 2400 if (nr != b.rows ())
2438 (*current_liboctave_error_handler) 2401 (*current_liboctave_error_handler)
2439 ("matrix dimension mismatch solution of linear equations"); 2402 ("matrix dimension mismatch solution of linear equations");
2440 else if (nr == 0 || nc == 0 || b.cols () == 0) 2403
2404 if (nr == 0 || nc == 0 || b.cols () == 0)
2441 retval = SparseComplexMatrix (nc, b.cols ()); 2405 retval = SparseComplexMatrix (nc, b.cols ());
2442 else 2406 else
2443 { 2407 {
2444 // Print spparms("spumoni") info if requested 2408 // Print spparms("spumoni") info if requested
2445 int typ = mattype.type (); 2409 int typ = mattype.type ();
2446 mattype.info (); 2410 mattype.info ();
2447 2411
2448 if (typ == MatrixType::Permuted_Upper || typ == MatrixType::Upper) 2412 if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper)
2449 { 2413 (*current_liboctave_error_handler) ("incorrect matrix type");
2450 double anorm = 0.; 2414
2451 double ainvnorm = 0.; 2415 double anorm = 0.;
2452 rcond = 1.; 2416 double ainvnorm = 0.;
2417 rcond = 1.;
2418
2419 if (calc_cond)
2420 {
2421 // Calculate the 1-norm of matrix for rcond calculation
2422 for (octave_idx_type j = 0; j < nc; j++)
2423 {
2424 double atmp = 0.;
2425 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2426 atmp += fabs (data (i));
2427 if (atmp > anorm)
2428 anorm = atmp;
2429 }
2430 }
2431
2432 octave_idx_type b_nc = b.cols ();
2433 octave_idx_type b_nz = b.nnz ();
2434 retval = SparseComplexMatrix (nc, b_nc, b_nz);
2435 retval.xcidx (0) = 0;
2436 octave_idx_type ii = 0;
2437 octave_idx_type x_nz = b_nz;
2438
2439 if (typ == MatrixType::Permuted_Upper)
2440 {
2441 octave_idx_type *perm = mattype.triangular_perm ();
2442 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2443
2444 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc);
2445 for (octave_idx_type i = 0; i < nc; i++)
2446 rperm[perm[i]] = i;
2447
2448 for (octave_idx_type j = 0; j < b_nc; j++)
2449 {
2450 for (octave_idx_type i = 0; i < nm; i++)
2451 cwork[i] = 0.;
2452 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2453 cwork[b.ridx (i)] = b.data (i);
2454
2455 for (octave_idx_type k = nc-1; k >= 0; k--)
2456 {
2457 octave_idx_type kidx = perm[k];
2458
2459 if (cwork[k] != 0.)
2460 {
2461 if (ridx (cidx (kidx+1)-1) != k
2462 || data (cidx (kidx+1)-1) == 0.)
2463 {
2464 err = -2;
2465 goto triangular_error;
2466 }
2467
2468 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2469 cwork[k] = tmp;
2470 for (octave_idx_type i = cidx (kidx);
2471 i < cidx (kidx+1)-1; i++)
2472 {
2473 octave_idx_type iidx = ridx (i);
2474 cwork[iidx] = cwork[iidx] - tmp * data (i);
2475 }
2476 }
2477 }
2478
2479 // Count nonzeros in work vector and adjust space in
2480 // retval if needed
2481 octave_idx_type new_nnz = 0;
2482 for (octave_idx_type i = 0; i < nc; i++)
2483 if (cwork[i] != 0.)
2484 new_nnz++;
2485
2486 if (ii + new_nnz > x_nz)
2487 {
2488 // Resize the sparse matrix
2489 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2490 retval.change_capacity (sz);
2491 x_nz = sz;
2492 }
2493
2494 for (octave_idx_type i = 0; i < nc; i++)
2495 if (cwork[rperm[i]] != 0.)
2496 {
2497 retval.xridx (ii) = i;
2498 retval.xdata (ii++) = cwork[rperm[i]];
2499 }
2500 retval.xcidx (j+1) = ii;
2501 }
2502
2503 retval.maybe_compress ();
2453 2504
2454 if (calc_cond) 2505 if (calc_cond)
2455 { 2506 {
2456 // Calculate the 1-norm of matrix for rcond calculation 2507 // Calculation of 1-norm of inv(*this)
2457 for (octave_idx_type j = 0; j < nc; j++) 2508 OCTAVE_LOCAL_BUFFER (double, work, nm);
2458 { 2509 for (octave_idx_type i = 0; i < nm; i++)
2459 double atmp = 0.; 2510 work[i] = 0.;
2460 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 2511
2461 atmp += fabs (data (i)); 2512 for (octave_idx_type j = 0; j < nr; j++)
2462 if (atmp > anorm) 2513 {
2463 anorm = atmp; 2514 work[j] = 1.;
2464 } 2515
2465 } 2516 for (octave_idx_type k = j; k >= 0; k--)
2466 2517 {
2467 octave_idx_type b_nc = b.cols (); 2518 octave_idx_type iidx = perm[k];
2468 octave_idx_type b_nz = b.nnz (); 2519
2469 retval = SparseComplexMatrix (nc, b_nc, b_nz); 2520 if (work[k] != 0.)
2470 retval.xcidx (0) = 0; 2521 {
2471 octave_idx_type ii = 0; 2522 double tmp = work[k] / data (cidx (iidx+1)-1);
2472 octave_idx_type x_nz = b_nz; 2523 work[k] = tmp;
2473 2524 for (octave_idx_type i = cidx (iidx);
2474 if (typ == MatrixType::Permuted_Upper) 2525 i < cidx (iidx+1)-1; i++)
2475 { 2526 {
2476 octave_idx_type *perm = mattype.triangular_perm (); 2527 octave_idx_type idx2 = ridx (i);
2477 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 2528 work[idx2] = work[idx2] - tmp * data (i);
2478 2529 }
2479 OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); 2530 }
2531 }
2532 double atmp = 0;
2533 for (octave_idx_type i = 0; i < j+1; i++)
2534 {
2535 atmp += fabs (work[i]);
2536 work[i] = 0.;
2537 }
2538 if (atmp > ainvnorm)
2539 ainvnorm = atmp;
2540 }
2541 rcond = 1. / ainvnorm / anorm;
2542 }
2543 }
2544 else
2545 {
2546 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
2547
2548 for (octave_idx_type j = 0; j < b_nc; j++)
2549 {
2550 for (octave_idx_type i = 0; i < nm; i++)
2551 cwork[i] = 0.;
2552 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2553 cwork[b.ridx (i)] = b.data (i);
2554
2555 for (octave_idx_type k = nc-1; k >= 0; k--)
2556 {
2557 if (cwork[k] != 0.)
2558 {
2559 if (ridx (cidx (k+1)-1) != k
2560 || data (cidx (k+1)-1) == 0.)
2561 {
2562 err = -2;
2563 goto triangular_error;
2564 }
2565
2566 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2567 cwork[k] = tmp;
2568 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2569 {
2570 octave_idx_type iidx = ridx (i);
2571 cwork[iidx] = cwork[iidx] - tmp * data (i);
2572 }
2573 }
2574 }
2575
2576 // Count nonzeros in work vector and adjust space in
2577 // retval if needed
2578 octave_idx_type new_nnz = 0;
2480 for (octave_idx_type i = 0; i < nc; i++) 2579 for (octave_idx_type i = 0; i < nc; i++)
2481 rperm[perm[i]] = i; 2580 if (cwork[i] != 0.)
2482 2581 new_nnz++;
2483 for (octave_idx_type j = 0; j < b_nc; j++) 2582
2484 { 2583 if (ii + new_nnz > x_nz)
2485 for (octave_idx_type i = 0; i < nm; i++) 2584 {
2486 cwork[i] = 0.; 2585 // Resize the sparse matrix
2487 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 2586 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2488 cwork[b.ridx (i)] = b.data (i); 2587 retval.change_capacity (sz);
2489 2588 x_nz = sz;
2490 for (octave_idx_type k = nc-1; k >= 0; k--) 2589 }
2491 { 2590
2492 octave_idx_type kidx = perm[k]; 2591 for (octave_idx_type i = 0; i < nc; i++)
2493 2592 if (cwork[i] != 0.)
2494 if (cwork[k] != 0.) 2593 {
2594 retval.xridx (ii) = i;
2595 retval.xdata (ii++) = cwork[i];
2596 }
2597 retval.xcidx (j+1) = ii;
2598 }
2599
2600 retval.maybe_compress ();
2601
2602 if (calc_cond)
2603 {
2604 // Calculation of 1-norm of inv(*this)
2605 OCTAVE_LOCAL_BUFFER (double, work, nm);
2606 for (octave_idx_type i = 0; i < nm; i++)
2607 work[i] = 0.;
2608
2609 for (octave_idx_type j = 0; j < nr; j++)
2610 {
2611 work[j] = 1.;
2612
2613 for (octave_idx_type k = j; k >= 0; k--)
2614 {
2615 if (work[k] != 0.)
2495 { 2616 {
2496 if (ridx (cidx (kidx+1)-1) != k 2617 double tmp = work[k] / data (cidx (k+1)-1);
2497 || data (cidx (kidx+1)-1) == 0.) 2618 work[k] = tmp;
2498 { 2619 for (octave_idx_type i = cidx (k);
2499 err = -2; 2620 i < cidx (k+1)-1; i++)
2500 goto triangular_error;
2501 }
2502
2503 Complex tmp = cwork[k] / data (cidx (kidx+1)-1);
2504 cwork[k] = tmp;
2505 for (octave_idx_type i = cidx (kidx);
2506 i < cidx (kidx+1)-1; i++)
2507 { 2621 {
2508 octave_idx_type iidx = ridx (i); 2622 octave_idx_type iidx = ridx (i);
2509 cwork[iidx] = cwork[iidx] - tmp * data (i); 2623 work[iidx] = work[iidx] - tmp * data (i);
2510 } 2624 }
2511 } 2625 }
2512 } 2626 }
2513 2627 double atmp = 0;
2514 // Count nonzeros in work vector and adjust space in 2628 for (octave_idx_type i = 0; i < j+1; i++)
2515 // retval if needed 2629 {
2516 octave_idx_type new_nnz = 0; 2630 atmp += fabs (work[i]);
2517 for (octave_idx_type i = 0; i < nc; i++) 2631 work[i] = 0.;
2518 if (cwork[i] != 0.) 2632 }
2519 new_nnz++; 2633 if (atmp > ainvnorm)
2520 2634 ainvnorm = atmp;
2521 if (ii + new_nnz > x_nz) 2635 }
2522 { 2636 rcond = 1. / ainvnorm / anorm;
2523 // Resize the sparse matrix 2637 }
2524 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 2638 }
2525 retval.change_capacity (sz); 2639
2526 x_nz = sz; 2640 triangular_error:
2527 } 2641 if (err != 0)
2528 2642 {
2529 for (octave_idx_type i = 0; i < nc; i++) 2643 if (sing_handler)
2530 if (cwork[rperm[i]] != 0.) 2644 {
2531 { 2645 sing_handler (rcond);
2532 retval.xridx (ii) = i; 2646 mattype.mark_as_rectangular ();
2533 retval.xdata (ii++) = cwork[rperm[i]];
2534 }
2535 retval.xcidx (j+1) = ii;
2536 }
2537
2538 retval.maybe_compress ();
2539
2540 if (calc_cond)
2541 {
2542 // Calculation of 1-norm of inv(*this)
2543 OCTAVE_LOCAL_BUFFER (double, work, nm);
2544 for (octave_idx_type i = 0; i < nm; i++)
2545 work[i] = 0.;
2546
2547 for (octave_idx_type j = 0; j < nr; j++)
2548 {
2549 work[j] = 1.;
2550
2551 for (octave_idx_type k = j; k >= 0; k--)
2552 {
2553 octave_idx_type iidx = perm[k];
2554
2555 if (work[k] != 0.)
2556 {
2557 double tmp = work[k] / data (cidx (iidx+1)-1);
2558 work[k] = tmp;
2559 for (octave_idx_type i = cidx (iidx);
2560 i < cidx (iidx+1)-1; i++)
2561 {
2562 octave_idx_type idx2 = ridx (i);
2563 work[idx2] = work[idx2] - tmp * data (i);
2564 }
2565 }
2566 }
2567 double atmp = 0;
2568 for (octave_idx_type i = 0; i < j+1; i++)
2569 {
2570 atmp += fabs (work[i]);
2571 work[i] = 0.;
2572 }
2573 if (atmp > ainvnorm)
2574 ainvnorm = atmp;
2575 }
2576 rcond = 1. / ainvnorm / anorm;
2577 }
2578 } 2647 }
2579 else 2648 else
2580 { 2649 warn_singular_matrix (rcond);
2581 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 2650 }
2582 2651
2583 for (octave_idx_type j = 0; j < b_nc; j++) 2652 volatile double rcond_plus_one = rcond + 1.0;
2584 { 2653
2585 for (octave_idx_type i = 0; i < nm; i++) 2654 if (rcond_plus_one == 1.0 || xisnan (rcond))
2586 cwork[i] = 0.; 2655 {
2587 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 2656 err = -2;
2588 cwork[b.ridx (i)] = b.data (i); 2657
2589 2658 if (sing_handler)
2590 for (octave_idx_type k = nc-1; k >= 0; k--) 2659 {
2591 { 2660 sing_handler (rcond);
2592 if (cwork[k] != 0.) 2661 mattype.mark_as_rectangular ();
2593 { 2662 }
2594 if (ridx (cidx (k+1)-1) != k 2663 else
2595 || data (cidx (k+1)-1) == 0.) 2664 warn_singular_matrix (rcond);
2596 { 2665 }
2597 err = -2;
2598 goto triangular_error;
2599 }
2600
2601 Complex tmp = cwork[k] / data (cidx (k+1)-1);
2602 cwork[k] = tmp;
2603 for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++)
2604 {
2605 octave_idx_type iidx = ridx (i);
2606 cwork[iidx] = cwork[iidx] - tmp * data (i);
2607 }
2608 }
2609 }
2610
2611 // Count nonzeros in work vector and adjust space in
2612 // retval if needed
2613 octave_idx_type new_nnz = 0;
2614 for (octave_idx_type i = 0; i < nc; i++)
2615 if (cwork[i] != 0.)
2616 new_nnz++;
2617
2618 if (ii + new_nnz > x_nz)
2619 {
2620 // Resize the sparse matrix
2621 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
2622 retval.change_capacity (sz);
2623 x_nz = sz;
2624 }
2625
2626 for (octave_idx_type i = 0; i < nc; i++)
2627 if (cwork[i] != 0.)
2628 {
2629 retval.xridx (ii) = i;
2630 retval.xdata (ii++) = cwork[i];
2631 }
2632 retval.xcidx (j+1) = ii;
2633 }
2634
2635 retval.maybe_compress ();
2636
2637 if (calc_cond)
2638 {
2639 // Calculation of 1-norm of inv(*this)
2640 OCTAVE_LOCAL_BUFFER (double, work, nm);
2641 for (octave_idx_type i = 0; i < nm; i++)
2642 work[i] = 0.;
2643
2644 for (octave_idx_type j = 0; j < nr; j++)
2645 {
2646 work[j] = 1.;
2647
2648 for (octave_idx_type k = j; k >= 0; k--)
2649 {
2650 if (work[k] != 0.)
2651 {
2652 double tmp = work[k] / data (cidx (k+1)-1);
2653 work[k] = tmp;
2654 for (octave_idx_type i = cidx (k);
2655 i < cidx (k+1)-1; i++)
2656 {
2657 octave_idx_type iidx = ridx (i);
2658 work[iidx] = work[iidx] - tmp * data (i);
2659 }
2660 }
2661 }
2662 double atmp = 0;
2663 for (octave_idx_type i = 0; i < j+1; i++)
2664 {
2665 atmp += fabs (work[i]);
2666 work[i] = 0.;
2667 }
2668 if (atmp > ainvnorm)
2669 ainvnorm = atmp;
2670 }
2671 rcond = 1. / ainvnorm / anorm;
2672 }
2673 }
2674
2675 triangular_error:
2676 if (err != 0)
2677 {
2678 if (sing_handler)
2679 {
2680 sing_handler (rcond);
2681 mattype.mark_as_rectangular ();
2682 }
2683 else
2684 warn_singular_matrix (rcond);
2685 }
2686
2687 volatile double rcond_plus_one = rcond + 1.0;
2688
2689 if (rcond_plus_one == 1.0 || xisnan (rcond))
2690 {
2691 err = -2;
2692
2693 if (sing_handler)
2694 {
2695 sing_handler (rcond);
2696 mattype.mark_as_rectangular ();
2697 }
2698 else
2699 warn_singular_matrix (rcond);
2700 }
2701 }
2702 else
2703 (*current_liboctave_error_handler) ("incorrect matrix type");
2704 } 2666 }
2705 2667
2706 return retval; 2668 return retval;
2707 } 2669 }
2708 2670
2720 err = 0; 2682 err = 0;
2721 2683
2722 if (nr != b.rows ()) 2684 if (nr != b.rows ())
2723 (*current_liboctave_error_handler) 2685 (*current_liboctave_error_handler)
2724 ("matrix dimension mismatch solution of linear equations"); 2686 ("matrix dimension mismatch solution of linear equations");
2725 else if (nr == 0 || nc == 0 || b.cols () == 0) 2687
2688 if (nr == 0 || nc == 0 || b.cols () == 0)
2726 retval = Matrix (nc, b.cols (), 0.0); 2689 retval = Matrix (nc, b.cols (), 0.0);
2727 else 2690 else
2728 { 2691 {
2729 // Print spparms("spumoni") info if requested 2692 // Print spparms("spumoni") info if requested
2730 int typ = mattype.type (); 2693 int typ = mattype.type ();
2731 mattype.info (); 2694 mattype.info ();
2732 2695
2733 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower) 2696 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2734 { 2697 (*current_liboctave_error_handler) ("incorrect matrix type");
2735 double anorm = 0.; 2698
2736 double ainvnorm = 0.; 2699 double anorm = 0.;
2737 octave_idx_type b_nc = b.cols (); 2700 double ainvnorm = 0.;
2738 rcond = 1.; 2701 octave_idx_type b_nc = b.cols ();
2702 rcond = 1.;
2703
2704 if (calc_cond)
2705 {
2706 // Calculate the 1-norm of matrix for rcond calculation
2707 for (octave_idx_type j = 0; j < nc; j++)
2708 {
2709 double atmp = 0.;
2710 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2711 atmp += fabs (data (i));
2712 if (atmp > anorm)
2713 anorm = atmp;
2714 }
2715 }
2716
2717 if (typ == MatrixType::Permuted_Lower)
2718 {
2719 retval.resize (nc, b_nc);
2720 OCTAVE_LOCAL_BUFFER (double, work, nm);
2721 octave_idx_type *perm = mattype.triangular_perm ();
2722
2723 for (octave_idx_type j = 0; j < b_nc; j++)
2724 {
2725 if (nc > nr)
2726 for (octave_idx_type i = 0; i < nm; i++)
2727 work[i] = 0.;
2728 for (octave_idx_type i = 0; i < nr; i++)
2729 work[perm[i]] = b(i,j);
2730
2731 for (octave_idx_type k = 0; k < nc; k++)
2732 {
2733 if (work[k] != 0.)
2734 {
2735 octave_idx_type minr = nr;
2736 octave_idx_type mini = 0;
2737
2738 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2739 if (perm[ridx (i)] < minr)
2740 {
2741 minr = perm[ridx (i)];
2742 mini = i;
2743 }
2744
2745 if (minr != k || data (mini) == 0)
2746 {
2747 err = -2;
2748 goto triangular_error;
2749 }
2750
2751 double tmp = work[k] / data (mini);
2752 work[k] = tmp;
2753 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2754 {
2755 if (i == mini)
2756 continue;
2757
2758 octave_idx_type iidx = perm[ridx (i)];
2759 work[iidx] = work[iidx] - tmp * data (i);
2760 }
2761 }
2762 }
2763
2764 for (octave_idx_type i = 0; i < nc; i++)
2765 retval(i, j) = work[i];
2766 }
2739 2767
2740 if (calc_cond) 2768 if (calc_cond)
2741 { 2769 {
2742 // Calculate the 1-norm of matrix for rcond calculation 2770 // Calculation of 1-norm of inv(*this)
2743 for (octave_idx_type j = 0; j < nc; j++) 2771 for (octave_idx_type i = 0; i < nm; i++)
2744 { 2772 work[i] = 0.;
2745 double atmp = 0.; 2773
2746 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 2774 for (octave_idx_type j = 0; j < nr; j++)
2747 atmp += fabs (data (i)); 2775 {
2748 if (atmp > anorm) 2776 work[j] = 1.;
2749 anorm = atmp;
2750 }
2751 }
2752
2753 if (typ == MatrixType::Permuted_Lower)
2754 {
2755 retval.resize (nc, b_nc);
2756 OCTAVE_LOCAL_BUFFER (double, work, nm);
2757 octave_idx_type *perm = mattype.triangular_perm ();
2758
2759 for (octave_idx_type j = 0; j < b_nc; j++)
2760 {
2761 if (nc > nr)
2762 for (octave_idx_type i = 0; i < nm; i++)
2763 work[i] = 0.;
2764 for (octave_idx_type i = 0; i < nr; i++)
2765 work[perm[i]] = b(i,j);
2766 2777
2767 for (octave_idx_type k = 0; k < nc; k++) 2778 for (octave_idx_type k = 0; k < nc; k++)
2768 { 2779 {
2769 if (work[k] != 0.) 2780 if (work[k] != 0.)
2770 { 2781 {
2771 octave_idx_type minr = nr; 2782 octave_idx_type minr = nr;
2772 octave_idx_type mini = 0; 2783 octave_idx_type mini = 0;
2773 2784
2774 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 2785 for (octave_idx_type i = cidx (k);
2786 i < cidx (k+1); i++)
2775 if (perm[ridx (i)] < minr) 2787 if (perm[ridx (i)] < minr)
2776 { 2788 {
2777 minr = perm[ridx (i)]; 2789 minr = perm[ridx (i)];
2778 mini = i; 2790 mini = i;
2779 } 2791 }
2780 2792
2781 if (minr != k || data (mini) == 0)
2782 {
2783 err = -2;
2784 goto triangular_error;
2785 }
2786
2787 double tmp = work[k] / data (mini); 2793 double tmp = work[k] / data (mini);
2788 work[k] = tmp; 2794 work[k] = tmp;
2789 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 2795 for (octave_idx_type i = cidx (k);
2796 i < cidx (k+1); i++)
2790 { 2797 {
2791 if (i == mini) 2798 if (i == mini)
2792 continue; 2799 continue;
2793 2800
2794 octave_idx_type iidx = perm[ridx (i)]; 2801 octave_idx_type iidx = perm[ridx (i)];
2795 work[iidx] = work[iidx] - tmp * data (i); 2802 work[iidx] = work[iidx] - tmp * data (i);
2796 } 2803 }
2797 } 2804 }
2798 } 2805 }
2799 2806
2800 for (octave_idx_type i = 0; i < nc; i++) 2807 double atmp = 0;
2801 retval(i, j) = work[i]; 2808 for (octave_idx_type i = j; i < nc; i++)
2802 } 2809 {
2803 2810 atmp += fabs (work[i]);
2804 if (calc_cond) 2811 work[i] = 0.;
2805 { 2812 }
2806 // Calculation of 1-norm of inv(*this) 2813 if (atmp > ainvnorm)
2807 for (octave_idx_type i = 0; i < nm; i++) 2814 ainvnorm = atmp;
2808 work[i] = 0.; 2815 }
2809 2816 rcond = 1. / ainvnorm / anorm;
2810 for (octave_idx_type j = 0; j < nr; j++) 2817 }
2811 { 2818 }
2812 work[j] = 1.; 2819 else
2813 2820 {
2814 for (octave_idx_type k = 0; k < nc; k++) 2821 OCTAVE_LOCAL_BUFFER (double, work, nm);
2822 retval.resize (nc, b_nc, 0.);
2823
2824 for (octave_idx_type j = 0; j < b_nc; j++)
2825 {
2826 for (octave_idx_type i = 0; i < nr; i++)
2827 work[i] = b(i,j);
2828 for (octave_idx_type i = nr; i < nc; i++)
2829 work[i] = 0.;
2830 for (octave_idx_type k = 0; k < nc; k++)
2831 {
2832 if (work[k] != 0.)
2833 {
2834 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2815 { 2835 {
2816 if (work[k] != 0.) 2836 err = -2;
2817 { 2837 goto triangular_error;
2818 octave_idx_type minr = nr;
2819 octave_idx_type mini = 0;
2820
2821 for (octave_idx_type i = cidx (k);
2822 i < cidx (k+1); i++)
2823 if (perm[ridx (i)] < minr)
2824 {
2825 minr = perm[ridx (i)];
2826 mini = i;
2827 }
2828
2829 double tmp = work[k] / data (mini);
2830 work[k] = tmp;
2831 for (octave_idx_type i = cidx (k);
2832 i < cidx (k+1); i++)
2833 {
2834 if (i == mini)
2835 continue;
2836
2837 octave_idx_type iidx = perm[ridx (i)];
2838 work[iidx] = work[iidx] - tmp * data (i);
2839 }
2840 }
2841 } 2838 }
2842 2839
2843 double atmp = 0; 2840 double tmp = work[k] / data (cidx (k));
2844 for (octave_idx_type i = j; i < nc; i++) 2841 work[k] = tmp;
2842 for (octave_idx_type i = cidx (k)+1;
2843 i < cidx (k+1); i++)
2845 { 2844 {
2846 atmp += fabs (work[i]); 2845 octave_idx_type iidx = ridx (i);
2847 work[i] = 0.; 2846 work[iidx] = work[iidx] - tmp * data (i);
2848 } 2847 }
2849 if (atmp > ainvnorm) 2848 }
2850 ainvnorm = atmp; 2849 }
2851 } 2850
2852 rcond = 1. / ainvnorm / anorm; 2851 for (octave_idx_type i = 0; i < nc; i++)
2853 } 2852 retval.xelem (i, j) = work[i];
2854 } 2853 }
2855 else 2854
2856 { 2855 if (calc_cond)
2857 OCTAVE_LOCAL_BUFFER (double, work, nm); 2856 {
2858 retval.resize (nc, b_nc, 0.); 2857 // Calculation of 1-norm of inv(*this)
2859 2858 for (octave_idx_type i = 0; i < nm; i++)
2860 for (octave_idx_type j = 0; j < b_nc; j++) 2859 work[i] = 0.;
2861 { 2860
2862 for (octave_idx_type i = 0; i < nr; i++) 2861 for (octave_idx_type j = 0; j < nr; j++)
2863 work[i] = b(i,j); 2862 {
2864 for (octave_idx_type i = nr; i < nc; i++) 2863 work[j] = 1.;
2865 work[i] = 0.; 2864
2866 for (octave_idx_type k = 0; k < nc; k++) 2865 for (octave_idx_type k = j; k < nc; k++)
2867 { 2866 {
2867
2868 if (work[k] != 0.) 2868 if (work[k] != 0.)
2869 { 2869 {
2870 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
2871 {
2872 err = -2;
2873 goto triangular_error;
2874 }
2875
2876 double tmp = work[k] / data (cidx (k)); 2870 double tmp = work[k] / data (cidx (k));
2877 work[k] = tmp; 2871 work[k] = tmp;
2878 for (octave_idx_type i = cidx (k)+1; 2872 for (octave_idx_type i = cidx (k)+1;
2879 i < cidx (k+1); i++) 2873 i < cidx (k+1); i++)
2880 { 2874 {
2881 octave_idx_type iidx = ridx (i); 2875 octave_idx_type iidx = ridx (i);
2882 work[iidx] = work[iidx] - tmp * data (i); 2876 work[iidx] = work[iidx] - tmp * data (i);
2883 } 2877 }
2884 } 2878 }
2885 } 2879 }
2886 2880 double atmp = 0;
2887 for (octave_idx_type i = 0; i < nc; i++) 2881 for (octave_idx_type i = j; i < nc; i++)
2888 retval.xelem (i, j) = work[i]; 2882 {
2889 } 2883 atmp += fabs (work[i]);
2890 2884 work[i] = 0.;
2891 if (calc_cond) 2885 }
2892 { 2886 if (atmp > ainvnorm)
2893 // Calculation of 1-norm of inv(*this) 2887 ainvnorm = atmp;
2894 for (octave_idx_type i = 0; i < nm; i++) 2888 }
2895 work[i] = 0.; 2889 rcond = 1. / ainvnorm / anorm;
2896 2890 }
2897 for (octave_idx_type j = 0; j < nr; j++) 2891 }
2898 { 2892
2899 work[j] = 1.; 2893 triangular_error:
2900 2894 if (err != 0)
2901 for (octave_idx_type k = j; k < nc; k++) 2895 {
2902 { 2896 if (sing_handler)
2903 2897 {
2904 if (work[k] != 0.) 2898 sing_handler (rcond);
2905 { 2899 mattype.mark_as_rectangular ();
2906 double tmp = work[k] / data (cidx (k)); 2900 }
2907 work[k] = tmp; 2901 else
2908 for (octave_idx_type i = cidx (k)+1; 2902 warn_singular_matrix (rcond);
2909 i < cidx (k+1); i++) 2903 }
2910 { 2904
2911 octave_idx_type iidx = ridx (i); 2905 volatile double rcond_plus_one = rcond + 1.0;
2912 work[iidx] = work[iidx] - tmp * data (i); 2906
2913 } 2907 if (rcond_plus_one == 1.0 || xisnan (rcond))
2914 } 2908 {
2915 } 2909 err = -2;
2916 double atmp = 0; 2910
2917 for (octave_idx_type i = j; i < nc; i++) 2911 if (sing_handler)
2918 { 2912 {
2919 atmp += fabs (work[i]); 2913 sing_handler (rcond);
2920 work[i] = 0.; 2914 mattype.mark_as_rectangular ();
2921 } 2915 }
2922 if (atmp > ainvnorm) 2916 else
2923 ainvnorm = atmp; 2917 warn_singular_matrix (rcond);
2924 } 2918 }
2925 rcond = 1. / ainvnorm / anorm;
2926 }
2927 }
2928
2929 triangular_error:
2930 if (err != 0)
2931 {
2932 if (sing_handler)
2933 {
2934 sing_handler (rcond);
2935 mattype.mark_as_rectangular ();
2936 }
2937 else
2938 warn_singular_matrix (rcond);
2939 }
2940
2941 volatile double rcond_plus_one = rcond + 1.0;
2942
2943 if (rcond_plus_one == 1.0 || xisnan (rcond))
2944 {
2945 err = -2;
2946
2947 if (sing_handler)
2948 {
2949 sing_handler (rcond);
2950 mattype.mark_as_rectangular ();
2951 }
2952 else
2953 warn_singular_matrix (rcond);
2954 }
2955 }
2956 else
2957 (*current_liboctave_error_handler) ("incorrect matrix type");
2958 } 2919 }
2959 2920
2960 return retval; 2921 return retval;
2961 } 2922 }
2962 2923
2974 err = 0; 2935 err = 0;
2975 2936
2976 if (nr != b.rows ()) 2937 if (nr != b.rows ())
2977 (*current_liboctave_error_handler) 2938 (*current_liboctave_error_handler)
2978 ("matrix dimension mismatch solution of linear equations"); 2939 ("matrix dimension mismatch solution of linear equations");
2979 else if (nr == 0 || nc == 0 || b.cols () == 0) 2940
2941 if (nr == 0 || nc == 0 || b.cols () == 0)
2980 retval = SparseMatrix (nc, b.cols ()); 2942 retval = SparseMatrix (nc, b.cols ());
2981 else 2943 else
2982 { 2944 {
2983 // Print spparms("spumoni") info if requested 2945 // Print spparms("spumoni") info if requested
2984 int typ = mattype.type (); 2946 int typ = mattype.type ();
2985 mattype.info (); 2947 mattype.info ();
2986 2948
2987 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower) 2949 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
2988 { 2950 (*current_liboctave_error_handler) ("incorrect matrix type");
2989 double anorm = 0.; 2951
2990 double ainvnorm = 0.; 2952 double anorm = 0.;
2991 rcond = 1.; 2953 double ainvnorm = 0.;
2954 rcond = 1.;
2955
2956 if (calc_cond)
2957 {
2958 // Calculate the 1-norm of matrix for rcond calculation
2959 for (octave_idx_type j = 0; j < nc; j++)
2960 {
2961 double atmp = 0.;
2962 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
2963 atmp += fabs (data (i));
2964 if (atmp > anorm)
2965 anorm = atmp;
2966 }
2967 }
2968
2969 octave_idx_type b_nc = b.cols ();
2970 octave_idx_type b_nz = b.nnz ();
2971 retval = SparseMatrix (nc, b_nc, b_nz);
2972 retval.xcidx (0) = 0;
2973 octave_idx_type ii = 0;
2974 octave_idx_type x_nz = b_nz;
2975
2976 if (typ == MatrixType::Permuted_Lower)
2977 {
2978 OCTAVE_LOCAL_BUFFER (double, work, nm);
2979 octave_idx_type *perm = mattype.triangular_perm ();
2980
2981 for (octave_idx_type j = 0; j < b_nc; j++)
2982 {
2983 for (octave_idx_type i = 0; i < nm; i++)
2984 work[i] = 0.;
2985 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
2986 work[perm[b.ridx (i)]] = b.data (i);
2987
2988 for (octave_idx_type k = 0; k < nc; k++)
2989 {
2990 if (work[k] != 0.)
2991 {
2992 octave_idx_type minr = nr;
2993 octave_idx_type mini = 0;
2994
2995 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
2996 if (perm[ridx (i)] < minr)
2997 {
2998 minr = perm[ridx (i)];
2999 mini = i;
3000 }
3001
3002 if (minr != k || data (mini) == 0)
3003 {
3004 err = -2;
3005 goto triangular_error;
3006 }
3007
3008 double tmp = work[k] / data (mini);
3009 work[k] = tmp;
3010 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3011 {
3012 if (i == mini)
3013 continue;
3014
3015 octave_idx_type iidx = perm[ridx (i)];
3016 work[iidx] = work[iidx] - tmp * data (i);
3017 }
3018 }
3019 }
3020
3021 // Count nonzeros in work vector and adjust space in
3022 // retval if needed
3023 octave_idx_type new_nnz = 0;
3024 for (octave_idx_type i = 0; i < nc; i++)
3025 if (work[i] != 0.)
3026 new_nnz++;
3027
3028 if (ii + new_nnz > x_nz)
3029 {
3030 // Resize the sparse matrix
3031 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3032 retval.change_capacity (sz);
3033 x_nz = sz;
3034 }
3035
3036 for (octave_idx_type i = 0; i < nc; i++)
3037 if (work[i] != 0.)
3038 {
3039 retval.xridx (ii) = i;
3040 retval.xdata (ii++) = work[i];
3041 }
3042 retval.xcidx (j+1) = ii;
3043 }
3044
3045 retval.maybe_compress ();
2992 3046
2993 if (calc_cond) 3047 if (calc_cond)
2994 { 3048 {
2995 // Calculate the 1-norm of matrix for rcond calculation 3049 // Calculation of 1-norm of inv(*this)
2996 for (octave_idx_type j = 0; j < nc; j++) 3050 for (octave_idx_type i = 0; i < nm; i++)
2997 { 3051 work[i] = 0.;
2998 double atmp = 0.; 3052
2999 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 3053 for (octave_idx_type j = 0; j < nr; j++)
3000 atmp += fabs (data (i)); 3054 {
3001 if (atmp > anorm) 3055 work[j] = 1.;
3002 anorm = atmp;
3003 }
3004 }
3005
3006 octave_idx_type b_nc = b.cols ();
3007 octave_idx_type b_nz = b.nnz ();
3008 retval = SparseMatrix (nc, b_nc, b_nz);
3009 retval.xcidx (0) = 0;
3010 octave_idx_type ii = 0;
3011 octave_idx_type x_nz = b_nz;
3012
3013 if (typ == MatrixType::Permuted_Lower)
3014 {
3015 OCTAVE_LOCAL_BUFFER (double, work, nm);
3016 octave_idx_type *perm = mattype.triangular_perm ();
3017
3018 for (octave_idx_type j = 0; j < b_nc; j++)
3019 {
3020 for (octave_idx_type i = 0; i < nm; i++)
3021 work[i] = 0.;
3022 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3023 work[perm[b.ridx (i)]] = b.data (i);
3024 3056
3025 for (octave_idx_type k = 0; k < nc; k++) 3057 for (octave_idx_type k = 0; k < nc; k++)
3026 { 3058 {
3027 if (work[k] != 0.) 3059 if (work[k] != 0.)
3028 { 3060 {
3029 octave_idx_type minr = nr; 3061 octave_idx_type minr = nr;
3030 octave_idx_type mini = 0; 3062 octave_idx_type mini = 0;
3031 3063
3032 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 3064 for (octave_idx_type i = cidx (k);
3065 i < cidx (k+1); i++)
3033 if (perm[ridx (i)] < minr) 3066 if (perm[ridx (i)] < minr)
3034 { 3067 {
3035 minr = perm[ridx (i)]; 3068 minr = perm[ridx (i)];
3036 mini = i; 3069 mini = i;
3037 } 3070 }
3038 3071
3039 if (minr != k || data (mini) == 0)
3040 {
3041 err = -2;
3042 goto triangular_error;
3043 }
3044
3045 double tmp = work[k] / data (mini); 3072 double tmp = work[k] / data (mini);
3046 work[k] = tmp; 3073 work[k] = tmp;
3047 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 3074 for (octave_idx_type i = cidx (k);
3075 i < cidx (k+1); i++)
3048 { 3076 {
3049 if (i == mini) 3077 if (i == mini)
3050 continue; 3078 continue;
3051 3079
3052 octave_idx_type iidx = perm[ridx (i)]; 3080 octave_idx_type iidx = perm[ridx (i)];
3053 work[iidx] = work[iidx] - tmp * data (i); 3081 work[iidx] = work[iidx] - tmp * data (i);
3054 } 3082 }
3055 } 3083 }
3056 } 3084 }
3057 3085
3058 // Count nonzeros in work vector and adjust space in 3086 double atmp = 0;
3059 // retval if needed 3087 for (octave_idx_type i = j; i < nr; i++)
3060 octave_idx_type new_nnz = 0; 3088 {
3061 for (octave_idx_type i = 0; i < nc; i++) 3089 atmp += fabs (work[i]);
3062 if (work[i] != 0.) 3090 work[i] = 0.;
3063 new_nnz++; 3091 }
3064 3092 if (atmp > ainvnorm)
3065 if (ii + new_nnz > x_nz) 3093 ainvnorm = atmp;
3066 { 3094 }
3067 // Resize the sparse matrix 3095 rcond = 1. / ainvnorm / anorm;
3068 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 3096 }
3069 retval.change_capacity (sz); 3097 }
3070 x_nz = sz; 3098 else
3071 } 3099 {
3072 3100 OCTAVE_LOCAL_BUFFER (double, work, nm);
3073 for (octave_idx_type i = 0; i < nc; i++) 3101
3074 if (work[i] != 0.) 3102 for (octave_idx_type j = 0; j < b_nc; j++)
3075 { 3103 {
3076 retval.xridx (ii) = i; 3104 for (octave_idx_type i = 0; i < nm; i++)
3077 retval.xdata (ii++) = work[i]; 3105 work[i] = 0.;
3078 } 3106 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3079 retval.xcidx (j+1) = ii; 3107 work[b.ridx (i)] = b.data (i);
3080 } 3108
3081 3109 for (octave_idx_type k = 0; k < nc; k++)
3082 retval.maybe_compress (); 3110 {
3083 3111 if (work[k] != 0.)
3084 if (calc_cond) 3112 {
3085 { 3113 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3086 // Calculation of 1-norm of inv(*this)
3087 for (octave_idx_type i = 0; i < nm; i++)
3088 work[i] = 0.;
3089
3090 for (octave_idx_type j = 0; j < nr; j++)
3091 {
3092 work[j] = 1.;
3093
3094 for (octave_idx_type k = 0; k < nc; k++)
3095 { 3114 {
3096 if (work[k] != 0.) 3115 err = -2;
3097 { 3116 goto triangular_error;
3098 octave_idx_type minr = nr;
3099 octave_idx_type mini = 0;
3100
3101 for (octave_idx_type i = cidx (k);
3102 i < cidx (k+1); i++)
3103 if (perm[ridx (i)] < minr)
3104 {
3105 minr = perm[ridx (i)];
3106 mini = i;
3107 }
3108
3109 double tmp = work[k] / data (mini);
3110 work[k] = tmp;
3111 for (octave_idx_type i = cidx (k);
3112 i < cidx (k+1); i++)
3113 {
3114 if (i == mini)
3115 continue;
3116
3117 octave_idx_type iidx = perm[ridx (i)];
3118 work[iidx] = work[iidx] - tmp * data (i);
3119 }
3120 }
3121 } 3117 }
3122 3118
3123 double atmp = 0; 3119 double tmp = work[k] / data (cidx (k));
3124 for (octave_idx_type i = j; i < nr; i++) 3120 work[k] = tmp;
3121 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3125 { 3122 {
3126 atmp += fabs (work[i]); 3123 octave_idx_type iidx = ridx (i);
3127 work[i] = 0.; 3124 work[iidx] = work[iidx] - tmp * data (i);
3128 } 3125 }
3129 if (atmp > ainvnorm) 3126 }
3130 ainvnorm = atmp; 3127 }
3131 } 3128
3132 rcond = 1. / ainvnorm / anorm; 3129 // Count nonzeros in work vector and adjust space in
3133 } 3130 // retval if needed
3134 } 3131 octave_idx_type new_nnz = 0;
3135 else 3132 for (octave_idx_type i = 0; i < nc; i++)
3136 { 3133 if (work[i] != 0.)
3137 OCTAVE_LOCAL_BUFFER (double, work, nm); 3134 new_nnz++;
3138 3135
3139 for (octave_idx_type j = 0; j < b_nc; j++) 3136 if (ii + new_nnz > x_nz)
3140 { 3137 {
3141 for (octave_idx_type i = 0; i < nm; i++) 3138 // Resize the sparse matrix
3142 work[i] = 0.; 3139 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3143 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 3140 retval.change_capacity (sz);
3144 work[b.ridx (i)] = b.data (i); 3141 x_nz = sz;
3145 3142 }
3146 for (octave_idx_type k = 0; k < nc; k++) 3143
3147 { 3144 for (octave_idx_type i = 0; i < nc; i++)
3145 if (work[i] != 0.)
3146 {
3147 retval.xridx (ii) = i;
3148 retval.xdata (ii++) = work[i];
3149 }
3150 retval.xcidx (j+1) = ii;
3151 }
3152
3153 retval.maybe_compress ();
3154
3155 if (calc_cond)
3156 {
3157 // Calculation of 1-norm of inv(*this)
3158 for (octave_idx_type i = 0; i < nm; i++)
3159 work[i] = 0.;
3160
3161 for (octave_idx_type j = 0; j < nr; j++)
3162 {
3163 work[j] = 1.;
3164
3165 for (octave_idx_type k = j; k < nc; k++)
3166 {
3167
3148 if (work[k] != 0.) 3168 if (work[k] != 0.)
3149 { 3169 {
3150 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3151 {
3152 err = -2;
3153 goto triangular_error;
3154 }
3155
3156 double tmp = work[k] / data (cidx (k)); 3170 double tmp = work[k] / data (cidx (k));
3157 work[k] = tmp; 3171 work[k] = tmp;
3158 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) 3172 for (octave_idx_type i = cidx (k)+1;
3173 i < cidx (k+1); i++)
3159 { 3174 {
3160 octave_idx_type iidx = ridx (i); 3175 octave_idx_type iidx = ridx (i);
3161 work[iidx] = work[iidx] - tmp * data (i); 3176 work[iidx] = work[iidx] - tmp * data (i);
3162 } 3177 }
3163 } 3178 }
3164 } 3179 }
3165 3180 double atmp = 0;
3166 // Count nonzeros in work vector and adjust space in 3181 for (octave_idx_type i = j; i < nc; i++)
3167 // retval if needed 3182 {
3168 octave_idx_type new_nnz = 0; 3183 atmp += fabs (work[i]);
3169 for (octave_idx_type i = 0; i < nc; i++) 3184 work[i] = 0.;
3170 if (work[i] != 0.) 3185 }
3171 new_nnz++; 3186 if (atmp > ainvnorm)
3172 3187 ainvnorm = atmp;
3173 if (ii + new_nnz > x_nz) 3188 }
3174 { 3189 rcond = 1. / ainvnorm / anorm;
3175 // Resize the sparse matrix 3190 }
3176 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 3191 }
3177 retval.change_capacity (sz); 3192
3178 x_nz = sz; 3193 triangular_error:
3179 } 3194 if (err != 0)
3180 3195 {
3181 for (octave_idx_type i = 0; i < nc; i++) 3196 if (sing_handler)
3182 if (work[i] != 0.) 3197 {
3183 { 3198 sing_handler (rcond);
3184 retval.xridx (ii) = i; 3199 mattype.mark_as_rectangular ();
3185 retval.xdata (ii++) = work[i]; 3200 }
3186 } 3201 else
3187 retval.xcidx (j+1) = ii; 3202 warn_singular_matrix (rcond);
3188 } 3203 }
3189 3204
3190 retval.maybe_compress (); 3205 volatile double rcond_plus_one = rcond + 1.0;
3191 3206
3192 if (calc_cond) 3207 if (rcond_plus_one == 1.0 || xisnan (rcond))
3193 { 3208 {
3194 // Calculation of 1-norm of inv(*this) 3209 err = -2;
3195 for (octave_idx_type i = 0; i < nm; i++) 3210
3196 work[i] = 0.; 3211 if (sing_handler)
3197 3212 {
3198 for (octave_idx_type j = 0; j < nr; j++) 3213 sing_handler (rcond);
3199 { 3214 mattype.mark_as_rectangular ();
3200 work[j] = 1.; 3215 }
3201 3216 else
3202 for (octave_idx_type k = j; k < nc; k++) 3217 warn_singular_matrix (rcond);
3203 { 3218 }
3204
3205 if (work[k] != 0.)
3206 {
3207 double tmp = work[k] / data (cidx (k));
3208 work[k] = tmp;
3209 for (octave_idx_type i = cidx (k)+1;
3210 i < cidx (k+1); i++)
3211 {
3212 octave_idx_type iidx = ridx (i);
3213 work[iidx] = work[iidx] - tmp * data (i);
3214 }
3215 }
3216 }
3217 double atmp = 0;
3218 for (octave_idx_type i = j; i < nc; i++)
3219 {
3220 atmp += fabs (work[i]);
3221 work[i] = 0.;
3222 }
3223 if (atmp > ainvnorm)
3224 ainvnorm = atmp;
3225 }
3226 rcond = 1. / ainvnorm / anorm;
3227 }
3228 }
3229
3230 triangular_error:
3231 if (err != 0)
3232 {
3233 if (sing_handler)
3234 {
3235 sing_handler (rcond);
3236 mattype.mark_as_rectangular ();
3237 }
3238 else
3239 warn_singular_matrix (rcond);
3240 }
3241
3242 volatile double rcond_plus_one = rcond + 1.0;
3243
3244 if (rcond_plus_one == 1.0 || xisnan (rcond))
3245 {
3246 err = -2;
3247
3248 if (sing_handler)
3249 {
3250 sing_handler (rcond);
3251 mattype.mark_as_rectangular ();
3252 }
3253 else
3254 warn_singular_matrix (rcond);
3255 }
3256 }
3257 else
3258 (*current_liboctave_error_handler) ("incorrect matrix type");
3259 } 3219 }
3260 3220
3261 return retval; 3221 return retval;
3262 } 3222 }
3263 3223
3275 err = 0; 3235 err = 0;
3276 3236
3277 if (nr != b.rows ()) 3237 if (nr != b.rows ())
3278 (*current_liboctave_error_handler) 3238 (*current_liboctave_error_handler)
3279 ("matrix dimension mismatch solution of linear equations"); 3239 ("matrix dimension mismatch solution of linear equations");
3280 else if (nr == 0 || nc == 0 || b.cols () == 0) 3240
3241 if (nr == 0 || nc == 0 || b.cols () == 0)
3281 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 3242 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
3282 else 3243 else
3283 { 3244 {
3284 // Print spparms("spumoni") info if requested 3245 // Print spparms("spumoni") info if requested
3285 int typ = mattype.type (); 3246 int typ = mattype.type ();
3286 mattype.info (); 3247 mattype.info ();
3287 3248
3288 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower) 3249 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3289 { 3250 (*current_liboctave_error_handler) ("incorrect matrix type");
3290 double anorm = 0.; 3251
3291 double ainvnorm = 0.; 3252 double anorm = 0.;
3292 octave_idx_type b_nc = b.cols (); 3253 double ainvnorm = 0.;
3293 rcond = 1.; 3254 octave_idx_type b_nc = b.cols ();
3255 rcond = 1.;
3256
3257 if (calc_cond)
3258 {
3259 // Calculate the 1-norm of matrix for rcond calculation
3260 for (octave_idx_type j = 0; j < nc; j++)
3261 {
3262 double atmp = 0.;
3263 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3264 atmp += fabs (data (i));
3265 if (atmp > anorm)
3266 anorm = atmp;
3267 }
3268 }
3269
3270 if (typ == MatrixType::Permuted_Lower)
3271 {
3272 retval.resize (nc, b_nc);
3273 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3274 octave_idx_type *perm = mattype.triangular_perm ();
3275
3276 for (octave_idx_type j = 0; j < b_nc; j++)
3277 {
3278 for (octave_idx_type i = 0; i < nm; i++)
3279 cwork[i] = 0.;
3280 for (octave_idx_type i = 0; i < nr; i++)
3281 cwork[perm[i]] = b(i,j);
3282
3283 for (octave_idx_type k = 0; k < nc; k++)
3284 {
3285 if (cwork[k] != 0.)
3286 {
3287 octave_idx_type minr = nr;
3288 octave_idx_type mini = 0;
3289
3290 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3291 if (perm[ridx (i)] < minr)
3292 {
3293 minr = perm[ridx (i)];
3294 mini = i;
3295 }
3296
3297 if (minr != k || data (mini) == 0)
3298 {
3299 err = -2;
3300 goto triangular_error;
3301 }
3302
3303 Complex tmp = cwork[k] / data (mini);
3304 cwork[k] = tmp;
3305 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3306 {
3307 if (i == mini)
3308 continue;
3309
3310 octave_idx_type iidx = perm[ridx (i)];
3311 cwork[iidx] = cwork[iidx] - tmp * data (i);
3312 }
3313 }
3314 }
3315
3316 for (octave_idx_type i = 0; i < nc; i++)
3317 retval(i, j) = cwork[i];
3318 }
3294 3319
3295 if (calc_cond) 3320 if (calc_cond)
3296 { 3321 {
3297 // Calculate the 1-norm of matrix for rcond calculation 3322 // Calculation of 1-norm of inv(*this)
3298 for (octave_idx_type j = 0; j < nc; j++) 3323 OCTAVE_LOCAL_BUFFER (double, work, nm);
3299 { 3324 for (octave_idx_type i = 0; i < nm; i++)
3300 double atmp = 0.; 3325 work[i] = 0.;
3301 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 3326
3302 atmp += fabs (data (i)); 3327 for (octave_idx_type j = 0; j < nr; j++)
3303 if (atmp > anorm) 3328 {
3304 anorm = atmp; 3329 work[j] = 1.;
3305 }
3306 }
3307
3308 if (typ == MatrixType::Permuted_Lower)
3309 {
3310 retval.resize (nc, b_nc);
3311 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3312 octave_idx_type *perm = mattype.triangular_perm ();
3313
3314 for (octave_idx_type j = 0; j < b_nc; j++)
3315 {
3316 for (octave_idx_type i = 0; i < nm; i++)
3317 cwork[i] = 0.;
3318 for (octave_idx_type i = 0; i < nr; i++)
3319 cwork[perm[i]] = b(i,j);
3320 3330
3321 for (octave_idx_type k = 0; k < nc; k++) 3331 for (octave_idx_type k = 0; k < nc; k++)
3322 { 3332 {
3323 if (cwork[k] != 0.) 3333 if (work[k] != 0.)
3324 { 3334 {
3325 octave_idx_type minr = nr; 3335 octave_idx_type minr = nr;
3326 octave_idx_type mini = 0; 3336 octave_idx_type mini = 0;
3327 3337
3328 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 3338 for (octave_idx_type i = cidx (k);
3339 i < cidx (k+1); i++)
3329 if (perm[ridx (i)] < minr) 3340 if (perm[ridx (i)] < minr)
3330 { 3341 {
3331 minr = perm[ridx (i)]; 3342 minr = perm[ridx (i)];
3332 mini = i; 3343 mini = i;
3333 } 3344 }
3334 3345
3335 if (minr != k || data (mini) == 0) 3346 double tmp = work[k] / data (mini);
3336 { 3347 work[k] = tmp;
3337 err = -2; 3348 for (octave_idx_type i = cidx (k);
3338 goto triangular_error; 3349 i < cidx (k+1); i++)
3339 }
3340
3341 Complex tmp = cwork[k] / data (mini);
3342 cwork[k] = tmp;
3343 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3344 { 3350 {
3345 if (i == mini) 3351 if (i == mini)
3346 continue; 3352 continue;
3347 3353
3348 octave_idx_type iidx = perm[ridx (i)]; 3354 octave_idx_type iidx = perm[ridx (i)];
3349 cwork[iidx] = cwork[iidx] - tmp * data (i); 3355 work[iidx] = work[iidx] - tmp * data (i);
3350 } 3356 }
3351 } 3357 }
3352 } 3358 }
3353 3359
3354 for (octave_idx_type i = 0; i < nc; i++) 3360 double atmp = 0;
3355 retval(i, j) = cwork[i]; 3361 for (octave_idx_type i = j; i < nc; i++)
3356 } 3362 {
3357 3363 atmp += fabs (work[i]);
3358 if (calc_cond) 3364 work[i] = 0.;
3359 { 3365 }
3360 // Calculation of 1-norm of inv(*this) 3366 if (atmp > ainvnorm)
3361 OCTAVE_LOCAL_BUFFER (double, work, nm); 3367 ainvnorm = atmp;
3362 for (octave_idx_type i = 0; i < nm; i++) 3368 }
3363 work[i] = 0.; 3369 rcond = 1. / ainvnorm / anorm;
3364 3370 }
3365 for (octave_idx_type j = 0; j < nr; j++) 3371 }
3366 { 3372 else
3367 work[j] = 1.; 3373 {
3368 3374 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3369 for (octave_idx_type k = 0; k < nc; k++) 3375 retval.resize (nc, b_nc, 0.);
3376
3377 for (octave_idx_type j = 0; j < b_nc; j++)
3378 {
3379 for (octave_idx_type i = 0; i < nr; i++)
3380 cwork[i] = b(i,j);
3381 for (octave_idx_type i = nr; i < nc; i++)
3382 cwork[i] = 0.;
3383
3384 for (octave_idx_type k = 0; k < nc; k++)
3385 {
3386 if (cwork[k] != 0.)
3387 {
3388 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3370 { 3389 {
3371 if (work[k] != 0.) 3390 err = -2;
3391 goto triangular_error;
3392 }
3393
3394 Complex tmp = cwork[k] / data (cidx (k));
3395 cwork[k] = tmp;
3396 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3397 {
3398 octave_idx_type iidx = ridx (i);
3399 cwork[iidx] = cwork[iidx] - tmp * data (i);
3400 }
3401 }
3402 }
3403
3404 for (octave_idx_type i = 0; i < nc; i++)
3405 retval.xelem (i, j) = cwork[i];
3406 }
3407
3408 if (calc_cond)
3409 {
3410 // Calculation of 1-norm of inv(*this)
3411 OCTAVE_LOCAL_BUFFER (double, work, nm);
3412 for (octave_idx_type i = 0; i < nm; i++)
3413 work[i] = 0.;
3414
3415 for (octave_idx_type j = 0; j < nr; j++)
3416 {
3417 work[j] = 1.;
3418
3419 for (octave_idx_type k = j; k < nc; k++)
3420 {
3421
3422 if (work[k] != 0.)
3423 {
3424 double tmp = work[k] / data (cidx (k));
3425 work[k] = tmp;
3426 for (octave_idx_type i = cidx (k)+1;
3427 i < cidx (k+1); i++)
3372 { 3428 {
3373 octave_idx_type minr = nr; 3429 octave_idx_type iidx = ridx (i);
3374 octave_idx_type mini = 0; 3430 work[iidx] = work[iidx] - tmp * data (i);
3375
3376 for (octave_idx_type i = cidx (k);
3377 i < cidx (k+1); i++)
3378 if (perm[ridx (i)] < minr)
3379 {
3380 minr = perm[ridx (i)];
3381 mini = i;
3382 }
3383
3384 double tmp = work[k] / data (mini);
3385 work[k] = tmp;
3386 for (octave_idx_type i = cidx (k);
3387 i < cidx (k+1); i++)
3388 {
3389 if (i == mini)
3390 continue;
3391
3392 octave_idx_type iidx = perm[ridx (i)];
3393 work[iidx] = work[iidx] - tmp * data (i);
3394 }
3395 } 3431 }
3396 } 3432 }
3397 3433 }
3398 double atmp = 0; 3434 double atmp = 0;
3399 for (octave_idx_type i = j; i < nc; i++) 3435 for (octave_idx_type i = j; i < nc; i++)
3400 { 3436 {
3401 atmp += fabs (work[i]); 3437 atmp += fabs (work[i]);
3402 work[i] = 0.; 3438 work[i] = 0.;
3403 } 3439 }
3404 if (atmp > ainvnorm) 3440 if (atmp > ainvnorm)
3405 ainvnorm = atmp; 3441 ainvnorm = atmp;
3406 } 3442 }
3407 rcond = 1. / ainvnorm / anorm; 3443 rcond = 1. / ainvnorm / anorm;
3408 } 3444 }
3445 }
3446
3447 triangular_error:
3448 if (err != 0)
3449 {
3450 if (sing_handler)
3451 {
3452 sing_handler (rcond);
3453 mattype.mark_as_rectangular ();
3409 } 3454 }
3410 else 3455 else
3411 { 3456 warn_singular_matrix (rcond);
3412 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 3457 }
3413 retval.resize (nc, b_nc, 0.); 3458
3414 3459 volatile double rcond_plus_one = rcond + 1.0;
3415 for (octave_idx_type j = 0; j < b_nc; j++) 3460
3416 { 3461 if (rcond_plus_one == 1.0 || xisnan (rcond))
3417 for (octave_idx_type i = 0; i < nr; i++) 3462 {
3418 cwork[i] = b(i,j); 3463 err = -2;
3419 for (octave_idx_type i = nr; i < nc; i++) 3464
3420 cwork[i] = 0.; 3465 if (sing_handler)
3421 3466 {
3422 for (octave_idx_type k = 0; k < nc; k++) 3467 sing_handler (rcond);
3423 { 3468 mattype.mark_as_rectangular ();
3424 if (cwork[k] != 0.) 3469 }
3425 { 3470 else
3426 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) 3471 warn_singular_matrix (rcond);
3427 { 3472 }
3428 err = -2;
3429 goto triangular_error;
3430 }
3431
3432 Complex tmp = cwork[k] / data (cidx (k));
3433 cwork[k] = tmp;
3434 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3435 {
3436 octave_idx_type iidx = ridx (i);
3437 cwork[iidx] = cwork[iidx] - tmp * data (i);
3438 }
3439 }
3440 }
3441
3442 for (octave_idx_type i = 0; i < nc; i++)
3443 retval.xelem (i, j) = cwork[i];
3444 }
3445
3446 if (calc_cond)
3447 {
3448 // Calculation of 1-norm of inv(*this)
3449 OCTAVE_LOCAL_BUFFER (double, work, nm);
3450 for (octave_idx_type i = 0; i < nm; i++)
3451 work[i] = 0.;
3452
3453 for (octave_idx_type j = 0; j < nr; j++)
3454 {
3455 work[j] = 1.;
3456
3457 for (octave_idx_type k = j; k < nc; k++)
3458 {
3459
3460 if (work[k] != 0.)
3461 {
3462 double tmp = work[k] / data (cidx (k));
3463 work[k] = tmp;
3464 for (octave_idx_type i = cidx (k)+1;
3465 i < cidx (k+1); i++)
3466 {
3467 octave_idx_type iidx = ridx (i);
3468 work[iidx] = work[iidx] - tmp * data (i);
3469 }
3470 }
3471 }
3472 double atmp = 0;
3473 for (octave_idx_type i = j; i < nc; i++)
3474 {
3475 atmp += fabs (work[i]);
3476 work[i] = 0.;
3477 }
3478 if (atmp > ainvnorm)
3479 ainvnorm = atmp;
3480 }
3481 rcond = 1. / ainvnorm / anorm;
3482 }
3483 }
3484
3485 triangular_error:
3486 if (err != 0)
3487 {
3488 if (sing_handler)
3489 {
3490 sing_handler (rcond);
3491 mattype.mark_as_rectangular ();
3492 }
3493 else
3494 warn_singular_matrix (rcond);
3495 }
3496
3497 volatile double rcond_plus_one = rcond + 1.0;
3498
3499 if (rcond_plus_one == 1.0 || xisnan (rcond))
3500 {
3501 err = -2;
3502
3503 if (sing_handler)
3504 {
3505 sing_handler (rcond);
3506 mattype.mark_as_rectangular ();
3507 }
3508 else
3509 warn_singular_matrix (rcond);
3510 }
3511 }
3512 else
3513 (*current_liboctave_error_handler) ("incorrect matrix type");
3514 } 3473 }
3515 3474
3516 return retval; 3475 return retval;
3517 } 3476 }
3518 3477
3530 err = 0; 3489 err = 0;
3531 3490
3532 if (nr != b.rows ()) 3491 if (nr != b.rows ())
3533 (*current_liboctave_error_handler) 3492 (*current_liboctave_error_handler)
3534 ("matrix dimension mismatch solution of linear equations"); 3493 ("matrix dimension mismatch solution of linear equations");
3535 else if (nr == 0 || nc == 0 || b.cols () == 0) 3494
3495 if (nr == 0 || nc == 0 || b.cols () == 0)
3536 retval = SparseComplexMatrix (nc, b.cols ()); 3496 retval = SparseComplexMatrix (nc, b.cols ());
3537 else 3497 else
3538 { 3498 {
3539 // Print spparms("spumoni") info if requested 3499 // Print spparms("spumoni") info if requested
3540 int typ = mattype.type (); 3500 int typ = mattype.type ();
3541 mattype.info (); 3501 mattype.info ();
3542 3502
3543 if (typ == MatrixType::Permuted_Lower || typ == MatrixType::Lower) 3503 if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower)
3544 { 3504 (*current_liboctave_error_handler) ("incorrect matrix type");
3545 double anorm = 0.; 3505
3546 double ainvnorm = 0.; 3506 double anorm = 0.;
3547 rcond = 1.; 3507 double ainvnorm = 0.;
3508 rcond = 1.;
3509
3510 if (calc_cond)
3511 {
3512 // Calculate the 1-norm of matrix for rcond calculation
3513 for (octave_idx_type j = 0; j < nc; j++)
3514 {
3515 double atmp = 0.;
3516 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++)
3517 atmp += fabs (data (i));
3518 if (atmp > anorm)
3519 anorm = atmp;
3520 }
3521 }
3522
3523 octave_idx_type b_nc = b.cols ();
3524 octave_idx_type b_nz = b.nnz ();
3525 retval = SparseComplexMatrix (nc, b_nc, b_nz);
3526 retval.xcidx (0) = 0;
3527 octave_idx_type ii = 0;
3528 octave_idx_type x_nz = b_nz;
3529
3530 if (typ == MatrixType::Permuted_Lower)
3531 {
3532 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3533 octave_idx_type *perm = mattype.triangular_perm ();
3534
3535 for (octave_idx_type j = 0; j < b_nc; j++)
3536 {
3537 for (octave_idx_type i = 0; i < nm; i++)
3538 cwork[i] = 0.;
3539 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3540 cwork[perm[b.ridx (i)]] = b.data (i);
3541
3542 for (octave_idx_type k = 0; k < nc; k++)
3543 {
3544 if (cwork[k] != 0.)
3545 {
3546 octave_idx_type minr = nr;
3547 octave_idx_type mini = 0;
3548
3549 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3550 if (perm[ridx (i)] < minr)
3551 {
3552 minr = perm[ridx (i)];
3553 mini = i;
3554 }
3555
3556 if (minr != k || data (mini) == 0)
3557 {
3558 err = -2;
3559 goto triangular_error;
3560 }
3561
3562 Complex tmp = cwork[k] / data (mini);
3563 cwork[k] = tmp;
3564 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3565 {
3566 if (i == mini)
3567 continue;
3568
3569 octave_idx_type iidx = perm[ridx (i)];
3570 cwork[iidx] = cwork[iidx] - tmp * data (i);
3571 }
3572 }
3573 }
3574
3575 // Count nonzeros in work vector and adjust space in
3576 // retval if needed
3577 octave_idx_type new_nnz = 0;
3578 for (octave_idx_type i = 0; i < nc; i++)
3579 if (cwork[i] != 0.)
3580 new_nnz++;
3581
3582 if (ii + new_nnz > x_nz)
3583 {
3584 // Resize the sparse matrix
3585 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3586 retval.change_capacity (sz);
3587 x_nz = sz;
3588 }
3589
3590 for (octave_idx_type i = 0; i < nc; i++)
3591 if (cwork[i] != 0.)
3592 {
3593 retval.xridx (ii) = i;
3594 retval.xdata (ii++) = cwork[i];
3595 }
3596 retval.xcidx (j+1) = ii;
3597 }
3598
3599 retval.maybe_compress ();
3548 3600
3549 if (calc_cond) 3601 if (calc_cond)
3550 { 3602 {
3551 // Calculate the 1-norm of matrix for rcond calculation 3603 // Calculation of 1-norm of inv(*this)
3552 for (octave_idx_type j = 0; j < nc; j++) 3604 OCTAVE_LOCAL_BUFFER (double, work, nm);
3553 { 3605 for (octave_idx_type i = 0; i < nm; i++)
3554 double atmp = 0.; 3606 work[i] = 0.;
3555 for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) 3607
3556 atmp += fabs (data (i)); 3608 for (octave_idx_type j = 0; j < nr; j++)
3557 if (atmp > anorm) 3609 {
3558 anorm = atmp; 3610 work[j] = 1.;
3559 }
3560 }
3561
3562 octave_idx_type b_nc = b.cols ();
3563 octave_idx_type b_nz = b.nnz ();
3564 retval = SparseComplexMatrix (nc, b_nc, b_nz);
3565 retval.xcidx (0) = 0;
3566 octave_idx_type ii = 0;
3567 octave_idx_type x_nz = b_nz;
3568
3569 if (typ == MatrixType::Permuted_Lower)
3570 {
3571 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3572 octave_idx_type *perm = mattype.triangular_perm ();
3573
3574 for (octave_idx_type j = 0; j < b_nc; j++)
3575 {
3576 for (octave_idx_type i = 0; i < nm; i++)
3577 cwork[i] = 0.;
3578 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3579 cwork[perm[b.ridx (i)]] = b.data (i);
3580 3611
3581 for (octave_idx_type k = 0; k < nc; k++) 3612 for (octave_idx_type k = 0; k < nc; k++)
3582 { 3613 {
3583 if (cwork[k] != 0.) 3614 if (work[k] != 0.)
3584 { 3615 {
3585 octave_idx_type minr = nr; 3616 octave_idx_type minr = nr;
3586 octave_idx_type mini = 0; 3617 octave_idx_type mini = 0;
3587 3618
3588 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 3619 for (octave_idx_type i = cidx (k);
3620 i < cidx (k+1); i++)
3589 if (perm[ridx (i)] < minr) 3621 if (perm[ridx (i)] < minr)
3590 { 3622 {
3591 minr = perm[ridx (i)]; 3623 minr = perm[ridx (i)];
3592 mini = i; 3624 mini = i;
3593 } 3625 }
3594 3626
3595 if (minr != k || data (mini) == 0) 3627 double tmp = work[k] / data (mini);
3596 { 3628 work[k] = tmp;
3597 err = -2; 3629 for (octave_idx_type i = cidx (k);
3598 goto triangular_error; 3630 i < cidx (k+1); i++)
3599 }
3600
3601 Complex tmp = cwork[k] / data (mini);
3602 cwork[k] = tmp;
3603 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
3604 { 3631 {
3605 if (i == mini) 3632 if (i == mini)
3606 continue; 3633 continue;
3607 3634
3608 octave_idx_type iidx = perm[ridx (i)]; 3635 octave_idx_type iidx = perm[ridx (i)];
3609 cwork[iidx] = cwork[iidx] - tmp * data (i); 3636 work[iidx] = work[iidx] - tmp * data (i);
3610 } 3637 }
3611 } 3638 }
3612 } 3639 }
3613 3640
3614 // Count nonzeros in work vector and adjust space in 3641 double atmp = 0;
3615 // retval if needed 3642 for (octave_idx_type i = j; i < nc; i++)
3616 octave_idx_type new_nnz = 0; 3643 {
3617 for (octave_idx_type i = 0; i < nc; i++) 3644 atmp += fabs (work[i]);
3618 if (cwork[i] != 0.) 3645 work[i] = 0.;
3619 new_nnz++; 3646 }
3620 3647 if (atmp > ainvnorm)
3621 if (ii + new_nnz > x_nz) 3648 ainvnorm = atmp;
3622 { 3649 }
3623 // Resize the sparse matrix 3650 rcond = 1. / ainvnorm / anorm;
3624 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; 3651 }
3625 retval.change_capacity (sz); 3652 }
3626 x_nz = sz; 3653 else
3627 } 3654 {
3628 3655 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm);
3629 for (octave_idx_type i = 0; i < nc; i++) 3656
3630 if (cwork[i] != 0.) 3657 for (octave_idx_type j = 0; j < b_nc; j++)
3631 { 3658 {
3632 retval.xridx (ii) = i; 3659 for (octave_idx_type i = 0; i < nm; i++)
3633 retval.xdata (ii++) = cwork[i]; 3660 cwork[i] = 0.;
3634 } 3661 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++)
3635 retval.xcidx (j+1) = ii; 3662 cwork[b.ridx (i)] = b.data (i);
3636 } 3663
3637 3664 for (octave_idx_type k = 0; k < nc; k++)
3638 retval.maybe_compress (); 3665 {
3639 3666 if (cwork[k] != 0.)
3640 if (calc_cond) 3667 {
3641 { 3668 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.)
3642 // Calculation of 1-norm of inv(*this)
3643 OCTAVE_LOCAL_BUFFER (double, work, nm);
3644 for (octave_idx_type i = 0; i < nm; i++)
3645 work[i] = 0.;
3646
3647 for (octave_idx_type j = 0; j < nr; j++)
3648 {
3649 work[j] = 1.;
3650
3651 for (octave_idx_type k = 0; k < nc; k++)
3652 { 3669 {
3653 if (work[k] != 0.) 3670 err = -2;
3671 goto triangular_error;
3672 }
3673
3674 Complex tmp = cwork[k] / data (cidx (k));
3675 cwork[k] = tmp;
3676 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3677 {
3678 octave_idx_type iidx = ridx (i);
3679 cwork[iidx] = cwork[iidx] - tmp * data (i);
3680 }
3681 }
3682 }
3683
3684 // Count nonzeros in work vector and adjust space in
3685 // retval if needed
3686 octave_idx_type new_nnz = 0;
3687 for (octave_idx_type i = 0; i < nc; i++)
3688 if (cwork[i] != 0.)
3689 new_nnz++;
3690
3691 if (ii + new_nnz > x_nz)
3692 {
3693 // Resize the sparse matrix
3694 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3695 retval.change_capacity (sz);
3696 x_nz = sz;
3697 }
3698
3699 for (octave_idx_type i = 0; i < nc; i++)
3700 if (cwork[i] != 0.)
3701 {
3702 retval.xridx (ii) = i;
3703 retval.xdata (ii++) = cwork[i];
3704 }
3705 retval.xcidx (j+1) = ii;
3706 }
3707
3708 retval.maybe_compress ();
3709
3710 if (calc_cond)
3711 {
3712 // Calculation of 1-norm of inv(*this)
3713 OCTAVE_LOCAL_BUFFER (double, work, nm);
3714 for (octave_idx_type i = 0; i < nm; i++)
3715 work[i] = 0.;
3716
3717 for (octave_idx_type j = 0; j < nr; j++)
3718 {
3719 work[j] = 1.;
3720
3721 for (octave_idx_type k = j; k < nc; k++)
3722 {
3723
3724 if (work[k] != 0.)
3725 {
3726 double tmp = work[k] / data (cidx (k));
3727 work[k] = tmp;
3728 for (octave_idx_type i = cidx (k)+1;
3729 i < cidx (k+1); i++)
3654 { 3730 {
3655 octave_idx_type minr = nr; 3731 octave_idx_type iidx = ridx (i);
3656 octave_idx_type mini = 0; 3732 work[iidx] = work[iidx] - tmp * data (i);
3657
3658 for (octave_idx_type i = cidx (k);
3659 i < cidx (k+1); i++)
3660 if (perm[ridx (i)] < minr)
3661 {
3662 minr = perm[ridx (i)];
3663 mini = i;
3664 }
3665
3666 double tmp = work[k] / data (mini);
3667 work[k] = tmp;
3668 for (octave_idx_type i = cidx (k);
3669 i < cidx (k+1); i++)
3670 {
3671 if (i == mini)
3672 continue;
3673
3674 octave_idx_type iidx = perm[ridx (i)];
3675 work[iidx] = work[iidx] - tmp * data (i);
3676 }
3677 } 3733 }
3678 } 3734 }
3679 3735 }
3680 double atmp = 0; 3736 double atmp = 0;
3681 for (octave_idx_type i = j; i < nc; i++) 3737 for (octave_idx_type i = j; i < nc; i++)
3682 { 3738 {
3683 atmp += fabs (work[i]); 3739 atmp += fabs (work[i]);
3684 work[i] = 0.; 3740 work[i] = 0.;
3685 } 3741 }
3686 if (atmp > ainvnorm) 3742 if (atmp > ainvnorm)
3687 ainvnorm = atmp; 3743 ainvnorm = atmp;
3688 } 3744 }
3689 rcond = 1. / ainvnorm / anorm; 3745 rcond = 1. / ainvnorm / anorm;
3690 } 3746 }
3747 }
3748
3749 triangular_error:
3750 if (err != 0)
3751 {
3752 if (sing_handler)
3753 {
3754 sing_handler (rcond);
3755 mattype.mark_as_rectangular ();
3691 } 3756 }
3692 else 3757 else
3693 { 3758 warn_singular_matrix (rcond);
3694 OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 3759 }
3695 3760
3696 for (octave_idx_type j = 0; j < b_nc; j++) 3761 volatile double rcond_plus_one = rcond + 1.0;
3697 { 3762
3698 for (octave_idx_type i = 0; i < nm; i++) 3763 if (rcond_plus_one == 1.0 || xisnan (rcond))
3699 cwork[i] = 0.; 3764 {
3700 for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) 3765 err = -2;
3701 cwork[b.ridx (i)] = b.data (i); 3766
3702 3767 if (sing_handler)
3703 for (octave_idx_type k = 0; k < nc; k++) 3768 {
3704 { 3769 sing_handler (rcond);
3705 if (cwork[k] != 0.) 3770 mattype.mark_as_rectangular ();
3706 { 3771 }
3707 if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) 3772 else
3708 { 3773 warn_singular_matrix (rcond);
3709 err = -2; 3774 }
3710 goto triangular_error;
3711 }
3712
3713 Complex tmp = cwork[k] / data (cidx (k));
3714 cwork[k] = tmp;
3715 for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++)
3716 {
3717 octave_idx_type iidx = ridx (i);
3718 cwork[iidx] = cwork[iidx] - tmp * data (i);
3719 }
3720 }
3721 }
3722
3723 // Count nonzeros in work vector and adjust space in
3724 // retval if needed
3725 octave_idx_type new_nnz = 0;
3726 for (octave_idx_type i = 0; i < nc; i++)
3727 if (cwork[i] != 0.)
3728 new_nnz++;
3729
3730 if (ii + new_nnz > x_nz)
3731 {
3732 // Resize the sparse matrix
3733 octave_idx_type sz = new_nnz * (b_nc - j) + x_nz;
3734 retval.change_capacity (sz);
3735 x_nz = sz;
3736 }
3737
3738 for (octave_idx_type i = 0; i < nc; i++)
3739 if (cwork[i] != 0.)
3740 {
3741 retval.xridx (ii) = i;
3742 retval.xdata (ii++) = cwork[i];
3743 }
3744 retval.xcidx (j+1) = ii;
3745 }
3746
3747 retval.maybe_compress ();
3748
3749 if (calc_cond)
3750 {
3751 // Calculation of 1-norm of inv(*this)
3752 OCTAVE_LOCAL_BUFFER (double, work, nm);
3753 for (octave_idx_type i = 0; i < nm; i++)
3754 work[i] = 0.;
3755
3756 for (octave_idx_type j = 0; j < nr; j++)
3757 {
3758 work[j] = 1.;
3759
3760 for (octave_idx_type k = j; k < nc; k++)
3761 {
3762
3763 if (work[k] != 0.)
3764 {
3765 double tmp = work[k] / data (cidx (k));
3766 work[k] = tmp;
3767 for (octave_idx_type i = cidx (k)+1;
3768 i < cidx (k+1); i++)
3769 {
3770 octave_idx_type iidx = ridx (i);
3771 work[iidx] = work[iidx] - tmp * data (i);
3772 }
3773 }
3774 }
3775 double atmp = 0;
3776 for (octave_idx_type i = j; i < nc; i++)
3777 {
3778 atmp += fabs (work[i]);
3779 work[i] = 0.;
3780 }
3781 if (atmp > ainvnorm)
3782 ainvnorm = atmp;
3783 }
3784 rcond = 1. / ainvnorm / anorm;
3785 }
3786 }
3787
3788 triangular_error:
3789 if (err != 0)
3790 {
3791 if (sing_handler)
3792 {
3793 sing_handler (rcond);
3794 mattype.mark_as_rectangular ();
3795 }
3796 else
3797 warn_singular_matrix (rcond);
3798 }
3799
3800 volatile double rcond_plus_one = rcond + 1.0;
3801
3802 if (rcond_plus_one == 1.0 || xisnan (rcond))
3803 {
3804 err = -2;
3805
3806 if (sing_handler)
3807 {
3808 sing_handler (rcond);
3809 mattype.mark_as_rectangular ();
3810 }
3811 else
3812 warn_singular_matrix (rcond);
3813 }
3814 }
3815 else
3816 (*current_liboctave_error_handler) ("incorrect matrix type");
3817 } 3775 }
3818 3776
3819 return retval; 3777 return retval;
3820 } 3778 }
3821 3779
3832 err = 0; 3790 err = 0;
3833 3791
3834 if (nr != nc || nr != b.rows ()) 3792 if (nr != nc || nr != b.rows ())
3835 (*current_liboctave_error_handler) 3793 (*current_liboctave_error_handler)
3836 ("matrix dimension mismatch solution of linear equations"); 3794 ("matrix dimension mismatch solution of linear equations");
3837 else if (nr == 0 || b.cols () == 0) 3795
3796 if (nr == 0 || b.cols () == 0)
3838 retval = Matrix (nc, b.cols (), 0.0); 3797 retval = Matrix (nc, b.cols (), 0.0);
3839 else if (calc_cond) 3798 else if (calc_cond)
3840 (*current_liboctave_error_handler) 3799 (*current_liboctave_error_handler)
3841 ("calculation of condition number not implemented"); 3800 ("calculation of condition number not implemented");
3842 else 3801 else
3981 err = 0; 3940 err = 0;
3982 3941
3983 if (nr != nc || nr != b.rows ()) 3942 if (nr != nc || nr != b.rows ())
3984 (*current_liboctave_error_handler) 3943 (*current_liboctave_error_handler)
3985 ("matrix dimension mismatch solution of linear equations"); 3944 ("matrix dimension mismatch solution of linear equations");
3986 else if (nr == 0 || b.cols () == 0) 3945
3946 if (nr == 0 || b.cols () == 0)
3987 retval = SparseMatrix (nc, b.cols ()); 3947 retval = SparseMatrix (nc, b.cols ());
3988 else if (calc_cond) 3948 else if (calc_cond)
3989 (*current_liboctave_error_handler) 3949 (*current_liboctave_error_handler)
3990 ("calculation of condition number not implemented"); 3950 ("calculation of condition number not implemented");
3991 else 3951 else
4126 err = 0; 4086 err = 0;
4127 4087
4128 if (nr != nc || nr != b.rows ()) 4088 if (nr != nc || nr != b.rows ())
4129 (*current_liboctave_error_handler) 4089 (*current_liboctave_error_handler)
4130 ("matrix dimension mismatch solution of linear equations"); 4090 ("matrix dimension mismatch solution of linear equations");
4131 else if (nr == 0 || b.cols () == 0) 4091
4092 if (nr == 0 || b.cols () == 0)
4132 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 4093 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
4133 else if (calc_cond) 4094 else if (calc_cond)
4134 (*current_liboctave_error_handler) 4095 (*current_liboctave_error_handler)
4135 ("calculation of condition number not implemented"); 4096 ("calculation of condition number not implemented");
4136 else 4097 else
4277 err = 0; 4238 err = 0;
4278 4239
4279 if (nr != nc || nr != b.rows ()) 4240 if (nr != nc || nr != b.rows ())
4280 (*current_liboctave_error_handler) 4241 (*current_liboctave_error_handler)
4281 ("matrix dimension mismatch solution of linear equations"); 4242 ("matrix dimension mismatch solution of linear equations");
4282 else if (nr == 0 || b.cols () == 0) 4243
4244 if (nr == 0 || b.cols () == 0)
4283 retval = SparseComplexMatrix (nc, b.cols ()); 4245 retval = SparseComplexMatrix (nc, b.cols ());
4284 else if (calc_cond) 4246 else if (calc_cond)
4285 (*current_liboctave_error_handler) 4247 (*current_liboctave_error_handler)
4286 ("calculation of condition number not implemented"); 4248 ("calculation of condition number not implemented");
4287 else 4249 else
4382 Bx, b_nr, err 4344 Bx, b_nr, err
4383 F77_CHAR_ARG_LEN (1))); 4345 F77_CHAR_ARG_LEN (1)));
4384 4346
4385 if (err != 0) 4347 if (err != 0)
4386 { 4348 {
4349 // FIXME: Should this be a warning?
4387 (*current_liboctave_error_handler) 4350 (*current_liboctave_error_handler)
4388 ("SparseMatrix::solve solve failed"); 4351 ("SparseMatrix::solve solve failed");
4389 4352
4390 err = -1; 4353 err = -1;
4391 break; 4354 break;
4397 Bz, b_nr, err 4360 Bz, b_nr, err
4398 F77_CHAR_ARG_LEN (1))); 4361 F77_CHAR_ARG_LEN (1)));
4399 4362
4400 if (err != 0) 4363 if (err != 0)
4401 { 4364 {
4365 // FIXME: Should this be a warning?
4402 (*current_liboctave_error_handler) 4366 (*current_liboctave_error_handler)
4403 ("SparseMatrix::solve solve failed"); 4367 ("SparseMatrix::solve solve failed");
4404 4368
4405 err = -1; 4369 err = -1;
4406 break; 4370 break;
4455 err = 0; 4419 err = 0;
4456 4420
4457 if (nr != nc || nr != b.rows ()) 4421 if (nr != nc || nr != b.rows ())
4458 (*current_liboctave_error_handler) 4422 (*current_liboctave_error_handler)
4459 ("matrix dimension mismatch solution of linear equations"); 4423 ("matrix dimension mismatch solution of linear equations");
4460 else if (nr == 0 || b.cols () == 0) 4424
4425 if (nr == 0 || b.cols () == 0)
4461 retval = Matrix (nc, b.cols (), 0.0); 4426 retval = Matrix (nc, b.cols (), 0.0);
4462 else 4427 else
4463 { 4428 {
4464 // Print spparms("spumoni") info if requested 4429 // Print spparms("spumoni") info if requested
4465 volatile int typ = mattype.type (); 4430 volatile int typ = mattype.type ();
4557 ldm, result, b.rows (), err 4522 ldm, result, b.rows (), err
4558 F77_CHAR_ARG_LEN (1))); 4523 F77_CHAR_ARG_LEN (1)));
4559 4524
4560 if (err != 0) 4525 if (err != 0)
4561 { 4526 {
4527 // FIXME: Should this be a warning?
4562 (*current_liboctave_error_handler) 4528 (*current_liboctave_error_handler)
4563 ("SparseMatrix::solve solve failed"); 4529 ("SparseMatrix::solve solve failed");
4564 err = -1; 4530 err = -1;
4565 } 4531 }
4566 } 4532 }
4698 err = 0; 4664 err = 0;
4699 4665
4700 if (nr != nc || nr != b.rows ()) 4666 if (nr != nc || nr != b.rows ())
4701 (*current_liboctave_error_handler) 4667 (*current_liboctave_error_handler)
4702 ("matrix dimension mismatch solution of linear equations"); 4668 ("matrix dimension mismatch solution of linear equations");
4703 else if (nr == 0 || b.cols () == 0) 4669
4670 if (nr == 0 || b.cols () == 0)
4704 retval = SparseMatrix (nc, b.cols ()); 4671 retval = SparseMatrix (nc, b.cols ());
4705 else 4672 else
4706 { 4673 {
4707 // Print spparms("spumoni") info if requested 4674 // Print spparms("spumoni") info if requested
4708 volatile int typ = mattype.type (); 4675 volatile int typ = mattype.type ();
4810 ldm, Bx, b_nr, err 4777 ldm, Bx, b_nr, err
4811 F77_CHAR_ARG_LEN (1))); 4778 F77_CHAR_ARG_LEN (1)));
4812 4779
4813 if (err != 0) 4780 if (err != 0)
4814 { 4781 {
4782 // FIXME: Should this be a warning?
4815 (*current_liboctave_error_handler) 4783 (*current_liboctave_error_handler)
4816 ("SparseMatrix::solve solve failed"); 4784 ("SparseMatrix::solve solve failed");
4817 err = -1; 4785 err = -1;
4818 break; 4786 break;
4819 } 4787 }
5010 err = 0; 4978 err = 0;
5011 4979
5012 if (nr != nc || nr != b.rows ()) 4980 if (nr != nc || nr != b.rows ())
5013 (*current_liboctave_error_handler) 4981 (*current_liboctave_error_handler)
5014 ("matrix dimension mismatch solution of linear equations"); 4982 ("matrix dimension mismatch solution of linear equations");
5015 else if (nr == 0 || b.cols () == 0) 4983
4984 if (nr == 0 || b.cols () == 0)
5016 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 4985 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
5017 else 4986 else
5018 { 4987 {
5019 // Print spparms("spumoni") info if requested 4988 // Print spparms("spumoni") info if requested
5020 volatile int typ = mattype.type (); 4989 volatile int typ = mattype.type ();
5125 ldm, Bx, b_nr, err 5094 ldm, Bx, b_nr, err
5126 F77_CHAR_ARG_LEN (1))); 5095 F77_CHAR_ARG_LEN (1)));
5127 5096
5128 if (err != 0) 5097 if (err != 0)
5129 { 5098 {
5099 // FIXME: Should this be a warning?
5130 (*current_liboctave_error_handler) 5100 (*current_liboctave_error_handler)
5131 ("SparseMatrix::solve solve failed"); 5101 ("SparseMatrix::solve solve failed");
5132 err = -1; 5102 err = -1;
5133 break; 5103 break;
5134 } 5104 }
5139 ldm, Bz, b.rows (), err 5109 ldm, Bz, b.rows (), err
5140 F77_CHAR_ARG_LEN (1))); 5110 F77_CHAR_ARG_LEN (1)));
5141 5111
5142 if (err != 0) 5112 if (err != 0)
5143 { 5113 {
5114 // FIXME: Should this be a warning?
5144 (*current_liboctave_error_handler) 5115 (*current_liboctave_error_handler)
5145 ("SparseMatrix::solve solve failed"); 5116 ("SparseMatrix::solve solve failed");
5146 err = -1; 5117 err = -1;
5147 break; 5118 break;
5148 } 5119 }
5303 err = 0; 5274 err = 0;
5304 5275
5305 if (nr != nc || nr != b.rows ()) 5276 if (nr != nc || nr != b.rows ())
5306 (*current_liboctave_error_handler) 5277 (*current_liboctave_error_handler)
5307 ("matrix dimension mismatch solution of linear equations"); 5278 ("matrix dimension mismatch solution of linear equations");
5308 else if (nr == 0 || b.cols () == 0) 5279
5280 if (nr == 0 || b.cols () == 0)
5309 retval = SparseComplexMatrix (nc, b.cols ()); 5281 retval = SparseComplexMatrix (nc, b.cols ());
5310 else 5282 else
5311 { 5283 {
5312 // Print spparms("spumoni") info if requested 5284 // Print spparms("spumoni") info if requested
5313 volatile int typ = mattype.type (); 5285 volatile int typ = mattype.type ();
5424 ldm, Bx, b_nr, err 5396 ldm, Bx, b_nr, err
5425 F77_CHAR_ARG_LEN (1))); 5397 F77_CHAR_ARG_LEN (1)));
5426 5398
5427 if (err != 0) 5399 if (err != 0)
5428 { 5400 {
5401 // FIXME: Should this be a warning?
5429 (*current_liboctave_error_handler) 5402 (*current_liboctave_error_handler)
5430 ("SparseMatrix::solve solve failed"); 5403 ("SparseMatrix::solve solve failed");
5431 err = -1; 5404 err = -1;
5432 break; 5405 break;
5433 } 5406 }
5438 ldm, Bz, b_nr, err 5411 ldm, Bz, b_nr, err
5439 F77_CHAR_ARG_LEN (1))); 5412 F77_CHAR_ARG_LEN (1)));
5440 5413
5441 if (err != 0) 5414 if (err != 0)
5442 { 5415 {
5416 // FIXME: Should this be a warning?
5443 (*current_liboctave_error_handler) 5417 (*current_liboctave_error_handler)
5444 ("SparseMatrix::solve solve failed"); 5418 ("SparseMatrix::solve solve failed");
5445 5419
5446 err = -1; 5420 err = -1;
5447 break; 5421 break;
5692 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0, 5666 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
5693 &Symbolic, control, info); 5667 &Symbolic, control, info);
5694 5668
5695 if (status < 0) 5669 if (status < 0)
5696 { 5670 {
5671 UMFPACK_DNAME (report_status) (control, status);
5672 UMFPACK_DNAME (report_info) (control, info);
5673
5674 UMFPACK_DNAME (free_symbolic) (&Symbolic);
5675
5676 // FIXME: Should this be a warning?
5697 (*current_liboctave_error_handler) 5677 (*current_liboctave_error_handler)
5698 ("SparseMatrix::solve symbolic factorization failed"); 5678 ("SparseMatrix::solve symbolic factorization failed");
5699 err = -1; 5679 err = -1;
5700
5701 UMFPACK_DNAME (report_status) (control, status);
5702 UMFPACK_DNAME (report_info) (control, info);
5703
5704 UMFPACK_DNAME (free_symbolic) (&Symbolic);
5705 } 5680 }
5706 else 5681 else
5707 { 5682 {
5708 UMFPACK_DNAME (report_symbolic) (Symbolic, control); 5683 UMFPACK_DNAME (report_symbolic) (Symbolic, control);
5709 5684
5729 else 5704 else
5730 warn_singular_matrix (rcond); 5705 warn_singular_matrix (rcond);
5731 } 5706 }
5732 else if (status < 0) 5707 else if (status < 0)
5733 { 5708 {
5709 UMFPACK_DNAME (report_status) (control, status);
5710 UMFPACK_DNAME (report_info) (control, info);
5711
5712 // FIXME: Should this be a warning?
5734 (*current_liboctave_error_handler) 5713 (*current_liboctave_error_handler)
5735 ("SparseMatrix::solve numeric factorization failed"); 5714 ("SparseMatrix::solve numeric factorization failed");
5736
5737 UMFPACK_DNAME (report_status) (control, status);
5738 UMFPACK_DNAME (report_info) (control, info);
5739 5715
5740 err = -1; 5716 err = -1;
5741 } 5717 }
5742 else 5718 else
5743 { 5719 {
5770 err = 0; 5746 err = 0;
5771 5747
5772 if (nr != nc || nr != b.rows ()) 5748 if (nr != nc || nr != b.rows ())
5773 (*current_liboctave_error_handler) 5749 (*current_liboctave_error_handler)
5774 ("matrix dimension mismatch solution of linear equations"); 5750 ("matrix dimension mismatch solution of linear equations");
5775 else if (nr == 0 || b.cols () == 0) 5751
5752 if (nr == 0 || b.cols () == 0)
5776 retval = Matrix (nc, b.cols (), 0.0); 5753 retval = Matrix (nc, b.cols (), 0.0);
5777 else 5754 else
5778 { 5755 {
5779 // Print spparms("spumoni") info if requested 5756 // Print spparms("spumoni") info if requested
5780 volatile int typ = mattype.type (); 5757 volatile int typ = mattype.type ();
5943 Ai, Ax, &result[iidx], 5920 Ai, Ax, &result[iidx],
5944 &Bx[iidx], Numeric, control, 5921 &Bx[iidx], Numeric, control,
5945 info); 5922 info);
5946 if (status < 0) 5923 if (status < 0)
5947 { 5924 {
5925 UMFPACK_DNAME (report_status) (control, status);
5926
5927 // FIXME: Should this be a warning?
5948 (*current_liboctave_error_handler) 5928 (*current_liboctave_error_handler)
5949 ("SparseMatrix::solve solve failed"); 5929 ("SparseMatrix::solve solve failed");
5950 5930
5951 UMFPACK_DNAME (report_status) (control, status);
5952
5953 err = -1; 5931 err = -1;
5954
5955 break; 5932 break;
5956 } 5933 }
5957 } 5934 }
5958 5935
5959 UMFPACK_DNAME (report_info) (control, info); 5936 UMFPACK_DNAME (report_info) (control, info);
5989 err = 0; 5966 err = 0;
5990 5967
5991 if (nr != nc || nr != b.rows ()) 5968 if (nr != nc || nr != b.rows ())
5992 (*current_liboctave_error_handler) 5969 (*current_liboctave_error_handler)
5993 ("matrix dimension mismatch solution of linear equations"); 5970 ("matrix dimension mismatch solution of linear equations");
5994 else if (nr == 0 || b.cols () == 0) 5971
5972 if (nr == 0 || b.cols () == 0)
5995 retval = SparseMatrix (nc, b.cols ()); 5973 retval = SparseMatrix (nc, b.cols ());
5996 else 5974 else
5997 { 5975 {
5998 // Print spparms("spumoni") info if requested 5976 // Print spparms("spumoni") info if requested
5999 volatile int typ = mattype.type (); 5977 volatile int typ = mattype.type ();
6186 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap, 6164 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
6187 Ai, Ax, Xx, Bx, Numeric, 6165 Ai, Ax, Xx, Bx, Numeric,
6188 control, info); 6166 control, info);
6189 if (status < 0) 6167 if (status < 0)
6190 { 6168 {
6169 UMFPACK_DNAME (report_status) (control, status);
6170
6171 // FIXME: Should this be a warning?
6191 (*current_liboctave_error_handler) 6172 (*current_liboctave_error_handler)
6192 ("SparseMatrix::solve solve failed"); 6173 ("SparseMatrix::solve solve failed");
6193 6174
6194 UMFPACK_DNAME (report_status) (control, status);
6195
6196 err = -1; 6175 err = -1;
6197
6198 break; 6176 break;
6199 } 6177 }
6200 6178
6201 for (octave_idx_type i = 0; i < b_nr; i++) 6179 for (octave_idx_type i = 0; i < b_nr; i++)
6202 { 6180 {
6253 err = 0; 6231 err = 0;
6254 6232
6255 if (nr != nc || nr != b.rows ()) 6233 if (nr != nc || nr != b.rows ())
6256 (*current_liboctave_error_handler) 6234 (*current_liboctave_error_handler)
6257 ("matrix dimension mismatch solution of linear equations"); 6235 ("matrix dimension mismatch solution of linear equations");
6258 else if (nr == 0 || b.cols () == 0) 6236
6237 if (nr == 0 || b.cols () == 0)
6259 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); 6238 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
6260 else 6239 else
6261 { 6240 {
6262 // Print spparms("spumoni") info if requested 6241 // Print spparms("spumoni") info if requested
6263 volatile int typ = mattype.type (); 6242 volatile int typ = mattype.type ();
6440 Ap, Ai, Ax, Xz, Bz, 6419 Ap, Ai, Ax, Xz, Bz,
6441 Numeric, control, info); 6420 Numeric, control, info);
6442 6421
6443 if (status < 0 || status2 < 0) 6422 if (status < 0 || status2 < 0)
6444 { 6423 {
6424 UMFPACK_DNAME (report_status) (control, status);
6425
6426 // FIXME: Should this be a warning?
6445 (*current_liboctave_error_handler) 6427 (*current_liboctave_error_handler)
6446 ("SparseMatrix::solve solve failed"); 6428 ("SparseMatrix::solve solve failed");
6447 6429
6448 UMFPACK_DNAME (report_status) (control, status);
6449
6450 err = -1; 6430 err = -1;
6451
6452 break; 6431 break;
6453 } 6432 }
6454 6433
6455 for (octave_idx_type i = 0; i < b_nr; i++) 6434 for (octave_idx_type i = 0; i < b_nr; i++)
6456 retval(i, j) = Complex (Xx[i], Xz[i]); 6435 retval(i, j) = Complex (Xx[i], Xz[i]);
6489 err = 0; 6468 err = 0;
6490 6469
6491 if (nr != nc || nr != b.rows ()) 6470 if (nr != nc || nr != b.rows ())
6492 (*current_liboctave_error_handler) 6471 (*current_liboctave_error_handler)
6493 ("matrix dimension mismatch solution of linear equations"); 6472 ("matrix dimension mismatch solution of linear equations");
6494 else if (nr == 0 || b.cols () == 0) 6473
6474 if (nr == 0 || b.cols () == 0)
6495 retval = SparseComplexMatrix (nc, b.cols ()); 6475 retval = SparseComplexMatrix (nc, b.cols ());
6496 else 6476 else
6497 { 6477 {
6498 // Print spparms("spumoni") info if requested 6478 // Print spparms("spumoni") info if requested
6499 volatile int typ = mattype.type (); 6479 volatile int typ = mattype.type ();
6697 Ap, Ai, Ax, Xz, Bz, 6677 Ap, Ai, Ax, Xz, Bz,
6698 Numeric, control, info); 6678 Numeric, control, info);
6699 6679
6700 if (status < 0 || status2 < 0) 6680 if (status < 0 || status2 < 0)
6701 { 6681 {
6682 UMFPACK_DNAME (report_status) (control, status);
6683
6684 // FIXME: Should this be a warning?
6702 (*current_liboctave_error_handler) 6685 (*current_liboctave_error_handler)
6703 ("SparseMatrix::solve solve failed"); 6686 ("SparseMatrix::solve solve failed");
6704 6687
6705 UMFPACK_DNAME (report_status) (control, status);
6706
6707 err = -1; 6688 err = -1;
6708
6709 break; 6689 break;
6710 } 6690 }
6711 6691
6712 for (octave_idx_type i = 0; i < b_nr; i++) 6692 for (octave_idx_type i = 0; i < b_nr; i++)
6713 { 6693 {
6797 || typ == MatrixType::Tridiagonal_Hermitian) 6777 || typ == MatrixType::Tridiagonal_Hermitian)
6798 retval = trisolve (mattype, b, err, rcond, sing_handler, false); 6778 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6799 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 6779 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6800 retval = fsolve (mattype, b, err, rcond, sing_handler, true); 6780 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6801 else if (typ != MatrixType::Rectangular) 6781 else if (typ != MatrixType::Rectangular)
6802 { 6782 (*current_liboctave_error_handler) ("unknown matrix type");
6803 (*current_liboctave_error_handler) ("unknown matrix type");
6804 return Matrix ();
6805 }
6806 6783
6807 // Rectangular or one of the above solvers flags a singular matrix 6784 // Rectangular or one of the above solvers flags a singular matrix
6808 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) 6785 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6809 { 6786 {
6810 rcond = 1.; 6787 rcond = 1.;
6865 || typ == MatrixType::Tridiagonal_Hermitian) 6842 || typ == MatrixType::Tridiagonal_Hermitian)
6866 retval = trisolve (mattype, b, err, rcond, sing_handler, false); 6843 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6867 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 6844 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6868 retval = fsolve (mattype, b, err, rcond, sing_handler, true); 6845 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6869 else if (typ != MatrixType::Rectangular) 6846 else if (typ != MatrixType::Rectangular)
6870 { 6847 (*current_liboctave_error_handler) ("unknown matrix type");
6871 (*current_liboctave_error_handler) ("unknown matrix type");
6872 return SparseMatrix ();
6873 }
6874 6848
6875 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) 6849 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6876 { 6850 {
6877 rcond = 1.; 6851 rcond = 1.;
6878 #ifdef USE_QRSOLVE 6852 #ifdef USE_QRSOLVE
6933 || typ == MatrixType::Tridiagonal_Hermitian) 6907 || typ == MatrixType::Tridiagonal_Hermitian)
6934 retval = trisolve (mattype, b, err, rcond, sing_handler, false); 6908 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
6935 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 6909 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
6936 retval = fsolve (mattype, b, err, rcond, sing_handler, true); 6910 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
6937 else if (typ != MatrixType::Rectangular) 6911 else if (typ != MatrixType::Rectangular)
6938 { 6912 (*current_liboctave_error_handler) ("unknown matrix type");
6939 (*current_liboctave_error_handler) ("unknown matrix type");
6940 return ComplexMatrix ();
6941 }
6942 6913
6943 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) 6914 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
6944 { 6915 {
6945 rcond = 1.; 6916 rcond = 1.;
6946 #ifdef USE_QRSOLVE 6917 #ifdef USE_QRSOLVE
7001 || typ == MatrixType::Tridiagonal_Hermitian) 6972 || typ == MatrixType::Tridiagonal_Hermitian)
7002 retval = trisolve (mattype, b, err, rcond, sing_handler, false); 6973 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
7003 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) 6974 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
7004 retval = fsolve (mattype, b, err, rcond, sing_handler, true); 6975 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
7005 else if (typ != MatrixType::Rectangular) 6976 else if (typ != MatrixType::Rectangular)
7006 { 6977 (*current_liboctave_error_handler) ("unknown matrix type");
7007 (*current_liboctave_error_handler) ("unknown matrix type");
7008 return SparseComplexMatrix ();
7009 }
7010 6978
7011 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) 6979 if (singular_fallback && mattype.type (false) == MatrixType::Rectangular)
7012 { 6980 {
7013 rcond = 1.; 6981 rcond = 1.;
7014 #ifdef USE_QRSOLVE 6982 #ifdef USE_QRSOLVE