anchor.diff:      INTEGER(KIND=4)            :: SH(2),J
   !
       SH = SHAPE(M)
  -    DO 110 J = 1, SH(1)
  -    110 IF (J.NE.COL) CALL MGEROW(M,J,COL,1)
  -    999 RETURN
  +    DO J = 1, SH(1)
  +      IF (J.NE.COL) CALL MGEROW(M,J,COL,1)
  +    ENDDO
     END SUBROUTINE MGJECOL
   !------------------------------------------------------------------------------+
   ! MGECOL: Matrix Gauss Eliminate COLumn
  
positive.code:       SH = SHAPE(M)
        IF (ABS(M(ROW,COL)).GE.TOL) THEN
          RATIO = M(ROW,COL) / M(COL,COL)
          DO 110 I = SCOL, SH(2)
          IF (I.EQ.COL) THEN
  !
  !  By definition this location must be zero.  Instead of doing the
  !  calculation, just set the value to zero.
  !
            M(ROW,I) = 0D0
          ELSE
            M(ROW,I) = M(ROW,I) - RATIO * M(COL,I)
          ENDIF
    110  CONTINUE
        ELSE
          M(ROW,COL) = 0D0
        END IF
    999 RETURN
    END SUBROUTINE MGEROW
  !------------------------------------------------------------------------------+
  ! MBSUB: Matrix Back SUBstitution
  
negative.code: !    M : Matrix
  !------------------------------------------------------------------------------+
    SUBROUTINE MBSUB(M)
      IMPLICIT NONE
      REAL(KIND=8),INTENT(INOUT) :: M(:,:)
      INTEGER I,J,K,SH(2)
      DOUBLE PRECISION SUM
      SH = SHAPE(M)
      DO I = SH(1)+1, SH(2)
        DO J = SH(1), 1, -1
          SUM = 0D0
          DO K = SH(1), J+1, -1
            SUM = SUM + M(J,K) * M(K,I)
          ENDDO
          M(J,I) = (M(J,I) - SUM) / M(J,J)
        ENDDO
      ENDDO
    END SUBROUTINE MBSUB
  END MODULE LINALG
  !------------------------------------------------------------------------------+
  