F6-1.FOR

c-----CSL/F6-1.FOR ELMINACION DE GAUSS
      DIMENSION A(10,11)
      PRINT *
      PRINT *, 'CSL/F6-1        ELMINACION DE GAUSS'
      DATA N/3/                                 !-- N es el orden de la matriz
      DATA (A(1,J),J=1,4)/ 0, -1, +2, 0/        !-- inicializa los elementos de matriz
      DATA (A(2,J),J=1,4)/-2, +2, -1, 0/        !-- inicializa los elementos de la matriz
      DATA (A(3,J),J=1,4)/-2, +4, +3, 1/        !-- inicializa los elementos de la matriz
      PRINT *
      PRINT *,'MATRIZ AUMENTADA'
      PRINT *
      DO I=1,N
        PRINT 61, (A(I,J),J=1,4)
61      FORMAT(1X, 1P6E12.4)
      END DO
      PRINT *
      CALL GAUSS(N,A)
65    PRINT *
68    PRINT *,'SOLUCION'
69    PRINT *,'---------------------------------------'
      PRINT *,'          I        X(I)'
70    PRINT *,'---------------------------------------'
    DO I=1,N
72        FORMAT(5X, I5, 1PE16.6)
          PRINT 72,I,A(I,N+1)
        END DO
75    PRINT *,'---------------------------------------'
80    PRINT *
      STOP
      END
C***************************
      SUBROUTINE GAUSS(N,A)                     ! Eliminación de Gauss
      INTEGER PV
      DIMENSION A(10,11)
      EPS=1.0
10    IF (1.0+EPS.GT.1.0) THEN
        EPS=EPS/2.0
        GOTO 10
      END IF
      EPS=EPS*2.
      PRINT *,'                 EPSILON DE LA MAQUINA=',EPS
      EPS2=EPS*2                !
1005  DET = 1                   ! Inicialización del determinante
      DO 1010 I=1,N-1
        PV=I
        DO J=I+1,N
          IF (ABS(A(PV,I)) .LT. ABS(A(J,I))) PV = J
        END DO
        IF (PV.EQ.I) GOTO 1050
        DO JC=1,N+1
          TM=A(I,JC)
          A(I,JC)=A(PV,JC)
          A(PV,JC)=TM
        END DO
1045    DET=-1*DET         ! Cada vez que se realice un pivoteo, cambia el signo DET
1050    IF (A(I,I).EQ.0) GOTO 1200 ! Una matriz singular si A(I,I) = 0
        DO JR=I+1, N                 ! Eliminación por debajo de la diagonal
          IF (A(JR,I).NE.0) THEN
            R=A(JR,I)/A(I,I)
            DO KC=I+1,N+1
              TEMP=A(JR,KC)
              A(JR,KC)=A(JR,KC)-R*A(I,KC)
              IF (ABS(A(JR,KC)).LT.EPS2*TEMP) A(JR,KC)=0.0
c               Si el resultado de la resta es menor que
c               el doble del épsilon de la máquina por el valor
c               original, se cambia su valor a cero.
            END DO
          END IF
1060    END DO
1010  CONTINUE
      DO I=1,N
        DET=DET*A(I,I)          ! Se calcula el determinante
      END DO
      PRINT *
      PRINT *,'DETERMINANTE = ', DET
      PRINT *
      IF (A(N,N).EQ.0) GOTO 1200
      A(N,N+1)=A(N,N+1)/A(N,N)
      DO NV=N-1,1,-1
        VA=A(NV,N+1)
        DO K=NV+1,N
          VA=VA-A(NV,K)*A(K,N+1)
        END DO
        A(NV,N+1)=VA/A(NV,NV)
      END DO
      RETURN
1200  PRINT *,'LA MATRIZ ES SINGUNLAR'
      STOP
      END