source: trunk/ADOL-C/examples/additional_examples/lufact/LUsolve.cpp @ 88

Last change on this file since 88 was 88, checked in by utke, 10 years ago

compiler message: ISO C++ does not support the `%le' printf format

  • Property svn:keywords set to Author Date Id Revision
File size: 5.4 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     LUsolve.cpp
4 Revision: $Id: LUsolve.cpp 88 2010-02-17 18:33:52Z utke $
5 Contents: example for
6             * Computation of the determinant of a matrix
7               by LU-decomposition of the system matrix without pivoting
8             * application of tapedoc to observe taping of
9               the new op_codes for the elementary operations
10                 
11                     y += x1 * x2;
12                     y -= x1 * x2;           
13
14 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
15               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
16 
17 This file is part of ADOL-C. This software is provided as open source.
18 Any use, reproduction, or distribution of the software constitutes
19 recipient's acceptance of the terms of the accompanying license file.
20 
21---------------------------------------------------------------------------*/
22
23/****************************************************************************/
24/*                                                                 INCLUDES */
25#include <../examples/additional_examples/lufact/LU.h>
26
27/****************************************************************************/
28/*                                                             MAIN PROGRAM */
29int main() { /*------------------------------------------------------------------------*/
30    /* variables */
31    const int tag   = 1;                       // tape tag
32    const int size  = 5;                       // system size
33    const int indep = size*size+size;          // # of indeps
34    const int depen = size;                    // # of deps
35
36    double  A[size][size], a1[size], a2[size], // passive variables
37    b[size], x[size];
38    adouble **AA, *AAp, *Abx;                  // active variables
39    double *args = myalloc1(indep);            // arguments
40    double **jac = myalloc2(depen,indep);      // the Jacobian
41    double *laghessvec = myalloc1(indep);      // Hessian-vector product
42
43    int i,j;
44
45
46    /*------------------------------------------------------------------------*/
47    /* Info */
48    fprintf(stdout,"LINEAR SYSTEM SOLVING by "
49            "LU-DECOMPOSITION (ADOL-C Example)\n\n");
50
51
52    /*------------------------------------------------------------------------*/
53    /* Allocation und initialization of the system matrix */
54    AA  = new adouble*[size];
55    AAp = new adouble[size*size];
56    for (i=0; i<size; i++) {
57        AA[i] = AAp;
58        AAp += size;
59    }
60    Abx = new adouble[size];
61    for(i=0; i<size; i++) {
62        a1[i] = i*0.25;
63        a2[i] = i*0.33;
64    }
65    for(i=0; i<size; i++) {
66        for(j=0; j<size; j++)
67            A[i][j] = a1[i]*a2[j];
68        A[i][i] += i+1;
69        b[i] = -i-1;
70    }
71
72
73    /*------------------------------------------------------------------------*/
74    /* Taping the computation of the determinant */
75    trace_on(tag);
76    /* marking indeps */
77    for(i=0; i<size; i++)
78        for(j=0; j<size; j++)
79            AA[i][j] <<= (args[i*size+j] = A[i][j]);
80    for(i=0; i<size; i++)
81        Abx[i] <<= (args[size*size+i] = b[i]);
82    /* LU-factorization and computation of solution */
83    LUfact(size,AA);
84    LUsolve(size,AA,Abx);
85    /* marking deps */
86    for (i=0; i<size; i++)
87        Abx[i] >>= x[i];
88    trace_off();
89    fprintf(stdout," x[0] (original):  %16.4E\n",x[0]);
90
91
92    /*------------------------------------------------------------------------*/
93    /* Recomputation  */
94    function(tag,depen,indep,args,x);
95    fprintf(stdout," x[0] (from tape): %16.4E\n",x[0]);
96
97
98    /*------------------------------------------------------------------------*/
99    /* Computation of Jacobian */
100    jacobian(tag,depen,indep,args,jac);
101    fprintf(stdout," Jacobian:\n");
102    for (i=0; i<depen; i++) {
103        for (j=0; j<indep; j++)
104            fprintf(stdout," %14.6E",jac[i][j]);
105        fprintf(stdout,"\n");
106    }
107
108    /*------------------------------------------------------------------------*/
109    /* Computation of Lagrange-Hessian-vector product */
110    lagra_hess_vec(tag,depen,indep,args,args,x,laghessvec);
111    fprintf(stdout," Part of Lagrange-Hessian-vector product:\n");
112    for (i=0; i<size; i++) {
113        for (j=0; j<size; j++)
114            fprintf(stdout," %14.6E",laghessvec[i*size+j]);
115        fprintf(stdout,"\n");
116    }
117
118
119    /*------------------------------------------------------------------------*/
120    /* Tape-documentation */
121    tape_doc(tag,depen,indep,args,x);
122
123
124    /*------------------------------------------------------------------------*/
125    /* Tape statistics */
126    int tape_stats[STAT_SIZE];
127    tapestats(tag,tape_stats);
128
129    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
130    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
131    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
132    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
133    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
134    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
135    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
136    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
137
138    /*------------------------------------------------------------------------*/
139    /* That's it */
140    return 1;
141}
142
143
144
145
146
147
148
149
Note: See TracBrowser for help on using the repository browser.