[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: dispatching for sparse
|
From: |
David Bateman |
|
Subject: |
Re: dispatching for sparse |
|
Date: |
Tue, 05 Feb 2008 11:33:08 +0100 |
|
User-agent: |
Thunderbird 2.0.0.6 (X11/20070914) |
John W. Eaton wrote:
> [I'm moving this discussion to the maintainers list. --jwe]
>
> On 4-Jan-2008, David Bateman wrote:
>
> | What we need to do to fix this is to use the newly introduced object
> | classes (ie the "@" directories in this case) and the superiorto
> | function to force mixed sparse/full function to be treated by the sparse
> | versions of the functions. As we are still missing the superiorto
> | function this isn't yet possible and so even converting to object
> | classes won't help here yet... I think this is one that will have to
> | wait at the moment.
>
> As it is currently implemented, the new class-based dispatching won't
> really help here since the class of a sparse object is always double.
>
> One possibility is to extend the class-based dispatch to be finer
> grained so that we can dispatch on "sparse", "sparse complex matrix",
> "sparse bool matrix", etc. (the string returned by typeinfo). But
> this adds even more complexity to the type dispatch code, and I think
> I'd prefer to avoid it.
>
> Another possibility is to merge the sparse code with the existing
> functions and dispatch internally. For example, we would eliminate
> spkron and just have kron, which then needs to do the dispatching
> itself for the various "internal" types.
>
>
Well, I'd thought to move kron the the Array<T> and Sparse<T> classes
like I did with sort. However, I'm not sure that is the right thing to
do.. What about the attached instead.
D.
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
*** ./src/DLD-FUNCTIONS/kron.cc.orig15 2008-02-05 11:21:26.165515769 +0100
--- ./src/DLD-FUNCTIONS/kron.cc 2008-02-05 11:21:04.114646257 +0100
***************
*** 69,74 ****
--- 69,119 ----
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,124 ****
{
print_usage ();
}
! else 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)
{
! ComplexMatrix c;
! kron (a, b, c);
! retval(0) = c;
}
}
! else
{
! Matrix a (args(0).matrix_value ());
! Matrix b (args(1).matrix_value ());
! if (! error_state)
{
! Matrix c;
! kron (a, b, c);
! retval (0) = c;
}
}
--- 142,199 ----
{
print_usage ();
}
! else if (args(0).is_sparse_type () || args(1).is_sparse_type ())
{
! 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;
! }
}
}
! else
{
! 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)
! {
! ComplexMatrix c;
! kron (a, b, c);
! retval(0) = c;
! }
! }
! else
{
! Matrix a (args(0).matrix_value ());
! Matrix b (args(1).matrix_value ());
!
! if (! error_state)
! {
! Matrix c;
! kron (a, b, c);
! retval (0) = c;
! }
}
}
*** ./src/DLD-FUNCTIONS/spkron.cc.orig15 2007-10-13 06:52:19.000000000
+0200
--- ./src/DLD-FUNCTIONS/spkron.cc 2008-02-05 11:21:38.703300097 +0100
***************
*** 1,161 ****
- /*
-
- 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 <address@hidden>
-
- #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: ***
- */
--- 0 ----
*** ./src/Makefile.in.orig15 2008-02-05 11:27:35.488486993 +0100
--- ./src/Makefile.in 2008-02-05 02:01:20.746090412 +0100
***************
*** 67,73 ****
gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
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 sort.cc sparse.cc \
spchol.cc spdet.cc spfind.cc spkron.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 \
--- 67,73 ----
gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \
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 \
sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \
urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \
*** doc/interpreter/sparse.txi.orig15 2008-02-05 11:30:21.978989515 +0100
--- doc/interpreter/sparse.txi 2008-02-05 11:30:55.131383617 +0100
***************
*** 489,495 ****
@item Linear algebra:
@dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv},
! @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, @dfn{spkron},
@dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, @dfn{condest},
@dfn{sprank}
@c @dfn{spaugment}
--- 489,495 ----
@item Linear algebra:
@dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv},
! @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv},
@dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, @dfn{condest},
@dfn{sprank}
@c @dfn{spaugment}
***************
*** 848,855 ****
@DOCSTRING(spinv)
- @DOCSTRING(spkron)
-
@DOCSTRING(splchol)
@DOCSTRING(splu)
--- 848,853 ----
2008-02-05 David Bateman <address@hidden>
* 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-05 David Bateman <address@hidden>
* interpreter/sparse.txi: Remove references to spkron.
- Re: dispatching for sparse,
David Bateman <=