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
25 BEGIN_C_DECLS
26
27 /****************************************************************************/
28 /* DRIVERS FOR PS FUNCTIONS */
29
30 /*--------------------------------------------------------------------------*/
31 /* abs-normal form */
32
abs_normal(short tag,int m,int n,int swchk,double * x,double * y,double * z,double * cz,double * cy,double ** Y,double ** J,double ** Z,double ** L)33 int 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 double *y, /* function value */
39 double *z, /* switching variables */
40 double *cz, /* first constant */
41 double *cy, /* second constant */
42 double **Y, /* m times n */
43 double **J, /* m times s */
44 double **Z, /* s times n */
45 double **L) /* s times s (lowtri) */
46 {
47
48 int i,j,s;
49 double *res, tmp;
50 s=get_num_switches(tag);
51
52 /* This check is required because the user is probably allocating his
53 * arrays sigma, cz, Z, L, Y, J according to swchk */
54 if (s != swchk) {
55 fprintf(DIAG_OUT, "ADOL-C error: Number of switches passed %d does not "
56 "match the one recorded on tape %d (%zu)\n", swchk, tag, s);
57 adolc_exit(-1,"",__func__,__FILE__,__LINE__);
58 }
59
60 res=(double*)myalloc1(n+s);
61
62 zos_pl_forward(tag,m,n,1,x,y,z);
63
64 for(i=0;i<m+s;i++){
65 int l = i - s;
66 fos_pl_reverse(tag,m,n,s,i,res);
67 if ( l < 0 ) {
68 cz[i]=z[i];
69 for(j=0;j<n;j++){
70 Z[i][j]=res[j];
71 }
72 for(j=0;j<s;j++) { /* L[i][i] .. L[i][s] are theoretically zero,
73 * we probably don't need to copy them */
74 L[i][j]=res[j+n];
75 if (j < i)
76 {
77 cz[i] = cz[i]-L[i][j]*fabs(z[j]);
78 }
79 }
80 } else {
81 cy[l]=y[l];
82 for(j=0;j<n;j++){
83 Y[l][j]=res[j];
84 }
85 for(j=0;j<s;j++){
86 J[l][j]=res[j+n];
87 cy[l] = cy[l]-J[l][j]*fabs(z[j]);
88 }
89 }
90 }
91
92 myfree1(res);
93 return 0;
94 }
95
96
97 /*--------------------------------------------------------------------------*/
98 /* directional_active_gradient */
99 /* */
directional_active_gradient(short tag,int n,double * x,double * d,double * g,short * sigma_g)100 int directional_active_gradient(short tag, /* trace identifier */
101 int n, /* number of independents */
102 double* x, /* value of independents */
103 double* d, /* direction */
104 double* g, /* directional active gradient */
105 short *sigma_g /* sigma of g */
106 )
107 {
108 int i, j, p, k, s, max_dk, done, sum, keep;
109 double max_entry, y, by;
110 double *z;
111 double **E, **invE, **grad, **gradu;
112
113 keep = 1;
114 by = 1;
115
116 s=get_num_switches(tag);
117
118 z = myalloc1(s);
119
120 grad = (double**)myalloc2(1,n);
121 gradu = (double**)myalloc2(s,n);
122 E = (double**)myalloc2(n,n);
123
124 #if !defined(ADOLC_USE_CALLOC)
125 memset(&(E[0][0]), 0, n*n*sizeof(double));
126 #endif
127
128 max_dk=0;
129 max_entry = -1;
130 for(i=0;i<n;i++){
131 E[i][0] = d[i];
132 if(max_entry<fabs(d[i])){
133 max_dk=i;
134 max_entry = fabs(d[i]);
135 }
136 }
137
138 k = 1; done = 0;
139 j = 0;
140
141 while((k<6) && (done == 0))
142 {
143 fov_pl_forward(tag,1,n,k,x,E,&y,grad,z,gradu,sigma_g);
144
145 sum = 0;
146 for(i=0;i<s;i++)
147 {
148 sum += fabs(sigma_g[i]);
149 }
150
151 if (sum == s)
152 {
153
154 zos_pl_forward(tag,1,n,keep,x,&y,z);
155 fos_pl_sig_reverse(tag,1,n,s,sigma_g, &by ,g);
156 done = 1;
157 }
158 else
159 {
160 if(j==max_dk)
161 j++;
162 E[j][k]=1;
163 j++;
164 k++;
165 }
166 }
167
168 myfree1(z); myfree2(E); myfree2(grad); myfree2(gradu);
169
170 if (done == 0)
171 {
172 fprintf(DIAG_OUT," NOT ENOUGH DIRECTIONS !!!!\n");
173 adolc_exit(-1,"",__func__,__FILE__,__LINE__);
174 }
175
176 return 0;
177 }
178