solve.cpp

Go to the documentation of this file.

#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#include <cmath>
#include <GfxTL/MathHelper.h>
 
int dmat_solve(int n, int rhs_num, double a[])
//******************************************************************************
//
//  Purpose:
//
//    DMAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system.
//
//  Discussion:
//
//    The doubly dimensioned array A is treated as a one dimensional vector,
//    stored by COLUMNS.  Entry A(I,J) is stored as A[I+J*N]
//
//  Modified:
//
//    29 August 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//
//    Input, int RHS_NUM, the number of right hand sides.  RHS_NUM
//    must be at least 0.
//
//    Input/output, double A[N*(N+RHS_NUM)], contains in rows and columns 1
//    to N the coefficient matrix, and in columns N+1 through
//    N+RHS_NUM, the right hand sides.  On output, the coefficient matrix
//    area has been destroyed, while the right hand sides have
//    been overwritten with the corresponding solutions.
//
//    Output, int DMAT_SOLVE, singularity flag.
//    0, the matrix was not singular, the solutions were computed;
//    J, factorization failed on step J, and the solutions could not
//    be computed.
//
{
    double apivot;
    double factor;
    int i;
    int ipivot;
    int j;
    int k;
    double temp;
 
    for (j = 0; j < n; j++)
    {
        //
        //  Choose a pivot row.
        //
        ipivot = j;
        apivot = a[j + j * n];
 
        for (i = j; i < n; i++)
        {
            if (fabs(apivot) < fabs(a[i + j * n]))
            {
                apivot = a[i + j * n];
                ipivot = i;
            }
        }
 
        if (apivot == 0.0)
        {
            return j;
        }
        //
        //  Interchange.
        //
        for (i = 0; i < n + rhs_num; i++)
        {
            temp          = a[ipivot + i * n];
            a[ipivot + i * n] = a[j + i * n];
            a[j + i * n]      = temp;
        }
        //
        //  A(J,J) becomes 1.
        //
        a[j + j * n] = 1.0;
        for (k = j; k < n + rhs_num; k++)
        {
            a[j + k * n] = a[j + k * n] / apivot;
        }
        //
        //  A(I,J) becomes 0.
        //
        for (i = 0; i < n; i++)
        {
            if (i != j)
            {
                factor = a[i + j * n];
                a[i + j * n] = 0.0;
                for (k = j; k < n + rhs_num; k++)
                {
                    a[i + k * n] = a[i + k * n] - factor * a[j + k * n];
                }
            }
 
        }
 
    }
 
    return 0;
}
//******************************************************************************