/*---------------------------------------------------------------------------- ADOL-C -- Automatic Differentiation by Overloading in C++ File: hessmat.cpp Revision: \$Id: hessmat.cpp 171 2010-10-04 13:57:19Z kulshres \$ Contents: example for testing the routines: hov_wk_forward ( = Higher Order Vector forward With Keep ) hos_ov_reverse ( = Higher Order Scalar reverse over vectors) Copyright (c) Andrea Walther, Andreas Kowarz, Olaf Vogel This file is part of ADOL-C. This software is provided as open source. Any use, reproduction, or distribution of the software constitutes recipient's acceptance of the terms of the accompanying license file. ---------------------------------------------------------------------------*/ /****************************************************************************/ /* INCLUDES */ #include #include #include using namespace std; /****************************************************************************/ /* MAIN */ int main() { int i,j,l,m,n,d,q,bd, keep; /*--------------------------------------------------------------------------*/ /* inputs */ cout << "vector x Hessian x matrix for the function \n\n"; cout << " y[0] = cos(x[0])* ...*cos(x[n]) \n"; cout << " y[1] = x[0]^n \n"; cout << " y[2] = condassign(y[i],y[0]>y[1],y[1],y[0]) \n"; cout << " y[3] = sin(x[0])+ ...+sin(x[n]) \n"; cout << " y[4] = exp(x[0])- ...-exp(x[n]) \n"; cout << " y[5] = pow(y[1],3) \n"; cout << " y[6] += y[5]*y[4] \n"; cout << " y[7] -= y[6]*y[5] \n"; cout << " y[j] = 1/x[0]/ .../x[n], j > 3 \n\n"; cout << " Number of independents = ?\n "; cin >> n; cout << " Number of dependents = ?\n "; cin >> m; cout << " Degree d (for forward) = ?\n"; cin >> d; cout << " keep (degree of corresponding reverse = keep-1) = ?\n"; cout << " keep <= d+1 must be valid \n"; cin >> keep; cout << " Number of directions = ?\n "; cin >> q; /*--------------------------------------------------------------------------*/ /* allocations and inits */ double* xp = new double[n]; /* passive indeps */ double* yp = new double[m]; /* passive depends */ /* vector x Hessian x matrix = Upp x H x XPPP */ double* Up = myalloc(m); /* vector on left-hand side */ double** Upp = myalloc(m,d+1); /* vector on left-hand side */ double*** Xppp = myalloc(n,q,d); /* matrix on right-hand side */ double*** Zppp = myalloc(q,n,d+1); /* result of Up x H x XPPP */ double*** Yppp = myalloc(m,q,d); /* results of needed hos_wk_forward */ /* check results with usual lagra-Hess-vec */ double** Xpp = myalloc(n,d); double** V = myalloc(n,q); double** W = myalloc(q,n); double** H = myalloc(n,n); double** Ypp = myalloc(m,d); double** Zpp = myalloc(n,d+1); /* inits */ for (l=0; ly[1],y[1],y[0]); break; case 3 : y[i] -= sin(x[j]); break; case 4 : y[i] -= exp(x[j]); break; case 5 : y[5] = pow(y[1],3); case 6 : y[6] += y[5]*y[4]; case 7 : y[7] -= y[6]*y[5]; default : y[i] /= x[j]; } } for (i=0; i>= yp[i] ; trace_off(); /*--------------------------------------------------------------------------*/ /* work on the tape */ /* compute results of lagra_hess_vec */ /* the following is equal to calls inside of lagra_hess_vec(..) */ /* direct calls to the basic routines hos_forward and hos_reverse */ /* seem to be faster than call of lagra_hess_vec(..) */ /* at least in some of our test cases */ hos_forward(1,m,n,d,keep,xp,Xpp,yp,Ypp); hos_reverse(1,m,n,keep-1,Up,Zpp); printf("\n Results of hos_reverse:\n\n"); for (i=0; i<=d; i++) { printf(" d = %d \n",i); for (j=0;j