Função para Cálculo pelo Método SOR para Gretl

function [x,norma_r,k]=sor(A,b,x,omega,tol,k_max)

    // [sol,norma_r,k]=sor(A,b,x,omega,tol,k_max)

    // Calcula a solução para o sistema Ax=b usando o método SOR.

    // A solução é aceita sujeita a um resíduo menor do que "tol" desde que

    // não exceda a "k_max" iterações.

    // (C) Rudnei Dias da Cunha 2008

    mode(0)

    [n,m]=size(A);

    D=diag(diag(A));

    Q=inv(D);

    U=-inv(D)*(triu(A)-D);

    L=-inv(D)*(tril(A)-D);

    d=Q*b;

    for k=1:k_max

        // r=b-A*x

        r = b-A*x;

        norma_r=norm(r,'fro');

        mprintf('k=%g norma_r=%g\n',k,norma_r);

        if (norma_r<tol) break; end;

        // w=U*x+d

        w = U*x+d;

        // Resolve o sistema (I+A_L)x_gs=w

        x_gs = (eye(n,n)-L)\w;   

        // x=omega*x_gs+(1-omega)*x

        x = omega*x_gs+(1.0-omega)*x;

    end

endfunction