Mercurial > hg > octave-lyh

comparison liboctave/CMatrix.cc @ 740:d8295febb0df

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1994-09-30 14:42:37 by jwe]
author jwe
date Fri, 30 Sep 1994 14:42:37 +0000
parents 01da6806197b
children 714fd17fca28
comparison
equal deleted inserted replaced
739:d452e6f52d12 740:d8295febb0df
29 #pragma implementation 29 #pragma implementation
30 #endif 30 #endif
31 31
32 #include <sys/types.h> 32 #include <sys/types.h>
33 #include <iostream.h> 33 #include <iostream.h>
34 #include <float.h>
34 35
35 #include <Complex.h> 36 #include <Complex.h>
36 37
37 #include "mx-base.h" 38 #include "mx-base.h"
38 #include "CmplxDET.h" 39 #include "CmplxDET.h"
40 #include "CmplxSVD.h"
39 #include "mx-inlines.cc" 41 #include "mx-inlines.cc"
40 #include "lo-error.h" 42 #include "lo-error.h"
41 #include "f77-uscore.h" 43 #include "f77-uscore.h"
42 44
43 // Fortran functions we call. 45 // Fortran functions we call.
952 954
953 delete [] ipvt; 955 delete [] ipvt;
954 delete [] z; 956 delete [] z;
955 957
956 return ComplexMatrix (tmp_data, nr, nc); 958 return ComplexMatrix (tmp_data, nr, nc);
959 }
960
961 ComplexMatrix
962 ComplexMatrix::pseudo_inverse (double tol)
963 {
964 ComplexSVD result (*this);
965
966 DiagMatrix S = result.singular_values ();
967 ComplexMatrix U = result.left_singular_matrix ();
968 ComplexMatrix V = result.right_singular_matrix ();
969
970 ColumnVector sigma = S.diag ();
971
972 int r = sigma.length () - 1;
973 int nr = rows ();
974 int nc = cols ();
975
976 if (tol <= 0.0)
977 {
978 if (nr > nc)
979 tol = nr * sigma.elem (0) * DBL_EPSILON;
980 else
981 tol = nc * sigma.elem (0) * DBL_EPSILON;
982 }
983
984 while (r >= 0 && sigma.elem (r) < tol)
985 r--;
986
987 if (r < 0)
988 return ComplexMatrix (nc, nr, 0.0);
989 else
990 {
991 ComplexMatrix Ur = U.extract (0, 0, nr-1, r);
992 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
993 ComplexMatrix Vr = V.extract (0, 0, nc-1, r);
994 return Vr * D * Ur.hermitian ();
995 }
957 } 996 }
958 997
959 ComplexMatrix 998 ComplexMatrix
960 ComplexMatrix::fourier (void) const 999 ComplexMatrix::fourier (void) const
961 { 1000 {