1 //
2 // cgsolve_example.c
3 //
4 // Solve linear system of equations A*x = b using the conjugate-
5 // gradient method where A is a symmetric positive-definite matrix.
6 // Compare speed to matrixf_linsolve() for same system.
7 //
8
9 #include <stdio.h>
10 #include <stdlib.h>
11 #include <string.h>
12 #include <math.h>
13 #include "liquid.h"
14
main()15 int main() {
16 // options
17 unsigned int n = 8;
18
19 unsigned int i;
20
21 // allocate memory for arrays
22 float A[n*n];
23 float b[n];
24 float x[n];
25 float x_hat[n];
26
27 // generate symmetric positive-definite matrix by first generating
28 // lower triangular matrix L and computing A = L*L'
29 float L[n*n];
30 unsigned int j;
31 for (i=0; i<n; i++) {
32 for (j=0; j<n; j++) {
33 #if 0
34 // sparse matrix
35 if (j > i) matrix_access(L,n,n,i,j) = 0.0;
36 else if (j == i) matrix_access(L,n,n,i,j) = randnf();
37 else if ((rand()%4)==0) matrix_access(L,n,n,i,j) = randnf();
38 else matrix_access(L,n,n,i,j) = 0.0;
39 #else
40 // full matrix
41 matrix_access(L,n,n,i,j) = (j < i) ? 0.0 : randnf();
42 #endif
43 }
44 }
45 matrixf_mul_transpose(L, n, n, A);
46
47 // generate random solution
48 for (i=0; i<n; i++)
49 x[i] = randnf();
50
51 // compute b
52 matrixf_mul(A, n, n,
53 x, n, 1,
54 b, n, 1);
55
56 // solve symmetric positive-definite system of equations
57 matrixf_cgsolve(A, n, b, x_hat, NULL);
58 //matrixf_linsolve(A, n, b, x_hat, NULL);
59
60 // print results
61
62 printf("A:\n"); matrixf_print(A, n, n);
63 printf("b:\n"); matrixf_print(b, n, 1);
64 printf("x (original):\n"); matrixf_print(x, n, 1);
65 printf("x (estimate):\n"); matrixf_print(x_hat, n, 1);
66
67 // compute error norm
68 float e = 0.0;
69 for (i=0; i<n; i++)
70 e += (x[i] - x_hat[i])*(x[i] - x_hat[i]);
71 e = sqrt(e);
72 printf("error norm: %12.4e\n", e);
73
74 printf("done.\n");
75 return 0;
76 }
77
78