Changeset 3673


Ignore:
Timestamp:
Apr 18, 2015 3:03:50 PM (5 years ago)
Author:
bradbell
Message:

merge to branch: trunk
from repository: https://github.com/coin-or/CppAD
start hash code: 65943182bbcf72a51ad9814fd7c2d933f7ed586c
end hash code: fe0bc2b6886b5fb80dc912dfb22e6c7f6e8c0545

commit fe0bc2b6886b5fb80dc912dfb22e6c7f6e8c0545
Author: Brad Bell <bradbell@…>
Date: Sat Apr 18 12:01:52 2015 -0700

Remove trialing white space.

commit 6e8df3db4dac5d1ae85beae04922a737210b82b3
Author: Brad Bell <bradbell@…>
Date: Sat Apr 18 12:01:13 2015 -0700

Document and test fact that sparsity pattern is only needed
when computing work structure.


boost_lu.sh: delete trailing white space.

commit 3ad59e465ad64fc3cda81f48ac4ccd20b5990190
Author: Brad Bell <bradbell@…>
Date: Tue Apr 14 08:04:41 2015 -0700

boost_lu.sh: Bug report by Peter Caspers.

Location:
trunk
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed
  • trunk/cppad/local/sparse_hessian.hpp

    r3651 r3673  
    112112$codei SparseHessian$$, it should be faster to calculate $icode p$$ once and
    113113pass this argument to $codei SparseHessian$$.
     114Furthermore, if you specify $icode work$$ in the calling sequence,
     115it is not necessary to keep the sparsity pattern; see the heading
     116$cref/p/sparse_hessian/work/p/$$ under the $icode work$$ description.
     117$pre
     118
     119$$
    114120In addition,
    115121if you specify $icode p$$, CppAD will use the same
     
    212218$cref/cmake command/cmake/CMake Command/$$ line.
    213219It also takes advantage of the fact that the Hessian matrix is symmetric.
     220
     221$subhead p$$
     222If $icode work$$ is present, and it is not the first call after
     223its construction or a clear,
     224the sparsity pattern $icode p$$ is not used.
     225This enables one to free the sparsity pattern
     226and still compute corresponding sparse Hessians.
    214227
    215228$head n_sweep$$
  • trunk/cppad/local/sparse_jacobian.hpp

    r3552 r3673  
    44
    55/* --------------------------------------------------------------------------
    6 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell
     6CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
    77
    88CppAD is distributed under multiple licenses. This distribution is under
    9 the terms of the 
     9the terms of the
    1010                    Eclipse Public License Version 1.0.
    1111
     
    3232        jac
    3333        Jacobian
     34        Jacobians
    3435        const
    3536        Taylor
     
    5253We use $latex F : \B{R}^n \rightarrow \B{R}^m$$ do denote the
    5354$cref/AD function/glossary/AD Function/$$
    54 corresponding to $icode f$$. 
    55 The syntax above sets $icode jac$$ to the Jacobian 
     55corresponding to $icode f$$.
     56The syntax above sets $icode jac$$ to the Jacobian
    5657$latex \[
    57         jac = F^{(1)} (x) 
     58        jac = F^{(1)} (x)
    5859\] $$
    5960This routine takes advantage of the sparsity of the Jacobian
     
    8081%$$
    8182(see $cref/VectorBase/sparse_jacobian/VectorBase/$$ below)
    82 and its size 
     83and its size
    8384must be equal to $icode n$$, the dimension of the
    8485$cref/domain/seq_property/Domain/$$ space for $icode f$$.
     
    9798its size is $latex m$$ and all its set elements are between
    9899zero and $latex n - 1$$.
    99 It specifies a 
    100 $cref/sparsity pattern/glossary/Sparsity Pattern/$$ 
     100It specifies a
     101$cref/sparsity pattern/glossary/Sparsity Pattern/$$
    101102for the Jacobian $latex F^{(1)} (x)$$.
    102103$pre
    103104
    104105$$
    105 If this sparsity pattern does not change between calls to 
    106 $codei SparseJacobian$$, it should be faster to calculate $icode p$$ once 
     106If this sparsity pattern does not change between calls to
     107$codei SparseJacobian$$, it should be faster to calculate $icode p$$ once
    107108(using $cref ForSparseJac$$ or $cref RevSparseJac$$)
    108109and then pass $icode p$$ to $codei SparseJacobian$$.
     110Furthermore, if you specify $icode work$$ in the calling sequence,
     111it is not necessary to keep the sparsity pattern; see the heading
     112$cref/p/sparse_jacobian/work/p/$$ under the $icode work$$ description.
     113$pre
     114
     115$$
    109116In addition,
    110117if you specify $icode p$$, CppAD will use the same
    111 type of sparsity representation 
     118type of sparsity representation
    112119(vectors of $code bool$$ or vectors of $code std::set<size_t>$$)
    113120for its internal calculations.
     
    149156$$
    150157In the case where the arguments $icode row$$ and $icode col$$ are present,
    151 we use $latex K$$ to denote the size of $icode jac$$. 
     158we use $latex K$$ to denote the size of $icode jac$$.
    152159The input value of its elements does not matter.
    153160Upon return, for $latex k = 0 , \ldots , K - 1$$,
     
    166173        sparse_jacobian_work& %work%
    167174%$$
    168 This object can only be used with the routines 
     175This object can only be used with the routines
    169176$code SparseJacobianForward$$ and $code SparseJacobianReverse$$.
    170177During its the first use, information is stored in $icode work$$.
    171178This is used to reduce the work done by future calls to the same mode
    172 (forward or reverse), 
     179(forward or reverse),
    173180the same $icode f$$, $icode p$$, $icode row$$, and $icode col$$.
    174181If a future call is for a different mode,
     
    184191        std::string %work%.color_method
    185192%$$
    186 and its default value (after a constructor or $code clear()$$) 
     193and its default value (after a constructor or $code clear()$$)
    187194is $code "cppad"$$.
    188195If $cref colpack_prefix$$ is specified on the
     
    193200$icode%work%.clear()%$$.
    194201
     202$subhead p$$
     203If $icode work$$ is present, and it is not the first call after
     204its construction or a clear,
     205the sparsity pattern $icode p$$ is not used.
     206This enables one to free the sparsity pattern
     207and still compute corresponding sparse Jacobians.
     208
    195209$head n_sweep$$
    196210The return value $icode n_sweep$$ has prototype
     
    198212        size_t %n_sweep%
    199213%$$
    200 If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used, 
    201 $icode n_sweep$$ is the number of first order forward (reverse) sweeps 
    202 used to compute the requested Jacobian values. 
     214If $code SparseJacobianForward$$ ($code SparseJacobianReverse$$) is used,
     215$icode n_sweep$$ is the number of first order forward (reverse) sweeps
     216used to compute the requested Jacobian values.
    203217(This is also the number of colors determined by the coloring method
    204218mentioned above).
    205 This is proportional to the total work that $code SparseJacobian$$ does, 
    206 not counting the zero order forward sweep, 
     219This is proportional to the total work that $code SparseJacobian$$ does,
     220not counting the zero order forward sweep,
    207221or the work to combine multiple columns (rows) into a single sweep.
    208222
     
    225239$subhead Restrictions$$
    226240If $icode VectorSet$$ has elements of $code std::set<size_t>$$,
    227 then $icode%p%[%i%]%$$ must return a reference (not a copy) to the 
     241then $icode%p%[%i%]%$$ must return a reference (not a copy) to the
    228242corresponding set.
    229243According to section 26.3.2.3 of the 1998 C++ standard,
    230244$code std::valarray< std::set<size_t> >$$ does not satisfy
    231 this condition. 
     245this condition.
    232246
    233247$head VectorSize$$
     
    240254$head Uses Forward$$
    241255After each call to $cref Forward$$,
    242 the object $icode f$$ contains the corresponding 
     256the object $icode f$$ contains the corresponding
    243257$cref/Taylor coefficients/glossary/Taylor Coefficient/$$.
    244258After a call to any of the sparse Jacobian routines,
     
    272286// ===========================================================================
    273287/*!
    274 class used by SparseJacobian to hold information so it does not need to be 
     288class used by SparseJacobian to hold information so it does not need to be
    275289recomputed.
    276290*/
     
    280294                /// (this field is set by user)
    281295                std::string color_method;
    282                 /// indices that sort the user row and col arrays by color 
     296                /// indices that sort the user row and col arrays by color
    283297                CppAD::vector<size_t> order;
    284298                /// results of the coloring algorithm
    285299                CppAD::vector<size_t> color;
    286        
     300
    287301                /// constructor
    288302                sparse_jacobian_work(void) : color_method("cppad")
     
    317331
    318332\param p_transpose [in]
    319 If <code>work.color.size() != 0</code>, 
     333If <code>work.color.size() != 0</code>,
    320334then \c p_transpose is not used.
    321335Otherwise, it is a
     
    341355this structure contains information that is computed by \c SparseJacobainFor.
    342356If the sparsity pattern, \c row vector, or \c col vectors
    343 are not the same between calls to \c SparseJacobianFor, 
     357are not the same between calls to \c SparseJacobianFor,
    344358\c work.clear() must be called to reinitialize \c work.
    345359
     
    450464
    451465                // set the corresponding components of the result
    452                 while( k < K && color[ col[order[k]] ] == ell ) 
     466                while( k < K && color[ col[order[k]] ] == ell )
    453467                {       jac[ order[k] ] = dy[row[order[k]]];
    454468                        k++;
     
    470484
    471485                // loop over colors we will do this tme
    472                 for(ell = 0; ell < r; ell++)   
     486                for(ell = 0; ell < r; ell++)
    473487                {       // combine all columns with this color
    474488                        for(j = 0; j < n; j++)
     
    482496
    483497                // store results
    484                 for(ell = 0; ell < r; ell++)   
     498                for(ell = 0; ell < r; ell++)
    485499                {       // set the components of the result for this color
    486                         while( k < K && color[ col[order[k]] ] == ell + count_color ) 
     500                        while( k < K && color[ col[order[k]] ] == ell + count_color )
    487501                        {       jac[ order[k] ] = dy[ row[order[k]] * r + ell ];
    488502                                k++;
     
    529543\param jac [out]
    530544is the vector of Jacobian values.
    531 It must have the same size as \c row. 
     545It must have the same size as \c row.
    532546The return value <code>jac[k]</code> is the partial of the
    533547<code>row[k]</code> range component of the function with respect
     
    538552This structure contains information that is computed by \c SparseJacobainRev.
    539553If the sparsity pattern, \c row vector, or \c col vectors
    540 are not the same between calls to \c SparseJacobianRev, 
     554are not the same between calls to \c SparseJacobianRev,
    541555\c work.clear() must be called to reinitialize \c work.
    542556
     
    621635        }
    622636        size_t n_color = 1;
    623         for(i = 0; i < m; i++) if( color[i] < m ) 
     637        for(i = 0; i < m; i++) if( color[i] < m )
    624638                n_color = std::max(n_color, color[i] + 1);
    625639
     
    649663
    650664                // set the corresponding components of the result
    651                 while( k < K && color[ row[order[k]] ]  == ell ) 
     665                while( k < K && color[ row[order[k]] ]  == ell )
    652666                {       jac[ order[k] ] = dw[col[order[k]]];
    653667                        k++;
     
    675689
    676690\tparam VectorSet
    677 is a simple vector class with elements of type 
     691is a simple vector class with elements of type
    678692\c bool or \c std::set<size_t>.
    679693
     
    696710\param jac [out]
    697711is the vector of Jacobian values.
    698 It must have the same size as \c row. 
     712It must have the same size as \c row.
    699713The return value <code>jac[k]</code> is the partial of the
    700714<code>row[k]</code> range component of the function with respect
     
    702716
    703717\param work [in,out]
    704 this structure contains information that depends on the function object, 
     718this structure contains information that depends on the function object,
    705719sparsity pattern, \c row vector, and \c col vector.
    706 If they are not the same between calls to \c SparseJacobianForward, 
     720If they are not the same between calls to \c SparseJacobianForward,
    707721\c work.clear() must be called to reinitialize them.
    708722
     
    730744                size_t(x.size()) == n ,
    731745                "SparseJacobianForward: size of x not equal domain dimension for f."
    732         ); 
     746        );
    733747        CPPAD_ASSERT_KNOWN(
    734748                size_t(row.size()) == K && size_t(col.size()) == K ,
    735749                "SparseJacobianForward: either r or c does not have "
    736750                "the same size as jac."
    737         ); 
     751        );
    738752        CPPAD_ASSERT_KNOWN(
    739753                work.color.size() == 0 || work.color.size() == n,
     
    783797
    784798\tparam VectorSet
    785 is a simple vector class with elements of type 
     799is a simple vector class with elements of type
    786800\c bool or \c std::set<size_t>.
    787801
     
    804818\param jac [out]
    805819is the vector of Jacobian values.
    806 It must have the same size as \c row. 
     820It must have the same size as \c row.
    807821The return value <code>jac[k]</code> is the partial of the
    808822<code>row[k]</code> range component of the function with respect
     
    810824
    811825\param work [in,out]
    812 this structure contains information that depends on the function object, 
     826this structure contains information that depends on the function object,
    813827sparsity pattern, \c row vector, and \c col vector.
    814 If they are not the same between calls to \c SparseJacobianReverse, 
     828If they are not the same between calls to \c SparseJacobianReverse,
    815829\c work.clear() must be called to reinitialize them.
    816830
     
    838852                size_t(x.size()) == n ,
    839853                "SparseJacobianReverse: size of x not equal domain dimension for f."
    840         ); 
     854        );
    841855        CPPAD_ASSERT_KNOWN(
    842856                size_t(row.size()) == K && size_t(col.size()) == K ,
    843857                "SparseJacobianReverse: either r or c does not have "
    844858                "the same size as jac."
    845         ); 
     859        );
    846860        CPPAD_ASSERT_KNOWN(
    847861                work.color.size() == 0 || work.color.size() == m,
     
    864878        );
    865879# endif
    866  
     880
    867881        typedef typename VectorSet::value_type Set_type;
    868882        typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
     
    891905
    892906\tparam VectorSet
    893 is a simple vector class with elements of type 
     907is a simple vector class with elements of type
    894908\c bool or \c std::set<size_t>.
    895909
     
    950964                                i = s_transpose.next_element();
    951965                        }
    952                 } 
     966                }
    953967                size_t K = k;
    954968                VectorBase J(K);
    955        
     969
    956970                // now we have folded this into the following case
    957971                SparseJacobianFor(x, s_transpose, row, col, J, work);
     
    978992                                j = s.next_element();
    979993                        }
    980                 } 
     994                }
    981995                size_t K = k;
    982996                VectorBase J(K);
     
    9911005
    9921006        return jac;
    993 } 
     1007}
    9941008
    9951009/*!
     
    10291043        {       size_t j, k;
    10301044
    1031                 // use forward mode 
     1045                // use forward mode
    10321046                VectorBool r(n * n);
    10331047                for(j = 0; j < n; j++)
     
    10411055        {       size_t i, k;
    10421056
    1043                 // use reverse mode 
     1057                // use reverse mode
    10441058                VectorBool s(m * m);
    10451059                for(i = 0; i < m; i++)
  • trunk/example/sparse_hessian.cpp

    r2506 r3673  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
    6 the terms of the 
     6the terms of the
    77                    Eclipse Public License Version 1.0.
    88
     
    162162        ell        = 1 * n + 1;
    163163        check[ell] = 6.0 * x[1];
    164         n_sweep    = f.SparseHessian(x, w, p_set, row, col, hes, work);
     164        s_vector   not_used;
     165        n_sweep    = f.SparseHessian(x, w, not_used, row, col, hes, work);
    165166        for(k = 0; k < K; k++)
    166167        {       ell = row[k] * n + col[k];
     
    168169        }
    169170        ok &= n_sweep == 2;
    170        
     171
    171172
    172173
  • trunk/example/sparse_jacobian.cpp

    r2506 r3673  
    11/* $Id$ */
    22/* --------------------------------------------------------------------------
    3 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
     3CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-15 Bradley M. Bell
    44
    55CppAD is distributed under multiple licenses. This distribution is under
    6 the terms of the 
     6the terms of the
    77                    Eclipse Public License Version 1.0.
    88
     
    105105
    106106        // using row and column indices to compute non-zero in rows 1 and 2
    107         // (skip row 0). 
     107        // (skip row 0).
    108108        size_t K = 6;
    109109        i_vector row(K), col(K);
     
    120120                        }
    121121                }
    122         } 
     122        }
    123123        ok &= k == K;
    124124
     
    126126        CppAD::sparse_jacobian_work work;
    127127
    128         // could use p_b 
     128        // could use p_b
    129129        size_t n_sweep = f.SparseJacobianReverse(x, p_s, row, col, jac, work);
    130130        for(k = 0; k < K; k++)
     
    136136        // now recompute at a different x value (using work from previous call)
    137137        check[11] = x[3] = 10.;
    138         n_sweep = f.SparseJacobianReverse(x, p_s, row, col, jac, work);
     138        std::vector< std::set<size_t> > not_used;
     139        n_sweep = f.SparseJacobianReverse(x, not_used, row, col, jac, work);
    139140        for(k = 0; k < K; k++)
    140141        {       ell = row[k] * n + col[k];
     
    190191        */
    191192        d_vector check(m * n);
    192         check[0] = 1.; check[1]  = 0.; check[2]  = 1.; 
     193        check[0] = 1.; check[1]  = 0.; check[2]  = 1.;
    193194        check[3] = 1.; check[4]  = 0.; check[5]  = 1.;
    194         check[6] = 0.; check[7]  = 1.; check[8]  = 1.; 
     195        check[6] = 0.; check[7]  = 1.; check[8]  = 1.;
    195196        check[9] = 0.; check[10] = 1.; check[11] = x[2];
    196197        for(ell = 0; ell < size_t(check.size()); ell++)
     
    219220
    220221        // using row and column indices to compute non-zero elements excluding
    221         // row 0 and column 0. 
     222        // row 0 and column 0.
    222223        size_t K = 5;
    223224        i_vector row(K), col(K);
     
    234235                        }
    235236                }
    236         } 
     237        }
    237238        ok &= k == K;
    238239
     
    240241        CppAD::sparse_jacobian_work work;
    241242
    242         // could use p_s 
     243        // could use p_s
    243244        size_t n_sweep = f.SparseJacobianForward(x, p_b, row, col, jac, work);
    244245        for(k = 0; k < K; k++)
     
    259260        return ok;
    260261}
    261 } // End empty namespace 
     262} // End empty namespace
    262263
    263264bool sparse_jacobian(void)
  • trunk/omh/whats_new/whats_new_15.omh

    r3670 r3673  
    4949assist you in learning about changes between various versions of CppAD.
    5050
     51$head 04-18$$
     52In the sparse jacobian and sparse hessian calculations,
     53If $icode work$$ is present, and has already been computed,
     54the sparsity pattern $icode p$$ is not used.
     55This has been added to the documentation; see
     56$cref/sparse jacobian/sparse_jacobian/work/p/$$ and
     57$cref/sparse hessian/sparse_hessian/work/p/$$ documentation
     58for $icode work$$ and $icode p$$.
     59
    5160$head 03-13$$
    5261Remove the syntax
Note: See TracChangeset for help on using the changeset viewer.