| 28 | // in a way that will record pivot operations on the AD<double> tape |
| 29 | typedef CPPAD_TESTVECTOR(CppAD::AD<double>) Vector; |
| 30 | Vector Solve(const Vector &a , const Vector &b) |
| 31 | { using namespace CppAD; |
| 32 | assert(a.size() == 4 && b.size() == 2); |
| 33 | |
| 34 | // copy the vector b into the VecAD object B |
| 35 | VecAD<double> B(2); |
| 36 | AD<double> u; |
| 37 | for(u = 0; u < 2; u += 1.) |
| 38 | B[u] = b[ size_t( Integer(u) ) ]; |
| 39 | |
| 40 | // copy the matrix a into the VecAD object A |
| 41 | VecAD<double> A(4); |
| 42 | for(u = 0; u < 4; u += 1.) |
| 43 | A[u] = a [ size_t( Integer(u) ) ]; |
| 44 | |
| 45 | // tape AD operation sequence that determines the row of A |
| 46 | // with maximum absolute element in column zero |
| 47 | AD<double> zero(0), one(1); |
| 48 | AD<double> rmax = CondExpGt(fabs(a[0]), fabs(a[2]), zero, one); |
| 49 | |
| 50 | // divide row rmax by A(rmax, 0) |
| 51 | A[rmax * 2 + 1] = A[rmax * 2 + 1] / A[rmax * 2 + 0]; |
| 52 | B[rmax] = B[rmax] / A[rmax * 2 + 0]; |
| 53 | A[rmax * 2 + 0] = one; |
| 54 | |
| 55 | // subtract A(other,0) times row A(rmax, *) from row A(other,*) |
| 56 | AD<double> other = one - rmax; |
| 57 | A[other * 2 + 1] = A[other * 2 + 1] |
| 58 | - A[other * 2 + 0] * A[rmax * 2 + 1]; |
| 59 | B[other] = B[other] |
| 60 | - A[other * 2 + 0] * B[rmax]; |
| 61 | A[other * 2 + 0] = zero; |
| 62 | |
| 63 | // back substitute to compute the solution vector x. |
| 64 | // Note that the columns of A correspond to rows of x. |
| 65 | // Also note that A[rmax * 2 + 0] is equal to one. |
| 66 | CPPAD_TESTVECTOR(AD<double>) x(2); |
| 67 | x[1] = B[other] / A[other * 2 + 1]; |
| 68 | x[0] = B[rmax] - A[rmax * 2 + 1] * x[1]; |
| 69 | |
| 70 | return x; |
| 71 | } |
| 72 | } |
| 73 | |
| 74 | bool vec_ad(void) |