Mercurial > hg > octave-lyh
changeset 7455:fe332ce262b5
eliminate spkron.cc; dispatch in kron
author | David Bateman <dbateman@free.fr> |
---|---|
date | Thu, 07 Feb 2008 06:00:37 -0500 |
parents | d03ae525c894 |
children | fee437704c67 |
files | doc/ChangeLog doc/interpreter/sparse.txi scripts/deprecated/spkron.m src/ChangeLog src/DLD-FUNCTIONS/kron.cc src/DLD-FUNCTIONS/spkron.cc src/Makefile.in |
diffstat | 7 files changed, 132 insertions(+), 179 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,7 @@ +2008-02-07 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Remove references to spkron. + 2008-02-06 John W. Eaton <jwe@octave.org> * interpreter/Makefile.in (stmp-html, HTML): New targets.
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -489,7 +489,7 @@ @item Linear algebra: @dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv}, - @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron}, + @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, @dfn{condest}, @dfn{sprank} @c @dfn{spaugment} @@ -848,8 +848,6 @@ @DOCSTRING(spinv) -@DOCSTRING(spkron) - @DOCSTRING(splchol) @DOCSTRING(splu)
new file mode 100644 --- /dev/null +++ b/scripts/deprecated/spkron.m @@ -0,0 +1,30 @@ +## Copyright (C) 2008 John W. Eaton +## +## 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 +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {} spkron (@var{a}, @var{b}) +## This function has been deprecated. Use @code{kron} instead. +## @end deftypefn + +## Author: jwe + +function retval = spkron (varargin) + + retval = kron (varargin{:}); + +endfunction
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,10 @@ +2008-02-07 David Bateman <dbateman@free.fr> + + * Makefile.in (DLD_XSRC): Delete spkron.cc. + * DLD-FUNCTIONS/spkron.cc: Delete. + * DLD-FUNCTIONS/kron.cc: Include here and dispatch to the sparse + version if either argument is sparse. + 2008-02-06 John W. Eaton <jwe@octave.org> * Makefile.in (parse.cc): Also use --defines option for bison.
--- a/src/DLD-FUNCTIONS/kron.cc +++ b/src/DLD-FUNCTIONS/kron.cc @@ -69,6 +69,51 @@ template void kron (const Array2<Complex>&, const Array2<Complex>&, Array2<Complex>&); + +#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) +extern void +kron (const Sparse<double>&, const Sparse<double>&, Sparse<double>&); + +extern void +kron (const Sparse<Complex>&, const Sparse<Complex>&, Sparse<Complex>&); +#endif + +template <class T> +void +kron (const Sparse<T>& A, const Sparse<T>& B, Sparse<T>& C) +{ + octave_idx_type idx = 0; + C = Sparse<T> (A.rows () * B.rows (), A.columns () * B.columns (), + A.nzmax () * B.nzmax ()); + + C.cidx (0) = 0; + + for (octave_idx_type Aj = 0; Aj < A.columns (); Aj++) + for (octave_idx_type Bj = 0; Bj < B.columns (); Bj++) + { + for (octave_idx_type Ai = A.cidx (Aj); Ai < A.cidx (Aj+1); Ai++) + { + octave_idx_type Ci = A.ridx(Ai) * B.rows (); + const T v = A.data (Ai); + + for (octave_idx_type Bi = B.cidx (Bj); Bi < B.cidx (Bj+1); Bi++) + { + OCTAVE_QUIT; + C.data (idx) = v * B.data (Bi); + C.ridx (idx++) = Ci + B.ridx (Bi); + } + } + C.cidx (Aj * B.columns () + Bj + 1) = idx; + } +} + +template void +kron (const Sparse<double>&, const Sparse<double>&, Sparse<double>&); + +template void +kron (const Sparse<Complex>&, const Sparse<Complex>&, Sparse<Complex>&); + + DEFUN_DLD (kron, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} kron (@var{a}, @var{b})\n\ Form the kronecker product of two matrices, defined block by block as\n\ @@ -97,28 +142,58 @@ { print_usage (); } - else if (args(0).is_complex_type () || args(1).is_complex_type ()) + else if (args(0).is_sparse_type () || args(1).is_sparse_type ()) { - ComplexMatrix a (args(0).complex_matrix_value()); - ComplexMatrix b (args(1).complex_matrix_value()); + if (args(0).is_complex_type () || args(1).is_complex_type ()) + { + SparseComplexMatrix a (args(0).sparse_complex_matrix_value()); + SparseComplexMatrix b (args(1).sparse_complex_matrix_value()); - if (! error_state) + if (! error_state) + { + SparseComplexMatrix c; + kron (a, b, c); + retval(0) = c; + } + } + else { - ComplexMatrix c; - kron (a, b, c); - retval(0) = c; + SparseMatrix a (args(0).sparse_matrix_value ()); + SparseMatrix b (args(1).sparse_matrix_value ()); + + if (! error_state) + { + SparseMatrix c; + kron (a, b, c); + retval (0) = c; + } } } - else + else { - Matrix a (args(0).matrix_value ()); - Matrix b (args(1).matrix_value ()); + if (args(0).is_complex_type () || args(1).is_complex_type ()) + { + ComplexMatrix a (args(0).complex_matrix_value()); + ComplexMatrix b (args(1).complex_matrix_value()); - if (! error_state) + if (! error_state) + { + ComplexMatrix c; + kron (a, b, c); + retval(0) = c; + } + } + else { - Matrix c; - kron (a, b, c); - retval (0) = c; + Matrix a (args(0).matrix_value ()); + Matrix b (args(1).matrix_value ()); + + if (! error_state) + { + Matrix c; + kron (a, b, c); + retval (0) = c; + } } }
--- a/src/DLD-FUNCTIONS/spkron.cc +++ b/src/DLD-FUNCTIONS/spkron.cc @@ -1,161 +0,0 @@ -/* - -Copyright (C) 2002, 2005, 2006, 2007 John W. Eaton -Copyright (C) 2005 David Bateman - -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 -<http://www.gnu.org/licenses/>. - -*/ - -// Author: Paul Kienzle <pkienzle@users.sf.net> - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "dSparse.h" -#include "CSparse.h" -#include "quit.h" - -#include "defun-dld.h" -#include "error.h" -#include "oct-obj.h" - -#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) -extern void -kron (const Sparse<double>&, const Sparse<double>&, Sparse<double>&); - -extern void -kron (const Sparse<Complex>&, const Sparse<Complex>&, Sparse<Complex>&); -#endif - -template <class T> -void -kron (const Sparse<T>& A, const Sparse<T>& B, Sparse<T>& C) -{ - octave_idx_type idx = 0; - C = Sparse<T> (A.rows () * B.rows (), A.columns () * B.columns (), - A.nzmax () * B.nzmax ()); - - C.cidx (0) = 0; - - for (octave_idx_type Aj = 0; Aj < A.columns (); Aj++) - for (octave_idx_type Bj = 0; Bj < B.columns (); Bj++) - { - for (octave_idx_type Ai = A.cidx (Aj); Ai < A.cidx (Aj+1); Ai++) - { - octave_idx_type Ci = A.ridx(Ai) * B.rows (); - const T v = A.data (Ai); - - for (octave_idx_type Bi = B.cidx (Bj); Bi < B.cidx (Bj+1); Bi++) - { - OCTAVE_QUIT; - C.data (idx) = v * B.data (Bi); - C.ridx (idx++) = Ci + B.ridx (Bi); - } - } - C.cidx (Aj * B.columns () + Bj + 1) = idx; - } -} - -template void -kron (const Sparse<double>&, const Sparse<double>&, Sparse<double>&); - -template void -kron (const Sparse<Complex>&, const Sparse<Complex>&, Sparse<Complex>&); - -// PKG_ADD: dispatch ("kron", "spkron", "sparse matrix"); -// PKG_ADD: dispatch ("kron", "spkron", "sparse complex matrix"); -// PKG_ADD: dispatch ("kron", "spkron", "sparse bool matrix"); -DEFUN_DLD (spkron, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spkron (@var{a}, @var{b})\n\ -Form the kronecker product of two sparse matrices. This is defined\n\ -block by block as\n\ -\n\ -@example\n\ -x = [a(i, j) b]\n\ -@end example\n\ -\n\ -For example,\n\ -\n\ -@example\n\ -@group\n\ -kron(speye(3),spdiag([1,2,3]))\n\ -@result{}\n\ -Compressed Column Sparse (rows = 9, cols = 9, nnz = 9)\n\ -\n\ - (1, 1) -> 1\n\ - (2, 2) -> 2\n\ - (3, 3) -> 3\n\ - (4, 4) -> 1\n\ - (5, 5) -> 2\n\ - (6, 6) -> 3\n\ - (7, 7) -> 1\n\ - (8, 8) -> 2\n\ - (9, 9) -> 3\n\ -@end group\n\ -@end example\n\ -@end deftypefn") -{ - octave_value_list retval; - - int nargin = args.length (); - - if (nargin != 2 || nargout > 1) - { - print_usage (); - } - else if (args(0).is_complex_type () || args(1).is_complex_type ()) - { - SparseComplexMatrix a (args(0).sparse_complex_matrix_value()); - SparseComplexMatrix b (args(1).sparse_complex_matrix_value()); - - if (! error_state) - { - SparseComplexMatrix c; - kron (a, b, c); - retval(0) = c; - } - } - else - { - SparseMatrix a (args(0).sparse_matrix_value ()); - SparseMatrix b (args(1).sparse_matrix_value ()); - - if (! error_state) - { - SparseMatrix c; - kron (a, b, c); - retval (0) = c; - } - } - - return retval; -} - -/* - -%!assert(spkron(spdiag([1,2,3]),spdiag([1,2,3])),sparse(kron(diag([1,2,3]),diag([1,2,3])))) -%!assert(spkron(spdiag([1i,2,3]),spdiag([1i,2,3])),sparse(kron(diag([1i,2,3]),diag([1i,2,3])))) - -*/ - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -70,7 +70,7 @@ givens.cc hess.cc inv.cc kron.cc lsode.cc \ lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \ - spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \ + spchol.cc spdet.cc spfind.cc splu.cc spparms.cc spqr.cc \ sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \ urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \ __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \