Mercurial > octave
view libinterp/corefcn/sparse-xdiv.cc @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 28 Dec 2021 18:22:40 -0500 |
parents | a61e1a0f6024 |
children | e88a07dec498 |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1998-2022 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // This file is part of Octave. // // Octave is free software: you can redistribute it and/or modify it // under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // Octave is distributed in the hope that it will be useful, but // WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with Octave; see the file COPYING. If not, see // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <cassert> #include "Array-util.h" #include "lo-array-errwarn.h" #include "oct-cmplx.h" #include "quit.h" #include "error.h" #include "lo-ieee.h" #include "dSparse.h" #include "dDiagMatrix.h" #include "CSparse.h" #include "CDiagMatrix.h" #include "oct-spparms.h" #include "sparse-xdiv.h" OCTAVE_NAMESPACE_BEGIN static void solve_singularity_warning (double rcond) { octave::warn_singular_matrix (rcond); } template <typename T1, typename T2> bool mx_leftdiv_conform (const T1& a, const T2& b) { octave_idx_type a_nr = a.rows (); octave_idx_type b_nr = b.rows (); if (a_nr != b_nr) { octave_idx_type a_nc = a.cols (); octave_idx_type b_nc = b.cols (); octave::err_nonconformant (R"(operator \)", a_nr, a_nc, b_nr, b_nc); } return true; } #define INSTANTIATE_MX_LEFTDIV_CONFORM(T1, T2) \ template bool mx_leftdiv_conform (const T1&, const T2&) INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, SparseComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, Matrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseMatrix, ComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, Matrix); INSTANTIATE_MX_LEFTDIV_CONFORM (SparseComplexMatrix, ComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (DiagMatrix, SparseComplexMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseMatrix); INSTANTIATE_MX_LEFTDIV_CONFORM (ComplexDiagMatrix, SparseComplexMatrix); template <typename T1, typename T2> bool mx_div_conform (const T1& a, const T2& b) { octave_idx_type a_nc = a.cols (); octave_idx_type b_nc = b.cols (); if (a_nc != b_nc) { octave_idx_type a_nr = a.rows (); octave_idx_type b_nr = b.rows (); octave::err_nonconformant ("operator /", a_nr, a_nc, b_nr, b_nc); } return true; } #define INSTANTIATE_MX_DIV_CONFORM(T1, T2) \ template bool mx_div_conform (const T1&, const T2&) INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, SparseComplexMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, SparseComplexMatrix); INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseMatrix); INSTANTIATE_MX_DIV_CONFORM (Matrix, SparseComplexMatrix); INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseMatrix); INSTANTIATE_MX_DIV_CONFORM (ComplexMatrix, SparseComplexMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, DiagMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseMatrix, ComplexDiagMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, DiagMatrix); INSTANTIATE_MX_DIV_CONFORM (SparseComplexMatrix, ComplexDiagMatrix); // Right division functions. X / Y = X * inv (Y) = (inv (Y') * X')' // // Y / X: m cm sm scm // +-- +---+----+----+----+ // sparse matrix | 1 | 3 | 5 | 7 | // +---+----+----+----+ // sparse complex_matrix | 2 | 4 | 6 | 8 | // +---+----+----+----+ // diagonal matrix | 9 | 11 | // +----+----+ // complex diag. matrix | 10 | 12 | // +----+----+ // -*- 1 -*- Matrix xdiv (const Matrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return Matrix (); Matrix atmp = a.transpose (); SparseMatrix btmp = b.transpose (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; Matrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.transpose (); } // -*- 2 -*- ComplexMatrix xdiv (const Matrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); Matrix atmp = a.transpose (); SparseComplexMatrix btmp = b.hermitian (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; ComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } // -*- 3 -*- ComplexMatrix xdiv (const ComplexMatrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); ComplexMatrix atmp = a.hermitian (); SparseMatrix btmp = b.transpose (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; ComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } // -*- 4 -*- ComplexMatrix xdiv (const ComplexMatrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return ComplexMatrix (); ComplexMatrix atmp = a.hermitian (); SparseComplexMatrix btmp = b.hermitian (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; ComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } // -*- 5 -*- SparseMatrix xdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return SparseMatrix (); SparseMatrix atmp = a.transpose (); SparseMatrix btmp = b.transpose (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; SparseMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.transpose (); } // -*- 6 -*- SparseComplexMatrix xdiv (const SparseMatrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseMatrix atmp = a.transpose (); SparseComplexMatrix btmp = b.hermitian (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; SparseComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } // -*- 7 -*- SparseComplexMatrix xdiv (const SparseComplexMatrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseComplexMatrix atmp = a.hermitian (); SparseMatrix btmp = b.transpose (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; SparseComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } // -*- 8 -*- SparseComplexMatrix xdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_div_conform (a, b)) return SparseComplexMatrix (); SparseComplexMatrix atmp = a.hermitian (); SparseComplexMatrix btmp = b.hermitian (); MatrixType btyp = typ.transpose (); octave_idx_type info; double rcond = 0.0; SparseComplexMatrix result = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning); typ = btyp.transpose (); return result.hermitian (); } template <typename RT, typename SM, typename DM> RT do_rightdiv_sm_dm (const SM& a, const DM& d) { const octave_idx_type d_nr = d.rows (); const octave_idx_type a_nr = a.rows (); const octave_idx_type a_nc = a.cols (); using std::min; const octave_idx_type nc = min (d_nr, a_nc); if (! mx_div_conform (a, d)) return RT (); const octave_idx_type nz = a.nnz (); RT r (a_nr, nc, nz); typedef typename DM::element_type DM_elt_type; const DM_elt_type zero = DM_elt_type (); octave_idx_type k_result = 0; for (octave_idx_type j = 0; j < nc; ++j) { octave_quit (); const DM_elt_type s = d.dgelem (j); const octave_idx_type colend = a.cidx (j+1); r.xcidx (j) = k_result; if (s != zero) for (octave_idx_type k = a.cidx (j); k < colend; ++k) { r.xdata (k_result) = a.data (k) / s; r.xridx (k_result) = a.ridx (k); ++k_result; } } r.xcidx (nc) = k_result; r.maybe_compress (true); return r; } // -*- 9 -*- SparseMatrix xdiv (const SparseMatrix& a, const DiagMatrix& b, MatrixType&) { return do_rightdiv_sm_dm<SparseMatrix> (a, b); } // -*- 10 -*- SparseComplexMatrix xdiv (const SparseMatrix& a, const ComplexDiagMatrix& b, MatrixType&) { return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b); } // -*- 11 -*- SparseComplexMatrix xdiv (const SparseComplexMatrix& a, const DiagMatrix& b, MatrixType&) { return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b); } // -*- 12 -*- SparseComplexMatrix xdiv (const SparseComplexMatrix& a, const ComplexDiagMatrix& b, MatrixType&) { return do_rightdiv_sm_dm<SparseComplexMatrix> (a, b); } // Funny element by element division operations. // // op2 \ op1: s cs // +-- +---+----+ // matrix | 1 | 3 | // +---+----+ // complex_matrix | 2 | 4 | // +---+----+ Matrix elem_xdiv (double a, const SparseMatrix& b) { octave_idx_type nr = b.rows (); octave_idx_type nc = b.cols (); Matrix result; if (a == 0.) result = Matrix (nr, nc, octave::numeric_limits<double>::NaN ()); else if (a > 0.) result = Matrix (nr, nc, octave::numeric_limits<double>::Inf ()); else result = Matrix (nr, nc, -octave::numeric_limits<double>::Inf ()); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { octave_quit (); result.elem (b.ridx (i), j) = a / b.data (i); } return result; } ComplexMatrix elem_xdiv (double a, const SparseComplexMatrix& b) { octave_idx_type nr = b.rows (); octave_idx_type nc = b.cols (); ComplexMatrix result (nr, nc, Complex (octave::numeric_limits<double>::NaN (), octave::numeric_limits<double>::NaN ())); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { octave_quit (); result.elem (b.ridx (i), j) = a / b.data (i); } return result; } ComplexMatrix elem_xdiv (const Complex& a, const SparseMatrix& b) { octave_idx_type nr = b.rows (); octave_idx_type nc = b.cols (); ComplexMatrix result (nr, nc, (a / 0.0)); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { octave_quit (); result.elem (b.ridx (i), j) = a / b.data (i); } return result; } ComplexMatrix elem_xdiv (const Complex& a, const SparseComplexMatrix& b) { octave_idx_type nr = b.rows (); octave_idx_type nc = b.cols (); ComplexMatrix result (nr, nc, (a / 0.0)); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { octave_quit (); result.elem (b.ridx (i), j) = a / b.data (i); } return result; } // Left division functions. X \ Y = inv (X) * Y // // Y \ X : sm scm dm dcm // +-- +---+----+ // matrix | 1 | 5 | // +---+----+ // complex_matrix | 2 | 6 | // +---+----+----+----+ // sparse matrix | 3 | 7 | 9 | 11 | // +---+----+----+----+ // sparse complex_matrix | 4 | 8 | 10 | 12 | // +---+----+----+----+ // -*- 1 -*- Matrix xleftdiv (const SparseMatrix& a, const Matrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return Matrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 2 -*- ComplexMatrix xleftdiv (const SparseMatrix& a, const ComplexMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 3 -*- SparseMatrix xleftdiv (const SparseMatrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return SparseMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 4 -*- SparseComplexMatrix xleftdiv (const SparseMatrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 5 -*- ComplexMatrix xleftdiv (const SparseComplexMatrix& a, const Matrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 6 -*- ComplexMatrix xleftdiv (const SparseComplexMatrix& a, const ComplexMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return ComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 7 -*- SparseComplexMatrix xleftdiv (const SparseComplexMatrix& a, const SparseMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } // -*- 8 -*- SparseComplexMatrix xleftdiv (const SparseComplexMatrix& a, const SparseComplexMatrix& b, MatrixType& typ) { if (! mx_leftdiv_conform (a, b)) return SparseComplexMatrix (); octave_idx_type info; double rcond = 0.0; return a.solve (typ, b, info, rcond, solve_singularity_warning); } template <typename RT, typename DM, typename SM> RT do_leftdiv_dm_sm (const DM& d, const SM& a) { const octave_idx_type a_nr = a.rows (); const octave_idx_type a_nc = a.cols (); const octave_idx_type d_nc = d.cols (); using std::min; const octave_idx_type nr = min (d_nc, a_nr); if (! mx_leftdiv_conform (d, a)) return RT (); const octave_idx_type nz = a.nnz (); RT r (nr, a_nc, nz); typedef typename DM::element_type DM_elt_type; const DM_elt_type zero = DM_elt_type (); octave_idx_type k_result = 0; for (octave_idx_type j = 0; j < a_nc; ++j) { octave_quit (); const octave_idx_type colend = a.cidx (j+1); r.xcidx (j) = k_result; for (octave_idx_type k = a.cidx (j); k < colend; ++k) { const octave_idx_type i = a.ridx (k); if (i < nr) { const DM_elt_type s = d.dgelem (i); if (s != zero) { r.xdata (k_result) = a.data (k) / s; r.xridx (k_result) = i; ++k_result; } } } } r.xcidx (a_nc) = k_result; r.maybe_compress (true); return r; } // -*- 9 -*- SparseMatrix xleftdiv (const DiagMatrix& d, const SparseMatrix& a, MatrixType&) { return do_leftdiv_dm_sm<SparseMatrix> (d, a); } // -*- 10 -*- SparseComplexMatrix xleftdiv (const DiagMatrix& d, const SparseComplexMatrix& a, MatrixType&) { return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a); } // -*- 11 -*- SparseComplexMatrix xleftdiv (const ComplexDiagMatrix& d, const SparseMatrix& a, MatrixType&) { return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a); } // -*- 12 -*- SparseComplexMatrix xleftdiv (const ComplexDiagMatrix& d, const SparseComplexMatrix& a, MatrixType&) { return do_leftdiv_dm_sm<SparseComplexMatrix> (d, a); } OCTAVE_NAMESPACE_END