octave-nkf: liboctave/dSparse.cc comparison

comparison liboctave/dSparse.cc @ 11586:12df7854fa7c

strip trailing whitespace from source files

author	John W. Eaton <jwe@octave.org>
date	Thu, 20 Jan 2011 17:24:59 -0500
parents	a83bad07f7e3
children	b646413c3d0e

comparison

equal deleted inserted replaced

-:1473d0cf86d2
+:12df7854fa7c
 SparseMatrix::SparseMatrix (const PermMatrix& a)
 : MSparse<double> (a.rows (), a.cols (), a.rows ())
 {
 octave_idx_type n = a.rows ();
 for (octave_idx_type i = 0; i <= n; i++)
 cidx (i) = i;
 const Array<octave_idx_type> pv = a.pvec ();
 if (a.is_row_perm ())
 {
 SparseMatrix result;
 dim_vector dv = dims ();
 if (dv.numel () == 0 || dim >= dv.length ())
 return result;
 if (dim < 0)
 dim = dv.first_non_singleton ();
 octave_idx_type nr = dv(0);
 octave_idx_type nc = dv(1);
 if (ridx(k) == i)
 {
 found = true;
 break;
 }
 if (!found)
 idx_arg.elem(i) = j;
 }
 SparseMatrix result;
 dim_vector dv = dims ();
 if (dv.numel () == 0 || dim >= dv.length ())
 return result;
 if (dim < 0)
 dim = dv.first_non_singleton ();
 octave_idx_type nr = dv(0);
 octave_idx_type nc = dv(1);
 if (ridx(k) == i)
 {
 found = true;
 break;
 }
 if (!found)
 idx_arg.elem(i) = j;
 }
 }
 return result;
 }
 RowVector
 SparseMatrix::row (octave_idx_type i) const
 {
 octave_idx_type nc = columns ();
 RowVector retval (nc, 0);
 }
 return retval;
 }
 ColumnVector
 SparseMatrix::column (octave_idx_type i) const
 {
 octave_idx_type nr = rows ();
 ColumnVector retval (nr);
 }
 return r;
 }
 SparseMatrix
 atan2 (const double& x, const SparseMatrix& y)
 {
 octave_idx_type nr = y.rows ();
 octave_idx_type nc = y.cols ();
 return SparseMatrix (tmp);
 }
 }
 SparseMatrix
 atan2 (const SparseMatrix& x, const double& y)
 {
 octave_idx_type nr = x.rows ();
 octave_idx_type nc = x.cols ();
 octave_idx_type nz = x.nnz ();
 }
 else
 return retval;
 }
 SparseMatrix
 atan2 (const SparseMatrix& x, const SparseMatrix& y)
 {
 SparseMatrix r;
 if ((x.rows() == y.rows()) && (x.cols() == y.cols()))
 {
 octave_idx_type x_nr = x.rows ();
 octave_idx_type x_nc = x.cols ();
 octave_idx_type y_nr = y.rows ();
 if (x_nr != y_nr || x_nc != y_nc)
 gripe_nonconformant ("atan2", x_nr, x_nc, y_nr, y_nc);
 else
 {
 r = SparseMatrix (x_nr, x_nc, (x.nnz () + y.nnz ()));
 octave_idx_type jx = 0;
 r.cidx (0) = 0;
 for (octave_idx_type i = 0 ; i < x_nc ; i++)
 {
 octave_idx_type  ja = x.cidx(i);
 octave_idx_type  ja_max = x.cidx(i+1);
 bool ja_lt_max= ja < ja_max;
 octave_idx_type  jb = y.cidx(i);
 octave_idx_type  jb_max = y.cidx(i+1);
 bool jb_lt_max = jb < jb_max;
 while (ja_lt_max || jb_lt_max )
 {
 octave_quit ();
 if ((! jb_lt_max) ||
 (ja_lt_max && (x.ridx(ja) < y.ridx(jb))))
 jb_lt_max= jb < jb_max;
 }
 }
 r.cidx(i+1) = jx;
 }
 r.maybe_compress ();
 }
 }
 else
 (*current_liboctave_error_handler) ("matrix size mismatch");
 {
 double rcond;
 return inverse (mattype, info, rcond, 0, 0);
 }
 SparseMatrix
 SparseMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info,
 double& rcond, const bool,
 const bool calccond) const
 {
 SparseMatrix retval;
 octave_idx_type nr = rows ();
 {
 if (typ == MatrixType::Permuted_Diagonal)
 retval = transpose();
 else
 retval = *this;
 // Force make_unique to be called
 double *v = retval.data();
 if (calccond)
 {
 double dmax = 0., dmin = octave_Inf;
 for (octave_idx_type i = 0; i < nr; i++)
 {
 double tmp = fabs(v[i]);
 if (tmp > dmax)
 dmax = tmp;
 }
 return retval;
 }
 SparseMatrix
 SparseMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info,
 double& rcond, const bool,
 const bool calccond) const
 {
 SparseMatrix retval;
 octave_idx_type nr = rows ();
 {
 // Print spparms("spumoni") info if requested
 int typ = mattyp.type ();
 mattyp.info ();
 if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper ||
 typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 for (octave_idx_type i = 0; i < nr; i++)
 {
 octave_quit ();
 // place the 1 in the identity position
 octave_idx_type cx_colstart = cx;
 if (cx == nz2)
 {
 nz2 *= 2;
 retval.change_capacity (nz2);
 }
 retval.xridx(cx) = i;
 retval.xdata(cx) = 1.0;
 cx++;
 // iterate accross columns of input matrix
 for (octave_idx_type j = i+1; j < nr; j++)
 {
 double v = 0.;
 // iterate to calculate sum
 octave_idx_type colXp = retval.xcidx(i);
 octave_idx_type colUp = cidx(j);
 octave_idx_type rpX, rpU;
 if (cidx(j) == cidx(j+1))
 {
 (*current_liboctave_error_handler)
 ("division by zero");
 goto inverse_singular;
 }
 do
 {
 octave_quit ();
 rpX = retval.xridx(colXp);
 rpU = ridx(colUp);
 if (rpX < rpU)
 colXp++;
 else if (rpX > rpU)
 colUp++;
 else
 {
 v -= retval.xdata(colXp) * data(colUp);
 colXp++;
 colUp++;
 }
 } while ((rpX<j) && (rpU<j) &&
 (colXp<cx) && (colUp<nz));
 // get A(m,m)
 if (typ == MatrixType::Upper)
 colUp = cidx(j+1) - 1;
 else
 colUp = cidx(j);
 double pivot = data(colUp);
 if (pivot == 0. || ridx(colUp) != j)
 {
 (*current_liboctave_error_handler)
 ("division by zero");
 goto inverse_singular;
 }
 if (v != 0.)
 if (typ == MatrixType::Upper)
 colUp = cidx(i+1) - 1;
 else
 colUp = cidx(i);
 double pivot = data(colUp);
 if (pivot == 0. || ridx(colUp) != i)
 {
 (*current_liboctave_error_handler) ("division by zero");
 goto inverse_singular;
 }
 // place the 1 in the identity position
 work[iidx] = 1.0;
 // iterate accross columns of input matrix
 for (octave_idx_type j = iidx+1; j < nr; j++)
 {
 double v = 0.;
 octave_idx_type jidx = perm[j];
 // iterate to calculate sum
 for (octave_idx_type k = cidx(jidx);
 k < cidx(jidx+1); k++)
 {
 octave_quit ();
 v -= work[ridx(k)] * data(k);
 }
 double pivot;
 if (typ == MatrixType::Permuted_Upper)
 pivot = data(cidx(jidx+1) - 1);
 else
 pivot = data(cidx(jidx));
 if (pivot == 0.)
 {
 (*current_liboctave_error_handler)
 ("division by zero");
 goto inverse_singular;
 }
 work[j] = v / pivot;
 colUp = cidx(perm[iidx]);
 double pivot = data(colUp);
 if (pivot == 0.)
 {
 (*current_liboctave_error_handler)
 ("division by zero");
 goto inverse_singular;
 }
 octave_idx_type new_cx = cx;
 {
 // Calculate the 1-norm of inverse matrix for rcond calculation
 for (octave_idx_type j = 0; j < nr; j++)
 {
 double atmp = 0.;
 for (octave_idx_type i = retval.cidx(j);
 i < retval.cidx(j+1); i++)
 atmp += fabs(retval.data(i));
 if (atmp > ainvnorm)
 ainvnorm = atmp;
 }
 rcond = 1. / ainvnorm / anorm;
 }
 }
 else
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 inverse_singular:
 return SparseMatrix();
 }
 SparseMatrix
 SparseMatrix::inverse (MatrixType &mattype, octave_idx_type& info,
 double& rcond, int, int calc_cond) const
 {
 int typ = mattype.type (false);
 SparseMatrix ret;
 MatrixType tmp_typ (MatrixType::Upper);
 SparseLU fact (*this, Qinit, Matrix(), false, false);
 rcond = fact.rcond();
 double rcond2;
 SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ,
 info, rcond2, true, false);
 SparseMatrix InvU = fact.U().tinverse(tmp_typ, info, rcond2,
 true, false).transpose();
 ret = fact.Pc().transpose() * InvU * InvL * fact.Pr();
 }
 // Set whether we are allowed to modify Q or not
 tmp = octave_sparse_params::get_key ("autoamd");
 if (!xisnan (tmp))
 Control (UMFPACK_FIXQ) = tmp;
 // Turn-off UMFPACK scaling for LU
 Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE;
 UMFPACK_DNAME (report_control) (control);
 const octave_idx_type *Ap = cidx ();
 UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control);
 void *Symbolic;
 Matrix Info (1, UMFPACK_INFO);
 double *info = Info.fortran_vec ();
 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai,
 Ax, 0, &Symbolic, control, info);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::determinant symbolic factorization failed");
 UMFPACK_DNAME (report_status) (control, status);
 UMFPACK_DNAME (report_info) (control, info);
 else
 {
 UMFPACK_DNAME (report_symbolic) (Symbolic, control);
 void *Numeric;
 status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic,
 &Numeric, control, info) ;
 UMFPACK_DNAME (free_symbolic) (&Symbolic) ;
 rcond = Info (UMFPACK_RCOND);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::determinant numeric factorization failed");
 UMFPACK_DNAME (report_status) (control, status);
 UMFPACK_DNAME (report_info) (control, info);
 status = UMFPACK_DNAME (get_determinant) (&c10, &e10, Numeric, info);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::determinant error calculating determinant");
 UMFPACK_DNAME (report_status) (control, status);
 UMFPACK_DNAME (report_info) (control, info);
 }
 else
 retval = DET (c10, e10, 10);
 return retval;
 }
 Matrix
 SparseMatrix::dsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err,
 double& rcond, solve_singularity_handler,
 bool calc_cond) const
 {
 Matrix retval;
 octave_idx_type nr = rows ();
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 retval(k,j) = b(ridx(i),j) / data (i);
 if (calc_cond)
 {
 double dmax = 0., dmin = octave_Inf;
 for (octave_idx_type i = 0; i < nm; i++)
 {
 double tmp = fabs(data(i));
 if (tmp > dmax)
 dmax = tmp;
 return retval;
 }
 SparseMatrix
 SparseMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler, bool calc_cond) const
 {
 SparseMatrix retval;
 octave_idx_type nr = rows ();
 retval.xcidx(j+1) = ii;
 }
 if (calc_cond)
 {
 double dmax = 0., dmin = octave_Inf;
 for (octave_idx_type i = 0; i < nm; i++)
 {
 double tmp = fabs(data(i));
 if (tmp > dmax)
 dmax = tmp;
 else
 for (octave_idx_type j = 0; j < b.cols(); j++)
 for (octave_idx_type k = 0; k < nc; k++)
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 retval(k,j) = b(ridx(i),j) / data (i);
 if (calc_cond)
 {
 double dmax = 0., dmin = octave_Inf;
 for (octave_idx_type i = 0; i < nm; i++)
 {
 double tmp = fabs(data(i));
 if (tmp > dmax)
 dmax = tmp;
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler, bool calc_cond) const
 {
 SparseComplexMatrix retval;
 octave_idx_type nr = rows ();
 retval.xdata (ii++) = b.data(k) / data (i);
 }
 }
 retval.xcidx(j+1) = ii;
 }
 if (calc_cond)
 {
 double dmax = 0., dmin = octave_Inf;
 for (octave_idx_type i = 0; i < nm; i++)
 {
 double tmp = fabs(data(i));
 if (tmp > dmax)
 dmax = tmp;
 }
 Matrix
 SparseMatrix::utsolve (MatrixType &mattype, const Matrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 Matrix retval;
 octave_idx_type nr = rows ();
 if (ridx(cidx(kidx+1)-1) != k ||
 data(cidx(kidx+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(kidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(kidx);
 i < cidx(kidx+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(iidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(iidx);
 i < cidx(iidx+1)-1; i++)
 {
 octave_idx_type idx2 = ridx(i);
 work[idx2] = work[idx2] - tmp * data(i);
 }
 if (ridx(cidx(k+1)-1) != k ||
 data(cidx(k+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(k+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 {
 return retval;
 }
 SparseMatrix
 SparseMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseMatrix retval;
 if (ridx(cidx(kidx+1)-1) != k ||
 data(cidx(kidx+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(kidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(kidx);
 i < cidx(kidx+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(iidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(iidx);
 i < cidx(iidx+1)-1; i++)
 {
 octave_idx_type idx2 = ridx(i);
 work[idx2] = work[idx2] - tmp * data(i);
 }
 if (ridx(cidx(k+1)-1) != k ||
 data(cidx(k+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(k+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 {
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 }
 return retval;
 }
 ComplexMatrix
 SparseMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 ComplexMatrix retval;
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Upper ||
 typ == MatrixType::Upper)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (ridx(cidx(kidx+1)-1) != k ||
 data(cidx(kidx+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(kidx+1)-1);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(kidx);
 i < cidx(kidx+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 cwork[iidx] = cwork[iidx] - tmp * data(i);
 }
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(iidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(iidx);
 i < cidx(iidx+1)-1; i++)
 {
 octave_idx_type idx2 = ridx(i);
 work[idx2] = work[idx2] - tmp * data(i);
 }
 if (ridx(cidx(k+1)-1) != k ||
 data(cidx(k+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(k+1)-1);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 {
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseComplexMatrix retval;
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Upper ||
 typ == MatrixType::Upper)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (ridx(cidx(kidx+1)-1) != k ||
 data(cidx(kidx+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(kidx+1)-1);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(kidx);
 i < cidx(kidx+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 cwork[iidx] = cwork[iidx] - tmp * data(i);
 }
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(iidx+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(iidx);
 i < cidx(iidx+1)-1; i++)
 {
 octave_idx_type idx2 = ridx(i);
 work[idx2] = work[idx2] - tmp * data(i);
 }
 if (ridx(cidx(k+1)-1) != k ||
 data(cidx(k+1)-1) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(k+1)-1);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++)
 {
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k+1)-1);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1)-1; i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Lower ||
 typ == MatrixType::Lower)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (minr != k || data(mini) == 0)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 octave_idx_type minr = nr;
 octave_idx_type mini = 0;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 if (perm[ridx(i)] < minr)
 {
 minr = perm[ridx(i)];
 mini = i;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 {
 if (i == mini)
 continue;
 if (ridx(cidx(k)) != k ||
 data(cidx(k)) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1;
 i < cidx(k+1); i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1;
 i < cidx(k+1); i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 return retval;
 }
 SparseMatrix
 SparseMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseMatrix retval;
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Lower ||
 typ == MatrixType::Lower)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (minr != k || data(mini) == 0)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 octave_idx_type minr = nr;
 octave_idx_type mini = 0;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 if (perm[ridx(i)] < minr)
 {
 minr = perm[ridx(i)];
 mini = i;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 {
 if (i == mini)
 continue;
 if (ridx(cidx(k)) != k ||
 data(cidx(k)) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1;
 i < cidx(k+1); i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 return retval;
 }
 ComplexMatrix
 SparseMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 ComplexMatrix retval;
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Lower ||
 typ == MatrixType::Lower)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (minr != k || data(mini) == 0)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(mini);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 octave_idx_type minr = nr;
 octave_idx_type mini = 0;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 if (perm[ridx(i)] < minr)
 {
 minr = perm[ridx(i)];
 mini = i;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 {
 if (i == mini)
 continue;
 if (ridx(cidx(k)) != k ||
 data(cidx(k)) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(k));
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1;
 i < cidx(k+1); i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::ltsolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseComplexMatrix retval;
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Permuted_Lower ||
 typ == MatrixType::Lower)
 {
 double anorm = 0.;
 double ainvnorm = 0.;
 if (minr != k || data(mini) == 0)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(mini);
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k); i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 octave_idx_type minr = nr;
 octave_idx_type mini = 0;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 if (perm[ridx(i)] < minr)
 {
 minr = perm[ridx(i)];
 mini = i;
 }
 double tmp = work[k] / data(mini);
 work[k] = tmp;
 for (octave_idx_type i = cidx(k);
 i < cidx(k+1); i++)
 {
 if (i == mini)
 continue;
 if (ridx(cidx(k)) != k ||
 data(cidx(k)) == 0.)
 {
 err = -2;
 goto triangular_error;
 }
 Complex tmp = cwork[k] / data(cidx(k));
 cwork[k] = tmp;
 for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++)
 {
 if (work[k] != 0.)
 {
 double tmp = work[k] / data(cidx(k));
 work[k] = tmp;
 for (octave_idx_type i = cidx(k)+1;
 i < cidx(k+1); i++)
 {
 octave_idx_type iidx = ridx(i);
 work[iidx] = work[iidx] - tmp * data(i);
 }
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
 else if (nr == 0 || b.cols () == 0)
 retval = Matrix (nc, b.cols (), 0.0);
 else if (calc_cond)
 (*current_liboctave_error_handler)
 ("calculation of condition number not implemented");
 else
 {
 // Print spparms("spumoni") info if requested
 volatile int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Tridiagonal_Hermitian)
 {
 OCTAVE_LOCAL_BUFFER (double, D, nr);
 OCTAVE_LOCAL_BUFFER (double, DL, nr - 1);
 D[j] = data(i);
 else if (ridx(i) == j + 1)
 DL[j] = data(i);
 }
 }
 octave_idx_type b_nc = b.cols();
 retval = b;
 double *result = retval.fortran_vec ();
 F77_XFCN (dptsv, DPTSV, (nr, b_nc, D, DL, result,
 b.rows(), err));
 if (err != 0)
 {
 err = 0;
 mattype.mark_as_unsymmetric ();
 typ = MatrixType::Tridiagonal;
 }
 else
 rcond = 1.;
 }
 if (typ == MatrixType::Tridiagonal)
 {
 octave_idx_type b_nc = b.cols();
 retval = b;
 double *result = retval.fortran_vec ();
 F77_XFCN (dgtsv, DGTSV, (nr, b_nc, DL, D, DU, result,
 b.rows(), err));
 if (err != 0)
 {
 rcond = 0.;
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 rcond = 1.;
 }
 else if (typ != MatrixType::Tridiagonal_Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 SparseMatrix
 SparseMatrix::trisolve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseMatrix retval;
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
 else if (nr == 0 || b.cols () == 0)
 retval = SparseMatrix (nc, b.cols ());
 else if (calc_cond)
 (*current_liboctave_error_handler)
 ("calculation of condition number not implemented");
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 // Note can't treat symmetric case as there is no dpttrf function
 if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 {
 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
 }
 }
 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
 if (err != 0)
 {
 rcond = 0.0;
 err = -2;
 if (sing_handler)
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 rcond = 1.0;
 char job = 'N';
 volatile octave_idx_type x_nz = b.nnz ();
 octave_idx_type b_nc = b.cols ();
 for (octave_idx_type i = 0; i < nr; i++)
 work[i] = 0.;
 for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
 work[b.ridx(i)] = b.data(i);
 F77_XFCN (dgttrs, DGTTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, 1, DL, D, DU, DU2, pipvt,
 work, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 // Count non-zeros in work vector and adjust
 // space in retval if needed
 octave_idx_type new_nnz = 0;
 for (octave_idx_type i = 0; i < nr; i++)
 if (work[i] != 0.)
 new_nnz++;
 return retval;
 }
 ComplexMatrix
 SparseMatrix::trisolve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 ComplexMatrix retval;
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
 else if (nr == 0 || b.cols () == 0)
 retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0));
 else if (calc_cond)
 (*current_liboctave_error_handler)
 ("calculation of condition number not implemented");
 else
 {
 // Print spparms("spumoni") info if requested
 volatile int typ = mattype.type ();
 mattype.info ();
 if (typ == MatrixType::Tridiagonal_Hermitian)
 {
 OCTAVE_LOCAL_BUFFER (double, D, nr);
 OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1);
 octave_idx_type b_nc = b.cols();
 rcond = 1.;
 retval = b;
 Complex *result = retval.fortran_vec ();
 F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, DL, result,
 b_nr, err));
 if (err != 0)
 {
 err = 0;
 octave_idx_type b_nc = b.cols();
 rcond = 1.;
 retval = b;
 Complex *result = retval.fortran_vec ();
 F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, DL, D, DU, result,
 b_nr, err));
 if (err != 0)
 {
 rcond = 0.;
 err = -2;
 if (sing_handler)
 {
 sing_handler (rcond);
 mattype.mark_as_rectangular ();
 }
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::trisolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseComplexMatrix retval;
 (*current_liboctave_error_handler)
 ("matrix dimension mismatch solution of linear equations");
 else if (nr == 0 || b.cols () == 0)
 retval = SparseComplexMatrix (nc, b.cols ());
 else if (calc_cond)
 (*current_liboctave_error_handler)
 ("calculation of condition number not implemented");
 else
 {
 // Print spparms("spumoni") info if requested
 int typ = mattype.type ();
 mattype.info ();
 // Note can't treat symmetric case as there is no dpttrf function
 if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 {
 OCTAVE_LOCAL_BUFFER (double, DU2, nr - 2);
 }
 }
 F77_XFCN (dgttrf, DGTTRF, (nr, DL, D, DU, DU2, pipvt, err));
 if (err != 0)
 {
 rcond = 0.0;
 err = -2;
 if (sing_handler)
 mattype.mark_as_rectangular ();
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 rcond = 1.;
 char job = 'N';
 octave_idx_type b_nr = b.rows ();
 octave_idx_type b_nc = b.cols ();
 OCTAVE_LOCAL_BUFFER (double, Bx, b_nr);
 Complex c = b (i,j);
 Bx[i] = std::real (c);
 Bz[i] = std::imag (c);
 }
 F77_XFCN (dgttrs, DGTTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, 1, DL, D, DU, DU2, pipvt,
 Bx, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 err = -1;
 break;
 }
 F77_XFCN (dgttrs, DGTTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, 1, DL, D, DU, DU2, pipvt,
 Bz, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 err = -1;
 break;
 }
 // Count non-zeros in work vector and adjust
 // space in retval if needed
 octave_idx_type new_nnz = 0;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 new_nnz++;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 {
 retval.xridx(ii) = i;
 retval.xdata(ii++) =
 Complex (Bx[i], Bz[i]);
 }
 retval.xcidx(j+1) = ii;
 }
 {
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 char job = 'L';
 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 // Matrix is not positive definite!! Fall through to
 // unsymmetric banded solver.
 mattype.mark_as_unsymmetric ();
 typ = MatrixType::Banded;
 rcond = 0.0;
 err = 0;
 }
 else
 {
 if (calc_cond)
 {
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dpbcon, DPBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 retval = b;
 double *result = retval.fortran_vec ();
 octave_idx_type b_nc = b.cols ();
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, b_nc, tmp_data,
 ldm, result, b.rows(), err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 err = -1;
 }
 }
 }
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 }
 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
 octave_idx_type *pipvt = ipvt.fortran_vec ();
 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 ldm, pipvt, err));
 // Throw-away extra info LAPACK gives so as to not
 // change output.
 if (err != 0)
 {
 err = -2;
 rcond = 0.0;
 if (sing_handler)
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 if (calc_cond)
 {
 char job = '1';
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dgbcon, DGBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 double *result = retval.fortran_vec ();
 octave_idx_type b_nc = b.cols ();
 char job = 'N';
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, b_nc, tmp_data,
 ldm, pipvt, result, b.rows(), err
 F77_CHAR_ARG_LEN (1)));
 }
 return retval;
 }
 SparseMatrix
 SparseMatrix::bsolve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseMatrix retval;
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 char job = 'L';
 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 mattype.mark_as_unsymmetric ();
 typ = MatrixType::Banded;
 rcond = 0.0;
 err = 0;
 }
 else
 {
 if (calc_cond)
 {
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dpbcon, DPBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 for (volatile octave_idx_type j = 0; j < b_nc; j++)
 {
 for (octave_idx_type i = 0; i < b_nr; i++)
 Bx[i] = b.elem (i, j);
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, 1, tmp_data,
 ldm, Bx, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 err = -1;
 break;
 }
 if (tmp != 0.0)
 {
 if (ii == x_nz)
 {
 // Resize the sparse matrix
 octave_idx_type sz = x_nz *
 (b_nc - j) / b_nc;
 sz = (sz > 10 ? sz : 10) + x_nz;
 retval.change_capacity (sz);
 x_nz = sz;
 }
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 }
 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
 octave_idx_type *pipvt = ipvt.fortran_vec ();
 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 ldm, pipvt, err));
 if (err != 0)
 {
 err = -2;
 rcond = 0.0;
 if (sing_handler)
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 if (calc_cond)
 {
 char job = '1';
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dgbcon, DGBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 for (volatile octave_idx_type j = 0; j < b_nc; j++)
 {
 for (octave_idx_type i = 0; i < nr; i++)
 work[i] = 0.;
 for (octave_idx_type i = b.cidx(j);
 i < b.cidx(j+1); i++)
 work[b.ridx(i)] = b.data(i);
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, 1, tmp_data,
 ldm, pipvt, work, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 // Count non-zeros in work vector and adjust
 // space in retval if needed
 octave_idx_type new_nnz = 0;
 for (octave_idx_type i = 0; i < nr; i++)
 if (work[i] != 0.)
 new_nnz++;
 return retval;
 }
 ComplexMatrix
 SparseMatrix::bsolve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 ComplexMatrix retval;
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 char job = 'L';
 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 // Matrix is not positive definite!! Fall through to
 // unsymmetric banded solver.
 mattype.mark_as_unsymmetric ();
 typ = MatrixType::Banded;
 rcond = 0.0;
 err = 0;
 }
 else
 {
 if (calc_cond)
 {
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dpbcon, DPBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 Complex c = b (i,j);
 Bx[i] = std::real (c);
 Bz[i] = std::imag (c);
 }
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, 1, tmp_data,
 ldm, Bx, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 err = -1;
 break;
 }
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, 1, tmp_data,
 ldm, Bz, b.rows(), err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 err = -1;
 break;
 }
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 }
 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
 octave_idx_type *pipvt = ipvt.fortran_vec ();
 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 ldm, pipvt, err));
 if (err != 0)
 {
 err = -2;
 rcond = 0.0;
 if (sing_handler)
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 if (calc_cond)
 {
 char job = '1';
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dpbcon, DPBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 Complex c = b (i, j);
 Bx[i] = std::real (c);
 Bz[i] = std::imag  (c);
 }
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, 1, tmp_data,
 ldm, pipvt, Bx, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, 1, tmp_data,
 ldm, pipvt, Bz, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::bsolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseComplexMatrix retval;
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 char job = 'L';
 F77_XFCN (dpbtrf, DPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 // Matrix is not positive definite!! Fall through to
 // unsymmetric banded solver.
 mattype.mark_as_unsymmetric ();
 typ = MatrixType::Banded;
 rcond = 0.0;
 err = 0;
 }
 else
 {
 if (calc_cond)
 {
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dpbcon, DPBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, tmp_data, ldm,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 Complex c = b (i,j);
 Bx[i] = std::real (c);
 Bz[i] = std::imag (c);
 }
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, 1, tmp_data,
 ldm, Bx, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 err = -1;
 break;
 }
 F77_XFCN (dpbtrs, DPBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, 1, tmp_data,
 ldm, Bz, b_nr, err
 F77_CHAR_ARG_LEN (1)));
 err = -1;
 break;
 }
 // Count non-zeros in work vector and adjust
 // space in retval if needed
 octave_idx_type new_nnz = 0;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 new_nnz++;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 {
 retval.xridx(ii) = i;
 retval.xdata(ii++) =
 Complex (Bx[i], Bz[i]);
 }
 retval.xcidx(j+1) = ii;
 }
 octave_idx_type n_lower = mattype.nlower ();
 octave_idx_type ldm = n_upper + 2 * n_lower + 1;
 Matrix m_band (ldm, nc);
 double *tmp_data = m_band.fortran_vec ();
 if (! mattype.is_dense ())
 {
 octave_idx_type ii = 0;
 for (octave_idx_type j = 0; j < ldm; j++)
 for (octave_idx_type i = 0; i < nc; i++)
 }
 Array<octave_idx_type> ipvt (dim_vector (nr, 1));
 octave_idx_type *pipvt = ipvt.fortran_vec ();
 F77_XFCN (dgbtrf, DGBTRF, (nr, nr, n_lower, n_upper, tmp_data,
 ldm, pipvt, err));
 if (err != 0)
 {
 err = -2;
 rcond = 0.0;
 if (sing_handler)
 }
 else
 (*current_liboctave_error_handler)
 ("matrix singular to machine precision");
 }
 else
 {
 if (calc_cond)
 {
 char job = '1';
 Array<double> z (dim_vector (3 * nr, 1));
 double *pz = z.fortran_vec ();
 Array<octave_idx_type> iz (dim_vector (nr, 1));
 octave_idx_type *piz = iz.fortran_vec ();
 F77_XFCN (dgbcon, DGBCON,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nc, n_lower, n_upper, tmp_data, ldm, pipvt,
 anorm, rcond, pz, piz, err
 F77_CHAR_ARG_LEN (1)));
 if (err != 0)
 err = -2;
 volatile double rcond_plus_one = rcond + 1.0;
 if (rcond_plus_one == 1.0 || xisnan (rcond))
 for (octave_idx_type i = 0; i < nr; i++)
 {
 Bx[i] = 0.;
 Bz[i] = 0.;
 }
 for (octave_idx_type i = b.cidx(j);
 i < b.cidx(j+1); i++)
 {
 Complex c = b.data(i);
 Bx[b.ridx(i)] = std::real (c);
 Bz[b.ridx(i)] = std::imag (c);
 }
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, 1, tmp_data,
 ldm, pipvt, Bx, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 F77_XFCN (dgbtrs, DGBTRS,
 (F77_CONST_CHAR_ARG2 (&job, 1),
 nr, n_lower, n_upper, 1, tmp_data,
 ldm, pipvt, Bz, b.rows (), err
 F77_CHAR_ARG_LEN (1)));
 // Count non-zeros in work vector and adjust
 // space in retval if needed
 octave_idx_type new_nnz = 0;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 new_nnz++;
 for (octave_idx_type i = 0; i < nr; i++)
 if (Bx[i] != 0. || Bz[i] != 0.)
 {
 retval.xridx(ii) = i;
 retval.xdata(ii++) =
 Complex (Bx[i], Bz[i]);
 }
 retval.xcidx(j+1) = ii;
 }
 }
 }
 else if (typ != MatrixType::Banded_Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 void *
 SparseMatrix::factorize (octave_idx_type& err, double &rcond, Matrix &Control,
 int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, 0,
 &Symbolic, control, info);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve symbolic factorization failed");
 err = -1;
 UMFPACK_DNAME (report_status) (control, status);
 UMFPACK_DNAME (report_info) (control, info);
 rcond = Info (UMFPACK_RCOND);
 else
 rcond = 1.;
 volatile double rcond_plus_one = rcond + 1.0;
 if (status == UMFPACK_WARNING_singular_matrix ||
 rcond_plus_one == 1.0 || xisnan (rcond))
 {
 UMFPACK_DNAME (report_numeric) (Numeric, control);
 err = -2;
 rcond);
 }
 else if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve numeric factorization failed");
 UMFPACK_DNAME (report_status) (control, status);
 UMFPACK_DNAME (report_info) (control, info);
 err = -1;
 }
 else
 {
 UMFPACK_DNAME (report_numeric) (Numeric, control);
 return Numeric;
 }
 Matrix
 SparseMatrix::fsolve (MatrixType &mattype, const Matrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 Matrix retval;
 }
 else
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
 rcond);
 return retval;
 }
 cholmod_dense *X;
 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 if (typ == MatrixType::Full)
 {
 #ifdef HAVE_UMFPACK
 Matrix Control, Info;
 void *Numeric =
 factorize (err, rcond, Control, Info, sing_handler, calc_cond);
 if (err == 0)
 {
 const double *Bx = b.fortran_vec ();
 const octave_idx_type *Ai = ridx ();
 const double *Ax = data ();
 for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr)
 {
 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
 Ai, Ax, &result[iidx], &Bx[iidx],
 Numeric, control, info);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 UMFPACK_DNAME (report_status) (control, status);
 err = -1;
 break;
 }
 }
 UMFPACK_DNAME (report_info) (control, info);
 UMFPACK_DNAME (free_numeric) (&Numeric);
 }
 else
 mattype.mark_as_rectangular ();
 #endif
 }
 else if (typ != MatrixType::Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 SparseMatrix
 SparseMatrix::fsolve (MatrixType &mattype, const SparseMatrix& b,
 }
 else
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
 rcond);
 return retval;
 }
 cholmod_sparse *X;
 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 retval = SparseMatrix (static_cast<octave_idx_type>(X->nrow),
 static_cast<octave_idx_type>(X->ncol),
 static_cast<octave_idx_type>(X->nzmax));
 for (octave_idx_type j = 0;
 j <= static_cast<octave_idx_type>(X->ncol); j++)
 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
 for (octave_idx_type j = 0;
 j < static_cast<octave_idx_type>(X->nzmax); j++)
 {
 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
 retval.xdata(j) = static_cast<double *>(X->x)[j];
 }
 if (typ == MatrixType::Full)
 {
 #ifdef HAVE_UMFPACK
 Matrix Control, Info;
 void *Numeric = factorize (err, rcond, Control, Info,
 sing_handler, calc_cond);
 if (err == 0)
 {
 octave_idx_type b_nr = b.rows ();
 {
 for (octave_idx_type i = 0; i < b_nr; i++)
 Bx[i] = b.elem (i, j);
 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
 Ai, Ax, Xx, Bx, Numeric, control,
 info);
 if (status < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 UMFPACK_DNAME (report_status) (control, status);
 err = -1;
 break;
 }
 for (octave_idx_type i = 0; i < b_nr; i++)
 {
 double tmp = Xx[i];
 if (tmp != 0.0)
 {
 #endif
 }
 else if (typ != MatrixType::Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 ComplexMatrix
 SparseMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 ComplexMatrix retval;
 }
 else
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
 rcond);
 return retval;
 }
 cholmod_dense *X;
 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 if (typ == MatrixType::Full)
 {
 #ifdef HAVE_UMFPACK
 Matrix Control, Info;
 void *Numeric = factorize (err, rcond, Control, Info,
 sing_handler, calc_cond);
 if (err == 0)
 {
 octave_idx_type b_nr = b.rows ();
 retval.resize (b_nr, b_nc);
 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
 for (octave_idx_type j = 0; j < b_nc; j++)
 {
 for (octave_idx_type i = 0; i < b_nr; i++)
 {
 Complex c = b (i,j);
 Bx[i] = std::real (c);
 Bz[i] = std::imag (c);
 }
 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
 Ai, Ax, Xx, Bx, Numeric, control,
 info);
 int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
 Ap, Ai, Ax, Xz, Bz, Numeric,
 control, info) ;
 if (status < 0 || status2 < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 UMFPACK_DNAME (report_status) (control, status);
 err = -1;
 break;
 }
 #endif
 }
 else if (typ != MatrixType::Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 SparseComplexMatrix
 SparseMatrix::fsolve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool calc_cond) const
 {
 SparseComplexMatrix retval;
 }
 else
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve matrix singular to machine precision, rcond = %g",
 rcond);
 return retval;
 }
 cholmod_sparse *X;
 BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm);
 END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
 retval = SparseComplexMatrix
 (static_cast<octave_idx_type>(X->nrow),
 static_cast<octave_idx_type>(X->ncol),
 static_cast<octave_idx_type>(X->nzmax));
 for (octave_idx_type j = 0;
 j <= static_cast<octave_idx_type>(X->ncol); j++)
 retval.xcidx(j) = static_cast<octave_idx_type *>(X->p)[j];
 for (octave_idx_type j = 0;
 j < static_cast<octave_idx_type>(X->nzmax); j++)
 {
 retval.xridx(j) = static_cast<octave_idx_type *>(X->i)[j];
 retval.xdata(j) = static_cast<Complex *>(X->x)[j];
 }
 if (typ == MatrixType::Full)
 {
 #ifdef HAVE_UMFPACK
 Matrix Control, Info;
 void *Numeric = factorize (err, rcond, Control, Info,
 sing_handler, calc_cond);
 if (err == 0)
 {
 octave_idx_type b_nr = b.rows ();
 octave_idx_type ii = 0;
 retval = SparseComplexMatrix (b_nr, b_nc, x_nz);
 OCTAVE_LOCAL_BUFFER (double, Xx, b_nr);
 OCTAVE_LOCAL_BUFFER (double, Xz, b_nr);
 retval.xcidx(0) = 0;
 for (octave_idx_type j = 0; j < b_nc; j++)
 {
 for (octave_idx_type i = 0; i < b_nr; i++)
 {
 Bx[i] = std::real (c);
 Bz[i] = std::imag (c);
 }
 status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap,
 Ai, Ax, Xx, Bx, Numeric, control,
 info);
 int status2 = UMFPACK_DNAME (solve) (UMFPACK_A,
 Ap, Ai, Ax, Xz, Bz, Numeric,
 control, info) ;
 if (status < 0 || status2 < 0)
 {
 (*current_liboctave_error_handler)
 ("SparseMatrix::solve solve failed");
 UMFPACK_DNAME (report_status) (control, status);
 err = -1;
 break;
 }
 #endif
 }
 else if (typ != MatrixType::Hermitian)
 (*current_liboctave_error_handler) ("incorrect matrix type");
 }
 return retval;
 }
 Matrix
 SparseMatrix::solve (MatrixType &mattype, const Matrix& b) const
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 Matrix
 SparseMatrix::solve (MatrixType &mattype, const Matrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 Matrix
 SparseMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
 double& rcond) const
 {
 return solve (mattype, b, info, rcond, 0);
 }
 Matrix
 SparseMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& err,
 double& rcond, solve_singularity_handler sing_handler,
 bool singular_fallback) const
 {
 Matrix retval;
 int typ = mattype.type (false);
 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
 else if (typ != MatrixType::Rectangular)
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 SparseMatrix
 SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 {
 return solve (mattype, b, info, rcond, 0);
 }
 SparseMatrix
 SparseMatrix::solve (MatrixType &mattype, const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool singular_fallback) const
 {
 SparseMatrix retval;
 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
 else if (typ != MatrixType::Rectangular)
 {
 rcond = 1.;
 #ifdef USE_QRSOLVE
 retval = qrsolve (*this, b, err);
 #else
 retval = dmsolve<SparseMatrix, SparseMatrix,
 SparseMatrix> (*this, b, err);
 #endif
 }
 return retval;
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 ComplexMatrix
 SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 ComplexMatrix
 SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& info, double& rcond) const
 {
 return solve (mattype, b, info, rcond, 0);
 }
 ComplexMatrix
 SparseMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool singular_fallback) const
 {
 ComplexMatrix retval;
 int typ = mattype.type (false);
 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
 else if (typ != MatrixType::Rectangular)
 {
 rcond = 1.;
 #ifdef USE_QRSOLVE
 retval = qrsolve (*this, b, err);
 #else
 retval = dmsolve<ComplexMatrix, SparseMatrix,
 ComplexMatrix> (*this, b, err);
 #endif
 }
 return retval;
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 SparseComplexMatrix
 SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 {
 return solve (mattype, b, info, rcond, 0);
 }
 SparseComplexMatrix
 SparseMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler,
 bool singular_fallback) const
 {
 SparseComplexMatrix retval;
 retval = utsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
 retval = ltsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian)
 retval = bsolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Tridiagonal ||
 typ == MatrixType::Tridiagonal_Hermitian)
 retval = trisolve (mattype, b, err, rcond, sing_handler, false);
 else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
 retval = fsolve (mattype, b, err, rcond, sing_handler, true);
 else if (typ != MatrixType::Rectangular)
 {
 rcond = 1.;
 #ifdef USE_QRSOLVE
 retval = qrsolve (*this, b, err);
 #else
 retval = dmsolve<SparseComplexMatrix, SparseMatrix,
 SparseComplexMatrix> (*this, b, err);
 #endif
 }
 return retval;
 double rcond;
 return solve (mattype, b, info, rcond, 0);
 }
 ComplexColumnVector
 SparseMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info,
 double& rcond) const
 {
 return solve (mattype, b, info, rcond, 0);
 }
 double rcond;
 return solve (b, info, rcond, 0);
 }
 Matrix
 SparseMatrix::solve (const Matrix& b, octave_idx_type& info,
 double& rcond) const
 {
 return solve (b, info, rcond, 0);
 }
 Matrix
 SparseMatrix::solve (const Matrix& b, octave_idx_type& err,
 double& rcond,
 solve_singularity_handler sing_handler) const
 {
 MatrixType mattype (*this);
 return solve (mattype, b, err, rcond, sing_handler);
 }
 double rcond;
 return solve (b, info, rcond, 0);
 }
 SparseMatrix
 SparseMatrix::solve (const SparseMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (b, info, rcond, 0);
 }
 {
 return solve (b, info, rcond, 0);
 }
 SparseMatrix
 SparseMatrix::solve (const SparseMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler) const
 {
 MatrixType mattype (*this);
 return solve (mattype, b, err, rcond, sing_handler);
 }
 ComplexMatrix
 SparseMatrix::solve (const ComplexMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (b, info, rcond, 0);
 }
 ComplexMatrix
 SparseMatrix::solve (const ComplexMatrix& b,
 octave_idx_type& info, double& rcond) const
 {
 return solve (b, info, rcond, 0);
 }
 ComplexMatrix
 SparseMatrix::solve (const ComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler) const
 {
 MatrixType mattype (*this);
 return solve (mattype, b, err, rcond, sing_handler);
 }
 double rcond;
 return solve (b, info, rcond, 0);
 }
 SparseComplexMatrix
 SparseMatrix::solve (const SparseComplexMatrix& b,
 octave_idx_type& info) const
 {
 double rcond;
 return solve (b, info, rcond, 0);
 }
 {
 return solve (b, info, rcond, 0);
 }
 SparseComplexMatrix
 SparseMatrix::solve (const SparseComplexMatrix& b,
 octave_idx_type& err, double& rcond,
 solve_singularity_handler sing_handler) const
 {
 MatrixType mattype (*this);
 return solve (mattype, b, err, rcond, sing_handler);
 double rcond;
 return solve (b, info, rcond, 0);
 }
 ComplexColumnVector
 SparseMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info,
 double& rcond) const
 {
 return solve (b, info, rcond, 0);
 }
 }
 return false;
 }
 SparseBoolMatrix
 SparseMatrix::operator ! (void) const
 {
 if (any_element_is_nan ())
 gripe_nan_to_logical_conversion ();
 octave_idx_type nr = rows ();
 octave_idx_type nc = cols ();
 octave_idx_type nz1 = nnz ();
 octave_idx_type nz2 = nr*nc - nz1;
 SparseBoolMatrix r (nr, nc, nz2);
 octave_idx_type ii = 0;
 octave_idx_type jj = 0;
 r.cidx (0) = 0;
 for (octave_idx_type i = 0; i < nc; i++)
 {
 {
 if ((rows() == 1 && dim == -1) || dim == 1)
 return transpose (). prod (0). transpose();
 else
 {
 SPARSE_REDUCTION_OP (SparseMatrix, double, *=,
 (cidx(j+1) - cidx(j) < nr ? 0.0 : 1.0), 1.0);
 }
 }
 SparseMatrix
 #define COL_EXPR \
 double d = data (i); \
 tmp[j] += d * d
 SPARSE_BASE_REDUCTION_OP (SparseMatrix, double, ROW_EXPR, COL_EXPR,
 0.0, 0.0);
 #undef ROW_EXPR
 #undef COL_EXPR
 }
 return read_sparse_matrix<elt_type> (is, a, octave_read_value<double>);
 }
 SparseMatrix
 SparseMatrix::squeeze (void) const
 {
 return MSparse<double>::squeeze ();
 }
 SparseMatrix
 SparseMatrix::reshape (const dim_vector& new_dims) const
 {
 SparseMatrix
 min (const SparseMatrix& a, const SparseMatrix& b)
 {
 SparseMatrix r;
 if ((a.rows() == b.rows()) && (a.cols() == b.cols()))
 {
 octave_idx_type a_nr = a.rows ();
 octave_idx_type a_nc = a.cols ();
 octave_idx_type b_nr = b.rows ();
 if (a_nr != b_nr || a_nc != b_nc)
 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
 else
 {
 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
 octave_idx_type jx = 0;
 r.cidx (0) = 0;
 for (octave_idx_type i = 0 ; i < a_nc ; i++)
 {
 octave_idx_type  ja = a.cidx(i);
 octave_idx_type  ja_max = a.cidx(i+1);
 bool ja_lt_max= ja < ja_max;
 octave_idx_type  jb = b.cidx(i);
 octave_idx_type  jb_max = b.cidx(i+1);
 bool jb_lt_max = jb < jb_max;
 while (ja_lt_max || jb_lt_max )
 {
 octave_quit ();
 if ((! jb_lt_max) ||
 (ja_lt_max && (a.ridx(ja) < b.ridx(jb))))
 jb_lt_max= jb < jb_max;
 }
 }
 r.cidx(i+1) = jx;
 }
 r.maybe_compress ();
 }
 }
 else
 (*current_liboctave_error_handler) ("matrix size mismatch");
 SparseMatrix
 max (const SparseMatrix& a, const SparseMatrix& b)
 {
 SparseMatrix r;
 if ((a.rows() == b.rows()) && (a.cols() == b.cols()))
 {
 octave_idx_type a_nr = a.rows ();
 octave_idx_type a_nc = a.cols ();
 octave_idx_type b_nr = b.rows ();
 if (a_nr != b_nr || a_nc != b_nc)
 gripe_nonconformant ("min", a_nr, a_nc, b_nr, b_nc);
 else
 {
 r = SparseMatrix (a_nr, a_nc, (a.nnz () + b.nnz ()));
 octave_idx_type jx = 0;
 r.cidx (0) = 0;
 for (octave_idx_type i = 0 ; i < a_nc ; i++)
 {
 octave_idx_type  ja = a.cidx(i);
 octave_idx_type  ja_max = a.cidx(i+1);
 bool ja_lt_max= ja < ja_max;
 octave_idx_type  jb = b.cidx(i);
 octave_idx_type  jb_max = b.cidx(i+1);
 bool jb_lt_max = jb < jb_max;
 while (ja_lt_max || jb_lt_max )
 {
 octave_quit ();
 if ((! jb_lt_max) ||
 (ja_lt_max && (a.ridx(ja) < b.ridx(jb))))
 jb_lt_max= jb < jb_max;
 }
 }
 r.cidx(i+1) = jx;
 }
 r.maybe_compress ();
 }
 }
 else
 (*current_liboctave_error_handler) ("matrix size mismatch");

Mercurial > octave-nkf

comparison liboctave/dSparse.cc @ 11586:12df7854fa7c