BibliotecaMatematica/gaussseidel.hpp at master · EngineeringLibrary/BibliotecaMatematica · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include "gaussseidel.h"

template<typename Type>
bool LinAlg::CritSassenfeld(const LinAlg::Matrix<Type> &A)
{
    LinAlg::Matrix<Type> beta(A.getNumberOfColumns(),1);
    Type BetaTemp = 0;
    for(unsigned i = 1; i <= A.getNumberOfColumns(); ++i )
    {
        for(unsigned j = 1; j <= i-1; ++j)
            beta(i,1) += Type (fabs(A(i,j)*beta(j,1)/(A(i,i))));

        for(unsigned j = i+1; j <= A.getNumberOfColumns(); ++j)
            beta(i,1) += Type (fabs(A(i,j)/(A(i,i))));

        if(BetaTemp < beta(i,1))
            BetaTemp = beta(i,1);
    }
    std::cout<< beta;
    if(BetaTemp >= 1)
    {
        return 0;
    }
    else
        return 1;
}

template<typename Type>
LinAlg::Matrix<Type> LinAlg::GaussSeidel(const LinAlg::Matrix<Type> &A, const LinAlg::Matrix<Type> &B,
                                         unsigned nIteracoes, Type tolerance)
{
    LinAlg::Matrix<Type> X = B, XA = B;

    for(unsigned i = 1; i <= nIteracoes; ++i)
    {
        Type MaxValue = 0;
        for(unsigned j = 1; j <= A.getNumberOfRows(); ++j)
        {
            X(j,1) = B(j,1)/A(j,j);

            for(unsigned k = 1; k <= A.getNumberOfColumns(); ++k)
                if(k != j)
                    X(j,1) -= A(j,k)*X(k,1)/A(j,j);

            if(fabs(Type(X(j,1)-XA(j,1))) > MaxValue)
                    MaxValue = fabs(Type(X(j,1)-XA(j,1)));
        }
        if(MaxValue < tolerance)
            break;
        XA = X;
    }
    return X;
}