| 72 | } |
| 73 | |
| 74 | bool vec_ad(void) |
| 75 | { bool ok = true; |
| 76 | |
| 77 | using CppAD::AD; |
| 78 | using CppAD::NearEqual; |
| 79 | double eps99 = 99.0 * std::numeric_limits<double>::epsilon(); |
| 80 | |
| 81 | // domain space vector |
| 82 | size_t n = 4; |
| 83 | CPPAD_TESTVECTOR(double) x(n); |
| 84 | CPPAD_TESTVECTOR(AD<double>) X(n); |
| 85 | // 2 * identity matrix (rmax in Solve will be 0) |
| 86 | X[0] = x[0] = 2.; X[1] = x[1] = 0.; |
| 87 | X[2] = x[2] = 0.; X[3] = x[3] = 2.; |
| 88 | |
| 89 | // declare independent variables and start tape recording |
| 90 | CppAD::Independent(X); |
| 91 | |
| 92 | // define the vector b |
| 93 | CPPAD_TESTVECTOR(double) b(2); |
| 94 | CPPAD_TESTVECTOR(AD<double>) B(2); |
| 95 | B[0] = b[0] = 0.; |
| 96 | B[1] = b[1] = 1.; |
| 97 | |
| 98 | // range space vector solves X * Y = b |
| 99 | size_t m = 2; |
| 100 | CPPAD_TESTVECTOR(AD<double>) Y(m); |
| 101 | Y = Solve(X, B); |
| 102 | |
| 103 | // create f: X -> Y and stop tape recording |
| 104 | CppAD::ADFun<double> f(X, Y); |
| 105 | |
| 106 | // By Cramer's rule: |
| 107 | // y[0] = [ b[0] * x[3] - x[1] * b[1] ] / [ x[0] * x[3] - x[1] * x[2] ] |
| 108 | // y[1] = [ x[0] * b[1] - b[0] * x[2] ] / [ x[0] * x[3] - x[1] * x[2] ] |
| 109 | |
| 110 | double den = x[0] * x[3] - x[1] * x[2]; |
| 111 | double dsq = den * den; |
| 112 | double num0 = b[0] * x[3] - x[1] * b[1]; |
| 113 | double num1 = x[0] * b[1] - b[0] * x[2]; |
| 114 | |
| 115 | // check value |
| 116 | ok &= NearEqual(Y[0] , num0 / den, eps99, eps99); |
| 117 | ok &= NearEqual(Y[1] , num1 / den, eps99, eps99); |
| 118 | |
| 119 | // forward computation of partials w.r.t. x[0] |
| 120 | CPPAD_TESTVECTOR(double) dx(n); |
| 121 | CPPAD_TESTVECTOR(double) dy(m); |
| 122 | dx[0] = 1.; dx[1] = 0.; |
| 123 | dx[2] = 0.; dx[3] = 0.; |
| 124 | dy = f.Forward(1, dx); |
| 125 | ok &= NearEqual(dy[0], 0. - num0 * x[3] / dsq, eps99, eps99); |
| 126 | ok &= NearEqual(dy[1], b[1] / den - num1 * x[3] / dsq, eps99, eps99); |
| 127 | |
| 128 | // compute the solution for a new x matrix such that pivioting |
| 129 | // on the original rmax row would divide by zero |
| 130 | CPPAD_TESTVECTOR(double) y(m); |
| 131 | x[0] = 0.; x[1] = 2.; |