Mercurial > octave-nkf

changeset 9698:7c6d5d8c8d37

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix diag*diag multiplication
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 06 Oct 2009 08:43:01 +0200
parents 51c17bd18563
children d896036d975b
files liboctave/CDiagMatrix.cc liboctave/ChangeLog liboctave/DiagArray2.h liboctave/PermMatrix.h liboctave/dDiagMatrix.cc liboctave/fCDiagMatrix.cc liboctave/fDiagMatrix.cc
diffstat 7 files changed, 42 insertions(+), 61 deletions(-) [+]
line wrap: on
line diff
--- 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;
 }
--- 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.
--- 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 (); }
--- 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 (); }
 
--- 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;
 }
--- 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;
 }
--- 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;
 }