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

Last change on this file since 71 was 42, checked in by awalther, 10 years ago

set svn keywords property

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