dispatching for sparse (was: Re: kron with a single ...)

[John W. Eaton]

[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.

Comments?

jwe

[dbateman]

John W. Eaton wrote:

[I'm moving this discussion to the maintainers list.  --jwe]

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.

Comments?

jwe

Sure, then what I would suggest doing is to create a kron method of the Array and Sparse classes and the octave_value classes and use them rather than testing for the type. Exactly like what you say you want for sort, min, max etc for 3.1. This is on my todo list so I'll add kron as well

D.

[David Bateman]

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                                David.Bateman@...
Motorola Labs - Paris                        [phone] (Ph)
Parc Les Algorithmes, Commune de St Aubin    [phone] (Mob)
91193 Gif-Sur-Yvette FRANCE                  [phone] (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 <pkienzle@...>
-
- #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  <dbateman@...>

        * 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  <dbateman@...>

        * interpreter/sparse.txi: Remove references to spkron.

[John W. Eaton]

On  5-Feb-2008, David Bateman wrote:

| 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.

I applied this patch and committed it.

I also added a scripts/deprecated/spkron.m to help ease the transition
for anyone who has been using spkron directly.

Thanks,

jwe