Prev Next Index-> contents reference index search external Up-> CppAD utility LuDetAndSolve LuSolve lu_solve.cpp utility-> ErrorHandler NearEqual speed_test SpeedTest time_test test_boolofvoid NumericType CheckNumericType SimpleVector CheckSimpleVector nan pow_int Poly LuDetAndSolve RombergOne RombergMul Runge45 Rosen34 OdeErrControl OdeGear OdeGearControl CppAD_vector thread_alloc index_sort to_string set_union sparse_rc sparse_rcv LuDetAndSolve-> LuSolve LuFactor LuInvert LuSolve-> lu_solve.cpp lu_solve.hpp lu_solve.cpp Headings

$\newcommand{\W}[1]{ \; #1 \; } \newcommand{\R}[1]{ {\rm #1} } \newcommand{\B}[1]{ {\bf #1} } \newcommand{\D}[2]{ \frac{\partial #1}{\partial #2} } \newcommand{\DD}[3]{ \frac{\partial^2 #1}{\partial #2 \partial #3} } \newcommand{\Dpow}[2]{ \frac{\partial^{#1}}{\partial {#2}^{#1}} } \newcommand{\dpow}[2]{ \frac{ {\rm d}^{#1}}{{\rm d}\, {#2}^{#1}} }$
LuSolve With Complex Arguments: Example and Test

# include <complex>               // for std::complex

typedef std::complex<double> Complex;    // define the Complex type
bool LuSolve(void)
{     bool  ok = true;

size_t   n = 3;           // number rows in A and B
size_t   m = 2;           // number columns in B, X and S

// A is an n by n matrix, B, X, and S are n by m matrices
CppAD::vector<Complex> A(n * n), B(n * m), X(n * m) , S(n * m);

Complex  logdet;          // log of determinant of A
int      signdet;         // zero if A is singular
Complex  det;             // determinant of A
size_t   i, j, k;         // some temporary indices

// set A equal to the n by n Hilbert Matrix
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
A[i * n + j] = 1. / (double) (i + j + 1);

// set S to the solution of the equation we will solve
for(j = 0; j < n; j++)
for(k = 0; k < m; k++)
S[ j * m + k ] = Complex(double(j), double(j + k));

// set B = A * S
size_t ik;
Complex sum;
for(k = 0; k < m; k++)
{     for(i = 0; i < n; i++)
{     sum = 0.;
for(j = 0; j < n; j++)
sum += A[i * n + j] * S[j * m + k];
B[i * m + k] = sum;
}
}

// solve the equation A * X = B and compute determinant of A
signdet = CppAD::LuSolve(n, m, A, B, X, logdet);
det     = Complex( signdet ) * exp( logdet );

double cond  = 4.62963e-4;       // condition number of A when n = 3
double determinant = 1. / 2160.; // determinant of A when n = 3
double delta = 1e-14 / cond;     // accuracy expected in X

// check determinant
ok &= CppAD::NearEqual(det, determinant, delta, delta);

// check solution
for(ik = 0; ik < n * m; ik++)
ok &= CppAD::NearEqual(X[ik], S[ik], delta, delta);

return ok;
}

Input File: example/utility/lu_solve.cpp