Função para Fatoração LUP para Gretl

PROGRAM TESTA_FATORACAO_LUP

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 FATORACAO_LUP(A,B,PIV,RSLT)

    PRINT *,'FATORACAO_LUP: RSLT=',RSLT

    PRINT *,'FATORACAO_LUP: PIV=',PIV

    IF (RSLT==0) THEN

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

            PRINT *,'RESOLVE_SISTEMA_LU: X=',X       

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

    END IF    

CONTAINS

SUBROUTINE FATORACAO_LUP(A,B,P,INFO)

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

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

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)

                    A(P(I),K) = Z

                    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

        ! DESALOCA S

        DEALLOCATE(S)

    ELSE

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

    END IF

END SUBROUTINE FATORACAO_LUP

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

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

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

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

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

INTEGER, INTENT(OUT) :: INFO

INTEGER :: I,J,N

REAL, DIMENSION(:), ALLOCATABLE :: Y

    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

        ALLOCATE(Y(N))

        ! RESOLVE LY=PB

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

        DO J=1,N-1

            Y(J+1:N) = Y(J+1:N)-A(P(J+1:N),J)*Y(J)

        END DO

        ! RESOLVE UX=Y

        X = Y

        DO J=N,2,-1

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

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

        END DO

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

        DEALLOCATE(Y)

    ELSE

        INFO = 1

    END IF

END SUBROUTINE RESOLVE_SISTEMA_LU

END PROGRAM TESTA_FATORACAO_LUP