[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Operations on sparse matricies (in 2.9.2)

From: David Bateman
Subject: Re: Operations on sparse matricies (in 2.9.2)
Date: Thu, 28 Apr 2005 22:03:10 +0200
User-agent: Mozilla Thunderbird 0.8 (X11/20040923)

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



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


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

#include "dSparse.h"
#include "CSparse.h"
#include "quit.h"

#include "defun-dld.h"
#include "error.h"
#include "oct-obj.h"

extern void
kron (const Sparse<double>&, const Sparse<double>&, Sparse<double>&);

extern void
kron (const Sparse<Complex>&, const Sparse<Complex>&, Sparse<Complex>&);

template <class T>
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 = (Ai);

            for (octave_idx_type Bi = B.cidx (Bj); Bi < B.cidx (Bj+1); Bi++)
       (idx) = v * (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\
x = [a(i, j) b]\n\
@end example\n\
For example,\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;
      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;

reply via email to

[Prev in Thread] Current Thread [Next in Thread]