# HG changeset patch
# User jwe
# Date 780936157 0
# Node ID d8295febb0df5c737e5f1a07c422b6114dfba6df
# Parent  d452e6f52d129286851fd33e92355c955d423454
[project @ 1994-09-30 14:42:37 by jwe]

diff --git a/liboctave/CMatrix.cc b/liboctave/CMatrix.cc
--- a/liboctave/CMatrix.cc
+++ b/liboctave/CMatrix.cc
@@ -31,11 +31,13 @@
 
 #include <sys/types.h>
 #include <iostream.h>
+#include <float.h>
 
 #include <Complex.h>
 
 #include "mx-base.h"
 #include "CmplxDET.h"
+#include "CmplxSVD.h"
 #include "mx-inlines.cc"
 #include "lo-error.h"
 #include "f77-uscore.h"
@@ -957,6 +959,43 @@
 }
 
 ComplexMatrix
+ComplexMatrix::pseudo_inverse (double tol)
+{
+  ComplexSVD result (*this);
+
+  DiagMatrix S = result.singular_values ();
+  ComplexMatrix U = result.left_singular_matrix ();
+  ComplexMatrix V = result.right_singular_matrix ();
+
+  ColumnVector sigma = S.diag ();
+
+  int r = sigma.length () - 1;
+  int nr = rows ();
+  int nc = cols ();
+
+  if (tol <= 0.0)
+    {
+      if (nr > nc)
+	tol = nr * sigma.elem (0) * DBL_EPSILON;
+      else
+	tol = nc * sigma.elem (0) * DBL_EPSILON;
+    }
+
+  while (r >= 0 && sigma.elem (r) < tol)
+    r--;
+
+  if (r < 0)
+    return ComplexMatrix (nc, nr, 0.0);
+  else
+    {
+      ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
+      DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
+      ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
+      return Vr * D * Ur.hermitian ();
+    }
+}
+
+ComplexMatrix
 ComplexMatrix::fourier (void) const
 {
   int nr = rows ();
diff --git a/liboctave/CMatrix.h b/liboctave/CMatrix.h
--- a/liboctave/CMatrix.h
+++ b/liboctave/CMatrix.h
@@ -134,6 +134,8 @@
   ComplexMatrix inverse (int& info) const;
   ComplexMatrix inverse (int& info, double& rcond) const;
 
+  ComplexMatrix pseudo_inverse (double tol = 0.0);
+
   ComplexMatrix fourier (void) const;
   ComplexMatrix ifourier (void) const;
 
diff --git a/liboctave/dMatrix.cc b/liboctave/dMatrix.cc
--- a/liboctave/dMatrix.cc
+++ b/liboctave/dMatrix.cc
@@ -32,11 +32,13 @@
 #include <sys/types.h>
 #include <iostream.h>
 #include <stdio.h>
+#include <float.h>
 
 #include <Complex.h>
 
 #include "mx-base.h"
 #include "dbleDET.h"
+#include "dbleSVD.h"
 #include "mx-inlines.cc"
 #include "lo-error.h"
 #include "f77-uscore.h"
@@ -605,6 +607,43 @@
   return Matrix (tmp_data, nr, nc);
 }
 
+Matrix
+Matrix::pseudo_inverse (double tol)
+{
+  SVD result (*this);
+
+  DiagMatrix S = result.singular_values ();
+  Matrix U = result.left_singular_matrix ();
+  Matrix V = result.right_singular_matrix ();
+
+  ColumnVector sigma = S.diag ();
+
+  int r = sigma.length () - 1;
+  int nr = rows ();
+  int nc = cols ();
+
+  if (tol <= 0.0)
+    {
+      if (nr > nc)
+	tol = nr * sigma.elem (0) * DBL_EPSILON;
+      else
+	tol = nc * sigma.elem (0) * DBL_EPSILON;
+    }
+
+  while (r >= 0 && sigma.elem (r) < tol)
+    r--;
+
+  if (r < 0)
+    return Matrix (nc, nr, 0.0);
+  else
+    {
+      Matrix Ur = U.extract (0, 0, nr-1, r);
+      DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
+      Matrix Vr = V.extract (0, 0, nc-1, r);
+      return Vr * D * Ur.transpose ();
+    }
+}
+
 ComplexMatrix
 Matrix::fourier (void) const
 {
diff --git a/liboctave/dMatrix.h b/liboctave/dMatrix.h
--- a/liboctave/dMatrix.h
+++ b/liboctave/dMatrix.h
@@ -109,6 +109,8 @@
   Matrix inverse (int& info) const;
   Matrix inverse (int& info, double& rcond) const;
 
+  Matrix pseudo_inverse (double tol = 0.0);
+
   ComplexMatrix fourier (void) const;
   ComplexMatrix ifourier (void) const;