source: trunk/ADOL-C/examples/additional_examples/lufact/LUdet.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.2 KB
RevLine 
[40]1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     LUdet.cpp
[42]4 Revision: $Id: LUdet.cpp 88 2010-02-17 18:33:52Z utke $
[40]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/****************************************************************************/
29/*                                                             MAIN PROGRAM */
30int main() { /*------------------------------------------------------------------------*/
31    /* variables */
32    const int tag   = 1;                       // tape tag
33    const int size  = 5;                       // system size
34    const int indep = size*size;               // # of indeps
35    const int depen = 1;                       // # of deps
36
37    double  A[size][size], a1[size], a2[size], det; // passive variables
38    adouble **AA, *AAp,     Adet;                   // active variables
39    double *args  = myalloc1(indep);                // arguments
40    double *grad  = myalloc1(indep);                // the gradient
41    double **hess = myalloc2(indep,indep);          // the hessian
42
43    int i,j;
44
45
46    /*------------------------------------------------------------------------*/
47    /* Info */
48    fprintf(stdout,"DETERMINANT by LU-DECOMPOSITION (ADOL-C Example)\n\n");
49
50
51    /*------------------------------------------------------------------------*/
52    /* Allcoation und initialization of the system matrix */
53    AA  = new adouble*[size];
54    AAp = new adouble[size*size];
55    for (i=0; i<size; i++) {
56        AA[i] = AAp;
57        AAp += size;
58    }
59    for(i=0; i<size; i++) {
60        a1[i] = i*0.25;
61        a2[i] = i*0.33;
62    }
63    for(i=0; i<size; i++) {
64        for(j=0; j<size; j++)
65            A[i][j] = a1[i]*a2[j];
66        A[i][i] += i+1;
67    }
68
69
70    /*------------------------------------------------------------------------*/
71    /* Taping the computation of the determinant */
72    trace_on(tag);
73    /* marking indeps */
74    for(i=0; i<size; i++)
75        for(j=0; j<size; j++)
76            AA[i][j] <<= (args[i*size+j] = A[i][j]);     // indep. vars
77    /* LU-factorization and computation of determinant */
78    LUfact(size,AA);
79    Adet = AA[0][0];
80    for (i=1; i<size; i++)
81        Adet *= AA[i][i];
82    /* marking deps */
83    Adet >>= det;
84    trace_off();
[88]85    fprintf(stdout," Determinant (original):  %16.4E\n",det);
[40]86
87
88    /*------------------------------------------------------------------------*/
89    /* Recomputation of determinant */
90    function(tag,depen,indep,args,&det);
[88]91    fprintf(stdout," Determinant (from tape): %16.4E\n",det);
[40]92
93
94    /*------------------------------------------------------------------------*/
95    /* Computation of gradient */
96    gradient(tag,indep,args,grad);
97    fprintf(stdout," Gradient:\n");
98    for (i=0; i<size; i++) {
99        for (j=0; j<size; j++)
[88]100            fprintf(stdout," %14.6E",grad[i*size+j]);
[40]101        fprintf(stdout,"\n");
102    }
103
104    /*------------------------------------------------------------------------*/
105    /* Computation of hessian */
106    hessian(tag,indep,args,hess);
107    fprintf(stdout," Part of Hessian:\n");
108    for (i=0; i<size; i++) {
109        for (j=0; j<size; j++)
[88]110            fprintf(stdout," %14.6E",hess[0][i*size+j]);
[40]111        fprintf(stdout,"\n");
112    }
113
114    /*------------------------------------------------------------------------*/
115    /* Tape-documentation */
116    tape_doc(tag,depen,indep,args,&det);
117
118    /*------------------------------------------------------------------------*/
119    /* Tape statistics */
120    int tape_stats[STAT_SIZE];
121    tapestats(tag,tape_stats);
122
123    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
124    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
125    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
126    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
127    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
128    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
129    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
130    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
131
132    /*------------------------------------------------------------------------*/
133    /* That's it */
134    return 1;
135}
136
137
138
139
140
141
142
143
Note: See TracBrowser for help on using the repository browser.