changeset 6927:50c0aa589d29 octave-forge

Add preliminary 'kronprod' class (still not quite functional)
author hauberg
date Thu, 25 Mar 2010 03:25:13 +0000
parents 56984eafa09d
children 14ac01869011
files main/linear-algebra/INDEX main/linear-algebra/inst/@kronprod/columns.m main/linear-algebra/inst/@kronprod/ctranspose.m main/linear-algebra/inst/@kronprod/det.m main/linear-algebra/inst/@kronprod/disp.m main/linear-algebra/inst/@kronprod/display.m main/linear-algebra/inst/@kronprod/eig.m main/linear-algebra/inst/@kronprod/full.m main/linear-algebra/inst/@kronprod/inv.m main/linear-algebra/inst/@kronprod/iscomplex.m main/linear-algebra/inst/@kronprod/ismatrix.m main/linear-algebra/inst/@kronprod/isreal.m main/linear-algebra/inst/@kronprod/issparse.m main/linear-algebra/inst/@kronprod/issquare.m main/linear-algebra/inst/@kronprod/kronprod.m main/linear-algebra/inst/@kronprod/minus.m main/linear-algebra/inst/@kronprod/mldivide.m main/linear-algebra/inst/@kronprod/mpower.m main/linear-algebra/inst/@kronprod/mtimes.m main/linear-algebra/inst/@kronprod/numel.m main/linear-algebra/inst/@kronprod/plus.m main/linear-algebra/inst/@kronprod/rank.m main/linear-algebra/inst/@kronprod/rdivide.m main/linear-algebra/inst/@kronprod/rows.m main/linear-algebra/inst/@kronprod/size.m main/linear-algebra/inst/@kronprod/sparse.m main/linear-algebra/inst/@kronprod/svd.m main/linear-algebra/inst/@kronprod/times.m main/linear-algebra/inst/@kronprod/trace.m main/linear-algebra/inst/@kronprod/transpose.m main/linear-algebra/inst/@kronprod/uminus.m main/linear-algebra/inst/@kronprod/uplus.m
diffstat 32 files changed, 1383 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- a/main/linear-algebra/INDEX	Wed Mar 24 20:46:28 2010 +0000
+++ b/main/linear-algebra/INDEX	Thu Mar 25 03:25:13 2010 +0000
@@ -17,4 +17,36 @@
 Iterative techniques
  bicg pgmres
 
+Kronecker Products
+ @kronprod/kronprod
+ @kronprod/rows
+ @kronprod/columns
+ @kronprod/numel
+ @kronprod/size
+ @kronprod/display
+ @kronprod/disp
+ @kronprod/issquare
+ @kronprod/issparse
+ @kronprod/isreal
+ @kronprod/iscomplex
+ @kronprod/ismatrix
+ @kronprod/sparse
+ @kronprod/full
+ @kronprod/transpose
+ @kronprod/ctranspose
+ @kronprod/uplus
+ @kronprod/uminus
+ @kronprod/plus
+ @kronprod/minus
+ @kronprod/times
+ @kronprod/rdivide
+ @kronprod/mldivide
+ @kronprod/mtimes
+ @kronprod/mpower
+ @kronprod/inv
+ @kronprod/det
+ @kronprod/rank
+ @kronprod/trace
+ @kronprod/eig
+ @kronprod/svd
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/columns.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,32 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} columns (@var{KP})
+## Return the number of columns in the Kronecker product @var{KP}.
+## @seealso{@@kronprod/rows, @@kronprod/size, @@kronprod/numel}
+## @end deftypefn
+
+function retval = columns (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("columns: input argument must be of class 'kronprod'");
+  endif
+
+  retval = columns (KP.A) * columns (KP.B);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/ctranspose.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,37 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} ctranspose (@var{KP})
+## The complex conjugate transpose of a Kronecker product. This is equivalent
+## to
+## 
+## @example
+## @var{KP}'
+## @end example
+## @seealso{ctranspose, @@kronprod/transpose}
+## @end deftypefn
+
+function retval = ctranspose (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("ctranspose: input argument must be of class 'kronprod'");
+  endif
+
+  retval = kronprod (ctranspose (KP.A), ctranspose (KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/det.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,56 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} det (@var{KP})
+## Compute the determinant of a Kronecker product.
+##
+## If @var{KP} is the Kronecker product of the @var{n}-by-@var{n} matrix @var{A}
+## and the @var{q}-by-@var{q} matrix @var{B}, then the determinant is computed
+## as
+##
+## @example
+## det (@var{A})^q * det (@var{B})^n
+## @end example
+##
+## If @var{KP} is not a Kronecker product of square matrices the determinant is
+## computed by forming the full matrix and then computing the determinant.
+## @seealso{det, @@kronprod/trace, @@kronprod/rank, @@kronprod/full}
+## @end deftypefn
+
+function retval = det (KP)
+  ## Check input
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("det: input argument must be of class 'kronprod'");
+  endif
+
+  if (!issquare (KP))
+    error ("det: argument must be a square matrix");
+  endif
+
+  ## Take action
+  [n, m] = size (KP.A);
+  [q, r] = size (KP.B);
+  if (n == m && q == r) # A and B are both square
+    retval = (det (KP.A)^q) * (det (KP.B)^n);
+  elseif (n*q == m*r) # kron (A, B) is square
+    ## XXX: Can we do something smarter here? We should be able to use the SVD...
+    retval = det (full (KP));
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/disp.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,28 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} disp (@var{KP})
+## Show the content of the Kronecker product @var{KP}. To avoid evaluating the
+## Kronecker product, this function displays the two matrices defining the product.
+## To display the actual values of @var{KP}, use @code{disp (full (@var{KP}))}.
+##
+## This function is equivalent to @code{@@kronprod/display}.
+## @seealso{@@kronprod/display, @@kronprod/full}
+## @end deftypefn
+
+function disp (KP)
+  display (KP);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/display.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,38 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} display (@var{KP})
+## Show the content of the Kronecker product @var{KP}. To avoid evaluating the
+## Kronecker product, this function displays the two matrices defining the product.
+## To display the actual values of @var{KP}, use @code{display (full (@var{KP}))}.
+## @seealso{@@kronprod/displ, @@kronprod/full}
+## @end deftypefn
+
+function display (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("display: input argument must be of class 'kronprod'");
+  endif
+  
+  disp ("Kronecker Product of A and B with");
+  disp ("A = ");
+  disp (KP.A);
+  disp ("B = ");
+  disp (KP.B);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/eig.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,71 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{lambda} =} eig (@var{KP})
+## @deftypefnx{Function File} {[var{V}, @var{lambda}] =} eig (@var{KP})
+## XXX: Write help text
+## @seealso{eig, @kronprod/svd}
+## @end deftypefn
+
+function [V, lambda] = eig (KP, A)
+  ## XXX: This implementation provides a different permutation of eigenvalues and
+  ## eigenvectors compared to 'eig (full (KP))'
+
+  ## Check input
+  if (nargin == 0 || nargin > 2)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("eig: first input argument must be of class 'kronprod'");
+  endif
+  
+  if (!issquare (KP))
+    error ("eig: first input must be a square matrix");
+  endif
+  
+  ## Take action
+  if (nargin == 1)
+    if (nargout <= 1)
+      ## Only eigenvalues were requested
+      if (issquare (KP.A) && issquare (KP.B))
+        lambda_A = eig (KP.A);
+        lambda_B = eig (KP.B);
+        V = kronprod (lambda_A, lambda_B);
+      else
+        ## We should be able to do this using SVD
+        error ("eig not implemented (yet) for Kronecker products of non-square matrices");
+      endif
+
+    elseif (nargout == 2)
+      ## Both eigenvectors and eigenvalues were requested
+      if (issquare (KP.A) && issquare (KP.B))
+        [V_A, lambda_A] = eig (KP.A);
+        [V_B, lambda_B] = eig (KP.B);
+        V = kronprod (V_A, V_B);
+        lambda = kronprod (lambda_A, lambda_B);
+      else
+        ## We should be able to do this using SVD
+        error ("eig not implemented (yet) for Kronecker products of non-square matrices");
+      endif
+    endif
+    
+  elseif (nargin == 2)
+    ## Solve generalised eigenvalue problem
+    ## XXX: Is there a fancy way of doing this?
+    [V, lambda] = eig (full (KP), full (A));
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/full.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,37 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} full (@var{KP})
+## Return the full matrix representation of the Kronecker product @var{KP}.
+##
+## If @var{KP} is the Kronecker product of an @var{n}-by-@var{m} matrix and a
+## @var{q}-by-@var{r} matrix, then the result is a @var{n}@var{q}-by-@var{m}@var{r}
+## matrix. Thus, the result can require vast amount of memory, so this function
+## should be avoided whenever possible.
+## @seealso{full, @kronprod/sparse}
+## @end deftypefn
+
+function retval = full (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("full: input argument must be of class 'kronprod'");
+  endif
+  
+  retval = full (kron (KP.A, KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/inv.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,52 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} inv (@var{KP})
+## Return the inverse of the Kronecker product @var{KP}.
+##
+## If @var{KP} is the Kronecker product of two square matrices @var{A} and @var{B},
+## the inverse will be computed as the Kronecker product of the inverse of
+## @var{A} and @var{B}.
+##
+## If @var{KP} is square but not a Kronecker product of square matrices, the
+## inverse will be computed using the SVD
+## @seealso{@kronprod/sparse}
+## @end deftypefn
+
+function retval = inv (KP)
+  ## Check input
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("inv: input argument must be of class 'kronprod'");
+  endif
+  
+  ## Do the computations
+  [n, m] = size (KP.A);
+  [q, r] = size (KP.B);
+  if (n == m && q == r) # A and B are both square
+    retval = kronprod (inv (KP.A), inv (KP.B));
+  elseif (n*q == m*r) # kron (A, B) is square
+    ## We use the SVD to compute the inverse.
+    ## XXX: Should we use 'eig' instead?
+    [U, S, V] = svd (KP);
+    retval = U * (1./S) * V';
+  else
+    error ("inv: argument must be a square matrix");
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/iscomplex.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,32 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} iscomplex (@var{KP})
+## Return @t{true} if the Kronecker product @var{KP} contains any complex values.
+## @seealso{iscomplex, @kronprod/isreal}
+## @end deftypefn
+
+function retval = iscomplex (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("iscomplex: input argument must be of class 'kronprod'");
+  endif
+  
+  retval = (iscomplex (KP.A) || iscomplex (KP.B);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/ismatrix.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,24 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} ismatrix (@var{KP})
+## Return @t{true} to indicate that the Kronecker product @var{KP} always is a
+## matrix.
+## @end deftypefn
+
+function retval = ismatrix (KP)
+  retval = true;
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/isreal.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,33 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} isreal (@var{KP})
+## Return @t{true} if the Kronecker product @var{KP} is real, i.e. has no
+## imaginary components.
+## @seealso{isreal, @@kronprod/iscomplex}
+## @end deftypefn
+
+function retval = isreal (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("isreal: input argument must be of class 'kronprod'");
+  endif
+  
+  retval = (isreal (KP.A) & isreal (KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/issparse.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,34 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} issparse (@var{KP})
+## Return @t{true} if one of the matrices in the Kronecker product @var{KP}
+## is sparse.
+## @seealso{@@kronprod/sparse}
+## @end deftypefn
+
+function  retval = issparse (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("issparse: input argument must be of class 'kronprod'");
+  endif
+
+  retval = (issparse(KP.A) || issparse(KP.B));
+endfunction
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/issquare.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,33 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} issquare (@var{KP})
+## Return @t{true} if the Kronecker product @var{KP} is a square matrix.
+## @seealso{@@kronprod/size}
+## @end deftypefn
+
+function retval = issquare (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("issquare: input argument must be of class 'kronprod'");
+  endif
+  
+  s = size (KP);
+  retval = (s (1) == s (2));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/kronprod.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,34 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} kronprod (@var{KP})
+## Construct a Kronecker product object.
+##
+## XXX: Write documentation
+## @end deftypefn
+
+function retval = kronprod (A, B)
+  if (nargin == 0)
+    KP.A = KP.B = KP.full = [];
+  elseif (nargin == 2 && ismatrix (A) && ismatrix (B))
+    KP.A = A;
+    KP.B = B;
+  else
+    print_usage ();
+  endif
+  
+  retval = class (KP, "kronprod");
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/minus.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,52 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} minus (@var{a}, @var{a})
+## Return the difference between a Kronecker product and another matrix. This is performed
+## by forming the full matrix of both inputs and is as such a potential expensive
+## operation.
+## @seealso{minus, @@kronprod/plus}
+## @end deftypefn
+
+function retval = minus (M1, M2)
+  if (nargin != 2)
+    print_usage ();
+  endif
+  
+  if (!ismatrix (M1) || !ismatrix (M2))
+    error ("minus: both input arguments must be matrices");
+  endif
+  
+  if (!size_equal (M1, M2))
+    error ("minus: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (M1), columns (M1), rows (M2), columns (M2));
+  endif
+
+  ## XXX: Can we do something smarter here?
+  if (issparse (M1))
+    M1 = sparse (M1);
+  else
+    M1 = full (M1);
+  endif
+  
+  if (issparse (M2))
+    M2 = sparse (M2);
+  else
+    M2 = full (M2);
+  endif
+  
+  retval = M1 - M2;
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/mldivide.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,59 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} mldivide (@var{M1}, @var{M2})
+## XXX: Write documentation
+## @end deftypefn
+
+function retval = mldivide (M1, M2)
+  ## Check input
+  if (nargin != 2)
+    print_usage ();
+  endif
+  
+  if (!ismatrix (M1) || !ismatrix (M2))
+    error ("mldivide: both input arguments must be matrices");
+  endif
+  
+  if (rows (M1) != rows (M2)
+    error ("mldivide: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (M1), columns (M1), rows (M2), columns (M2));
+  endif
+
+  ## Take action depending on types
+  M1_is_KP = isa (M1, "kronprod");
+  M2_is_KP = isa (M2, "kronprod");
+  
+  if (M1_is_KP && M2_is_KP) # Left division of Kronecker Products
+    error ("mldividide: this part not yet implemented as I'm lazy...");
+
+  elseif (M1_is_KP) # Left division of Kronecker Product and Matrix
+    ## XXX: Does this give the same minimum-norm solution as when using
+    ## XXX:   full (M1) \ M2
+    ## XXX: ? It is the same when M1 is invertible.
+    retval = zeros (columns (M1), columns (M2));
+    for n = 1:columns (M2)
+      M = reshape (M2 (:, n), [rows(M1.B), rows(M1.A)]);
+      retval (:, n) = vec ((M1.A \ (M1.B \ M)')');
+    endfor
+
+  elseif (M2_is_KP) # Left division of Matrix and Kronecker Product
+    error ("mldividide: this part not yet implemented as I'm lazy...");
+      
+  else
+    error ("mldivide: internal error for 'kronprod'");
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/mpower.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,44 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} mpower (@var{KP}, @var{k})
+## XXX: Write documentation
+## @end deftypefn
+
+function retval = mpower (KP, k)
+  ## Check input
+  if (nargin != 2)
+    print_usage ();
+  endif
+  
+  if (!ismatrix (KP))
+    error ("mpower: first input argument must be a matrix");
+  endif
+  
+  if (!isscalar (k))
+    error ("mpower: second input argument must be a scalar");
+  endif
+
+  ## Do the actual computation
+  if (issquare (KP.A) && issquare (KP.B) && k == round (k))
+    retval = kronprod (KP.A^k, KP.B^k);
+  elseif (issquare (KP))
+    ## XXX: Can we do something smarter here?
+    retval = full (KP)^k;
+  else
+    error ("for A^b, A must be square");
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/mtimes.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,92 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} mtimes (@var{KP})
+## XXX: Write documentation
+## @end deftypefn
+
+function retval = mtimes (M1, M2)
+  ## Check input
+  if (nargin == 0)
+    print_usage ();
+  elseif (nargin == 1)
+    ## This seems to be what happens for full and sparse matrices, so we copy this behaviour
+    retval = M1;
+    return;
+  endif
+  
+  if (!ismatrix (M1) || !ismatrix (M2))
+    error ("mtimes: input arguments must be matrices");
+  endif
+  
+  if (columns (M1) != rows (M2))
+    error ("mtimes: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (M1), columns (M1), rows (M2), columns (M2));
+  endif
+
+  ## Take action depending on input types
+  M1_is_KP = isa (M1, "kronprod");
+  M2_is_KP = isa (M2, "kronprod");
+  
+  if (M1_is_KP && M2_is_KP) # Product of Kronecker Products
+    ## Check if the size match such that the result is a Kronecker Product
+    if (columns (M1.A) == rows (M2.A) && columns (M1.B) == rows (M2.B))
+      retval = kronprod (M1.A * M2.A, M1.B * M2.B);
+    else
+      ## Form the full matrix of the smallest matrix and use that to compute the
+      ## final product
+      ## XXX: Can we do something smarter here?
+      numel1 = numel (M1);
+      numel2 = numel (M2);
+      if (numel1 < numel2)
+        retval = full (M1) * M2;
+      else
+        retval = M1 * full (M2);
+      endif
+    endif
+    
+  elseif (M1_is_KP && isscalar (M2)) # Product of Kronecker Product and scalar
+    if (numel (M1.A) < numel (M1.B))
+      retval = kronprod (M2 * M1.A, M1.B);
+    else
+      retval = kronprod (M1.A, M2 * M1.B);
+    endif
+    
+  elseif (M1_is_KP && ismatrix (M2)) # Product of Kronecker Product and Matrix
+    retval = zeros (rows (M1), columns (M2));
+    for n = 1:columns (M2)
+      M = reshape (M2 (:, n), [columns(M1.B), columns(M1.A)]);
+      retval (:, n) = vec (M1.B * M * M1.A');
+    endfor
+  
+  elseif (isscalar (M1) && M2_is_KP) # Product of scalar and Kronecker Product
+    if (numel (M2.A) < numel (M2.B))
+      retval = kronprod (M1 * M2.A, M2.B);
+    else
+      retval = kronprod (M2.A, M1 * M2.B);
+    endif
+    
+  elseif (ismatrix (M1) && M2_is_KP) # Product of Matrix and Kronecker Product
+    retval = zeros (rows (M1), columns (M2));
+    for n = 1:rows (M1)
+      M = reshape (M1 (n, :), [rows(M2.B), rows(M2.A)]);
+      retval (n, :) = vec (M2.B' * M * M2.A);
+    endfor
+      
+  else
+    error ("mtimes: internal error for 'kronprod'");
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/numel.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,32 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} numel (@var{KP})
+## Return the number of elements in the Kronecker product @var{KP}.
+## @seealso{numel, @@kronprod/rows, @@kronprod/columns, @@kronprod/size}
+## @end deftypefn
+
+function retval = numel (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("numel: input must be of class 'kronprod'");
+  endif
+  
+  retval = prod (size (KP.A) .* size (KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/plus.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,56 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} plus (@var{a}, @var{a})
+## Return the sum of a Kronecker product and another matrix. This is performed
+## by forming the full matrix of both inputs and is as such a potential expensive
+## operation.
+## @seealso{plus, @@kronprod/minus}
+## @end deftypefn
+
+function retval = plus (M1, M2)
+  if (nargin == 0 || nargin > 2)
+    print_usage ();
+  elseif (nargin == 1)
+    ## This seems to be the behaviour for the built-in types so we copy that
+    retval = M1;
+    return;
+  endif
+  
+  if (!ismatrix (M1) || !ismatrix (M2))
+    error ("plus: input arguments must be matrics");
+  endif
+  
+  if (!size_equal (M1, M2))
+    error ("plus: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (M1), columns (M1), rows (M2), columns (M2));
+  endif
+  
+  ## XXX: Can we do something smarter here?
+  if (issparse (M1))
+    M1 = sparse (M1);
+  else
+    M1 = full (M1);
+  endif
+  
+  if (issparse (M2))
+    M2 = sparse (M2);
+  else
+    M2 = full (M2);
+  endif
+  
+  retval = M1 + M2;
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/rank.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,33 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} rank (@var{KP})
+## Return the rank of the Kronecker product @var{KP}. This is computed as the
+## product of the ranks of the matrices forming the product.
+## @seealso{rank, @@kronprod/det, @@kronprod/trace}
+## @end deftypefn
+
+function retval = rank (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("rank: input must be of class 'kronprod'");
+  endif
+  
+  retval = rank (KP.A) * rank (KP.B);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/rdivide.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,53 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} rdivide (@var{a}, @var{b})
+## XXX: Write help text.
+## @end deftypefn
+
+function retval = rdivide (a, b)
+  ## Check input
+  if (nargin < 2)
+    print_usage ();
+  endif
+  
+  if (!ismatrix (a) || !ismatrix (a))
+    error ("rdivide: input arguments must be scalars or matrices");
+  endif
+  
+  if (!size_equal (a, b) || !isscalar (b))
+    error ("times: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (a), columns (a), rows (b), columns (b));
+  endif
+
+  ## Take action depending on input
+  if (isscalar (a) && isa (b, "kronprod"))
+    retval = kronprod (a ./ b.A, 1 ./ b.B);
+  elseif (isa (a, "kronprod") && isscalar (b))
+    if (numel (a.A) < numel (a.B))
+      retval = kronprod (a.A ./ b, a.B);
+    else
+      retval = kronprod (a.A, a.B ./ b);
+    endif
+  elseif (isa (a, "kronprod") && isa (b, "kronprod"))
+    ## XXX: Can we do something smarter here?
+    retval = full (a) ./ full (b);
+  else
+    ## XXX: We should probably handle sparse cases and all sorts of other
+    ## situations better here
+    retval = full (a) ./ full (b);
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/rows.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,32 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} rows (@var{KP})
+## Return the number of rows in the Kronecker product @var{KP}.
+## @seealso{rows, @@kronprod/size, @@kronprod/columns, @@kronprod/numel}
+## @end deftypefn
+
+function retval = rows (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("rows: input must be of class 'kronprod'");
+  endif
+  
+  retval = rows (KP.A) * rows (KP.B);
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/size.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,39 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} size (@var{KP})
+## @deftypefnx{Function File} size (@var{KP}, @var{dim})
+## Return the size of the Kronecker product @var{KP} as a vector.
+## @seealso{size, @@kronprod/rows, @@kronprod/columns, @@kronprod/numel}
+## @end deftypefn
+
+function retval = size (KP, dim)
+  if (nargin < 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("size: input must be of class 'kronprod'");
+  endif
+  if (nargin > 1 && !(isscalar (dim) && dim == round (dim) && dim > 0))
+    error ("size: optional second input must be a positive integer");
+  endif
+  
+  retval = size (KP.A) .* size (KP.B);
+  if (nargin > 1)
+    retval = retval (dim);
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/sparse.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,33 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} sparse (@var{KP})
+## Return the Kronecker product @var{KP} represented as a sparse matrix.
+## @seealso{sparse, @@kronprod/issparse, @@kronprod/full}
+## @end deftypefn
+
+function retval = sparse (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("sparse: input argument must be of class 'kronprod'");
+  endif
+
+  ## XXX: Would this be better? kron (sparse (KP.A), sparse (KP.B)))
+  retval = sparse (kron (KP.A, KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/svd.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,47 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} svd (@var{KP})
+## XXX: Write documentation
+## @end deftypefn
+
+function [U, S, V] = svd (KP)
+  if (nargin < 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("svd: input must be of class 'kronprod'");
+  endif
+  
+  ## XXX: I don't think this works properly for non-square A and B
+  
+  if (nargout <= 1)
+    ## Only singular values were requested
+    S_A = svd (KP.A);
+    S_B = svd (KP.B);
+    U = kronprod (S_A, S_B);
+  elseif (nargout == 3)
+    ## The full SVD was requested
+    [U_A, S_A, V_A] = svd (KP.A);
+    [U_B, S_B, V_B] = svd (KP.B);
+    U = kronprod (U_A, U_B);
+    S = kronprod (S_A, S_B);
+    V = kronprod (V_A, V_B);
+  else
+    print_usage ();
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/times.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,97 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} times (@var{KP})
+## XXX: Write documentation
+## @end deftypefn
+
+function retval = times (M1, M2)
+  ## Check input
+  if (nargin == 0)
+    print_usage ();
+  elseif (nargin == 1)
+    ## This seems to be what happens for full and sparse matrices, so we copy this behaviour
+    retval = M1;
+    return;
+  endif
+  
+  if (!ismatrix (M1) || !ismatrix (M2))
+    error ("times: input arguments must be matrices");
+  endif
+  
+  if (!size_equal (M1, M2))
+    error ("times: nonconformant arguments (op1 is %dx%d, op2 is %dx%d)",
+           rows (M1), columns (M1), rows (M2), columns (M2));
+  endif
+
+  ## Take action depending on input types
+  M1_is_KP = isa (M1, "kronprod");
+  M2_is_KP = isa (M2, "kronprod");
+  
+  if (M1_is_KP && M2_is_KP) # Product of Kronecker Products
+    ## Check if the size match such that the result is a Kronecker Product
+    if (size_equal (M1.A, M2.A) && size_equal (M1.B, M2.B))
+      retval = kronprod (M1.A .* M2.A, M1.B .* M2.B);
+    else
+      ## Form the full matrix or sparse matrix of both matrices
+      ## XXX: Can we do something smarter here?
+      if (issparse (M1))
+        M1 = sparse (M1);
+      else
+        M1 = full (M1);
+      endif
+      
+      if (issparse (M2))
+        M2 = sparse (M2);
+      else
+        M2 = full (M2);
+      endif
+      
+      retval = M1 .* M2;
+    endif
+    
+  elseif (M1_is_KP && isscalar (M2)) # Product of Kronecker Product and scalar
+    if (numel (M1.A) < numel (M1.B))
+      retval = kronprod (M2 * M1.A, M1.B);
+    else
+      retval = kronprod (M1.A, M2 * M1.B);
+    endif
+    
+  elseif (M1_is_KP && ismatrix (M2)) # Product of Kronecker Product and Matrix
+    retval = zeros (rows (M1), columns (M2));
+    for n = 1:columns (M2)
+      M = reshape (M2 (:, n), [columns(M1.B), columns(M1.A)]);
+      retval (:, n) = vec (M1.B * M * M1.A');
+    endfor
+  
+  elseif (isscalar (M1) && M2_is_KP) # Product of scalar and Kronecker Product
+    if (numel (M2.A) < numel (M2.B))
+      retval = kronprod (M1 * M2.A, M2.B);
+    else
+      retval = kronprod (M2.A, M1 * M2.B);
+    endif
+    
+  elseif (ismatrix (M1) && M2_is_KP) # Product of Matrix and Kronecker Product
+    retval = zeros (rows (M1), columns (M2));
+    for n = 1:rows (M1)
+      M = reshape (M1 (n, :), [rows(M2.B), rows(M2.A)]);
+      retval (n, :) = vec (M2.B' * M * M2.A);
+    endfor
+      
+  else
+    error ("mtimes: internal error for 'kronprod'");
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/trace.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,41 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} trace (@var{KP})
+## Returns the trace of the Kronecker product @var{KP}.
+##
+## If @var{KP} is a Kronecker product of two square matrices, the trace is
+## computed as the product of the trace of these two matrices. Otherwise the
+## trace is computed by forming the full matrix.
+## @seealso{@@kronprod/det, @@kronprod/rank, @@kronprod/full}
+## @end deftypefn
+
+function retval = trace (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("trace: input must be of class 'kronprod'");
+  endif
+  
+  if (issquare (KP.A) && issquare (KP.B))
+    retval = trace (KP.A) * trace (KP.B);
+  else
+    ## XXX: Can we do something smarter here? Using 'eig' or 'svd'?
+    retval = trace (full (KP));
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/transpose.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,37 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} transpose (@var{KP})
+## Returns the transpose of the Kronecker product @var{KP}. This is equivalent
+## to
+## 
+## @example
+## @var{KP}.'
+## @end example
+## @seealso{transpose, @@kronprod/ctranspose}
+## @end deftypefn
+
+function retval = transpose (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("transpose: input must be of class 'kronprod'");
+  endif
+  
+  retval = kronprod (transpose (KP.A), transpose (KP.B));
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/uminus.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,38 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} uminus (@var{KP})
+## Returns the unary minus operator working on the Kronecker product @var{KP}.
+## This corresponds to @code{-@var{KP}} and simply returns the Kronecker
+## product with the sign of the smallest matrix in the product reversed.
+## @seealso{@@kronprod/uminus}
+## @end deftypefn
+
+function KP = uplus (KP)
+  if (nargin != 1)
+    print_usage ();
+  endif
+  
+  if (!isa (KP, "kronprod"))
+    error ("uplus: input must be of class 'kronprod'");
+  endif
+  
+  if (numel (KP.A) < numel (KP.B))
+    KP.A *= -1;
+  else
+    KP.B *= -1;
+  endif
+endfunction
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/linear-algebra/inst/@kronprod/uplus.m	Thu Mar 25 03:25:13 2010 +0000
@@ -0,0 +1,25 @@
+## Copyright (C) 2010  Soren Hauberg
+## 
+## This program 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, or (at your option)
+## any later version.
+## 
+## This program 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 file.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} uplus (@var{KP})
+## Returns the unary plus operator working on the Kronecker product @var{KP}.
+## This corresponds to @code{+@var{KP}} and simply returns the Kronecker
+## product unchanged.
+## @seealso{@@kronprod/uminus}
+## @end deftypefn
+
+function KP = uplus (KP)
+endfunction