source: trunk/ADOL-C/src/lie/adolc_lie_c.c

Last change on this file was 708, checked in by kulshres, 3 years ago

Merge branch 'master' of 'gitclone' into 'svn'

The following changes have been merged:

commit e2291bde44a282a133894b0db350aeb0b92a87db
Author: Mladen Banovic <mladenbanovic2705@…>
Date: Fri Jul 8 10:15:51 2016 +0200

Add methods getNumLiveVar and getNumDir in adtl.h, change counter type in FOR_I_EQ_0_LT_NUMDIR macro to size_t (instead of int). Update chunk size of BOOST pool in adouble_tl.cpp according to adouble::numDir.

commit 2ffb294465b973bfd4bf1f73d84478f8233c0d2f
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Thu Jun 23 12:32:14 2016 +0200

implement missing ref_eq_mult_p und ref_eq_min_p in ho_rev.c

somehow these were left out when parameters were being implemented.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 8cf0e5c1bd36f1dcf3be72cd67de631b2e1d0ee6
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Thu Jun 23 12:31:04 2016 +0200

make sure the result is the last locint written in trace for each operation

since we're trying to generate ascii traces in the future, we'll need this
convention that the last location is the result, and previous locations
are arguments. This has been the case for almost all operations anyway
except for a few new one's that I wrote without keeping this in mind.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 9ae0ff220f37463f2ed85cafc8a626c24e472f2f
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue Jun 21 14:16:27 2016 +0200

on some compilers newer boost interferes with AC_FUNC_MALLOC test

so do AC_FUNC_MALLOC and AC_FUNC_REALLOC as usual and check for boost
library later.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit b746f620772cc8cce53e8f350adc6281279caf72
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Jun 20 15:32:22 2016 +0200

make Klaus Röbenack's name UTF-8 instead of ISO-8859-1

These are the only places where we're not simple ASCII or UTF-8 already

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 1171aa3961b5eb46a5d2ee64751c02a393a8a6f5
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri Jun 17 10:42:39 2016 +0200

correct short_ref document about include file

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 2c6b2aac2ef04431ece2c6ff80e574aa2e58814b
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri Jun 17 10:40:34 2016 +0200

correct error message to new semantics

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 506cde73451740bf0a15eff7d4abb158ee719ab0
Author: mflehmig <martin.flehmig@…>
Date: Fri Jun 17 10:14:26 2016 +0200

Fixed include of ADOL-C header.

ADOL-C header was included in old fashion (without adolc directory) for this example.

commit 2a023d3281d3d6d9824bad724a5768e3ee2fff94
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Thu Jun 16 13:50:39 2016 +0200

Try to use boost::pool for allocating advals in traceless vector mode

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

commit 80f1e2019ac1faab96fe06f3e9da47efcc1bcd23
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon May 23 15:13:22 2016 +0200

correct a reference in doc and rebuild

commit d7ab5283afe58bacb2e8739d72ede4e17f4c8081
Author: Mladen Banovic <mladenbanovic2705@…>
Date: Fri May 20 16:42:13 2016 +0200

Update section 7 of adolc-manual related to the Traceless forward differentiation.

commit bedb8e36f959c5272e4610fe504acc83208e5e9d
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue May 17 16:09:36 2016 +0200

macro name correction

commit 92ff596a0331776901df7f172ca347572e3daafd
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue May 17 15:56:17 2016 +0200

Add a warning about using static build of ADOL-C

static build of ADOL-C does not call constructors
for internal global objects, thereby causing
segmentation faults.

Signed-off-by: Kshitij Kulshreshtha <kshitij@…>

File size: 12.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     lie/adolc_lie_c.c
4 Revision: $Id$
5 Contents: Implementation of functions for computation of Lie derivatives
6 
7
8 Copyright (c) Siquian Wang, Klaus Röbenack, Jan Winkler, Mirko Franke
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#include <adolc/adolc.h>
16#include <adolc/adouble.h>
17#include <adolc/lie/drivers.h>
18#include "dvlparms.h"
19
20
21
22
23/** Computes the total derivative of the output
24 *
25 *  @param p Number of rows of C (number of outputs)
26 *  @param n Number of columns of B (number of inputs)
27 *  @param deg Order of derivative (d)
28 *  @retval ***B Total derivative dx(j+1)/dx0
29 *  @param ***C Partial derivative of the output mapping
30 *  @param ***D Total derivative of the output mapping
31 */
32void accodeout (int p, int n, int deg, double ***B,     double ***C, double ***D)
33{
34        int i, j, k, l, ip, jp, kp;
35
36        // D_k:=C_k (1. summation step)
37        for (k=0; k <= deg; k++)
38        {
39                for (j=0; j < p; j++)
40                        for (i=0; i < n; i++)
41                                D[j][i][k]=C[j][i][k];
42               
43                // add sum if necessary
44                if (k >= 1)
45                {
46                        for (l=1; l<=k; l++)
47                                for (jp=0; jp<p; jp++)
48                                        for (ip=0; ip<n; ip++)
49                                        {
50                                                double x=0.0;
51                                                for (kp=0; kp<n; kp++)
52                                                {
53                                                        x+=C[jp][kp][k-l]*B[kp][ip][l-1];
54                                                };
55                                                D[jp][ip][k]+=x;
56                                        };
57                };
58        };
59};
60
61       
62
63/** Helper function for calculation of the Lie derivatives of a covector field
64 *
65 *  @param n Number of rows = number of columns
66 *  @param d Order of derivative (d)
67 *  @param ***B Total derivative dx(j+1)/dx0 (B[n][n][d])
68 *  @param **C Taylor coefficients of the output mapping (C[n][d])
69 *  @retval **D Lie derivative along covector field (D[n][d])
70 */
71void acccov(int n, int d, double ***B, double **C, double **D)
72{
73        int i, k, l, ip, kp;
74
75        // factorial     
76        int Fak = 1;
77       
78        // D_k:=Fak*C_k (1. summation step)
79        for (k = 0; k <= d; k++)
80        {
81                if(k == 0)
82                        Fak = 1;
83                else
84                        Fak = Fak*k;
85
86                for (i = 0; i < n; i++)
87                        D[i][k]=Fak*C[i][k];
88
89                // add sum if necessary
90                if (k>=1)
91                {
92                        double x;
93                        for (l = 1; l <= k; l++)
94                        {
95                                for (ip = 0; ip < n; ip++)
96                                {
97                                        x = 0.0;
98                                        for (kp = 0; kp < n; kp++)
99                                        {
100                                                x+=C[kp][k-l]*B[kp][ip][l-1];
101                                        };
102                                        D[ip][k]+= Fak*x;
103                                };
104                        };
105                };
106        };
107};
108
109       
110
111
112/** Helper function for calculation of Lie-brackets (solution of the adjoint variational equation)
113 *
114 *  @param n Number of columns of B (number of inputs)
115 *  @param deg Order of derivative
116 *  @param ***A Total derivative of A
117 *  @retval ***Bs Solution of adjoint variational equation
118 */
119void accadj(int n, int deg, double ***A, double ***Bs)
120{
121        int i, j, k, l, ip, jp, kp;
122
123        // (1. summation step)
124        for (k = 0; k <= deg; k++)
125        {
126                for (j = 0; j < n; j++)
127                {
128                        for (i = 0; i < n; i++)
129                        {
130                                Bs[j][i][k] = -A[i][j][k]/(k+1);
131                        }
132                }
133                       
134                // add sum if necessary
135                if (k >= 1)
136                {
137                        double x = 0.0;
138                        for (l = 1; l <= k; l++)
139                        {
140                                for (jp = 0; jp < n; jp++)
141                                {
142                                        for (ip = 0; ip < n; ip++)
143                                        {
144                                                x = 0.0;
145                                                for (kp = 0; kp < n; kp++)
146                                                {
147                                                        x+=A[kp][jp][k-l]*Bs[kp][ip][l-1];
148                                                };
149                                                Bs[jp][ip][k] -= x/(k+1);
150                                        };
151                                };
152                        };
153                };
154        };
155};
156
157
158
159/** Calculates the Lie-derivative along a co vector field
160 *
161 * @param n Number of rows and columns
162 * @param d Order of derivative (d)
163 * @param Bs Solution of adjoint variational equation (Bs[n][n][d])
164 * @param b Taylor-coefficients of output mapping (b[n][d])
165 * @retval result Lie derivative along co-vector field (result[n][d])
166 */
167void accbrac(int n,     int d, double ***Bs, double **b, double **result) 
168{
169        int i, j, k, l, jp, kp;
170
171        // factorial     
172        int Fak = 1;
173
174        // formula 3.58
175        for (k = 0; k <= d; k++)
176        {
177                if(k == 0)
178                        Fak = 1;
179                else
180                        Fak = Fak*k;
181
182                for (j = 0; j < n; j++)
183                {
184                        for (i = 0; i < n; i++)
185                        {
186                                result[i][k] = Fak*b[i][k];
187                        }
188
189                        if(k >= 1)
190                        {
191                                double x;
192                                for (l = 1; l <= k; l++)
193                                {
194                                        for (jp = 0; jp < n; jp++)
195                                        {
196                                                x = 0.0;
197                                                for (kp = 0; kp < n; kp++)
198                                                {
199                                                        x+=Bs[kp][jp][l-1]*b[kp][k-l];
200                                                }
201                                                result[jp][k]+=Fak*x;
202                                        };
203                                };
204                        };
205                };
206        };
207};
208
209       
210
211
212
213
214/** Computes the Lie derivative of smooth map h : D -> R^m along a vector field f : D -> R^n
215 * \param Tape_F tape identification of vector field f
216 * \param Tape_H tape identification of vector field h
217 * \param n      number of independent variables n
218 * \param m          number of dependent variables m
219 * \param x0     values of independent variables x0 (dimension [n])
220 * \param d      highest derivative degree d
221 * \param result resulting Lie derivatives of vectorial scalar fields (dimension [m][d+1])
222 */
223int lie_scalarcv(short Tape_F, short Tape_H, short n, short m, double* x0, short d,     double** result) 
224{
225        double** X = myalloc2(n, d+1);   // Taylorcoeff. expansion x(t)
226        double** Y = myalloc2(m, d+1);   // Taylorcoeff. expansion y(t)
227        double*  x = myalloc1(n);
228        double*  y = myalloc1(m);
229        int i=0, j=0, k=0;
230        double Fak = 1.0;
231
232        for (i = 0; i < n; i++) 
233        {
234                X[i][0] = x0[i];
235        };
236       
237        //see odedrivers
238        forodec(Tape_F, n, 1.0, 0, d, X);
239
240        //prepare for input
241        for (i = 0; i < n; i++) 
242        {
243                x[i] = X[i][0];
244                for (k = 0; k < d; k++)
245                        X[i][k] = X[i][k+1];
246        }
247
248        hos_forward(Tape_H, m, n, d, 0, x, X, y, Y);
249
250        //postprocess output of hos_forward
251        for (i=0; i<m; i++) 
252        {
253                for (k = d; k > 0; k--)
254                {
255                        Y[i][k] = Y[i][k-1];
256                };
257                Y[i][0] = y[i];
258        };
259
260        // prepare output for lie_Scalar
261        for(j=0;j<m;j++)
262        {
263                for (i = 0; i <= d; i++)
264                {
265                        result[j][i] = Fak*Y[j][i];
266                        Fak = Fak*(i+1);
267                }
268                Fak=1.0;
269        }
270
271        myfree2(X);
272        myfree2(Y);
273        myfree1(x);
274        myfree1(y);
275
276        return -1;
277}
278
279
280/** Lie derivative of scalar field h : D -> R^m along vector field f : D -> R^n
281 *  \param Tape_F tape identification of vector field f
282 *  \param Tape_H tape identification of scalar field h
283 *  \param n      number of independent variables n and m = 1
284 *  \param x0     values of independent variables x0 (dimension [n])
285 *  \param d      highest derivative degree d
286 *  \retval result resulting Lie derivatives of a scalar field (dimension [d+1])
287 */
288int lie_scalarc(short Tape_F, short Tape_H, short n, double* x0, short d, double* result) 
289{
290        int rc= -1;
291        short m = 1, i=0;
292        double** Temp = myalloc2(m, d+1);
293
294        rc = lie_scalarcv(Tape_F, Tape_H, n, m, x0, d, Temp);
295
296        for (i = 0; i <= d; i++)
297        {
298                result[i] = Temp[0][i];
299        }
300
301        myfree2(Temp);     
302
303        return rc;
304}
305
306
307/** Calculates the jacobians of the Lie derivatives of scalar fields h : D -> R^m
308 *  \param Tape_F tape identification of vector field f
309 *  \param Tape_H tape identification of vector field h
310 *  \param n      number of independent variables n
311 *  \param m      number of dependent variables m
312 *  \param x0     values of independent variables x0 (dimension [n])
313 *  \param d      highest derivative degree d
314 *  \retval result resulting jacobians of Lie derivatives of vectorial scalar fields (dimension [m][n][d+1])
315 */
316int lie_gradientcv(short Tape_F, short Tape_H, short n, short m, double* x0, short d, double*** result) 
317{
318        double **X=myalloc2(n,d+1);
319        double **Y=myalloc2(m,d+1);
320        double ***Pc=myalloc3(m,n,d+1);
321        double ***A=myalloc3(n,n,d);
322        double ***B=myalloc3(n,n,d);
323        double ***D=myalloc3(m,n,d+1);
324
325        double* x = myalloc1(n);
326        double* y = myalloc1(m);
327        double* xp = myalloc1(n);
328        double* yp = myalloc1(m);
329
330        static int depax_m,depax_n;
331        static double** In;
332        static double** Im;
333
334        int i=0, j=0, l=0, k=0, rc=-1;
335        double Fak=1.0;
336
337        for (i = 0; i < n; i++) 
338                X[i][0] = x0[i];
339
340        forodec(Tape_F, n, 1.0, 0, d, X);
341
342        if (n adolc_compsize depax_n) {
343                if (depax_n)
344                        myfreeI2(depax_n,In);
345                In = myallocI2(depax_n = n);
346        }
347        if (m adolc_compsize depax_m) {
348                if (depax_m)
349                        myfreeI2(depax_m,Im);
350                Im = myallocI2(depax_m = m);
351        }
352
353        hov_reverse(Tape_F,n,n,d-1,n,In,A,0);// explanation in interfaces.cpp
354        accodec(n, 1.0, d-1, A, B, 0);           // explanation in odedrivers.c
355
356        //prepare for input
357        for (i=0; i<n; i++) 
358        {
359                x[i] = X[i][0];
360                if (d == 1)
361                {
362                        xp[i] = X[i][1];
363                }
364                else
365                {
366                        for (k=0; k<d; k++)
367                        {
368                                X[i][k] = X[i][k+1];
369                        }
370                }
371        }
372
373        hos_forward(Tape_H, m, n, d, d+1, x, X, y, Y);
374
375        hov_reverse(Tape_H, m, n, d, m, Im, Pc, 0); 
376        accodeout(m, n, d, B, Pc, D);
377
378        for(l=0; l<m; l++)
379        {
380                for(i=0; i<n; i++)
381                {
382                        Fak=1.0;
383                        for(j=0; j<=d; j++)
384                        {
385                                result[l][i][j]=Fak*D[l][i][j];
386                                Fak=Fak*(j+1);
387                        }
388                }
389        }
390
391        myfree2(X);
392        myfree2(Y);
393        myfree3(Pc);
394        myfree3(A);
395        myfree3(B);
396        myfree3(D);
397
398        myfree1(x);
399        myfree1(y);
400        myfree1(xp);
401        myfree1(yp);
402
403        return rc;
404}
405
406
407
408
409/** Computes the gradients of the Lie derivatives of a scalar field h : D -> R
410 *
411 *  \param Tape_F tape identification of vector field f
412 *  \param Tape_H tape identification of vector field h
413 *  \param n      number of independent variables n
414 *  \param x0     values of independent variables x0 (dimension [n])
415 *  \param d      highest derivative degree d
416 *  \retval result resulting jacobians of Lie derivatives of a scalar field (dimension [n][d+1])
417 */
418int lie_gradientc(short Tape_F, short Tape_H, short n, double* x0, short d,     double** result) 
419{
420        int rc= -1;
421        short m = 1, i=0, j=0;
422        double*** Temp = myalloc3(m, n, d+1);
423
424        rc = lie_gradientcv(Tape_F, Tape_H, n, m, x0, d, Temp);
425
426        for(i=0; i<n; i++)
427                for(j=0; j<=d; j++)
428                {
429                        result[i][j]=Temp[0][i][j];
430
431                }
432
433                myfree3(Temp);
434
435                return rc;
436}
437
438       
439
440/** Computes the Lie derivatives of the covector field w : D -> (R^m)* along the vector field f : D -> R^n
441 *
442 *  \param Tape_F tape identification of vector field f
443 *  \param Tape_W tape identification of covector field h
444 *  \param n      number of independent variables n
445 *  \param x0     values of independent variables x0 (dimension [n])
446 *  \param d      highest derivative degree d
447 *  \retval result resulting Lie derivatives of a covector field (dimension [n][d+1])
448 */
449int lie_covector( short int Tape_F, short int Tape_W, short int n, double* x0, short int d, double** result)
450{
451        int m=n;                   
452        double** X = myalloc2(n, d+1);   // Taylorcoeff. expansion x(t)
453        double** Y = myalloc2(m, d+1);   // Taylorcoeff. expansion y(t)
454
455        double***A = myalloc3(n, n, d);   
456        double***B = myalloc3(n, n, d+1);
457
458        double* x = myalloc1(n);
459        double* y = myalloc1(m);
460        double* xp = myalloc1(n);
461        double* yp = myalloc1(m);
462
463        int i=0, k=0;
464
465        static int depax_m,depax_n;
466        static double** In;
467        static double** Im;
468
469        for (i = 0; i < n; i++) X[i][0] = x0[i];
470
471        forodec(Tape_F, n, 1.0, 0, d, X);
472
473        if (n adolc_compsize depax_n) {
474                if (depax_n)
475                        myfreeI2(depax_n,In);
476                In = myallocI2(depax_n = n);
477        }
478        if (m adolc_compsize depax_m) {
479                if (depax_m)
480                        myfreeI2(depax_m,Im);
481                Im = myallocI2(depax_m = m);
482        }
483
484        hov_reverse(Tape_F,n,n,d-1,n,In,A,0);   // explanation in interfaces.cpp
485
486        //prepare for input
487        for (i=0; i<n; i++) {
488                x[i] = X[i][0];
489                if (d == 1)
490                        xp[i] = X[i][1];
491                else
492                        for (k=0; k<d; k++)
493                                X[i][k] = X[i][k+1];
494        }
495
496        hos_forward(Tape_W, m, n, d, d+1, x, X, y, Y);
497
498        //postprocess output of hos_forward
499        for (i=0; i<m; i++) {
500                if (d == 1)
501                        Y[i][1] = yp[i];
502                else
503                        for (k=d; k>0; k--)
504                                Y[i][k] = Y[i][k-1];
505                Y[i][0] = y[i];
506        }
507
508        accodec(n, 1.0, d-1, A, B, 0);      // explanation in odedrivers.c
509        acccov(n, d, B, Y, result); 
510
511        myfree2(X);
512        myfree2(Y);
513        myfree3(A);
514        myfree3(B);
515
516        myfree1(x);
517        myfree1(y);
518        myfree1(xp);
519        myfree1(yp);
520
521        return -1;
522}
523
524
525
526
527/** Calculates the iterated Lie derivatives (Lie brackets) of the vector field g : D -> R^n along the vector field f : D -> R^n
528 *
529 *  \param Tape_F tape identification of vector field f
530 *  \param Tape_G tape identification of vector field g
531 *  \param n      number of independent variables n
532 *  \param x0     values of independent variables x0 (dimension [n])
533 *  \param d      highest derivative degree d
534 */
535int lie_bracket(short int Tape_F, short int Tape_G, short int n, double* x0, short int d, double** result)
536{
537        int m = n;                     
538        double  **= myalloc2(n, d+2);   // Taylorcoeff. expansion x(t)
539        double  **= myalloc2(m, d+2);   // Taylorcoeff. expansion y(t)
540        double ***= myalloc3(n, n, d+1);   
541        double ***Xs = myalloc3(n, n, d+1);
542
543        double* y = myalloc1(m);
544
545        int i, k;
546
547        //static identity matrix for hov_reverse
548        static int      depax_n = 0;
549        static double** In      = NULL;
550
551
552        for (i = 0; i < n; i++) 
553        {
554                X[i][0] = x0[i];
555        };
556        forodec(Tape_F, n, 1.0, 0, d+1, X);
557
558        // for hov_reverse
559        if (n > depax_n) 
560        {
561                if (depax_n)
562                {
563                        myfreeI2(depax_n, In);
564                };
565                In = myallocI2(depax_n = n);
566        }
567        hov_reverse(Tape_F, n, n, d, n, In, A, 0);
568
569        //prepare for input
570        for (i = 0; i < n; i++) 
571        {
572                for (k = 0; k < d; k++)
573                {
574                        X[i][k] = X[i][k+1];
575                }
576        }
577        hos_forward(Tape_G, m, n, d+1, d+2, x0, X, y, Y);
578
579        //postprocess output of hos_forward
580        for (i = 0; i < m; i++) 
581        {
582                for (k = d; k > 0; k--)
583                {
584                        Y[i][k] = Y[i][k-1];
585                }
586                Y[i][0] = y[i];
587        }
588
589        accadj(n, d, A, Xs); 
590        accbrac(n, d, Xs, Y, result); 
591
592        myfree1(y);
593        myfree2(X);
594        myfree2(Y);
595        myfree3(A);
596        myfree3(Xs);
597       
598        return -1;
599};
600
601
Note: See TracBrowser for help on using the repository browser.