1 /*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File: luexam.cpp
4 Revision: $Id: luexam.cpp 431 2013-06-18 20:14:07Z kulshres $
5 Contents: computation of LU factorization with pivoting
6
7 Copyright (c) Kshitij Kulshreshtha
8
9 This file is part of ADOL-C. This software is provided as open source.
10 Any use, reproduction, or distribution of the software constitutes
11 recipient's acceptance of the terms of the accompanying license file.
12
13 ---------------------------------------------------------------------------*/
14
15 #include <adolc/adouble.h>
16 #include <adolc/advector.h>
17 #include <adolc/taping.h>
18
19 #include <adolc/interfaces.h>
20 #include <adolc/drivers/drivers.h>
21
22 #include <iostream>
23 #include <fstream>
24 #include <cstring>
25 #include <iomanip>
26 #include <sstream>
27
28 using namespace std;
29
findmaxindex(const size_t n,const advector & col,const size_t k)30 adouble findmaxindex(const size_t n, const advector& col, const size_t k) {
31 adouble idx = k;
32 for (size_t j = k + 1; j < n; j++ )
33 condassign(idx,(fabs(col[j]) - fabs(col[idx])), adouble((double)j));
34 return idx;
35 }
36
37 // Assuming A is stored row-wise in the vector
38
lufactorize(const size_t n,advector & A,advector & p)39 void lufactorize(const size_t n, advector& A, advector& p) {
40 adouble idx, tmp;
41 for (size_t j = 0; j < n; j++)
42 p[j] = j;
43 for (size_t k = 0; k < n; k++) {
44 advector column(n);
45 for(size_t j = 0; j < n; j++)
46 condassign(column[j], adouble(double(j - k + 1)), A[j*n + k]);
47 idx = findmaxindex(n, column, k);
48 tmp = p[k];
49 p[k] = p[idx];
50 p[idx] = tmp;
51 for(size_t j = 0; j < n; j++) {
52 tmp = A[k*n + j];
53 A[k*n + j] = A[idx*n + j];
54 A[idx*n + j] = tmp;
55 }
56 tmp = 1.0/A[k*n + k];
57 for (size_t i = k + 1; i < n; i++) {
58 A[i*n + k] *= tmp;
59 for (size_t j = k + 1; j < n; j++) {
60 A[i*n + j] -= A[i*n + k] * A[k*n+j];
61 }
62 }
63 }
64 }
65
Lsolve(const size_t n,const advector & A,const advector & p,advector & b,advector & x)66 void Lsolve(const size_t n, const advector& A, const advector& p, advector& b, advector& x) {
67 for (size_t j = 0; j < n; j++) {
68 x[j] = b[p[j]];
69 for (size_t k = j+1; k <n; k++) {
70 b[p[k]] -= A[k*n+j]*x[j];
71 }
72 }
73 }
74
Rsolve(const size_t n,const advector & A,advector & x)75 void Rsolve(const size_t n, const advector& A, advector& x) {
76 for (size_t j = 1; j <= n; j++) {
77 x[n-j] *= 1.0/A[(n-j)*n + n-j];
78 for (size_t k = 0; k < n-j; k++) {
79 x[k] -= A[k*n + n-j]*x[n-j];
80 }
81 }
82 }
83
printL(const size_t n,const advector & A,ostream & outf=std::cout)84 void printL(const size_t n, const advector& A, ostream &outf = std::cout) {
85 for (size_t i = 0; i < n; i++) {
86 for (size_t j = 0; j < n; j++)
87 if (j < i)
88 outf << setw(8) << A[i*n + j].value() << " ";
89 else if (j == i)
90 outf << setw(8) << 1.0 << " ";
91 else
92 outf << setw(8) << 0.0 << " ";
93 outf << "\n";
94 }
95 }
96
printR(const size_t n,const advector & A,ostream & outf=std::cout)97 void printR(const size_t n, const advector& A, ostream &outf = std::cout) {
98 for (size_t i = 0; i < n; i++) {
99 for (size_t j = 0; j < n; j++)
100 if (j >= i)
101 outf << setw(8) << A[i*n + j].value() << " ";
102 else
103 outf << setw(8) << 0.0 << " ";
104 outf << "\n";
105 }
106 }
107
norm2(const double * const v,const size_t n)108 double norm2(const double *const v, const size_t n)
109 {
110 size_t j;
111 double abs,scale,sum;
112 scale=0.0;
113 sum=0.0;
114 for (j=0; j<n; j++) {
115 if (v[j] != 0.0) {
116 abs = fabs(v[j]);
117 if (scale <= abs) {
118 sum = 1.0 + sum * (scale/abs)*(scale/abs);
119 scale = abs;
120 } else
121 sum += (abs/scale)*(abs/scale);
122 }
123 }
124 sum = sqrt(sum)*scale;
125 return sum;
126 }
127
scalar(double * x,double * y,size_t n)128 double scalar(double *x, double *y, size_t n)
129 {
130 size_t j;
131 int8_t sign;
132 double abs,scale,sum,prod;
133 scale = 0.0;
134 sum = 0.0;
135 for(j=0; j<n; j++) {
136 sign = 1;
137 prod = x[j]*y[j];
138 if( prod != 0.0) {
139 abs = prod;
140 if (abs < 0.0) {
141 sign = -1;
142 abs = -abs;
143 }
144 if( scale <= fabs(abs)) {
145 sum = sign * 1.0 + sum *(scale/abs);
146 scale = abs;
147 } else
148 sum+=sign*(abs/scale);
149 }
150 }
151 sum = sum*scale;
152 return sum;
153 }
154
matvec(const double * const A,const size_t m,const double * const b,const size_t n,double * const ret)155 void matvec(const double *const A, const size_t m, const double *const b, const size_t n, double *const ret)
156 {
157 size_t i,j;
158 memset(ret,0,n*sizeof(double));
159 for (i=0; i<m; i++) {
160 double abs, scale = 0.0, prod, sum = 0.0;
161 int8_t sign;
162 for (j=0; j<n; j++) {
163 sign = 1;
164 prod = A[i*n+j]*b[j];
165 if (prod != 0.0) {
166 abs = prod;
167 if (abs < 0.0) {
168 sign = -1;
169 abs = -abs;
170 }
171 if (scale <=fabs(abs)) {
172 sum = sign * 1.0 + sum * (scale/abs);
173 scale = abs;
174 } else
175 sum += sign*(abs/scale);
176 }
177 }
178 ret[i] = sum*scale;
179 }
180 }
181
residue(const double * const A,const size_t m,const double * const b,const size_t n,const double * const x,double * const ret)182 void residue(const double *const A, const size_t m, const double *const b, const size_t n, const double *const x, double *const ret) {
183 double *b2 = new double[n];
184 matvec(A,m,x,n,b2);
185 for (size_t i = 0; i < n; i++)
186 ret[i] = b[i] - b2[i];
187 delete[] b2;
188 }
189
normresidue(const double * const A,const size_t m,const double * const b,const size_t n,const double * const x)190 double normresidue(const double *const A, const size_t m, const double *const b, const size_t n, const double*const x) {
191 double *res = new double[n];
192 residue(A,m,b,n,x,res);
193 double ans = norm2(res, n);
194 delete[] res;
195 return ans;
196 }
197
main()198 int main() {
199 int tag = 1;
200 int keep = 1;
201 int n;
202 string matrixname, rhsname;
203 ifstream matf, rhsf;
204 double *mat, *rhs, *ans, err, normx, normb;
205
206 cout << "COMPUTATION OF LU-Factorization with pivoting (ADOL-C Documented Example)\n\n";
207 cout << "order of matrix = ? \n"; // select matrix size
208 cin >> n;
209
210 cout << "---------------------------------\nNow tracing:\n";
211 rhs = new double[n*n + n];
212 mat = rhs + n;
213 ans = new double[n];
214 cout << "file name for matrix = ?\n";
215 cin >> matrixname;
216 cout << "file name for rhs = ?\n";
217 cin >> rhsname;
218
219
220 matf.open(matrixname.c_str());
221 for (size_t i = 0; i < n*n; i++)
222 matf >> mat[i];
223 matf.close();
224
225 rhsf.open(rhsname.c_str());
226 for (size_t i = 0; i < n; i++)
227 rhsf >> rhs[i];
228 rhsf.close();
229
230 {
231 trace_on(tag,keep); // tag=1=keep
232 advector A(n*n), b(n), x(n), p(n);
233 for(size_t i = 0; i < n; i++)
234 b[i] <<= rhs[i];
235 for(size_t i = 0; i < n*n; i++)
236 A[i] <<= mat[i];
237 lufactorize(n, A, p);
238 Lsolve(n, A, p, b, x);
239 Rsolve(n, A, x);
240 for(size_t i = 0; i < n; i++)
241 x[i] >>= ans[i];
242 trace_off();
243 }
244
245 err = normresidue(mat, n, rhs, n, ans);
246 normb = norm2(rhs, n);
247 normx = norm2(ans, n);
248 cout << "Norm rhs = " << normb <<"\n";
249 cout << "Norm solution = " << normx <<"\n";
250 cout << "Norm of residue = " << err <<"\t relative error = " << err/normx << "\n";
251
252 cout << "---------------------------------\nNow computing from trace:\n";
253 cout << "file name for matrix = ?\n";
254 cin >> matrixname;
255 cout << "file name for rhs = ?\n";
256 cin >> rhsname;
257
258
259 matf.open(matrixname.c_str());
260 for (size_t i = 0; i < n*n; i++)
261 matf >> mat[i];
262 matf.close();
263
264 rhsf.open(rhsname.c_str());
265 for (size_t i = 0; i < n; i++)
266 rhsf >> rhs[i];
267 rhsf.close();
268
269 zos_forward(tag, n, n*n + n, keep, rhs, ans);
270
271 err = normresidue(mat, n, rhs, n, ans);
272 normb = norm2(rhs, n);
273 normx = norm2(ans, n);
274 cout << "Norm rhs = " << normb <<"\n";
275 cout << "Norm solution = " << normx <<"\n";
276 cout << "Norm of residue = " << err <<"\t relative error = " << err/normx <<"\n";
277 double *ansbar = new double[n];
278 double *matcol = new double[n];
279 double *rhsbar = new double[n*n+n];
280 double *matbar = rhsbar + n;
281 double scprod = 0.0;
282
283 memset(rhsbar, 0, (n*n+n)*sizeof(double));
284 memset(ansbar, 0, n*sizeof(double));
285 for (size_t k = 0; k < n; k++) {
286 cout << "computing gradient of element " << k + 1 << " of solution w.r.t. matrix elements and rhs\n";
287 ansbar[k] = 1.0;
288
289 fos_reverse(tag, n, n*n+n, ansbar, rhsbar);
290
291 for (size_t i = 0; i < n*n + n; i++)
292 cout << "bar[" << i << "] = " << rhsbar[i] << "\n";
293
294 for (size_t j = 0; j < n; j++) {
295 for (size_t i = 0; i < n; i++)
296 matcol[i] = mat[i*n + j];
297 scprod = scalar(rhsbar, matcol, n);
298 cout << "gradient w.r.t. rhs times column " << j + 1 << " of matrix = " << scprod << "\n";
299 }
300 ansbar[k] = 0.0;
301 }
302 delete[] ansbar;
303 delete[] matcol;
304 delete[] rhsbar;
305
306 delete[] rhs;
307 delete[] ans;
308 }
309