-
Notifications
You must be signed in to change notification settings - Fork 0
/
LUtest.c
129 lines (122 loc) · 3.31 KB
/
LUtest.c
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
121
122
123
124
125
126
127
128
129
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "LUdecomp.h"
double A5_[5][5] = {
{-1, -2, -3, 4, 5},
{ 2, 3, -4, -5, -6},
{-3, 4, -5, 6, -7},
{-4, 5, 6, -7, 8},
{ 5, -6, 7, -8, 9}
};
double A10_[10][10] = {
{ 8, 0, 8,-10,-10, -8, 9, 2,-10, -5},
{ 3, 5, -7, -5, -3, 2, -6, 1, -6, 4},
{ 1, 3, 9, 2, -1, -9, 0, 10, -9, 10},
{ 3, 3, -9, -5, 1, -2, 0, -2, 0, -6},
{ 2, 5, 10, -8, -4, -7, 4, -8, -4, 7},
{ 7, 2, -8, 1, 6, -3, 2, -2, 10, 5},
{-2, 3, -7, 9, -4, 7,-10, 3, 5, -3},
{ 1, -1, 8, 0, -8, -6, 3, -5, -6, -2},
{-7, -5, 0, 6, 6, -7, 8, -6, -2, 8},
{-7, 8, -6, 6, -8, -1, 8, -8, 9, 7}
};
int main() {
const double *A5[5] = {
A5_[0], A5_[1], A5_[2], A5_[3], A5_[4],
};
const double b5[] = {1, 2, 3, 4, 5};
//
// Solution from wolfram-alpha
// http://goo.gl/Nwe6f5
//
printf("solving 5x5 system...\n");
LUdecomp *LU = LUdecompose(5, A5);
double x5[5];
LUsolve(LU, b5, x5);
double x5_soln[5] = {
263.0/12, 107.0/6, 61.0/20, 139.0/15, 92.0/15
};
for (int i = 0; i < 5; i++) {
const double err = x5[i] - x5_soln[i];
printf("x[i] = %11.7f (%11.7f, error=%0.7e)\n",
x5[i], x5_soln[i], err);
}
LUdestroy(LU);
//
// compute inverse of A10.
// Solve A*X = I where I is 10x10 identity matrix.
//
printf("computing 10x10 inverse...\n");
const double *A10[10] = {
A10_[0], A10_[1], A10_[2], A10_[3], A10_[4],
A10_[5], A10_[6], A10_[7], A10_[8], A10_[9],
};
LU = LUdecompose(10, A10);
double X10_[10][10];
double *X10[10] = { // will hold transpose of inverse
X10_[0], X10_[1], X10_[2], X10_[3], X10_[4],
X10_[5], X10_[6], X10_[7], X10_[8], X10_[9],
};
for (int i = 0; i < 10; i++) {
double b[10];
for (int j = 0; j < 10; j++)
b[j] = (i == j) ? 1.0 : 0.0;
LUsolve(LU, b, X10[i]);
}
printf("checking solution (A*A^-1 = I)...\n");
double I[10][10];
double maxError = 0.0;
for (int i = 0; i < 10; i++)
for (int j = 0; j < 10; j++) {
double sum = 0.0;
for (int k = 0; k < 10; k++)
sum += A10[i][k]*X10[j][k];
I[i][j] = sum;
const double error = fabs(sum - ((i == j) ? 1.0 : 0.0));
if (error > maxError)
maxError = error;
}
printf("max error = %0.12e\n", maxError);
/* XXX
for (int i = 0; i < 10; i++) {
for (int j = 0; j < 10; j++) {
printf("%7.3f ", I[i][j]);
}
printf("\n");
}
*/
LUdestroy(LU);
//
// Big NxN system
//
const int N = 600;
printf("solving %dx%d system...\n", N, N);
double **S = (double **) malloc(N*sizeof(double*));
for (int i = 0; i < N; i++) {
S[i] = (double *) malloc(N*sizeof(double));
const double f = i + 1;
const double s = sin(f);
for (int j = 0; j < N; j++)
S[i][j] = sin(s*exp(j));
}
LU = LUdecompose(N, (const double **) S);
double *B = (double *) malloc(N*sizeof(double));
for (int i = 0; i < N; i++)
B[i] = i;
double *X = (double *) malloc(N*sizeof(double));
LUsolve(LU, B, X);
printf("checking solution...\n");
maxError = 0.0;
for (int i = 0; i < N; i++) {
double sum = 0.0;
for (int j = 0; j < N; j++)
sum += S[i][j]*X[j];
const double error = fabs(B[i] - sum);
if (error > maxError)
maxError = error;
}
LUdestroy(LU);
printf("max error = %0.12e\n", maxError);
return 0;
}