# HG changeset patch
# User dbateman
# Date 1124972484 0
# Node ID 2042301733ced43a33c51bac0db96b46ccf8ff52
# Parent  2a16423e4aa042b7feb1f03328b9d3b58019ee28
[project @ 2005-08-25 12:21:24 by dbateman]

diff --git a/liboctave/CSparse.cc b/liboctave/CSparse.cc
--- a/liboctave/CSparse.cc
+++ b/liboctave/CSparse.cc
@@ -6361,6 +6361,54 @@
 #endif
 }
 
+ComplexMatrix
+operator * (const ComplexMatrix& m, const SparseMatrix& a)
+{
+  SparseComplexMatrix tmp (a);
+  return m * tmp;
+}
+
+ComplexMatrix
+operator * (const Matrix& m, const SparseComplexMatrix& a)
+{
+  ComplexMatrix tmp (m);
+  return tmp * a;
+}
+
+ComplexMatrix
+operator * (const ComplexMatrix& m, const SparseComplexMatrix& a)
+{
+#ifdef HAVE_SPARSE_BLAS
+  // XXX FIXME XXX Isn't there a sparse BLAS ??
+#else
+  FULL_SPARSE_MUL (ComplexMatrix, Complex);
+#endif
+}
+
+ComplexMatrix
+operator * (const SparseComplexMatrix& m, const Matrix& a)
+{
+  ComplexMatrix tmp (a);
+  return m * tmp;
+}
+
+ComplexMatrix
+operator * (const SparseMatrix& m, const ComplexMatrix& a)
+{
+  SparseComplexMatrix tmp (m);
+  return tmp * a;
+}
+
+ComplexMatrix
+operator * (const SparseComplexMatrix& m, const ComplexMatrix& a)
+{
+#ifdef HAVE_SPARSE_BLAS
+  // XXX FIXME XXX Isn't there a sparse BLAS ??
+#else
+  SPARSE_FULL_MUL (ComplexMatrix, Complex);
+#endif
+}
+
 // XXX FIXME XXX -- it would be nice to share code among the min/max
 // functions below.
 
diff --git a/liboctave/CSparse.h b/liboctave/CSparse.h
--- a/liboctave/CSparse.h
+++ b/liboctave/CSparse.h
@@ -398,6 +398,20 @@
 extern SparseComplexMatrix operator * (const SparseComplexMatrix&, 
 				       const SparseComplexMatrix&);
 
+extern ComplexMatrix operator * (const Matrix&,        
+				       const SparseComplexMatrix&);
+extern ComplexMatrix operator * (const ComplexMatrix&, 
+				       const SparseMatrix&);
+extern ComplexMatrix operator * (const ComplexMatrix&, 
+				       const SparseComplexMatrix&);
+
+extern ComplexMatrix operator * (const SparseMatrix&,        
+				       const ComplexMatrix&);
+extern ComplexMatrix operator * (const SparseComplexMatrix&, 
+				       const Matrix&);
+extern ComplexMatrix operator * (const SparseComplexMatrix&, 
+				       const ComplexMatrix&);
+
 extern SparseComplexMatrix min (const Complex& c, 
 				const SparseComplexMatrix& m);
 extern SparseComplexMatrix min (const SparseComplexMatrix& m, 
diff --git a/liboctave/ChangeLog b/liboctave/ChangeLog
--- a/liboctave/ChangeLog
+++ b/liboctave/ChangeLog
@@ -1,3 +1,12 @@
+2005-08-25  David Bateman <dbateman@free.fr>
+
+	* Sparse-op-defs.h (FULL_SPARSE_MUL, SPARSE_FULL_MUL): Macro for
+	mixed sparse/full multiply.
+	* dSparse.cc (operator *), CSparse.cc (operator *): New operators for
+	mixed sparse/full multiply.
+	* dSparse.h (operator *), CSparse.h (operator *): Declaration of
+	mixed sparse/full multiply operators.
+
 2005-07-25   Erik de Castro Lopo  <erikd@zip.com.au>
 
 	* oct-inttypes.h (OCTAVE_S_US_FTR): Compare <= 0 instead of < 0 to
diff --git a/liboctave/Sparse-op-defs.h b/liboctave/Sparse-op-defs.h
--- a/liboctave/Sparse-op-defs.h
+++ b/liboctave/Sparse-op-defs.h
@@ -1603,6 +1603,67 @@
 	} \
     }
 
+#define SPARSE_FULL_MUL( RET_TYPE, EL_TYPE ) \
+  octave_idx_type nr = m.rows (); \
+  octave_idx_type nc = m.cols (); \
+  \
+  octave_idx_type a_nr = a.rows (); \
+  octave_idx_type a_nc = a.cols (); \
+  \
+  if (nc != a_nr) \
+    { \
+      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
+      return RET_TYPE (); \
+    } \
+  else \
+    { \
+      RET_TYPE retval (nr, a_nc, EL_TYPE ()); \
+      \
+      for (octave_idx_type i = 0; i < a_nc ; i++) \
+	{ \
+	  for (octave_idx_type j = 0; j < a_nr; j++) \
+	    { \
+              OCTAVE_QUIT; \
+	      \
+              EL_TYPE tmpval = a.elem(j,i); \
+	      for (octave_idx_type k = m.cidx(j) ; k < m.cidx(j+1); k++) \
+	        retval.elem (m.ridx(k),i) += tmpval * m.data(k); \
+	    } \
+        } \
+      return retval; \
+    }
+
+#define FULL_SPARSE_MUL( RET_TYPE, EL_TYPE ) \
+  octave_idx_type nr = m.rows (); \
+  octave_idx_type nc = m.cols (); \
+  \
+  octave_idx_type a_nr = a.rows (); \
+  octave_idx_type a_nc = a.cols (); \
+  \
+  if (nc != a_nr) \
+    { \
+      gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
+      return RET_TYPE (); \
+    } \
+  else \
+    { \
+      RET_TYPE retval (nr, a_nc, EL_TYPE ()); \
+      \
+      for (octave_idx_type i = 0; i < a_nc ; i++) \
+	{ \
+	   for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
+	     { \
+	        octave_idx_type col = a.ridx(j); \
+	        EL_TYPE tmpval = a.data(j); \
+                OCTAVE_QUIT; \
+	        \
+	        for (octave_idx_type k = 0 ; k < nr; k++) \
+	          retval.elem (k,i) += tmpval * m.elem(k,col); \
+	    } \
+        } \
+      return retval; \
+    }
+
 #endif
 
 /*
diff --git a/liboctave/dSparse.cc b/liboctave/dSparse.cc
--- a/liboctave/dSparse.cc
+++ b/liboctave/dSparse.cc
@@ -6478,13 +6478,33 @@
 operator * (const SparseMatrix& m, const SparseMatrix& a)
 {
 #ifdef HAVE_SPARSE_BLAS
-  // XXX FIXME XXX Isn't there a sparse BLAS ??
+  // XXX FIXME XXX Isn't there a sparse BLAS ?? Is it faster??
 #else
   // Use Andy's sparse matrix multiply function
   SPARSE_SPARSE_MUL (SparseMatrix, double);
 #endif
 }
 
+Matrix
+operator * (const Matrix& m, const SparseMatrix& a)
+{
+#ifdef HAVE_SPARSE_BLAS
+  // XXX FIXME XXX Isn't there a sparse BLAS ?? Is it faster??
+#else
+  FULL_SPARSE_MUL (Matrix, double);
+#endif
+}
+
+Matrix
+operator * (const SparseMatrix& m, const Matrix& a)
+{
+#ifdef HAVE_SPARSE_BLAS
+  // XXX FIXME XXX Isn't there a sparse BLAS ?? Is it faster??
+#else
+  SPARSE_FULL_MUL (Matrix, double);
+#endif
+}
+
 // XXX FIXME XXX -- it would be nice to share code among the min/max
 // functions below.
 
diff --git a/liboctave/dSparse.h b/liboctave/dSparse.h
--- a/liboctave/dSparse.h
+++ b/liboctave/dSparse.h
@@ -384,6 +384,10 @@
 
 extern SparseMatrix operator * (const SparseMatrix& a, 
 				const SparseMatrix& b);
+extern Matrix operator * (const Matrix& a, 
+				const SparseMatrix& b);
+extern Matrix operator * (const SparseMatrix& a, 
+				const Matrix& b);
 
 extern SparseMatrix min (double d, const SparseMatrix& m);
 extern SparseMatrix min (const SparseMatrix& m, double d);
diff --git a/src/ChangeLog b/src/ChangeLog
--- a/src/ChangeLog
+++ b/src/ChangeLog
@@ -1,3 +1,10 @@
+2005-08-25  David Bateman <dbateman@free.fr>
+
+        * OPERATORS/op-sm-m.cc, OPERATORS/op-sm-cm.cc, OPERATORS/op-scm-m.cc,
+	OPERATORS/op-scm-cm.cc, OPERATORS/op-m-sm.cc, OPERATORS/op-m-scm.cc,
+	OPERATORS/op-cm-sm.cc, OPERATORS/op-cm-scm.cc: Use mixed matrix/sparse
+	multiply operator rather than casting sparse to matrix.
+
 2005-07-18  John W. Eaton  <jwe@octave.org>
 
 	* strfns.cc (Fstrcmp): New function from Soren Hauberg
diff --git a/src/OPERATORS/op-cm-scm.cc b/src/OPERATORS/op-cm-scm.cc
--- a/src/OPERATORS/op-cm-scm.cc
+++ b/src/OPERATORS/op-cm-scm.cc
@@ -43,15 +43,7 @@
 DEFBINOP_OP (add, complex_matrix, sparse_complex_matrix,+)
 DEFBINOP_OP (sub, complex_matrix, sparse_complex_matrix,-)
 
-DEFBINOP (mul, complex_matrix, sparse_complex_matrix)
-{
-  CAST_BINOP_ARGS (const octave_complex_matrix&, 
-		   const octave_sparse_complex_matrix&);
-  
-  ComplexMatrix tmp (v2.complex_matrix_value ());
-
-  return octave_value ( v1.complex_matrix_value() * tmp);
-}
+DEFBINOP_OP (mul, complex_matrix, sparse_complex_matrix, *)
 
 DEFBINOP (div, complex_matrix, sparse_complex_matrix)
 {
diff --git a/src/OPERATORS/op-cm-sm.cc b/src/OPERATORS/op-cm-sm.cc
--- a/src/OPERATORS/op-cm-sm.cc
+++ b/src/OPERATORS/op-cm-sm.cc
@@ -43,15 +43,7 @@
 DEFBINOP_OP (add, complex_matrix, sparse_matrix, +)
 DEFBINOP_OP (sub, complex_matrix, sparse_matrix, -)
 
-DEFBINOP (mul, complex_matrix, sparse_matrix)
-{
-  CAST_BINOP_ARGS (const octave_complex_matrix&, 
-		   const octave_sparse_matrix&);
-  
-  Matrix tmp (v2.matrix_value ());
-
-  return octave_value (v1.complex_matrix_value() * tmp);
-}
+DEFBINOP_OP (mul, complex_matrix, sparse_matrix, *)
 
 DEFBINOP (div, complex_matrix, sparse_matrix)
 {
diff --git a/src/OPERATORS/op-m-scm.cc b/src/OPERATORS/op-m-scm.cc
--- a/src/OPERATORS/op-m-scm.cc
+++ b/src/OPERATORS/op-m-scm.cc
@@ -44,15 +44,7 @@
 DEFBINOP_OP (add, matrix, sparse_complex_matrix, +)
 DEFBINOP_OP (sub, matrix, sparse_complex_matrix, -)
 
-DEFBINOP (mul, matrix, sparse_complex_matrix)
-{
-  CAST_BINOP_ARGS (const octave_matrix&, 
-		   const octave_sparse_complex_matrix&);
-  
-  ComplexMatrix tmp (v2.complex_matrix_value ());
-
-  return octave_value ( v1.matrix_value() * tmp);
-}
+DEFBINOP_OP (mul, matrix, sparse_complex_matrix, *)
 
 DEFBINOP (div, matrix, sparse_complex_matrix)
 {
diff --git a/src/OPERATORS/op-m-sm.cc b/src/OPERATORS/op-m-sm.cc
--- a/src/OPERATORS/op-m-sm.cc
+++ b/src/OPERATORS/op-m-sm.cc
@@ -43,14 +43,7 @@
 DEFBINOP_OP (add, matrix, sparse_matrix, +)
 DEFBINOP_OP (sub, matrix, sparse_matrix, -)
 
-DEFBINOP (mul, matrix, sparse_matrix)
-{
-  CAST_BINOP_ARGS (const octave_matrix&, const octave_sparse_matrix&);
-  
-  Matrix tmp (v2.matrix_value ());
-
-  return octave_value ( v1.matrix_value() * tmp);
-}
+DEFBINOP_OP (mul, matrix, sparse_matrix, *)
 
 DEFBINOP (div, matrix, sparse_matrix)
 {
diff --git a/src/OPERATORS/op-scm-cm.cc b/src/OPERATORS/op-scm-cm.cc
--- a/src/OPERATORS/op-scm-cm.cc
+++ b/src/OPERATORS/op-scm-cm.cc
@@ -43,15 +43,7 @@
 DEFBINOP_OP (add, sparse_complex_matrix, complex_matrix, +)
 DEFBINOP_OP (sub, sparse_complex_matrix, complex_matrix, -)
 
-DEFBINOP (mul, sparse_complex_matrix, complex_matrix)
-{
-  CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, 
-		   const octave_complex_matrix&);
-  
-  ComplexMatrix tmp (v1.complex_matrix_value ());
-
-  return octave_value ( tmp * v2.complex_matrix_value());
-}
+DEFBINOP_OP (mul, sparse_complex_matrix, complex_matrix, *)
 
 DEFBINOP (div, sparse_complex_matrix, complex_matrix)
 {
diff --git a/src/OPERATORS/op-scm-m.cc b/src/OPERATORS/op-scm-m.cc
--- a/src/OPERATORS/op-scm-m.cc
+++ b/src/OPERATORS/op-scm-m.cc
@@ -44,15 +44,7 @@
 DEFBINOP_OP (add, sparse_complex_matrix, matrix, +)
 DEFBINOP_OP (sub, sparse_complex_matrix, matrix, -)
 
-DEFBINOP (mul, sparse_complex_matrix, matrix)
-{
-  CAST_BINOP_ARGS (const octave_sparse_complex_matrix&, 
-		   const octave_matrix&);
-  
-  ComplexMatrix tmp (v1.complex_matrix_value ());
-
-  return octave_value ( tmp * v2.matrix_value());
-}
+DEFBINOP_OP (mul, sparse_complex_matrix, matrix, *)
 
 DEFBINOP (div, sparse_complex_matrix, matrix)
 {
diff --git a/src/OPERATORS/op-sm-cm.cc b/src/OPERATORS/op-sm-cm.cc
--- a/src/OPERATORS/op-sm-cm.cc
+++ b/src/OPERATORS/op-sm-cm.cc
@@ -43,15 +43,7 @@
 DEFBINOP_OP (add, sparse_matrix, complex_matrix, +)
 DEFBINOP_OP (sub, sparse_matrix, complex_matrix, -)
 
-DEFBINOP (mul, sparse_matrix, complex_matrix)
-{
-  CAST_BINOP_ARGS (const octave_sparse_matrix&, 
-		   const octave_complex_matrix&);
-  
-  Matrix tmp (v1.matrix_value ());
-
-  return octave_value ( tmp * v2.complex_matrix_value());
-}
+DEFBINOP_OP (mul, sparse_matrix, complex_matrix, *)
 
 DEFBINOP (div, sparse_matrix, complex_matrix)
 {
diff --git a/src/OPERATORS/op-sm-m.cc b/src/OPERATORS/op-sm-m.cc
--- a/src/OPERATORS/op-sm-m.cc
+++ b/src/OPERATORS/op-sm-m.cc
@@ -43,14 +43,7 @@
 DEFBINOP_OP (add, sparse_matrix, matrix, +)
 DEFBINOP_OP (sub, sparse_matrix, matrix, -)
 
-DEFBINOP (mul, sparse_matrix, matrix)
-{
-  CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_matrix&);
-  
-  Matrix tmp (v1.matrix_value ());
-
-  return octave_value ( tmp * v2.matrix_value());
-}
+DEFBINOP_OP (mul, sparse_matrix, matrix, *)
 
 DEFBINOP (div, sparse_matrix, matrix)
 {