# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1254811381 -7200
# Node ID 7c6d5d8c8d375b85e2a0f60272276a1205915852
# Parent  51c17bd185632c6ba8a60a150ddeeab3c6ce42d7
fix diag*diag multiplication

diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/CDiagMatrix.cc
--- a/liboctave/CDiagMatrix.cc	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/CDiagMatrix.cc	Tue Oct 06 08:43:01 2009 +0200
@@ -460,25 +460,16 @@
   octave_idx_type b_nc = b.cols ();
 
   if (a_nc != b_nr)
-    {
-      gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
-      return ComplexDiagMatrix ();
-    }
-
-  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
-    return ComplexDiagMatrix (a_nr, a_nc, 0.0);
+    gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
 
   ComplexDiagMatrix c (a_nr, b_nc);
 
-  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
+  octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
-  for (octave_idx_type i = 0; i < len; i++)
-    {
-      Complex a_element = a.elem (i, i);
-      double b_element = b.elem (i, i);
-
-      c.elem (i, i) = a_element * b_element;
-    }
+  for (octave_idx_type i = 0; i < lenm; i++)
+    c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+  for (octave_idx_type i = lenm; i < len; i++)
+    c.dgxelem (i) = 0.0;
 
   return c;
 }
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/ChangeLog	Tue Oct 06 08:43:01 2009 +0200
@@ -1,3 +1,18 @@
+2009-10-06  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dDiagMatrix.cc (operator *(const DiagMatrix&, const DiagMatrix&)):
+	Rewrite.
+	* fDiagMatrix.cc (operator *(const FloatDiagMatrix&, const FloatDiagMatrix&)):
+	Rewrite.
+	* CDiagMatrix.cc (operator *(const ComplexDiagMatrix&, const ComplexDiagMatrix&)):
+	Rewrite.
+	* fCDiagMatrix.cc (operator *(const FloatComplexDiagMatrix&, const
+	FloatComplexDiagMatrix&)):
+	Rewrite.
+	* DiagArray2.h (DiagArray2::diag_length): New method.
+	* PermMatrix.h (PermMatrix::length): Make consistent with
+	DiagArray2::length.
+
 2009-10-05  Jaroslav Hajek  <highegg@gmail.com>
 
 	* base-lu.cc (base_lu::unpack): Unpack getp as well.
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/DiagArray2.h
--- a/liboctave/DiagArray2.h	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/DiagArray2.h	Tue Oct 06 08:43:01 2009 +0200
@@ -143,6 +143,7 @@
   octave_idx_type cols (void) const { return dim2 (); }
   octave_idx_type columns (void) const { return dim2 (); }
 
+  octave_idx_type diag_length (void) const { return Array<T>::length (); }
   // FIXME: a dangerous ambiguity?
   octave_idx_type length (void) const { return Array<T>::length (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/PermMatrix.h
--- a/liboctave/PermMatrix.h	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/PermMatrix.h	Tue Oct 06 08:43:01 2009 +0200
@@ -57,8 +57,9 @@
 
   octave_idx_type perm_length (void) const 
     { return Array<octave_idx_type>::length (); }
+  // FIXME: a dangerous ambiguity?
   octave_idx_type length (void) const 
-    { return dim1 () * dim2 (); }
+    { return perm_length (); }
   octave_idx_type nelem (void) const { return dim1 () * dim2 (); }
   octave_idx_type numel (void) const { return nelem (); }
 
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/dDiagMatrix.cc
--- a/liboctave/dDiagMatrix.cc	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/dDiagMatrix.cc	Tue Oct 06 08:43:01 2009 +0200
@@ -338,25 +338,16 @@
   octave_idx_type b_nc = b.cols ();
 
   if (a_nc != b_nr)
-    {
-      gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
-      return DiagMatrix ();
-    }
-
-  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
-    return DiagMatrix (a_nr, a_nc, 0.0);
+    gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
 
   DiagMatrix c (a_nr, b_nc);
 
-  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
+  octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
-  for (octave_idx_type i = 0; i < len; i++)
-    {
-      double a_element = a.elem (i, i);
-      double b_element = b.elem (i, i);
-
-      c.elem (i, i) = a_element * b_element;
-    }
+  for (octave_idx_type i = 0; i < lenm; i++)
+    c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+  for (octave_idx_type i = lenm; i < len; i++)
+    c.dgxelem (i) = 0.0;
 
   return c;
 }
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/fCDiagMatrix.cc
--- a/liboctave/fCDiagMatrix.cc	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/fCDiagMatrix.cc	Tue Oct 06 08:43:01 2009 +0200
@@ -460,25 +460,16 @@
   octave_idx_type b_nc = b.cols ();
 
   if (a_nc != b_nr)
-    {
-      gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
-      return FloatComplexDiagMatrix ();
-    }
-
-  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
-    return FloatComplexDiagMatrix (a_nr, a_nc, 0.0);
+    gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
 
   FloatComplexDiagMatrix c (a_nr, b_nc);
 
-  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
+  octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
-  for (octave_idx_type i = 0; i < len; i++)
-    {
-      FloatComplex a_element = a.elem (i, i);
-      float b_element = b.elem (i, i);
-
-      c.elem (i, i) = a_element * b_element;
-    }
+  for (octave_idx_type i = 0; i < lenm; i++)
+    c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+  for (octave_idx_type i = lenm; i < len; i++)
+    c.dgxelem (i) = 0.0f;
 
   return c;
 }
diff -r 51c17bd18563 -r 7c6d5d8c8d37 liboctave/fDiagMatrix.cc
--- a/liboctave/fDiagMatrix.cc	Tue Oct 06 08:22:56 2009 +0200
+++ b/liboctave/fDiagMatrix.cc	Tue Oct 06 08:43:01 2009 +0200
@@ -338,25 +338,16 @@
   octave_idx_type b_nc = b.cols ();
 
   if (a_nc != b_nr)
-    {
-      gripe_nonconformant ("operaotr *", a_nr, a_nc, b_nr, b_nc);
-      return FloatDiagMatrix ();
-    }
-
-  if (a_nr == 0 || a_nc == 0 || b_nc == 0)
-    return FloatDiagMatrix (a_nr, a_nc, 0.0);
+    gripe_nonconformant ("operator *", a_nr, a_nc, b_nr, b_nc);
 
   FloatDiagMatrix c (a_nr, b_nc);
 
-  octave_idx_type len = a_nr < b_nc ? a_nr : b_nc;
+  octave_idx_type len = c.length (), lenm = len < a_nc ? len : a_nc;
 
-  for (octave_idx_type i = 0; i < len; i++)
-    {
-      float a_element = a.elem (i, i);
-      float b_element = b.elem (i, i);
-
-      c.elem (i, i) = a_element * b_element;
-    }
+  for (octave_idx_type i = 0; i < lenm; i++)
+    c.dgxelem (i) = a.dgelem (i) * b.dgelem (i);
+  for (octave_idx_type i = lenm; i < len; i++)
+    c.dgxelem (i) = 0.0f;
 
   return c;
 }