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

Last change on this file since 638 was 638, checked in by kulshres, 4 years ago

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

The following commits were merged:

commit dd665e6cf8a03b18c4ff44aa5cf7bced713af86a
Author: franke <mirko.franke@…>
Date: Wed Nov 11 11:15:25 2015 +0100

free memory in lie_gradientcv and lie_covector

commit 122d7f6eaba1127af9ba79f08a7e9cdeee6a486e
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue Oct 20 12:53:40 2015 +0200

skip_tapefile_cleanup() is a C function not C++

commit 516f5daf18d4e2825c2c927be92e9982debfcb68
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri Sep 18 15:10:11 2015 +0200

correct pointer free'ing in fortran interfaces and drivers

it seems no one uses this, otherwise people would have seen
segmentation faults all over the place

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

File size: 13.4 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        //prepare output for hos_forward
251        for (i = 0; i < n; i++)
252        {
253                if (d > 1) 
254                {
255                        for (k = d; k > 0; k--)
256                        {
257                                X[i][k] = X[i][k-1];
258                        };
259                        X[i][0] = x[i];
260                };
261        };
262
263        for (i=0; i<m; i++) 
264        {
265                for (k = d; k > 0; k--)
266                {
267                        Y[i][k] = Y[i][k-1];
268                };
269                Y[i][0] = y[i];
270        };
271
272        // prepare output for lie_Scalar
273        for(j=0;j<m;j++)
274        {
275                for (i = 0; i <= d; i++)
276                {
277                        result[j][i] = Fak*Y[j][i];
278                        Fak = Fak*(i+1);
279                }
280                Fak=1.0;
281        }
282
283        myfree2(X);
284        myfree2(Y);
285        myfree1(x);
286        myfree1(y);
287
288        return -1;
289}
290
291
292/** Lie derivative of scalar field h : D -> R^m along vector field f : D -> R^n
293 *  \param Tape_F tape identification of vector field f
294 *  \param Tape_H tape identification of scalar field h
295 *  \param n      number of independent variables n and m = 1
296 *  \param x0     values of independent variables x0 (dimension [n])
297 *  \param d      highest derivative degree d
298 *  \retval result resulting Lie derivatives of a scalar field (dimension [d+1])
299 */
300int lie_scalarc(short Tape_F, short Tape_H, short n, double* x0, short d, double* result) 
301{
302        int rc= -1;
303        short m = 1, i=0;
304        double** Temp = myalloc2(m, d+1);
305
306        rc = lie_scalarcv(Tape_F, Tape_H, n, m, x0, d, Temp);
307
308        for (i = 0; i <= d; i++)
309        {
310                result[i] = Temp[0][i];
311        }
312
313        myfree2(Temp);     
314
315        return rc;
316}
317
318
319/** Calculates the jacobians of the Lie derivatives of scalar fields h : D -> R^m
320 *  \param Tape_F tape identification of vector field f
321 *  \param Tape_H tape identification of vector field h
322 *  \param n      number of independent variables n
323 *  \param m      number of dependent variables m
324 *  \param x0     values of independent variables x0 (dimension [n])
325 *  \param d      highest derivative degree d
326 *  \retval result resulting jacobians of Lie derivatives of vectorial scalar fields (dimension [m][n][d+1])
327 */
328int lie_gradientcv(short Tape_F, short Tape_H, short n, short m, double* x0, short d, double*** result) 
329{
330        double **X=myalloc2(n,d+1);
331        double **Y=myalloc2(m,d+1);
332        double ***Pc=myalloc3(m,n,d+1);
333        double ***A=myalloc3(n,n,d);
334        double ***B=myalloc3(n,n,d);
335        double ***D=myalloc3(m,n,d+1);
336
337        double* x = myalloc1(n);
338        double* y = myalloc1(m);
339        double* xp = myalloc1(n);
340        double* yp = myalloc1(m);
341
342        static int depax_m,depax_n;
343        static double** In;
344        static double** Im;
345
346        int i=0, j=0, l=0, k=0, rc=-1;
347        double Fak=1.0;
348
349        for (i = 0; i < n; i++) 
350                X[i][0] = x0[i];
351
352        forodec(Tape_F, n, 1.0, 0, d, X);
353
354        if (n adolc_compsize depax_n) {
355                if (depax_n)
356                        myfreeI2(depax_n,In);
357                In = myallocI2(depax_n = n);
358        }
359        if (m adolc_compsize depax_m) {
360                if (depax_m)
361                        myfreeI2(depax_m,Im);
362                Im = myallocI2(depax_m = m);
363        }
364
365        hov_reverse(Tape_F,n,n,d-1,n,In,A,0);// explanation in interfaces.cpp
366        accodec(n, 1.0, d-1, A, B, 0);           // explanation in odedrivers.c
367
368        //prepare for input
369        for (i=0; i<n; i++) 
370        {
371                x[i] = X[i][0];
372                if (d == 1)
373                {
374                        xp[i] = X[i][1];
375                }
376                else
377                {
378                        for (k=0; k<d; k++)
379                        {
380                                X[i][k] = X[i][k+1];
381                        }
382                }
383        }
384
385        hos_forward(Tape_H, m, n, d, d+1, x, X, y, Y);
386
387        //prepare output for hos_forward
388        for (i=0; i<n; i++)
389        {
390                if (d > 1) 
391                {
392                        for (k=d; k>0; k--)
393                        {
394                                X[i][k] = X[i][k-1];
395                        }
396                        X[i][0] = x[i];
397                }
398        }
399
400        for (i=0; i<m; i++) 
401        {
402                if (d == 1)
403                {
404                        Y[i][1] = yp[i];
405                }
406                else
407                {
408                        for (k=d; k>0; k--)
409                        {
410                                Y[i][k] = Y[i][k-1];
411                        }
412                        Y[i][0] = y[i];
413                }
414        }
415
416        hov_reverse(Tape_H, m, n, d, m, Im, Pc, 0); 
417        accodeout(m, n, d, B, Pc, D);
418
419        for(l=0; l<m; l++)
420        {
421                for(i=0; i<n; i++)
422                {
423                        Fak=1.0;
424                        for(j=0; j<=d; j++)
425                        {
426                                result[l][i][j]=Fak*D[l][i][j];
427                                Fak=Fak*(j+1);
428                        }
429                }
430        }
431
432        myfree2(X);
433        myfree2(Y);
434        myfree3(Pc);
435        myfree3(A);
436        myfree3(B);
437        myfree3(D);
438
439        myfree1(x);
440        myfree1(y);
441        myfree1(xp);
442        myfree1(yp);
443
444        return rc;
445}
446
447
448
449
450/** Computes the gradients of the Lie derivatives of a scalar field h : D -> R
451 *
452 *  \param Tape_F tape identification of vector field f
453 *  \param Tape_H tape identification of vector field h
454 *  \param n      number of independent variables n
455 *  \param x0     values of independent variables x0 (dimension [n])
456 *  \param d      highest derivative degree d
457 *  \retval result resulting jacobians of Lie derivatives of a scalar field (dimension [n][d+1])
458 */
459int lie_gradientc(short Tape_F, short Tape_H, short n, double* x0, short d,     double** result) 
460{
461        int rc= -1;
462        short m = 1, i=0, j=0;
463        double*** Temp = myalloc3(m, n, d+1);
464
465        rc = lie_gradientcv(Tape_F, Tape_H, n, m, x0, d, Temp);
466
467        for(i=0; i<n; i++)
468                for(j=0; j<=d; j++)
469                {
470                        result[i][j]=Temp[0][i][j];
471
472                }
473
474                myfree3(Temp);
475
476                return rc;
477}
478
479       
480
481/** Computes the Lie derivatives of the covector field w : D -> (R^m)* along the vector field f : D -> R^n
482 *
483 *  \param Tape_F tape identification of vector field f
484 *  \param Tape_W tape identification of covector field h
485 *  \param n      number of independent variables n
486 *  \param x0     values of independent variables x0 (dimension [n])
487 *  \param d      highest derivative degree d
488 *  \retval result resulting Lie derivatives of a covector field (dimension [n][d+1])
489 */
490int lie_covector( short int Tape_F, short int Tape_W, short int n, double* x0, short int d, double** result)
491{
492        int m=n;                   
493        double** X = myalloc2(n, d+1);   // Taylorcoeff. expansion x(t)
494        double** Y = myalloc2(m, d+1);   // Taylorcoeff. expansion y(t)
495
496        double***A = myalloc3(n, n, d);   
497        double***B = myalloc3(n, n, d+1);
498
499        double* x = myalloc1(n);
500        double* y = myalloc1(m);
501        double* xp = myalloc1(n);
502        double* yp = myalloc1(m);
503
504        int i=0, k=0;
505
506        static int depax_m,depax_n;
507        static double** In;
508        static double** Im;
509
510        for (i = 0; i < n; i++) X[i][0] = x0[i];
511
512        forodec(Tape_F, n, 1.0, 0, d, X);
513
514        if (n adolc_compsize depax_n) {
515                if (depax_n)
516                        myfreeI2(depax_n,In);
517                In = myallocI2(depax_n = n);
518        }
519        if (m adolc_compsize depax_m) {
520                if (depax_m)
521                        myfreeI2(depax_m,Im);
522                Im = myallocI2(depax_m = m);
523        }
524
525        hov_reverse(Tape_F,n,n,d-1,n,In,A,0);   // explanation in interfaces.cpp
526
527        //prepare for input
528        for (i=0; i<n; i++) {
529                x[i] = X[i][0];
530                if (d == 1)
531                        xp[i] = X[i][1];
532                else
533                        for (k=0; k<d; k++)
534                                X[i][k] = X[i][k+1];
535        }
536
537        hos_forward(Tape_W, m, n, d, d+1, x, X, y, Y);
538
539        //prepare output for hos_forward
540        for (i=0; i<n; i++)
541                if (d > 1) {
542                        for (k=d; k>0; k--)
543                                X[i][k] = X[i][k-1];
544                        X[i][0] = x[i];
545                }
546
547        for (i=0; i<m; i++) {
548                if (d == 1)
549                        Y[i][1] = yp[i];
550                else
551                        for (k=d; k>0; k--)
552                                Y[i][k] = Y[i][k-1];
553                Y[i][0] = y[i];
554        }
555
556        accodec(n, 1.0, d-1, A, B, 0);      // explanation in odedrivers.c
557        acccov(n, d, B, Y, result); 
558
559        myfree2(X);
560        myfree2(Y);
561        myfree3(A);
562        myfree3(B);
563
564        myfree1(x);
565        myfree1(y);
566        myfree1(xp);
567        myfree1(yp);
568
569        return -1;
570}
571
572
573
574
575/** Calculates the iterated Lie derivatives (Lie brackets) of the vector field g : D -> R^n along the vector field f : D -> R^n
576 *
577 *  \param Tape_F tape identification of vector field f
578 *  \param Tape_G tape identification of vector field g
579 *  \param n      number of independent variables n
580 *  \param x0     values of independent variables x0 (dimension [n])
581 *  \param d      highest derivative degree d
582 */
583int lie_bracket(short int Tape_F, short int Tape_G, short int n, double* x0, short int d, double** result)
584{
585        int m = n;                     
586        double  **= myalloc2(n, d+2);   // Taylorcoeff. expansion x(t)
587        double  **= myalloc2(m, d+2);   // Taylorcoeff. expansion y(t)
588        double ***= myalloc3(n, n, d+1);   
589        double ***Xs = myalloc3(n, n, d+1);
590
591        double* y = myalloc1(m);
592
593        int i, k;
594
595        //static identity matrix for hov_reverse
596        static int      depax_n = 0;
597        static double** In      = NULL;
598
599
600        for (i = 0; i < n; i++) 
601        {
602                X[i][0] = x0[i];
603        };
604        forodec(Tape_F, n, 1.0, 0, d+1, X);
605
606        // for hov_reverse
607        if (n > depax_n) 
608        {
609                if (depax_n)
610                {
611                        myfreeI2(depax_n, In);
612                };
613                In = myallocI2(depax_n = n);
614        }
615        hov_reverse(Tape_F, n, n, d, n, In, A, 0);
616
617        //prepare for input
618        for (i = 0; i < n; i++) 
619        {
620                for (k = 0; k < d; k++)
621                {
622                        X[i][k] = X[i][k+1];
623                }
624        }
625        hos_forward(Tape_G, m, n, d+1, d+2, x0, X, y, Y);
626
627        //postprocess output of hos_forward
628        for (i = 0; i < m; i++) 
629        {
630                for (k = d; k > 0; k--)
631                {
632                        Y[i][k] = Y[i][k-1];
633                }
634                Y[i][0] = y[i];
635        }
636
637        accadj(n, d, A, Xs); 
638        accbrac(n, d, Xs, Y, result); 
639
640        myfree1(y);
641        myfree2(X);
642        myfree2(Y);
643        myfree3(A);
644        myfree3(Xs);
645       
646        return -1;
647};
648
649
Note: See TracBrowser for help on using the repository browser.