main/sc/source/core/tool/interpr5.cxx
/* Matrix LUP decomposition according to the pseudocode of "Introduction to
 * Algorithms" by Cormen, Leiserson, Rivest, Stein.
 *
 * Added scaling for numeric stability.
 *
 * Given an n x n nonsingular matrix A, find a permutation matrix P, a unit
 * lower-triangular matrix L, and an upper-triangular matrix U such that PA=LU.
 * Compute L and U "in place" in the matrix A, the original content is
 * destroyed. Note that the diagonal elements of the U triangular matrix
 * replace the diagonal elements of the L-unit matrix (that are each ==1). The
 * permutation matrix P is an array, where P[i]=j means that the i-th row of P
 * contains a 1 in column j. Additionally keep track of the number of
 * permutations (row exchanges).
 *
 * Returns 0 if a singular matrix is encountered, else +1 if an even number of
 * permutations occurred, or -1 if odd, which is the sign of the determinant.
 * This may be used to calculate the determinant by multiplying the sign with
 * the product of the diagonal elements of the LU matrix.
 */
Matrix LUP decomposition according to the pseudocode of "Introduction to
Algorithms" by Cormen, Leiserson, Rivest, Stein.

Added scaling for numeric stability.

Given an n x n nonsingular matrix A, find a permutation matrix P, a unit
lower-triangular matrix L, and an upper-triangular matrix U such that PA=LU.
Compute L and U "in place" in the matrix A, the original content is
destroyed. Note that the diagonal elements of the U triangular matrix
replace the diagonal elements of the L-unit matrix (that are each ==1). The
permutation matrix P is an array, where P[i]=j means that the i-th row of P
contains a 1 in column j. Additionally keep track of the number of
permutations (row exchanges).

Returns 0 if a singular matrix is encountered, else +1 if an even number of
permutations occurred, or -1 if odd, which is the sign of the determinant.
This may be used to calculate the determinant by multiplying the sign with
the product of the diagonal elements of the LU matrix.
Matrix LUP decomposition according to the pseudocode of " Introduction to Algorithms " by Cormen Leiserson Rivest Stein . Added scaling for numeric stability . Given an n x n nonsingular matrix A find a permutation matrix P a unit lower - triangular matrix L and an upper - triangular matrix U such that PA = LU . Compute L and U " in place " in the matrix A the original content is destroyed . Note that the diagonal elements of the U triangular matrix replace the diagonal elements of the L - unit matrix ( that are each == 1 ) . The permutation matrix P is an array where P [ i ] = j means that the i - th row of P contains a 1 in column j . Additionally keep track of the number of permutations ( row exchanges ) . Returns 0 if a singular matrix is encountered else + 1 if an even number of permutations occurred or - 1 if odd which is the sign of the determinant . This may be used to calculate the determinant by multiplying the sign with the product of the diagonal elements of the LU matrix .
static int lcl_LUP_decompose( ScMatrix* mA, const SCSIZE n,
        ::std::vector< SCSIZE> & P )
{
    int nSign = 1;
    ::std::vector< double> aScale(n);
    for (SCSIZE i=0; i < n; ++i)
    {
        double fMax = 0.0;
        for (SCSIZE j=0; j < n; ++j)
        {
            double fTmp = fabs( mA->GetDouble( j, i));
            if (fMax < fTmp)
                fMax = fTmp;
        }
        if (fMax == 0.0)
            return 0;       
        aScale[i] = 1.0 / fMax;
    }
    for (SCSIZE i=0; i < n; ++i)
        P[i] = i;
    SCSIZE l = n - 1;
    for (SCSIZE k=0; k < l; ++k)
    {
        double fMax = 0.0;
        double fScale = aScale[k];
        SCSIZE kp = k;
        for (SCSIZE i = k; i < n; ++i)
        {
            double fTmp = fScale * fabs( mA->GetDouble( k, i));
            if (fMax < fTmp)
            {
                fMax = fTmp;
                kp = i;
            }
        }
        if (fMax == 0.0)
            return 0;       
        if (k != kp)
        {
            SCSIZE nTmp = P[k];
            P[k]        = P[kp];
            P[kp]       = nTmp;
            nSign       = -nSign;
            double fTmp = aScale[k];
            aScale[k]   = aScale[kp];
            aScale[kp]  = fTmp;
            for (SCSIZE i=0; i < n; ++i)
            {
                double fMatTmp = mA->GetDouble( i, k);
                mA->PutDouble( mA->GetDouble( i, kp), i, k);
                mA->PutDouble( fMatTmp, i, kp);
            }
        }
        for (SCSIZE i = k+1; i < n; ++i)
        {
            double fTmp = mA->GetDouble( k, i) / mA->GetDouble( k, k);
            mA->PutDouble( fTmp, k, i);
            for (SCSIZE j = k+1; j < n; ++j)
                mA->PutDouble( mA->GetDouble( j, i) - fTmp * mA->GetDouble( j,
                            k), j, i);
        }
    }
#if OSL_DEBUG_LEVEL > 1
    fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): LU");
    for (SCSIZE i=0; i < n; ++i)
    {
        for (SCSIZE j=0; j < n; ++j)
            fprintf( stderr, "%8.2g  ", mA->GetDouble( j, i));
        fprintf( stderr, "\n%s\n", "");
    }
    fprintf( stderr, "\n%s\n", "lcl_LUP_decompose(): P");
    for (SCSIZE j=0; j < n; ++j)
        fprintf( stderr, "%5u ", (unsigned)P[j]);
    fprintf( stderr, "\n%s\n", "");
#endif
    bool bSingular=false;
    for (SCSIZE i=0; i < n && !bSingular; i++)
        bSingular = (mA->GetDouble(i,i) == 0.0);
    if (bSingular)
        nSign = 0;
    return nSign;
}
static int lcl_LUP_decompose ( ScMatrix * mA , const SCSIZE n , :: std :: vector < SCSIZE > & P ) { int nSign = 1 ; :: std :: vector < double > aScale ( n ) ; for ( SCSIZE i = 0 ; i < n ; ++ i ) { double fMax = 0.0 ; for ( SCSIZE j = 0 ; j < n ; ++ j ) { double fTmp = fabs ( mA -> GetDouble ( j , i ) ) ; if ( fMax < fTmp ) fMax = fTmp ; } if ( fMax == 0.0 ) return 0 ; aScale [ i ] = 1.0 / fMax ; } for ( SCSIZE i = 0 ; i < n ; ++ i ) P [ i ] = i ; SCSIZE l = n - 1 ; for ( SCSIZE k = 0 ; k < l ; ++ k ) { double fMax = 0.0 ; double fScale = aScale [ k ] ; SCSIZE kp = k ; for ( SCSIZE i = k ; i < n ; ++ i ) { double fTmp = fScale * fabs ( mA -> GetDouble ( k , i ) ) ; if ( fMax < fTmp ) { fMax = fTmp ; kp = i ; } } if ( fMax == 0.0 ) return 0 ; if ( k != kp ) { SCSIZE nTmp = P [ k ] ; P [ k ] = P [ kp ] ; P [ kp ] = nTmp ; nSign = - nSign ; double fTmp = aScale [ k ] ; aScale [ k ] = aScale [ kp ] ; aScale [ kp ] = fTmp ; for ( SCSIZE i = 0 ; i < n ; ++ i ) { double fMatTmp = mA -> GetDouble ( i , k ) ; mA -> PutDouble ( mA -> GetDouble ( i , kp ) , i , k ) ; mA -> PutDouble ( fMatTmp , i , kp ) ; } } for ( SCSIZE i = k + 1 ; i < n ; ++ i ) { double fTmp = mA -> GetDouble ( k , i ) / mA -> GetDouble ( k , k ) ; mA -> PutDouble ( fTmp , k , i ) ; for ( SCSIZE j = k + 1 ; j < n ; ++ j ) mA -> PutDouble ( mA -> GetDouble ( j , i ) - fTmp * mA -> GetDouble ( j , k ) , j , i ) ; } } #if OSL_DEBUG_LEVEL > 1 fprintf ( stderr , " \n \n " , " " ) ; for ( SCSIZE i = 0 ; i < n ; ++ i ) { for ( SCSIZE j = 0 ; j < n ; ++ j ) fprintf ( stderr , " " , mA -> GetDouble ( j , i ) ) ; fprintf ( stderr , " \n \n " , " " ) ; } fprintf ( stderr , " \n \n " , " " ) ; for ( SCSIZE j = 0 ; j < n ; ++ j ) fprintf ( stderr , " " , ( unsigned ) P [ j ] ) ; fprintf ( stderr , " \n \n " , " " ) ; #endif bool bSingular = false ; for ( SCSIZE i = 0 ; i < n && ! bSingular ; i ++ ) bSingular = ( mA -> GetDouble ( i , i ) == 0.0 ) ; if ( bSingular ) nSign = 0 ; return nSign ; }
Matrix LUP decomposition according to the pseudocode of "Introduction to
 Algorithms" by Cormen, Leiserson, Rivest, Stein.
Matrix LUP decomposition according to the pseudocode of " Introduction to Algorithms " by Cormen Leiserson Rivest Stein .
// Find scale of each row. // singular matrix // Represent identity permutation, P[i]=i // "Recursion" on the diagonale. // Implicit pivoting. With the scale found for a row, compare values of // a column and pick largest. // singular matrix // Swap rows. The pivot element will be at mA[k,kp] (row,col notation) // permutations // scales // elements // Compute Schur complement.