# Changeset 2397

Ignore:
Timestamp:
May 23, 2012 1:10:12 AM (22 months ago)
Message:

Change row and colum indices from (r,c) to (row, col).

File:
1 edited

### Legend:

Unmodified
 r2394 $icode%hes% = %f%.SparseHessian(%x%, %w%) %hes% = %f%.SparseHessian(%x%, %w%, %p%) %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %r%, %c%, %hes%, %work%) %n_sweep% = %f%.SparseHessian(%x%, %w%, %p%, %row%, %col%, %hes%, %work%) %$$This routine takes advantage of the sparsity of the Hessian in order to reduce the amount of computation necessary. If icode r$$ and$icode c$$are present, it also takes If icode row$$ and $icode col$$are present, it also takes advantage of the reduced set of elements of the Hessian that need to be computed. for the internal calculations is unspecified. head r, c$$ The arguments$icode r$$and icode c$$ are optional and have prototype $head row, col$$The arguments icode row$$ and$icode col$$are optional and have prototype codei% const %VectorSize%& %r% const %VectorSize%& %c% const %VectorSize%& %row% const %VectorSize%& %col% %$$ (see $cref/VectorSize/sparse_hessian/VectorSize/$$below). returned and in what order. We use latex K$$ to denote the value$icode%hes%.size()%$$which must also equal the size of icode r$$ and $icode c$$. which must also equal the size of icode row$$ and$icode col$$. Furthermore, for latex k = 0 , \ldots , K-1$$, it must hold that $latex r[k] < n$$and latex c[k] < n$$.$latex row[k] < n$$and latex col[k] < n$$. In addition, all of the $latex (r[k], c[k])$$pairs must correspond to a true value all of the latex (row[k], col[k])$$ pairs must correspond to a true value in the sparsity pattern$icode p$$. %VectorBase% %hes% %$$ In the case where $icode r$$and icode c$$ are not present, In the case where$icode row$$and icode col$$ are not present, the size of $icode hes$$is latex n * n$$ and its size is$latex n * n$$.$$ In the case where the arguments $icode r$$and icode c$$ are present, In the case where the arguments$icode row$$and icode col$$ are present, we use $latex K$$to denote the size of icode hes$$. The input value of its elements does not matter. \; , \; \; {\rm where} \; j = r[k] j = row[k] \; {\rm and } \; \ell = c[k] \ell = col[k] \] $$During its the first use, information is stored in icode work$$. This is used to reduce the work done by future calls to$code SparseHessian$$with the same icode f$$, $icode p$$, icode r$$, and$icode c$$. with the same icode f$$, $icode p$$, icode row$$, and$icode col$$. If a future call is make where any of these values have changed, you must first call icode%work%.clear()%$$ \tparam Base See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorBase See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorSet \param x See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param w See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param sparsity \param hes See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param work See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \return See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). */ template size_t i, j, k, ell; CppAD::vector& r(work.r_sort); CppAD::vector& c(work.c_sort); CppAD::vector& row(work.r_sort); CppAD::vector& col(work.c_sort); CppAD::vector& user_k(work.k_sort); CppAD::vector& color(work.color); // number of components of Hessian that are required size_t K = hes.size(); CPPAD_ASSERT_UNKNOWN( r.size() == K+1 ); CPPAD_ASSERT_UNKNOWN( c.size() == K ); CPPAD_ASSERT_UNKNOWN( r[K] == n ); CPPAD_ASSERT_UNKNOWN( row.size() == K+1 ); CPPAD_ASSERT_UNKNOWN( col.size() == K ); CPPAD_ASSERT_UNKNOWN( row[K] == n ); // Point at which we are evaluating the Hessian k = 0; while( k < K ) {       CPPAD_ASSERT_UNKNOWN( r[k] < n && c[k] < n ); CPPAD_ASSERT_UNKNOWN( k == 0 || r[k-1] <= r[k] ); {       CPPAD_ASSERT_UNKNOWN( row[k] < n && col[k] < n ); CPPAD_ASSERT_UNKNOWN( k == 0 || row[k-1] <= row[k] ); CPPAD_ASSERT_KNOWN( sparsity.is_element(r[k], c[k]) , "SparseHessian: an (r,c) pair is not in sparsity pattern." sparsity.is_element(row[k], col[k]) , "SparseHessian: an (row, col) pair is not in sparsity pattern." ); r_used.add_element(c[k], r[k]); c_used.add_element(r[k], c[k]); r_used.add_element(col[k], row[k]); c_used.add_element(row[k], col[k]); k++; } k = 0; for(i = 0; i < n; i++) if( color[i] == ell ) {       // find first k such that r[k] has color ell {       // find first k such that row[k] has color ell if( ! any ) {       while( r[k] < i ) {       while( row[k] < i ) k++; any = r[k] == i; any = row[k] == i; } } for(i = 0; i < n; i++) if( color[i] == ell ) {       // find first index in c for this column while( r[k] < i ) while( row[k] < i ) k++; // extract the results for this row while( r[k] == i ) {       hes[ user_k[k] ] = ddw[ c[k] * 2 + 1 ]; while( row[k] == i ) {       hes[ user_k[k] ] = ddw[ col[k] * 2 + 1 ]; k++; } \tparam Base See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorBase See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorSet \param x See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param w See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param p \param hes See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param work See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \return See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). */ template \tparam Base See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorBase See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \tparam VectorSet \param x See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param w See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param p \param hes See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \param work See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). \return See \c SparseHessian(x, w, p, r, c, hes, work). See \c SparseHessian(x, w, p, row, col, hes, work). */ template ); sparse_hessian_work work; CppAD::vector& r(work.r_sort); CppAD::vector& c(work.c_sort); CppAD::vector& row(work.r_sort); CppAD::vector& col(work.c_sort); CppAD::vector& user_k(work.k_sort); size_t K = k; VectorBase H(K); r.resize(K+1); c.resize(K); row.resize(K+1); col.resize(K); user_k.resize(K); k = 0; {       for(j = 0; j < n; j++) if( p[i * n + j] ) {       r[k]      = i; c[k]      = j; {       row[k]      = i; col[k]      = j; user_k[k] = k; k++; } } r[K] = n; row[K] = n; SparseHessianCase(set_type, x, w, p, H, work); } for(k = 0; k < K; k++) hes[r[k] * n + c[k]] = H[k]; hes[row[k] * n + col[k]] = H[k]; } /*! ); sparse_hessian_work work; CppAD::vector& r(work.r_sort); CppAD::vector& c(work.c_sort); CppAD::vector& row(work.r_sort); CppAD::vector& col(work.c_sort); CppAD::vector& user_k(work.k_sort); size_t K = k; VectorBase H(K); r.resize(K+1); c.resize(K); row.resize(K+1); col.resize(K); user_k.resize(K); k = 0; {       itr = p[i].begin(); while( itr != p[i].end() ) {       r[k]      = i; c[k]      = *itr; {       row[k]      = i; col[k]      = *itr; user_k[k] = k; itr++; } } r[K] = n; row[K] = n; SparseHessianCase(set_type, x, w, p, H, work); } for(k = 0; k < K; k++) hes[r[k] * n + c[k]] = H[k]; hes[row[k] * n + col[k]] = H[k]; } // =========================================================================== The C++ source code corresponding to this operation is \verbatim SparceHessian(x, w, p, r, c, hes, work) SparceHessian(x, w, p, row, col, hes, work) \endverbatim The return value hes[k] is the second partial of \f$w^{\rm T} F(x)\f$ with respect to the r[k] and c[k] component of \f$x\f$. row[k] and col[k] component of \f$x\f$. \param work contains information that depends on the function object, sparsity pattern, \c r, and \c c vector. \c row, and \c col vector. If these values are the same, \c work does not need to be recomputed. To be more specific, \c r_sort is sorted copy of \c r , \c c_sort is sorted copy of \c c , \c r_sort is sorted copy of \c row , \c c_sort is sorted copy of \c col , k_sort[k] is the original index corresponding to the values r_sort[k] and c_sort[k]. The order for the sort is by columns. Let \c n the domain dimension, and \c K the size of \c r , \c c , and \c hes. and \c K the size of \c row , \c col , and \c hes. There is one extra entry in the sorted row array and it has value r_sort[K]=n. const VectorBase&     w    , const VectorSet&      p    , const VectorSize&     r    , const VectorSize&     c    , const VectorSize&     row  , const VectorSize&     col  , VectorBase&           hes  , sparse_hessian_work&  work ) size_t k, K = hes.size(); if( work.r_sort.size() == 0 ) {       // create version of (r,c,k) sorted by rows {       // create version of (row, col,k) sorted by rows size_t min_bytes = 3 * K * sizeof(size_t); size_t cap_bytes; for(k = 0; k < K; k++) {       // must put column first so it is used for sorting rck[3 * k + 0] = r[k]; rck[3 * k + 1] = c[k]; rck[3 * k + 0] = row[k]; rck[3 * k + 1] = col[k]; rck[3 * k + 2] = k; } ); CPPAD_ASSERT_KNOWN( r.size() == K && c.size() == K , row.size() == K && col.size() == K , "SparseHessian: either r or c does not have the same size as ehs." ); for(k = 0; k < K; k++) {       CPPAD_ASSERT_KNOWN( r[k] < n, row[k] < n, "SparseHessian: invalid value in r." ); CPPAD_ASSERT_KNOWN( c[k] < n, col[k] < n, "SparseHessian: invalid value in c." ); ); CPPAD_ASSERT_KNOWN( work.r_sort[k] == r[ work.k_sort[k] ] , work.r_sort[k] == row[ work.k_sort[k] ] , "SparseHessian: invalid value in work." ); CPPAD_ASSERT_KNOWN( work.c_sort[k] == c[ work.k_sort[k] ], work.c_sort[k] == col[ work.k_sort[k] ], "SparseHessian: invalid value in work." );