1 /*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File: detexam.cpp
4 Revision: $Id: detexam.cpp 511 2014-05-14 13:01:01Z kulshres $
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/adolc.h>
19 #include "../clock/myclock.h"
20
21
22 /****************************************************************************/
23 /* DOUBLE ROUTINE */
24 int n,it;
25 double** PA;
pdet(int k,int m)26 double 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 */
63 adouble** A;
64 adouble zero = 0;
det(int k,int m)65 adouble 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 */
main()102 int 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 size_t tape_stats[STAT_SIZE];
152 tapestats(tag,tape_stats);
153
154 fprintf(stdout,"\n independents %zu\n",tape_stats[NUM_INDEPENDENTS]);
155 fprintf(stdout," dependents %zu\n",tape_stats[NUM_DEPENDENTS]);
156 fprintf(stdout," operations %zu\n",tape_stats[NUM_OPERATIONS]);
157 fprintf(stdout," operations buffer size %zu\n",tape_stats[OP_BUFFER_SIZE]);
158 fprintf(stdout," locations buffer size %zu\n",tape_stats[LOC_BUFFER_SIZE]);
159 fprintf(stdout," constants buffer size %zu\n",tape_stats[VAL_BUFFER_SIZE]);
160 fprintf(stdout," maxlive %zu\n",tape_stats[NUM_MAX_LIVES]);
161 fprintf(stdout," valstack size %zu\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 }
226