Mercurial > hg > octave-lyh
changeset 5876:565d0cd4d9d0
[project @ 2006-07-01 19:42:06 by dbateman]
author | dbateman |
---|---|
date | Sat, 01 Jul 2006 19:42:07 +0000 |
parents | f6ddc0ee2315 |
children | 50d43cdbec80 |
files | liboctave/CSparse.cc liboctave/ChangeLog liboctave/Sparse-op-defs.h liboctave/dSparse.cc |
diffstat | 4 files changed, 118 insertions(+), 30 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -750,6 +750,14 @@ octave_idx_type colXp = retval.xcidx(i); octave_idx_type colUp = cidx(j); octave_idx_type rpX, rpU; + + if (cidx(j) == cidx(j+1)) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } + do { OCTAVE_QUIT; @@ -769,12 +777,19 @@ } while ((rpX<j) && (rpU<j) && (colXp<cx) && (colUp<nz)); + // get A(m,m) - colUp = cidx(j+1) - 1; + if (typ == MatrixType::Upper) + colUp = cidx(j+1) - 1; + else + colUp = cidx(j) - 1; Complex pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + if (pivot == 0. || colUp != j) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } if (v != 0.) { @@ -791,10 +806,17 @@ } // get A(m,m) - octave_idx_type colUp = cidx(i+1) - 1; + octave_idx_type colUp; + if (typ == MatrixType::Upper) + colUp = cidx(i+1) - 1; + else + colUp = cidx(i) - 1; Complex pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) ("division by zero"); + if (pivot == 0. || colUp != i) + { + (*current_liboctave_error_handler) ("division by zero"); + goto inverse_singular; + } if (pivot != 1.0) for (octave_idx_type j = cx_colstart; j < cx; j++) @@ -852,20 +874,35 @@ } // get A(m,m) - Complex pivot = data(cidx(jidx+1) - 1); + Complex pivot; + if (typ == MatrixType::Permuted_Upper) + pivot = data(cidx(jidx+1) - 1); + else + pivot = data(cidx(jidx) - 1); if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } work[j] = v / pivot; } // get A(m,m) - octave_idx_type colUp = cidx(perm[iidx]+1) - 1; - Complex pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + octave_idx_type colUp; + if (typ == MatrixType::Permuted_Upper) + colUp = cidx(perm[iidx]+1) - 1; + else + colUp = cidx(perm[iidx]) - 1; + + Complex pivot = data(colUp); + if (pivot == 0.) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } octave_idx_type new_cx = cx; for (octave_idx_type j = iidx; j < nr; j++) @@ -916,6 +953,9 @@ } return retval; + + inverse_singular: + return SparseComplexMatrix(); } SparseComplexMatrix
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,12 @@ +2006-07-01 David Bateman <dbateman@free.fr> + + * dSparse.cc (tinverse): Check for rows with no elements and zero + elements on the diagonal. Allow both Upper and Lower triangular + matrices to be treated. + * CSparse.cc (tinverse): ditto. + * Sparse-op-defs.h (SPARSE_SPARSE_MUL): Take into account 64-bit + constant assignment. + 2006-06-30 John W. Eaton <jwe@octave.org> * lo-sysdep.cc (octave_chdir): Perform tilde expansion here.
--- a/liboctave/Sparse-op-defs.h +++ b/liboctave/Sparse-op-defs.h @@ -1548,7 +1548,7 @@ else \ { \ OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \ - RET_TYPE retval (nr, a_nc, 0); \ + RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \ for (octave_idx_type i = 0; i < nr; i++) \ w[i] = 0; \ retval.xcidx(0) = 0; \
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -827,6 +827,14 @@ octave_idx_type colXp = retval.xcidx(i); octave_idx_type colUp = cidx(j); octave_idx_type rpX, rpU; + + if (cidx(j) == cidx(j+1)) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } + do { OCTAVE_QUIT; @@ -847,11 +855,17 @@ (colXp<cx) && (colUp<nz)); // get A(m,m) - colUp = cidx(j+1) - 1; + if (typ == MatrixType::Upper) + colUp = cidx(j+1) - 1; + else + colUp = cidx(j) - 1; double pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + if (pivot == 0. || colUp != j) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } if (v != 0.) { @@ -868,10 +882,17 @@ } // get A(m,m) - octave_idx_type colUp = cidx(i+1) - 1; + octave_idx_type colUp; + if (typ == MatrixType::Upper) + colUp = cidx(i+1) - 1; + else + colUp = cidx(i) - 1; double pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) ("division by zero"); + if (pivot == 0. || colUp != i) + { + (*current_liboctave_error_handler) ("division by zero"); + goto inverse_singular; + } if (pivot != 1.0) for (octave_idx_type j = cx_colstart; j < cx; j++) @@ -929,20 +950,35 @@ } // get A(m,m) - double pivot = data(cidx(jidx+1) - 1); + double pivot; + if (typ == MatrixType::Permuted_Upper) + pivot = data(cidx(jidx+1) - 1); + else + pivot = data(cidx(jidx) - 1); if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } work[j] = v / pivot; } // get A(m,m) - octave_idx_type colUp = cidx(perm[iidx]+1) - 1; + octave_idx_type colUp; + if (typ == MatrixType::Permuted_Upper) + colUp = cidx(perm[iidx]+1) - 1; + else + colUp = cidx(perm[iidx]) - 1; + double pivot = data(colUp); - if (pivot == 0.) - (*current_liboctave_error_handler) - ("division by zero"); + if (pivot == 0.) + { + (*current_liboctave_error_handler) + ("division by zero"); + goto inverse_singular; + } octave_idx_type new_cx = cx; for (octave_idx_type j = iidx; j < nr; j++) @@ -993,6 +1029,9 @@ } return retval; + + inverse_singular: + return SparseMatrix(); } SparseMatrix