anchor.diff:  !
       SH = SHAPE(M)
       IF (ABS(M(ROW,ROW)).LE.TOL) THEN
  -    DO 110 J = ROW+1, SH(1)
  -    110 IF (ABS(M(J,ROW)).GE.TOL) EXIT
  -    DO 120 K = ROW,SH(2)
  -    TMP = M(ROW,K)
  -    M(ROW,K) = M(J,K)
  -    120 M(J,K) = TMP
  +    DO J = ROW+1, SH(1)
  +      IF (ABS(M(J,ROW)).GE.TOL) EXIT
  +    ENDDO
  +    DO K = ROW,SH(2)
  +      TMP = M(ROW,K)
  +      M(ROW,K) = M(J,K)
  +      M(J,K) = TMP
  +    ENDDO
       END IF
  -    999 RETURN
     END SUBROUTINE MPPVT
   !------------------------------------------------------------------------------+
   ! MGJECOL: Matrix Gauss-Jordan Eliminate COLumn
  
positive.code:     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
    END SUBROUTINE MGJECOL
  !------------------------------------------------------------------------------+
  ! MGECOL: Matrix Gauss Eliminate COLumn
  
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
  !------------------------------------------------------------------------------+
  