diff liboctave/CMatrix.cc @ 740:d8295febb0df

[project @ 1994-09-30 14:42:37 by jwe]
author jwe
date Fri, 30 Sep 1994 14:42:37 +0000
parents 01da6806197b
children 714fd17fca28
line wrap: on
line diff
--- 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 ();