CppAD: A C++ Algorithmic Differentiation Package 20110419
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-08 Bradley M. Bell
00007 
00008 CppAD is distributed under multiple licenses. This distribution is under
00009 the terms of the 
00010                     Common 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 $syntax%%sign% = LuRatio(%ip%, %jp%, %LU%, %ratio%)%$$
00044 
00045 
00046 $head Description$$
00047 Computes an LU factorization of the matrix $italic A$$ 
00048 where $italic A$$ is a square matrix.
00049 A measure of the numerical stability called $italic ratio$$ is calculated.
00050 This ratio is useful when the results of $code LuRatio$$ are
00051 used as part of an $xref/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 $italic sign$$ has prototype
00070 $syntax%
00071         int %sign%
00072 %$$
00073 If $italic A$$ is invertible, $italic sign$$ is plus or minus one
00074 and is the sign of the permutation corresponding to the row ordering
00075 $italic ip$$ and column ordering $italic jp$$.
00076 If $italic A$$ is not invertible, $italic sign$$ is zero.
00077 
00078 $head ip$$
00079 The argument $italic ip$$ has prototype
00080 $syntax%
00081         %SizeVector% &%ip%
00082 %$$
00083 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below).
00084 The size of $italic ip$$ is referred to as $italic n$$ in the
00085 specifications below.
00086 The input value of the elements of $italic ip$$ does not matter.
00087 The output value of the elements of $italic ip$$ determine
00088 the order of the rows in the permuted matrix.
00089 
00090 $head jp$$
00091 The argument $italic jp$$ has prototype
00092 $syntax%
00093         %SizeVector% &%jp%
00094 %$$
00095 (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below).
00096 The size of $italic jp$$ must be equal to $italic n$$.
00097 The input value of the elements of $italic jp$$ does not matter.
00098 The output value of the elements of $italic jp$$ determine
00099 the order of the columns in the permuted matrix.
00100 
00101 $head LU$$
00102 The argument $italic LU$$ has the prototype
00103 $syntax%
00104         %ADvector% &%LU%
00105 %$$
00106 and the size of $italic LU$$ must equal $latex n * n$$
00107 (see description of $xref/LuRatio/ADvector/ADvector/$$ below).
00108 
00109 $subhead A$$
00110 We define $italic A$$ as the matrix corresponding to the input 
00111 value of $italic LU$$.
00112 
00113 $subhead P$$
00114 We define the permuted matrix $italic P$$ in terms of $italic A$$ by
00115 $syntax%
00116         %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ]
00117 %$$
00118 
00119 $subhead L$$
00120 We define the lower triangular matrix $italic L$$ in terms of the 
00121 output value of $italic LU$$.
00122 The matrix $italic L$$ is zero above the diagonal
00123 and the rest of the elements are defined by
00124 $syntax%
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 $italic U$$ in terms of the
00131 output value of $italic LU$$.
00132 The matrix $italic U$$ is zero below the diagonal,
00133 one on the diagonal,
00134 and the rest of the elements are defined by
00135 $syntax%
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 $italic sign$$ is non-zero,
00142 $syntax%
00143         %L% * %U% = %P%
00144 %$$
00145 If the return value of $italic sign$$ is zero,
00146 the contents of $italic L$$ and $italic U$$ are not defined. 
00147 
00148 $subhead Determinant$$
00149 $index determinant$$
00150 If the return value $italic sign$$ is zero,
00151 the determinant of $italic A$$ is zero.
00152 If $italic sign$$ is non-zero,
00153 using the output value of $italic LU$$
00154 the determinant of the matrix $italic A$$ is equal to
00155 $syntax%
00156 %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] 
00157 %$$
00158 
00159 $head ratio$$
00160 The argument $italic ratio$$ has prototype
00161 $syntax%
00162         AD<%Base%> &%ratio%
00163 %$$
00164 On input, the value of $italic 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 $italic 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 $syntax%ADFun<%Base%>%$$ object $italic F$$.
00177 Then a call to $xref/Forward/$$ of the form
00178 $syntax%
00179         %F%.Forward(%k%, %xk%)
00180 %$$
00181 with $italic k$$ equal to zero will revaluate this Lu factorization
00182 with the same pivots and a new value for $italic A$$.
00183 In this case, the resulting $italic ratio$$ may not be one.
00184 If $italic 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 $italic A$$.
00186 A better choice for the pivots (for this value of $italic A$$)
00187 will be made if you recreate the $code ADFun$$ object
00188 starting with the $xref/Independent/$$ variable values
00189 that correspond to the vector $italic xk$$.
00190 
00191 $head SizeVector$$
00192 The type $italic SizeVector$$ must be a $xref/SimpleVector/$$ class with
00193 $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$.
00194 The routine $xref/CheckSimpleVector/$$ will generate an error message
00195 if this is not the case.
00196 
00197 $head ADvector$$
00198 The type $italic ADvector$$ must be a 
00199 $xref/SimpleVector//simple vector class/$$ with elements of type
00200 $syntax%AD<%Base%>%$$.
00201 The routine $xref/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 $xref/LuRatio.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 = ip.size();
00244         CPPAD_ASSERT_KNOWN(
00245                 jp.size() == n,
00246                 "Error in LuFactor: jp must have size equal to n"
00247         );
00248         CPPAD_ASSERT_KNOWN(
00249                 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