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