CppAD: A C++ Algorithmic Differentiation Package  20130918
lu_ratio.hpp
Go to the documentation of this file.
00001 /* $Id$ */
00002 # ifndef CPPAD_LU_RATIO_INCLUDED
00003 # define CPPAD_LU_RATIO_INCLUDED
00004 
00005 /* --------------------------------------------------------------------------
00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-12 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Eclipse Public License Version 1.0.
00011 
00012 A copy of this license is included in the COPYING file of this distribution.
00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
00014 -------------------------------------------------------------------------- */
00015 
00016 /*
00017 $begin LuRatio$$
00018 $spell
00019      cppad.hpp
00020      xk
00021      Cpp
00022      Lu
00023      bool
00024      const
00025      ip
00026      jp
00027      std
00028      ADvector
00029 $$
00030 
00031 $index LuRatio$$
00032 $index linear, Lu factor equation$$
00033 $index equation, Lu factor$$
00034 $index determinant, Lu factor$$
00035 $index solve, Lu factor$$
00036 
00037 $section LU Factorization of A Square Matrix and Stability Calculation$$
00038 
00039 $head Syntax$$
00040 $code# include <cppad/cppad.hpp>$$
00041 $pre
00042 $$
00043 $icode%sign% = LuRatio(%ip%, %jp%, %LU%, %ratio%)%$$
00044 
00045 
00046 $head Description$$
00047 Computes an LU factorization of the matrix $icode A$$ 
00048 where $icode A$$ is a square matrix.
00049 A measure of the numerical stability called $icode ratio$$ is calculated.
00050 This ratio is useful when the results of $code LuRatio$$ are
00051 used as part of an $cref ADFun$$ object.
00052 
00053 $head Include$$
00054 This routine is designed to be used with AD objects and
00055 requires the $code cppad/cppad.hpp$$ file to be included.
00056 
00057 $head Matrix Storage$$
00058 All matrices are stored in row major order.
00059 To be specific, if $latex Y$$ is a vector
00060 that contains a $latex p$$ by $latex q$$ matrix,
00061 the size of $latex Y$$ must be equal to $latex  p * q $$ and for
00062 $latex i = 0 , \ldots , p-1$$,
00063 $latex j = 0 , \ldots , q-1$$,
00064 $latex \[
00065      Y_{i,j} = Y[ i * q + j ]
00066 \] $$
00067 
00068 $head sign$$
00069 The return value $icode sign$$ has prototype
00070 $codei%
00071      int %sign%
00072 %$$
00073 If $icode A$$ is invertible, $icode sign$$ is plus or minus one
00074 and is the sign of the permutation corresponding to the row ordering
00075 $icode ip$$ and column ordering $icode jp$$.
00076 If $icode A$$ is not invertible, $icode sign$$ is zero.
00077 
00078 $head ip$$
00079 The argument $icode ip$$ has prototype
00080 $codei%
00081      %SizeVector% &%ip%
00082 %$$
00083 (see description of $cref/SizeVector/LuFactor/SizeVector/$$ below).
00084 The size of $icode ip$$ is referred to as $icode n$$ in the
00085 specifications below.
00086 The input value of the elements of $icode ip$$ does not matter.
00087 The output value of the elements of $icode ip$$ determine
00088 the order of the rows in the permuted matrix.
00089 
00090 $head jp$$
00091 The argument $icode jp$$ has prototype
00092 $codei%
00093      %SizeVector% &%jp%
00094 %$$
00095 (see description of $cref/SizeVector/LuFactor/SizeVector/$$ below).
00096 The size of $icode jp$$ must be equal to $icode n$$.
00097 The input value of the elements of $icode jp$$ does not matter.
00098 The output value of the elements of $icode jp$$ determine
00099 the order of the columns in the permuted matrix.
00100 
00101 $head LU$$
00102 The argument $icode LU$$ has the prototype
00103 $codei%
00104      %ADvector% &%LU%
00105 %$$
00106 and the size of $icode LU$$ must equal $latex n * n$$
00107 (see description of $cref/ADvector/LuRatio/ADvector/$$ below).
00108 
00109 $subhead A$$
00110 We define $icode A$$ as the matrix corresponding to the input 
00111 value of $icode LU$$.
00112 
00113 $subhead P$$
00114 We define the permuted matrix $icode P$$ in terms of $icode A$$ by
00115 $codei%
00116      %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00117 %$$
00118 
00119 $subhead L$$
00120 We define the lower triangular matrix $icode L$$ in terms of the 
00121 output value of $icode LU$$.
00122 The matrix $icode L$$ is zero above the diagonal
00123 and the rest of the elements are defined by
00124 $codei%
00125      %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00126 %$$
00127 for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$.
00128 
00129 $subhead U$$
00130 We define the upper triangular matrix $icode U$$ in terms of the
00131 output value of $icode LU$$.
00132 The matrix $icode U$$ is zero below the diagonal,
00133 one on the diagonal,
00134 and the rest of the elements are defined by
00135 $codei%
00136      %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00137 %$$
00138 for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$.
00139 
00140 $subhead Factor$$
00141 If the return value $icode sign$$ is non-zero,
00142 $codei%
00143      %L% * %U% = %P%
00144 %$$
00145 If the return value of $icode sign$$ is zero,
00146 the contents of $icode L$$ and $icode U$$ are not defined. 
00147 
00148 $subhead Determinant$$
00149 $index determinant$$
00150 If the return value $icode sign$$ is zero,
00151 the determinant of $icode A$$ is zero.
00152 If $icode sign$$ is non-zero,
00153 using the output value of $icode LU$$
00154 the determinant of the matrix $icode A$$ is equal to
00155 $codei%
00156 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 
00157 %$$
00158 
00159 $head ratio$$
00160 The argument $icode ratio$$ has prototype
00161 $codei%
00162         AD<%Base%> &%ratio%
00163 %$$
00164 On input, the value of $icode ratio$$ does not matter.
00165 On output it is a measure of how good the choice of pivots is.
00166 For $latex p = 0 , \ldots , n-1$$, 
00167 the $th p$$ pivot element is the element of maximum absolute value of a 
00168 $latex (n-p) \times (n-p)$$ sub-matrix.
00169 The ratio of each element of sub-matrix divided by the pivot element
00170 is computed.
00171 The return value of $icode ratio$$ is the maximum absolute value of
00172 such ratios over with respect to all elements and all the pivots.
00173 
00174 $subhead Purpose$$
00175 Suppose that the execution of a call to $code LuRatio$$ 
00176 is recorded in the $codei%ADFun<%Base%>%$$ object $icode F$$.
00177 Then a call to $cref Forward$$ of the form
00178 $codei%
00179      %F%.Forward(%k%, %xk%)
00180 %$$
00181 with $icode k$$ equal to zero will revaluate this Lu factorization
00182 with the same pivots and a new value for $icode A$$.
00183 In this case, the resulting $icode ratio$$ may not be one.
00184 If $icode ratio$$ is too large (the meaning of too large is up to you), 
00185 the current pivots do not yield a stable LU factorization of $icode A$$.
00186 A better choice for the pivots (for this value of $icode A$$)
00187 will be made if you recreate the $code ADFun$$ object
00188 starting with the $cref Independent$$ variable values
00189 that correspond to the vector $icode xk$$.
00190 
00191 $head SizeVector$$
00192 The type $icode SizeVector$$ must be a $cref SimpleVector$$ class with
00193 $cref/elements of type size_t/SimpleVector/Elements of Specified Type/$$.
00194 The routine $cref CheckSimpleVector$$ will generate an error message
00195 if this is not the case.
00196 
00197 $head ADvector$$
00198 The type $icode ADvector$$ must be a 
00199 $cref/simple vector class/SimpleVector/$$ with elements of type
00200 $codei%AD<%Base%>%$$.
00201 The routine $cref CheckSimpleVector$$ will generate an error message
00202 if this is not the case.
00203 
00204 
00205 $head Example$$
00206 $children%
00207      example/lu_ratio.cpp
00208 %$$
00209 The file $cref lu_ratio.cpp$$
00210 contains an example and test of using $code LuRatio$$.
00211 It returns true if it succeeds and false otherwise.
00212 
00213 $end
00214 --------------------------------------------------------------------------
00215 */
00216 namespace CppAD { // BEGIN CppAD namespace
00217 
00218 // Lines different from the code in cppad/lu_factor.hpp end with           //
00219 template <class SizeVector, class ADvector, class Base>                    //
00220 int LuRatio(SizeVector &ip, SizeVector &jp, ADvector &LU, AD<Base> &ratio) //
00221 {    
00222      typedef ADvector FloatVector;                                       //
00223      typedef AD<Base>       Float;                                       //
00224 
00225      // check numeric type specifications
00226      CheckNumericType<Float>();
00227 
00228      // check simple vector class specifications
00229      CheckSimpleVector<Float, FloatVector>();
00230      CheckSimpleVector<size_t, SizeVector>();
00231 
00232      size_t  i, j;          // some temporary indices
00233      const Float zero( 0 ); // the value zero as a Float object
00234      size_t  imax;          // row index of maximum element
00235      size_t  jmax;          // column indx of maximum element
00236      Float    emax;         // maximum absolute value
00237      size_t  p;             // count pivots
00238      int     sign;          // sign of the permutation
00239      Float   etmp;          // temporary element
00240      Float   pivot;         // pivot element
00241 
00242      // -------------------------------------------------------
00243      size_t n = size_t(ip.size());
00244      CPPAD_ASSERT_KNOWN(
00245           size_t(jp.size()) == n,
00246           "Error in LuFactor: jp must have size equal to n"
00247      );
00248      CPPAD_ASSERT_KNOWN(
00249           size_t(LU.size()) == n * n,
00250           "Error in LuFactor: LU must have size equal to n * m"
00251      );
00252      // -------------------------------------------------------
00253 
00254      // initialize row and column order in matrix not yet pivoted
00255      for(i = 0; i < n; i++)
00256      {    ip[i] = i;
00257           jp[i] = i;
00258      }
00259      // initialize the sign of the permutation
00260      sign = 1;
00261      // initialize the ratio                                             //
00262      ratio = Float(1);                                                   //
00263      // ---------------------------------------------------------
00264 
00265      // Reduce the matrix P to L * U using n pivots
00266      for(p = 0; p < n; p++)
00267      {    // determine row and column corresponding to element of 
00268           // maximum absolute value in remaining part of P
00269           imax = jmax = n;
00270           emax = zero;
00271           for(i = p; i < n; i++)
00272           {    for(j = p; j < n; j++)
00273                {    CPPAD_ASSERT_UNKNOWN(
00274                          (ip[i] < n) & (jp[j] < n)
00275                     );
00276                     etmp = LU[ ip[i] * n + jp[j] ];
00277 
00278                     // check if maximum absolute value so far
00279                     if( AbsGeq (etmp, emax) )
00280                     {    imax = i;
00281                          jmax = j;
00282                          emax = etmp;
00283                     }
00284                }
00285           }
00286           for(i = p; i < n; i++)                                       //
00287           {    for(j = p; j < n; j++)                               //
00288                {    etmp  = abs(LU[ ip[i] * n + jp[j] ] / emax); //
00289                     ratio =                                      //
00290                     CondExpGt(etmp, ratio, etmp, ratio);         //
00291                }                                                    //
00292           }                                                            //
00293           CPPAD_ASSERT_KNOWN( 
00294                (imax < n) & (jmax < n) ,
00295                "AbsGeq must return true when second argument is zero"
00296           );
00297           if( imax != p )
00298           {    // switch rows so max absolute element is in row p
00299                i        = ip[p];
00300                ip[p]    = ip[imax];
00301                ip[imax] = i;
00302                sign     = -sign;
00303           }
00304           if( jmax != p )
00305           {    // switch columns so max absolute element is in column p
00306                j        = jp[p];
00307                jp[p]    = jp[jmax];
00308                jp[jmax] = j;
00309                sign     = -sign;
00310           }
00311           // pivot using the max absolute element
00312           pivot   = LU[ ip[p] * n + jp[p] ];
00313 
00314           // check for determinant equal to zero
00315           if( pivot == zero )
00316           {    // abort the mission
00317                return   0;
00318           }
00319 
00320           // Reduce U by the elementary transformations that maps 
00321           // LU( ip[p], jp[p] ) to one.  Only need transform elements
00322           // above the diagonal in U and LU( ip[p] , jp[p] ) is
00323           // corresponding value below diagonal in L.
00324           for(j = p+1; j < n; j++)
00325                LU[ ip[p] * n + jp[j] ] /= pivot;
00326 
00327           // Reduce U by the elementary transformations that maps 
00328           // LU( ip[i], jp[p] ) to zero. Only need transform elements 
00329           // above the diagonal in U and LU( ip[i], jp[p] ) is 
00330           // corresponding value below diagonal in L.
00331           for(i = p+1; i < n; i++ )
00332           {    etmp = LU[ ip[i] * n + jp[p] ];
00333                for(j = p+1; j < n; j++)
00334                {    LU[ ip[i] * n + jp[j] ] -= 
00335                          etmp * LU[ ip[p] * n + jp[j] ];
00336                } 
00337           }
00338      }
00339      return sign;
00340 }
00341 } // END CppAD namespace 
00342 
00343 # endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines