// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#ifndef MGCLINEARSYSTEM_H
#define MGCLINEARSYSTEM_H

//#include "MgcMath.h"


class MgcLinearSystem
{
public:
    MgcLinearSystem () { /**/ }

    // 2x2 and 3x3 systems (avoids overhead of Gaussian elimination)
    static float& Tolerance ();
    static bool Solve2 (const float aafA[2][2], const float afB[2],
        float afX[2]);
    static bool Solve3 (const float aafA[3][3], const float afB[3],
        float afX[3]);

    // convenience for allocation, memory is zeroed out
    static float** NewMatrix (int iSize);
    static void DeleteMatrix (int iSize, float** aafA);
    static float* NewVector (int iSize);
    static void DeleteVector (int iSize, float* afB);

    static bool Inverse (int iSize, float** aafA);
    // Input:
    //     A[iSize][iSize], entries are A[row][col]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     A[iSize][iSize], inverse matrix

    static bool Solve (int iSize, float** aafA, float* afB);
    // Input:
    //     A[iSize][iSize] coefficient matrix, entries are A[row][col]
    //     B[iSize] vector, entries are B[row]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     A[iSize][iSize] is inverse matrix
    //     B[iSize] is solution x to Ax = B

    static bool SolveTri (int iSize, float* afA, float* afB, float* afC,
        float* afR, float* afU);
    // Input:
    //     Matrix is tridiagonal.
    //     Lower diagonal A[iSize-1]
    //     Main  diagonal B[iSize]
    //     Upper diagonal C[iSize-1]
    //     Right-hand side R[iSize]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     U[iSize] is solution

    static bool SolveConstTri (int iSize, float fA, float fB, float fC,
        float* afR, float* afU);
    // Input:
    //     Matrix is tridiagonal.
    //     Lower diagonal is constant, A
    //     Main  diagonal is constant, B
    //     Upper diagonal is constant, C
    //     Right-hand side Rr[iSize]
    // Output:
    //     return value is TRUE if successful, FALSE if pivoting failed
    //     U[iSize] is solution

    static bool SolveSymmetric (int iSize, float** aafA, float* afB);
    // Input:
    //     A[iSize][iSize] symmetric matrix, entries are A[row][col]
    //     B[iSize] vector, entries are B[row]
    // Output:
    //     return value is TRUE if successful, FALSE if (nearly) singular
    //     decomposition A = L D L^t (diagonal terms of L are all 1)
    //         A[i][i] = entries of diagonal D
    //         A[i][j] for i > j = lower triangular part of L
    //     B[iSize] is solution to x to Ax = B

    static bool SymmetricInverse (int iSize, float** aafA,
        float** aafAInv);
    // Input:
    //     A[iSize][iSize], entries are A[row][col]
    // Output:
    //     return value is TRUE if successful, FALSE if algorithm failed
    //     AInv[iSize][iSize], inverse matrix

protected:
    // tolerance for 2x2 and 3x3 system solving
    static float ms_fTolerance;
};

//----------------------------------------------------------------------------
inline float& MgcLinearSystem::Tolerance ()
{
    return ms_fTolerance;
}
//----------------------------------------------------------------------------


#endif
