00001 00033 #include <itpp/base/det.h> 00034 #include <itpp/base/lu.h> 00035 00036 00037 namespace itpp { 00038 00039 00040 /* Determinant of square matrix. 00041 Calculate determinant of inmatrix (Uses LU-factorisation) 00042 (See Theorem 3.2.1 p. 97 in Golub & van Loan, "Matrix Computations"). 00043 00044 det(X) = det(P')*det(L)*det(U) = det(P')*1*prod(diag(U)) 00045 */ 00046 double det(const mat &X) 00047 { 00048 it_assert1(X.rows()==X.cols(),"det : Only square matrices"); 00049 00050 mat L, U; 00051 ivec p; 00052 double s=1.0; 00053 int i; 00054 00055 lu(X,L,U,p); // calculate LU-factorisation 00056 00057 double temp=U(0,0); 00058 for (i=1;i<X.rows();i++) { 00059 temp*=U(i,i); 00060 } 00061 00062 // Calculate det(P'). Equal to (-1)^(no row changes) 00063 for (i=0; i<p.size(); i++) 00064 if (i != p(i)) 00065 s *=-1.0; 00066 00067 return temp*s; 00068 } 00069 00070 00071 /* Determinant of complex square matrix. 00072 Calculate determinant of inmatrix (Uses LU-factorisation) 00073 (See Theorem 3.2.1 p. 97 in Golub & van Loan, "Matrix Computations"). 00074 00075 det(X) = det(P')*det(L)*det(U) = det(P')*1*prod(diag(U)) 00076 00077 Needs LU-factorization of complex matrices (LAPACK) 00078 */ 00079 std::complex<double> det(const cmat &X) 00080 { 00081 it_assert1(X.rows()==X.cols(),"det : Only square matrices"); 00082 00083 int i; 00084 cmat L, U; 00085 ivec p; 00086 double s=1.0; 00087 00088 lu(X,L,U,p); // calculate LU-factorisation 00089 00090 std::complex<double> temp=U(0,0); 00091 for (i=1;i<X.rows();i++) { 00092 temp*=U(i,i); 00093 } 00094 00095 // Calculate det(P'). Equal to (-1)^(no row changes) 00096 for (i=0; i<p.size(); i++) 00097 if (i != p(i)) 00098 s *=-1.0; 00099 00100 return temp*s; 00101 } 00102 00103 00104 } // namespace itpp
Generated on Fri Jan 11 08:51:36 2008 for IT++ by Doxygen 1.3.9.1