1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin lu_solve.cpp}
LuSolve With Complex Arguments: Example and Test
################################################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end lu_solve.cpp}
*/
// BEGIN C++
# include <cppad/utility/lu_solve.hpp> // for CppAD::LuSolve
# include <cppad/utility/near_equal.hpp> // for CppAD::NearEqual
# include <cppad/utility/vector.hpp> // for CppAD::vector
# include <complex> // for std::complex
typedef std::complex<double> Complex; // define the Complex type
bool LuSolve(void)
{ bool ok = true;
using namespace CppAD;
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;
}
// END C++
|
