Mercurial > hg > octave-nkf
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 ();