CppAD: A C++ Algorithmic Differentiation Package 20110419
jac_g_map.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 "jac_g_map.hpp"
00014 CPPAD_BEGIN_NAMESPACE
00015 /*!
00016 \file jac_g_map.cpp
00017 \brief Creates a mapping between two representations for Jacobian of g. 
00018 */
00019 
00020 /*!
00021 Create mapping from CppAD to Ipopt sparse representations of Jacobian of g.
00022 
00023 The functions 
00024 \f$ f : {\bf R}^n \rightarrow {\bf R} \f$ and
00025 \f$ g : {\bf R}^n \rightarrow {\bf R}^m \f$ are defined by
00026 the \ref Users_Representation.
00027 
00028 \param fg_info
00029 For <tt>k = 0 , ... , K-1</tt>,
00030 for <tt>ell = 0 , ... , L[k]</tt>,
00031 the function call 
00032 \verbatim
00033         fg_info->index(k, ell, I, J);   
00034 \endverbatim
00035 is made by \c jac_g_map. 
00036 The values \c k and \c ell are inputs.
00037 The input size of \c I ( \c J ) 
00038 is greater than or equal <tt>p[k] ( q[k] )</tt>
00039 and this size is not changed.
00040 The input values of the elements of \c I and \c J are not specified.
00041 The output value of the elements of \c I define
00042 \f[
00043 I_{k, \ell} = ( {\rm I[0]} , \cdots , {\rm I[p[k]-1]} )
00044 \f]
00045 The output value of the elements of \c J define
00046 \f[
00047 J_{k, \ell} = ( {\rm J[0]} , \cdots , {\rm J[q[k]-1]} )
00048 \f]
00049 
00050 \param m
00051 is the dimension of the range space for \f$ g(x) \f$; i.e.,
00052 \f$ g(x) \in {\bf R}^m \f$.
00053 
00054 \param n
00055 is the dimension of the domain space for \f$ f(x) \f$ and \f$ g(x) \f$;
00056 i.e., \f$ x \in {\bf R}^n \f$.
00057 
00058 \param K
00059 is the number of functions \f$ r_k ( u ) \f$ used for the representation of 
00060 \f$ f(x) \f$ and \f$ g(x) \f$.
00061 
00062 \param L
00063 is a vector with size \c K.
00064 For <tt>k = 0 , ... , K-1, L[k]</tt>
00065 is the number of terms that use \f$ r_k (u) \f$
00066 in the representation of \f$ f(x) \f$ and \f$ g(x) \f$.
00067 
00068 \param p
00069 is a vector with size \c K.
00070 For <tt>k = 0 , ... , K-1, p[k]</tt>
00071 is dimension of the range space for \f$ r_k (u) \f$; i.e.,
00072 \f$ r_k (u) \in {\bf R}^{p(k)} \f$.
00073 
00074 \param q
00075 is a vector with size \c K.
00076 For <tt>k = 0 , ... , K-1, q[k]</tt>
00077 is dimension of the domain space for \f$ r_k (u) \f$; i.e.,
00078 \f$ u \in {\bf R}^{q(k)} \f$.
00079 
00080 \param pattern_jac_r
00081 is a vector with size \c K.
00082 For <tt>k = 0 , ... , K-1, pattern_jac_r[k]</tt>
00083 is a CppAD sparsity pattern for the Jacobian of the function 
00084 \f$ r_k : {\bf R}^{q(k)} \rightarrow {\bf R}^{p(k)} \f$.
00085 As such, <tt>pattern_jac_r[k].size() == p[k] * q[k]</tt>.
00086 
00087 \param I
00088 is a work vector of length greater than or equal <tt>p[k]</tt> for all \c k.
00089 The input and output value of its elements are unspecified.
00090 The size of \c I is not changed.
00091 
00092 \param J
00093 is a work vector of length greater than or equal <tt>q[k]</tt> for all \c k.
00094 The input and output value of its elements are unspecified.
00095 The size of \c J is not changed.
00096 
00097 \param index_jac_g:
00098 On input, this is empty; i.e., <tt>index_jac_g.size() == 0</tt>.
00099 On output, it is the index mapping from \f$ (i, j) \f$ in the Jacobian of 
00100 \f$ g(x) \f$ to the corresponding index value used by Ipopt to represent
00101 the Jacobian.
00102 Furthermore, if <tt>index_jac_g[i].find(j) == index_jac_g[i].end()</tt>,
00103 then the \f$ (i, j)\f$ entry in the Jacobian of \f$ g(x) \f$ is always zero.
00104 */
00105 void jac_g_map(
00106         cppad_ipopt_fg_info*  fg_info                                  , 
00107         size_t                                          m              ,
00108         size_t                                          n              ,
00109         size_t                                          K              ,
00110         const CppAD::vector<size_t>&                    L              ,
00111         const CppAD::vector<size_t>&                    p              ,
00112         const CppAD::vector<size_t>&                    q              ,
00113         const CppAD::vector<CppAD::vectorBool>&         pattern_jac_r  ,
00114         CppAD::vector<size_t>&                          I              ,
00115         CppAD::vector<size_t>&                          J              ,
00116         CppAD::vector< std::map<size_t,size_t> >&       index_jac_g    )
00117 {
00118         using CppAD::vectorBool;
00119         size_t i, j, ij, k, ell;
00120 
00121         CPPAD_ASSERT_UNKNOWN( K == L.size() );
00122         CPPAD_ASSERT_UNKNOWN( K == p.size() );
00123         CPPAD_ASSERT_UNKNOWN( K == q.size() );
00124         CPPAD_ASSERT_UNKNOWN( K == pattern_jac_r.size() );
00125 # ifndef NDEBUG
00126         for(k = 0; k < K; k++)
00127         {       CPPAD_ASSERT_UNKNOWN( p[k] <= I.size() );
00128                 CPPAD_ASSERT_UNKNOWN( q[k] <= J.size() );
00129                 CPPAD_ASSERT_UNKNOWN( p[k]*q[k] == pattern_jac_r[k].size() );
00130         }
00131 # endif
00132         // Now compute pattern for g
00133         // (use standard set representation because can be huge).
00134         CppAD::vector< std::set<size_t> > pattern_jac_g(m);
00135         for(k = 0; k < K; k++) for(ell = 0; ell < L[k]; ell++)
00136         {       fg_info->index(k, ell, I, J);   
00137                 for(i = 0; i < p[k]; i++) if( I[i] != 0 )
00138                 {       for(j = 0; j < q[k]; j++)
00139                         {       ij  = i * q[k] + j;
00140                                 if( pattern_jac_r[k][ij] )
00141                                         pattern_jac_g[I[i]-1].insert(J[j]);
00142                         }
00143                 }
00144         }
00145 
00146         // Now compute the mapping from (i, j) in the Jacobian of g to the 
00147         // corresponding index value used by Ipopt to represent the Jacobian.
00148         CPPAD_ASSERT_UNKNOWN( index_jac_g.size() == 0 );
00149         index_jac_g.resize(m);
00150         std::set<size_t>::const_iterator itr;
00151         ell = 0;
00152         for(i = 0; i < m; i++)
00153         {       for(    itr = pattern_jac_g[i].begin(); 
00154                         itr != pattern_jac_g[i].end(); itr++)
00155                 {
00156                         index_jac_g[i][*itr] = ell++;
00157                 }
00158         }
00159         return;
00160 }
00161 
00162 CPPAD_END_NAMESPACE