source: trunk/ADOL-C/examples/additional_examples/detexam/detexam.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: 7.4 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     detexam.cpp
4 Revision: $Id: detexam.cpp 88 2010-02-17 18:33:52Z utke $
5 Contents: modified computation of determinants
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 <adolc.h>
19#include <../examples/additional_examples/clock/myclock.h>
20
21
22/****************************************************************************/
23/*                                                           DOUBLE ROUTINE */
24int n,it;
25double** PA;
26double pdet( int k, int m ) {
27    if (m == 0)
28        return 1.0 ;
29    else {
30        double* pt = PA[k-1];
31        double t = 0;
32        int p = 1;
33        int s;
34        if (k%2)
35            s = 1;
36        else
37            s = -1;
38        for (int i=0; i<n; i++) {
39            int p1 = 2*p;
40            if (m%p1 >= p) {
41                if (m == p) {
42                    if (s>0)
43                        t += *pt;
44                    else
45                        t -= *pt;
46                } else {
47                    if (s>0)
48                        t += *pt*pdet(k-1, m-p);
49                    else
50                        t -= *pt*pdet(k-1, m-p);
51                }
52                s = -s;
53            }
54            ++pt;
55            p = p1;
56        }
57        return t;
58    }
59}
60
61/****************************************************************************/
62/*                                                          ADOUBLE ROUTINE */
63adouble** A;
64adouble zero = 0;
65adouble det( int k, int m ) {
66    if (m == 0)
67        return 1.0;
68    else {
69        adouble* pt = A[k-1];
70        adouble t = zero;
71        int p = 1;
72        int s;
73        if (k%2)
74            s = 1;
75        else
76            s = -1;
77        for (int i=0; i<n; i++) {
78            int p1 = 2*p;
79            if (m%p1 >= p) {
80                if (m == p) {
81                    if (s>0)
82                        t += *pt;
83                    else
84                        t -= *pt;
85                } else {
86                    if (s>0)
87                        t += *pt*det(k-1, m-p);
88                    else
89                        t -= *pt*det(k-1, m-p);
90                }
91                s = -s;
92            }
93            ++pt;
94            p = p1;
95        }
96        return t;
97    }
98}
99
100/****************************************************************************/
101/*                                                             MAIN PROGRAM */
102int main() {
103    int i, j;
104    int tag = 1;
105    fprintf(stdout,"COMPUTATION OF DETERMINANTS Type 1 (ADOL-C Example)\n\n");
106    fprintf(stdout,"order of matrix = ? \n");
107    scanf("%d",&n);
108    A  = new adouble*[n];
109    PA = new double*[n];
110    int n2 = n*n;
111    double* a = new double[n2];
112
113    /*--------------------------------------------------------------------------*/
114    /* Preparation */
115    double diag = 0;
116    int m = 1;
117    double* pa = a;
118    for (i=0; i<n; i++) {
119        m *= 2;
120        PA[i] = new double[n];
121        double* ppt = PA[i];
122        for (j=0; j<n; j++) {
123            *ppt++ = j/(1.0+i);
124            *pa++  = j/(1.0+i);
125        }
126        diag += PA[i][i];   // val corrected to value 2/23/91
127        PA[i][i] += 1.0;
128        a[i*n+i] += 1.0;
129    }
130    diag += 1;
131
132    /*--------------------------------------------------------------------------*/
133    double t00 = myclock();                               /* 0. time (taping) */
134    trace_on(tag);
135    for (i=0; i<n; i++) {
136        A[i] = new adouble[n];
137        adouble* pt = A[i];
138        double* ppt = PA[i];
139        for (j=0; j<n; j++)
140            *pt++ <<= *ppt++;
141    }
142    adouble deter;
143    deter = det(n,m-1);
144    double detout = 0.0;
145    deter >>= detout;
146    trace_off();
147    double t01 = myclock();
148    fprintf(stdout,"\n %f =? %f should be the same \n",detout,diag);
149
150    /*--------------------------------------------------------------------------*/
151    int tape_stats[STAT_SIZE];
152    tapestats(tag,tape_stats);
153
154    fprintf(stdout,"\n    independents            %d\n",tape_stats[NUM_INDEPENDENTS]);
155    fprintf(stdout,"    dependents              %d\n",tape_stats[NUM_DEPENDENTS]);
156    fprintf(stdout,"    operations              %d\n",tape_stats[NUM_OPERATIONS]);
157    fprintf(stdout,"    operations buffer size  %d\n",tape_stats[OP_BUFFER_SIZE]);
158    fprintf(stdout,"    locations buffer size   %d\n",tape_stats[LOC_BUFFER_SIZE]);
159    fprintf(stdout,"    constants buffer size   %d\n",tape_stats[VAL_BUFFER_SIZE]);
160    fprintf(stdout,"    maxlive                 %d\n",tape_stats[NUM_MAX_LIVES]);
161    fprintf(stdout,"    valstack size           %d\n\n",tape_stats[TAY_STACK_SIZE]);
162
163    /*--------------------------------------------------------------------------*/
164    int itu = 8-n;
165    itu = itu*itu*itu*itu;
166    itu = itu > 0 ? itu : 1;
167    double raus;
168
169    /*--------------------------------------------------------------------------*/
170    double t10 = myclock();                             /* 1. time (original) */
171    for (it = 0; it < itu; it++)
172        raus = pdet(n,m-1);
173    double t11 = myclock();
174    double rtu = itu/(t11-t10);
175
176    double* B = new double[n2];
177    double* detaut = new double[1];
178
179    /*--------------------------------------------------------------------------*/
180    double t40 = myclock();                      /* 4. time (forward no keep) */
181    for (it = 0; it < itu; it++)
182        forward(tag,1,n2,0,a,detaut);
183    double t41 = myclock();
184
185    /*--------------------------------------------------------------------------*/
186    double t20 = myclock();                         /* 2. time (forward+keep) */
187    for(it = 0; it < itu; it++)
188        forward(tag,1,n2,1,a,detaut);
189    double t21 = myclock();
190    // fprintf(stdout,"\n %f =? %f should be the same \n",detout,*detaut);
191
192    double u[1];
193    u[0] = 1.0;
194
195    /*--------------------------------------------------------------------------*/
196    double t30 = myclock();                              /* 3. time (reverse) */
197    for (it = 0; it < itu; it++)
198        reverse(tag,1,n2,0,u,B);
199    double t31 = myclock();
200
201    /*--------------------------------------------------------------------------*/
202    /* output of results */
203    // optional generation of tape_doc.tex
204    // tape_doc(tag,1,n2,a,detaut);
205    fprintf(stdout,"\n first base? :   \n");
206    for (i=0; i<n; i++) {
207        adouble sum = 0;
208        adouble* pt;
209        pt = A[i];
210        for (j=0; j<n; j++)
211            sum += (*pt++)*B[j];
212        fprintf(stdout,"%E ",sum.value());
213    }
214    fprintf(stdout,"\n\n times for ");
215    fprintf(stdout,"\n tracing          : \t%E",(t01-t00)*rtu);
216    fprintf(stdout," units \t%E    seconds",(t01-t00));
217    fprintf(stdout,"\n forward (no keep): \t%E",(t41-t40)*rtu/itu);
218    fprintf(stdout," units \t%E    seconds",(t41-t40)/itu);
219    fprintf(stdout,"\n forward + keep   : \t%E",(t21-t20)*rtu/itu);
220    fprintf(stdout," units \t%E    seconds",(t21-t20)/itu);
221    fprintf(stdout,"\n reverse          : \t%E",(t31-t30)*rtu/itu);
222    fprintf(stdout," units \t%E    seconds\n",(t31-t30)/itu);
223
224    return 1;
225}
Note: See TracBrowser for help on using the repository browser.