Função para Eliminação Gaussiana para Gretl

PROGRAM TESTA_ELIMINACAO_GAUSSIANA

REAL, DIMENSION(2,2) :: A

REAL, DIMENSION(2) :: B,X

INTEGER, DIMENSION(2) :: PIV

INTEGER :: RSLT

REAL :: EPS

        PRINT *,'EPS=?'

    READ *,EPS

   

        A(:,1) = (/ EPS, 1. /)

    A(:,2) = (/ 1., 1. /)

    B = (/ 1., 2. /)

        CALL ELIMINACAO_GAUSSIANA(A,B,X,RSLT)

    PRINT *,'ELIMINACAO_GAUSSIANA: RSLT=',RSLT

    PRINT *,'ELIMINACAO_GAUSSIANA: X=',X

    IF (RSLT==0) THEN

        A(:,1) = (/ EPS, 1. /)

        A(:,2) = (/ 1., 1. /)

        B = (/ 1., 2. /)     

            PRINT *,'A*X-B=',MATMUL(A,X)-B

    END IF

        A(:,1) = (/ EPS, 1. /)

    A(:,2) = (/ 1., 1. /)

    B = (/ 1., 2. /)

    CALL ELIMINACAO_GAUSSIANA_PIVOTAMENTO(A,B,X,PIV,RSLT)

    PRINT *,'ELIMINACAO_GAUSSIANA_PIVOTAMENTO: RSLT=',RSLT

    PRINT *,'ELIMINACAO_GAUSSIANA_PIVOTAMENTO: X=',X

    PRINT *,'ELIMINACAO_GAUSSIANA_PIVOTAMENTO: PIV=',PIV

    IF (RSLT==0) THEN

        A(:,1) = (/ EPS, 1. /)

        A(:,2) = (/ 1., 1. /)

        B = (/ 1., 2. /)

            PRINT *,'A*X-B=',MATMUL(A,X)-B

    END IF

CONTAINS

SUBROUTINE ELIMINACAO_GAUSSIANA(A,B,X,INFO)

REAL, DIMENSION(:,:), INTENT(INOUT) :: A

REAL, DIMENSION(:), INTENT(INOUT) :: B

REAL, DIMENSION(:), INTENT(OUT) :: X

INTEGER, INTENT(OUT) :: INFO

INTEGER :: I,J,K,N

    N = SIZE(A,DIM=1)

    IF ((N==SIZE(A,DIM=2)) .AND. (N==SIZE(B)) .AND. (N==SIZE(X))) THEN

        INFO = 0 ! INICIALIZA INFO COM O VALOR QUE INDICAR� SUCESSO NA EXECU��O DA SUBROTINA     

        DO K=1,N-1

            IF (A(K,K)==0) THEN

                INFO = -K

                RETURN

            END IF

            DO I=K+1,N

                Z = A(I,K) / A(K,K)

                A(I,K) = 0.0   ! APENAS PARA FINS ILUSTRATIVOS, NA PR�TICA N�O � NECESS�RIO

                DO J = K+1,N

                    A(I,J) = A(I,J)-Z*A(K,J)

                END DO

                B(I) = B(I)-Z*B(K)

            END DO

        END DO

        IF (A(N,N)==0) THEN

                INFO = -N

                RETURN

        END IF

        CALL RETRO_SUBSTITUICAO(A,B,X)

    ELSE

        INFO = 1  ! DIMENS�ES INCOMPAT�VEIS ENTRE A MATRIZ E/OU VETORES

    END IF

END SUBROUTINE ELIMINACAO_GAUSSIANA

SUBROUTINE RETRO_SUBSTITUICAO(A,B,X)

REAL, DIMENSION(:,:), INTENT(IN) :: A

REAL, DIMENSION(:), INTENT(IN) :: B

REAL, DIMENSION(:), INTENT(OUT) :: X

INTEGER :: I,N

    N = SIZE(A,DIM=1)

    IF ((N==SIZE(A,DIM=2)) .AND. (N==SIZE(B)) .AND. (N==SIZE(X))) THEN

        X = B

        DO I=N,2,-1

            X(I) = X(I)/A(I,I)  ! ASSUME QUE A(I,I)/=0.0

            X(1:I-1) = X(1:I-1)-X(I)*A(1:I-1,I)

        END DO

        X(1) = X(1)/A(1,1) ! ASSUME QUE A(1,1)/=0.0

    END IF

END SUBROUTINE RETRO_SUBSTITUICAO

SUBROUTINE ELIMINACAO_GAUSSIANA_PIVOTAMENTO(A,B,X,P,INFO)

REAL, DIMENSION(:,:), INTENT(INOUT) :: A

REAL, DIMENSION(:), INTENT(INOUT) :: B

REAL, DIMENSION(:), INTENT(OUT) :: X

INTEGER, DIMENSION(:), INTENT(OUT) :: P

INTEGER, INTENT(OUT) :: INFO

REAL, DIMENSION(:), ALLOCATABLE :: S

INTEGER :: I,J,K,N,Q

    N = SIZE(A,DIM=1)

    IF ((N==SIZE(A,DIM=2)) .AND. (N==SIZE(B)) .AND. (N==SIZE(X)) .AND. (N==SIZE(P))) THEN

        INFO = 0 ! INICIALIZA INFO COM O VALOR QUE INDICAR� SUCESSO NA EXECU��O DA SUBROTINA

        ALLOCATE(S(N))

        DO I=1,N

            P(I) = I ; S(I) = MAXVAL(ABS(A(I,:)))

        END DO

        DO K=1,N-1

            !  LOCALIZA O PIV�, INDICADO PELO VALOR DE J

            J = K

            DO I=K+1,N

                IF (ABS(A(P(I),K))/S(P(I)) >= ABS(A(P(J),K))/S(P(J))) THEN

                    J = I

                    EXIT

                END IF

            END DO

            ! FAZ A TROCA DA LINHA K COM A LINHA J, USANDO O VETOR P

            Q = P(K) ; P(K) = P(J) ; P(J) = Q

            IF (A(P(K),K)/=0.0) THEN

                DO I=K+1,N

                    Z = A(P(I),K)/A(P(K),K)

                    DO J=K+1,N

                        A(P(I),J) = A(P(I),J)-Z*A(P(K),J)

                    END DO

                    B(P(I)) = B(P(I))-Z*B(P(K))

                END DO

            ELSE

                INFO = -P(K)

            END IF

        END DO

        IF (A(P(N),N)==0.0) THEN

            INFO = -P(N)

            RETURN

        END IF

        ! RESOLVE O SISTEMA TRIANGULAR SUPERIOR

        X = B((/(P(I),I=1,N)/))

        DO I=N,2,-1

            X(I) = X(I)/A(P(I),I)

            X(1:I-1) = X(1:I-1)-X(I)*A(P(1:I-1),I)

        END DO

        X(1) = X(1)/A(P(1),1)

       

        ! DESALOCA S

        DEALLOCATE(S)

    ELSE

        INFO = 1  ! DIMENS�ES INCOMPAT�VEIS ENTRE A MATRIZ E/OU VETORES

    END IF

END SUBROUTINE ELIMINACAO_GAUSSIANA_PIVOTAMENTO

END PROGRAM TESTA_ELIMINACAO_GAUSSIANA