Mercurial > hg > octave-nkf

comparison liboctave/dMatrix.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
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 <stdio.h> 34 #include <stdio.h>
35 #include <float.h>
35 36
36 #include <Complex.h> 37 #include <Complex.h>
37 38
38 #include "mx-base.h" 39 #include "mx-base.h"
39 #include "dbleDET.h" 40 #include "dbleDET.h"
41 #include "dbleSVD.h"
40 #include "mx-inlines.cc" 42 #include "mx-inlines.cc"
41 #include "lo-error.h" 43 #include "lo-error.h"
42 #include "f77-uscore.h" 44 #include "f77-uscore.h"
43 45
44 // Fortran functions we call. 46 // Fortran functions we call.
601 603
602 delete [] ipvt; 604 delete [] ipvt;
603 delete [] z; 605 delete [] z;
604 606
605 return Matrix (tmp_data, nr, nc); 607 return Matrix (tmp_data, nr, nc);
608 }
609
610 Matrix
611 Matrix::pseudo_inverse (double tol)
612 {
613 SVD result (*this);
614
615 DiagMatrix S = result.singular_values ();
616 Matrix U = result.left_singular_matrix ();
617 Matrix V = result.right_singular_matrix ();
618
619 ColumnVector sigma = S.diag ();
620
621 int r = sigma.length () - 1;
622 int nr = rows ();
623 int nc = cols ();
624
625 if (tol <= 0.0)
626 {
627 if (nr > nc)
628 tol = nr * sigma.elem (0) * DBL_EPSILON;
629 else
630 tol = nc * sigma.elem (0) * DBL_EPSILON;
631 }
632
633 while (r >= 0 && sigma.elem (r) < tol)
634 r--;
635
636 if (r < 0)
637 return Matrix (nc, nr, 0.0);
638 else
639 {
640 Matrix Ur = U.extract (0, 0, nr-1, r);
641 DiagMatrix D = DiagMatrix (sigma.extract (0, r)) . inverse ();
642 Matrix Vr = V.extract (0, 0, nc-1, r);
643 return Vr * D * Ur.transpose ();
644 }
606 } 645 }
607 646
608 ComplexMatrix 647 ComplexMatrix
609 Matrix::fourier (void) const 648 Matrix::fourier (void) const
610 { 649 {