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