CppAD: A C++ Algorithmic Differentiation Package 20110419
vec_fun_pattern.cpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 /* --------------------------------------------------------------------------
00003 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-10 Bradley M. Bell
00004 
00005 CppAD is distributed under multiple licenses. This distribution is under
00006 the terms of the 
00007                     Common Public License Version 1.0.
00008 
00009 A copy of this license is included in the COPYING file of this distribution.
00010 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00011 -------------------------------------------------------------------------- */
00012 # include "cppad_ipopt_nlp.hpp"
00013 # include "vec_fun_pattern.hpp"
00014 CPPAD_BEGIN_NAMESPACE
00015 /*!
00016 \file vec_fun_pattern.cpp
00017 \brief Determine a sparsity pattern for a vector of AD function objects.
00018 */
00019 
00020 /*!
00021 Determine a sparsity patterns for each function in a vector of functions.
00022 
00023 \param K
00024 is the number of functions that we are computing the sparsity pattern for.
00025 
00026 \param p
00027 is a vector with size \c K.
00028 For <tt>k = 0 , ... , K-1, p[k]</tt>
00029 is dimension of the range space for \f$ r_k (u) \f$; i.e.,
00030 \f$ r_k (u) \in {\bf R}^{p(k)} \f$.
00031 
00032 \param q
00033 is a vector with size \c K.
00034 For <tt>k = 0 , ... , K-1, q[k]</tt>
00035 is dimension of the domain space for \f$ r_k (u) \f$; i.e.,
00036 \f$ u \in {\bf R}^{q(k)} \f$.
00037 
00038 \param retape
00039 is a vector with size \c K.
00040 For <tt>k = 0 , ... , K-1</tt>,
00041 if <tt>retape[k]</tt> is true,
00042 the function object <tt>r[k]</tt> is a valid representation
00043 for \f$ r_k (u) \f$ for all \f$ u \in {\bf R}^{q(k)} \f$.
00044 Otherwise, the function object must be retaped for each
00045 value of \f$ u \f$.
00046 
00047 \param r_fun
00048 is the vector of AD function objects which has size size \c K.
00049 For <tt>k = 0 , ... , K-1</tt>,
00050 if <tt>retape[k]</tt> is true, <tt>r_fun[k]</tt> is not used.
00051 If <tt>retape[k]</tt> is false, <tt>r_fun[k]</tt> is not used.
00052 is a CppAD function object correspopnding to the function
00053 \f$ r_k : {\bf R}^{q[k]} \rightarrow {\bf R}^{p[k]} \f$.
00054 The following non-constant member functions will be called:
00055 \verbatim
00056         r_fun[k].ForSparseJac(q[k], pattern_domain)
00057         r_fun[k].RevSparseHes(p[k], pattern_range)
00058 \endverbatim
00059 The following \c const member functions <tt>r_fun[k].Range()</tt>
00060 and <tt>r_fun[k].Domain()</tt> may also be called.
00061 
00062 \param pattern_jac_r
00063 is a vector with size \c K.
00064 On input, For <tt>k = 0 , ... , K-1, pattern_jac_r[k]</tt>
00065 is a vector of length p[k] * q[k]
00066 and the value of its elements does not matter.
00067 On output it is a CppAD sparsity pattern for the Jacobian of 
00068 \f$ r_k (u) \f$.
00069 
00070 \param pattern_hes_r
00071 is a vector with size \c K.
00072 On input, For <tt>k = 0 , ... , K-1, pattern_hes_r[k]</tt>
00073 is a vector of length q[k] * q[k]
00074 and the value of its elements does not matter.
00075 On output it is a CppAD sparsity pattern for the Hessian of 
00076 \f$ R : {\bf R}^{q[k]} \rightarrow {\bf R} \f$ which is defined by
00077 \f[
00078         R(u) = \sum_{i=0}^{p[k]-1} r_k (u)_i
00079 \f]
00080 */
00081 void vec_fun_pattern(
00082         size_t                                          K              ,
00083         const CppAD::vector<size_t>&                    p              ,
00084         const CppAD::vector<size_t>&                    q              ,
00085         const CppAD::vectorBool&                        retape         ,
00086         CppAD::vector< CppAD::ADFun<Ipopt::Number> >&   r_fun          ,
00087         CppAD::vector<CppAD::vectorBool>&               pattern_jac_r  ,
00088         CppAD::vector<CppAD::vectorBool>&               pattern_hes_r  )
00089 {       // check some assumptions
00090         CPPAD_ASSERT_UNKNOWN( K == p.size() );
00091         CPPAD_ASSERT_UNKNOWN( K == q.size() );
00092         CPPAD_ASSERT_UNKNOWN( K == retape.size() );
00093         CPPAD_ASSERT_UNKNOWN( K == r_fun.size() );
00094         CPPAD_ASSERT_UNKNOWN( K == pattern_jac_r.size() );
00095         CPPAD_ASSERT_UNKNOWN( K == pattern_hes_r.size() );
00096 
00097         using CppAD::vectorBool;
00098         size_t i, j, k;
00099 
00100         for(k = 0; k < K; k++)
00101         {       // check some k specific assumptions
00102                 CPPAD_ASSERT_UNKNOWN( pattern_jac_r[k].size() == p[k] * q[k] );
00103                 CPPAD_ASSERT_UNKNOWN( pattern_hes_r[k].size() == q[k] * q[k] );
00104 
00105                 if( retape[k] )
00106                 {       for(i = 0; i < p[k]; i++)
00107                         {       for(j = 0; j < q[k]; j++)
00108                                         pattern_jac_r[k][i*q[k] + j] = true; 
00109                         }
00110                         for(i = 0; i < q[k]; i++)
00111                         {       for(j = 0; j < q[k]; j++)
00112                                         pattern_hes_r[k][i*q[k] + j] = true; 
00113                         }
00114                 }
00115                 else 
00116                 {       // check assumptions about r_k
00117                         CPPAD_ASSERT_UNKNOWN( r_fun[k].Range() == p[k] );
00118                         CPPAD_ASSERT_UNKNOWN( r_fun[k].Domain() == q[k] );
00119 
00120                         // pattern for the identity matrix
00121                         CppAD::vectorBool pattern_domain(q[k] * q[k]);
00122                         for(i = 0; i < q[k]; i++)
00123                         {       for(j = 0; j < q[k]; j++) 
00124                                         pattern_domain[i*q[k] + j] = (i == j);
00125                         }
00126                         // use forward mode to compute Jacobian sparsity
00127                         pattern_jac_r[k] = 
00128                                 r_fun[k].ForSparseJac(q[k], pattern_domain);
00129                         // user reverse mode to compute Hessian sparsity
00130                         CppAD::vectorBool pattern_ones(p[k]);
00131                         for(i = 0; i < p[k]; i++)
00132                                 pattern_ones[i] = true;
00133                         pattern_hes_r[k] =
00134                                 r_fun[k].RevSparseHes(q[k], pattern_ones);
00135                 }
00136         }
00137 }
00138 CPPAD_END_NAMESPACE