source: trunk/ADOL-C/src/fixpoint.cpp @ 708

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

merge branch 'master' of 'gitclone' into 'svn'

The following commits have been merged:

commit 48aee4916d2ed907b772dbd1c1d6ce46cb273651
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Mon Aug 10 21:28:49 2015 +0200

modernise configure.ac

disable static library building. this causes more
problems than it solves.

commit 47332811a4c5c27cb884f75792c910c813378ef4
Merge: 0ee77fd 0d4eeec
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Thu Aug 6 22:33:46 2015 +0200

Merge branch 'edf-memory'

This is to reduce memory allocation and copying in ext_diff_fct

commit 0ee77fd33a1d6d55fcc67ad419937b2cb777ed4e
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Wed Aug 5 15:49:33 2015 +0200

Remove empty file from dist

this can be created during compilation

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

commit 51505c34571aa61b4b21ebce6cdf1728ff56ddaa
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Tue Aug 4 17:36:26 2015 +0200

adouble(const adub&) should match operator=(const adub&)

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

commit 03e49097aa0455337647d280cda530064987e6b9
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri Jul 3 11:17:53 2015 +0200

make a define for default contiguous locations

this is not needed during compilation of the library
only during compilation of user-code, if the user
wants to have all adouble* allocations to have
contiguous locations.

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

commit f00cfb5d0dc8a8993581fd8c08dd8c6c5cd23248
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Wed Jul 1 11:56:39 2015 +0200

rename adolc_lie.h to drivers.h

the old name led to an include <adolc/lie/adolc_lie.h>
which looks highly redundant.
new name makes for include <adolc/lie/drivers.h>

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

commit fcf78bf8426a227750a0bcaa32ff65e57ef329b8
Author: franke <mirko.franke@…>
Date: Wed May 20 16:39:16 2015 +0200

added Lie drivers

commit 0d4eeec7b6212aa64c8997db8a511f81b604b3e1
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri Jun 26 14:19:41 2015 +0200

minimise extra memory requirement and copies in ext_diff

This should in theory reduce the amount of memory required
to run an external function with the old interface. It also
reduced some copying operations.

Fingers crossed that we've not broken checkpointing and/or
fixpoint iterations.

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

commit 8811f02a4d4a15946a18f7513ab17dada66509c3
Author: Kshitij Kulshreshtha <kshitij@…>
Date: Fri May 22 12:41:45 2015 +0200

try to streamline data copies in ext_diff_v2

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

  • Property svn:keywords set to Author Date Id Revision
File size: 7.9 KB
Line 
1
2/*----------------------------------------------------------------------------
3 ADOL-C -- Automatic Differentiation by Overloading in C++
4 File:     fixpoint.c
5 Revision: $Id: fixpoint.cpp 608 2015-08-10 20:06:55Z kulshres $
6 Contents: all C functions directly accessing at least one of the four tapes
7           (operations, locations, constants, value stack)
8 
9 Copyright (c) Andreas Kowarz, Sebastian Schlenkrich
10 
11 This file is part of ADOL-C. This software is provided as open source.
12 Any use, reproduction, or distribution of the software constitutes
13 recipient's acceptance of the terms of the accompanying license file.
14
15----------------------------------------------------------------------------*/
16
17#include "taping_p.h"
18#include <adolc/adolc.h>
19#include <adolc/fixpoint.h>
20#include "dvlparms.h"
21
22#include <vector>
23
24using namespace std;
25
26
27/*--------------------------------------------------------------------------*/
28
29/* F(x,u,y,dim_x,dim_u) */
30/* norm(x,dim_x)        */
31typedef struct {
32    locint     edf_index;
33    int        sub_tape_num;
34    int      (*double_F)(double*, double* ,double*, int, int);
35    int      (*adouble_F)(adouble*, adouble*, adouble*, int, int);
36    double   (*norm)(double*, int);
37    double   (*norm_deriv)(double*, int);
38    double     epsilon;
39    double     epsilon_deriv;
40    int        N_max;
41    int        N_max_deriv;
42}
43fpi_data;
44
45
46static vector<fpi_data*> fpi_stack;
47
48
49static int iteration ( int dim_xu, double *xu, int dim_x, double *x_fix ) {
50    int i, k;
51    double err;
52    fpi_data *current = fpi_stack.back();
53    for (i=0; i<dim_x; i++) x_fix[i] = xu[i];
54    for (k=1; k<=current->N_max; k++) {
55        for (i=0; i<dim_x; i++) xu[i] = x_fix[i];
56        (*current->double_F)( xu, xu+dim_x, x_fix, dim_x, dim_xu-dim_x );
57        for (i=0; i<dim_x; i++) xu[i] = x_fix[i] - xu[i];
58        err = (*current->norm)(xu,dim_x);
59        if (err<current->epsilon) return k;
60    }
61    return -1;
62}
63
64
65static int fp_zos_forward ( int dim_xu, double *xu, int dim_x, double *x_fix ) {
66    int i, k;
67    double err;
68    ADOLC_OPENMP_THREAD_NUMBER;
69    ADOLC_OPENMP_GET_THREAD_NUMBER;
70    locint edf_index = ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index;
71    fpi_data *current=0;
72    vector<fpi_data*>::iterator fpi_stack_iterator;
73    for (fpi_stack_iterator=fpi_stack.begin();
74            fpi_stack_iterator!=fpi_stack.end();
75            ++fpi_stack_iterator) {
76        current = *fpi_stack_iterator;
77        if (edf_index==current->edf_index) break;
78    }
79    if (fpi_stack_iterator==fpi_stack.end()) {
80        fprintf(stderr,"ADOL-C Error! No edf found for fixpoint iteration.\n");
81        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
82    }
83    for (i=0; i<dim_x; i++) x_fix[i] = xu[i];
84    for (k=1; k<=current->N_max; k++) {
85        for (i=0; i<dim_x; i++) xu[i] = x_fix[i];
86        (*current->double_F)( xu, xu+dim_x, x_fix, dim_x, dim_xu-dim_x );
87        for (i=0; i<dim_x; i++) xu[i] = x_fix[i] - xu[i];
88        err = (*current->norm)(xu,dim_x);
89        if (err<current->epsilon) return k;
90    }
91    return -1;
92}
93
94static int fp_fos_forward ( int dim_xu, double *xu, double *xu_dot,
95                            int dim_x, double *x_fix, double *x_fix_dot) {
96    // Piggy back
97    int i, k;
98    double err, err_deriv;
99    ADOLC_OPENMP_THREAD_NUMBER;
100    ADOLC_OPENMP_GET_THREAD_NUMBER;
101    locint edf_index = ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index;
102    fpi_data *current=0;
103    vector<fpi_data*>::iterator fpi_stack_iterator;
104    for (fpi_stack_iterator=fpi_stack.begin();
105            fpi_stack_iterator!=fpi_stack.end();
106            ++fpi_stack_iterator) {
107        current = *fpi_stack_iterator;
108        if (edf_index==current->edf_index) break;
109    }
110    if (fpi_stack_iterator==fpi_stack.end()) {
111        fprintf(stderr,"ADOL-C Error! No edf found for fixpoint iteration.\n");
112        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
113    }
114    for (k=1; (k<current->N_max_deriv)|(k<current->N_max); k++) {
115        for (i=0; i<dim_x; i++) xu[i] = x_fix[i];
116        for (i=0; i<dim_x; i++) xu_dot[i] = x_fix_dot[i];
117        fos_forward ( current->sub_tape_num, dim_x, dim_xu, 0, xu, xu_dot, x_fix, x_fix_dot);
118        for (i=0; i<dim_x; i++)  xu[i] = x_fix[i] - xu[i];
119        err = (*current->norm)(xu,dim_x);
120        for (i=0; i<dim_x; i++) xu_dot[i] = x_fix_dot[i] -  xu_dot[i];
121        err_deriv = (*current->norm_deriv)(xu_dot,dim_x);
122        if ((err<current->epsilon)&(err_deriv<current->epsilon_deriv)) {
123            return k;
124        }
125    }
126    return -1;
127}
128
129static int fp_fos_reverse ( int dim_x, double *x_fix_bar, int dim_xu, double *xu_bar, double* /*unused*/, double* /*unused*/) {
130    // (d x_fix) / (d x_0) = 0 (!)
131    int i, k;
132    double err;
133    ADOLC_OPENMP_THREAD_NUMBER;
134    ADOLC_OPENMP_GET_THREAD_NUMBER;
135    locint edf_index = ADOLC_CURRENT_TAPE_INFOS.ext_diff_fct_index;
136    fpi_data *current=0;
137    vector<fpi_data*>::iterator fpi_stack_iterator;
138    for (fpi_stack_iterator=fpi_stack.begin();
139            fpi_stack_iterator!=fpi_stack.end();
140            ++fpi_stack_iterator) {
141        current = *fpi_stack_iterator;
142        if (edf_index==current->edf_index) break;
143    }
144    if (fpi_stack_iterator==fpi_stack.end()) {
145        fprintf(stderr,"ADOL-C Error! No edf found for fixpoint iteration.\n");
146        adolc_exit(-1,"",__func__,__FILE__,__LINE__);
147    }
148    double *U = new double[dim_xu];
149    double *xi = new double[dim_x];
150
151    for (k=1; k<current->N_max_deriv; k++) {
152        for (i=0; i<dim_x; i++) xi[i] = U[i];
153        fos_reverse ( current->sub_tape_num, dim_x, dim_xu, xi, U );
154        for (i=0; i<dim_x; i++) U[i] += x_fix_bar[i];
155        for (i=0; i<dim_x; i++) xi[i] = U[i] - xi[i];
156        err = (*current->norm_deriv)(xi,dim_x);
157        printf(" fp_fos_reverse: k = %d  err = %e\n",k,err);
158        if (err<current->epsilon_deriv) {
159            for (i=0; i<dim_xu-dim_x; i++) {
160                xu_bar[dim_x+i] += U[dim_x+i];
161            }
162
163            delete[] xi;
164            delete[] U;
165            return k;
166        }
167    }
168    for (i=0; i<dim_xu-dim_x; i++) xu_bar[dim_x+i] += U[dim_x+i];
169    delete[] xi;
170    delete[] U;
171    return -1;
172}
173
174
175int fp_iteration ( int        sub_tape_num,
176                   int      (*double_F)(double*, double* ,double*, int, int),
177                   int      (*adouble_F)(adouble*, adouble*, adouble*, int, int),
178                   double   (*norm)(double*, int),
179                   double   (*norm_deriv)(double*, int),
180                   double     epsilon,
181                   double     epsilon_deriv,
182                   int        N_max,
183                   int        N_max_deriv,
184                   adouble   *x_0,
185                   adouble   *u,
186                   adouble   *x_fix,
187                   int        dim_x,
188                   int        dim_u ) {
189    int i, k;
190    double dummy;
191    // add new fp information
192    fpi_data *data = new fpi_data;
193    data->sub_tape_num = sub_tape_num;
194    data->double_F     = double_F;
195    data->adouble_F    = adouble_F;
196    data->norm         = norm;
197    data->norm_deriv   = norm_deriv;
198    data->epsilon      = epsilon;
199    data->epsilon_deriv = epsilon_deriv;
200    data->N_max        = N_max;
201    data->N_max_deriv  = N_max_deriv;
202    fpi_stack.push_back(data);
203
204    // declare extern differentiated function and data
205    ext_diff_fct *edf_iteration = reg_ext_fct(&iteration);
206    data->edf_index = edf_iteration->index;
207    edf_iteration->zos_forward = &fp_zos_forward;
208    edf_iteration->fos_forward = &fp_fos_forward;
209    edf_iteration->fos_reverse = &fp_fos_reverse;
210
211    // put x and u together
212    adouble   *xu      = new adouble[dim_x+dim_u];
213    for (i=0; i<dim_x; i++) xu[i] = x_0[i];
214    for (i=0; i<dim_u; i++) xu[dim_x+i] = u[i];
215
216    k = call_ext_fct ( edf_iteration, dim_x+dim_u, xu,
217                       dim_x, x_fix );
218
219    // tape near solution
220    trace_on(sub_tape_num,1);
221    for (i=0; i<dim_x; i++) xu[i] <<= x_fix[i].getValue();
222    for (i=0; i<dim_u; i++) xu[dim_x+i] <<= u[i].getValue();
223    adouble_F(xu, xu+dim_x, x_fix, dim_x, dim_u);
224    for (i=0; i<dim_x; i++) x_fix[i] >>= dummy;
225    trace_off();
226
227    delete[] xu;
228    return k;
229}
Note: See TracBrowser for help on using the repository browser.