source: trunk/ADOL-C/src/drivers/psdrivers.c @ 706

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

fix number of arguments

File size: 4.9 KB
Line 
1/*----------------------------------------------------------------------------
2 ADOL-C -- Automatic Differentiation by Overloading in C++
3 File:     drivers/psdrivers.c
4 Revision: $Id$
5 Contents: Easy to use drivers for piecewise smooth functions
6           (with C and C++ callable interfaces including Fortran
7            callable versions).
8
9 Copyright (c) Andrea Walther, Sabrina Fiege
10
11 This file is part of ADOL-C. This software is provided as open source.
12 Any use, reproduct ion, or distribution of the software constitutes
13 recipient's acceptance of the terms of the accompanying license file.
14
15----------------------------------------------------------------------------*/
16#include <adolc/drivers/psdrivers.h>
17#include <adolc/interfaces.h>
18#include <adolc/adalloc.h>
19#include "taping_p.h"
20#include "dvlparms.h"
21
22#include <math.h>
23
24
25BEGIN_C_DECLS
26
27/****************************************************************************/
28/*                                                 DRIVERS FOR PS FUNCTIONS */
29
30/*--------------------------------------------------------------------------*/
31/*                                                          abs-normal form */
32
33int abs_normal(short tag,      /* tape identifier */
34               int m,          /* number od dependents   */
35               int n,          /* number of independents */
36               int swchk,      /* number of switches (check) */
37               double *x,      /* base point */
38               short *sigma,   /* sigma of x */
39               double *y,      /* function value */
40               double *z,      /* switching variables */
41               double *cz,     /* first constant */
42               double *cy,     /* second constant */
43               double **J,
44               double **Y,
45               double **Z,
46               double **L)
47{
48
49  int i,j,s;
50  double *res, tmp;
51  s=get_num_switches(tag);
52
53  /* This check is required because the user is probably allocating his
54   * arrays sigma, cz, Z, L, Y, J according to swchk */
55  if (s != swchk) {
56      fprintf(DIAG_OUT, "ADOL-C error: Number of switches passed %d does not "
57              "match the one recorded on tape %d (%zu)\n", swchk, tag, s);
58      adolc_exit(-1,"",__func__,__FILE__,__LINE__);
59  }
60
61  res=(double*)myalloc1(n+s);
62
63  zos_pl_forward(tag,m,n,1,x,y,z);
64
65  for(i=0;i<m+s;i++){
66    int l = i - s;
67    fos_pl_reverse(tag,m,n,s,i,res);
68    if ( l < 0 ) {
69        cz[i]=z[i];
70        for(j=0;j<n;j++){
71            Z[i][j]=res[j];
72        }
73        for(j=0;j<s;j++) { /* L[i][i] .. L[i][s] are theoretically zero,
74                            *  we probably don't need to copy them */
75            L[i][j]=res[j+n];
76            if (j < i)
77              {
78                cz[i] = cz[i]-L[i][j]*fabs(z[j]);
79              }
80        }
81    } else {
82        cy[l]=y[l];
83        for(j=0;j<n;j++){
84            J[l][j]=res[j];
85        }
86        for(j=0;j<s;j++){
87            Y[l][j]=res[j+n];
88            cy[l] = cy[l]-Y[l][j]*fabs(z[j]);
89        }
90    }
91  }
92
93  myfree1(res);
94  return 0;
95}
96
97
98/*--------------------------------------------------------------------------*/
99/*                                              directional_active_gradient */
100/*                                                                          */
101int directional_active_gradient(short tag,      /* trace identifier */
102                                int n,          /* number of independents */
103                                double* x,      /* value of independents */
104                                double* d,      /* direction */
105                                double* g,      /* directional active gradient */
106                                short *sigma_g  /* sigma of g */
107                                )
108{
109  int i, j, p, k, s, max_dk, done, sum, keep;
110  double max_entry, y, by;
111  double *z;
112  double **E, **invE, **grad, **gradu;
113
114  keep = 1;
115  by = 1;
116
117  s=get_num_switches(tag);
118
119  z = myalloc1(s);
120
121  grad = (double**)myalloc2(1,n);
122  gradu = (double**)myalloc2(s,n);
123  E = (double**)myalloc2(n,n);
124
125#if !defined(ADOLC_USE_CALLOC)
126  memset(&(E[0][0]), 0, n*n*sizeof(double));
127#endif
128
129  max_dk=0;
130  max_entry = -1;
131  for(i=0;i<n;i++){
132    E[i][0] = d[i];
133    if(max_entry<fabs(d[i])){
134      max_dk=i;
135      max_entry = fabs(d[i]);
136    }
137  }
138
139  k = 1; done = 0;
140  j = 0;
141
142  while((k<6) && (done == 0))
143    {
144      fov_pl_forward(tag,1,n,k,x,E,&y,grad,z,gradu,sigma_g);
145
146      sum = 0;
147      for(i=0;i<s;i++)
148        {
149          sum += fabs(sigma_g[i]);
150        }
151
152       if (sum == s)
153        {
154
155          zos_pl_forward(tag,1,n,keep,x,&y,z);
156          fos_pl_sig_reverse(tag,1,n,s,sigma_g, &by ,g);
157          done = 1;
158        }
159      else
160        {
161          if(j==max_dk)
162            j++;
163          E[j][k]=1;
164          j++;
165          k++;
166        }
167    }
168
169  myfree1(z); myfree2(E); myfree2(grad); myfree2(gradu);
170
171  if (done == 0)
172    {
173      fprintf(DIAG_OUT," NOT ENOUGH DIRECTIONS !!!!\n");
174      adolc_exit(-1,"",__func__,__FILE__,__LINE__);
175    }
176
177  return 0;
178}
Note: See TracBrowser for help on using the repository browser.