MODULE Prototypes
   IMPLICIT NONE
!
   INTERFACE
   !
   	  FUNCTION DASUM( N, X, INCX ) RESULT(S)
   	  	INTEGER, INTENT(IN) :: N, INCX
   	  	REAL(8), DIMENSION(N), INTENT(IN) :: X
   	  	REAL(8) :: S   	  
   	  END FUNCTION DASUM
   	  
   	  FUNCTION DNRM2( N, X, INCX ) RESULT(DN)
   	  	INTEGER, INTENT(IN) :: N, INCX
   	  	REAL(8), DIMENSION(N), INTENT(IN) :: X
   	  	REAL(8) :: DN  	  
   	  END FUNCTION DNRM2
   	  
      SUBROUTINE DGECON( &
    &		NORM,   &
    &   	N,   &
	&	  	A,   &
	&	  	LDA,   &
	&	  	ANORM,   &
	&	 	RCOND,   &
	&		WORK,   &
	&		IWORK,   &
	&	 	INFO   &
	&	) 	
	 	CHARACTER(1), INTENT(IN) ::	NORM
		INTEGER, INTENT(IN) :: LDA, N
		REAL(8), DIMENSION( LDA, N ), INTENT(IN) ::	A
		REAL(8), INTENT(IN) :: ANORM
		REAL(8), INTENT(OUT) :: RCOND
		REAL(8), DIMENSION( 4*N ), INTENT(IN OUT) :: WORK
		INTEGER, DIMENSION( N ), INTENT(IN OUT) :: IWORK
		INTEGER, INTENT(OUT) :: INFO 
	  END SUBROUTINE DGECON
	  
	  subroutine DGETRF( &
	&		M,   &
	&	  	N,   &
	&	  	A,   &
	&	  	LDA,   &
	&	  	IPIV,   &
	&	  	INFO    &
	&	) 	
	 	INTEGER, INTENT(IN) :: M, N, LDA
		REAL(8), dimension( LDA, N ), INTENT(IN OUT) ::	A
		INTEGER, dimension( * ), INTENT(OUT) ::IPIV
		INTEGER, INTENT(OUT) ::	INFO 
	  END SUBROUTINE DGETRF
	  
	  SUBROUTINE DGETRS( &
	&    	TRANS,   &
	&	  	N,   &
	&	  	NRHS,   &
	&	  	A,   &
	&	  	LDA,   &
	&	  	IPIV,   &
	&	  	B,   &
	&	  	LDB,   &
	&	  	INFO   &
	&	) 	
		CHARACTER(1), INTENT(IN) :: TRANS
		INTEGER, INTENT(IN) :: N, LDA, LDB, NRHS
		REAL(8), dimension( LDA,N), INTENT(IN) :: A
		INTEGER, dimension(N), INTENT(IN) :: IPIV
		REAL(8), dimension( LDB,NRHS), INTENT(IN OUT) :: B
		INTEGER, INTENT(OUT) :: INFO
	  END SUBROUTINE DGETRS
	  !
   END INTERFACE
!
CONTAINS
	!
      FUNCTION MatrixNorm( Norm, N, A, LDA, Info) RESULT(anorm)
      	! Computes 1, Infinite or 2 norm of a matrix A(LDA,N)
      	! Info < 0 if input is not consistent or troubles arise
	 	CHARACTER(1), INTENT(IN) ::	Norm
		INTEGER, INTENT(IN) :: LDA, N
		INTEGER, INTENT(OUT) :: Info
		REAL(8), DIMENSION( LDA, N ), INTENT(IN) ::	A
		REAL(8) :: ANORM
		!
		REAL(8), ALLOCATABLE, DIMENSION(:) :: TV
		INTEGER :: j, rc
		!
		INFO = 0
		anorm = 0.0D0
		SELECT CASE( Norm(1:1) )
		CASE("1")
			! Column wise
			ALLOCATE( TV(N), STAT=rc)
			IF ( rc /= 0 ) THEN
				INFO = -9
				RETURN
			END IF
			DO j = 1, N
				TV(j) = DASUM(LDA,A(1,j),1)
			END DO
			anorm = MAXVAL(TV(1:N))
			DEALLOCATE(TV)
		CASE("2")
			anorm = DNRM2(LDA*N,A,1)
		CASE("i","I")
			! Row wise
			ALLOCATE( TV(LDA), STAT=rc)
			IF ( rc /= 0 ) THEN
				INFO = -9
				RETURN
			END IF
			DO j = 1, LDA
				TV(j) = DASUM(N,A(j,1),LDA)
			END DO
			anorm = MAXVAL(TV(1:LDA))
			DEALLOCATE(TV)
		CASE DEFAULT
			INFO = -1
		END SELECT
		RETURN
	  END FUNCTION MatrixNorm

END MODULE Prototypes



PROGRAM CondNumber
   USE Prototypes
   IMPLICIT NONE
   !
   ! Computes the Condition Number of a matrix A(LA,N) 
   !
   REAL(8), ALLOCATABLE, DIMENSION(:,:) :: A
   REAL(8), ALLOCATABLE, DIMENSION(:) :: work
   INTEGER, ALLOCATABLE, DIMENSION(:) :: Ipiv, iwork
   REAL(8) :: x, y, xl, xr, yl, yu, dx, dy, Eps, Solution
   REAL(8) :: a1norm, a2norm, aInorm, rcond, rmin, rmax, v

   INTEGER :: N, Np1, Nm1, LA, np, pos
   INTEGER :: i, j, rc
   
   ! Input reading in the format:
   ! LA N
   ! A(1,1)    A(1,2)   A(1,3)   ... A(1,N)
   !		.	.	.
   ! A(LA,1)  A(LA,2) A(LA,3) ... A(LA,N)
   !
   PRINT*," LA, N ?"   
   READ(*,*) LA, N
   IF ( LA < 0 .OR. LA > 1000 ) LA = 1000
   IF ( N < 0 .OR. N > 1000 ) N = 1000
   PRINT*," LA, N = ",LA,N
!   
   ALLOCATE(A(LA,N),STAT=rc)
   IF ( rc /= 0 ) THEN
   	   WRITE(*,*) "Error allocating A matrix"
   	   STOP
   END IF
!    
   ALLOCATE(Ipiv(LA), work(4*LA), iwork(LA), STAT=rc)
   IF ( rc /= 0 ) THEN
   	   WRITE(*,*) "Error allocating Ipiv vector"
   	   STOP
   END IF
!
!  Read A matrix 
!
   A = 0.0D0
   DO i = 1, LA
      READ(*,*) A(i,:)
   END DO
   
   IF ( LA <= 5 ) THEN
      print*," A = "
      DO i = 1, LA
         PRINT*,A(i,:)
      END DO
      print*,""
   END IF

   WRITE(*,*) "Computing matrix norm..." 
   a1norm = MatrixNorm("1",LA,A,LA,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After MatrixNorm('1') INFO = ",RC
   END IF
   a2norm = MatrixNorm("2",LA,A,LA,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After MatrixNorm('2') INFO = ",RC
   END IF
   aInorm = MatrixNorm("I",LA,A,LA,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After MatrixNorm('I') INFO = ",RC
   END IF
   
   WRITE(*,*) "Computing matrix factorization..."   
   Ipiv = 0.0D0
   CALL DGETRF( LA, LA, A, LA, Ipiv, rc )
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGETRF INFO = ",RC
   END IF
   
   IF ( LA <= 5 ) THEN
      print*," Factorized A = "
      DO i = 1, LA
         PRINT*,A(i,:)
      END DO
      print*,""
   END IF
   
   WRITE(*,*) "Computing Condition number:"  
   work = 0.0D0; iwork = 0.0D0
   CALL DGECON("1",LA,A,LA,a1norm,rcond,work,iwork,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGECON INFO = ",RC
   END IF
   WRITE(*,*) " A1norm, Rcond = ", a1norm, rcond 
   
   WRITE(*,*) " A2norm        = ", a2norm 
   
   work = 0.0D0; iwork = 0.0D0
   CALL DGECON("I",LA,A,LA,aInorm,rcond,work,iwork,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGECON INFO = ",RC
   END IF
   WRITE(*,*) " AInorm, Rcond = ", aInorm, rcond  
   
   STOP
END PROGRAM CondNumber