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/*----------------------------------------------------------------------------
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
21#include <taping.h>           // use of taping
22
23#include <iostream>
24using namespace std;
25
26/****************************************************************************/
28int n;
29adouble **A;                        // A is an n x n matrix
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    {
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
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;
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.
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++) {