source: trunk/ADOL-C/examples/detexam.cpp @ 42

Last change on this file since 42 was 42, checked in by awalther, 10 years ago

set svn keywords property

  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     detexam.cpp
4 Revision: $Id: detexam.cpp 42 2009-07-15 18:37:17Z awalther $
5 Contents: computation of determinants, described in the manual
6
7 Copyright (c) Andrea Walther, Andreas Griewank, Andreas Kowarz,
8               Hristo Mitev, Sebastian Schlenkrich, Jean Utke, Olaf Vogel
9 
10 This file is part of ADOL-C. This software is provided as open source.
11 Any use, reproduction, or distribution of the software constitutes
12 recipient's acceptance of the terms of the accompanying license file.
13 
14---------------------------------------------------------------------------*/
15
16/****************************************************************************/
17/*                                                                 INCLUDES */
18#include <adouble.h>          // use of active doubles
19#include <interfaces.h>       // use of basic forward/reverse
20// interfaces of ADOL-C
21#include <taping.h>           // use of taping
22
23#include <iostream>
24using namespace std;
25
26/****************************************************************************/
27/*                                                          ADOUBLE ROUTINE */
28int n;
29adouble **A;                        // A is an n x n matrix
30adouble zero = 0;
31
32adouble det(int k, int m)           // k <= n is the order of the submatrix
33{ if (m == 0)                       // its column indices
34        return 1.0;
35    else                              // are encoded in m
36    {
37        adouble *pt = A[k-1];
38        adouble   t = zero;
39        int p = 1;
40        int s;
41        if (k%2)
42            s = 1;
43        else
44            s = -1;
45        for(int i=0; i<n; i++) {
46            int p1 = 2*p;
47            if (m%p1 >= p) {
48                if (m == p) {
49                    if (s>0)
50                        t += *pt;
51                    else
52                        t -= *pt;
53                } else {
54                    if (s>0)
55                        t += *pt*det(k-1, m-p); // recursive call to det
56                    else
57                        t -= *pt*det(k-1, m-p); // recursive call to det
58                }
59                s = -s;
60            }
61            ++pt;
62            p = p1;
63        }
64        return t;
65    }
66}
67
68/****************************************************************************/
69/*                                                             MAIN PROGRAM */
70int main() {
71    int i,j, m = 1;
72    int tag = 1;
73    int keep = 1;
74
75    cout << "COMPUTATION OF DETERMINANTS (ADOL-C Documented Example)\n\n";
76    cout << "order of matrix = ? \n"; // select matrix size
77    cin >> n;
78
79    A = new adouble*[n];
80    adouble ad;
81
82    trace_on(tag,keep);               // tag=1=keep
83    double detout = 0.0, diag = 1.0;// here keep the intermediates for
84    for (i=0; i<n; i++)             // the subsequent call to reverse
85    { m *= 2;
86        A[i] = new adouble[n];
87        for (j=0; j<n; j++)
88            A[i][j] <<= j/(1.0+i);      // make all elements of A independent
89        diag += A[i][i].value();       // value(adouble) converts to double
90        A[i][i] += 1.0;
91    }
92    ad = det(n,m-1);                // actual function call.
93    ad >>= detout;
94    printf("\n %f - %f = %f  (should be 0)\n",detout,diag,detout-diag);
95    trace_off();
96
97    double u[1];
98    u[0] = 1.0;
99    double* B = new double[n*n];
100
101    reverse(tag,1,n*n,0,u,B);         // call reverse to calculate the gradient
102
103    cout << " \n first base? : ";
104    for (i=0; i<n; i++) {
105        adouble sum = 0;
106        for (j=0; j<n; j++)             // the matrix A times the first n
107            sum += A[i][j]*B[j];          // components of the gradient B
108        cout << sum.value() << " ";      // must be a Cartesian basis vector
109    }
110    cout << "\n";
111
112    return 1;
113}
114
Note: See TracBrowser for help on using the repository browser.