# HG changeset patch # User dbateman # Date 1114780236 0 # Node ID e5a68648db9ce229b9c77ae4ebcdda605d7cd53a # Parent 22994a5730f9d67861a6e067d3e130997c11c548 [project @ 2005-04-29 13:10:36 by dbateman] diff -r 22994a5730f9 -r e5a68648db9c src/DLD-FUNCTIONS/matrix_type.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/matrix_type.cc Fri Apr 29 13:10:36 2005 +0000 @@ -0,0 +1,288 @@ +/* + +Copyright (C) 2005 David Bateman + +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 2, 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 this program; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include "ov.h" +#include "defun-dld.h" +#include "error.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" +#include "SparseType.h" + +DEFUN_DLD (matrix_type, args, , + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{type} =} matrix_type (@var{a})\n\ +@deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, @var{type})\n\ +@deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, 'upper', @var{perm})\n\ +@deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, 'lower', @var{perm})\n\ +@deftypefnx {Loadable Function} {@var{a} =} matrix_type (@var{a}, 'banded', @var{nl}, @var{nu})\n\ +Identify the matrix type or mark a matrix as a particular type. This allows rapid\n\ +for solutions of linear equations involving @var{a} to be performed. Called with a\n\ +single argument, @code{matrix_type} returns the type of the matrix and caches it for\n\ +future use. Called with more than one argument, @code{matrix_type} allows the type\n\ +of the matrix to be defined.\n\ +\n\ +The possible matrix types depend on whether the matrix is full or sparse, and can be\n\ +one of the following\n\ +\n\ +@table @asis\n\ +@item 'unknown'\n\ +Remove any previously cached matrix type, and mark type as unknown\n\ +\n\ +@item 'full'\n\ +Mark the matrix as full.\n\ +\n\ +@item 'positive definite'\n\ +Full positive definite matrix.\n\ +\n\ +@item 'diagonal'\n\ +Diagonal Matrix. (Sparse matrices only)\n\ +\n\ +@item 'permuted diagonal'\n\ +Permuted Diagonal matrix. The permutation does not need to be specifically\n\ +indicated, as the structure of the matrix explicitly gives this. (Sparse matrices\n\ +only)\n\ +\n\ +@item 'upper'\n\ +Upper triangular. If the optional third argument @var{perm} is given, the matrix is\n\ +assumed to be a permuted upper triangular with the permutations defined by the\n\ +vector @var{perm}.\n\ +\n\ +@item 'lower'\n\ +Lower triangular. If the optional third argument @var{perm} is given, the matrix is\n\ +assumed to be a permuted lower triangular with the permutations defined by the\n\ +vector @var{perm}.\n\ +\n\ +@item 'banded'\n\ +@itemx 'banded positive definite'\n\ +Banded matrix with the band size of @var{nl} below the diagonal and @var{nu} above\n\ +it. If @var{nl} and @var{nu} are 1, then the matrix is tridiagonal and treated\n\ +with specialized code. In addition the matrix can be marked as positive definite\n\ +(Sparse matrices only)\n\ +\n\ +@item 'singular'\n\ +The matrix is assumed to be singular and will be treated with a minimum norm solution\n\ +\n\ +@end table\n\ +\n\ +Note that the matrix type will be discovered automatically on the first attempt to\n\ +solve a linear equation involving @var{a}. Therefore @code{matrix_type} is only\n\ +useful to give @sc{Octave} hints of the matrix type. Incorrectly defining the\n\ +matrix type will result in incorrect results from solutions of linear equations,\n\ +and so it is entirely the responsibility of the user to correctly indentify the\n\ +matrix type.\n\ +@end deftypefn") +{ + int nargin = args.length (); + octave_value retval; + + if (nargin == 0) + print_usage ("matrix_type"); + else if (nargin > 4) + error ("matrix_type: incorrect number of arguments"); + else + { + if (args(0).class_name () == "sparse") + { + if (nargin == 1) + { + SparseType mattyp; + const octave_value& rep = args(0).get_rep (); + + if (args(0).type_name () == "sparse complex matrix" ) + { + mattyp = + ((const octave_sparse_complex_matrix &)rep).sparse_type (); + + if (mattyp.is_unknown ()) + { + mattyp = SparseType (args(0).sparse_complex_matrix_value ()); + ((octave_sparse_complex_matrix &)rep).sparse_type (mattyp); + } + } + else + { + mattyp = ((const octave_sparse_matrix &)rep).sparse_type (); + + if (mattyp.is_unknown ()) + { + mattyp = SparseType (args(0).sparse_matrix_value ()); + ((octave_sparse_matrix &)rep).sparse_type (mattyp); + } + } + + int typ = mattyp.type (); + + if (typ == SparseType::Diagonal) + retval = octave_value ("Diagonal"); + else if (typ == SparseType::Permuted_Diagonal) + retval = octave_value ("Permuted Diagonal"); + else if (typ == SparseType::Upper) + retval = octave_value ("Upper"); + else if (typ == SparseType::Permuted_Upper) + retval = octave_value ("Permuted Upper"); + else if (typ == SparseType::Lower) + retval = octave_value ("Lower"); + else if (typ == SparseType::Permuted_Lower) + retval = octave_value ("Permuted Lower"); + else if (typ == SparseType::Banded) + retval = octave_value ("Banded"); + else if (typ == SparseType::Banded_Hermitian) + retval = octave_value ("Banded Positive Definite"); + else if (typ == SparseType::Tridiagonal) + retval = octave_value ("Tridiagonal"); + else if (typ == SparseType::Tridiagonal_Hermitian) + retval = octave_value ("Tridiagonal Positive Definite"); + else if (typ == SparseType::Hermitian) + retval = octave_value ("Positive Definite"); + else if (typ == SparseType::Full) + retval = octave_value ("Full"); + else + // This should never happen!!! + retval = octave_value ("Unknown"); + } + else + { + // Ok, we're changing the matrix type + std::string str_typ = args(1).string_value (); + + // XXX FIXME, why do I have to explicitly call the constructor? + SparseType mattyp = SparseType (); + + int nl = 0; + int nu = 0; + + if (error_state) + error ("Matrix type must be a string"); + else + { + // Use STL function to convert to lower case + transform (str_typ.begin (), str_typ.end (), str_typ.begin (), + tolower); + + if (str_typ == "diagonal") + mattyp.mark_as_diagonal (); + if (str_typ == "permuted diagonal") + mattyp.mark_as_permuted_diagonal (); + else if (str_typ == "upper") + mattyp.mark_as_upper_triangular (); + else if (str_typ == "lower") + mattyp.mark_as_lower_triangular (); + else if (str_typ == "banded" || str_typ == "banded positive definite") + { + if (nargin != 4) + error ("matrix_type: banded matrix type requires 4 arguments"); + else + { + nl = args(2).nint_value (); + nu = args(3).nint_value (); + + if (error_state) + error ("matrix_type: band size must be integer"); + else + { + if (nl == 1 && nu == 1) + mattyp.mark_as_tridiagonal (); + else + mattyp.mark_as_banded (nu, nl); + + if (str_typ == "banded positive definite") + mattyp.mark_as_symmetric (); + } + } + } + else if (str_typ == "positive definite") + { + mattyp.mark_as_full (); + mattyp.mark_as_symmetric (); + } + else if (str_typ == "singular") + mattyp.mark_as_rectangular (); + else if (str_typ == "full") + mattyp.mark_as_full (); + else if (str_typ == "unknown") + mattyp.invalidate_type (); + else + error ("matrix_type: Unknown matrix type %s", str_typ.c_str()); + + if (! error_state) + { + if (nargin == 3 && (str_typ == "upper" || str_typ == "lower")) + { + const ColumnVector perm = + ColumnVector (args (2).vector_value ()); + + if (error_state) + error ("matrix_type: Invalid permutation vector"); + else + { + int len = perm.length (); + dim_vector dv = args(0).dims (); + + if (len != dv(0)) + error ("matrix_type: Invalid permutation vector"); + else + { + OCTAVE_LOCAL_BUFFER (octave_idx_type, p, len); + + for (int i = 0; i < len; i++) + p[i] = (int) (perm (i)); + + if (str_typ == "upper") + mattyp.mark_as_permuted (len, p); + else + mattyp.mark_as_permuted (len, p); + } + } + } + else if (nargin != 2 && str_typ != "banded positive definite" && + str_typ != "banded") + error ("matrix_type: Invalid number of arguments"); + + if (! error_state) + { + // Set the matrix type + if (args(0).type_name () == "sparse complex matrix" ) + retval = + octave_value (args(0).sparse_complex_matrix_value (), + mattyp); + else + retval = octave_value (args(0).sparse_matrix_value (), + mattyp); + } + } + } + } + } + else + error ("matrix_type: Only sparse matrices treated at the moment"); + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ diff -r 22994a5730f9 -r e5a68648db9c src/DLD-FUNCTIONS/spkron.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/spkron.cc Fri Apr 29 13:10:36 2005 +0000 @@ -0,0 +1,140 @@ +/* + +Copyright (C) 2002 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 2, 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, write to the Free +Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +02110-1301, USA. + +*/ + +// Author: Paul Kienzle + +#ifdef HAVE_CONFIG_H +#include +#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&, const Sparse&, Sparse&); + +extern void +kron (const Sparse&, const Sparse&, Sparse&); +#endif + +template +void +kron (const Sparse& A, const Sparse& B, Sparse& C) +{ + octave_idx_type idx = 0; + C = Sparse (A.rows () * B.rows (), A.columns () * B.columns (), + A.nnz () * B.nnz ()); + + 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&, const Sparse&, Sparse&); + +template void +kron (const Sparse&, const Sparse&, Sparse&); + +// 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 {Function File} {} 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 (1:4, ones (3, 1))\n\ + @result{} 1 2 3 4\n\ + 1 2 3 4\n\ + 1 2 3 4\n\ +@end group\n\ +@end example\n\ +@end deftypefn") +{ + octave_value_list retval; + + int nargin = args.length (); + + if (nargin != 2 || nargout > 1) + { + print_usage ("kron"); + } + 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; +}