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