main/sc/source/core/tool/interpr5.cxx
/* Solve a LUP decomposed equation Ax=b. LU is a combined matrix of L and U
 * triangulars and P the permutation vector as obtained from
 * lcl_LUP_decompose(). B is the right-hand side input vector, X is used to
 * return the solution vector.
 */
Solve a LUP decomposed equation Ax=b. LU is a combined matrix of L and U
triangulars and P the permutation vector as obtained from
lcl_LUP_decompose(). B is the right-hand side input vector, X is used to
return the solution vector.
Solve a LUP decomposed equation Ax = b . LU is a combined matrix of L and U triangulars and P the permutation vector as obtained from lcl_LUP_decompose () . B is the right - hand side input vector X is used to return the solution vector .
static void lcl_LUP_solve( const ScMatrix* mLU, const SCSIZE n,
        const ::std::vector< SCSIZE> & P, const ::std::vector< double> & B,
        ::std::vector< double> & X )
{
    SCSIZE nFirst = SCSIZE_MAX;
    for (SCSIZE i=0; i < n; ++i)
    {
        double fSum = B[P[i]];
        if (nFirst != SCSIZE_MAX)
        {
            for (SCSIZE j = nFirst; j < i; ++j)
                fSum -= mLU->GetDouble( j, i) * X[j];   
        }
        else if (fSum)
            nFirst = i;
        X[i] = fSum;                                    
    }
    for (SCSIZE i = n; i--; )
    {
        double fSum = X[i];                             
        for (SCSIZE j = i+1; j < n; ++j)
            fSum -= mLU->GetDouble( j, i) * X[j];       
        X[i] = fSum / mLU->GetDouble( i, i);            
    }
#if OSL_DEBUG_LEVEL >1
    fprintf( stderr, "\n%s\n", "lcl_LUP_solve():");
    for (SCSIZE i=0; i < n; ++i)
        fprintf( stderr, "%8.2g  ", X[i]);
    fprintf( stderr, "%s\n", "");
#endif
}
static void lcl_LUP_solve ( const ScMatrix * mLU , const SCSIZE n , const :: std :: vector < SCSIZE > & P , const :: std :: vector < double > & B , :: std :: vector < double > & X ) { SCSIZE nFirst = SCSIZE_MAX ; for ( SCSIZE i = 0 ; i < n ; ++ i ) { double fSum = B [ P [ i ] ] ; if ( nFirst != SCSIZE_MAX ) { for ( SCSIZE j = nFirst ; j < i ; ++ j ) fSum -= mLU -> GetDouble ( j , i ) * X [ j ] ; } else if ( fSum ) nFirst = i ; X [ i ] = fSum ; } for ( SCSIZE i = n ; i -- ; ) { double fSum = X [ i ] ; for ( SCSIZE j = i + 1 ; j < n ; ++ j ) fSum -= mLU -> GetDouble ( j , i ) * X [ j ] ; X [ i ] = fSum / mLU -> GetDouble ( i , i ) ; } #if OSL_DEBUG_LEVEL > 1 fprintf ( stderr , " \n \n " , " " ) ; for ( SCSIZE i = 0 ; i < n ; ++ i ) fprintf ( stderr , " " , X [ i ] ) ; fprintf ( stderr , " \n " , " " ) ; #endif }
Solve a LUP decomposed equation Ax=b.
Solve a LUP decomposed equation Ax = b .
// Ax=b => PAx=Pb, with decomposition LUx=Pb. // Define y=Ux and solve for y in Ly=Pb using forward substitution. // Matrix inversion comes with a lot of zeros in the B vectors, we // don't have to do all the computing with results multiplied by zero. // Until then, simply lookout for the position of the first nonzero // value. // X[j] === y[j] // X[i] === y[i] // Solve for x in Ux=y using back substitution. // X[i] === y[i] // X[j] === x[j] // X[i] === x[i]