# ifndef CPPAD_LU_RATIO_INCLUDED # define CPPAD_LU_RATIO_INCLUDED /* -------------------------------------------------------------------------- CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-07 Bradley M. Bell CppAD is distributed under multiple licenses. This distribution is under the terms of the Common Public License Version 1.0. A copy of this license is included in the COPYING file of this distribution. Please visit http://www.coin-or.org/CppAD/ for information on other licenses. -------------------------------------------------------------------------- */ /* $begin LuRatio$$ $spell cppad.hpp xk Cpp Lu bool const ip jp std ADvector $$ $index LuRatio$$ $index linear, Lu factor equation$$ $index equation, Lu factor$$ $index determinant, Lu factor$$ $index solve, Lu factor$$ $section LU Factorization of A Square Matrix and Stability Calculation$$ $head Syntax$$ $code# include $$ $pre $$ $syntax%%sign% = LuRatio(%ip%, %jp%, %LU%, %ratio%)%$$ $head Description$$ Computes an LU factorization of the matrix $italic A$$ where $italic A$$ is a square matrix. A measure of the numerical stability called $italic ratio$$ is calculated. This ratio is useful when the results of $code LuRatio$$ are used as part of an $xref/ADFun/$$ object. $head Include$$ This routine is designed to be used with AD objects and requires the $code cppad/cppad.hpp$$ file to be included. $head Matrix Storage$$ All matrices are stored in row major order. To be specific, if $latex Y$$ is a vector that contains a $latex p$$ by $latex q$$ matrix, the size of $latex Y$$ must be equal to $latex p * q $$ and for $latex i = 0 , \ldots , p-1$$, $latex j = 0 , \ldots , q-1$$, $latex \[ Y_{i,j} = Y[ i * q + j ] \] $$ $head sign$$ The return value $italic sign$$ has prototype $syntax% int %sign% %$$ If $italic A$$ is invertible, $italic sign$$ is plus or minus one and is the sign of the permutation corresponding to the row ordering $italic ip$$ and column ordering $italic jp$$. If $italic A$$ is not invertible, $italic sign$$ is zero. $head ip$$ The argument $italic ip$$ has prototype $syntax% %SizeVector% &%ip% %$$ (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). The size of $italic ip$$ is referred to as $italic n$$ in the specifications below. The input value of the elements of $italic ip$$ does not matter. The output value of the elements of $italic ip$$ determine the order of the rows in the permuted matrix. $head jp$$ The argument $italic jp$$ has prototype $syntax% %SizeVector% &%jp% %$$ (see description of $xref/LuFactor/SizeVector/SizeVector/$$ below). The size of $italic jp$$ must be equal to $italic n$$. The input value of the elements of $italic jp$$ does not matter. The output value of the elements of $italic jp$$ determine the order of the columns in the permuted matrix. $head LU$$ The argument $italic LU$$ has the prototype $syntax% %ADvector% &%LU% %$$ and the size of $italic LU$$ must equal $latex n * n$$ (see description of $xref/LuRatio/ADvector/ADvector/$$ below). $subhead A$$ We define $italic A$$ as the matrix corresponding to the input value of $italic LU$$. $subhead P$$ We define the permuted matrix $italic P$$ in terms of $italic A$$ by $syntax% %P%(%i%, %j%) = %A%[ %ip%[%i%] * %n% + %jp%[%j%] ] %$$ $subhead L$$ We define the lower triangular matrix $italic L$$ in terms of the output value of $italic LU$$. The matrix $italic L$$ is zero above the diagonal and the rest of the elements are defined by $syntax% %L%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] %$$ for $latex i = 0 , \ldots , n-1$$ and $latex j = 0 , \ldots , i$$. $subhead U$$ We define the upper triangular matrix $italic U$$ in terms of the output value of $italic LU$$. The matrix $italic U$$ is zero below the diagonal, one on the diagonal, and the rest of the elements are defined by $syntax% %U%(%i%, %j%) = %LU%[ %ip%[%i%] * %n% + %jp%[%j%] ] %$$ for $latex i = 0 , \ldots , n-2$$ and $latex j = i+1 , \ldots , n-1$$. $subhead Factor$$ If the return value $italic sign$$ is non-zero, $syntax% %L% * %U% = %P% %$$ If the return value of $italic sign$$ is zero, the contents of $italic L$$ and $italic U$$ are not defined. $subhead Determinant$$ $index determinant$$ If the return value $italic sign$$ is zero, the determinant of $italic A$$ is zero. If $italic sign$$ is non-zero, using the output value of $italic LU$$ the determinant of the matrix $italic A$$ is equal to $syntax% %sign% * %LU%[%ip%[0], %jp%[0]] * %...% * %LU%[%ip%[%n%-1], %jp%[%n%-1]] %$$ $head ratio$$ The argument $italic ratio$$ has prototype $syntax% AD<%Base%> &%ratio% %$$ On input, the value of $italic ratio$$ does not matter. On output it is a measure of how good the choice of pivots is. For $latex p = 0 , \ldots , n-1$$, the $th p$$ pivot element is the element of maximum absolute value of a $latex (n-p) \times (n-p)$$ sub-matrix. The ratio of each element of sub-matrix divided by the pivot element is computed. The return value of $italic ratio$$ is the maximum absolute value of such ratios over with respect to all elements and all the pivots. $subhead Purpose$$ Suppose that the execution of a call to $code LuRatio$$ is recorded in the $syntax%ADFun<%Base%>%$$ object $italic F$$. Then a call to $xref/Forward/$$ of the form $syntax% %F%.Forward(%k%, %xk%) %$$ with $italic k$$ equal to zero will revaluate this Lu factorization with the same pivots and a new value for $italic A$$. In this case, the resulting $italic ratio$$ may not be one. If $italic ratio$$ is too large (the meaning of too large is up to you), the current pivots do not yield a stable LU factorization of $italic A$$. A better choice for the pivots (for this value of $italic A$$) will be made if you recreate the $code ADFun$$ object starting with the $xref/Independent/$$ variable values that correspond to the vector $italic xk$$. $head SizeVector$$ The type $italic SizeVector$$ must be a $xref/SimpleVector/$$ class with $xref/SimpleVector/Elements of Specified Type/elements of type size_t/$$. The routine $xref/CheckSimpleVector/$$ will generate an error message if this is not the case. $head Base$$ This type is determined by the argument $italic ratio$$ which has type $syntax%AD<%Base%>%$$. The type $italic Base$$ must be $code float$$, $code double$$, or any type in the $cref/AD levels above/glossary/AD Levels Above Base/$$ above $code float$$ or $code double$$. $head ADvector$$ The type $italic ADvector$$ must be a $xref/SimpleVector//simple vector class/$$ with elements of type $syntax%AD<%Base%>%$$. The routine $xref/CheckSimpleVector/$$ will generate an error message if this is not the case. $head Example$$ $children% example/lu_ratio.cpp %$$ The file $xref/LuRatio.cpp/$$ contains an example and test of using $code LuRatio$$. It returns true if it succeeds and false otherwise. $end -------------------------------------------------------------------------- */ namespace CppAD { // BEGIN CppAD namespace // Lines different from the code in cppad/lu_factor.hpp end with // template // int LuRatio(SizeVector &ip, SizeVector &jp, ADvector &LU, AD &ratio) // { typedef ADvector FloatVector; // typedef AD Float; // // check numeric type specifications CheckNumericType(); // check simple vector class specifications CheckSimpleVector(); CheckSimpleVector(); size_t i, j; // some temporary indices const Float zero( 0 ); // the value zero as a Float object size_t imax; // row index of maximum element size_t jmax; // column indx of maximum element Float emax; // maximum absolute value size_t p; // count pivots int sign; // sign of the permutation Float etmp; // temporary element Float pivot; // pivot element // ------------------------------------------------------- size_t n = ip.size(); CPPAD_ASSERT_KNOWN( jp.size() == n, "Error in LuFactor: jp must have size equal to n" ); CPPAD_ASSERT_KNOWN( LU.size() == n * n, "Error in LuFactor: LU must have size equal to n * m" ); // ------------------------------------------------------- // initialize row and column order in matrix not yet pivoted for(i = 0; i < n; i++) { ip[i] = i; jp[i] = i; } // initialize the sign of the permutation sign = 1; // initialize the ratio // ratio = Float(1); // // --------------------------------------------------------- // Reduce the matrix P to L * U using n pivots for(p = 0; p < n; p++) { // determine row and column corresponding to element of // maximum absolute value in remaining part of P imax = jmax = n; emax = zero; for(i = p; i < n; i++) { for(j = p; j < n; j++) { CPPAD_ASSERT_UNKNOWN( (ip[i] < n) & (jp[j] < n) ); etmp = LU[ ip[i] * n + jp[j] ]; // check if maximum absolute value so far if( AbsGeq (etmp, emax) ) { imax = i; jmax = j; emax = etmp; } } } for(i = p; i < n; i++) // { for(j = p; j < n; j++) // { etmp = abs(LU[ ip[i] * n + jp[j] ] / emax); // ratio = // CondExpGt(etmp, ratio, etmp, ratio); // } // } // CPPAD_ASSERT_KNOWN( (imax < n) & (jmax < n) , "AbsGeq must return true when second argument is zero" ); if( imax != p ) { // switch rows so max absolute element is in row p i = ip[p]; ip[p] = ip[imax]; ip[imax] = i; sign = -sign; } if( jmax != p ) { // switch columns so max absolute element is in column p j = jp[p]; jp[p] = jp[jmax]; jp[jmax] = j; sign = -sign; } // pivot using the max absolute element pivot = LU[ ip[p] * n + jp[p] ]; // check for determinant equal to zero if( pivot == zero ) { // abort the mission return 0; } // Reduce U by the elementary transformations that maps // LU( ip[p], jp[p] ) to one. Only need transform elements // above the diagonal in U and LU( ip[p] , jp[p] ) is // corresponding value below diagonal in L. for(j = p+1; j < n; j++) LU[ ip[p] * n + jp[j] ] /= pivot; // Reduce U by the elementary transformations that maps // LU( ip[i], jp[p] ) to zero. Only need transform elements // above the diagonal in U and LU( ip[i], jp[p] ) is // corresponding value below diagonal in L. for(i = p+1; i < n; i++ ) { etmp = LU[ ip[i] * n + jp[p] ]; for(j = p+1; j < n; j++) { LU[ ip[i] * n + jp[j] ] -= etmp * LU[ ip[p] * n + jp[j] ]; } } } return sign; } } // END CppAD namespace # endif