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

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

add the new files for the lie driver and example

File size: 13.3 KB
RevLine 
[609]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        return rc;
440}
441
442
443
444
445/** Computes the gradients of the Lie derivatives of a scalar field h : D -> R
446 *
447 *  \param Tape_F tape identification of vector field f
448 *  \param Tape_H tape identification of vector field h
449 *  \param n      number of independent variables n
450 *  \param x0     values of independent variables x0 (dimension [n])
451 *  \param d      highest derivative degree d
452 *  \retval result resulting jacobians of Lie derivatives of a scalar field (dimension [n][d+1])
453 */
454int lie_gradientc(short Tape_F, short Tape_H, short n, double* x0, short d,     double** result) 
455{
456        int rc= -1;
457        short m = 1, i=0, j=0;
458        double*** Temp = myalloc3(m, n, d+1);
459
460        rc = lie_gradientcv(Tape_F, Tape_H, n, m, x0, d, Temp);
461
462        for(i=0; i<n; i++)
463                for(j=0; j<=d; j++)
464                {
465                        result[i][j]=Temp[0][i][j];
466
467                }
468
469                myfree3(Temp);
470
471                return rc;
472}
473
474       
475
476/** Computes the Lie derivatives of the covector field w : D -> (R^m)* along the vector field f : D -> R^n
477 *
478 *  \param Tape_F tape identification of vector field f
479 *  \param Tape_W tape identification of covector field h
480 *  \param n      number of independent variables n
481 *  \param x0     values of independent variables x0 (dimension [n])
482 *  \param d      highest derivative degree d
483 *  \retval result resulting Lie derivatives of a covector field (dimension [n][d+1])
484 */
485int lie_covector( short int Tape_F, short int Tape_W, short int n, double* x0, short int d, double** result)
486{
487        int m=n;                   
488        double** X = myalloc2(n, d+1);   // Taylorcoeff. expansion x(t)
489        double** Y = myalloc2(m, d+1);   // Taylorcoeff. expansion y(t)
490
491        double***A = myalloc3(n, n, d);   
492        double***B = myalloc3(n, n, d+1);
493
494        double* x = myalloc1(n);
495        double* y = myalloc1(m);
496        double* xp = myalloc1(n);
497        double* yp = myalloc1(m);
498
499        int i=0, k=0;
500
501        static int depax_m,depax_n;
502        static double** In;
503        static double** Im;
504
505        for (i = 0; i < n; i++) X[i][0] = x0[i];
506
507        forodec(Tape_F, n, 1.0, 0, d, X);
508
509        if (n adolc_compsize depax_n) {
510                if (depax_n)
511                        myfreeI2(depax_n,In);
512                In = myallocI2(depax_n = n);
513        }
514        if (m adolc_compsize depax_m) {
515                if (depax_m)
516                        myfreeI2(depax_m,Im);
517                Im = myallocI2(depax_m = m);
518        }
519
520        hov_reverse(Tape_F,n,n,d-1,n,In,A,0);   // explanation in interfaces.cpp
521
522        //prepare for input
523        for (i=0; i<n; i++) {
524                x[i] = X[i][0];
525                if (d == 1)
526                        xp[i] = X[i][1];
527                else
528                        for (k=0; k<d; k++)
529                                X[i][k] = X[i][k+1];
530        }
531
532        hos_forward(Tape_W, m, n, d, d+1, x, X, y, Y);
533
534        //prepare output for hos_forward
535        for (i=0; i<n; i++)
536                if (d > 1) {
537                        for (k=d; k>0; k--)
538                                X[i][k] = X[i][k-1];
539                        X[i][0] = x[i];
540                }
541
542        for (i=0; i<m; i++) {
543                if (d == 1)
544                        Y[i][1] = yp[i];
545                else
546                        for (k=d; k>0; k--)
547                                Y[i][k] = Y[i][k-1];
548                Y[i][0] = y[i];
549        }
550
551        accodec(n, 1.0, d-1, A, B, 0);      // explanation in odedrivers.c
552        acccov(n, d, B, Y, result); 
553
554        return -1;
555}
556
557
558
559
560/** Calculates the iterated Lie derivatives (Lie brackets) of the vector field g : D -> R^n along the vector field f : D -> R^n
561 *
562 *  \param Tape_F tape identification of vector field f
563 *  \param Tape_G tape identification of vector field g
564 *  \param n      number of independent variables n
565 *  \param x0     values of independent variables x0 (dimension [n])
566 *  \param d      highest derivative degree d
567 */
568int lie_bracket(short int Tape_F, short int Tape_G, short int n, double* x0, short int d, double** result)
569{
570        int m = n;                     
571        double  **= myalloc2(n, d+2);   // Taylorcoeff. expansion x(t)
572        double  **= myalloc2(m, d+2);   // Taylorcoeff. expansion y(t)
573        double ***= myalloc3(n, n, d+1);   
574        double ***Xs = myalloc3(n, n, d+1);
575
576        double* y = myalloc1(m);
577
578        int i, k;
579
580        //static identity matrix for hov_reverse
581        static int      depax_n = 0;
582        static double** In      = NULL;
583
584
585        for (i = 0; i < n; i++) 
586        {
587                X[i][0] = x0[i];
588        };
589        forodec(Tape_F, n, 1.0, 0, d+1, X);
590
591        // for hov_reverse
592        if (n > depax_n) 
593        {
594                if (depax_n)
595                {
596                        myfreeI2(depax_n, In);
597                };
598                In = myallocI2(depax_n = n);
599        }
600        hov_reverse(Tape_F, n, n, d, n, In, A, 0);
601
602        //prepare for input
603        for (i = 0; i < n; i++) 
604        {
605                for (k = 0; k < d; k++)
606                {
607                        X[i][k] = X[i][k+1];
608                }
609        }
610        hos_forward(Tape_G, m, n, d+1, d+2, x0, X, y, Y);
611
612        //postprocess output of hos_forward
613        for (i = 0; i < m; i++) 
614        {
615                for (k = d; k > 0; k--)
616                {
617                        Y[i][k] = Y[i][k-1];
618                }
619                Y[i][0] = y[i];
620        }
621
622        accadj(n, d, A, Xs); 
623        accbrac(n, d, Xs, Y, result); 
624
625        myfree1(y);
626        myfree2(X);
627        myfree2(Y);
628        myfree3(A);
629        myfree3(Xs);
630       
631        return -1;
632};
633
634
Note: See TracBrowser for help on using the repository browser.