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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
// 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_invert.cpp}
LuInvert: Example and Test
##########################
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end lu_invert.cpp}
*/
// BEGIN C++
# include <cstdlib> // for rand function
# include <cppad/utility/lu_invert.hpp> // for CppAD::LuInvert
# include <cppad/utility/near_equal.hpp> // for CppAD::NearEqual
# include <cppad/utility/vector.hpp> // for CppAD::vector
bool LuInvert(void)
{ bool ok = true;
# ifndef _MSC_VER
using std::rand;
using std::srand;
# endif
double eps200 = 200.0 * std::numeric_limits<double>::epsilon();
size_t n = 7; // number rows in A
size_t m = 3; // number columns in B
double rand_max = double(RAND_MAX); // maximum rand value
double sum; // element of L * U
size_t i, j, k; // temporary indices
// dimension matrices
CppAD::vector<double>
A(n*n), X(n*m), B(n*m), LU(n*n), L(n*n), U(n*n);
// seed the random number generator
srand(123);
// pivot vectors
CppAD::vector<size_t> ip(n);
CppAD::vector<size_t> jp(n);
// set pivot vectors
for(i = 0; i < n; i++)
{ ip[i] = (i + 2) % n; // ip = 2 , 3, ... , n-1, 0, 1
jp[i] = (n + 2 - i) % n; // jp = 2 , 1, n-1, n-2, ... , 3
}
// chose L, a random lower triangular matrix
for(i = 0; i < n; i++)
{ for(j = 0; j <= i; j++)
L [i * n + j] = rand() / rand_max;
for(j = i+1; j < n; j++)
L [i * n + j] = 0.;
}
// chose U, a random upper triangular matrix with ones on diagonal
for(i = 0; i < n; i++)
{ for(j = 0; j < i; j++)
U [i * n + j] = 0.;
U[ i * n + i ] = 1.;
for(j = i+1; j < n; j++)
U [i * n + j] = rand() / rand_max;
}
// chose X, a random matrix
for(i = 0; i < n; i++)
{ for(k = 0; k < m; k++)
X[i * m + k] = rand() / rand_max;
}
// set LU to a permuted combination of both L and U
for(i = 0; i < n; i++)
{ for(j = 0; j <= i; j++)
LU [ ip[i] * n + jp[j] ] = L[i * n + j];
for(j = i+1; j < n; j++)
LU [ ip[i] * n + jp[j] ] = U[i * n + j];
}
// set A to a permuted version of L * U
for(i = 0; i < n; i++)
{ for(j = 0; j < n; j++)
{ // compute (i,j) entry in permuted matrix
sum = 0.;
for(k = 0; k < n; k++)
sum += L[i * n + k] * U[k * n + j];
A[ ip[i] * n + jp[j] ] = sum;
}
}
// set B to A * X
for(i = 0; i < n; i++)
{ for(k = 0; k < m; k++)
{ // compute (i,k) entry of B
sum = 0.;
for(j = 0; j < n; j++)
sum += A[i * n + j] * X[j * m + k];
B[i * m + k] = sum;
}
}
// solve for X
CppAD::LuInvert(ip, jp, LU, B);
// check result
for(i = 0; i < n; i++)
{ for(k = 0; k < m; k++)
{ ok &= CppAD::NearEqual(
X[i * m + k], B[i * m + k], eps200, eps200
);
}
}
return ok;
}
// END C++
|
