Mercurial > hg > octave-nkf

comparison liboctave/array/CSparse.cc @ 17769:49a5a4be04a1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use GNU style coding conventions for code in liboctave/ * liboctave/array/Array-C.cc, liboctave/array/Array-b.cc, liboctave/array/Array-ch.cc, liboctave/array/Array-d.cc, liboctave/array/Array-f.cc, liboctave/array/Array-fC.cc, liboctave/array/Array-util.cc, liboctave/array/Array-util.h, liboctave/array/Array.cc, liboctave/array/Array.h, liboctave/array/Array3.h, liboctave/array/CColVector.cc, liboctave/array/CColVector.h, liboctave/array/CDiagMatrix.cc, liboctave/array/CDiagMatrix.h, liboctave/array/CMatrix.cc, liboctave/array/CMatrix.h, liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h, liboctave/array/CRowVector.cc, liboctave/array/CRowVector.h, liboctave/array/CSparse.cc, liboctave/array/CSparse.h, liboctave/array/DiagArray2.h, liboctave/array/MArray.cc, liboctave/array/MArray.h, liboctave/array/MDiagArray2.cc, liboctave/array/MDiagArray2.h, liboctave/array/MSparse.cc, liboctave/array/MSparse.h, liboctave/array/MatrixType.cc, liboctave/array/MatrixType.h, liboctave/array/PermMatrix.h, liboctave/array/Range.cc, liboctave/array/Range.h, liboctave/array/Sparse.cc, liboctave/array/Sparse.h, liboctave/array/boolMatrix.cc, liboctave/array/boolMatrix.h, liboctave/array/boolNDArray.cc, liboctave/array/boolNDArray.h, liboctave/array/boolSparse.cc, liboctave/array/boolSparse.h, liboctave/array/chMatrix.cc, liboctave/array/chMatrix.h, liboctave/array/chNDArray.cc, liboctave/array/chNDArray.h, liboctave/array/dColVector.h, liboctave/array/dDiagMatrix.cc, liboctave/array/dDiagMatrix.h, liboctave/array/dMatrix.cc, liboctave/array/dMatrix.h, liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h, liboctave/array/dRowVector.h, liboctave/array/dSparse.cc, liboctave/array/dSparse.h, liboctave/array/dim-vector.cc, liboctave/array/dim-vector.h, liboctave/array/fCColVector.cc, liboctave/array/fCColVector.h, liboctave/array/fCDiagMatrix.cc, liboctave/array/fCDiagMatrix.h, liboctave/array/fCMatrix.cc, liboctave/array/fCMatrix.h, liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h, liboctave/array/fCRowVector.cc, liboctave/array/fCRowVector.h, liboctave/array/fColVector.h, liboctave/array/fDiagMatrix.cc, liboctave/array/fDiagMatrix.h, liboctave/array/fMatrix.cc, liboctave/array/fMatrix.h, liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h, liboctave/array/fRowVector.h, liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h, liboctave/array/intNDArray.cc, liboctave/array/intNDArray.h, liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/quit.cc, liboctave/numeric/CmplxCHOL.cc, liboctave/numeric/CmplxCHOL.h, liboctave/numeric/CmplxGEPBAL.cc, liboctave/numeric/CmplxGEPBAL.h, liboctave/numeric/CmplxHESS.h, liboctave/numeric/CmplxLU.cc, liboctave/numeric/CmplxLU.h, liboctave/numeric/CmplxQR.cc, liboctave/numeric/CmplxQRP.cc, liboctave/numeric/CmplxQRP.h, liboctave/numeric/CmplxSCHUR.h, liboctave/numeric/CmplxSVD.cc, liboctave/numeric/CmplxSVD.h, liboctave/numeric/CollocWt.h, liboctave/numeric/DAE.h, liboctave/numeric/DAEFunc.h, liboctave/numeric/DAERT.h, liboctave/numeric/DAERTFunc.h, liboctave/numeric/DASPK.cc, liboctave/numeric/DASRT.cc, liboctave/numeric/DASRT.h, liboctave/numeric/DASSL.cc, liboctave/numeric/DET.h, liboctave/numeric/EIG.cc, liboctave/numeric/EIG.h, liboctave/numeric/LSODE.cc, liboctave/numeric/ODE.h, liboctave/numeric/ODEFunc.h, liboctave/numeric/ODES.h, liboctave/numeric/ODESFunc.h, liboctave/numeric/Quad.cc, liboctave/numeric/Quad.h, liboctave/numeric/SparseCmplxCHOL.h, liboctave/numeric/SparseCmplxLU.cc, liboctave/numeric/SparseCmplxLU.h, liboctave/numeric/SparseCmplxQR.cc, liboctave/numeric/SparseCmplxQR.h, liboctave/numeric/SparseQR.cc, liboctave/numeric/SparseQR.h, liboctave/numeric/SparsedbleCHOL.h, liboctave/numeric/SparsedbleLU.cc, liboctave/numeric/SparsedbleLU.h, liboctave/numeric/base-aepbal.h, liboctave/numeric/base-dae.h, liboctave/numeric/base-de.h, liboctave/numeric/base-lu.cc, liboctave/numeric/base-lu.h, liboctave/numeric/base-min.h, liboctave/numeric/base-qr.h, liboctave/numeric/bsxfun.h, liboctave/numeric/dbleCHOL.cc, liboctave/numeric/dbleCHOL.h, liboctave/numeric/dbleGEPBAL.h, liboctave/numeric/dbleHESS.h, liboctave/numeric/dbleLU.cc, liboctave/numeric/dbleLU.h, liboctave/numeric/dbleQR.cc, liboctave/numeric/dbleQRP.cc, liboctave/numeric/dbleQRP.h, liboctave/numeric/dbleSCHUR.cc, liboctave/numeric/dbleSCHUR.h, liboctave/numeric/dbleSVD.cc, liboctave/numeric/dbleSVD.h, liboctave/numeric/eigs-base.cc, liboctave/numeric/fCmplxAEPBAL.cc, liboctave/numeric/fCmplxAEPBAL.h, liboctave/numeric/fCmplxCHOL.cc, liboctave/numeric/fCmplxCHOL.h, liboctave/numeric/fCmplxGEPBAL.cc, liboctave/numeric/fCmplxGEPBAL.h, liboctave/numeric/fCmplxHESS.h, liboctave/numeric/fCmplxLU.cc, liboctave/numeric/fCmplxLU.h, liboctave/numeric/fCmplxQR.cc, liboctave/numeric/fCmplxQR.h, liboctave/numeric/fCmplxQRP.cc, liboctave/numeric/fCmplxQRP.h, liboctave/numeric/fCmplxSCHUR.cc, liboctave/numeric/fCmplxSCHUR.h, liboctave/numeric/fCmplxSVD.h, liboctave/numeric/fEIG.cc, liboctave/numeric/fEIG.h, liboctave/numeric/floatCHOL.cc, liboctave/numeric/floatCHOL.h, liboctave/numeric/floatGEPBAL.cc, liboctave/numeric/floatGEPBAL.h, liboctave/numeric/floatHESS.h, liboctave/numeric/floatLU.cc, liboctave/numeric/floatLU.h, liboctave/numeric/floatQR.cc, liboctave/numeric/floatQRP.cc, liboctave/numeric/floatQRP.h, liboctave/numeric/floatSCHUR.cc, liboctave/numeric/floatSCHUR.h, liboctave/numeric/floatSVD.cc, liboctave/numeric/floatSVD.h, liboctave/numeric/lo-mappers.cc, liboctave/numeric/lo-mappers.h, liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h, liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-fftw.cc, liboctave/numeric/oct-fftw.h, liboctave/numeric/oct-norm.cc, liboctave/numeric/oct-rand.cc, liboctave/numeric/oct-rand.h, liboctave/numeric/randgamma.c, liboctave/numeric/randgamma.h, liboctave/numeric/randmtzig.c, liboctave/numeric/randpoisson.c, liboctave/numeric/randpoisson.h, liboctave/numeric/sparse-base-chol.h, liboctave/numeric/sparse-base-lu.h, liboctave/numeric/sparse-dmsolve.cc, liboctave/operators/Sparse-diag-op-defs.h, liboctave/operators/Sparse-op-defs.h, liboctave/operators/mx-inlines.cc, liboctave/system/dir-ops.h, liboctave/system/file-ops.cc, liboctave/system/file-stat.cc, liboctave/system/file-stat.h, liboctave/system/lo-sysdep.cc, liboctave/system/lo-sysdep.h, liboctave/system/mach-info.cc, liboctave/system/mach-info.h, liboctave/system/oct-env.cc, liboctave/system/oct-group.cc, liboctave/system/oct-syscalls.cc, liboctave/system/oct-syscalls.h, liboctave/system/oct-time.h, liboctave/system/tempname.c, liboctave/util/action-container.h, liboctave/util/base-list.h, liboctave/util/cmd-edit.cc, liboctave/util/cmd-edit.h, liboctave/util/cmd-hist.cc, liboctave/util/cmd-hist.h, liboctave/util/data-conv.cc, liboctave/util/data-conv.h, liboctave/util/kpse.cc, liboctave/util/lo-array-gripes.cc, liboctave/util/lo-cieee.c, liboctave/util/lo-regexp.cc, liboctave/util/lo-utils.cc, liboctave/util/oct-alloc.cc, liboctave/util/oct-base64.cc, liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h, liboctave/util/oct-glob.cc, liboctave/util/oct-inttypes.cc, liboctave/util/oct-inttypes.h, liboctave/util/oct-locbuf.cc, liboctave/util/oct-locbuf.h, liboctave/util/oct-mem.h, liboctave/util/oct-mutex.cc, liboctave/util/oct-refcount.h, liboctave/util/oct-shlib.cc, liboctave/util/oct-shlib.h, liboctave/util/oct-sort.cc, liboctave/util/oct-sort.h, liboctave/util/pathsearch.cc, liboctave/util/pathsearch.h, liboctave/util/sparse-util.cc, liboctave/util/str-vec.cc, liboctave/util/str-vec.h, liboctave/util/unwind-prot.h, liboctave/util/url-transfer.cc, liboctave/util/url-transfer.h: Use GNU style coding conventions.
author Rik <rik@octave.org>
date Sat, 26 Oct 2013 18:57:05 -0700
parents d63878346099
children fffd0c0ca2dc
comparison
equal deleted inserted replaced
17768:271c0cce0f64 17769:49a5a4be04a1
187 if (nr != nr_a || nc != nc_a || nz != nz_a) 187 if (nr != nr_a || nc != nc_a || nz != nz_a)
188 return false; 188 return false;
189 189
190 for (octave_idx_type i = 0; i < nc + 1; i++) 190 for (octave_idx_type i = 0; i < nc + 1; i++)
191 if (cidx (i) != a.cidx (i)) 191 if (cidx (i) != a.cidx (i))
192 return false; 192 return false;
193 193
194 for (octave_idx_type i = 0; i < nz; i++) 194 for (octave_idx_type i = 0; i < nz; i++)
195 if (data (i) != a.data (i) || ridx (i) != a.ridx (i)) 195 if (data (i) != a.data (i) || ridx (i) != a.ridx (i))
196 return false; 196 return false;
197 197
349 349
350 for (octave_idx_type j = 0; j < nc; j++) 350 for (octave_idx_type j = 0; j < nc; j++)
351 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 351 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
352 if (found [ridx (i)] == -j) 352 if (found [ridx (i)] == -j)
353 found [ridx (i)] = -j - 1; 353 found [ridx (i)] = -j - 1;
354 354
355 for (octave_idx_type i = 0; i < nr; i++) 355 for (octave_idx_type i = 0; i < nr; i++)
356 if (found [i] > -nc && found [i] < 0) 356 if (found [i] > -nc && found [i] < 0)
357 idx_arg.elem (i) = -found [i]; 357 idx_arg.elem (i) = -found [i];
358 358
359 for (octave_idx_type j = 0; j < nc; j++) 359 for (octave_idx_type j = 0; j < nc; j++)
506 506
507 for (octave_idx_type j = 0; j < nc; j++) 507 for (octave_idx_type j = 0; j < nc; j++)
508 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) 508 for (octave_idx_type i = cidx(j); i < cidx(j+1); i++)
509 if (found [ridx (i)] == -j) 509 if (found [ridx (i)] == -j)
510 found [ridx (i)] = -j - 1; 510 found [ridx (i)] = -j - 1;
511 511
512 for (octave_idx_type i = 0; i < nr; i++) 512 for (octave_idx_type i = 0; i < nr; i++)
513 if (found [i] > -nc && found [i] < 0) 513 if (found [i] > -nc && found [i] < 0)
514 idx_arg.elem (i) = -found [i]; 514 idx_arg.elem (i) = -found [i];
515 515
516 for (octave_idx_type j = 0; j < nc; j++) 516 for (octave_idx_type j = 0; j < nc; j++)
561 return result; 561 return result;
562 } 562 }
563 563
564 /* 564 /*
565 565
566 %!assert (max (max (speye (65536) * 1i)), sparse (1i)) 566 %!assert (max (max (speye (65536) * 1i)), sparse (1i))
567 %!assert (min (min (speye (65536) * 1i)), sparse (0)) 567 %!assert (min (min (speye (65536) * 1i)), sparse (0))
568 %!assert (size (max (sparse (8, 0), [], 1)), [1, 0]) 568 %!assert (size (max (sparse (8, 0), [], 1)), [1, 0])
569 %!assert (size (max (sparse (8, 0), [], 2)), [8, 0]) 569 %!assert (size (max (sparse (8, 0), [], 2)), [8, 0])
570 %!assert (size (max (sparse (0, 8), [], 1)), [0, 8]) 570 %!assert (size (max (sparse (0, 8), [], 1)), [0, 8])
571 %!assert (size (max (sparse (0, 8), [], 2)), [0, 1]) 571 %!assert (size (max (sparse (0, 8), [], 2)), [0, 1])
572 %!assert (size (min (sparse (8, 0), [], 1)), [1, 0]) 572 %!assert (size (min (sparse (8, 0), [], 1)), [1, 0])
573 %!assert (size (min (sparse (8, 0), [], 2)), [8, 0]) 573 %!assert (size (min (sparse (8, 0), [], 2)), [8, 0])
574 %!assert (size (min (sparse (0, 8), [], 1)), [0, 8]) 574 %!assert (size (min (sparse (0, 8), [], 1)), [0, 8])
575 %!assert (size (min (sparse (0, 8), [], 2)), [0, 1]) 575 %!assert (size (min (sparse (0, 8), [], 2)), [0, 1])
576 576
577 */ 577 */
578 578
579 ComplexRowVector 579 ComplexRowVector
580 SparseComplexMatrix::row (octave_idx_type i) const 580 SparseComplexMatrix::row (octave_idx_type i) const
608 } 608 }
609 609
610 // destructive insert/delete/reorder operations 610 // destructive insert/delete/reorder operations
611 611
612 SparseComplexMatrix& 612 SparseComplexMatrix&
613 SparseComplexMatrix::insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c) 613 SparseComplexMatrix::insert (const SparseMatrix& a,
614 octave_idx_type r, octave_idx_type c)
614 { 615 {
615 SparseComplexMatrix tmp (a); 616 SparseComplexMatrix tmp (a);
616 return insert (tmp /*a*/, r, c); 617 return insert (tmp /*a*/, r, c);
617 } 618 }
618 619
619 SparseComplexMatrix& 620 SparseComplexMatrix&
620 SparseComplexMatrix::insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c) 621 SparseComplexMatrix::insert (const SparseComplexMatrix& a,
622 octave_idx_type r, octave_idx_type c)
621 { 623 {
622 MSparse<Complex>::insert (a, r, c); 624 MSparse<Complex>::insert (a, r, c);
623 return *this; 625 return *this;
624 } 626 }
625 627
626 SparseComplexMatrix& 628 SparseComplexMatrix&
627 SparseComplexMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx) 629 SparseComplexMatrix::insert (const SparseMatrix& a,
630 const Array<octave_idx_type>& indx)
628 { 631 {
629 SparseComplexMatrix tmp (a); 632 SparseComplexMatrix tmp (a);
630 return insert (tmp /*a*/, indx); 633 return insert (tmp /*a*/, indx);
631 } 634 }
632 635
633 SparseComplexMatrix& 636 SparseComplexMatrix&
634 SparseComplexMatrix::insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx) 637 SparseComplexMatrix::insert (const SparseComplexMatrix& a,
638 const Array<octave_idx_type>& indx)
635 { 639 {
636 MSparse<Complex>::insert (a, indx); 640 MSparse<Complex>::insert (a, indx);
637 return *this; 641 return *this;
638 } 642 }
639 643
646 insert (rb, ra_idx(0), ra_idx(1)); 650 insert (rb, ra_idx(0), ra_idx(1));
647 return *this; 651 return *this;
648 } 652 }
649 653
650 SparseComplexMatrix 654 SparseComplexMatrix
651 SparseComplexMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx) 655 SparseComplexMatrix::concat (const SparseMatrix& rb,
656 const Array<octave_idx_type>& ra_idx)
652 { 657 {
653 SparseComplexMatrix tmp (rb); 658 SparseComplexMatrix tmp (rb);
654 if (rb.rows () > 0 && rb.cols () > 0) 659 if (rb.rows () > 0 && rb.cols () > 0)
655 insert (tmp, ra_idx(0), ra_idx(1)); 660 insert (tmp, ra_idx(0), ra_idx(1));
656 return *this; 661 return *this;
740 return inverse (mattype, info, rcond, 0, 0); 745 return inverse (mattype, info, rcond, 0, 0);
741 } 746 }
742 747
743 SparseComplexMatrix 748 SparseComplexMatrix
744 SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info, 749 SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info,
745 double& rcond, const bool, 750 double& rcond, const bool,
746 const bool calccond) const 751 const bool calccond) const
747 { 752 {
748 SparseComplexMatrix retval; 753 SparseComplexMatrix retval;
749 754
750 octave_idx_type nr = rows (); 755 octave_idx_type nr = rows ();
751 octave_idx_type nc = cols (); 756 octave_idx_type nc = cols ();
1066 (*current_liboctave_error_handler) ("incorrect matrix type"); 1071 (*current_liboctave_error_handler) ("incorrect matrix type");
1067 } 1072 }
1068 1073
1069 return retval; 1074 return retval;
1070 1075
1071 inverse_singular: 1076 inverse_singular:
1072 return SparseComplexMatrix (); 1077 return SparseComplexMatrix ();
1073 } 1078 }
1074 1079
1075 SparseComplexMatrix 1080 SparseComplexMatrix
1076 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info, 1081 SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info,
1101 if (info == 0) 1106 if (info == 0)
1102 { 1107 {
1103 double rcond2; 1108 double rcond2;
1104 SparseMatrix Q = fact.Q (); 1109 SparseMatrix Q = fact.Q ();
1105 SparseComplexMatrix InvL = fact.L ().transpose (). 1110 SparseComplexMatrix InvL = fact.L ().transpose ().
1106 tinverse (tmp_typ, info, rcond2, true, false); 1111 tinverse (tmp_typ, info, rcond2,
1112 true, false);
1107 ret = Q * InvL.hermitian () * InvL * Q.transpose (); 1113 ret = Q * InvL.hermitian () * InvL * Q.transpose ();
1108 } 1114 }
1109 else 1115 else
1110 { 1116 {
1111 // Matrix is either singular or not positive definite 1117 // Matrix is either singular or not positive definite
1124 MatrixType tmp_typ (MatrixType::Upper); 1130 MatrixType tmp_typ (MatrixType::Upper);
1125 SparseComplexLU fact (*this, Qinit, Matrix (), false, false); 1131 SparseComplexLU fact (*this, Qinit, Matrix (), false, false);
1126 rcond = fact.rcond (); 1132 rcond = fact.rcond ();
1127 double rcond2; 1133 double rcond2;
1128 SparseComplexMatrix InvL = fact.L ().transpose (). 1134 SparseComplexMatrix InvL = fact.L ().transpose ().
1129 tinverse (tmp_typ, info, rcond2, true, false); 1135 tinverse (tmp_typ, info, rcond2,
1136 true, false);
1130 SparseComplexMatrix InvU = fact.U (). 1137 SparseComplexMatrix InvU = fact.U ().
1131 tinverse (tmp_typ, info, rcond2, true, false).transpose (); 1138 tinverse (tmp_typ, info, rcond2,
1139 true, false).transpose ();
1132 ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr (); 1140 ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr ();
1133 } 1141 }
1134 } 1142 }
1135 1143
1136 return ret; 1144 return ret;
1150 double rcond; 1158 double rcond;
1151 return determinant (info, rcond, 0); 1159 return determinant (info, rcond, 0);
1152 } 1160 }
1153 1161
1154 ComplexDET 1162 ComplexDET
1155 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int) const 1163 SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond,
1164 int) const
1156 { 1165 {
1157 ComplexDET retval; 1166 ComplexDET retval;
1158 #ifdef HAVE_UMFPACK 1167 #ifdef HAVE_UMFPACK
1159 1168
1160 octave_idx_type nr = rows (); 1169 octave_idx_type nr = rows ();
1204 1213
1205 void *Symbolic; 1214 void *Symbolic;
1206 Matrix Info (1, UMFPACK_INFO); 1215 Matrix Info (1, UMFPACK_INFO);
1207 double *info = Info.fortran_vec (); 1216 double *info = Info.fortran_vec ();
1208 int status = UMFPACK_ZNAME (qsymbolic) 1217 int status = UMFPACK_ZNAME (qsymbolic)
1209 (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, 1218 (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0,
1210 0, &Symbolic, control, info); 1219 0, &Symbolic, control, info);
1211 1220
1212 if (status < 0) 1221 if (status < 0)
1213 { 1222 {
1214 (*current_liboctave_error_handler) 1223 (*current_liboctave_error_handler)
1215 ("SparseComplexMatrix::determinant symbolic factorization failed"); 1224 ("SparseComplexMatrix::determinant symbolic factorization failed");
1216 1225
1217 UMFPACK_ZNAME (report_status) (control, status); 1226 UMFPACK_ZNAME (report_status) (control, status);
1218 UMFPACK_ZNAME (report_info) (control, info); 1227 UMFPACK_ZNAME (report_info) (control, info);
1219 1228
1220 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 1229 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1221 } 1230 }
1222 else 1231 else
1223 { 1232 {
1224 UMFPACK_ZNAME (report_symbolic) (Symbolic, control); 1233 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
1225 1234
1226 void *Numeric; 1235 void *Numeric;
1227 status 1236 status
1228 = UMFPACK_ZNAME (numeric) (Ap, Ai, 1237 = UMFPACK_ZNAME (numeric) (Ap, Ai,
1229 reinterpret_cast<const double *> (Ax), 1238 reinterpret_cast<const double *> (Ax),
1230 0, Symbolic, &Numeric, control, info) ; 1239 0, Symbolic, &Numeric, control, info);
1231 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 1240 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
1232 1241
1233 rcond = Info (UMFPACK_RCOND); 1242 rcond = Info (UMFPACK_RCOND);
1234 1243
1235 if (status < 0) 1244 if (status < 0)
1236 { 1245 {
1300 typ == MatrixType::Permuted_Diagonal) 1309 typ == MatrixType::Permuted_Diagonal)
1301 { 1310 {
1302 retval.resize (nc, b.cols (), Complex (0.,0.)); 1311 retval.resize (nc, b.cols (), Complex (0.,0.));
1303 if (typ == MatrixType::Diagonal) 1312 if (typ == MatrixType::Diagonal)
1304 for (octave_idx_type j = 0; j < b.cols (); j++) 1313 for (octave_idx_type j = 0; j < b.cols (); j++)
1305 for (octave_idx_type i = 0; i < nm; i++) 1314 for (octave_idx_type i = 0; i < nm; i++)
1306 retval(i,j) = b(i,j) / data (i); 1315 retval(i,j) = b(i,j) / data (i);
1307 else 1316 else
1308 for (octave_idx_type j = 0; j < b.cols (); j++) 1317 for (octave_idx_type j = 0; j < b.cols (); j++)
1309 for (octave_idx_type k = 0; k < nc; k++) 1318 for (octave_idx_type k = 0; k < nc; k++)
1310 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) 1319 for (octave_idx_type i = cidx (k); i < cidx (k+1); i++)
1311 retval(k,j) = b(ridx (i),j) / data (i); 1320 retval(k,j) = b(ridx (i),j) / data (i);
3897 } 3906 }
3898 else 3907 else
3899 rcond = 1.; 3908 rcond = 1.;
3900 } 3909 }
3901 else if (typ != MatrixType::Tridiagonal_Hermitian) 3910 else if (typ != MatrixType::Tridiagonal_Hermitian)
3902 (*current_liboctave_error_handler) ("incorrect matrix type"); 3911 (*current_liboctave_error_handler) ("incorrect matrix type");
3903 } 3912 }
3904 3913
3905 return retval; 3914 return retval;
3906 } 3915 }
3907 3916
4437 Complex *pz = z.fortran_vec (); 4446 Complex *pz = z.fortran_vec ();
4438 Array<double> iz (dim_vector (nr, 1)); 4447 Array<double> iz (dim_vector (nr, 1));
4439 double *piz = iz.fortran_vec (); 4448 double *piz = iz.fortran_vec ();
4440 4449
4441 F77_XFCN (zpbcon, ZPBCON, 4450 F77_XFCN (zpbcon, ZPBCON,
4442 (F77_CONST_CHAR_ARG2 (&job, 1), 4451 (F77_CONST_CHAR_ARG2 (&job, 1),
4443 nr, n_lower, tmp_data, ldm, 4452 nr, n_lower, tmp_data, ldm,
4444 anorm, rcond, pz, piz, err 4453 anorm, rcond, pz, piz, err
4445 F77_CHAR_ARG_LEN (1))); 4454 F77_CHAR_ARG_LEN (1)));
4446 4455
4447 if (err != 0) 4456 if (err != 0)
4448 err = -2; 4457 err = -2;
4449 4458
4450 volatile double rcond_plus_one = rcond + 1.0; 4459 volatile double rcond_plus_one = rcond + 1.0;
4558 Complex *pz = z.fortran_vec (); 4567 Complex *pz = z.fortran_vec ();
4559 Array<double> iz (dim_vector (nr, 1)); 4568 Array<double> iz (dim_vector (nr, 1));
4560 double *piz = iz.fortran_vec (); 4569 double *piz = iz.fortran_vec ();
4561 4570
4562 F77_XFCN (zgbcon, ZGBCON, 4571 F77_XFCN (zgbcon, ZGBCON,
4563 (F77_CONST_CHAR_ARG2 (&job, 1), 4572 (F77_CONST_CHAR_ARG2 (&job, 1),
4564 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 4573 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4565 anorm, rcond, pz, piz, err 4574 anorm, rcond, pz, piz, err
4566 F77_CHAR_ARG_LEN (1))); 4575 F77_CHAR_ARG_LEN (1)));
4567 4576
4568 if (err != 0) 4577 if (err != 0)
4569 err = -2; 4578 err = -2;
4570 4579
4571 volatile double rcond_plus_one = rcond + 1.0; 4580 volatile double rcond_plus_one = rcond + 1.0;
4572 4581
4573 if (rcond_plus_one == 1.0 || xisnan (rcond)) 4582 if (rcond_plus_one == 1.0 || xisnan (rcond))
4684 Complex *pz = z.fortran_vec (); 4693 Complex *pz = z.fortran_vec ();
4685 Array<double> iz (dim_vector (nr, 1)); 4694 Array<double> iz (dim_vector (nr, 1));
4686 double *piz = iz.fortran_vec (); 4695 double *piz = iz.fortran_vec ();
4687 4696
4688 F77_XFCN (zpbcon, ZPBCON, 4697 F77_XFCN (zpbcon, ZPBCON,
4689 (F77_CONST_CHAR_ARG2 (&job, 1), 4698 (F77_CONST_CHAR_ARG2 (&job, 1),
4690 nr, n_lower, tmp_data, ldm, 4699 nr, n_lower, tmp_data, ldm,
4691 anorm, rcond, pz, piz, err 4700 anorm, rcond, pz, piz, err
4692 F77_CHAR_ARG_LEN (1))); 4701 F77_CHAR_ARG_LEN (1)));
4693 4702
4694 if (err != 0) 4703 if (err != 0)
4695 err = -2; 4704 err = -2;
4696 4705
4697 volatile double rcond_plus_one = rcond + 1.0; 4706 volatile double rcond_plus_one = rcond + 1.0;
4753 { 4762 {
4754 if (ii == x_nz) 4763 if (ii == x_nz)
4755 { 4764 {
4756 // Resize the sparse matrix 4765 // Resize the sparse matrix
4757 octave_idx_type sz = x_nz * 4766 octave_idx_type sz = x_nz *
4758 (b_nc - j) / b_nc; 4767 (b_nc - j) / b_nc;
4759 sz = (sz > 10 ? sz : 10) + x_nz; 4768 sz = (sz > 10 ? sz : 10) + x_nz;
4760 retval.change_capacity (sz); 4769 retval.change_capacity (sz);
4761 x_nz = sz; 4770 x_nz = sz;
4762 } 4771 }
4763 retval.xdata (ii) = tmp; 4772 retval.xdata (ii) = tmp;
4820 rcond = 0.0; 4829 rcond = 0.0;
4821 err = -2; 4830 err = -2;
4822 4831
4823 if (sing_handler) 4832 if (sing_handler)
4824 { 4833 {
4825 sing_handler (rcond); 4834 sing_handler (rcond);
4826 mattype.mark_as_rectangular (); 4835 mattype.mark_as_rectangular ();
4827 } 4836 }
4828 else 4837 else
4829 (*current_liboctave_error_handler) 4838 (*current_liboctave_error_handler)
4830 ("matrix singular to machine precision"); 4839 ("matrix singular to machine precision");
4831 4840
4839 Complex *pz = z.fortran_vec (); 4848 Complex *pz = z.fortran_vec ();
4840 Array<double> iz (dim_vector (nr, 1)); 4849 Array<double> iz (dim_vector (nr, 1));
4841 double *piz = iz.fortran_vec (); 4850 double *piz = iz.fortran_vec ();
4842 4851
4843 F77_XFCN (zgbcon, ZGBCON, 4852 F77_XFCN (zgbcon, ZGBCON,
4844 (F77_CONST_CHAR_ARG2 (&job, 1), 4853 (F77_CONST_CHAR_ARG2 (&job, 1),
4845 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 4854 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
4846 anorm, rcond, pz, piz, err 4855 anorm, rcond, pz, piz, err
4847 F77_CHAR_ARG_LEN (1))); 4856 F77_CHAR_ARG_LEN (1)));
4848 4857
4849 if (err != 0) 4858 if (err != 0)
4850 err = -2; 4859 err = -2;
4851 4860
4852 volatile double rcond_plus_one = rcond + 1.0; 4861 volatile double rcond_plus_one = rcond + 1.0;
4853 4862
4854 if (rcond_plus_one == 1.0 || xisnan (rcond)) 4863 if (rcond_plus_one == 1.0 || xisnan (rcond))
5004 Complex *pz = z.fortran_vec (); 5013 Complex *pz = z.fortran_vec ();
5005 Array<double> iz (dim_vector (nr, 1)); 5014 Array<double> iz (dim_vector (nr, 1));
5006 double *piz = iz.fortran_vec (); 5015 double *piz = iz.fortran_vec ();
5007 5016
5008 F77_XFCN (zpbcon, ZPBCON, 5017 F77_XFCN (zpbcon, ZPBCON,
5009 (F77_CONST_CHAR_ARG2 (&job, 1), 5018 (F77_CONST_CHAR_ARG2 (&job, 1),
5010 nr, n_lower, tmp_data, ldm, 5019 nr, n_lower, tmp_data, ldm,
5011 anorm, rcond, pz, piz, err 5020 anorm, rcond, pz, piz, err
5012 F77_CHAR_ARG_LEN (1))); 5021 F77_CHAR_ARG_LEN (1)));
5013 5022
5014 if (err != 0) 5023 if (err != 0)
5015 err = -2; 5024 err = -2;
5016 5025
5017 volatile double rcond_plus_one = rcond + 1.0; 5026 volatile double rcond_plus_one = rcond + 1.0;
5123 Complex *pz = z.fortran_vec (); 5132 Complex *pz = z.fortran_vec ();
5124 Array<double> iz (dim_vector (nr, 1)); 5133 Array<double> iz (dim_vector (nr, 1));
5125 double *piz = iz.fortran_vec (); 5134 double *piz = iz.fortran_vec ();
5126 5135
5127 F77_XFCN (zgbcon, ZGBCON, 5136 F77_XFCN (zgbcon, ZGBCON,
5128 (F77_CONST_CHAR_ARG2 (&job, 1), 5137 (F77_CONST_CHAR_ARG2 (&job, 1),
5129 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 5138 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5130 anorm, rcond, pz, piz, err 5139 anorm, rcond, pz, piz, err
5131 F77_CHAR_ARG_LEN (1))); 5140 F77_CHAR_ARG_LEN (1)));
5132 5141
5133 if (err != 0) 5142 if (err != 0)
5134 err = -2; 5143 err = -2;
5135 5144
5136 volatile double rcond_plus_one = rcond + 1.0; 5145 volatile double rcond_plus_one = rcond + 1.0;
5137 5146
5138 if (rcond_plus_one == 1.0 || xisnan (rcond)) 5147 if (rcond_plus_one == 1.0 || xisnan (rcond))
5251 Complex *pz = z.fortran_vec (); 5260 Complex *pz = z.fortran_vec ();
5252 Array<double> iz (dim_vector (nr, 1)); 5261 Array<double> iz (dim_vector (nr, 1));
5253 double *piz = iz.fortran_vec (); 5262 double *piz = iz.fortran_vec ();
5254 5263
5255 F77_XFCN (zpbcon, ZPBCON, 5264 F77_XFCN (zpbcon, ZPBCON,
5256 (F77_CONST_CHAR_ARG2 (&job, 1), 5265 (F77_CONST_CHAR_ARG2 (&job, 1),
5257 nr, n_lower, tmp_data, ldm, 5266 nr, n_lower, tmp_data, ldm,
5258 anorm, rcond, pz, piz, err 5267 anorm, rcond, pz, piz, err
5259 F77_CHAR_ARG_LEN (1))); 5268 F77_CHAR_ARG_LEN (1)));
5260 5269
5261 if (err != 0) 5270 if (err != 0)
5262 err = -2; 5271 err = -2;
5263 5272
5264 volatile double rcond_plus_one = rcond + 1.0; 5273 volatile double rcond_plus_one = rcond + 1.0;
5411 Complex *pz = z.fortran_vec (); 5420 Complex *pz = z.fortran_vec ();
5412 Array<double> iz (dim_vector (nr, 1)); 5421 Array<double> iz (dim_vector (nr, 1));
5413 double *piz = iz.fortran_vec (); 5422 double *piz = iz.fortran_vec ();
5414 5423
5415 F77_XFCN (zgbcon, ZGBCON, 5424 F77_XFCN (zgbcon, ZGBCON,
5416 (F77_CONST_CHAR_ARG2 (&job, 1), 5425 (F77_CONST_CHAR_ARG2 (&job, 1),
5417 nc, n_lower, n_upper, tmp_data, ldm, pipvt, 5426 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
5418 anorm, rcond, pz, piz, err 5427 anorm, rcond, pz, piz, err
5419 F77_CHAR_ARG_LEN (1))); 5428 F77_CHAR_ARG_LEN (1)));
5420 5429
5421 if (err != 0) 5430 if (err != 0)
5422 err = -2; 5431 err = -2;
5423 5432
5424 volatile double rcond_plus_one = rcond + 1.0; 5433 volatile double rcond_plus_one = rcond + 1.0;
5425 5434
5426 if (rcond_plus_one == 1.0 || xisnan (rcond)) 5435 if (rcond_plus_one == 1.0 || xisnan (rcond))
5547 5556
5548 void *Symbolic; 5557 void *Symbolic;
5549 Info = Matrix (1, UMFPACK_INFO); 5558 Info = Matrix (1, UMFPACK_INFO);
5550 double *info = Info.fortran_vec (); 5559 double *info = Info.fortran_vec ();
5551 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, 5560 int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai,
5552 reinterpret_cast<const double *> (Ax), 5561 reinterpret_cast<const double *> (Ax),
5553 0, 0, &Symbolic, control, info); 5562 0, 0, &Symbolic, control, info);
5554 5563
5555 if (status < 0) 5564 if (status < 0)
5556 { 5565 {
5557 (*current_liboctave_error_handler) 5566 (*current_liboctave_error_handler)
5558 ("SparseComplexMatrix::solve symbolic factorization failed"); 5567 ("SparseComplexMatrix::solve symbolic factorization failed");
5559 err = -1; 5568 err = -1;
5560 5569
5561 UMFPACK_ZNAME (report_status) (control, status); 5570 UMFPACK_ZNAME (report_status) (control, status);
5562 UMFPACK_ZNAME (report_info) (control, info); 5571 UMFPACK_ZNAME (report_info) (control, info);
5563 5572
5564 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 5573 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5565 } 5574 }
5566 else 5575 else
5567 { 5576 {
5568 UMFPACK_ZNAME (report_symbolic) (Symbolic, control); 5577 UMFPACK_ZNAME (report_symbolic) (Symbolic, control);
5569 5578
5570 status = UMFPACK_ZNAME (numeric) (Ap, Ai, 5579 status = UMFPACK_ZNAME (numeric) (Ap, Ai,
5571 reinterpret_cast<const double *> (Ax), 0, 5580 reinterpret_cast<const double *> (Ax),
5572 Symbolic, &Numeric, control, info) ; 5581 0, Symbolic, &Numeric, control, info);
5573 UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; 5582 UMFPACK_ZNAME (free_symbolic) (&Symbolic);
5574 5583
5575 if (calc_cond) 5584 if (calc_cond)
5576 rcond = Info (UMFPACK_RCOND); 5585 rcond = Info (UMFPACK_RCOND);
5577 else 5586 else
5578 rcond = 1.; 5587 rcond = 1.;
5592 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", 5601 ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g",
5593 rcond); 5602 rcond);
5594 5603
5595 } 5604 }
5596 else if (status < 0) 5605 else if (status < 0)
5597 { 5606 {
5598 (*current_liboctave_error_handler) 5607 (*current_liboctave_error_handler)
5599 ("SparseComplexMatrix::solve numeric factorization failed"); 5608 ("SparseComplexMatrix::solve numeric factorization failed");
5600 5609
5601 UMFPACK_ZNAME (report_status) (control, status); 5610 UMFPACK_ZNAME (report_status) (control, status);
5602 UMFPACK_ZNAME (report_info) (control, info); 5611 UMFPACK_ZNAME (report_info) (control, info);
5603 5612
5604 err = -1; 5613 err = -1;
5605 } 5614 }
5606 else 5615 else
5607 { 5616 {
5608 UMFPACK_ZNAME (report_numeric) (Numeric, control); 5617 UMFPACK_ZNAME (report_numeric) (Numeric, control);
5609 } 5618 }
5610 } 5619 }
5611 5620
5612 if (err != 0) 5621 if (err != 0)
5613 UMFPACK_ZNAME (free_numeric) (&Numeric); 5622 UMFPACK_ZNAME (free_numeric) (&Numeric);
5614 #else 5623 #else
5806 5815
5807 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) 5816 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
5808 { 5817 {
5809 #ifdef UMFPACK_SEPARATE_SPLIT 5818 #ifdef UMFPACK_SEPARATE_SPLIT
5810 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, 5819 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5811 Ai, 5820 Ai,
5812 reinterpret_cast<const double *> (Ax), 5821 reinterpret_cast<const double *> (Ax),
5813 0, 5822 0,
5814 reinterpret_cast<double *> (&Xx[iidx]), 5823 reinterpret_cast<double *> (&Xx[iidx]),
5815 0, 5824 0,
5816 &Bx[iidx], Bz, Numeric, 5825 &Bx[iidx], Bz, Numeric,
5817 control, info); 5826 control, info);
5818 #else 5827 #else
5819 for (octave_idx_type i = 0; i < b_nr; i++) 5828 for (octave_idx_type i = 0; i < b_nr; i++)
5820 Bz[i] = b.elem (i, j); 5829 Bz[i] = b.elem (i, j);
5821 5830
5822 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, 5831 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
5823 Ai, 5832 Ai,
5824 reinterpret_cast<const double *> (Ax), 5833 reinterpret_cast<const double *> (Ax),
5825 0, 5834 0,
5826 reinterpret_cast<double *> (&Xx[iidx]), 5835 reinterpret_cast<double *> (&Xx[iidx]),
5827 0, 5836 0,
5828 reinterpret_cast<const double *> (Bz), 5837 reinterpret_cast<const double *> (Bz),
5829 0, Numeric, 5838 0, Numeric,
5830 control, info); 5839 control, info);
5831 #endif 5840 #endif
5832 5841
5833 if (status < 0) 5842 if (status < 0)
5834 { 5843 {
5835 (*current_liboctave_error_handler) 5844 (*current_liboctave_error_handler)
6003 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6012 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6004 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); 6013 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6005 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6014 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6006 6015
6007 retval = SparseComplexMatrix 6016 retval = SparseComplexMatrix
6008 (static_cast<octave_idx_type>(X->nrow), 6017 (static_cast<octave_idx_type>(X->nrow),
6009 static_cast<octave_idx_type>(X->ncol), 6018 static_cast<octave_idx_type>(X->ncol),
6010 static_cast<octave_idx_type>(X->nzmax)); 6019 static_cast<octave_idx_type>(X->nzmax));
6011 for (octave_idx_type j = 0; 6020 for (octave_idx_type j = 0;
6012 j <= static_cast<octave_idx_type>(X->ncol); j++) 6021 j <= static_cast<octave_idx_type>(X->ncol); j++)
6013 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j]; 6022 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6014 for (octave_idx_type j = 0; 6023 for (octave_idx_type j = 0;
6015 j < static_cast<octave_idx_type>(X->nzmax); j++) 6024 j < static_cast<octave_idx_type>(X->nzmax); j++)
6077 #ifdef UMFPACK_SEPARATE_SPLIT 6086 #ifdef UMFPACK_SEPARATE_SPLIT
6078 for (octave_idx_type i = 0; i < b_nr; i++) 6087 for (octave_idx_type i = 0; i < b_nr; i++)
6079 Bx[i] = b.elem (i, j); 6088 Bx[i] = b.elem (i, j);
6080 6089
6081 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, 6090 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
6082 Ai, 6091 Ai,
6083 reinterpret_cast<const double *> (Ax), 6092 reinterpret_cast<const double *> (Ax),
6084 0, 6093 0,
6085 reinterpret_cast<double *> (Xx), 6094 reinterpret_cast<double *> (Xx),
6086 0, 6095 0,
6087 Bx, Bz, Numeric, control, 6096 Bx, Bz, Numeric, control,
6088 info); 6097 info);
6089 #else 6098 #else
6090 for (octave_idx_type i = 0; i < b_nr; i++) 6099 for (octave_idx_type i = 0; i < b_nr; i++)
6091 Bz[i] = b.elem (i, j); 6100 Bz[i] = b.elem (i, j);
6092 6101
6093 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, 6102 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
6094 reinterpret_cast<const double *> (Ax), 6103 reinterpret_cast<const double *> (Ax),
6095 0, 6104 0,
6096 reinterpret_cast<double *> (Xx), 6105 reinterpret_cast<double *> (Xx),
6097 0, 6106 0,
6098 reinterpret_cast<double *> (Bz), 6107 reinterpret_cast<double *> (Bz),
6099 0, 6108 0,
6100 Numeric, control, 6109 Numeric, control,
6101 info); 6110 info);
6102 #endif 6111 #endif
6103 if (status < 0) 6112 if (status < 0)
6104 { 6113 {
6105 (*current_liboctave_error_handler) 6114 (*current_liboctave_error_handler)
6106 ("SparseComplexMatrix::solve solve failed"); 6115 ("SparseComplexMatrix::solve solve failed");
6334 6343
6335 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) 6344 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
6336 { 6345 {
6337 status = 6346 status =
6338 UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, 6347 UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai,
6339 reinterpret_cast<const double *> (Ax), 6348 reinterpret_cast<const double *> (Ax),
6340 0, 6349 0,
6341 reinterpret_cast<double *> (&Xx[iidx]), 6350 reinterpret_cast<double *> (&Xx[iidx]),
6342 0, 6351 0,
6343 reinterpret_cast<const double *> (&Bx[iidx]), 6352 reinterpret_cast<const double *> (&Bx[iidx]),
6344 0, Numeric, control, info); 6353 0, Numeric, control, info);
6345 6354
6346 if (status < 0) 6355 if (status < 0)
6347 { 6356 {
6348 (*current_liboctave_error_handler) 6357 (*current_liboctave_error_handler)
6349 ("SparseComplexMatrix::solve solve failed"); 6358 ("SparseComplexMatrix::solve solve failed");
6516 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6525 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6517 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); 6526 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
6518 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; 6527 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
6519 6528
6520 retval = SparseComplexMatrix 6529 retval = SparseComplexMatrix
6521 (static_cast<octave_idx_type>(X->nrow), 6530 (static_cast<octave_idx_type>(X->nrow),
6522 static_cast<octave_idx_type>(X->ncol), 6531 static_cast<octave_idx_type>(X->ncol),
6523 static_cast<octave_idx_type>(X->nzmax)); 6532 static_cast<octave_idx_type>(X->nzmax));
6524 for (octave_idx_type j = 0; 6533 for (octave_idx_type j = 0;
6525 j <= static_cast<octave_idx_type>(X->ncol); j++) 6534 j <= static_cast<octave_idx_type>(X->ncol); j++)
6526 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j]; 6535 retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j];
6527 for (octave_idx_type j = 0; 6536 for (octave_idx_type j = 0;
6528 j < static_cast<octave_idx_type>(X->nzmax); j++) 6537 j < static_cast<octave_idx_type>(X->nzmax); j++)
6581 { 6590 {
6582 for (octave_idx_type i = 0; i < b_nr; i++) 6591 for (octave_idx_type i = 0; i < b_nr; i++)
6583 Bx[i] = b (i,j); 6592 Bx[i] = b (i,j);
6584 6593
6585 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, 6594 status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap,
6586 Ai, 6595 Ai,
6587 reinterpret_cast<const double *> (Ax), 6596 reinterpret_cast<const double *> (Ax),
6588 0, 6597 0,
6589 reinterpret_cast<double *> (Xx), 6598 reinterpret_cast<double *> (Xx),
6590 0, 6599 0,
6591 reinterpret_cast<double *> (Bx), 6600 reinterpret_cast<double *> (Bx),
6592 0, Numeric, control, info); 6601 0, Numeric, control, info);
6593 6602
6594 if (status < 0) 6603 if (status < 0)
6595 { 6604 {
6596 (*current_liboctave_error_handler) 6605 (*current_liboctave_error_handler)
6597 ("SparseComplexMatrix::solve solve failed"); 6606 ("SparseComplexMatrix::solve solve failed");
6718 { 6727 {
6719 rcond = 1.; 6728 rcond = 1.;
6720 #ifdef USE_QRSOLVE 6729 #ifdef USE_QRSOLVE
6721 retval = qrsolve (*this, b, err); 6730 retval = qrsolve (*this, b, err);
6722 #else 6731 #else
6723 retval = dmsolve<ComplexMatrix, SparseComplexMatrix, 6732 retval = dmsolve<ComplexMatrix, SparseComplexMatrix, Matrix>
6724 Matrix> (*this, b, err); 6733 (*this, b, err);
6725 #endif 6734 #endif
6726 } 6735 }
6727 6736
6728 return retval; 6737 return retval;
6729 } 6738 }
6736 return solve (mattype, b, info, rcond, 0); 6745 return solve (mattype, b, info, rcond, 0);
6737 } 6746 }
6738 6747
6739 SparseComplexMatrix 6748 SparseComplexMatrix
6740 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b, 6749 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
6741 octave_idx_type& info) const 6750 octave_idx_type& info) const
6742 { 6751 {
6743 double rcond; 6752 double rcond;
6744 return solve (mattype, b, info, rcond, 0); 6753 return solve (mattype, b, info, rcond, 0);
6745 } 6754 }
6746 6755
6747 SparseComplexMatrix 6756 SparseComplexMatrix
6748 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b, 6757 SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
6749 octave_idx_type& info, double& rcond) const 6758 octave_idx_type& info, double& rcond) const
6750 { 6759 {
6751 return solve (mattype, b, info, rcond, 0); 6760 return solve (mattype, b, info, rcond, 0);
6752 } 6761 }
6753 6762
6754 SparseComplexMatrix 6763 SparseComplexMatrix
6786 { 6795 {
6787 rcond = 1.; 6796 rcond = 1.;
6788 #ifdef USE_QRSOLVE 6797 #ifdef USE_QRSOLVE
6789 retval = qrsolve (*this, b, err); 6798 retval = qrsolve (*this, b, err);
6790 #else 6799 #else
6791 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, 6800 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, SparseMatrix>
6792 SparseMatrix> (*this, b, err); 6801 (*this, b, err);
6793 #endif 6802 #endif
6794 } 6803 }
6795 6804
6796 return retval; 6805 return retval;
6797 } 6806 }
6854 { 6863 {
6855 rcond = 1.; 6864 rcond = 1.;
6856 #ifdef USE_QRSOLVE 6865 #ifdef USE_QRSOLVE
6857 retval = qrsolve (*this, b, err); 6866 retval = qrsolve (*this, b, err);
6858 #else 6867 #else
6859 retval = dmsolve<ComplexMatrix, SparseComplexMatrix, 6868 retval = dmsolve<ComplexMatrix, SparseComplexMatrix, ComplexMatrix>
6860 ComplexMatrix> (*this, b, err); 6869 (*this, b, err);
6861 #endif 6870 #endif
6862 } 6871 }
6863 6872
6864 return retval; 6873 return retval;
6865 } 6874 }
6924 rcond = 1.; 6933 rcond = 1.;
6925 #ifdef USE_QRSOLVE 6934 #ifdef USE_QRSOLVE
6926 retval = qrsolve (*this, b, err); 6935 retval = qrsolve (*this, b, err);
6927 #else 6936 #else
6928 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, 6937 retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix,
6929 SparseComplexMatrix> (*this, b, err); 6938 SparseComplexMatrix> (*this, b, err);
6930 #endif 6939 #endif
6931 } 6940 }
6932 6941
6933 return retval; 6942 return retval;
6934 } 6943 }
6959 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b, 6968 SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b,
6960 octave_idx_type& info, double& rcond, 6969 octave_idx_type& info, double& rcond,
6961 solve_singularity_handler sing_handler) const 6970 solve_singularity_handler sing_handler) const
6962 { 6971 {
6963 Matrix tmp (b); 6972 Matrix tmp (b);
6964 return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); 6973 return solve (mattype, tmp, info, rcond,
6974 sing_handler).column (static_cast<octave_idx_type> (0));
6965 } 6975 }
6966 6976
6967 ComplexColumnVector 6977 ComplexColumnVector
6968 SparseComplexMatrix::solve (MatrixType &mattype, 6978 SparseComplexMatrix::solve (MatrixType &mattype,
6969 const ComplexColumnVector& b) const 6979 const ComplexColumnVector& b) const
6992 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, 7002 SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b,
6993 octave_idx_type& info, double& rcond, 7003 octave_idx_type& info, double& rcond,
6994 solve_singularity_handler sing_handler) const 7004 solve_singularity_handler sing_handler) const
6995 { 7005 {
6996 ComplexMatrix tmp (b); 7006 ComplexMatrix tmp (b);
6997 return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); 7007 return solve (mattype, tmp, info, rcond,
7008 sing_handler).column (static_cast<octave_idx_type> (0));
6998 } 7009 }
6999 7010
7000 ComplexMatrix 7011 ComplexMatrix
7001 SparseComplexMatrix::solve (const Matrix& b) const 7012 SparseComplexMatrix::solve (const Matrix& b) const
7002 { 7013 {
7012 return solve (b, info, rcond, 0); 7023 return solve (b, info, rcond, 0);
7013 } 7024 }
7014 7025
7015 ComplexMatrix 7026 ComplexMatrix
7016 SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info, 7027 SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info,
7017 double& rcond) const 7028 double& rcond) const
7018 { 7029 {
7019 return solve (b, info, rcond, 0); 7030 return solve (b, info, rcond, 0);
7020 } 7031 }
7021 7032
7022 ComplexMatrix 7033 ComplexMatrix
7036 return solve (b, info, rcond, 0); 7047 return solve (b, info, rcond, 0);
7037 } 7048 }
7038 7049
7039 SparseComplexMatrix 7050 SparseComplexMatrix
7040 SparseComplexMatrix::solve (const SparseMatrix& b, 7051 SparseComplexMatrix::solve (const SparseMatrix& b,
7041 octave_idx_type& info) const 7052 octave_idx_type& info) const
7042 { 7053 {
7043 double rcond; 7054 double rcond;
7044 return solve (b, info, rcond, 0); 7055 return solve (b, info, rcond, 0);
7045 } 7056 }
7046 7057
7047 SparseComplexMatrix 7058 SparseComplexMatrix
7048 SparseComplexMatrix::solve (const SparseMatrix& b, 7059 SparseComplexMatrix::solve (const SparseMatrix& b,
7049 octave_idx_type& info, double& rcond) const 7060 octave_idx_type& info, double& rcond) const
7050 { 7061 {
7051 return solve (b, info, rcond, 0); 7062 return solve (b, info, rcond, 0);
7052 } 7063 }
7053 7064
7054 SparseComplexMatrix 7065 SparseComplexMatrix
7055 SparseComplexMatrix::solve (const SparseMatrix& b, 7066 SparseComplexMatrix::solve (const SparseMatrix& b,
7056 octave_idx_type& err, double& rcond, 7067 octave_idx_type& err, double& rcond,
7057 solve_singularity_handler sing_handler) const 7068 solve_singularity_handler sing_handler) const
7058 { 7069 {
7059 MatrixType mattype (*this); 7070 MatrixType mattype (*this);
7060 return solve (mattype, b, err, rcond, sing_handler); 7071 return solve (mattype, b, err, rcond, sing_handler);
7061 } 7072 }
7062 7073
7068 return solve (b, info, rcond, 0); 7079 return solve (b, info, rcond, 0);
7069 } 7080 }
7070 7081
7071 ComplexMatrix 7082 ComplexMatrix
7072 SparseComplexMatrix::solve (const ComplexMatrix& b, 7083 SparseComplexMatrix::solve (const ComplexMatrix& b,
7073 octave_idx_type& info, double& rcond) const 7084 octave_idx_type& info, double& rcond) const
7074 { 7085 {
7075 return solve (b, info, rcond, 0); 7086 return solve (b, info, rcond, 0);
7076 } 7087 }
7077 7088
7078 ComplexMatrix 7089 ComplexMatrix
7079 SparseComplexMatrix::solve (const ComplexMatrix& b, 7090 SparseComplexMatrix::solve (const ComplexMatrix& b,
7080 octave_idx_type& err, double& rcond, 7091 octave_idx_type& err, double& rcond,
7081 solve_singularity_handler sing_handler) const 7092 solve_singularity_handler sing_handler) const
7082 { 7093 {
7083 MatrixType mattype (*this); 7094 MatrixType mattype (*this);
7084 return solve (mattype, b, err, rcond, sing_handler); 7095 return solve (mattype, b, err, rcond, sing_handler);
7085 } 7096 }
7086 7097
7092 return solve (b, info, rcond, 0); 7103 return solve (b, info, rcond, 0);
7093 } 7104 }
7094 7105
7095 SparseComplexMatrix 7106 SparseComplexMatrix
7096 SparseComplexMatrix::solve (const SparseComplexMatrix& b, 7107 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7097 octave_idx_type& info) const 7108 octave_idx_type& info) const
7098 { 7109 {
7099 double rcond; 7110 double rcond;
7100 return solve (b, info, rcond, 0); 7111 return solve (b, info, rcond, 0);
7101 } 7112 }
7102 7113
7103 SparseComplexMatrix 7114 SparseComplexMatrix
7104 SparseComplexMatrix::solve (const SparseComplexMatrix& b, 7115 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7105 octave_idx_type& info, double& rcond) const 7116 octave_idx_type& info, double& rcond) const
7106 { 7117 {
7107 return solve (b, info, rcond, 0); 7118 return solve (b, info, rcond, 0);
7108 } 7119 }
7109 7120
7110 SparseComplexMatrix 7121 SparseComplexMatrix
7111 SparseComplexMatrix::solve (const SparseComplexMatrix& b, 7122 SparseComplexMatrix::solve (const SparseComplexMatrix& b,
7112 octave_idx_type& err, double& rcond, 7123 octave_idx_type& err, double& rcond,
7113 solve_singularity_handler sing_handler) const 7124 solve_singularity_handler sing_handler) const
7114 { 7125 {
7115 MatrixType mattype (*this); 7126 MatrixType mattype (*this);
7116 return solve (mattype, b, err, rcond, sing_handler); 7127 return solve (mattype, b, err, rcond, sing_handler);
7117 } 7128 }
7118 7129
7136 { 7147 {
7137 return solve (b, info, rcond, 0); 7148 return solve (b, info, rcond, 0);
7138 } 7149 }
7139 7150
7140 ComplexColumnVector 7151 ComplexColumnVector
7141 SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond, 7152 SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
7153 double& rcond,
7142 solve_singularity_handler sing_handler) const 7154 solve_singularity_handler sing_handler) const
7143 { 7155 {
7144 Matrix tmp (b); 7156 Matrix tmp (b);
7145 return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); 7157 return solve (tmp, info, rcond,
7158 sing_handler).column (static_cast<octave_idx_type> (0));
7146 } 7159 }
7147 7160
7148 ComplexColumnVector 7161 ComplexColumnVector
7149 SparseComplexMatrix::solve (const ComplexColumnVector& b) const 7162 SparseComplexMatrix::solve (const ComplexColumnVector& b) const
7150 { 7163 {
7152 double rcond; 7165 double rcond;
7153 return solve (b, info, rcond, 0); 7166 return solve (b, info, rcond, 0);
7154 } 7167 }
7155 7168
7156 ComplexColumnVector 7169 ComplexColumnVector
7157 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const 7170 SparseComplexMatrix::solve (const ComplexColumnVector& b,
7171 octave_idx_type& info) const
7158 { 7172 {
7159 double rcond; 7173 double rcond;
7160 return solve (b, info, rcond, 0); 7174 return solve (b, info, rcond, 0);
7161 } 7175 }
7162 7176
7163 ComplexColumnVector 7177 ComplexColumnVector
7164 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, 7178 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
7165 double& rcond) const 7179 double& rcond) const
7166 { 7180 {
7167 return solve (b, info, rcond, 0); 7181 return solve (b, info, rcond, 0);
7168 } 7182 }
7169 7183
7170 ComplexColumnVector 7184 ComplexColumnVector
7171 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, 7185 SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
7172 double& rcond, 7186 double& rcond,
7173 solve_singularity_handler sing_handler) const 7187 solve_singularity_handler sing_handler) const
7174 { 7188 {
7175 ComplexMatrix tmp (b); 7189 ComplexMatrix tmp (b);
7176 return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); 7190 return solve (tmp, info, rcond,
7191 sing_handler).column (static_cast<octave_idx_type> (0));
7177 } 7192 }
7178 7193
7179 // unary operations 7194 // unary operations
7180 SparseBoolMatrix 7195 SparseBoolMatrix
7181 SparseComplexMatrix::operator ! (void) const 7196 SparseComplexMatrix::operator ! (void) const
7290 max_val = std::real (data (0)); 7305 max_val = std::real (data (0));
7291 min_val = std::real (data (0)); 7306 min_val = std::real (data (0));
7292 7307
7293 for (octave_idx_type i = 0; i < nel; i++) 7308 for (octave_idx_type i = 0; i < nel; i++)
7294 { 7309 {
7295 Complex val = data (i); 7310 Complex val = data (i);
7296 7311
7297 double r_val = std::real (val); 7312 double r_val = std::real (val);
7298 double i_val = std::imag (val); 7313 double i_val = std::imag (val);
7299 7314
7300 if (r_val > max_val) 7315 if (r_val > max_val)
7301 max_val = r_val; 7316 max_val = r_val;
7302 7317
7303 if (i_val > max_val) 7318 if (i_val > max_val)
7304 max_val = i_val; 7319 max_val = i_val;
7305 7320
7306 if (r_val < min_val) 7321 if (r_val < min_val)
7307 min_val = r_val; 7322 min_val = r_val;
7308 7323
7309 if (i_val < min_val) 7324 if (i_val < min_val)
7310 min_val = i_val; 7325 min_val = i_val;
7311 7326
7312 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val) 7327 if (D_NINT (r_val) != r_val || D_NINT (i_val) != i_val)
7313 return false; 7328 return false;
7314 } 7329 }
7315 7330
7316 return true; 7331 return true;
7317 } 7332 }
7318 7333
7320 SparseComplexMatrix::too_large_for_float (void) const 7335 SparseComplexMatrix::too_large_for_float (void) const
7321 { 7336 {
7322 return test_any (xtoo_large_for_float); 7337 return test_any (xtoo_large_for_float);
7323 } 7338 }
7324 7339
7325 // FIXME Do these really belong here? Maybe they should be 7340 // FIXME: Do these really belong here? Maybe they should be in a base class?
7326 // in a base class?
7327 7341
7328 SparseBoolMatrix 7342 SparseBoolMatrix
7329 SparseComplexMatrix::all (int dim) const 7343 SparseComplexMatrix::all (int dim) const
7330 { 7344 {
7331 SPARSE_ALL_OP (dim); 7345 SPARSE_ALL_OP (dim);
7413 std::ostream& 7427 std::ostream&
7414 operator << (std::ostream& os, const SparseComplexMatrix& a) 7428 operator << (std::ostream& os, const SparseComplexMatrix& a)
7415 { 7429 {
7416 octave_idx_type nc = a.cols (); 7430 octave_idx_type nc = a.cols ();
7417 7431
7418 // add one to the printed indices to go from 7432 // add one to the printed indices to go from
7419 // zero-based to one-based arrays 7433 // zero-based to one-based arrays
7420 for (octave_idx_type j = 0; j < nc; j++) 7434 for (octave_idx_type j = 0; j < nc; j++)
7421 { 7435 {
7422 octave_quit (); 7436 octave_quit ();
7423 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) 7437 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
7424 { 7438 {
7425 os << a.ridx (i) + 1 << " " << j + 1 << " "; 7439 os << a.ridx (i) + 1 << " " << j + 1 << " ";
7426 octave_write_complex (os, a.data (i)); 7440 octave_write_complex (os, a.data (i));
7427 os << "\n"; 7441 os << "\n";
7428 } 7442 }
7429 } 7443 }
7430 7444
7431 return os; 7445 return os;
7432 } 7446 }
7433 7447
7434 std::istream& 7448 std::istream&
7625 operator * (const SparseComplexMatrix& a, const PermMatrix& p) 7639 operator * (const SparseComplexMatrix& a, const PermMatrix& p)
7626 { 7640 {
7627 return octinternal_do_mul_sm_pm (a, p); 7641 return octinternal_do_mul_sm_pm (a, p);
7628 } 7642 }
7629 7643
7630 // FIXME -- it would be nice to share code among the min/max 7644 // FIXME: it would be nice to share code among the min/max functions below.
7631 // functions below.
7632 7645
7633 #define EMPTY_RETURN_CHECK(T) \ 7646 #define EMPTY_RETURN_CHECK(T) \
7634 if (nr == 0 || nc == 0) \ 7647 if (nr == 0 || nc == 0) \
7635 return T (nr, nc); 7648 return T (nr, nc);
7636 7649
7877 7890
7878 return r; 7891 return r;
7879 } 7892 }
7880 7893
7881 SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex, 7894 SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex,
7882 0.0, real) 7895 0.0, real)
7883 SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0) 7896 SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0)
7884 7897
7885 SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix, 7898 SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix,
7886 0.0, real) 7899 0.0, real)
7887 SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0) 7900 SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0)
7888 7901
7889 SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix, 7902 SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix,
7890 0.0, real) 7903 0.0, real)
7891 SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0) 7904 SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0)