|
1 | 1 | #include "qdldl.h"
|
2 |
| -#include <stdlib.h> |
3 | 2 | #include <stdio.h>
|
| 3 | +#include <stdlib.h> |
4 | 4 |
|
5 |
| -void print_arrayi(const QDLDL_int* data, QDLDL_int n,char* varName); |
| 5 | +void print_arrayi(const QDLDL_int* data, QDLDL_int n, char* varName); |
6 | 6 | void print_arrayf(const QDLDL_float* data, QDLDL_int n, char* varName);
|
7 | 7 | void print_line(void);
|
8 | 8 |
|
9 |
| -const QDLDL_int An = 10; |
10 |
| -const QDLDL_int Ap[] = {0, 1, 2, 4, 5, 6, 8, 10, 12, 14, 17}; |
11 |
| -const QDLDL_int Ai[] = {0, 1, 1, 2, 3, 4, 1, 5, 0, 6, 3, 7, 6, 8, 1, 2, 9}; |
12 |
| -const QDLDL_float Ax[] = {1.0, 0.460641, -0.121189, 0.417928, 0.177828, 0.1, |
13 |
| - -0.0290058, -1.0, 0.350321, -0.441092, -0.0845395, |
14 |
| - -0.316228, 0.178663, -0.299077, 0.182452, -1.56506, -0.1}; |
15 |
| -const QDLDL_float b[] = {1,2,3,4,5,6,7,8,9,10}; |
| 9 | +const QDLDL_int An = 10; |
| 10 | +const QDLDL_int Ap[] = { 0, 1, 2, 4, 5, 6, 8, 10, 12, 14, 17 }; |
| 11 | +const QDLDL_int Ai[] = { 0, 1, 1, 2, 3, 4, 1, 5, 0, 6, 3, 7, 6, 8, 1, 2, 9 }; |
| 12 | +const QDLDL_float Ax[] = { 1.0, 0.460641, -0.121189, 0.417928, 0.177828, 0.1, |
| 13 | + -0.0290058, -1.0, 0.350321, -0.441092, -0.0845395, -0.316228, |
| 14 | + 0.178663, -0.299077, 0.182452, -1.56506, -0.1 }; |
| 15 | +const QDLDL_float b[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; |
16 | 16 |
|
17 | 17 | int main()
|
18 |
| - |
19 | 18 | {
|
20 |
| - QDLDL_int i; // Counter |
21 |
| - |
22 |
| - //data for L and D factors |
23 |
| - QDLDL_int Ln = An; |
24 |
| - QDLDL_int *Lp; |
25 |
| - QDLDL_int *Li; |
26 |
| - QDLDL_float *Lx; |
27 |
| - QDLDL_float *D; |
28 |
| - QDLDL_float *Dinv; |
29 |
| - |
30 |
| - //data for elim tree calculation |
31 |
| - QDLDL_int *etree; |
32 |
| - QDLDL_int *Lnz; |
33 |
| - QDLDL_int sumLnz; |
34 |
| - |
35 |
| - //working data for factorisation |
36 |
| - QDLDL_int *iwork; |
37 |
| - QDLDL_bool *bwork; |
38 |
| - QDLDL_float *fwork; |
39 |
| - |
40 |
| - //Data for results of A\b |
41 |
| - QDLDL_float *x; |
42 |
| - |
43 |
| - |
44 |
| - /*-------------------------------- |
45 |
| - * pre-factorisation memory allocations |
46 |
| - *---------------------------------*/ |
47 |
| - |
48 |
| - //These can happen *before* the etree is calculated |
49 |
| - //since the sizes are not sparsity pattern specific |
50 |
| - |
51 |
| - //For the elimination tree |
52 |
| - etree = (QDLDL_int*)malloc(sizeof(QDLDL_int)*An); |
53 |
| - Lnz = (QDLDL_int*)malloc(sizeof(QDLDL_int)*An); |
54 |
| - |
55 |
| - //For the L factors. Li and Lx are sparsity dependent |
56 |
| - //so must be done after the etree is constructed |
57 |
| - Lp = (QDLDL_int*)malloc(sizeof(QDLDL_int)*(An+1)); |
58 |
| - D = (QDLDL_float*)malloc(sizeof(QDLDL_float)*An); |
59 |
| - Dinv = (QDLDL_float*)malloc(sizeof(QDLDL_float)*An); |
60 |
| - |
61 |
| - //Working memory. Note that both the etree and factor |
62 |
| - //calls requires a working vector of QDLDL_int, with |
63 |
| - //the factor function requiring 3*An elements and the |
64 |
| - //etree only An elements. Just allocate the larger |
65 |
| - //amount here and use it in both places |
66 |
| - iwork = (QDLDL_int*)malloc(sizeof(QDLDL_int)*(3*An)); |
67 |
| - bwork = (QDLDL_bool*)malloc(sizeof(QDLDL_bool)*An); |
68 |
| - fwork = (QDLDL_float*)malloc(sizeof(QDLDL_float)*An); |
69 |
| - |
70 |
| - /*-------------------------------- |
71 |
| - * elimination tree calculation |
72 |
| - *---------------------------------*/ |
73 |
| - sumLnz = QDLDL_etree(An,Ap,Ai,iwork,Lnz,etree); |
74 |
| - |
75 |
| - /*-------------------------------- |
76 |
| - * LDL factorisation |
77 |
| - *---------------------------------*/ |
78 |
| - |
79 |
| - //First allocate memory for Li and Lx |
80 |
| - Li = (QDLDL_int*)malloc(sizeof(QDLDL_int)*sumLnz); |
81 |
| - Lx = (QDLDL_float*)malloc(sizeof(QDLDL_float)*sumLnz); |
82 |
| - |
83 |
| - //now factor |
84 |
| - QDLDL_factor(An,Ap,Ai,Ax,Lp,Li,Lx,D,Dinv,Lnz,etree,bwork,iwork,fwork); |
85 |
| - |
86 |
| - /*-------------------------------- |
87 |
| - * solve |
88 |
| - *---------------------------------*/ |
89 |
| - x = (QDLDL_float*)malloc(sizeof(QDLDL_float)*An); |
90 |
| - |
91 |
| - //when solving A\b, start with x = b |
92 |
| - for(i=0;i < Ln; i++) x[i] = b[i]; |
93 |
| - QDLDL_solve(Ln,Lp,Li,Lx,Dinv,x); |
94 |
| - |
95 |
| - /*-------------------------------- |
96 |
| - * print factors and solution |
97 |
| - *---------------------------------*/ |
98 |
| - printf("\n"); |
99 |
| - printf("A (CSC format):\n"); |
100 |
| - print_line(); |
101 |
| - print_arrayi(Ap, An + 1, "A.p"); |
102 |
| - print_arrayi(Ai, Ap[An], "A.i"); |
103 |
| - print_arrayf(Ax, Ap[An], "A.x"); |
104 |
| - printf("\n\n"); |
105 |
| - |
106 |
| - printf("elimination tree:\n"); |
107 |
| - print_line(); |
108 |
| - print_arrayi(etree, Ln, "etree"); |
109 |
| - print_arrayi(Lnz, Ln, "Lnz"); |
110 |
| - printf("\n\n"); |
111 |
| - |
112 |
| - printf("L (CSC format):\n"); |
113 |
| - print_line(); |
114 |
| - print_arrayi(Lp, Ln + 1, "L.p"); |
115 |
| - print_arrayi(Li, Lp[Ln], "L.i"); |
116 |
| - print_arrayf(Lx, Lp[Ln], "L.x"); |
117 |
| - printf("\n\n"); |
118 |
| - |
119 |
| - printf("D:\n"); |
120 |
| - print_line(); |
121 |
| - print_arrayf(D, An, "diag(D) "); |
122 |
| - print_arrayf(Dinv, An, "diag(D^{-1})"); |
123 |
| - printf("\n\n"); |
124 |
| - |
125 |
| - printf("solve results:\n"); |
126 |
| - print_line(); |
127 |
| - print_arrayf(b, An, "b"); |
128 |
| - print_arrayf(x, An, "A\\b"); |
129 |
| - printf("\n\n"); |
130 |
| - |
131 |
| - |
132 |
| - /*-------------------------------- |
133 |
| - * clean up |
134 |
| - *---------------------------------*/ |
135 |
| - free(Lp); |
136 |
| - free(Li); |
137 |
| - free(Lx); |
138 |
| - free(D); |
139 |
| - free(Dinv); |
140 |
| - free(etree); |
141 |
| - free(Lnz); |
142 |
| - free(iwork); |
143 |
| - free(bwork); |
144 |
| - free(fwork); |
145 |
| - free(x); |
146 |
| - |
147 |
| -return 0 ; |
148 |
| - |
149 |
| - |
| 19 | + QDLDL_int i; // Counter |
| 20 | + |
| 21 | + // data for L and D factors |
| 22 | + QDLDL_int Ln = An; |
| 23 | + QDLDL_int* Lp; |
| 24 | + QDLDL_int* Li; |
| 25 | + QDLDL_float* Lx; |
| 26 | + QDLDL_float* D; |
| 27 | + QDLDL_float* Dinv; |
| 28 | + |
| 29 | + // data for elim tree calculation |
| 30 | + QDLDL_int* etree; |
| 31 | + QDLDL_int* Lnz; |
| 32 | + QDLDL_int sumLnz; |
| 33 | + |
| 34 | + // working data for factorisation |
| 35 | + QDLDL_int* iwork; |
| 36 | + QDLDL_bool* bwork; |
| 37 | + QDLDL_float* fwork; |
| 38 | + |
| 39 | + // Data for results of A\b |
| 40 | + QDLDL_float* x; |
| 41 | + |
| 42 | + /*-------------------------------- |
| 43 | + * pre-factorisation memory allocations |
| 44 | + *---------------------------------*/ |
| 45 | + |
| 46 | + // These can happen *before* the etree is calculated |
| 47 | + // since the sizes are not sparsity pattern specific |
| 48 | + |
| 49 | + // For the elimination tree |
| 50 | + etree = (QDLDL_int*) malloc(sizeof(QDLDL_int) * An); |
| 51 | + Lnz = (QDLDL_int*) malloc(sizeof(QDLDL_int) * An); |
| 52 | + |
| 53 | + // For the L factors. Li and Lx are sparsity dependent |
| 54 | + // so must be done after the etree is constructed |
| 55 | + Lp = (QDLDL_int*) malloc(sizeof(QDLDL_int) * (An + 1)); |
| 56 | + D = (QDLDL_float*) malloc(sizeof(QDLDL_float) * An); |
| 57 | + Dinv = (QDLDL_float*) malloc(sizeof(QDLDL_float) * An); |
| 58 | + |
| 59 | + // Working memory. Note that both the etree and factor |
| 60 | + // calls requires a working vector of QDLDL_int, with |
| 61 | + // the factor function requiring 3*An elements and the |
| 62 | + // etree only An elements. Just allocate the larger |
| 63 | + // amount here and use it in both places |
| 64 | + iwork = (QDLDL_int*) malloc(sizeof(QDLDL_int) * (3 * An)); |
| 65 | + bwork = (QDLDL_bool*) malloc(sizeof(QDLDL_bool) * An); |
| 66 | + fwork = (QDLDL_float*) malloc(sizeof(QDLDL_float) * An); |
| 67 | + |
| 68 | + /*-------------------------------- |
| 69 | + * elimination tree calculation |
| 70 | + *---------------------------------*/ |
| 71 | + sumLnz = QDLDL_etree(An, Ap, Ai, iwork, Lnz, etree); |
| 72 | + |
| 73 | + /*-------------------------------- |
| 74 | + * LDL factorisation |
| 75 | + *---------------------------------*/ |
| 76 | + |
| 77 | + // First allocate memory for Li and Lx |
| 78 | + Li = (QDLDL_int*) malloc(sizeof(QDLDL_int) * sumLnz); |
| 79 | + Lx = (QDLDL_float*) malloc(sizeof(QDLDL_float) * sumLnz); |
| 80 | + |
| 81 | + // now factor |
| 82 | + QDLDL_factor(An, Ap, Ai, Ax, Lp, Li, Lx, D, Dinv, Lnz, etree, bwork, iwork, fwork); |
| 83 | + |
| 84 | + /*-------------------------------- |
| 85 | + * solve |
| 86 | + *---------------------------------*/ |
| 87 | + x = (QDLDL_float*) malloc(sizeof(QDLDL_float) * An); |
| 88 | + |
| 89 | + // when solving A\b, start with x = b |
| 90 | + for(i = 0; i < Ln; i++) { |
| 91 | + x[i] = b[i]; |
| 92 | + } |
| 93 | + QDLDL_solve(Ln, Lp, Li, Lx, Dinv, x); |
| 94 | + |
| 95 | + /*-------------------------------- |
| 96 | + * print factors and solution |
| 97 | + *---------------------------------*/ |
| 98 | + printf("\n"); |
| 99 | + printf("A (CSC format):\n"); |
| 100 | + print_line(); |
| 101 | + print_arrayi(Ap, An + 1, "A.p"); |
| 102 | + print_arrayi(Ai, Ap[An], "A.i"); |
| 103 | + print_arrayf(Ax, Ap[An], "A.x"); |
| 104 | + printf("\n\n"); |
| 105 | + |
| 106 | + printf("elimination tree:\n"); |
| 107 | + print_line(); |
| 108 | + print_arrayi(etree, Ln, "etree"); |
| 109 | + print_arrayi(Lnz, Ln, "Lnz"); |
| 110 | + printf("\n\n"); |
| 111 | + |
| 112 | + printf("L (CSC format):\n"); |
| 113 | + print_line(); |
| 114 | + print_arrayi(Lp, Ln + 1, "L.p"); |
| 115 | + print_arrayi(Li, Lp[Ln], "L.i"); |
| 116 | + print_arrayf(Lx, Lp[Ln], "L.x"); |
| 117 | + printf("\n\n"); |
| 118 | + |
| 119 | + printf("D:\n"); |
| 120 | + print_line(); |
| 121 | + print_arrayf(D, An, "diag(D) "); |
| 122 | + print_arrayf(Dinv, An, "diag(D^{-1})"); |
| 123 | + printf("\n\n"); |
| 124 | + |
| 125 | + printf("solve results:\n"); |
| 126 | + print_line(); |
| 127 | + print_arrayf(b, An, "b"); |
| 128 | + print_arrayf(x, An, "A\\b"); |
| 129 | + printf("\n\n"); |
| 130 | + |
| 131 | + /*-------------------------------- |
| 132 | + * clean up |
| 133 | + *---------------------------------*/ |
| 134 | + free(Lp); |
| 135 | + free(Li); |
| 136 | + free(Lx); |
| 137 | + free(D); |
| 138 | + free(Dinv); |
| 139 | + free(etree); |
| 140 | + free(Lnz); |
| 141 | + free(iwork); |
| 142 | + free(bwork); |
| 143 | + free(fwork); |
| 144 | + free(x); |
| 145 | + |
| 146 | + return 0; |
150 | 147 | }
|
151 | 148 |
|
152 |
| - |
153 |
| -void print_line(void){ |
154 |
| - printf("--------------------------\n"); |
| 149 | +void print_line(void) { |
| 150 | + printf("--------------------------\n"); |
155 | 151 | }
|
156 | 152 |
|
157 |
| -void print_arrayi(const QDLDL_int* data, QDLDL_int n,char* varName){ |
158 |
| - |
159 |
| - QDLDL_int i; |
160 |
| - printf("%s = [",varName); |
161 |
| - for(i=0; i< n; i++){ |
162 |
| - printf("%i,",(int)data[i]); |
163 |
| - } |
164 |
| - printf("]\n"); |
| 153 | +void print_arrayi(const QDLDL_int* data, QDLDL_int n, char* varName) { |
| 154 | + QDLDL_int i; |
| 155 | + printf("%s = [", varName); |
165 | 156 |
|
| 157 | + for(i = 0; i < n; i++) { |
| 158 | + printf("%i,", (int) data[i]); |
| 159 | + } |
| 160 | + printf("]\n"); |
166 | 161 | }
|
167 | 162 |
|
168 |
| -void print_arrayf(const QDLDL_float* data, QDLDL_int n, char* varName){ |
169 |
| - |
170 |
| - QDLDL_int i; |
171 |
| - printf("%s = [",varName); |
172 |
| - for(i=0; i< n; i++){ |
173 |
| - printf("%.3g,",data[i]); |
174 |
| - } |
175 |
| - printf("]\n"); |
| 163 | +void print_arrayf(const QDLDL_float* data, QDLDL_int n, char* varName) { |
| 164 | + QDLDL_int i; |
| 165 | + printf("%s = [", varName); |
176 | 166 |
|
| 167 | + for(i = 0; i < n; i++) { |
| 168 | + printf("%.3g,", data[i]); |
| 169 | + } |
| 170 | + printf("]\n"); |
177 | 171 | }
|
0 commit comments