# HG changeset patch
# User Jaroslav Hajek <highegg@gmail.com>
# Date 1235397296 -3600
# Node ID c690e37725839344fb2579ccf5dc68ef9f04a1d3
# Parent  fcba62cc454939ea068fbfc8a581640d85bb66e3
support diagonal matrices in pinv

diff -r fcba62cc4549 -r c690e3772583 liboctave/CDiagMatrix.cc
--- a/liboctave/CDiagMatrix.cc	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/CDiagMatrix.cc	Mon Feb 23 14:54:56 2009 +0100
@@ -387,6 +387,26 @@
   return retval;
 }
 
+ComplexDiagMatrix
+ComplexDiagMatrix::pseudo_inverse (void) const
+{
+  octave_idx_type r = rows ();
+  octave_idx_type c = cols ();
+  octave_idx_type len = length ();
+
+  ComplexDiagMatrix retval (c, r);
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      if (elem (i, i) != 0.0)
+        retval.elem (i, i) = 1.0 / elem (i, i);
+      else
+        retval.elem (i, i) = 0.0;
+    }
+
+  return retval;
+}
+
 bool
 ComplexDiagMatrix::all_elements_are_real (void) const
 {
diff -r fcba62cc4549 -r c690e3772583 liboctave/CDiagMatrix.h
--- a/liboctave/CDiagMatrix.h	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/CDiagMatrix.h	Mon Feb 23 14:54:56 2009 +0100
@@ -112,6 +112,7 @@
 
   ComplexDiagMatrix inverse (octave_idx_type& info) const;
   ComplexDiagMatrix inverse (void) const;
+  ComplexDiagMatrix pseudo_inverse (void) const;
 
   bool all_elements_are_real (void) const;
 
diff -r fcba62cc4549 -r c690e3772583 liboctave/ChangeLog
--- a/liboctave/ChangeLog	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/ChangeLog	Mon Feb 23 14:54:56 2009 +0100
@@ -1,3 +1,14 @@
+2009-02-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* dDiagMatrix.cc (DiagMatrix::pseudo_inverse): New method.
+	* dDiagMatrix.h: Declare it.
+	* fDiagMatrix.cc (FloatDiagMatrix::pseudo_inverse): New method.
+	* fDiagMatrix.h: Declare it.
+	* CDiagMatrix.cc (ComplexDiagMatrix::pseudo_inverse): New method.
+	* CDiagMatrix.h: Declare it.
+	* fCDiagMatrix.cc (FloatComplexDiagMatrix::pseudo_inverse): New method.
+	* fCDiagMatrix.h: Declare it.
+
 2009-02-20  Jaroslav Hajek  <highegg@gmail.com>
 
 	* oct-sort.h (octave_sort<T>::MergeState::MergeState): New
diff -r fcba62cc4549 -r c690e3772583 liboctave/dDiagMatrix.cc
--- a/liboctave/dDiagMatrix.cc	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/dDiagMatrix.cc	Mon Feb 23 14:54:56 2009 +0100
@@ -303,6 +303,26 @@
   return retval;
 }
 
+DiagMatrix
+DiagMatrix::pseudo_inverse (void) const
+{
+  octave_idx_type r = rows ();
+  octave_idx_type c = cols ();
+  octave_idx_type len = length ();
+
+  DiagMatrix retval (c, r);
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      if (elem (i, i) != 0.0)
+        retval.elem (i, i) = 1.0 / elem (i, i);
+      else
+        retval.elem (i, i) = 0.0;
+    }
+
+  return retval;
+}
+
 // diagonal matrix by diagonal matrix -> diagonal matrix operations
 
 // diagonal matrix by diagonal matrix -> diagonal matrix operations
diff -r fcba62cc4549 -r c690e3772583 liboctave/dDiagMatrix.h
--- a/liboctave/dDiagMatrix.h	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/dDiagMatrix.h	Mon Feb 23 14:54:56 2009 +0100
@@ -94,6 +94,7 @@
 
   DiagMatrix inverse (void) const;
   DiagMatrix inverse (octave_idx_type& info) const;
+  DiagMatrix pseudo_inverse (void) const;
 
   // other operations
 
diff -r fcba62cc4549 -r c690e3772583 liboctave/fCDiagMatrix.cc
--- a/liboctave/fCDiagMatrix.cc	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/fCDiagMatrix.cc	Mon Feb 23 14:54:56 2009 +0100
@@ -387,6 +387,26 @@
   return retval;
 }
 
+FloatComplexDiagMatrix
+FloatComplexDiagMatrix::pseudo_inverse (void) const
+{
+  octave_idx_type r = rows ();
+  octave_idx_type c = cols ();
+  octave_idx_type len = length ();
+
+  FloatComplexDiagMatrix retval (c, r);
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      if (elem (i, i) != 0.0f)
+        retval.elem (i, i) = 1.0f / elem (i, i);
+      else
+        retval.elem (i, i) = 0.0f;
+    }
+
+  return retval;
+}
+
 bool
 FloatComplexDiagMatrix::all_elements_are_real (void) const
 {
diff -r fcba62cc4549 -r c690e3772583 liboctave/fCDiagMatrix.h
--- a/liboctave/fCDiagMatrix.h	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/fCDiagMatrix.h	Mon Feb 23 14:54:56 2009 +0100
@@ -112,6 +112,7 @@
 
   FloatComplexDiagMatrix inverse (octave_idx_type& info) const;
   FloatComplexDiagMatrix inverse (void) const;
+  FloatComplexDiagMatrix pseudo_inverse (void) const;
 
   bool all_elements_are_real (void) const;
 
diff -r fcba62cc4549 -r c690e3772583 liboctave/fDiagMatrix.cc
--- a/liboctave/fDiagMatrix.cc	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/fDiagMatrix.cc	Mon Feb 23 14:54:56 2009 +0100
@@ -303,6 +303,26 @@
   return retval;
 }
 
+FloatDiagMatrix
+FloatDiagMatrix::pseudo_inverse (void) const
+{
+  octave_idx_type r = rows ();
+  octave_idx_type c = cols ();
+  octave_idx_type len = length ();
+
+  FloatDiagMatrix retval (c, r);
+
+  for (octave_idx_type i = 0; i < len; i++)
+    {
+      if (elem (i, i) != 0.0f)
+        retval.elem (i, i) = 1.0f / elem (i, i);
+      else
+        retval.elem (i, i) = 0.0f;
+    }
+
+  return retval;
+}
+
 // diagonal matrix by diagonal matrix -> diagonal matrix operations
 
 // diagonal matrix by diagonal matrix -> diagonal matrix operations
diff -r fcba62cc4549 -r c690e3772583 liboctave/fDiagMatrix.h
--- a/liboctave/fDiagMatrix.h	Mon Feb 23 13:55:44 2009 +0100
+++ b/liboctave/fDiagMatrix.h	Mon Feb 23 14:54:56 2009 +0100
@@ -94,6 +94,7 @@
 
   FloatDiagMatrix inverse (void) const;
   FloatDiagMatrix inverse (octave_idx_type& info) const;
+  FloatDiagMatrix pseudo_inverse (void) const;
 
   // other operations
 
diff -r fcba62cc4549 -r c690e3772583 src/ChangeLog
--- a/src/ChangeLog	Mon Feb 23 13:55:44 2009 +0100
+++ b/src/ChangeLog	Mon Feb 23 14:54:56 2009 +0100
@@ -1,3 +1,7 @@
+2009-02-23  Jaroslav Hajek  <highegg@gmail.com>
+
+	* DLD-FUNCTIONS/pinv.cc: Support diagonal and permutation matrices.
+
 2008-02-21  David Bateman  <dbateman@free.fr>
 
 	* OPERATORS/op-cs-scm.cc, OPERATORS/op-cs-sm.cc, OPERATORS/op-s-scm.cc,
diff -r fcba62cc4549 -r c690e3772583 src/DLD-FUNCTIONS/pinv.cc
--- a/src/DLD-FUNCTIONS/pinv.cc	Mon Feb 23 13:55:44 2009 +0100
+++ b/src/DLD-FUNCTIONS/pinv.cc	Mon Feb 23 14:54:56 2009 +0100
@@ -29,6 +29,13 @@
 #include "gripes.h"
 #include "oct-obj.h"
 #include "utils.h"
+#include "ops.h"
+#include "ov-re-diag.h"
+#include "ov-cx-diag.h"
+#include "ov-flt-re-diag.h"
+#include "ov-flt-cx-diag.h"
+#include "ov-perm.h"
+#include "ov-flt-perm.h"
 
 DEFUN_DLD (pinv, args, ,
   "-*- texinfo -*-\n\
@@ -65,7 +72,56 @@
   else if (arg_is_empty > 0)
     return octave_value (Matrix ());
 
-  if (arg.is_single_type ())
+  bool isfloat = arg.is_single_type ();
+
+  if (arg.is_diag_matrix ())
+    {
+      if (nargin == 2)
+        warning ("pinv: tol is ignored for diagonal matrices");
+
+      const octave_base_value& a = arg.get_rep ();
+      if (arg.is_complex_type ())
+        {
+          if (isfloat)
+            {
+              CAST_CONV_ARG (const octave_float_complex_diag_matrix&);
+              retval = v.float_complex_diag_matrix_value ().pseudo_inverse ();
+            }
+          else
+            {
+              CAST_CONV_ARG (const octave_complex_diag_matrix&);
+              retval = v.complex_diag_matrix_value ().pseudo_inverse ();
+            }
+        }
+      else
+        {
+          if (isfloat)
+            {
+              CAST_CONV_ARG (const octave_float_diag_matrix&);
+              retval = v.float_diag_matrix_value ().pseudo_inverse ();
+            }
+          else
+            {
+              CAST_CONV_ARG (const octave_diag_matrix&);
+              retval = v.diag_matrix_value ().pseudo_inverse ();
+            }
+        }
+    }
+  else if (arg.is_perm_matrix ())
+    {
+      const octave_base_value& a = arg.get_rep ();
+      if (isfloat)
+        {
+          CAST_CONV_ARG (const octave_float_perm_matrix&);
+          retval = v.perm_matrix_value ().inverse ();
+        }
+      else
+        {
+          CAST_CONV_ARG (const octave_perm_matrix&);
+          retval = v.perm_matrix_value ().inverse ();
+        }
+    }
+  else if (isfloat)
     {
       float tol = 0.0;
       if (nargin == 2)