Changeset 3941 for trunk/example/sparse/colpack_jac.cpp
 Timestamp:
 Jun 2, 2017 1:36:10 AM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/example/sparse/colpack_jac.cpp
r3875 r3941 17 17 $$ 18 18 19 $section Using ColPack: Example and Test$$ 20 $mindex colpack jacobian sparse$$ 19 $section ColPack: Sparse Jacobian Example and Test$$ 21 20 22 21 … … 34 33 using CppAD::AD; 35 34 using CppAD::NearEqual; 36 typedef CPPAD_TESTVECTOR(AD<double>) a_vector;37 typedef CPPAD_TESTVECTOR(double) d_vector;38 typedef CppAD::vector<size_t> i_vector;39 size_t i, j, k, ell;40 double eps = 10. * CppAD::numeric_limits<double>::epsilon();35 typedef CPPAD_TESTVECTOR(AD<double>) a_vector; 36 typedef CPPAD_TESTVECTOR(double) d_vector; 37 typedef CppAD::vector<size_t> i_vector; 38 typedef CppAD::sparse_rc<i_vector> sparsity; 39 typedef CppAD::sparse_rcv<i_vector, d_vector> sparse_matrix; 41 40 42 41 // domain space vector 43 42 size_t n = 4; 44 43 a_vector a_x(n); 45 for( j = 0; j < n; j++)44 for(size_t j = 0; j < n; j++) 46 45 a_x[j] = AD<double> (0); 47 46 … … 60 59 // new value for the independent variable vector 61 60 d_vector x(n); 62 for( j = 0; j < n; j++)61 for(size_t j = 0; j < n; j++) 63 62 x[j] = double(j); 64 63 … … 68 67 [ 1 1 1 x_3] 69 68 */ 70 d_vector check(m * n); 71 check[0] = 1.; check[1] = 1.; check[2] = 0.; check[3] = 0.; 72 check[4] = 0.; check[5] = 0.; check[6] = 1.; check[7] = 1.; 73 check[8] = 1.; check[9] = 1.; check[10] = 1.; check[11] = x[3]; 74 75 // Normally one would use f.ForSparseJac or f.RevSparseJac to compute 76 // sparsity pattern, but for this example we extract it from check. 77 std::vector< std::set<size_t> > p(m); 69 // Normally one would use CppAD to compute sparsity pattern, but for this 70 // example we set it directly 71 size_t nr = m; 72 size_t nc = n; 73 size_t nnz = 8; 74 sparsity pattern(nr, nc, nnz); 75 d_vector check(nnz); 76 for(size_t k = 0; k < nnz; k++) 77 { size_t r, c; 78 if( k < 2 ) 79 { r = 0; 80 c = k; 81 } 82 else if( k < 4 ) 83 { r = 1; 84 c = k; 85 } 86 else 87 { r = 2; 88 c = k  4; 89 } 90 pattern.set(k, r, c); 91 if( k == nnz  1 ) 92 check[k] = x[3]; 93 else 94 check[k] = 1.0; 95 } 78 96 79 97 // using row and column indices to compute nonzero in rows 1 and 2 80 i_vector row, col; 81 for(i = 0; i < m; i++) 82 { for(j = 0; j < n; j++) 83 { ell = i * n + j; 84 if( check[ell] != 0. ) 85 { row.push_back(i); 86 col.push_back(j); 87 p[i].insert(j); 88 } 98 sparse_matrix subset( pattern ); 99 100 // check results for both CppAD and Colpack 101 for(size_t i_method = 0; i_method < 4; i_method++) 102 { // coloring method 103 std::string coloring; 104 if( i_method % 2 == 0 ) 105 coloring = "cppad"; 106 else 107 coloring = "colpack"; 108 // 109 CppAD::sparse_jac_work work; 110 size_t group_max = 1; 111 if( i_method / 2 == 0 ) 112 { size_t n_sweep = f.sparse_jac_for( 113 group_max, x, subset, pattern, coloring, work 114 ); 115 ok &= n_sweep == 4; 89 116 } 117 else 118 { size_t n_sweep = f.sparse_jac_rev( 119 x, subset, pattern, coloring, work 120 ); 121 ok &= n_sweep == 2; 122 } 123 const d_vector& hes( subset.val() ); 124 for(size_t k = 0; k < nnz; k++) 125 ok &= check[k] == hes[k]; 90 126 } 91 size_t K = row.size();92 d_vector jac(K);93 94 // empty work structure95 CppAD::sparse_jacobian_work work;96 ok &= work.color_method == "cppad";97 98 // choose to use ColPack99 work.color_method = "colpack";100 101 // forward mode102 size_t n_sweep = f.SparseJacobianForward(x, p, row, col, jac, work);103 for(k = 0; k < K; k++)104 { ell = row[k] * n + col[k];105 ok &= NearEqual(check[ell], jac[k], eps, eps);106 }107 ok &= n_sweep == 4;108 109 // reverse mode110 work.clear();111 work.color_method = "colpack";112 n_sweep = f.SparseJacobianReverse(x, p, row, col, jac, work);113 for(k = 0; k < K; k++)114 { ell = row[k] * n + col[k];115 ok &= NearEqual(check[ell], jac[k], eps, eps);116 }117 ok &= n_sweep == 2;118 119 127 return ok; 120 128 }
Note: See TracChangeset
for help on using the changeset viewer.