Mercurial > hg > octave-thorsten
changeset 5630:512d0d11ae39
[project @ 2006-02-20 22:05:30 by dbateman]
| author | dbateman |
|---|---|
| date | Mon, 20 Feb 2006 22:05:32 +0000 |
| parents | 489a475073d7 |
| children | 7171d19706df |
| files | liboctave/CSparse.cc liboctave/ChangeLog liboctave/Sparse.cc liboctave/SparseType.cc liboctave/dSparse.cc src/ChangeLog src/DLD-FUNCTIONS/matrix_type.cc src/pt-mat.cc test/ChangeLog test/build_sparse_tests.sh |
| diffstat | 10 files changed, 1159 insertions(+), 886 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/CSparse.cc +++ b/liboctave/CSparse.cc @@ -1107,16 +1107,18 @@ } ComplexMatrix -SparseComplexMatrix::dsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, +SparseComplexMatrix::dsolve (SparseType &mattype, const Matrix& b, + octave_idx_type& err, double& rcond, solve_singularity_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1128,18 +1130,19 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - retval.resize (b.rows (), b.cols()); + retval.resize (nc, b.cols(), Complex(0.,0.)); if (typ == SparseType::Diagonal) for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(i,j) / data (i); + for (octave_idx_type i = 0; i < nm; i++) + retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(ridx(i),j) / data (i); + for (octave_idx_type k = 0; k < nc; k++) + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) + retval(k,j) = b(ridx(i),j) / data (i); double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs(data(i)); if (tmp > dmax) @@ -1158,15 +1161,17 @@ SparseComplexMatrix SparseComplexMatrix::dsolve (SparseType &mattype, const SparseMatrix& b, - octave_idx_type& err, double& rcond, solve_singularity_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1178,10 +1183,9 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; @@ -1198,27 +1202,28 @@ else for (octave_idx_type j = 0; j < b.cols(); j++) { - for (octave_idx_type i = 0; i < nr; i++) - { - bool found = false; - octave_idx_type k; - for (k = b.cidx(j); k < b.cidx(j+1); k++) - if (ridx(i) == b.ridx(k)) + for (octave_idx_type l = 0; l < nc; l++) + for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) + { + bool found = false; + octave_idx_type k; + for (k = b.cidx(j); k < b.cidx(j+1); k++) + if (ridx(i) == b.ridx(k)) + { + found = true; + break; + } + if (found) { - found = true; - break; + retval.xridx (ii) = l; + retval.xdata (ii++) = b.data(k) / data (i); } - if (found) - { - retval.xridx (ii) = i; - retval.xdata (ii++) = b.data(k) / data (i); - } - } + } retval.xcidx(j+1) = ii; } double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs(data(i)); if (tmp > dmax) @@ -1237,15 +1242,17 @@ ComplexMatrix SparseComplexMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b, - octave_idx_type& err, double& rcond, solve_singularity_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1257,15 +1264,16 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - retval.resize (b.rows (), b.cols()); + retval.resize (nc, b.cols(), Complex(0.,0.)); if (typ == SparseType::Diagonal) for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(ridx(i),j) / data (i); + for (octave_idx_type k = 0; k < nc; k++) + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) + retval(k,j) = b(ridx(i),j) / data (i); double dmax = 0., dmin = octave_Inf; for (octave_idx_type i = 0; i < nr; i++) @@ -1287,16 +1295,17 @@ SparseComplexMatrix SparseComplexMatrix::dsolve (SparseType &mattype, const SparseComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1308,10 +1317,9 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; @@ -1328,27 +1336,28 @@ else for (octave_idx_type j = 0; j < b.cols(); j++) { - for (octave_idx_type i = 0; i < nr; i++) - { - bool found = false; - octave_idx_type k; - for (k = b.cidx(j); k < b.cidx(j+1); k++) - if (ridx(i) == b.ridx(k)) + for (octave_idx_type l = 0; l < nc; l++) + for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) + { + bool found = false; + octave_idx_type k; + for (k = b.cidx(j); k < b.cidx(j+1); k++) + if (ridx(i) == b.ridx(k)) + { + found = true; + break; + } + if (found) { - found = true; - break; + retval.xridx (ii) = l; + retval.xdata (ii++) = b.data(k) / data (i); } - if (found) - { - retval.xridx (ii) = i; - retval.xdata (ii++) = b.data(k) / data (i); - } - } + } retval.xcidx(j+1) = ii; } double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs(data(i)); if (tmp > dmax) @@ -1366,17 +1375,18 @@ } ComplexMatrix -SparseComplexMatrix::utsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, - double& rcond, - solve_singularity_handler sing_handler) const +SparseComplexMatrix::utsolve (SparseType &mattype, const Matrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1390,11 +1400,11 @@ { double anorm = 0.; double ainvnorm = 0.; - octave_idx_type b_cols = b.cols (); + octave_idx_type b_nc = b.cols (); rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1405,16 +1415,18 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (nr, b_cols); + retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nr); - for (octave_idx_type j = 0; j < b_cols; j++) + for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); - - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1437,12 +1449,12 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1477,15 +1489,19 @@ } else { - retval = ComplexMatrix (b); - Complex *x_vec = retval.fortran_vec (); - - for (octave_idx_type j = 0; j < b_cols; j++) + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + retval.resize (nc, b_nc); + + for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k+1)-1) != k) { @@ -1493,22 +1509,22 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k+1)-1); - x_vec[k+offset] = tmp; + Complex tmp = work[k] / data(cidx(k+1)-1); + work[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1575,16 +1591,17 @@ SparseComplexMatrix SparseComplexMatrix::utsolve (SparseType &mattype, const SparseMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1601,7 +1618,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1610,10 +1627,9 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; @@ -1621,20 +1637,20 @@ if (typ == SparseType::Permuted_Upper) { octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - - OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); + for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1660,7 +1676,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -1672,7 +1688,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; @@ -1684,7 +1700,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1719,16 +1735,16 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { if (work[k] != 0.) { @@ -1751,7 +1767,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -1763,7 +1779,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -1775,7 +1791,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1841,16 +1857,17 @@ ComplexMatrix SparseComplexMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1868,7 +1885,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1879,16 +1896,18 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (nr, b_nc); + retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); - - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1911,12 +1930,12 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1951,15 +1970,19 @@ } else { - retval = b; - Complex *x_vec = retval.fortran_vec (); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + retval.resize (nc, b_nc); for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k+1)-1) != k) { @@ -1967,22 +1990,22 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k+1)-1); - x_vec[k+offset] = tmp; + Complex tmp = work[k] / data(cidx(k+1)-1); + work[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2049,16 +2072,17 @@ SparseComplexMatrix SparseComplexMatrix::utsolve (SparseType &mattype, const SparseComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2075,7 +2099,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2084,10 +2108,9 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; @@ -2095,20 +2118,20 @@ if (typ == SparseType::Permuted_Upper) { octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - - OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); + for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -2134,7 +2157,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2146,7 +2169,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; @@ -2158,7 +2181,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2193,11 +2216,11 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); @@ -2225,7 +2248,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2237,7 +2260,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -2249,7 +2272,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2315,16 +2338,18 @@ } ComplexMatrix -SparseComplexMatrix::ltsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler sing_handler) const +SparseComplexMatrix::ltsolve (SparseType &mattype, const Matrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2338,11 +2363,11 @@ { double anorm = 0.; double ainvnorm = 0.; - octave_idx_type b_cols = b.cols (); + octave_idx_type b_nc = b.cols (); rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2353,16 +2378,18 @@ if (typ == SparseType::Permuted_Lower) { - retval.resize (b.rows (), b.cols ()); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + retval.resize (nc, b_nc); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); - for (octave_idx_type j = 0; j < b_cols; j++) + for (octave_idx_type j = 0; j < b_nc; j++) { + for (octave_idx_type i = 0; i < nm; i++) + work[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) work[perm[i]] = b(i,j); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2395,19 +2422,19 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2435,7 +2462,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2446,15 +2473,18 @@ } else { - retval = ComplexMatrix (b); - Complex *x_vec = retval.fortran_vec (); - - for (octave_idx_type j = 0; j < b_cols; j++) + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + retval.resize (nc, b_nc, 0.); + + for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + for (octave_idx_type k = 0; k < nc; k++) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k)) != k) { @@ -2462,29 +2492,28 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k)); - x_vec[k+offset] = tmp; + Complex tmp = work[k] / data(cidx(k)); + work[k] = tmp; for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -2499,7 +2528,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2545,16 +2574,18 @@ SparseComplexMatrix SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); + err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2571,7 +2602,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2580,27 +2611,26 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == SparseType::Permuted_Lower) { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); 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++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) work[perm[b.ridx(i)]] = b.data(i); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2636,7 +2666,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2648,7 +2678,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -2660,14 +2690,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2695,7 +2725,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2706,16 +2736,16 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2738,7 +2768,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2750,7 +2780,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -2762,14 +2792,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -2784,7 +2814,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2831,16 +2861,17 @@ ComplexMatrix SparseComplexMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2858,7 +2889,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2869,16 +2900,18 @@ if (typ == SparseType::Permuted_Lower) { - retval.resize (b.rows (), b.cols ()); - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + retval.resize (nc, b_nc); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { + for (octave_idx_type i = 0; i < nm; i++) + work[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) work[perm[i]] = b(i,j); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2911,19 +2944,19 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2951,7 +2984,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -2962,15 +2995,20 @@ } else { - retval = b; - Complex *x_vec = retval.fortran_vec (); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); + retval.resize (nc, b_nc, 0.); + for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = 0; k < nc; k++) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k)) != k) { @@ -2978,29 +3016,29 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k)); - x_vec[k+offset] = tmp; + Complex tmp = work[k] / data(cidx(k)); + work[k] = tmp; for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (Complex, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -3015,7 +3053,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -3062,16 +3100,17 @@ SparseComplexMatrix SparseComplexMatrix::ltsolve (SparseType &mattype, const SparseComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -3088,7 +3127,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -3097,27 +3136,26 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == SparseType::Permuted_Lower) { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); 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++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) work[perm[b.ridx(i)]] = b.data(i); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -3153,7 +3191,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -3165,7 +3203,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -3177,14 +3215,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -3212,7 +3250,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -3223,16 +3261,16 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -3255,7 +3293,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -3267,7 +3305,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -3279,14 +3317,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -3301,7 +3339,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs(work[i]); work[i] = 0.; @@ -4122,7 +4160,7 @@ Array<octave_idx_type> ipvt (nr); octave_idx_type *pipvt = ipvt.fortran_vec (); - F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, tmp_data, + F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, tmp_data, ldm, pipvt, err)); if (f77_exception_encountered)
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,17 @@ +2006-02-20 David Bateman <dbateman@free.fr> + + * dSparse.cc (dsolve, utsolve, ltsolve): Remove restriction that + matrix must be square in diagonal, permuted diagonal, triangular + and permuted triangular back/forward substitution code. Change + ambiguous use of no. rows and columns. + * CSParse.cc (dsolve, utsolve, ltsolve): ditto. + * SparseType.cc (SparseType::SparseType(const SparseMatrix&), + SparseType::SparseType(const SparseComplexMatrix&)): Recognize + rectangular diagonal, permuted diagonal, triangular and permuted + triangular matrices. + * Sparse.cc (Sparse<T>::Sparse (octave_idx_type, octave_idx_type, T)): + Treat case where third argument is zero. + 2006-02-15 John W. Eaton <jwe@octave.org> * kpse.cc: Do define ST_NLINK_TRICK for Cygwin systems.
--- a/liboctave/Sparse.cc +++ b/liboctave/Sparse.cc @@ -222,20 +222,29 @@ template <class T> Sparse<T>::Sparse (octave_idx_type nr, octave_idx_type nc, T val) - : rep (new typename Sparse<T>::SparseRep (nr, nc, nr*nc)), - dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) + : dimensions (dim_vector (nr, nc)), idx (0), idx_count (0) { - - octave_idx_type ii = 0; - xcidx (0) = 0; - for (octave_idx_type j = 0; j < nc; j++) + if (val != T ()) { - for (octave_idx_type i = 0; i < nr; i++) + rep = new typename Sparse<T>::SparseRep (nr, nc, nr*nc); + + octave_idx_type ii = 0; + xcidx (0) = 0; + for (octave_idx_type j = 0; j < nc; j++) { - xdata (ii) = val; - xridx (ii++) = i; - } - xcidx (j+1) = ii; + for (octave_idx_type i = 0; i < nr; i++) + { + xdata (ii) = val; + xridx (ii++) = i; + } + xcidx (j+1) = ii; + } + } + else + { + rep = new typename Sparse<T>::SparseRep (nr, nc, 0); + for (octave_idx_type j = 0; j < nc+1; j++) + xcidx(j) = 0; } }
--- a/liboctave/SparseType.cc +++ b/liboctave/SparseType.cc @@ -55,6 +55,7 @@ { octave_idx_type nrows = a.rows (); octave_idx_type ncols = a.cols (); + octave_idx_type nm = (ncols < nrows ? ncols : nrows); octave_idx_type nnz = a.nzmax (); if (Voctave_sparse_controls.get_key ("spumoni") != 0.) @@ -63,202 +64,211 @@ nperm = 0; - if (nrows != ncols) - typ = SparseType::Rectangular; - else + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); + bool maybe_hermitian = false; + typ = SparseType::Full; + + if (nnz == nm) { - sp_bandden = Voctave_sparse_controls.get_key ("bandden"); - bool maybe_hermitian = false; - typ = SparseType::Full; - - if (nnz == ncols) + matrix_type tmp_typ = SparseType::Diagonal; + octave_idx_type i; + // Maybe the matrix is diagonal + for (i = 0; i < nm; i++) { - matrix_type tmp_typ = SparseType::Diagonal; - octave_idx_type i; - // Maybe the matrix is diagonal - for (i = 0; i < ncols; i++) + if (a.cidx(i+1) != a.cidx(i) + 1) + { + tmp_typ = SparseType::Full; + break; + } + if (a.ridx(i) != i) { - if (a.cidx(i+1) != a.cidx(i) + 1) + tmp_typ = SparseType::Permuted_Diagonal; + break; + } + } + + if (tmp_typ == SparseType::Permuted_Diagonal) + { + bool found [nrows]; + + for (octave_idx_type j = 0; j < i; j++) + found [j] = true; + for (octave_idx_type j = i; j < nrows; j++) + found [j] = false; + + for (octave_idx_type j = i; j < nm; j++) + { + if ((a.cidx(j+1) > a.cidx(j) + 1) || + ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)])) { - tmp_typ = Full; - break; - } - if (a.ridx(i) != i) - { - tmp_typ = SparseType::Permuted_Diagonal; + tmp_typ = SparseType::Full; break; } + found [a.ridx(j)] = true; } - - if (tmp_typ == SparseType::Permuted_Diagonal) - { - bool found [ncols]; + } + typ = tmp_typ; + } - for (octave_idx_type j = 0; j < i; j++) - found [j] = true; - for (octave_idx_type j = i; j < ncols; j++) - found [j] = false; - - for (octave_idx_type j = i; j < ncols; j++) - { - if ((a.cidx(j+1) != a.cidx(j) + 1) || found [a.ridx(j)]) - { - tmp_typ = Full; - break; - } - found [a.ridx(j)] = true; - } - } - typ = tmp_typ; - } - - if (typ == SparseType::Full) + if (typ == SparseType::Full) + { + // Search for banded, upper and lower triangular matrices + bool singular = false; + upper_band = 0; + lower_band = 0; + for (octave_idx_type j = 0; j < ncols; j++) { - // Search for banded, upper and lower triangular matrices - bool singular = false; - upper_band = 0; - lower_band = 0; - for (octave_idx_type j = 0; j < ncols; j++) + bool zero_on_diagonal = false; + if (j < nrows) { - bool zero_on_diagonal = true; + zero_on_diagonal = true; for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) if (a.ridx(i) == j) { zero_on_diagonal = false; break; } - - if (zero_on_diagonal) - { - singular = true; - break; - } + } - if (a.cidx(j+1) - a.cidx(j) > 0) - { - octave_idx_type ru = a.ridx(a.cidx(j)); - octave_idx_type rl = a.ridx(a.cidx(j+1)-1); - - if (j - ru > upper_band) - upper_band = j - ru; - - if (rl - j > lower_band) - lower_band = rl - j; - } + if (zero_on_diagonal) + { + singular = true; + break; } - if (!singular) + if (a.cidx(j+1) != a.cidx(j)) { - bandden = double (nnz) / - (double (ncols) * (double (lower_band) + - double (upper_band)) - - 0.5 * double (upper_band + 1) * double (upper_band) - - 0.5 * double (lower_band + 1) * double (lower_band)); + octave_idx_type ru = a.ridx(a.cidx(j)); + octave_idx_type rl = a.ridx(a.cidx(j+1)-1); + + if (j - ru > upper_band) + upper_band = j - ru; + + if (rl - j > lower_band) + lower_band = rl - j; + } + } - if (sp_bandden != 1. && bandden > sp_bandden) - { - if (upper_band == 1 && lower_band == 1) - typ = SparseType::Tridiagonal; - else - typ = SparseType::Banded; + if (!singular) + { + bandden = double (nnz) / + (double (ncols) * (double (lower_band) + + double (upper_band)) - + 0.5 * double (upper_band + 1) * double (upper_band) - + 0.5 * double (lower_band + 1) * double (lower_band)); + + if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden) + { + if (upper_band == 1 && lower_band == 1) + typ = SparseType::Tridiagonal; + else + typ = SparseType::Banded; - octave_idx_type nnz_in_band = (upper_band + lower_band + 1) * nrows - - (1 + upper_band) * upper_band / 2 - - (1 + lower_band) * lower_band / 2; - if (nnz_in_band == nnz) - dense = true; - else - dense = false; + octave_idx_type nnz_in_band = + (upper_band + lower_band + 1) * nrows - + (1 + upper_band) * upper_band / 2 - + (1 + lower_band) * lower_band / 2; + if (nnz_in_band == nnz) + dense = true; + else + dense = false; + } + else if (upper_band == 0) + typ = SparseType::Lower; + else if (lower_band == 0) + typ = SparseType::Upper; + + if (upper_band == lower_band && nrows == ncols) + maybe_hermitian = true; + } + + if (typ == SparseType::Full) + { + // Search for a permuted triangular matrix, and test if + // permutation is singular + + // XXX FIXME XXX + // Perhaps this should be based on a dmperm algorithm + bool found = false; + + nperm = ncols; + perm = new octave_idx_type [nperm]; + + for (octave_idx_type i = 0; i < nm; 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; + } } - else if (upper_band == 0) - typ = SparseType::Lower; - else if (lower_band == 0) - typ = SparseType::Upper; - if (upper_band == lower_band) - maybe_hermitian = true; + if (!found) + break; } - if (typ == SparseType::Full) + if (found) + typ = SparseType::Permuted_Upper; + else { - // Search for a permuted triangular matrix, and test if - // permutation is singular - - // XXX FIXME XXX Perhaps this should be based on a dmperm algorithm - bool found = false; - + delete [] perm; nperm = nrows; perm = new octave_idx_type [nperm]; + OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, 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; + 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) - typ = SparseType::Permuted_Upper; - else - { - OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nperm); - - for (octave_idx_type i = 0; i < nperm; i++) + for (octave_idx_type i = 0; i < nm; i++) + if (tmp[i] == -1) { - perm [i] = -1; - tmp [i] = -1; + found = false; + break; } - 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; - } + if (found) + typ = SparseType::Permuted_Lower; + else + { + delete [] perm; + nperm = 0; } } } - if (maybe_hermitian && (typ == Full || typ == Tridiagonal || - typ == Banded)) + if (typ == SparseType::Full && ncols != nrows) + typ = SparseType::Rectangular; + + if (maybe_hermitian && (typ == SparseType::Full || + typ == SparseType::Tridiagonal || + typ == SparseType::Banded)) { // Check for symmetry, with positive real diagonal, which // has a very good chance of being symmetric positive @@ -311,12 +321,12 @@ if (is_herm) { - if (typ == Full) - typ = Hermitian; - else if (typ == Banded) - typ = Banded_Hermitian; + if (typ == SparseType::Full) + typ = SparseType::Hermitian; + else if (typ == SparseType::Banded) + typ = SparseType::Banded_Hermitian; else - typ = Tridiagonal_Hermitian; + typ = SparseType::Tridiagonal_Hermitian; } } } @@ -326,6 +336,7 @@ { octave_idx_type nrows = a.rows (); octave_idx_type ncols = a.cols (); + octave_idx_type nm = (ncols < nrows ? ncols : nrows); octave_idx_type nnz = a.nzmax (); if (Voctave_sparse_controls.get_key ("spumoni") != 0.) @@ -334,202 +345,211 @@ nperm = 0; - if (nrows != ncols) - typ = SparseType::Rectangular; - else + sp_bandden = Voctave_sparse_controls.get_key ("bandden"); + bool maybe_hermitian = false; + typ = SparseType::Full; + + if (nnz == nm) { - sp_bandden = Voctave_sparse_controls.get_key ("bandden"); - bool maybe_hermitian = false; - typ = SparseType::Full; - - if (nnz == ncols) + matrix_type tmp_typ = SparseType::Diagonal; + octave_idx_type i; + // Maybe the matrix is diagonal + for (i = 0; i < nm; i++) { - matrix_type tmp_typ = SparseType::Diagonal; - octave_idx_type i; - // Maybe the matrix is diagonal - for (i = 0; i < ncols; i++) + if (a.cidx(i+1) != a.cidx(i) + 1) + { + tmp_typ = SparseType::Full; + break; + } + if (a.ridx(i) != i) { - if (a.cidx(i+1) != a.cidx(i) + 1) + tmp_typ = SparseType::Permuted_Diagonal; + break; + } + } + + if (tmp_typ == SparseType::Permuted_Diagonal) + { + bool found [nrows]; + + for (octave_idx_type j = 0; j < i; j++) + found [j] = true; + for (octave_idx_type j = i; j < nrows; j++) + found [j] = false; + + for (octave_idx_type j = i; j < nm; j++) + { + if ((a.cidx(j+1) > a.cidx(j) + 1) || + ((a.cidx(j+1) == a.cidx(j) + 1) && found [a.ridx(j)])) { - tmp_typ = Full; - break; - } - if (a.ridx(i) != i) - { - tmp_typ = SparseType::Permuted_Diagonal; + tmp_typ = SparseType::Full; break; } + found [a.ridx(j)] = true; } - - if (tmp_typ == SparseType::Permuted_Diagonal) - { - bool found [ncols]; + } + typ = tmp_typ; + } - for (octave_idx_type j = 0; j < i; j++) - found [j] = true; - for (octave_idx_type j = i; j < ncols; j++) - found [j] = false; - - for (octave_idx_type j = i; j < ncols; j++) - { - if ((a.cidx(j+1) != a.cidx(j) + 1) || found [a.ridx(j)]) - { - tmp_typ = Full; - break; - } - found [a.ridx(j)] = true; - } - } - typ = tmp_typ; - } - - if (typ == Full) + if (typ == SparseType::Full) + { + // Search for banded, upper and lower triangular matrices + bool singular = false; + upper_band = 0; + lower_band = 0; + for (octave_idx_type j = 0; j < ncols; j++) { - // Search for banded, upper and lower triangular matrices - bool singular = false; - upper_band = 0; - lower_band = 0; - for (octave_idx_type j = 0; j < ncols; j++) + bool zero_on_diagonal = false; + if (j < nrows) { - bool zero_on_diagonal = true; + zero_on_diagonal = true; for (octave_idx_type i = a.cidx(j); i < a.cidx(j+1); i++) if (a.ridx(i) == j) { zero_on_diagonal = false; break; } - - if (zero_on_diagonal) - { - singular = true; - break; - } + } - if (a.cidx(j+1) - a.cidx(j) > 0) - { - octave_idx_type ru = a.ridx(a.cidx(j)); - octave_idx_type rl = a.ridx(a.cidx(j+1)-1); - - if (j - ru > upper_band) - upper_band = j - ru; - - if (rl - j > lower_band) - lower_band = rl - j; - } + if (zero_on_diagonal) + { + singular = true; + break; } - if (!singular) + if (a.cidx(j+1) != a.cidx(j)) { - bandden = double (nnz) / - (double (ncols) * (double (lower_band) + - double (upper_band)) - - 0.5 * double (upper_band + 1) * double (upper_band) - - 0.5 * double (lower_band + 1) * double (lower_band)); + octave_idx_type ru = a.ridx(a.cidx(j)); + octave_idx_type rl = a.ridx(a.cidx(j+1)-1); + + if (j - ru > upper_band) + upper_band = j - ru; + + if (rl - j > lower_band) + lower_band = rl - j; + } + } - if (sp_bandden != 1. && bandden > sp_bandden) - { - if (upper_band == 1 && lower_band == 1) - typ = SparseType::Tridiagonal; - else - typ = SparseType::Banded; + if (!singular) + { + bandden = double (nnz) / + (double (ncols) * (double (lower_band) + + double (upper_band)) - + 0.5 * double (upper_band + 1) * double (upper_band) - + 0.5 * double (lower_band + 1) * double (lower_band)); + + if (nrows == ncols && sp_bandden != 1. && bandden > sp_bandden) + { + if (upper_band == 1 && lower_band == 1) + typ = SparseType::Tridiagonal; + else + typ = SparseType::Banded; - octave_idx_type nnz_in_band = (upper_band + lower_band + 1) * nrows - - (1 + upper_band) * upper_band / 2 - - (1 + lower_band) * lower_band / 2; - if (nnz_in_band == nnz) - dense = true; - else - dense = false; + octave_idx_type nnz_in_band = + (upper_band + lower_band + 1) * nrows - + (1 + upper_band) * upper_band / 2 - + (1 + lower_band) * lower_band / 2; + if (nnz_in_band == nnz) + dense = true; + else + dense = false; + } + else if (upper_band == 0) + typ = SparseType::Lower; + else if (lower_band == 0) + typ = SparseType::Upper; + + if (upper_band == lower_band && nrows == ncols) + maybe_hermitian = true; + } + + if (typ == SparseType::Full) + { + // Search for a permuted triangular matrix, and test if + // permutation is singular + + // XXX FIXME XXX + // Perhaps this should be based on a dmperm algorithm + bool found = false; + + nperm = ncols; + perm = new octave_idx_type [nperm]; + + for (octave_idx_type i = 0; i < nm; 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; + } } - else if (upper_band == 0) - typ = SparseType::Lower; - else if (lower_band == 0) - typ = SparseType::Upper; - if (upper_band == lower_band) - maybe_hermitian = true; + if (!found) + break; } - if (typ == Full) + if (found) + typ = SparseType::Permuted_Upper; + else { - // Search for a permuted triangular matrix, and test if - // permutation is singular - - // XXX FIXME XXX Perhaps this should be based on a dmperm algorithm - bool found = false; - + delete [] perm; nperm = nrows; perm = new octave_idx_type [nperm]; + OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, 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; + 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) - typ = SparseType::Permuted_Upper; - else - { - OCTAVE_LOCAL_BUFFER (octave_idx_type, tmp, nperm); - - for (octave_idx_type i = 0; i < nperm; i++) + for (octave_idx_type i = 0; i < nm; i++) + if (tmp[i] == -1) { - perm [i] = -1; - tmp [i] = -1; + found = false; + break; } - 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; - } + if (found) + typ = SparseType::Permuted_Lower; + else + { + delete [] perm; + nperm = 0; } } } - if (maybe_hermitian && (typ == Full || typ == Tridiagonal || - typ == Banded)) + if (typ == SparseType::Full && ncols != nrows) + typ = SparseType::Rectangular; + + if (maybe_hermitian && (typ == SparseType::Full || + typ == SparseType::Tridiagonal || + typ == SparseType::Banded)) { // Check for symmetry, with positive real diagonal, which // has a very good chance of being symmetric positive @@ -582,12 +602,12 @@ if (is_herm) { - if (typ == Full) - typ = Hermitian; - else if (typ == Banded) - typ = Banded_Hermitian; + if (typ == SparseType::Full) + typ = SparseType::Hermitian; + else if (typ == SparseType::Banded) + typ = SparseType::Banded_Hermitian; else - typ = Tridiagonal_Hermitian; + typ = SparseType::Tridiagonal_Hermitian; } } } @@ -889,9 +909,9 @@ else if (typ == SparseType::Permuted_Upper) retval.typ = SparseType::Permuted_Lower; else if (typ == SparseType::Lower) - retval.typ = Upper; + retval.typ = SparseType::Upper; else if (typ == SparseType::Permuted_Lower) - retval.typ = Permuted_Upper; + retval.typ = SparseType::Permuted_Upper; else if (typ == SparseType::Banded) { retval.upper_band = lower_band;
--- a/liboctave/dSparse.cc +++ b/liboctave/dSparse.cc @@ -1194,9 +1194,10 @@ octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1208,18 +1209,19 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - retval.resize (b.rows (), b.cols()); + retval.resize (nc, b.cols(), 0.); if (typ == SparseType::Diagonal) for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(ridx(i),j) / data (i); - + for (octave_idx_type k = 0; k < nc; k++) + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) + retval(k,j) = b(ridx(i),j) / data (i); + double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = fabs(data(i)); if (tmp > dmax) @@ -1237,16 +1239,18 @@ } SparseMatrix -SparseMatrix::dsolve (SparseType &mattype, const SparseMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler) const +SparseMatrix::dsolve (SparseType &mattype, const SparseMatrix& b, + octave_idx_type& err, + double& rcond, solve_singularity_handler) const { SparseMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1258,10 +1262,9 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseMatrix (b_nr, b_nc, b_nz); + retval = SparseMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; @@ -1278,27 +1281,28 @@ else for (octave_idx_type j = 0; j < b.cols(); j++) { - for (octave_idx_type i = 0; i < nr; i++) - { - bool found = false; - octave_idx_type k; - for (k = b.cidx(j); k < b.cidx(j+1); k++) - if (ridx(i) == b.ridx(k)) + for (octave_idx_type l = 0; l < nc; l++) + for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) + { + bool found = false; + octave_idx_type k; + for (k = b.cidx(j); k < b.cidx(j+1); k++) + if (ridx(i) == b.ridx(k)) + { + found = true; + break; + } + if (found) { - found = true; - break; + retval.xridx (ii) = l; + retval.xdata (ii++) = b.data(k) / data (i); } - if (found) - { - retval.xridx (ii) = i; - retval.xdata (ii++) = b.data(k) / data (i); - } - } + } retval.xcidx(j+1) = ii; } - + double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = fabs(data(i)); if (tmp > dmax) @@ -1316,16 +1320,18 @@ } ComplexMatrix -SparseMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler) const +SparseMatrix::dsolve (SparseType &mattype, const ComplexMatrix& b, + octave_idx_type& err, + double& rcond, solve_singularity_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1337,18 +1343,19 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - retval.resize (b.rows (), b.cols()); + retval.resize (nc, b.cols(), 0); if (typ == SparseType::Diagonal) for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(i,j) / data (i); + for (octave_idx_type i = 0; i < nm; i++) + retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols(); j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = b(ridx(i),j) / data (i); + for (octave_idx_type k = 0; k < nc; k++) + for (octave_idx_type i = cidx(k); i < cidx(k+1); i++) + retval(k,j) = b(ridx(i),j) / data (i); double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = fabs(data(i)); if (tmp > dmax) @@ -1374,9 +1381,10 @@ octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc < nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1388,10 +1396,9 @@ if (typ == SparseType::Diagonal || typ == SparseType::Permuted_Diagonal) { - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; @@ -1408,27 +1415,28 @@ else for (octave_idx_type j = 0; j < b.cols(); j++) { - for (octave_idx_type i = 0; i < nr; i++) - { - bool found = false; - octave_idx_type k; - for (k = b.cidx(j); k < b.cidx(j+1); k++) - if (ridx(i) == b.ridx(k)) + for (octave_idx_type l = 0; l < nc; l++) + for (octave_idx_type i = cidx(l); i < cidx(l+1); i++) + { + bool found = false; + octave_idx_type k; + for (k = b.cidx(j); k < b.cidx(j+1); k++) + if (ridx(i) == b.ridx(k)) + { + found = true; + break; + } + if (found) { - found = true; - break; + retval.xridx (ii) = l; + retval.xdata (ii++) = b.data(k) / data (i); } - if (found) - { - retval.xridx (ii) = i; - retval.xdata (ii++) = b.data(k) / data (i); - } - } + } retval.xcidx(j+1) = ii; } double dmax = 0., dmin = octave_Inf; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) { double tmp = fabs(data(i)); if (tmp > dmax) @@ -1446,17 +1454,18 @@ } Matrix -SparseMatrix::utsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, - double& rcond, +SparseMatrix::utsolve (SparseType &mattype, const Matrix& b, + octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { Matrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1470,11 +1479,11 @@ { double anorm = 0.; double ainvnorm = 0.; - octave_idx_type b_cols = b.cols (); + octave_idx_type b_nc = b.cols (); rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1485,16 +1494,18 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (nr, b_cols); + retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (double, work, nr); - - for (octave_idx_type j = 0; j < b_cols; j++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + + for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); - - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1517,12 +1528,12 @@ } } - for (octave_idx_type i = 0; i < nr; i++) - retval (perm[i], j) = work[i]; + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (perm[i], j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1556,15 +1567,19 @@ } else { - retval = b; - double *x_vec = retval.fortran_vec (); - - for (octave_idx_type j = 0; j < b_cols; j++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + retval.resize (nc, b_nc); + + for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k+1)-1) != k) { @@ -1572,22 +1587,22 @@ goto triangular_error; } - double tmp = x_vec[k+offset] / - data(cidx(k+1)-1); - x_vec[k+offset] = tmp; + double tmp = work[k] / data(cidx(k+1)-1); + work[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1653,16 +1668,18 @@ } SparseMatrix -SparseMatrix::utsolve (SparseType &mattype, const SparseMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler sing_handler) const +SparseMatrix::utsolve (SparseType &mattype, const SparseMatrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1679,7 +1696,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1688,10 +1705,9 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseMatrix (b_nr, b_nc, b_nz); + retval = SparseMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; @@ -1699,20 +1715,20 @@ if (typ == SparseType::Permuted_Upper) { octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (double, work, nr); - - OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); + for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1738,7 +1754,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -1750,7 +1766,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[rperm[i]] != 0.) { retval.xridx(ii) = i; @@ -1762,7 +1778,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1796,16 +1812,16 @@ } else { - OCTAVE_LOCAL_BUFFER (double, work, nr); + OCTAVE_LOCAL_BUFFER (double, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { if (work[k] != 0.) { @@ -1828,7 +1844,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -1840,7 +1856,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -1852,7 +1868,7 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -1917,16 +1933,18 @@ } ComplexMatrix -SparseMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler sing_handler) const +SparseMatrix::utsolve (SparseType &mattype, const ComplexMatrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -1944,7 +1962,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -1955,16 +1973,18 @@ if (typ == SparseType::Permuted_Upper) { - retval.resize (nr, b_nc); + retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); - OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) cwork[i] = b(i,j); - - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = nr; i < nc; i++) + cwork[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -1987,13 +2007,13 @@ } } - for (octave_idx_type i = 0; i < nr; i++) - retval (perm[i], j) = cwork[i]; + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (perm[i], j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2027,15 +2047,19 @@ } else { - retval = b; - Complex *x_vec = retval.fortran_vec (); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); + retval.resize (nc, b_nc); for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type i = 0; i < nr; i++) + cwork[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + cwork[i] = 0.; + + for (octave_idx_type k = nc-1; k >= 0; k--) { - if (x_vec[k+offset] != 0.) + if (cwork[k] != 0.) { if (ridx(cidx(k+1)-1) != k) { @@ -2043,22 +2067,23 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k+1)-1); - x_vec[k+offset] = tmp; + Complex tmp = cwork[k] / data(cidx(k+1)-1); + cwork[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2125,16 +2150,17 @@ SparseComplexMatrix SparseMatrix::utsolve (SparseType &mattype, const SparseComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2151,7 +2177,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2160,10 +2186,9 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; @@ -2171,20 +2196,20 @@ if (typ == SparseType::Permuted_Upper) { 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++) + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); + + OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); + for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) cwork[b.ridx(i)] = b.data(i); - for (octave_idx_type k = nr-1; k >= 0; k--) + for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; @@ -2210,7 +2235,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (cwork[i] != 0.) new_nnz++; @@ -2222,7 +2247,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (cwork[rperm[i]] != 0.) { retval.xridx(ii) = i; @@ -2234,8 +2259,8 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) @@ -2269,18 +2294,18 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 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 = 0; i < nm; i++) + cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - for (octave_idx_type k = nr-1; k >= 0; k--) + cwork[b.ridx(i)] = b.data(i); + + for (octave_idx_type k = nc-1; k >= 0; k--) { - if (work[k] != 0.) + if (cwork[k] != 0.) { if (ridx(cidx(k+1)-1) != k) { @@ -2288,12 +2313,12 @@ goto triangular_error; } - Complex tmp = work[k] / data(cidx(k+1)-1); - work[k] = tmp; + Complex tmp = cwork[k] / data(cidx(k+1)-1); + cwork[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - work[iidx] = work[iidx] - tmp * data(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } @@ -2301,8 +2326,8 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + for (octave_idx_type i = 0; i < nc; i++) + if (cwork[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) @@ -2313,11 +2338,11 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + for (octave_idx_type i = 0; i < nc; i++) + if (cwork[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; + retval.xdata(ii++) = cwork[i]; } retval.xcidx(j+1) = ii; } @@ -2325,32 +2350,32 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work2, nr); - for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[j] = 1.; + work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { - if (work2[k] != 0.) + if (work[k] != 0.) { - double tmp = work2[k] / data(cidx(k+1)-1); - work2[k] = tmp; + double tmp = work[k] / data(cidx(k+1)-1); + work[k] = tmp; for (octave_idx_type i = cidx(k); i < cidx(k+1)-1; i++) { octave_idx_type iidx = ridx(i); - work2[iidx] = work2[iidx] - tmp * data(i); + work[iidx] = work[iidx] - 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; @@ -2392,17 +2417,18 @@ } Matrix -SparseMatrix::ltsolve (SparseType &mattype, const Matrix& b, octave_idx_type& err, - double& rcond, +SparseMatrix::ltsolve (SparseType &mattype, const Matrix& b, + octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { Matrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2416,11 +2442,11 @@ { double anorm = 0.; double ainvnorm = 0.; - octave_idx_type b_cols = b.cols (); + octave_idx_type b_nc = b.cols (); rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2431,16 +2457,19 @@ if (typ == SparseType::Permuted_Lower) { - retval.resize (b.rows (), b.cols ()); - OCTAVE_LOCAL_BUFFER (double, work, nr); + retval.resize (nc, b_nc); + OCTAVE_LOCAL_BUFFER (double, work, nm); octave_idx_type *perm = mattype.triangular_perm (); - for (octave_idx_type j = 0; j < b_cols; j++) + for (octave_idx_type j = 0; j < b_nc; j++) { + if (nc > nr) + for (octave_idx_type i = 0; i < nm; i++) + work[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) work[perm[i]] = b(i,j); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2473,19 +2502,19 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2513,7 +2542,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -2524,15 +2553,18 @@ } else { - retval = b; - double *x_vec = retval.fortran_vec (); - - for (octave_idx_type j = 0; j < b_cols; j++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + retval.resize (nc, b_nc, 0.); + + for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type i = 0; i < nr; i++) + work[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + work[i] = 0.; + for (octave_idx_type k = 0; k < nc; k++) { - if (x_vec[k+offset] != 0.) + if (work[k] != 0.) { if (ridx(cidx(k)) != k) { @@ -2540,29 +2572,30 @@ goto triangular_error; } - double tmp = x_vec[k+offset] / - data(cidx(k)); - x_vec[k+offset] = tmp; - for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) + double tmp = work[k] / data(cidx(k)); + work[k] = tmp; + for (octave_idx_type i = cidx(k)+1; + i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = work[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -2577,7 +2610,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -2623,16 +2656,18 @@ } SparseMatrix -SparseMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler sing_handler) const +SparseMatrix::ltsolve (SparseType &mattype, const SparseMatrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2649,7 +2684,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2673,12 +2708,12 @@ for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) work[perm[b.ridx(i)]] = b.data(i); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2714,7 +2749,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2726,7 +2761,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -2738,14 +2773,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2788,12 +2823,12 @@ for (octave_idx_type j = 0; j < b_nc; j++) { - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; 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); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -2816,7 +2851,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; @@ -2828,7 +2863,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx(ii) = i; @@ -2840,14 +2875,14 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -2862,7 +2897,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -2908,16 +2943,18 @@ } ComplexMatrix -SparseMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b, octave_idx_type& err, - double& rcond, solve_singularity_handler sing_handler) const +SparseMatrix::ltsolve (SparseType &mattype, const ComplexMatrix& b, + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -2935,7 +2972,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -2946,16 +2983,18 @@ if (typ == SparseType::Permuted_Lower) { - retval.resize (b.rows (), b.cols ()); + retval.resize (nc, b_nc); 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 < nm; i++) + cwork[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) cwork[perm[i]] = b(i,j); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (cwork[k] != 0.) { @@ -2988,20 +3027,20 @@ } } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) retval (i, j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -3029,7 +3068,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -3040,15 +3079,19 @@ } else { - retval = b; - Complex *x_vec = retval.fortran_vec (); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); + retval.resize (nc, b_nc, 0.); for (octave_idx_type j = 0; j < b_nc; j++) { - octave_idx_type offset = j * nr; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type i = 0; i < nr; i++) + cwork[i] = b(i,j); + for (octave_idx_type i = nr; i < nc; i++) + cwork[i] = 0.; + + for (octave_idx_type k = 0; k < nc; k++) { - if (x_vec[k+offset] != 0.) + if (cwork[k] != 0.) { if (ridx(cidx(k)) != k) { @@ -3056,29 +3099,30 @@ goto triangular_error; } - Complex tmp = x_vec[k+offset] / - data(cidx(k)); - x_vec[k+offset] = tmp; + Complex tmp = cwork[k] / data(cidx(k)); + cwork[k] = tmp; for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - x_vec[iidx+offset] = - x_vec[iidx+offset] - tmp * data(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } + + for (octave_idx_type i = 0; i < nc; i++) + retval.xelem (i, j) = cwork[i]; } // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = j; k < nr; k++) + for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) @@ -3093,7 +3137,7 @@ } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -3140,16 +3184,17 @@ SparseComplexMatrix SparseMatrix::ltsolve (SparseType &mattype, const SparseComplexMatrix& b, - octave_idx_type& err, double& rcond, - solve_singularity_handler sing_handler) const + octave_idx_type& err, double& rcond, + solve_singularity_handler sing_handler) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); + octave_idx_type nm = (nc > nr ? nc : nr); err = 0; - if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ()) + if (nr == 0 || nc == 0 || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); else @@ -3166,7 +3211,7 @@ rcond = 0.; // Calculate the 1-norm of matrix for rcond calculation - for (octave_idx_type j = 0; j < nr; j++) + for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx(j); i < cidx(j+1); i++) @@ -3175,27 +3220,26 @@ anorm = atmp; } - octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nzmax (); - retval = SparseComplexMatrix (b_nr, b_nc, b_nz); + retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx(0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == SparseType::Permuted_Lower) { - OCTAVE_LOCAL_BUFFER (Complex, cwork, nr); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 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++) + for (octave_idx_type i = 0; i < nm; i++) cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) cwork[perm[b.ridx(i)]] = b.data(i); - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (cwork[k] != 0.) { @@ -3231,7 +3275,7 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (cwork[i] != 0.) new_nnz++; @@ -3243,7 +3287,7 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) + for (octave_idx_type i = 0; i < nc; i++) if (cwork[i] != 0.) { retval.xridx(ii) = i; @@ -3255,15 +3299,15 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work, nr); - for (octave_idx_type i = 0; i < nr; i++) + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; - for (octave_idx_type k = 0; k < nr; k++) + for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { @@ -3291,7 +3335,7 @@ } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { atmp += fabs(work[i]); work[i] = 0.; @@ -3302,18 +3346,18 @@ } else { - OCTAVE_LOCAL_BUFFER (Complex, work, nr); + OCTAVE_LOCAL_BUFFER (Complex, cwork, nm); 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 = 0; i < nm; i++) + cwork[i] = 0.; for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++) - work[b.ridx(i)] = b.data(i); - - for (octave_idx_type k = 0; k < nr; k++) + cwork[b.ridx(i)] = b.data(i); + + for (octave_idx_type k = 0; k < nc; k++) { - if (work[k] != 0.) + if (cwork[k] != 0.) { if (ridx(cidx(k)) != k) { @@ -3321,12 +3365,12 @@ goto triangular_error; } - Complex tmp = work[k] / data(cidx(k)); - work[k] = tmp; + Complex tmp = cwork[k] / data(cidx(k)); + cwork[k] = tmp; for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - work[iidx] = work[iidx] - tmp * data(i); + cwork[iidx] = cwork[iidx] - tmp * data(i); } } } @@ -3334,8 +3378,8 @@ // Count non-zeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + for (octave_idx_type i = 0; i < nc; i++) + if (cwork[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) @@ -3346,11 +3390,11 @@ x_nz = sz; } - for (octave_idx_type i = 0; i < nr; i++) - if (work[i] != 0.) + for (octave_idx_type i = 0; i < nc; i++) + if (cwork[i] != 0.) { retval.xridx(ii) = i; - retval.xdata(ii++) = work[i]; + retval.xdata(ii++) = cwork[i]; } retval.xcidx(j+1) = ii; } @@ -3358,33 +3402,33 @@ retval.maybe_compress (); // Calculation of 1-norm of inv(*this) - OCTAVE_LOCAL_BUFFER (double, work2, nr); - for (octave_idx_type i = 0; i < nr; i++) - work2[i] = 0.; + OCTAVE_LOCAL_BUFFER (double, work, nm); + for (octave_idx_type i = 0; i < nm; i++) + work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { - work2[j] = 1.; - - for (octave_idx_type k = j; k < nr; k++) + work[j] = 1.; + + for (octave_idx_type k = j; k < nc; k++) { - if (work2[k] != 0.) + if (work[k] != 0.) { - double tmp = work2[k] / data(cidx(k)); - work2[k] = tmp; + double tmp = work[k] / data(cidx(k)); + work[k] = tmp; for (octave_idx_type i = cidx(k)+1; i < cidx(k+1); i++) { octave_idx_type iidx = ridx(i); - work2[iidx] = work2[iidx] - tmp * data(i); + work[iidx] = work[iidx] - tmp * data(i); } } } double atmp = 0; - for (octave_idx_type i = j; i < nr; i++) + for (octave_idx_type i = j; i < nc; i++) { - atmp += fabs(work2[i]); - work2[i] = 0.; + atmp += fabs(work[i]); + work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp;
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,16 @@ +2006-02-20 David Bateman <dbateman@free.fr> + + * pt-mat.cc (class tm_row_const): Add any_sparse bool variable. + (tm_row_const::tm_row_const_rep::do_init_element): Initialize + any_sparse variable if any matrice is sparse. + (class tm_const): Add any_sparse bool variable. + (tm_const::init): Initialize any_sparse variable. + (tree_matrix::rvalue): If any matrix is sparse use sparse matrix + as initial matrix for concatenation + * DLD-FUNCTIONS/matrix_type.cc: Add tests for new rectangular + diagonal, permuted diagonal, triangular and permuted triangular + matrices + 2006-02-20 John W. Eaton <jwe@octave.org> * toplev.cc (__builtin_delete, __builtin_new): Use std::cerr for
--- a/src/DLD-FUNCTIONS/matrix_type.cc +++ b/src/DLD-FUNCTIONS/matrix_type.cc @@ -316,7 +316,19 @@ %! a=[speye(10,10),[sparse(9,1);1];-1,sparse(1,9),1]; %! assert(matrix_type(a),"Full"); %! assert(matrix_type(a'*a),"Positive Definite"); -%!assert(matrix_type(speye(10,11)),"Rectangular"); +%!assert(matrix_type(speye(10,11)),"Diagonal"); +%!assert(matrix_type(speye(10,11)([2:10,1],:)),"Permuted Diagonal"); +%!assert(matrix_type(speye(11,10)),"Diagonal"); +%!assert(matrix_type(speye(11,10)([2:11,1],:)),"Permuted Diagonal"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[[1,1];sparse(9,2);[1,1]]]),"Upper"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[[1,1];sparse(9,2);[1,1]]](:,[2,1,3:12])),"Permuted Upper"); +%!assert(matrix_type([speye(11,9),[1;sparse(8,1);1;0]]),"Upper"); +%!assert(matrix_type([speye(11,9),[1;sparse(8,1);1;0]](:,[2,1,3:10])),"Permuted Upper"); +%!assert(matrix_type([speye(10,10),sparse(10,1);[1;1],sparse(2,9),[1;1]]),"Lower"); +%!assert(matrix_type([speye(10,10),sparse(10,1);[1;1],sparse(2,9),[1;1]]([2,1,3:12],:)),"Permuted Lower"); +%!assert(matrix_type([speye(9,11);[1,sparse(1,8),1,0]]),"Lower"); +%!assert(matrix_type([speye(9,11);[1,sparse(1,8),1,0]]([2,1,3:10],:)),"Permuted Lower"); +%!assert(matrix_type(spdiags(randn(10,4),[-2:1],10,9)),"Rectangular") %!assert(matrix_type(1i*speye(10,10)),"Diagonal"); %!assert(matrix_type(1i*speye(10,10)([2:10,1],:)),"Permuted Diagonal"); @@ -342,7 +354,19 @@ %! a=[speye(10,10),[sparse(9,1);1i];-1,sparse(1,9),1]; %! assert(matrix_type(a),"Full"); %! assert(matrix_type(a'*a),"Positive Definite"); -%!assert(matrix_type(speye(10,11)),"Rectangular"); +%!assert(matrix_type(1i*speye(10,11)),"Diagonal"); +%!assert(matrix_type(1i*speye(10,11)([2:10,1],:)),"Permuted Diagonal"); +%!assert(matrix_type(1i*speye(11,10)),"Diagonal"); +%!assert(matrix_type(1i*speye(11,10)([2:11,1],:)),"Permuted Diagonal"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[[1i,1i];sparse(9,2);[1i,1i]]]),"Upper"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[[1i,1i];sparse(9,2);[1i,1i]]](:,[2,1,3:12])),"Permuted Upper"); +%!assert(matrix_type([speye(11,9),[1i;sparse(8,1);1i;0]]),"Upper"); +%!assert(matrix_type([speye(11,9),[1i;sparse(8,1);1i;0]](:,[2,1,3:10])),"Permuted Upper"); +%!assert(matrix_type([speye(10,10),sparse(10,1);[1i;1i],sparse(2,9),[1i;1i]]),"Lower"); +%!assert(matrix_type([speye(10,10),sparse(10,1);[1i;1i],sparse(2,9),[1i;1i]]([2,1,3:12],:)),"Permuted Lower"); +%!assert(matrix_type([speye(9,11);[1i,sparse(1,8),1i,0]]),"Lower"); +%!assert(matrix_type([speye(9,11);[1i,sparse(1,8),1i,0]]([2,1,3:10],:)),"Permuted Lower"); +%!assert(matrix_type(1i*spdiags(randn(10,4),[-2:1],10,9)),"Rectangular") %!test %! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular");
--- a/src/pt-mat.cc +++ b/src/pt-mat.cc @@ -41,6 +41,9 @@ #include "ov.h" #include "variables.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + // If TRUE, print a warning message for empty elements in a matrix list. static bool Vwarn_empty_list_elements; @@ -70,15 +73,15 @@ : count (1), dv (0, 0), all_str (false), all_sq_str (false), all_dq_str (false), some_str (false), all_real (false), all_cmplx (false), - all_mt (true), class_nm (octave_base_value::static_class_name ()), - ok (false) + all_mt (true), any_sparse (false), + class_nm (octave_base_value::static_class_name ()), ok (false) { } tm_row_const_rep (const tree_argument_list& row) : count (1), dv (0, 0), all_str (false), all_sq_str (false), some_str (false), all_real (false), all_cmplx (false), - all_mt (true), class_nm (octave_base_value::static_class_name ()), - ok (false) + all_mt (true), any_sparse (false), + class_nm (octave_base_value::static_class_name ()), ok (false) { init (row); } ~tm_row_const_rep (void) { } @@ -94,6 +97,7 @@ bool all_real; bool all_cmplx; bool all_mt; + bool any_sparse; std::string class_nm; @@ -167,6 +171,7 @@ bool all_real_p (void) const { return rep->all_real; } bool all_complex_p (void) const { return rep->all_cmplx; } bool all_empty_p (void) const { return rep->all_mt; } + bool any_sparse_p (void) const { return rep->any_sparse; } std::string class_name (void) const { return rep->class_nm; } @@ -341,6 +346,9 @@ if (all_cmplx && ! (val.is_complex_type () || val.is_real_type ())) all_cmplx = false; + if (!any_sparse && val.class_name() == "sparse") + any_sparse = true; + return true; } @@ -352,6 +360,7 @@ all_dq_str = true; all_real = true; all_cmplx = true; + any_sparse = false; bool first_elem = true; @@ -432,8 +441,8 @@ tm_const (const tree_matrix& tm) : dv (0, 0), all_str (false), all_sq_str (false), all_dq_str (false), some_str (false), all_real (false), all_cmplx (false), - all_mt (true), class_nm (octave_base_value::static_class_name ()), - ok (false) + all_mt (true), any_sparse (false), + class_nm (octave_base_value::static_class_name ()), ok (false) { init (tm); } ~tm_const (void) { } @@ -450,6 +459,7 @@ bool all_real_p (void) const { return all_real; } bool all_complex_p (void) const { return all_cmplx; } bool all_empty_p (void) const { return all_mt; } + bool any_sparse_p (void) const { return any_sparse; } std::string class_name (void) const { return class_nm; } @@ -466,6 +476,7 @@ bool all_real; bool all_cmplx; bool all_mt; + bool any_sparse; std::string class_nm; @@ -488,6 +499,7 @@ all_dq_str = true; all_real = true; all_cmplx = true; + any_sparse = false; bool first_elem = true; @@ -527,6 +539,9 @@ if (all_mt && ! tmp.all_empty_p ()) all_mt = false; + if (!any_sparse && tmp.any_sparse_p ()) + any_sparse = true; + append (tmp); } else @@ -754,6 +769,7 @@ bool all_empty_p = false; bool all_real_p = false; bool all_complex_p = false; + bool any_sparse_p = false; bool frc_str_conv = false; tm_const tmp (*this); @@ -767,6 +783,7 @@ all_empty_p = tmp.all_empty_p (); all_real_p = tmp.all_real_p (); all_complex_p = tmp.all_complex_p (); + any_sparse_p = tmp.any_sparse_p (); frc_str_conv = tmp.some_strings_p (); // Try to speed up the common cases. @@ -836,30 +853,42 @@ // Find the first non-empty object - for (tm_const::iterator p = tmp.begin (); p != tmp.end (); p++) + if (any_sparse_p) { - OCTAVE_QUIT; - - tm_row_const row = *p; - - for (tm_row_const::iterator q = row.begin (); - q != row.end (); q++) + // Start with sparse matrix to avoid issues memory issues + // with things like [ones(1,4),sprandn(1e8,4,1e-4)] + if (all_real_p) + ctmp = octave_sparse_matrix ().resize (dv); + else + ctmp = octave_sparse_complex_matrix ().resize (dv); + } + else + { + for (tm_const::iterator p = tmp.begin (); p != tmp.end (); p++) { OCTAVE_QUIT; - ctmp = *q; + tm_row_const row = *p; - if (! ctmp.all_zero_dims ()) - goto found_non_empty; - } - } + for (tm_row_const::iterator q = row.begin (); + q != row.end (); q++) + { + OCTAVE_QUIT; + + ctmp = *q; - ctmp = (*(tmp.begin() -> begin())); + if (! ctmp.all_zero_dims ()) + goto found_non_empty; + } + } - found_non_empty: + ctmp = (*(tmp.begin() -> begin())); + + found_non_empty: - if (! all_empty_p) - ctmp = ctmp.resize (dim_vector (0,0)).resize (dv); + if (! all_empty_p) + ctmp = ctmp.resize (dim_vector (0,0)).resize (dv); + } if (! error_state) {
--- a/test/ChangeLog +++ b/test/ChangeLog @@ -1,6 +1,12 @@ +2006-02-20 David Bateman <dbateman@free.fr> + + * build_spase_tests.sh: Add tests for ldiv tests for rectangular + diagonal, permuted diagonal, triangular and permuted triangular + matrices. + 2006-02-09 David Bateman <dbateman@free.fr> - * test/build_sparse_tests.sh: Add tests for sparse QR solvers. + * build_sparse_tests.sh: Add tests for sparse QR solvers. 2006-01-21 David Bateman <dbateman@free.fr>
--- a/test/build_sparse_tests.sh +++ b/test/build_sparse_tests.sh @@ -858,10 +858,6 @@ # ============================================================= # Specific solver tests for matrices that will test all of the solver # code. Uses alpha and beta -# -# Note these tests can still fail with a singular matrix, as sprandn -# is used and QR solvers not implemented. However, that should happen -# rarely gen_solver_tests() { if $preset; then @@ -929,10 +925,90 @@ %! assert (sparse(a * x), b, feps); %!test %! a = alpha*sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10); -%! ## Test this by forcing matrix_type +%! ## Test this by forcing matrix_type, as can't get a certain +%! ## result for over-determined systems. %! a = alpha*sprandn(10,10,0.2)+speye(10,10); matrix_type(a, "Singular"); %! f(a,[10,2],1e-10); +%% Rectanguar solver tests that don't use QR + +%!test +%! ds = alpha * spdiags([1:11]',0,10,11); +%! df = full(ds); +%! xf = beta * ones(10,1); +%! xs = speye(10,10); +%!assert(ds\xf,df\xf,100*eps) +%!assert(ds\xs,sparse(df\xs,true),100*eps) +%!test +%! pds = ds([2,1,3:10],:); +%! pdf = full(pds); +%!assert(pds\xf,pdf\xf,100*eps) +%!assert(pds\xs,sparse(pdf\xs,true),100*eps) +%!test +%! ds = alpha * spdiags([1:11]',0,11,10); +%! df = full(ds); +%! xf = beta * ones(11,1); +%! xs = speye(11,11); +%!assert(ds\xf,df\xf,100*eps) +%!assert(ds\xs,sparse(df\xs,true),100*eps) +%!test +%! pds = ds([2,1,3:11],:); +%! pdf = full(pds); +%!assert(pds\xf,pdf\xf,100*eps) +%!assert(pds\xs,sparse(pdf\xs,true),100*eps) +%!test +%! us = alpha*[[speye(10,10);sparse(1,10)],[[1,1];sparse(9,2);[1,1]]]; +%!assert(us*(us\xf),xf,100*eps) +%!assert(us*(us\xs),xs,100*eps) +%!test +%! pus = us(:,[2,1,3:12]); +%!assert(pus*(pus\xf),xf,100*eps) +%!assert(pus*(pus\xs),xs,100*eps) +%!test +%! us = alpha*[speye(11,9),[1;sparse(8,1);1;0]]; +%!test +%! [c,r] = spqr (us, xf); +%! assert(us\xf,r\c,100*eps) +%!test +%! [c,r] = spqr (us, xs); +%! assert(us\xs,r\c,100*eps) +%!test +%! pus = us(:,[1:8,10,9]); +%!test +%! [c,r] = spqr (pus, xf); +%! assert(pus\xf,r\c,100*eps) +%!test +%! [c,r] = spqr (pus, xs); +%! assert(pus\xs,r\c,100*eps) +%!test +%! ls = alpha*[speye(9,11);[1,sparse(1,8),1,0]]; +%! xf = beta * ones(10,1); +%! xs = speye(10,10); +%!assert(ls*(ls\xf),xf,100*eps) +%!assert(ls*(ls\xs),xs,100*eps) +%!test +%! pls = ls([1:8,10,9],:); +%!assert(pls*(pls\xf),xf,100*eps) +%!assert(pls*(pls\xs),xs,100*eps) +%!test +%! ls = alpha*[speye(10,10),sparse(10,1);[1;1],sparse(2,9),[1;1]]; +%! xf = beta * ones(12,1); +%! xs = speye(12,12); +%!test +%! [c,r] = spqr (ls, xf); +%! assert(ls\xf,r\c,100*eps) +%!test +%! [c,r] = spqr (ls, xs); +%! assert(ls\xs,r\c,100*eps) +%!test +%! pls = ls(:,[1:8,10,9]); +%!test +%! [c,r] = spqr (pls, xf); +%! assert(pls\xf,r\c,100*eps) +%!test +%! [c,r] = spqr (pls, xs); +%! assert(pls\xs,r\c,100*eps) + EOF }