Re: Operations on sparse matricies (in 2.9.2)

[David Bateman]

Shaun Cloherty wrote:

> I recently downloaded version 2.9.2 to explore its support for sparse
> matrices. In my script I have a number of operations like this;
>
> AM = kron(speye(N), D2M);
>
> where N = 130 and D2M is a 82 x 82 sparse matrix.
>
> I kind of expected the result (AM) to be a 10660 x 10660 sparse matrix
> but octave returns a 10660 x 10660 full matrix (occupying something like
> 900Mb) which, to my naive understanding, seems to negate to some extent
> the benefit of using sparse matrices. I am a complete novice on the
> subject of sparse matrices so I have little idea if expectation is
> reasonable or not. Am I missing something here or should I file a bug
> report?
>
> Regards,
>
> Shaun
>
>  
Kron is an oct-file, and as such is currently using full matrices only. 
Pass it sparse matrices and it'll give you full matrices back. I attach 
a sparse version, that just took a few minutes to adapt from the full 
version and I will check it in soon, but my CVS tree is a mess due to 
the sparse matrix type caching code I'm currently writing and am not 
ready to check in, so it'll be a couple of days till this turns up on 
the CVS.

The PKG_ADD commands will eventually add the alias to allow the attached 
version to be called as kron, but if you build it with mkoctfile you'll 
need to call it as spkron, or run the commands

dispatch ("kron", "spkron", "sparse matrix")
dispatch ("kron", "spkron", "sparse complex matrix")
dispatch ("kron", "spkron", "sparse bool matrix")

prior to calling your scripts to manually install the aliases.. In any 
case I can't imagine that the first thing you tried is a kron!!
You didn't try somthing like A*B or A\B first ? Oh well, you'll never 
know what people will want till you let them loose on your code...

Regards
David

--
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

/*

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 <address@hidden>


// XXX FIXME Comment out HAVE_CONFIG_H so Shaun Cloherty can temporarily build 
// outside of octave source tree with mkoctfile...
//#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.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<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 {Function File} {} spkron (@var{a}, @var{b})\n\
Form the kronecker product of two sparse matrices, defined 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;
}