Last change on this file was 338, checked in by kulshres, 7 years ago

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

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