If A is a square or symmetric matrix use
CroutMatrix X = A; // carries out LU decomposition
Matrix AP = X.i()*P; Matrix AQ = X.i()*Q;
LogAndSign ld = X.LogDeterminant();
rather than
Matrix AP = A.i()*P; Matrix AQ = A.i()*Q;
LogAndSign ld = A.LogDeterminant();
since each operation will repeat the LU decomposition.
If A is a BandMatrix or a SymmetricBandMatrix begin with
BandLUMatrix X = A; // carries out LU decomposition
A CroutMatrix or a BandLUMatrix can't be manipulated or
copied. Use
references as an alternative to copying.
Alternatively use
LinearEquationSolver X = A;
This will choose the most appropiate decomposition of A. That is, the
band form if A is banded; the Crout decomposition if A is
square or symmetric and no decomposition if A is triangular or
diagonal. If you want to use the LinearEquationSolver #include
newmatap.h.