Mercurial > octave-nkf
changeset 5322:22994a5730f9
[project @ 2005-04-29 13:04:24 by dbateman]
line wrap: on
line diff
--- a/ChangeLog Fri Apr 29 05:20:01 2005 +0000 +++ b/ChangeLog Fri Apr 29 13:04:25 2005 +0000 @@ -1,3 +1,7 @@ +2005-04-29 David Bateman <dbateman@free.fr> + + * configure.in: Add UMFPACK_LONG_IDX + 2005-04-21 John W. Eaton <jwe@octave.org> * configure.in (AC_CONFIG_FILES): Remove install-octave from the list.
--- a/PROJECTS Fri Apr 29 05:20:01 2005 +0000 +++ b/PROJECTS Fri Apr 29 13:04:25 2005 +0000 @@ -113,8 +113,8 @@ * QR factorization functions, also for use in lssolve functions. Write svds function based on this. Write sprank function based on svds. - * Once dmperm is implemented, use the technique to detect permuted - triangular matrices. Test the permuted triangular matrix solver code + * Once dmperm is implemented, perhaps use the technique to detect permuted + triangular matrices with permutations in both rows and cols. * Sparse inverse function, based on Andy's code from octave-forge. @@ -160,7 +160,8 @@ way to treat this but a complete rewrite of the sparse assignment functions... - * Port the sparse testing code into the DejaGNU testing code. + * Port the sparse testing code into the DejaGNU testing code. Add tests + for spkron, matrix_type, luinc.. * Treat the dispatching of the splu, spdet, functions, etc within octave itself. This either means that classes need implementing or at the
--- a/configure.in Fri Apr 29 05:20:01 2005 +0000 +++ b/configure.in Fri Apr 29 13:04:25 2005 +0000 @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.474 $) +AC_REVISION($Revision: 1.475 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -160,6 +160,7 @@ OCTAVE_IDX_TYPE=int elif test $ac_cv_sizeof_long -eq 8; then OCTAVE_IDX_TYPE=long + AC_DEFINE(UMFPACK_LONG_IDX, 1, [Define to 1 to use long int versions of UMFPACK]) else AC_MSG_WARN([no suitable type found for octave_idx_type so disabling 64-bit features]) USE_64_BIT_IDX_T=false
--- a/doc/ChangeLog Fri Apr 29 05:20:01 2005 +0000 +++ b/doc/ChangeLog Fri Apr 29 13:04:25 2005 +0000 @@ -1,3 +1,8 @@ +2005-04-29 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Add matrix_type, spkron, and document changes + in solve code + 2005-03-14 David Bateman <dbateman@free.fr> * interpreter/sparse.txi: Add luinc function.
--- a/doc/interpreter/sparse.txi Fri Apr 29 05:20:01 2005 +0000 +++ b/doc/interpreter/sparse.txi Fri Apr 29 13:04:25 2005 +0000 @@ -415,9 +415,10 @@ @sc{Octave} includes a poly-morphic solver for sparse matrices, where the exact solver used to factorize the matrix, depends on the properties -of the sparse matrix, itself. As the cost of determining the matrix type -is small relative to the cost of factorizing the matrix itself, the matrix -type is re-determined each time it is used in a linear equation. +of the sparse matrix, itself. The cost of determining the matrix type +is small relative to the cost of factorizing the matrix itself, but in any +case the matrix type is cached once it is calculated, so that it is not +re-determined each time it is used in a linear equation. The selection tree for how the linear equation is solve is @@ -456,11 +457,9 @@ @item If the matrix is upper or lower triangular perform a sparse forward or backward subsitution, and goto 9 -@item If the matrix is a permuted upper or lower triangular matrix, perform -a sparse forward or backward subsitution, and goto 9 - -FIXME: Detection of permuted triangular matrices not written yet, and so - the code itself is not tested either +@item If the matrix is a upper triangular matrix with column permutations +or lower triangular matrix with row permutations, perform a sparse forward +or backward subsitution, and goto 9 @item If the matrix is hermitian with a real positive diagonal, attempt sparse Cholesky factorization. @@ -494,6 +493,12 @@ In cases where, this might be a problem the user is recommended to disable the banded solvers as above, at a significant cost in terms of speed. +The user can force the type of the matrix with the @code{matrix_type} +function. This overcomes the cost of discovering the type of the matrix. +However, it should be noted incorrectly identifying the type of the matrix +will lead to unpredictable results, and so @code{matrix_type} should be +use dwith care. + @node Iterative Techniques, Oct-Files, Sparse Linear Algebra, Sparse Matrices @section Iterative Techniques applied to sparse matrices @@ -965,7 +970,7 @@ @item colperm Returns the column permutations such that the columns of `S (:, P)' are ordered in terms of increase number of non-zero elements. @item dmperm -Perfrom a Deulmage-Mendelsohn permutation on the sparse matrix S. +Perform a Deulmage-Mendelsohn permutation on the sparse matrix S. @item symamd For a symmetric positive definite matrix S, returns the permutation vector p such that `S (P, P)' tends to have a sparser Cholesky factor than S. @item symrcm @@ -979,12 +984,16 @@ @emph{Not implemented} @item eigs @emph{Not implemented} +@item matrix_type +Identify the matrix type or mark a matrix as a particular type. @item normest @emph{Not implemented} @item spdet Compute the determinant of sparse matrix A using UMFPACK. @item spinv Compute the inverse of the square matrix A. +@item spkron +Form the kronecker product of two sparse matrices. @item splu Compute the LU decomposition of the sparse matrix A, using subroutines from UMFPACK. @item sprank @@ -1066,6 +1075,8 @@ matrix. * luinc:: Produce the incomplete LU factorization of the sparse A. +* matrix_type:: Identify the matrix type or mark a matrix as a particular + type. * nnz:: returns number of non zero elements in SM See also: sparse * nonzeros:: Returns a vector of the non-zero values of the sparse matrix S @@ -1089,6 +1100,7 @@ * spfun:: Compute `f(X)' for the non-zero values of X This results in a sparse matrix with the same structure as X. * spinv:: Compute the inverse of the square matrix A. +* spkron:: Form the kronecker product of two sparse matrices. * splu:: Compute the LU decomposition of the sparse matrix A, using subroutines from UMFPACK. * spmax:: For a vector argument, return the maximum value. @@ -1147,12 +1159,17 @@ @DOCSTRING(issparse) -@node luinc, nnz, issparse, Function Reference +@node luinc, matrix_type, issparse, Function Reference @subsubsection luinc @DOCSTRING(luinc) -@node nnz, nonzeros, luinc, Function Reference +@node matrix_type, nnz, luinc, Function Reference +@subsubsection matrix_type + +@DOCSTRING(matrix_type) + +@node nnz, nonzeros, matrix_type, Function Reference @subsubsection nnz @DOCSTRING(nnz) @@ -1227,12 +1244,17 @@ @DOCSTRING(spfun) -@node spinv, splu, spfun, Function Reference +@node spinv, spkron, spfun, Function Reference @subsubsection spinv @DOCSTRING(spinv) -@node splu, spmax, spinv, Function Reference +@node spkron, splu, spinv, Function Reference +@subsubsection spkron + +@DOCSTRING(spkron) + +@node splu, spmax, spkron, Function Reference @subsubsection splu @DOCSTRING(splu)
--- a/liboctave/CSparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/CSparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -98,7 +98,7 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*, + F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, double*, Complex*, Complex*, const octave_idx_type&, octave_idx_type&); F77_RET_T @@ -649,7 +649,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_zi_defaults (control); + UMFPACK_ZNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -670,20 +670,20 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_zi_report_control (control); + UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); - umfpack_zi_report_matrix (nr, nc, Ap, Ai, - X_CAST (const double *, Ax), - NULL, 1, control); + UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, + X_CAST (const double *, Ax), + NULL, 1, control); void *Symbolic; Matrix Info (1, UMFPACK_INFO); double *info = Info.fortran_vec (); - int status = umfpack_zi_qsymbolic + int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, NULL, &Symbolic, control, info); @@ -692,20 +692,20 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant symbolic factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); - - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_zi_report_symbolic (Symbolic, control); + UMFPACK_ZNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), - NULL, Symbolic, &Numeric, - control, info) ; - umfpack_zi_free_symbolic (&Symbolic) ; + status = UMFPACK_ZNAME (numeric) (Ap, Ai, + X_CAST (const double *, Ax), NULL, + Symbolic, &Numeric, control, info) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; rcond = Info (UMFPACK_RCOND); @@ -714,19 +714,19 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant numeric factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); - - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_numeric) (&Numeric); } else { - umfpack_zi_report_numeric (Numeric, control); + UMFPACK_ZNAME (report_numeric) (Numeric, control); Complex d[2]; double d_exponent; - status = umfpack_zi_get_determinant + status = UMFPACK_ZNAME (get_determinant) (X_CAST (double *, &d[0]), NULL, &d_exponent, Numeric, info); d[1] = d_exponent; @@ -736,10 +736,10 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant error calculating determinant"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); } else retval = ComplexDET (d); @@ -1052,13 +1052,9 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (b.rows (), b.cols ()); + retval.resize (nr, b_cols); + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); for (octave_idx_type j = 0; j < b_cols; j++) { @@ -1067,28 +1063,29 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; + retval (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -1097,19 +1094,20 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -1269,12 +1267,12 @@ if (typ == SparseType::Permuted_Upper) { + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { @@ -1285,22 +1283,23 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -1321,10 +1320,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[rperm[i]]; } retval.xcidx(j+1) = ii; } @@ -1337,19 +1336,20 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -1526,13 +1526,9 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (b.rows (), b.cols ()); + retval.resize (nr, b_nc); + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); for (octave_idx_type j = 0; j < b_nc; j++) { @@ -1541,29 +1537,29 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -1572,19 +1568,20 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -1744,12 +1741,12 @@ if (typ == SparseType::Permuted_Upper) { + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { @@ -1760,22 +1757,23 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -1796,10 +1794,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[rperm[i]]; } retval.xcidx(j+1) = ii; } @@ -1812,19 +1810,20 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + Complex 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -2003,42 +2002,48 @@ { retval.resize (b.rows (), b.cols ()); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_cols; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = b(i,j); + work[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -2047,25 +2052,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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; + } + + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2221,37 +2238,44 @@ if (typ == SparseType::Permuted_Lower) { OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (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); + work[perm[b.ridx(i)]] = b.data(i); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -2272,10 +2296,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[i]; } retval.xcidx(j+1) = ii; } @@ -2288,25 +2312,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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; + } + + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2482,42 +2518,48 @@ { retval.resize (b.rows (), b.cols ()); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = b(i,j); + work[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -2526,25 +2568,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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; + } + + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2701,37 +2755,44 @@ if (typ == SparseType::Permuted_Lower) { OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseComplexMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (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); + work[perm[b.ridx(i)]] = b.data(i); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -2752,10 +2813,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[i]; } retval.xcidx(j+1) = ii; } @@ -2768,25 +2829,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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; + } + + Complex tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2942,7 +3015,7 @@ if (typ == SparseType::Tridiagonal_Hermitian) { - OCTAVE_LOCAL_BUFFER (Complex, D, nr); + OCTAVE_LOCAL_BUFFER (double, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) @@ -2951,11 +3024,11 @@ for (octave_idx_type j = 0; j < nc-1; j++) { - D[j] = data(ii++); + D[j] = std::real(data(ii++)); DL[j] = data(ii); ii += 2; } - D[nc-1] = data(ii); + D[nc-1] = std::real(data(ii)); } else { @@ -2970,7 +3043,7 @@ for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) { if (ridx(i) == j) - D[j] = data(i); + D[j] = std::real(data(i)); else if (ridx(i) == j + 1) DL[j] = data(i); } @@ -3032,7 +3105,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3129,7 +3202,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3238,10 +3311,9 @@ volatile int typ = mattype.type (); mattype.info (); - // Note can't treat symmetric case as there is no dpttrf function if (typ == SparseType::Tridiagonal_Hermitian) { - OCTAVE_LOCAL_BUFFER (Complex, D, nr); + OCTAVE_LOCAL_BUFFER (double, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) @@ -3250,11 +3322,11 @@ for (octave_idx_type j = 0; j < nc-1; j++) { - D[j] = data(ii++); + D[j] = std::real(data(ii++)); DL[j] = data(ii); ii += 2; } - D[nc-1] = data(ii); + D[nc-1] = std::real(data(ii)); } else { @@ -3269,7 +3341,7 @@ for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) { if (ridx(i) == j) - D[j] = data(i); + D[j] = std::real (data(i)); else if (ridx(i) == j + 1) DL[j] = data(i); } @@ -3335,7 +3407,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3435,7 +3507,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -4457,7 +4529,7 @@ // Setup the control parameters Control = Matrix (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_zi_defaults (control); + UMFPACK_ZNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -4474,7 +4546,7 @@ if (!xisnan (tmp)) Control (UMFPACK_FIXQ) = tmp; - umfpack_zi_report_control (control); + UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); @@ -4482,13 +4554,13 @@ octave_idx_type nr = rows (); octave_idx_type nc = cols (); - umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), - NULL, 1, control); + UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, + X_CAST (const double *, Ax), NULL, 1, control); void *Symbolic; Info = Matrix (1, UMFPACK_INFO); double *info = Info.fortran_vec (); - int status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, + int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, NULL, &Symbolic, control, info); @@ -4498,18 +4570,19 @@ ("SparseComplexMatrix::solve symbolic factorization failed"); err = -1; - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); - - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_zi_report_symbolic (Symbolic, control); - - status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL, + UMFPACK_ZNAME (report_symbolic) (Symbolic, control); + + status = UMFPACK_ZNAME (numeric) (Ap, Ai, + X_CAST (const double *, Ax), NULL, Symbolic, &Numeric, control, info) ; - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; #ifdef HAVE_LSSOLVE rcond = Info (UMFPACK_RCOND); @@ -4518,7 +4591,7 @@ if (status == UMFPACK_WARNING_singular_matrix || rcond_plus_one == 1.0 || xisnan (rcond)) { - umfpack_zi_report_numeric (Numeric, control); + UMFPACK_ZNAME (report_numeric) (Numeric, control); err = -2; @@ -4537,19 +4610,19 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::solve numeric factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); err = -1; } else { - umfpack_zi_report_numeric (Numeric, control); + UMFPACK_ZNAME (report_numeric) (Numeric, control); } } if (err != 0) - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); #else (*current_liboctave_error_handler) ("UMFPACK not installed"); #endif @@ -4620,8 +4693,8 @@ for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { #ifdef UMFPACK_SEPARATE_SPLIT - status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, - X_CAST (const double *, Ax), + status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, + Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, &Xx[iidx]), NULL, @@ -4631,8 +4704,8 @@ for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = b.elem (i, j); - status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, - X_CAST (const double *, Ax), + status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, + Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, &Xx[iidx]), NULL, @@ -4646,7 +4719,7 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, status); err = -1; @@ -4673,9 +4746,9 @@ } #endif - umfpack_zi_report_info (control, info); - - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -4762,8 +4835,8 @@ for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b.elem (i, j); - status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, - X_CAST (const double *, Ax), + status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, + Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, Xx), NULL, Bx, Bz, Numeric, control, @@ -4772,7 +4845,7 @@ for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = b.elem (i, j); - status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, + status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, Xx), NULL, @@ -4785,7 +4858,7 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, status); err = -1; @@ -4833,9 +4906,9 @@ } #endif - umfpack_zi_report_info (control, info); - - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -4904,7 +4977,7 @@ for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { status = - umfpack_zi_solve (UMFPACK_A, Ap, Ai, + UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, &Xx[iidx]), NULL, X_CAST (const double *, &Bx[iidx]), @@ -4915,7 +4988,7 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, status); err = -1; @@ -4942,9 +5015,9 @@ } #endif - umfpack_zi_report_info (control, info); - - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -5022,8 +5095,8 @@ for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b (i,j); - status = umfpack_zi_solve (UMFPACK_A, Ap, Ai, - X_CAST (const double *, Ax), + status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, + Ai, X_CAST (const double *, Ax), NULL, X_CAST (double *, Xx), NULL, X_CAST (double *, Bx), NULL, Numeric, control, info); @@ -5033,7 +5106,7 @@ (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, status); err = -1; @@ -5081,9 +5154,9 @@ } #endif - umfpack_zi_report_info (control, info); - - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (report_info) (control, info); + + UMFPACK_ZNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -5124,7 +5197,7 @@ double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5178,7 +5251,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5232,7 +5305,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5287,7 +5360,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this);
--- a/liboctave/CSparse.h Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/CSparse.h Fri Apr 29 13:04:25 2005 +0000 @@ -423,6 +423,12 @@ SPARSE_FORWARD_DEFS (MSparse, SparseComplexMatrix, ComplexMatrix, Complex) +#ifdef UMFPACK_LONG_IDX +#define UMFPACK_ZNAME(name) umfpack_zl_ ## name +#else +#define UMFPACK_ZNAME(name) umfpack_zi_ ## name +#endif + #endif /*
--- a/liboctave/ChangeLog Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/ChangeLog Fri Apr 29 13:04:25 2005 +0000 @@ -1,3 +1,43 @@ +2005-04-29 David Bateman <dbateman@free.fr> + + * dSparse.cc (trisolve): Diagonal passed to lapack zptsv is type double. + Correct indexing for upper diagonal elements for sparse tridiagonal. + * CSparse.cc (trisolve): ditto. + + * CSparse.h (UMFPACK_ZNAME): Define macro to pick version of UMFPACK for + 64-bit. + * CSparse.cc (UMFPACK_ZNAME): Replace all umfpack_zi_* with + UMFPACK_ZNAME(*). + * SparseCmplxLU.cc (UMFPACK_ZNAME): ditto + + * dSparse.h (UMFPACK_DNAME): Define macro to pick version of UMFPACK for + 64-bit. + * dSparse.cc (UMFPACK_DNAME): Replace all umfpack_di_* with + UMFPACK_DNAME(*). + * SparsedbleLU.cc (UMFPACK_DNAME): ditto + + * dSparse.cc (ltsolve, utsolve): Correct permuted upper/lower triangular + back/forward substitution code. + * CSparse.cc (ltsolve, utsolve): ditto. + + * dSparse.cc (solve): Use mattype.type (false) to force messaging from + spparms("spumoni",1). + * CSparse.cc (solve): ditto + + * SparseType.cc (SparseType(void)): Print info for spparms("spumoni",1). + (SparseType(const matrix_type), SparseType(const matrix_type, const + octave_idx_type, const octave_idx_type*), SparseType(const matrix_type, + const octave_idx_type, const octave_idx_type)): New constructors. + (SparseType (const SparseMatrix&), SparseType (SparseComplexMatrix&)): + Detect row permuted lower triangular and column permuted upper triangular + matrices. Remove one of the permutation vectors.. + + * SparseType.h: Simplify the permutation code. + (SparseType(const matrix_type), SparseType + (const matrix_type, const octave_idx_type, const octave_idx_type*), + SparseType(const matrix_type, const octave_idx_type, + const octave_idx_type)): Declarations. + 2005-04-25 John W. Eaton <jwe@octave.org> * str-vec.cc (string_vector::delete_c_str_vec): Correctly free
--- a/liboctave/SparseCmplxLU.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/SparseCmplxLU.cc Fri Apr 29 13:04:25 2005 +0000 @@ -55,7 +55,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_zi_defaults (control); + UMFPACK_ZNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -84,19 +84,19 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_zi_report_control (control); + UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = a.cidx (); const octave_idx_type *Ai = a.ridx (); const Complex *Ax = a.data (); - umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), - NULL, 1, control); + UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, + X_CAST (const double *, Ax), NULL, 1, control); void *Symbolic; Matrix Info (1, UMFPACK_INFO); double *info = Info.fortran_vec (); - int status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, + int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, NULL, &Symbolic, control, info); @@ -105,19 +105,20 @@ (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_zi_report_symbolic (Symbolic, control); + UMFPACK_ZNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL, + status = UMFPACK_ZNAME (numeric) (Ap, Ai, + X_CAST (const double *, Ax), NULL, Symbolic, &Numeric, control, info) ; - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; cond = Info (UMFPACK_RCOND); @@ -126,94 +127,93 @@ (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU numeric factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); } else { - umfpack_zi_report_numeric (Numeric, control); + UMFPACK_ZNAME (report_numeric) (Numeric, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, - &ignore3, Numeric) ; + octave_idx_type lnz, unz, ignore1, ignore2, ignore3; + status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, &ignore1, + &ignore2, &ignore3, Numeric) ; if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); } else { - int n_inner = (nr < nc ? nr : nc); + octave_idx_type n_inner = (nr < nc ? nr : nc); if (lnz < 1) - Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, + Lfact = SparseComplexMatrix (n_inner, nr, static_cast<octave_idx_type> (1)); else - Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); + Lfact = SparseComplexMatrix (n_inner, nr, lnz); octave_idx_type *Ltp = Lfact.cidx (); octave_idx_type *Ltj = Lfact.ridx (); Complex *Ltx = Lfact.data (); if (unz < 1) - Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, + Ufact = SparseComplexMatrix (n_inner, nc, static_cast<octave_idx_type> (1)); else - Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, unz); + Ufact = SparseComplexMatrix (n_inner, nc, unz); octave_idx_type *Up = Ufact.cidx (); octave_idx_type *Uj = Ufact.ridx (); Complex *Ux = Ufact.data (); P.resize (nr); - int *p = P.fortran_vec (); + octave_idx_type *p = P.fortran_vec (); Q.resize (nc); - int *q = Q.fortran_vec (); + octave_idx_type *q = Q.fortran_vec (); - int do_recip; - status = umfpack_zi_get_numeric (Ltp, Ltj, X_CAST (double *, Ltx), - NULL, Up, Uj, + octave_idx_type do_recip; + status = UMFPACK_ZNAME (get_numeric) (Ltp, Ltj, + X_CAST (double *, Ltx), NULL, Up, Uj, X_CAST (double *, Ux), NULL, p, q, NULL, NULL, &do_recip, NULL, Numeric) ; - umfpack_zi_free_numeric (&Numeric) ; + UMFPACK_ZNAME (free_numeric) (&Numeric) ; if (status < 0 || do_recip) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, status); } else { Lfact = Lfact.transpose (); - umfpack_zi_report_matrix (nr, n_inner, Lfact.cidx (), - Lfact.ridx (), + UMFPACK_ZNAME (report_matrix) (nr, n_inner, + Lfact.cidx (), Lfact.ridx (), X_CAST (double *, Lfact.data()), NULL, 1, control); - umfpack_zi_report_matrix (n_inner, nc, Ufact.cidx (), - Ufact.ridx (), + UMFPACK_ZNAME (report_matrix) (n_inner, nc, + Ufact.cidx (), Ufact.ridx (), X_CAST (double *, Ufact.data()), NULL, 1, control); - umfpack_zi_report_perm (nr, p, control); - umfpack_zi_report_perm (nc, q, control); + UMFPACK_ZNAME (report_perm) (nr, p, control); + UMFPACK_ZNAME (report_perm) (nc, q, control); } - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_info) (control, info); } } } @@ -239,7 +239,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_zi_defaults (control); + UMFPACK_ZNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -276,13 +276,13 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_zi_report_control (control); + UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = a.cidx (); const octave_idx_type *Ai = a.ridx (); const Complex *Ax = a.data (); - umfpack_zi_report_matrix (nr, nc, Ap, Ai, + UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, 1, control); @@ -294,12 +294,12 @@ // Null loop so that qinit is imediately deallocated when not // needed do { - OCTAVE_LOCAL_BUFFER (int, qinit, nc); + OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc); - for (int i = 0; i < nc; i++) - qinit [i] = static_cast<int> (Qinit (i)); + for (octave_idx_type i = 0; i < nc; i++) + qinit [i] = static_cast<octave_idx_type> (Qinit (i)); - status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, + status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, qinit, &Symbolic, control, info); @@ -310,20 +310,20 @@ (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_zi_report_symbolic (Symbolic, control); + UMFPACK_ZNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_zi_numeric (Ap, Ai, + status = UMFPACK_ZNAME (numeric) (Ap, Ai, X_CAST (const double *, Ax), NULL, Symbolic, &Numeric, control, info) ; - umfpack_zi_free_symbolic (&Symbolic) ; + UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; cond = Info (UMFPACK_RCOND); @@ -332,102 +332,98 @@ (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU numeric factorization failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); } else { - umfpack_zi_report_numeric (Numeric, control); + UMFPACK_ZNAME (report_numeric) (Numeric, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, - &ignore2, &ignore3, Numeric); + octave_idx_type lnz, unz, ignore1, ignore2, ignore3; + status = UMFPACK_ZNAME (get_lunz) (&lnz, &unz, + &ignore1, &ignore2, &ignore3, Numeric); if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_status) (control, status); + UMFPACK_ZNAME (report_info) (control, info); - umfpack_zi_free_numeric (&Numeric); + UMFPACK_ZNAME (free_numeric) (&Numeric); } else { - int n_inner = (nr < nc ? nr : nc); + octave_idx_type n_inner = (nr < nc ? nr : nc); if (lnz < 1) - Lfact = SparseComplexMatrix - (static_cast<octave_idx_type> (n_inner), nr, + Lfact = SparseComplexMatrix (n_inner, nr, static_cast<octave_idx_type> (1)); else - Lfact = SparseComplexMatrix - (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); + Lfact = SparseComplexMatrix (n_inner, nr, lnz); octave_idx_type *Ltp = Lfact.cidx (); octave_idx_type *Ltj = Lfact.ridx (); Complex *Ltx = Lfact.data (); if (unz < 1) - Ufact = SparseComplexMatrix - (static_cast<octave_idx_type> (n_inner), nc, + Ufact = SparseComplexMatrix (n_inner, nc, static_cast<octave_idx_type> (1)); else - Ufact = SparseComplexMatrix - (static_cast<octave_idx_type> (n_inner), nc, unz); + Ufact = SparseComplexMatrix (n_inner, nc, unz); octave_idx_type *Up = Ufact.cidx (); octave_idx_type *Uj = Ufact.ridx (); Complex *Ux = Ufact.data (); P.resize (nr); - int *p = P.fortran_vec (); + octave_idx_type *p = P.fortran_vec (); Q.resize (nc); - int *q = Q.fortran_vec (); + octave_idx_type *q = Q.fortran_vec (); - int do_recip; + octave_idx_type do_recip; status = - umfpack_zi_get_numeric (Ltp, Ltj, + UMFPACK_ZNAME (get_numeric) (Ltp, Ltj, X_CAST (double *, Ltx), NULL, Up, Uj, X_CAST (double *, Ux), NULL, p, q, NULL, NULL, &do_recip, NULL, Numeric) ; - umfpack_zi_free_numeric (&Numeric) ; + UMFPACK_ZNAME (free_numeric) (&Numeric) ; if (status < 0 || do_recip) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); - umfpack_zi_report_status (control, status); + UMFPACK_ZNAME (report_status) (control, + status); } else { Lfact = Lfact.transpose (); - umfpack_zi_report_matrix (nr, n_inner, + UMFPACK_ZNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (), X_CAST (double *, Lfact.data()), NULL, 1, control); - umfpack_zi_report_matrix (n_inner, nc, + UMFPACK_ZNAME (report_matrix) (n_inner, nc, Ufact.cidx (), Ufact.ridx (), X_CAST (double *, Ufact.data()), NULL, 1, control); - umfpack_zi_report_perm (nr, p, control); - umfpack_zi_report_perm (nc, q, control); + UMFPACK_ZNAME (report_perm) (nr, p, control); + UMFPACK_ZNAME (report_perm) (nc, q, control); } - umfpack_zi_report_info (control, info); + UMFPACK_ZNAME (report_info) (control, info); } } }
--- a/liboctave/SparseType.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/SparseType.cc Fri Apr 29 13:04:25 2005 +0000 @@ -24,6 +24,8 @@ #include <config.h> #endif +#include <vector> + #include "SparseType.h" #include "dSparse.h" #include "CSparse.h" @@ -31,6 +33,11 @@ // XXX FIXME XXX There is a large code duplication here +SparseType::SparseType (void) : typ (SparseType::Unknown), nperm (0) +{ + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); +} + SparseType::SparseType (const SparseType &a) : typ (a.typ), sp_bandden (a.sp_bandden), bandden (a.bandden), upper_band (a.upper_band), lower_band (a.lower_band), @@ -38,13 +45,9 @@ { if (nperm != 0) { - row_perm = new octave_idx_type [nperm]; - col_perm = new octave_idx_type [nperm]; + perm = new octave_idx_type [nperm]; for (octave_idx_type i = 0; i < nperm; i++) - { - row_perm[i] = a.row_perm[i]; - col_perm[i] = a.col_perm[i]; - } + perm[i] = a.perm[i]; } } @@ -54,6 +57,10 @@ octave_idx_type ncols = a.cols (); octave_idx_type nnz = a.nnz (); + if (Voctave_sparse_controls.get_key ("spumoni") != 0.) + (*current_liboctave_warning_handler) + ("Calculating Sparse Matrix Type"); + nperm = 0; if (nrows != ncols) @@ -177,8 +184,76 @@ // Search for a permuted triangular matrix, and test if // permutation is singular - // XXX FIXME XXX Write this test based on dmperm - + // XXX FIXME XXX Perhaps this should be based on a dmperm algorithm + bool found = false; + + nperm = nrows; + perm = new octave_idx_type [nperm]; + + for (octave_idx_type i = 0; i < nperm; i++) + { + found = false; + + for (octave_idx_type j = 0; j < ncols; j++) + { + + if ((a.cidx(j+1) - a.cidx(j)) > 0 && + a.ridx(a.cidx(j+1)-1) == i) + { + perm [i] = j; + found = true; + break; + } + } + + if (!found) + break; + } + + if (found) + typ = SparseType::Permuted_Upper; + else + { + OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nperm); + + for (octave_idx_type i = 0; i < nperm; i++) + { + perm [i] = -1; + tmp [i] = -1; + } + + for (octave_idx_type j = 0; j < ncols; j++) + for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) + perm [a.ridx(i)] = j; + + found = true; + for (octave_idx_type i = 0; i < nperm; i++) + if (perm[i] == -1) + { + found = false; + break; + } + else + { + tmp[perm[i]] = 1; + } + + if (found) + for (octave_idx_type i = 0; i < nperm; i++) + if (tmp[i] == -1) + { + found = false; + break; + } + + if (found) + typ = SparseType::Permuted_Lower; + else + { + delete [] perm; + nperm = 0; + } + } } } @@ -253,6 +328,10 @@ octave_idx_type ncols = a.cols (); octave_idx_type nnz = a.nnz (); + if (Voctave_sparse_controls.get_key ("spumoni") != 0.) + (*current_liboctave_warning_handler) + ("Calculating Sparse Matrix Type"); + nperm = 0; if (nrows != ncols) @@ -376,8 +455,76 @@ // Search for a permuted triangular matrix, and test if // permutation is singular - // XXX FIXME XXX Write this test based on dmperm - + // XXX FIXME XXX Perhaps this should be based on a dmperm algorithm + bool found = false; + + nperm = nrows; + perm = new octave_idx_type [nperm]; + + for (octave_idx_type i = 0; i < nperm; i++) + { + found = false; + + for (octave_idx_type j = 0; j < ncols; j++) + { + + if ((a.cidx(j+1) - a.cidx(j)) > 0 && + a.ridx(a.cidx(j+1)-1) == i) + { + perm [i] = j; + found = true; + break; + } + } + + if (!found) + break; + } + + if (found) + typ = SparseType::Permuted_Upper; + else + { + OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nperm); + + for (octave_idx_type i = 0; i < nperm; i++) + { + perm [i] = -1; + tmp [i] = -1; + } + + for (octave_idx_type j = 0; j < ncols; j++) + for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) + perm [a.ridx(i)] = j; + + found = true; + for (octave_idx_type i = 0; i < nperm; i++) + if (perm[i] == -1) + { + found = false; + break; + } + else + { + tmp[perm[i]] = 1; + } + + if (found) + for (octave_idx_type i = 0; i < nperm; i++) + if (tmp[i] == -1) + { + found = false; + break; + } + + if (found) + typ = SparseType::Permuted_Lower; + else + { + delete [] perm; + nperm = 0; + } + } } } @@ -446,12 +593,59 @@ } } +SparseType::SparseType (const matrix_type t) : typ (SparseType::Unknown), + nperm (0) +{ + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); + + if (t == SparseType::Full || t == SparseType::Diagonal || + t == SparseType::Permuted_Diagonal || t == SparseType::Upper || + t == SparseType::Lower || t == SparseType::Tridiagonal || + t == SparseType::Tridiagonal_Hermitian || t == SparseType::Rectangular) + typ = t; + else + (*current_liboctave_warning_handler) ("Invalid sparse matrix type"); +} + +SparseType::SparseType (const matrix_type t, const octave_idx_type np, + const octave_idx_type *p) : typ (SparseType::Unknown), + nperm (0) +{ + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); + + if (t == SparseType::Permuted_Upper || t == SparseType::Permuted_Lower) + { + typ = t; + nperm = np; + perm = new octave_idx_type [nperm]; + for (octave_idx_type i = 0; i < nperm; i++) + perm[i] = p[i]; + } + else + (*current_liboctave_warning_handler) ("Invalid sparse matrix type"); +} + +SparseType::SparseType (const matrix_type t, const octave_idx_type ku, + const octave_idx_type kl) : typ (SparseType::Unknown), + nperm (0) +{ + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); + + if (t == SparseType::Banded || t == SparseType::Banded_Hermitian) + { + typ = t; + upper_band = ku; + lower_band = kl; + } + else + (*current_liboctave_warning_handler) ("Invalid sparse matrix type"); +} + SparseType::~SparseType (void) { if (nperm != 0) { - delete [] row_perm; - delete [] col_perm; + delete [] perm; } } @@ -470,17 +664,37 @@ if (nperm != 0) { - row_perm = new octave_idx_type [nperm]; - col_perm = new octave_idx_type [nperm]; + perm = new octave_idx_type [nperm]; for (octave_idx_type i = 0; i < nperm; i++) - { - row_perm[i] = a.row_perm[i]; - col_perm[i] = a.col_perm[i]; - } + perm[i] = a.perm[i]; + } } + return *this; +} + +int +SparseType::type (bool quiet) +{ + if (typ != SparseType::Unknown && + sp_bandden == Voctave_sparse_controls.get_key ("bandden")) + { + if (!quiet && + Voctave_sparse_controls.get_key ("spumoni") != 0.) + (*current_liboctave_warning_handler) + ("Using Cached Sparse Matrix Type"); + + return typ; } - return *this; + + if (typ != SparseType::Unknown && + Voctave_sparse_controls.get_key ("spumoni") != 0.) + (*current_liboctave_warning_handler) + ("Invalidating Sparse Matrix Type"); + + typ = SparseType::Unknown; + + return typ; } int @@ -496,11 +710,6 @@ return typ; } - if (Voctave_sparse_controls.get_key ("spumoni") != 0.) - (*current_liboctave_warning_handler) - ("Calculating Sparse Matrix Type"); - - SparseType tmp_typ (a); typ = tmp_typ.typ; sp_bandden = tmp_typ.sp_bandden; @@ -512,13 +721,9 @@ if (nperm != 0) { - row_perm = new octave_idx_type [nperm]; - col_perm = new octave_idx_type [nperm]; + perm = new octave_idx_type [nperm]; for (octave_idx_type i = 0; i < nperm; i++) - { - row_perm[i] = tmp_typ.row_perm[i]; - col_perm[i] = tmp_typ.col_perm[i]; - } + perm[i] = tmp_typ.perm[i]; } return typ; @@ -537,11 +742,6 @@ return typ; } - if (Voctave_sparse_controls.get_key ("spumoni") != 0.) - (*current_liboctave_warning_handler) - ("Calculating Sparse Matrix Type"); - - SparseType tmp_typ (a); typ = tmp_typ.typ; sp_bandden = tmp_typ.sp_bandden; @@ -553,13 +753,9 @@ if (nperm != 0) { - row_perm = new octave_idx_type [nperm]; - col_perm = new octave_idx_type [nperm]; + perm = new octave_idx_type [nperm]; for (octave_idx_type i = 0; i < nperm; i++) - { - row_perm[i] = tmp_typ.row_perm[i]; - col_perm[i] = tmp_typ.col_perm[i]; - } + perm[i] = tmp_typ.perm[i]; } return typ; @@ -593,11 +789,11 @@ ("Permuted Lower Triangular Sparse Matrix"); else if (typ == SparseType::Banded) (*current_liboctave_warning_handler) - ("Banded Sparse Matrix %g-1-%g (Density %g)", lower_band, + ("Banded Sparse Matrix %d-1-%d (Density %f)", lower_band, upper_band, bandden); else if (typ == SparseType::Banded_Hermitian) (*current_liboctave_warning_handler) - ("Banded Hermitian/Symmetric Sparse Matrix %g-1-%g (Density %g)", + ("Banded Hermitian/Symmetric Sparse Matrix %d-1-%d (Density %f)", lower_band, upper_band, bandden); else if (typ == SparseType::Hermitian) (*current_liboctave_warning_handler) @@ -649,16 +845,12 @@ } void -SparseType::mark_as_permuted (const octave_idx_type np, const octave_idx_type *pr, const octave_idx_type *pc) +SparseType::mark_as_permuted (const octave_idx_type np, const octave_idx_type *p) { nperm = np; - row_perm = new octave_idx_type [nperm]; - col_perm = new octave_idx_type [nperm]; + perm = new octave_idx_type [nperm]; for (octave_idx_type i = 0; i < nperm; i++) - { - row_perm[i] = pr[i]; - col_perm[i] = pc[i]; - } + perm[i] = p[i]; if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) typ = SparseType::Permuted_Diagonal; @@ -677,8 +869,7 @@ if (nperm) { nperm = 0; - delete [] row_perm; - delete [] col_perm; + delete [] perm; } if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) @@ -694,23 +885,13 @@ { SparseType retval (*this); if (typ == SparseType::Upper) - retval.typ = Lower; + retval.typ = SparseType::Lower; else if (typ == SparseType::Permuted_Upper) - { - octave_idx_type *tmp = retval.row_perm; - retval.row_perm = retval.col_perm; - retval.col_perm = tmp; - retval.typ = Lower; - } + retval.typ = SparseType::Permuted_Lower; else if (typ == SparseType::Lower) retval.typ = Upper; - else if (typ == SparseType::Permuted_Upper) - { - octave_idx_type *tmp = retval.row_perm; - retval.row_perm = retval.col_perm; - retval.col_perm = tmp; - retval.typ = Upper; - } + else if (typ == SparseType::Permuted_Lower) + retval.typ = Permuted_Upper; else if (typ == SparseType::Banded) { retval.upper_band = lower_band;
--- a/liboctave/SparseType.h Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/SparseType.h Fri Apr 29 13:04:25 2005 +0000 @@ -47,7 +47,7 @@ Rectangular }; - SparseType (void) : typ (Unknown), nperm (0) { } + SparseType (void); SparseType (const SparseType &a); @@ -55,11 +55,19 @@ SparseType (const SparseComplexMatrix &a); + SparseType (const matrix_type t); + + SparseType (const matrix_type t, const octave_idx_type np, + const octave_idx_type *p); + + SparseType (const matrix_type t, const octave_idx_type ku, + const octave_idx_type kl); + ~SparseType (void); SparseType& operator = (const SparseType& a); - int type (void) const { return typ; } + int type (bool quiet = true); int type (const SparseMatrix &a); @@ -100,14 +108,14 @@ void info (void) const; - octave_idx_type * triangular_row_perm (void) const { return row_perm; } + octave_idx_type * triangular_perm (void) const { return perm; } - octave_idx_type * triangular_col_perm (void) const { return col_perm; } - - void invaldate_type (void) { typ = Unknown; } + void invalidate_type (void) { typ = Unknown; } void mark_as_diagonal (void) { typ = Diagonal; } + void mark_as_permuted_diagonal (void) { typ = Permuted_Diagonal; } + void mark_as_upper_triangular (void) { typ = Upper; } void mark_as_lower_triangular (void) { typ = Lower; } @@ -129,7 +137,7 @@ void mark_as_unsymmetric (void); - void mark_as_permuted (const octave_idx_type np, const octave_idx_type *pr, const octave_idx_type *pc); + void mark_as_permuted (const octave_idx_type np, const octave_idx_type *p); void mark_as_unpermuted (void); @@ -145,8 +153,7 @@ octave_idx_type lower_band; bool dense; octave_idx_type nperm; - octave_idx_type *row_perm; - octave_idx_type *col_perm; + octave_idx_type *perm; }; #endif
--- a/liboctave/SparsedbleLU.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/SparsedbleLU.cc Fri Apr 29 13:04:25 2005 +0000 @@ -54,7 +54,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_di_defaults (control); + UMFPACK_DNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -84,18 +84,18 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_di_report_control (control); + UMFPACK_DNAME (report_control) (control); const octave_idx_type *Ap = a.cidx (); const octave_idx_type *Ai = a.ridx (); const double *Ax = a.data (); - umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + 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_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL, + int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, NULL, &Symbolic, control, info); if (status < 0) @@ -103,19 +103,19 @@ (*current_liboctave_error_handler) ("SparseLU::SparseLU symbolic factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_symbolic (&Symbolic) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_di_report_symbolic (Symbolic, control); + UMFPACK_DNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, - control, info) ; - umfpack_di_free_symbolic (&Symbolic) ; + status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, + &Numeric, control, info) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; cond = Info (UMFPACK_RCOND); @@ -124,91 +124,89 @@ (*current_liboctave_error_handler) ("SparseLU::SparseLU numeric factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); } else { - umfpack_di_report_numeric (Numeric, control); + UMFPACK_DNAME (report_numeric) (Numeric, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, - &ignore3, Numeric) ; + octave_idx_type lnz, unz, ignore1, ignore2, ignore3; + status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1, + &ignore2, &ignore3, Numeric) ; if (status < 0) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); } else { - int n_inner = (nr < nc ? nr : nc); + octave_idx_type n_inner = (nr < nc ? nr : nc); if (lnz < 1) - Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, + Lfact = SparseMatrix (n_inner, nr, static_cast<octave_idx_type> (1)); else - Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); + Lfact = SparseMatrix (n_inner, nr, lnz); octave_idx_type *Ltp = Lfact.cidx (); octave_idx_type *Ltj = Lfact.ridx (); double *Ltx = Lfact.data (); if (unz < 1) - Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc, + Ufact = SparseMatrix (n_inner, nc, static_cast<octave_idx_type> (1)); else - Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (unz)); + Ufact = SparseMatrix (n_inner, nc, unz); octave_idx_type *Up = Ufact.cidx (); octave_idx_type *Uj = Ufact.ridx (); double *Ux = Ufact.data (); P.resize (nr); - int *p = P.fortran_vec (); + octave_idx_type *p = P.fortran_vec (); Q.resize (nc); - int *q = Q.fortran_vec (); + octave_idx_type *q = Q.fortran_vec (); - int do_recip; - status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj, - Ux, p, q, (double *) NULL, + octave_idx_type do_recip; + status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx, + Up, Uj, Ux, p, q, (double *) NULL, &do_recip, (double *) NULL, Numeric) ; - umfpack_di_free_numeric (&Numeric) ; + UMFPACK_DNAME (free_numeric) (&Numeric) ; if (status < 0 || do_recip) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); - umfpack_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); } else { Lfact = Lfact.transpose (); - umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (), - Lfact.ridx (), Lfact.data (), - 1, control); - umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (), - Ufact.ridx (), Ufact.data (), - 1, control); - umfpack_di_report_perm (nr, p, control); - umfpack_di_report_perm (nc, q, control); + UMFPACK_DNAME (report_matrix) (nr, n_inner, + Lfact.cidx (), Lfact.ridx (), + Lfact.data (), 1, control); + UMFPACK_DNAME (report_matrix) (n_inner, nc, + Ufact.cidx (), Ufact.ridx (), + Ufact.data (), 1, control); + UMFPACK_DNAME (report_perm) (nr, p, control); + UMFPACK_DNAME (report_perm) (nc, q, control); } - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_info) (control, info); } } } @@ -233,7 +231,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_di_defaults (control); + UMFPACK_DNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -271,13 +269,14 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_di_report_control (control); + UMFPACK_DNAME (report_control) (control); const octave_idx_type *Ap = a.cidx (); const octave_idx_type *Ai = a.ridx (); const double *Ax = a.data (); - umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, + control); void *Symbolic; Matrix Info (1, UMFPACK_INFO); @@ -286,13 +285,13 @@ // Null loop so that qinit is imediately deallocated when not needed do { - OCTAVE_LOCAL_BUFFER (int, qinit, nc); + OCTAVE_LOCAL_BUFFER (octave_idx_type, qinit, nc); - for (int i = 0; i < nc; i++) - qinit [i] = static_cast<int> (Qinit (i)); + for (octave_idx_type i = 0; i < nc; i++) + qinit [i] = static_cast<octave_idx_type> (Qinit (i)); - status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit, - &Symbolic, control, info); + status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, + qinit, &Symbolic, control, info); } while (0); if (status < 0) @@ -300,19 +299,19 @@ (*current_liboctave_error_handler) ("SparseLU::SparseLU symbolic factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_symbolic (&Symbolic) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_di_report_symbolic (Symbolic, control); + UMFPACK_DNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, - control, info) ; - umfpack_di_free_symbolic (&Symbolic) ; + status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, + &Numeric, control, info) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; cond = Info (UMFPACK_RCOND); @@ -321,100 +320,94 @@ (*current_liboctave_error_handler) ("SparseLU::SparseLU numeric factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); } else { - umfpack_di_report_numeric (Numeric, control); + UMFPACK_DNAME (report_numeric) (Numeric, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, - &ignore3, Numeric) ; + octave_idx_type lnz, unz, ignore1, ignore2, ignore3; + status = UMFPACK_DNAME (get_lunz) (&lnz, &unz, &ignore1, &ignore2, + &ignore3, Numeric) ; if (status < 0) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); } else { - int n_inner = (nr < nc ? nr : nc); + octave_idx_type n_inner = (nr < nc ? nr : nc); if (lnz < 1) - Lfact = SparseMatrix - (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (1)); + Lfact = SparseMatrix (n_inner, nr, + static_cast<octave_idx_type> (1)); else - Lfact = SparseMatrix - (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); + Lfact = SparseMatrix (n_inner, nr, lnz); octave_idx_type *Ltp = Lfact.cidx (); octave_idx_type *Ltj = Lfact.ridx (); double *Ltx = Lfact.data (); if (unz < 1) - Ufact = SparseMatrix - (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (1)); + Ufact = SparseMatrix (n_inner, nc, + static_cast<octave_idx_type> (1)); else - Ufact = SparseMatrix - (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (unz)); + Ufact = SparseMatrix (n_inner, nc, unz); octave_idx_type *Up = Ufact.cidx (); octave_idx_type *Uj = Ufact.ridx (); double *Ux = Ufact.data (); P.resize (nr); - int *p = P.fortran_vec (); + octave_idx_type *p = P.fortran_vec (); Q.resize (nc); - int *q = Q.fortran_vec (); + octave_idx_type *q = Q.fortran_vec (); - int do_recip; - status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj, - Ux, p, q, + octave_idx_type do_recip; + status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, + Ltx, Up, Uj, Ux, p, q, (double *) NULL, &do_recip, (double *) NULL, Numeric) ; - umfpack_di_free_numeric (&Numeric) ; + UMFPACK_DNAME (free_numeric) (&Numeric) ; if (status < 0 || do_recip) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); - umfpack_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); } else { Lfact = Lfact.transpose (); - umfpack_di_report_matrix (nr, n_inner, + UMFPACK_DNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), 1, control); - umfpack_di_report_matrix (n_inner, nc, + UMFPACK_DNAME (report_matrix) (n_inner, nc, Ufact.cidx (), Ufact.ridx (), Ufact.data (), 1, control); - umfpack_di_report_perm (nr, p, control); - umfpack_di_report_perm (nc, q, control); + UMFPACK_DNAME (report_perm) (nr, p, control); + UMFPACK_DNAME (report_perm) (nc, q, control); } - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_info) (control, info); } } }
--- a/liboctave/dSparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/dSparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -106,7 +106,7 @@ F77_CHAR_ARG_LEN_DECL); F77_RET_T - F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, Complex*, Complex*, + F77_FUNC (zptsv, ZPTSV) (const octave_idx_type&, const octave_idx_type&, double*, Complex*, Complex*, const octave_idx_type&, octave_idx_type&); F77_RET_T @@ -735,7 +735,7 @@ // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_di_defaults (control); + UMFPACK_DNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -756,38 +756,38 @@ // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - umfpack_di_report_control (control); + UMFPACK_DNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const double *Ax = data (); - umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + 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_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL, - &Symbolic, control, info); + int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, + Ax, NULL, &Symbolic, control, info); if (status < 0) { (*current_liboctave_error_handler) ("SparseMatrix::determinant symbolic factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); - - umfpack_di_free_symbolic (&Symbolic) ; + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_di_report_symbolic (Symbolic, control); + UMFPACK_DNAME (report_symbolic) (Symbolic, control); void *Numeric; - status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, - control, info) ; - umfpack_di_free_symbolic (&Symbolic) ; + status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, + &Numeric, control, info) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; rcond = Info (UMFPACK_RCOND); @@ -796,29 +796,29 @@ (*current_liboctave_error_handler) ("SparseMatrix::determinant numeric factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); - - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_numeric) (&Numeric); } else { - umfpack_di_report_numeric (Numeric, control); + UMFPACK_DNAME (report_numeric) (Numeric, control); double d[2]; - status = umfpack_di_get_determinant (&d[0], &d[1], Numeric, - info); + status = UMFPACK_DNAME (get_determinant) (&d[0], + &d[1], Numeric, info); if (status < 0) { (*current_liboctave_error_handler) ("SparseMatrix::determinant error calculating determinant"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); - - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_numeric) (&Numeric); } else retval = DET (d); @@ -1131,13 +1131,9 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (b.rows (), b.cols ()); + retval.resize (nr, b_cols); + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (double, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); for (octave_idx_type j = 0; j < b_cols; j++) { @@ -1146,28 +1142,29 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; + retval (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -1176,19 +1173,19 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; + 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -1347,12 +1344,12 @@ if (typ == SparseType::Permuted_Upper) { + octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (double, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { @@ -1363,22 +1360,23 @@ for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (work[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -1399,10 +1397,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[rperm[i]]; } retval.xcidx(j+1) = ii; } @@ -1415,19 +1413,19 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; + 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 = q_perm[ridx(i)]; + octave_idx_type idx2 = ridx(i); work[idx2] = work[idx2] - tmp * data(i); } } @@ -1603,75 +1601,71 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (b.rows (), b.cols ()); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + retval.resize (nr, b_nc); + octave_idx_type *perm = mattype.triangular_perm (); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = b(i,j); + cwork[i] = b(i,j); for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (cwork[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (perm[i], j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work2, nr); + OCTAVE_LOCAL_BUFFER (double, work, nr); for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work2[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - double tmp = work2[iidx] / data(cidx(iidx+1)-1); - work2[iidx] = tmp; + 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 = q_perm[ridx(i)]; - work2[idx2] = work2[idx2] - tmp * data(i); + octave_idx_type idx2 = ridx(i); + work[idx2] = work[idx2] - tmp * data(i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { - atmp += fabs(work2[i]); - work2[i] = 0.; + atmp += fabs(work[i]); + work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; @@ -1822,38 +1816,39 @@ if (typ == SparseType::Permuted_Upper) { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); + for (octave_idx_type i = 0; i < nr; i++) + rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; + cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); + cwork[b.ridx(i)] = b.data(i); for (octave_idx_type k = nr-1; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + octave_idx_type kidx = perm[k]; + + if (cwork[k] != 0.) { - if (ridx(cidx(iidx+1)-1) != iidx) + if (ridx(cidx(kidx+1)-1) != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx); i < cidx(iidx+1)-1; i++) + 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 idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + octave_idx_type iidx = ridx(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } @@ -1862,7 +1857,7 @@ // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + if (cwork[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) @@ -1874,45 +1869,45 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (cwork[rperm[i]] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = cwork[rperm[i]]; } retval.xcidx(j+1) = ii; } retval.maybe_compress (); - OCTAVE_LOCAL_BUFFER (double, work2, nr); // Calculation of 1-norm of inv(*this) + OCTAVE_LOCAL_BUFFER (double, work, nr); for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - octave_idx_type iidx = q_perm[k]; - - if (work2[iidx] != 0.) + octave_idx_type iidx = perm[k]; + + if (work[k] != 0.) { - double tmp = work2[iidx] / data(cidx(iidx+1)-1); - work2[iidx] = tmp; + 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 = q_perm[ridx(i)]; - work2[idx2] = work2[idx2] - tmp * data(i); + octave_idx_type idx2 = ridx(i); + work[idx2] = work[idx2] - tmp * data(i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { - atmp += fabs(work2[i]); - work2[i] = 0.; + atmp += fabs(work[i]); + work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; @@ -2084,42 +2079,48 @@ { retval.resize (b.rows (), b.cols ()); OCTAVE_LOCAL_BUFFER (double, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_cols; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = b(i,j); + work[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + double tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) @@ -2128,25 +2129,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -2302,37 +2315,44 @@ if (typ == SparseType::Permuted_Lower) { OCTAVE_LOCAL_BUFFER (double, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + octave_idx_type *perm = mattype.triangular_perm (); for (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); + work[perm[b.ridx(i)]] = b.data(i); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (work[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + double tmp = work[k] / data(mini); + work[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } @@ -2353,10 +2373,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (work[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = work[i]; } retval.xcidx(j+1) = ii; } @@ -2369,25 +2389,37 @@ for (octave_idx_type j = 0; j < nr; j++) { - work[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work[iidx] != 0.) + if (work[k] != 0.) { - double tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -2561,74 +2593,92 @@ if (typ == SparseType::Permuted_Lower) { retval.resize (b.rows (), b.cols ()); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); + octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = b(i,j); + cwork[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (cwork[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = cwork[k] / data(mini); + cwork[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } for (octave_idx_type i = 0; i < nr; i++) - retval (i, j) = work[p_perm[i]]; - + retval (i, j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work2, nr); + OCTAVE_LOCAL_BUFFER (double, work, nr); for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work2[iidx] != 0.) + if (work[k] != 0.) { - double tmp = work2[iidx] / data(cidx(iidx+1)-1); - work2[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work2[idx2] = work2[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { - atmp += fabs(work2[i]); - work2[i] = 0.; + atmp += fabs(work[i]); + work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; @@ -2781,38 +2831,45 @@ if (typ == SparseType::Permuted_Lower) { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - octave_idx_type *p_perm = mattype.triangular_row_perm (); - octave_idx_type *q_perm = mattype.triangular_col_perm (); - - (*current_liboctave_warning_handler) - ("SparseMatrix::solve XXX FIXME XXX permuted triangular code not tested"); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); + octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) - work[i] = 0.; + cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); + cwork[perm[b.ridx(i)]] = b.data(i); for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - if (work[iidx] != 0.) + if (cwork[k] != 0.) { - if (ridx(cidx(iidx)) != iidx) + 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; + } + + if (minr != k) { err = -2; goto triangular_error; } - Complex tmp = work[iidx] / data(cidx(iidx+1)-1); - work[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + Complex tmp = cwork[k] / data(mini); + cwork[k] = tmp; + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work[idx2] = - work[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } @@ -2821,7 +2878,7 @@ // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + if (cwork[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) @@ -2833,10 +2890,10 @@ } for (octave_idx_type i = 0; i < nr; i++) - if (work[p_perm[i]] != 0.) + if (cwork[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[p_perm[i]]; + retval.xdata(ii++) = cwork[i]; } retval.xcidx(j+1) = ii; } @@ -2844,34 +2901,46 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work2, nr); + OCTAVE_LOCAL_BUFFER (double, work, nr); for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[q_perm[j]] = 1.; + work[j] = 1.; for (octave_idx_type k = 0; k < nr; k++) { - octave_idx_type iidx = q_perm[k]; - - if (work2[iidx] != 0.) + if (work[k] != 0.) { - double tmp = work2[iidx] / data(cidx(iidx+1)-1); - work2[iidx] = tmp; - for (octave_idx_type i = cidx(iidx)+1; i < cidx(iidx+1); i++) + 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++) { - octave_idx_type idx2 = q_perm[ridx(i)]; - work2[idx2] = work2[idx2] - tmp * data(i); + if (i == mini) + continue; + + octave_idx_type iidx = perm[ridx(i)]; + work[iidx] = work[iidx] - tmp * data(i); } } } + double atmp = 0; - for (octave_idx_type i = 0; i < j+1; i++) + for (octave_idx_type i = j; i < nr; i++) { - atmp += fabs(work2[i]); - work2[i] = 0.; + atmp += fabs(work[i]); + work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; @@ -3115,7 +3184,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3211,7 +3280,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3319,10 +3388,9 @@ volatile int typ = mattype.type (); mattype.info (); - // Note can't treat symmetric case as there is no dpttrf function if (typ == SparseType::Tridiagonal_Hermitian) { - OCTAVE_LOCAL_BUFFER (Complex, D, nr); + OCTAVE_LOCAL_BUFFER (double, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) @@ -3416,7 +3484,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3516,7 +3584,7 @@ else if (ridx(i) == j + 1) DL[j] = data(i); else if (ridx(i) == j - 1) - DU[j] = data(i); + DU[j-1] = data(i); } } @@ -3564,7 +3632,6 @@ Bz[i] = std::imag (c); } - F77_XFCN (dgttrs, DGTTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 1, DL, D, DU, DU2, pipvt, @@ -4682,7 +4749,7 @@ // Setup the control parameters Control = Matrix (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); - umfpack_di_defaults (control); + UMFPACK_DNAME (defaults) (control); double tmp = Voctave_sparse_controls.get_key ("spumoni"); if (!xisnan (tmp)) @@ -4699,7 +4766,7 @@ if (!xisnan (tmp)) Control (UMFPACK_FIXQ) = tmp; - umfpack_di_report_control (control); + UMFPACK_DNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); @@ -4707,12 +4774,12 @@ octave_idx_type nr = rows (); octave_idx_type nc = cols (); - umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + UMFPACK_DNAME (report_matrix) (nr, nc, Ap, Ai, Ax, 1, control); void *Symbolic; Info = Matrix (1, UMFPACK_INFO); double *info = Info.fortran_vec (); - int status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, NULL, + int status = UMFPACK_DNAME (qsymbolic) (nr, nc, Ap, Ai, Ax, NULL, &Symbolic, control, info); if (status < 0) @@ -4721,18 +4788,18 @@ ("SparseMatrix::solve symbolic factorization failed"); err = -1; - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); - - umfpack_di_free_symbolic (&Symbolic) ; + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; } else { - umfpack_di_report_symbolic (Symbolic, control); - - status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, - control, info) ; - umfpack_di_free_symbolic (&Symbolic) ; + UMFPACK_DNAME (report_symbolic) (Symbolic, control); + + status = UMFPACK_DNAME (numeric) (Ap, Ai, Ax, Symbolic, + &Numeric, control, info) ; + UMFPACK_DNAME (free_symbolic) (&Symbolic) ; #ifdef HAVE_LSSOLVE rcond = Info (UMFPACK_RCOND); @@ -4741,7 +4808,7 @@ if (status == UMFPACK_WARNING_singular_matrix || rcond_plus_one == 1.0 || xisnan (rcond)) { - umfpack_di_report_numeric (Numeric, control); + UMFPACK_DNAME (report_numeric) (Numeric, control); err = -2; @@ -4760,19 +4827,19 @@ (*current_liboctave_error_handler) ("SparseMatrix::solve numeric factorization failed"); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_status) (control, status); + UMFPACK_DNAME (report_info) (control, info); err = -1; } else { - umfpack_di_report_numeric (Numeric, control); + UMFPACK_DNAME (report_numeric) (Numeric, control); } } if (err != 0) - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -4836,15 +4903,15 @@ for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { - status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, - &result[iidx], &Bx[iidx], + 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_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); err = -1; @@ -4871,9 +4938,9 @@ } #endif - umfpack_di_report_info (control, info); + UMFPACK_DNAME (report_info) (control, info); - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -4951,15 +5018,15 @@ for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b.elem (i, j); - status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, Xx, - Bx, Numeric, control, + 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_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); err = -1; @@ -5007,9 +5074,9 @@ } #endif - umfpack_di_report_info (control, info); - - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -5088,11 +5155,11 @@ Bz[i] = std::imag (c); } - status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, - Xx, Bx, Numeric, control, + status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap, + Ai, Ax, Xx, Bx, Numeric, control, info); - int status2 = umfpack_di_solve (UMFPACK_A, Ap, Ai, - Ax, Xz, Bz, Numeric, + int status2 = UMFPACK_DNAME (solve) (UMFPACK_A, + Ap, Ai, Ax, Xz, Bz, Numeric, control, info) ; if (status < 0 || status2 < 0) @@ -5100,7 +5167,7 @@ (*current_liboctave_error_handler) ("SparseMatrix::solve solve failed"); - umfpack_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); err = -1; @@ -5130,9 +5197,9 @@ } #endif - umfpack_di_report_info (control, info); - - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -5217,11 +5284,11 @@ Bz[i] = std::imag (c); } - status = umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, Xx, - Bx, Numeric, control, + status = UMFPACK_DNAME (solve) (UMFPACK_A, Ap, + Ai, Ax, Xx, Bx, Numeric, control, info); - int status2 = umfpack_di_solve (UMFPACK_A, Ap, Ai, - Ax, Xz, Bz, Numeric, + int status2 = UMFPACK_DNAME (solve) (UMFPACK_A, + Ap, Ai, Ax, Xz, Bz, Numeric, control, info) ; if (status < 0 || status2 < 0) @@ -5229,7 +5296,7 @@ (*current_liboctave_error_handler) ("SparseMatrix::solve solve failed"); - umfpack_di_report_status (control, status); + UMFPACK_DNAME (report_status) (control, status); err = -1; @@ -5277,9 +5344,9 @@ } #endif - umfpack_di_report_info (control, info); - - umfpack_di_free_numeric (&Numeric); + UMFPACK_DNAME (report_info) (control, info); + + UMFPACK_DNAME (free_numeric) (&Numeric); } #else (*current_liboctave_error_handler) ("UMFPACK not installed"); @@ -5319,7 +5386,7 @@ double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5373,7 +5440,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5427,7 +5494,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this); @@ -5481,7 +5548,7 @@ octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { - int typ = mattype.type (); + int typ = mattype.type (false); if (typ == SparseType::Unknown) typ = mattype.type (*this);
--- a/liboctave/dSparse.h Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/dSparse.h Fri Apr 29 13:04:25 2005 +0000 @@ -404,6 +404,12 @@ SPARSE_FORWARD_DEFS (MSparse, SparseMatrix, Matrix, double) +#ifdef UMFPACK_LONG_IDX +#define UMFPACK_DNAME(name) umfpack_dl_ ## name +#else +#define UMFPACK_DNAME(name) umfpack_di_ ## name +#endif + #endif /*
--- a/liboctave/sparse-base-lu.h Fri Apr 29 05:20:01 2005 +0000 +++ b/liboctave/sparse-base-lu.h Fri Apr 29 13:04:25 2005 +0000 @@ -60,9 +60,9 @@ p_type Pr (void) const; - MArray<int> row_perm (void) const { return P; } + const octave_idx_type * row_perm (void) const { return P.fortran_vec (); } - MArray<int> col_perm (void) const { return Q; } + const octave_idx_type * col_perm (void) const { return Q.fortran_vec (); } double rcond (void) const { return cond; } @@ -73,8 +73,8 @@ double cond; - MArray<int> P; - MArray<int> Q; + MArray<octave_idx_type> P; + MArray<octave_idx_type> Q; }; #endif
--- a/src/ChangeLog Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ChangeLog Fri Apr 29 13:04:25 2005 +0000 @@ -1,3 +1,44 @@ +2005-04-29 David Bateman <dbateman@free.fr> + + * Makefile.in: Add matrix_type.cc and spkron.cc to DLD_XSRC + + * ls.mat.cc (read_mat5_binary_element): Allow for endian change for + compressed data elements. + + * ov-base-sparse.cc (assign): Invalidate matrix type. + + * ov-base-sparse.cc (SparseType sparse_type (void), + SparseType sparse_type (const SparseType&): Functions to read and set + sparse matrix type. + + * ov-bool-sparse.cc (load_binary): Remove third argument + (load_hdf5): Cast hsize_t comparisions with int to avoid warning + * ov-cx-sparse.cc (load_hdf5): ditto + * ov-re-sparse.cc (load_hdf5): ditto + + * ov-re-sparse.cc (convert_to_str_internal): Add third argument for string + type. + * ov-re-sparse.h (convert_to_str_internal): Adject declaration. + + * sparse-xdiv.cc (xdiv, xleftdiv): Pass SparseType as third argument, use + it and return it to allow caching of type. + * sparse-xdiv.h (xdiv, xleftdiv): Change declarations for third argument + of type SparseType. + + * DLD-FUNCTIONS/luinc.cc (Fluinc): use type_name and not class_name to + test for real/complex sparse matrices. Set matrix type + + * DLD-FUNCTIONS/splu.cc (Fsplu): Set matrix type. + + * OPERATORS/op-cm-scm.cc, OPERATORS/op-cm-sm.cc, OPERATORS/op-cs-scm.cc, + OPERATORS/op-cs-sm.cc, OPERATORS/op-m-scm.cc, OPERATORS/op-m-sm.cc, + OPERATORS/op-s-scm.cc, OPERATORS/op-s-sm.cc, OPERATORS/op-scm-cm.cc, + OPERATORS/op-scm-cs.cc, OPERATORS/op-scm-m.cc, OPERATORS/op-scm-s.cc, + OPERATORS/op-scm-scm.cc, OPERATORS/op-scm-sm.cc, OPERATORS/op-sm-cm.cc, + OPERATORS/op-sm-cs.cc, OPERATORS/op-sm-m.cc, OPERATORS/op-sm-s.cc, + OPERATORS/op-sm-scm.cc, OPERATORS/op-sm-sm.cc (div, ldiv): Pass and + recache SparseType wirh xdiv/xleftdiv. + 2005-04-29 John W. Eaton <jwe@octave.org> * file-io.cc (maybe_warn_interface_change): Delete function.
--- a/src/DLD-FUNCTIONS/luinc.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/DLD-FUNCTIONS/luinc.cc Fri Apr 29 13:04:25 2005 +0000 @@ -30,6 +30,7 @@ #include "utils.h" #include "oct-map.h" +#include "SparseType.h" #include "SparseCmplxLU.h" #include "SparsedbleLU.h" #include "ov-re-sparse.h" @@ -146,9 +147,10 @@ if (!error_state) { - if (args(0).class_name () == "sparse") + if (args(0).type_name () == "sparse matrix") { SparseMatrix sm = args(0).sparse_matrix_value (); + octave_idx_type sm_nr = sm.rows (); octave_idx_type sm_nc = sm.cols (); ColumnVector Qinit (sm_nc); @@ -168,8 +170,11 @@ SparseMatrix P = fact.Pr (); SparseMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); - retval(0) = L; + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (L, SparseType + (SparseType::Permuted_Lower, + sm_nr, fact.row_perm ())); } break; @@ -179,8 +184,10 @@ milu, udiag); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; @@ -192,17 +199,20 @@ retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; } } } - else if (args(0).class_name () == "sparse complex") + else if (args(0).type_name () == "sparse complex matrix") { SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); + octave_idx_type sm_nr = sm.rows (); octave_idx_type sm_nc = sm.cols (); ColumnVector Qinit (sm_nc); @@ -222,8 +232,11 @@ SparseMatrix P = fact.Pr (); SparseComplexMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); - retval(0) = L; + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (L, SparseType + (SparseType::Permuted_Lower, + sm_nr, fact.row_perm ())); } break; @@ -233,8 +246,10 @@ droptol, milu, udiag); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; @@ -246,8 +261,10 @@ retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; }
--- a/src/DLD-FUNCTIONS/splu.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/DLD-FUNCTIONS/splu.cc Fri Apr 29 13:04:25 2005 +0000 @@ -239,11 +239,14 @@ SparseMatrix P = fact.Pr (); SparseMatrix L = P.transpose () * fact.L (); if (have_Qinit) - retval(1) = fact.U () * fact.Pc ().transpose (); + retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), + SparseType (SparseType::Permuted_Upper, nc, fact.col_perm ())); else - retval(1) = fact.U (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); - retval(0) = L; + retval(0) = octave_value (L, + SparseType (SparseType::Permuted_Lower, nr, fact.row_perm ())); } break; @@ -253,10 +256,14 @@ retval(2) = fact.Pr (); if (have_Qinit) - retval(1) = fact.U () * fact.Pc ().transpose (); + retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), + SparseType (SparseType::Permuted_Upper, nc, fact.col_perm ())); else - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; @@ -269,8 +276,10 @@ retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } else { @@ -278,8 +287,10 @@ retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } } break; @@ -312,10 +323,14 @@ SparseComplexMatrix L = P.transpose () * fact.L (); if (have_Qinit) - retval(1) = fact.U () * fact.Pc ().transpose (); + retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), + SparseType (SparseType::Permuted_Upper, nc, fact.col_perm ())); else - retval(1) = fact.U (); - retval(0) = L; + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + + retval(0) = octave_value (L, + SparseType (SparseType::Permuted_Lower, nr, fact.row_perm ())); } break; @@ -325,10 +340,14 @@ retval(2) = fact.Pr (); if (have_Qinit) - retval(1) = fact.U () * fact.Pc ().transpose (); + retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), + SparseType (SparseType::Permuted_Upper, nc, fact.col_perm ())); else - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } break; @@ -338,11 +357,13 @@ if (have_Qinit) { SparseComplexLU fact (m, Qinit, thres, false); - + retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } else { @@ -350,8 +371,10 @@ retval(3) = fact.Pc (); retval(2) = fact.Pr (); - retval(1) = fact.U (); - retval(0) = fact.L (); + retval(1) = octave_value (fact.U (), + SparseType (SparseType::Upper)); + retval(0) = octave_value (fact.L (), + SparseType (SparseType::Lower)); } } break;
--- a/src/Makefile.in Fri Apr 29 05:20:01 2005 +0000 +++ b/src/Makefile.in Fri Apr 29 13:04:25 2005 +0000 @@ -47,10 +47,10 @@ eig.cc expm.cc fft.cc fft2.cc fftn.cc fftw_wisdom.cc \ filter.cc find.cc fsolve.cc gammainc.cc gcd.cc getgrent.cc \ getpwent.cc getrusage.cc givens.cc hess.cc inv.cc kron.cc \ - lpsolve.cc lsode.cc lu.cc luinc.cc minmax.cc pinv.cc qr.cc \ - quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \ - splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l \ - __glpk__.cc __qp__.cc + lpsolve.cc lsode.cc lu.cc luinc.cc matrix_type.cc minmax.cc \ + pinv.cc qr.cc quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc \ + spdet.cc spkron.cc splu.cc spparms.cc sqrtm.cc svd.cc syl.cc \ + time.cc gplot.l __glpk__.cc __qp__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- a/src/OPERATORS/op-cm-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-cm-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -55,11 +55,15 @@ DEFBINOP (div, complex_matrix, sparse_complex_matrix) { - CAST_BINOP_ARGS (const octave_complex_matrix&, - const octave_sparse_complex_matrix&); + CAST_BINOP_ARGS (const octave_complex_matrix&, octave_sparse_complex_matrix&); - return xdiv (v1.complex_matrix_value (), - v2.sparse_complex_matrix_value ()); + SparseType typ = v2.sparse_type (); + + ComplexMatrix ret = xdiv (v1.complex_matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; } DEFBINOPX (pow, complex_matrix, sparse_complex_matrix)
--- a/src/OPERATORS/op-cm-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-cm-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -55,10 +55,15 @@ DEFBINOP (div, complex_matrix, sparse_matrix) { - CAST_BINOP_ARGS (const octave_complex_matrix&, - const octave_sparse_matrix&); + CAST_BINOP_ARGS (const octave_complex_matrix&, octave_sparse_matrix&); - return xdiv (v1.complex_matrix_value (), v2.sparse_matrix_value ()); + SparseType typ = v2.sparse_type (); + + ComplexMatrix ret = xdiv (v1.complex_matrix_value (), + v2.sparse_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; } DEFBINOPX (pow, complex_matrix, sparse_matrix)
--- a/src/OPERATORS/op-cs-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-cs-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -45,13 +45,15 @@ DEFBINOP (div, complex, sparse_complex_matrix) { - CAST_BINOP_ARGS (const octave_complex&, - const octave_sparse_complex_matrix&); + CAST_BINOP_ARGS (const octave_complex&, octave_sparse_complex_matrix&); + SparseType typ = v2.sparse_type (); ComplexMatrix m1 = ComplexMatrix (1, 1, v1.complex_value ()); SparseComplexMatrix m2 = v2.sparse_complex_matrix_value (); + ComplexMatrix ret = xdiv (m1, m2, typ); + v2.sparse_type (typ); - return xdiv (m1, m2); + return ret; } DEFBINOP (pow, complex, sparse_complex_matrix)
--- a/src/OPERATORS/op-cs-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-cs-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -47,12 +47,15 @@ DEFBINOP (div, complex, sparse_matrix) { - CAST_BINOP_ARGS (const octave_complex&, const octave_sparse_matrix&); + CAST_BINOP_ARGS (const octave_complex&, octave_sparse_matrix&); + SparseType typ = v2.sparse_type (); ComplexMatrix m1 = ComplexMatrix (1, 1, v1.complex_value ()); SparseMatrix m2 = v2.sparse_matrix_value (); + ComplexMatrix ret = xdiv (m1, m2, typ); + v2.sparse_type (typ); - return xdiv (m1, m2); + return ret; } DEFBINOP (pow, complex, sparse_matrix)
--- a/src/OPERATORS/op-m-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-m-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -56,10 +56,15 @@ DEFBINOP (div, matrix, sparse_complex_matrix) { - CAST_BINOP_ARGS (const octave_matrix&, - const octave_sparse_complex_matrix&); - - return xdiv (v1.matrix_value (), v2.sparse_complex_matrix_value ()); + CAST_BINOP_ARGS (const octave_matrix&, octave_sparse_complex_matrix&); + + SparseType typ = v2.sparse_type (); + + ComplexMatrix ret = xdiv (v1.matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; } DEFBINOPX (pow, matrix, sparse_complex_matrix)
--- a/src/OPERATORS/op-m-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-m-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -54,9 +54,13 @@ DEFBINOP (div, matrix, sparse_matrix) { - CAST_BINOP_ARGS (const octave_matrix&, const octave_sparse_matrix&); - - return xdiv (v1.matrix_value (), v2.sparse_matrix_value ()); + CAST_BINOP_ARGS (const octave_matrix&, octave_sparse_matrix&); + SparseType typ = v2.sparse_type (); + + Matrix ret = xdiv (v1.matrix_value (), v2.sparse_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; } DEFBINOPX (pow, matrix, sparse_matrix)
--- a/src/OPERATORS/op-s-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-s-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -48,13 +48,15 @@ DEFBINOP (div, scalar, sparse_complex_matrix) { - CAST_BINOP_ARGS (const octave_scalar&, - const octave_sparse_complex_matrix&); + CAST_BINOP_ARGS (const octave_scalar&, octave_sparse_complex_matrix&); + SparseType typ = v2.sparse_type (); Matrix m1 = Matrix (1, 1, v1.scalar_value ()); SparseComplexMatrix m2 = v2.sparse_complex_matrix_value (); + ComplexMatrix ret = xdiv (m1, m2, typ); + v2.sparse_type (typ); - return xdiv (m1, m2); + return ret; } DEFBINOP (pow, scalar, sparse_complex_matrix)
--- a/src/OPERATORS/op-s-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-s-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -44,12 +44,15 @@ DEFBINOP (div, scalar, sparse_matrix) { - CAST_BINOP_ARGS (const octave_scalar&, const octave_sparse_matrix&); + CAST_BINOP_ARGS (const octave_scalar&, octave_sparse_matrix&); + SparseType typ = v2.sparse_type (); Matrix m1 = Matrix (1, 1, v1.double_value ()); SparseMatrix m2 = v2.sparse_matrix_value (); + Matrix ret = xdiv (m1, m2, typ); + v2.sparse_type (typ); - return xdiv (m1, m2); + return ret; } DEFBINOP (pow, scalar, sparse_matrix)
--- a/src/OPERATORS/op-scm-cm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-cm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -69,11 +69,14 @@ DEFBINOP (ldiv, sparse_complex_matrix, complex_matrix) { - CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, - const octave_complex_matrix&); - - return xleftdiv (v1.sparse_complex_matrix_value (), - v2.complex_matrix_value ()); + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, const octave_complex_matrix&); + SparseType typ = v1.sparse_type (); + + ComplexMatrix ret = xleftdiv (v1.sparse_complex_matrix_value (), + v2.complex_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; } DEFBINOP_FN (lt, sparse_complex_matrix, complex_matrix, mx_el_lt)
--- a/src/OPERATORS/op-scm-cs.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-cs.cc Fri Apr 29 13:04:25 2005 +0000 @@ -72,13 +72,15 @@ DEFBINOP (ldiv, sparse_complex_matrix, complex) { - CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, - const octave_complex&); + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, const octave_complex&); + SparseType typ = v1.sparse_type (); SparseComplexMatrix m1 = v1.sparse_complex_matrix_value (); ComplexMatrix m2 = ComplexMatrix (1, 1, v2.complex_value ()); + ComplexMatrix ret = xleftdiv (m1, m2, typ); + v1.sparse_type (typ); - return xleftdiv (m1, m2); + return ret; } DEFBINOP_FN (lt, sparse_complex_matrix, complex, mx_el_lt)
--- a/src/OPERATORS/op-scm-m.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-m.cc Fri Apr 29 13:04:25 2005 +0000 @@ -57,7 +57,7 @@ DEFBINOP (div, sparse_complex_matrix, matrix) { CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, - const octave_matrix&); + const octave_matrix&); return xdiv (v1.complex_matrix_value (), v2.matrix_value ()); } @@ -70,10 +70,15 @@ DEFBINOP (ldiv, sparse_complex_matrix, matrix) { - CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, - const octave_matrix&); + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, const octave_matrix&); - return xleftdiv (v1.sparse_complex_matrix_value (), v2.matrix_value ()); + SparseType typ = v1.sparse_type (); + + ComplexMatrix ret = xleftdiv (v1.sparse_complex_matrix_value (), + v2.matrix_value (), typ); + + v1.sparse_type (typ); + return ret; } DEFBINOP_FN (lt, sparse_complex_matrix, matrix, mx_el_lt)
--- a/src/OPERATORS/op-scm-s.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-s.cc Fri Apr 29 13:04:25 2005 +0000 @@ -80,13 +80,15 @@ DEFBINOP (ldiv, sparse_complex_matrix, scalar) { - CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, - const octave_scalar&); + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, const octave_scalar&); + SparseType typ = v1.sparse_type (); SparseComplexMatrix m1 = v1.sparse_complex_matrix_value (); Matrix m2 = Matrix (1, 1, v2.scalar_value ()); + ComplexMatrix ret = xleftdiv (m1, m2, typ); + v1.sparse_type (typ); - return xleftdiv (m1, m2); + return ret; } DEFBINOP_FN (lt, sparse_complex_matrix, scalar, mx_el_lt)
--- a/src/OPERATORS/op-scm-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -57,15 +57,17 @@ DEFUNOP (transpose, sparse_complex_matrix) { CAST_UNOP_ARG (const octave_sparse_complex_matrix&); - - return octave_value (v.sparse_complex_matrix_value().transpose ()); + return octave_value + (v.sparse_complex_matrix_value().transpose (), + v.sparse_type ().transpose ()); } DEFUNOP (hermitian, sparse_complex_matrix) { CAST_UNOP_ARG (const octave_sparse_complex_matrix&); - - return octave_value (v.sparse_complex_matrix_value().hermitian ()); + return octave_value + (v.sparse_complex_matrix_value().hermitian (), + v.sparse_type ().transpose ()); } #if 0 @@ -91,7 +93,17 @@ DEFBINOP_OP (mul, sparse_complex_matrix, sparse_complex_matrix, *) -DEFBINOP_FN (div, sparse_complex_matrix, sparse_complex_matrix, xdiv) +DEFBINOP (div, sparse_complex_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, + octave_sparse_complex_matrix&); + SparseType typ = v2.sparse_type (); + SparseComplexMatrix ret = xdiv (v1.sparse_complex_matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; +} DEFBINOPX (pow, sparse_complex_matrix, sparse_complex_matrix) { @@ -99,7 +111,18 @@ return octave_value (); } -DEFBINOP_FN (ldiv, sparse_complex_matrix, sparse_complex_matrix, xleftdiv) +DEFBINOP (ldiv, sparse_complex_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, + const octave_sparse_complex_matrix&); + SparseType typ = v1.sparse_type (); + + SparseComplexMatrix ret = xleftdiv (v1.sparse_complex_matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; +} DEFBINOP_FN (lt, sparse_complex_matrix, sparse_complex_matrix, mx_el_lt) DEFBINOP_FN (le, sparse_complex_matrix, sparse_complex_matrix, mx_el_le)
--- a/src/OPERATORS/op-scm-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-scm-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -44,7 +44,16 @@ DEFBINOP_OP (mul, sparse_complex_matrix, sparse_matrix, *) -DEFBINOP_FN (div, sparse_complex_matrix, sparse_matrix, xdiv) +DEFBINOP (div, sparse_complex_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, octave_sparse_matrix&); + SparseType typ = v2.sparse_type (); + SparseComplexMatrix ret = xdiv (v1.sparse_complex_matrix_value (), + v2.sparse_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; +} DEFBINOPX (pow, sparse_complex_matrix, sparse_matrix) { @@ -52,7 +61,17 @@ return octave_value (); } -DEFBINOP_FN (ldiv, sparse_complex_matrix, sparse_matrix, xleftdiv) +DEFBINOP (ldiv, sparse_complex_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (octave_sparse_complex_matrix&, const octave_sparse_matrix&); + SparseType typ = v1.sparse_type (); + + SparseComplexMatrix ret = xleftdiv (v1.sparse_complex_matrix_value (), + v2.sparse_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; +} DEFBINOP_FN (lt, sparse_complex_matrix, sparse_matrix, mx_el_lt) DEFBINOP_FN (le, sparse_complex_matrix, sparse_matrix, mx_el_le)
--- a/src/OPERATORS/op-sm-cm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-cm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -69,10 +69,14 @@ DEFBINOP (ldiv, sparse_matrix, complex_matrix) { - CAST_BINOP_ARGS (const octave_sparse_matrix&, - const octave_complex_matrix&); - - return xleftdiv (v1.sparse_matrix_value (), v2.complex_matrix_value ()); + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_complex_matrix&); + SparseType typ = v1.sparse_type (); + + ComplexMatrix ret = xleftdiv (v1.sparse_matrix_value (), + v2.complex_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; } DEFBINOP_FN (lt, sparse_matrix, complex_matrix, mx_el_lt)
--- a/src/OPERATORS/op-sm-cs.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-cs.cc Fri Apr 29 13:04:25 2005 +0000 @@ -72,12 +72,15 @@ DEFBINOP (ldiv, sparse_matrix, complex) { - CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_complex&); + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_complex&); + SparseType typ = v1.sparse_type (); SparseMatrix m1 = v1.sparse_matrix_value (); ComplexMatrix m2 = ComplexMatrix (1, 1, v2.complex_value ()); + ComplexMatrix ret = xleftdiv (m1, m2, typ); + v1.sparse_type (typ); - return xleftdiv (m1, m2); + return ret; } DEFBINOP_FN (lt, sparse_matrix, complex, mx_el_lt)
--- a/src/OPERATORS/op-sm-m.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-m.cc Fri Apr 29 13:04:25 2005 +0000 @@ -67,12 +67,16 @@ DEFBINOP (ldiv, sparse_matrix, matrix) { - CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_matrix&); - - return xleftdiv (v1.sparse_matrix_value (), v2.matrix_value ()); + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_matrix&); + SparseType typ = v1.sparse_type (); + + Matrix ret = xleftdiv (v1.sparse_matrix_value (), + v2.matrix_value (), typ); + + v1.sparse_type (typ); + return ret; } - DEFBINOP_FN (lt, sparse_matrix, matrix, mx_el_lt) DEFBINOP_FN (le, sparse_matrix, matrix, mx_el_le) DEFBINOP_FN (eq, sparse_matrix, matrix, mx_el_eq)
--- a/src/OPERATORS/op-sm-s.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-s.cc Fri Apr 29 13:04:25 2005 +0000 @@ -74,12 +74,15 @@ DEFBINOP (ldiv, sparse_matrix, scalar) { - CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_scalar&); + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_scalar&); + SparseType typ = v1.sparse_type (); SparseMatrix m1 = v1.sparse_matrix_value (); Matrix m2 = Matrix (1, 1, v2.scalar_value ()); + Matrix ret = xleftdiv (m1, m2, typ); + v1.sparse_type (typ); - return xleftdiv (m1, m2); + return ret; } DEFBINOP_FN (lt, sparse_matrix, scalar, mx_el_lt)
--- a/src/OPERATORS/op-sm-scm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-scm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -44,7 +44,16 @@ DEFBINOP_OP (mul, sparse_matrix, sparse_complex_matrix, *) -DEFBINOP_FN (div, sparse_matrix, sparse_complex_matrix, xdiv) +DEFBINOP (div, sparse_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, octave_sparse_complex_matrix&); + SparseType typ = v2.sparse_type (); + SparseComplexMatrix ret = xdiv (v1.sparse_matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; +} DEFBINOPX (pow, sparse_matrix, sparse_complex_matrix) { @@ -52,7 +61,17 @@ return octave_value (); } -DEFBINOP_FN (ldiv, sparse_matrix, sparse_complex_matrix, xleftdiv) +DEFBINOP (ldiv, sparse_matrix, sparse_complex_matrix) +{ + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_sparse_complex_matrix&); + SparseType typ = v1.sparse_type (); + + SparseComplexMatrix ret = xleftdiv (v1.sparse_matrix_value (), + v2.sparse_complex_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; +} DEFBINOP_FN (lt, sparse_matrix, sparse_complex_matrix, mx_el_lt) DEFBINOP_FN (le, sparse_matrix, sparse_complex_matrix, mx_el_le)
--- a/src/OPERATORS/op-sm-sm.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/OPERATORS/op-sm-sm.cc Fri Apr 29 13:04:25 2005 +0000 @@ -44,7 +44,8 @@ DEFUNOP (transpose, sparse_matrix) { CAST_UNOP_ARG (const octave_sparse_matrix&); - return octave_value (v.sparse_matrix_value().transpose ()); + return octave_value (v.sparse_matrix_value().transpose (), + v.sparse_type ().transpose ()); } // sparse matrix by sparse matrix ops. @@ -53,7 +54,16 @@ DEFBINOP_OP (sub, sparse_matrix, sparse_matrix, -) DEFBINOP_OP (mul, sparse_matrix, sparse_matrix, *) -DEFBINOP_FN (div, sparse_matrix, sparse_matrix, xdiv) +DEFBINOP (div, sparse_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (const octave_sparse_matrix&, octave_sparse_matrix&); + SparseType typ = v2.sparse_type (); + SparseMatrix ret = xdiv (v1.sparse_matrix_value (), + v2.sparse_matrix_value (), typ); + + v2.sparse_type (typ); + return ret; +} DEFBINOPX (pow, sparse_matrix, sparse_matrix) { @@ -61,7 +71,17 @@ return octave_value (); } -DEFBINOP_FN (ldiv, sparse_matrix, sparse_matrix, xleftdiv) +DEFBINOP (ldiv, sparse_matrix, sparse_matrix) +{ + CAST_BINOP_ARGS (octave_sparse_matrix&, const octave_sparse_matrix&); + SparseType typ = v1.sparse_type (); + + SparseMatrix ret = xleftdiv (v1.sparse_matrix_value (), + v2.sparse_matrix_value (), typ); + + v1.sparse_type (typ); + return ret; +} DEFBINOP_FN (lt, sparse_matrix, sparse_matrix, mx_el_lt) DEFBINOP_FN (le, sparse_matrix, sparse_matrix, mx_el_le)
--- a/src/load-save.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/load-save.cc Fri Apr 29 13:04:25 2005 +0000 @@ -1151,6 +1151,10 @@ { use_zlib = true; } + else if (argv[i] == "-nozip" || argv[i] == "-nz") + { + use_zlib = false; + } #endif else break; @@ -1507,7 +1511,13 @@ @itemx -z\n\ Use the gzip algorithm to compress the file. This works equally on files that\n\ are compressed with gzip outside of octave, and gzip can equally be used to\n\ -convert the files for backward compatibility" +convert the files for backward compatibility." + +HAVE_ZLIB_HELP_STRING + +"@item -nozip\n\ +@itemx -nz\n\ +Disable the use of the file compression." HAVE_ZLIB_HELP_STRING
--- a/src/ls-mat5.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ls-mat5.cc Fri Apr 29 13:04:25 2005 +0000 @@ -441,6 +441,9 @@ X_CAST (Bytef *, inbuf), element_length) != Z_MEM_ERROR) { // Why should I have to initialize outbuf as I'll just overwrite!! + if (swap) + swap_bytes<4> (tmp, 2); + destLen = tmp[1] + 8; std::string outbuf (destLen, ' ');
--- a/src/ov-base-sparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-base-sparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -188,6 +188,9 @@ matrix.set_index (idx(i).index_vector ()); ::assign (matrix, rhs); + + // Invalidate matrix type. + typ.invalidate_type (); }
--- a/src/ov-base-sparse.h Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-base-sparse.h Fri Apr 29 13:04:25 2005 +0000 @@ -117,6 +117,10 @@ octave_value all (int dim = 0) const { return matrix.all (dim); } octave_value any (int dim = 0) const { return matrix.any (dim); } + SparseType sparse_type (void) const { return typ; } + SparseType sparse_type (const SparseType& _typ) + { SparseType ret = typ; typ = _typ; return ret; } + bool is_matrix_type (void) const { return true; } bool is_numeric_type (void) const { return true; }
--- a/src/ov-bool-sparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-bool-sparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -249,7 +249,7 @@ bool octave_sparse_bool_matrix::load_binary (std::istream& is, bool swap, - oct_mach_info::float_format fmt) + oct_mach_info::float_format /* fmt */) { FOUR_BYTE_INT nz, nc, nr, tmp; if (! is.read (X_CAST (char *, &tmp), 4)) @@ -596,7 +596,8 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nc + 1 || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nc + 1 || + static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -631,7 +632,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -666,7 +667,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid);
--- a/src/ov-cx-sparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-cx-sparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -650,7 +650,8 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nc + 1 || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nc + 1 || + static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -685,7 +686,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -732,7 +733,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid);
--- a/src/ov-re-sparse.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-re-sparse.cc Fri Apr 29 13:04:25 2005 +0000 @@ -184,7 +184,7 @@ } octave_value -octave_sparse_matrix::convert_to_str_internal (bool, bool) const +octave_sparse_matrix::convert_to_str_internal (bool, bool, char type) const { octave_value retval; dim_vector dv = dims (); @@ -193,7 +193,7 @@ if (nel == 0) { char s = '\0'; - retval = octave_value (&s); + retval = octave_value (&s, type); } else { @@ -238,7 +238,7 @@ static_cast<char> (ival); } } - retval = octave_value (chm, 1); + retval = octave_value (chm, true, type); } return retval; @@ -368,7 +368,8 @@ if (! is.read (X_CAST (char *, &tmp), 1)) return false; - read_doubles (is, m.xdata(), X_CAST (save_type, tmp), nz, swap, fmt); + double *re = m.xdata (); + read_doubles (is, re, X_CAST (save_type, tmp), nz, swap, fmt); if (error_state || ! is) return false; @@ -678,7 +679,8 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nc + 1 || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nc + 1 || + static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -713,7 +715,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid); @@ -748,7 +750,7 @@ H5Sget_simple_extent_dims (space_hid, hdims, maxdims); - if (hdims[0] != nz || hdims[1] != 1) + if (static_cast<int> (hdims[0]) != nz || static_cast<int> (hdims[1]) != 1) { H5Sclose (space_hid); H5Dclose (data_hid);
--- a/src/ov-re-sparse.h Fri Apr 29 05:20:01 2005 +0000 +++ b/src/ov-re-sparse.h Fri Apr 29 13:04:25 2005 +0000 @@ -113,7 +113,7 @@ streamoff_array streamoff_array_value (void) const; - octave_value convert_to_str_internal (bool pad, bool force) const; + octave_value convert_to_str_internal (bool pad, bool force, char type) const; #if 0 int write (octave_stream& os, int block_size,
--- a/src/sparse-xdiv.cc Fri Apr 29 05:20:01 2005 +0000 +++ b/src/sparse-xdiv.cc Fri Apr 29 13:04:25 2005 +0000 @@ -127,41 +127,47 @@ // -*- 1 -*- Matrix -xdiv (const Matrix& a, const SparseMatrix& b) +xdiv (const Matrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return Matrix (); Matrix atmp = a.transpose (); SparseMatrix btmp = b.transpose (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) { double rcond = 0.0; - Matrix result = btmp.solve (atmp, info, rcond, + Matrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return Matrix (result.transpose ()); + { + typ = btyp.transpose (); + return Matrix (result.transpose ()); + } } octave_idx_type rank; Matrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.transpose (); } // -*- 2 -*- ComplexMatrix -xdiv (const Matrix& a, const SparseComplexMatrix& b) +xdiv (const Matrix& a, const SparseComplexMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); Matrix atmp = a.transpose (); SparseComplexMatrix btmp = b.hermitian (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -169,27 +175,32 @@ double rcond = 0.0; ComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; ComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } // -*- 3 -*- ComplexMatrix -xdiv (const ComplexMatrix& a, const SparseMatrix& b) +xdiv (const ComplexMatrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); ComplexMatrix atmp = a.hermitian (); SparseMatrix btmp = b.transpose (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -197,27 +208,32 @@ double rcond = 0.0; ComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; ComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } // -*- 4 -*- ComplexMatrix -xdiv (const ComplexMatrix& a, const SparseComplexMatrix& b) +xdiv (const ComplexMatrix& a, const SparseComplexMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); ComplexMatrix atmp = a.hermitian (); SparseComplexMatrix btmp = b.hermitian (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -225,55 +241,65 @@ double rcond = 0.0; ComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; ComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } // -*- 5 -*- SparseMatrix -xdiv (const SparseMatrix& a, const SparseMatrix& b) +xdiv (const SparseMatrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return SparseMatrix (); SparseMatrix atmp = a.transpose (); SparseMatrix btmp = b.transpose (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) { double rcond = 0.0; - SparseMatrix result = btmp.solve (atmp, info, rcond, + SparseMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return SparseMatrix (result.transpose ()); + { + typ = btyp.transpose (); + return SparseMatrix (result.transpose ()); + } } octave_idx_type rank; SparseMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.transpose (); } // -*- 6 -*- SparseComplexMatrix -xdiv (const SparseMatrix& a, const SparseComplexMatrix& b) +xdiv (const SparseMatrix& a, const SparseComplexMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseMatrix atmp = a.transpose (); SparseComplexMatrix btmp = b.hermitian (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -281,27 +307,32 @@ double rcond = 0.0; SparseComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } // -*- 7 -*- SparseComplexMatrix -xdiv (const SparseComplexMatrix& a, const SparseMatrix& b) +xdiv (const SparseComplexMatrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseComplexMatrix atmp = a.hermitian (); SparseMatrix btmp = b.transpose (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -309,27 +340,32 @@ double rcond = 0.0; SparseComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } // -*- 8 -*- SparseComplexMatrix -xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b) +xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, SparseType &typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseComplexMatrix atmp = a.hermitian (); SparseComplexMatrix btmp = b.hermitian (); + SparseType btyp = typ.transpose (); octave_idx_type info; if (btmp.rows () == btmp.columns ()) @@ -337,14 +373,18 @@ double rcond = 0.0; SparseComplexMatrix result - = btmp.solve (atmp, info, rcond, solve_singularity_warning); + = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); if (result_ok (info)) - return result.hermitian (); + { + typ = btyp.transpose (); + return result.hermitian (); + } } octave_idx_type rank; SparseComplexMatrix result = btmp.lssolve (atmp, info, rank); + typ = btyp.transpose (); return result.hermitian (); } @@ -452,7 +492,7 @@ // -*- 1 -*- Matrix -xleftdiv (const SparseMatrix& a, const Matrix& b) +xleftdiv (const SparseMatrix& a, const Matrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return Matrix (); @@ -463,7 +503,7 @@ double rcond = 0.0; Matrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -475,7 +515,7 @@ // -*- 2 -*- ComplexMatrix -xleftdiv (const SparseMatrix& a, const ComplexMatrix& b) +xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); @@ -486,7 +526,7 @@ double rcond = 0.0; ComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -498,7 +538,7 @@ // -*- 3 -*- SparseMatrix -xleftdiv (const SparseMatrix& a, const SparseMatrix& b) +xleftdiv (const SparseMatrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return SparseMatrix (); @@ -509,7 +549,7 @@ double rcond = 0.0; SparseMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -521,7 +561,7 @@ // -*- 4 -*- SparseComplexMatrix -xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b) +xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); @@ -532,7 +572,7 @@ double rcond = 0.0; SparseComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -544,7 +584,7 @@ // -*- 5 -*- ComplexMatrix -xleftdiv (const SparseComplexMatrix& a, const Matrix& b) +xleftdiv (const SparseComplexMatrix& a, const Matrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); @@ -555,7 +595,7 @@ double rcond = 0.0; ComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -567,7 +607,7 @@ // -*- 6 -*- ComplexMatrix -xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b) +xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); @@ -578,7 +618,7 @@ double rcond = 0.0; ComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -590,7 +630,7 @@ // -*- 7 -*- SparseComplexMatrix -xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b) +xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b, SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); @@ -601,7 +641,7 @@ double rcond = 0.0; SparseComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result; @@ -613,7 +653,8 @@ // -*- 8 -*- SparseComplexMatrix -xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b) +xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, + SparseType &typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); @@ -624,7 +665,7 @@ double rcond = 0.0; SparseComplexMatrix result - = a.solve (b, info, rcond, solve_singularity_warning); + = a.solve (typ, b, info, rcond, solve_singularity_warning); if (result_ok (info)) return result;
--- a/src/sparse-xdiv.h Fri Apr 29 05:20:01 2005 +0000 +++ b/src/sparse-xdiv.h Fri Apr 29 13:04:25 2005 +0000 @@ -24,23 +24,27 @@ #define octave_sparse_xdiv_h 1 #include "oct-cmplx.h" +#include "SparseType.h" class SparseMatrix; class SparseComplexMatrix; -extern Matrix xdiv (const Matrix& a, const SparseMatrix& b); -extern ComplexMatrix xdiv (const Matrix& a, const SparseComplexMatrix& b); -extern ComplexMatrix xdiv (const ComplexMatrix& a, const SparseMatrix& b); +extern Matrix xdiv (const Matrix& a, const SparseMatrix& b, SparseType &typ); +extern ComplexMatrix xdiv (const Matrix& a, const SparseComplexMatrix& b, + SparseType &typ); +extern ComplexMatrix xdiv (const ComplexMatrix& a, const SparseMatrix& b, + SparseType &typ); extern ComplexMatrix xdiv (const ComplexMatrix& a, - const SparseComplexMatrix& b); + const SparseComplexMatrix& b, SparseType &typ); -extern SparseMatrix xdiv (const SparseMatrix& a, const SparseMatrix& b); +extern SparseMatrix xdiv (const SparseMatrix& a, const SparseMatrix& b, + SparseType &typ); extern SparseComplexMatrix xdiv (const SparseMatrix& a, - const SparseComplexMatrix& b); + const SparseComplexMatrix& b, SparseType &typ); extern SparseComplexMatrix xdiv (const SparseComplexMatrix& a, - const SparseMatrix& b); + const SparseMatrix& b, SparseType &typ); extern SparseComplexMatrix xdiv (const SparseComplexMatrix& a, - const SparseComplexMatrix& b); + const SparseComplexMatrix& b, SparseType &typ); extern Matrix x_el_div (double a, const SparseMatrix& b); extern ComplexMatrix x_el_div (double a, const SparseComplexMatrix& b); @@ -48,19 +52,23 @@ extern ComplexMatrix x_el_div (const Complex a, const SparseComplexMatrix& b); -extern Matrix xleftdiv (const SparseMatrix& a, const Matrix& b); -extern ComplexMatrix xleftdiv (const SparseMatrix& a, const ComplexMatrix& b); -extern ComplexMatrix xleftdiv (const SparseComplexMatrix& a, const Matrix& b); +extern Matrix xleftdiv (const SparseMatrix& a, const Matrix& b, + SparseType& typ); +extern ComplexMatrix xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, + SparseType &typ); +extern ComplexMatrix xleftdiv (const SparseComplexMatrix& a, const Matrix& b, + SparseType &typ); extern ComplexMatrix xleftdiv (const SparseComplexMatrix& a, - const ComplexMatrix& b); + const ComplexMatrix& b, SparseType &typ); -extern SparseMatrix xleftdiv (const SparseMatrix& a, const SparseMatrix& b); +extern SparseMatrix xleftdiv (const SparseMatrix& a, const SparseMatrix& b, + SparseType &typ); extern SparseComplexMatrix xleftdiv (const SparseMatrix& a, - const SparseComplexMatrix& b); + const SparseComplexMatrix& b, SparseType &typ); extern SparseComplexMatrix xleftdiv (const SparseComplexMatrix& a, - const SparseMatrix& b); + const SparseMatrix& b, SparseType &typ); extern SparseComplexMatrix xleftdiv (const SparseComplexMatrix& a, - const SparseComplexMatrix& b); + const SparseComplexMatrix& b, SparseType &typ); #endif