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
|
#include <iostream>
#include <vector>
using namespace std;
using matrix = vector< vector<double> >;
//======================================================================
void print( const char *s, const matrix &M )
{
cout << s << '\n';
for ( auto row : M )
{
for ( auto e : row ) cout << e << '\t';
cout << '\n';
}
}
//======================================================================
matrix matmulSquare( const matrix &A, const matrix &B )
{
int N = A.size();
matrix C( N, vector<double>( N, 0 ) );
for ( int i = 0; i < N; i++ )
{
for ( int j = 0; j < N; j++ )
{
for ( int k = 0; k < N; k++ ) C[i][j] += A[i][k] * B[k][j];
}
}
return C;
}
//======================================================================
matrix invertLT( const matrix &L )
{
int N = L.size();
matrix M( N, vector<double>( N, 0 ) );
for ( int i = 0; i < N; i++ )
{
if ( L[i][i] == 0.0 )
{
cout << "*** Singular matrix ***\n";
return M;
}
M[i][i] = 1.0 / L[i][i];
for ( int j = 0; j < i; j++ )
{
for ( int k = j; k < i; k++ ) M[i][j] += L[i][k] * M[k][j];
M[i][j] = -M[i][j] / L[i][i];
}
}
return M;
}
//======================================================================
int main()
{
matrix L = { { 1, 0, 0, 0 },
{ 4, 1, 0, 0 },
{ 2, 5, 1, 0 },
{ 6, 8, 4, 1 } };
matrix M = invertLT( L );
matrix I = matmulSquare( L, M );
print( "\nL:", L );
print( "\ninverse(L):", M );
print( "\nProduct:", I );
}
|
