gauss_jordan.h

Go to the documentation of this file.

#ifndef TBCI_SOLVER_GAUSS_JORDAN_H
#define TBCI_SOLVER_GAUSS_JORDAN_H  

#include "tbci/basics.h"
#include "tbci/matrix.h"

#ifdef PRAGMA_I
# pragma interface "gauss_jordan.h"
#endif

#ifndef GJ_MINVAL
# define GJ_MINVAL 1e-16
#endif

//include <huge_val.h>
#ifndef HVAL
# define HVAL 1e38
#endif

NAMESPACE_TBCI

INST(template <typename T> class Matrix friend char gaussj (Matrix<T>&, Matrix<T>&);)

template <typename T>
char gaussj(Matrix<T>& a, Matrix<T>& b)
{
        const int n = a.row, m = b.col;
        char err = 0;

        BVector<unsigned> indxc((unsigned)n);
        BVector<unsigned> indxr((unsigned)n);
        BVector<unsigned> ipiv ((unsigned)0, n);

        BCHK(a.col != b.row, MatErr, Wrong sizes in gaussj, b.row, 1);

        for (int i = 0; i < n; ++i) {
                double big = 0.0;
                int icol = 0, irow = 0; 
                for (int j = 0; j < n; ++j)
                        if (ipiv(j) != 1)
                                for (REGISTER int k = 0; k < n; ++k) {
                                        if (ipiv(k) == 0) {
                                                if (MATH__ fabs(a(j, k)) >= MATH__ fabs(big)) {
                                                        big = MATH__ fabs(a(j, k));
                                                        irow = j; icol = k;
                                                 }
                                        } else if (ipiv(k) > 1) { 
                                                STD__ cerr << "\n\n error: gaussj: singular Matrix-1 \n\n"; 
                                                err |= 1; 
                                        }
                                }
                ++(ipiv(icol));
                if (irow != icol) {
                        for (REGISTER int l = 0; l < n; ++l) 
                                SWAP(a(irow, l), a(icol, l));
                        for (REGISTER int z = 0; z < m; ++z) 
                                SWAP(b(irow, z), b(icol, z));
                }
                indxr(i) = irow;
                indxc(i) = icol;
                T pivinv;
                if (MATH__ fabs(a(icol, icol)) < GJ_MINVAL) {
                        STD__ cerr << "\n\n error: gaussj: singular Matrix-2 \n" << STD__ endl;
                        err |= 2; pivinv = (T)HVAL;
                } else 
                        pivinv = 1.0 / a(icol, icol);
                a(icol, icol) = (T)1.0;
                for (REGISTER int l = 0; l < n; ++l) 
                        a(icol, l) *= pivinv;
                for (REGISTER int z = 0; z < m; ++z) 
                        b(icol, z) *= pivinv;
                for (int ll = 0; ll < n; ++ll)
                        if (ll != icol) {
                                REGISTER T dum = a(ll, icol);
                                a(ll, icol) = (T)0.0;
                                for (REGISTER int v = 0; v < n; ++v)
                                        a(ll, v) -= a(icol, v) * dum;
                                for (REGISTER int w = 0; w < m; ++w)
                                        b(ll, w) -= b(icol, w) * dum;
                        }
        }
        for (int l = n-1; l >= 0; --l) {
                if (indxr(l) != indxc(l))
                        for(REGISTER int k = 0; k < n; ++k) 
                                SWAP(a(k, indxr(l)), a(k, indxc(l)));
        }
        return err;
}

NAMESPACE_END
#endif  /* TBCI_SOLVER_GAUSS_JORDAN_H */