MODULE Prototypes
   IMPLICIT NONE
!
   INTERFACE
   !
   	  FUNCTION DASUM( N, X, INCX ) RESULT(S)
   	  ! Sum of the absolute values.
   	  	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)
   	  ! DNRM2 returns the euclidean norm of a vector
   	  	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)
	    ! DGECON estimates the reciprocal of the condition number of a general
        ! real matrix A, in either the 1-norm or the infinity-norm, using
        ! the LU factorization computed by DGETRF.
	 	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 ) 	
	    ! DGETRF computes an LU factorization of a general M-by-N matrix A
        ! using partial pivoting with row interchanges.
        ! The factorization has the form
        !    A = P * L * U
	 	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
	  
	  FUNCTION DLANGE(NORM, M, N, A, LDA, WORK) RESULT(VALUE)	
	  ! DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or 
	  ! the largest absolute value of any element of a general rectangular matrix.
	  ! DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm'
	     IMPLICIT NONE
	     CHARACTER, INTENT(IN) :: NORM
	     INTEGER, INTENT(IN) :: M, N, LDA
	     REAL(8), DIMENSION(LDA,*), INTENT(IN) :: A
	     REAL(8), DIMENSION(*), INTENT(OUT) :: WORK
	     REAL(8) :: VALUE
	  END FUNCTION DLANGE
	  !
   END INTERFACE
!
CONTAINS
	!
      FUNCTION MatrixNorm( Norm, N, A, LDA, Info) RESULT(anorm)
      	! Computes 1, Infinite or 2 (Frobenius) 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
	        !! call DASUM to compute the absolute values of A[][]
	        !! use TV for storing the results:
				.	.	.
			anorm = MAXVAL(TV(1:N))
			DEALLOCATE(TV)
		CASE("2")
			!! call DNRM2 to compute this norm:
				.	.	.
		CASE("i","I")
			! Row wise
			ALLOCATE( TV(LDA), STAT=rc)
			IF ( rc /= 0 ) THEN
				INFO = -9
				RETURN
			END IF
	        !! call DASUM to compute the absolute values of A[][]
	        !! use TV for storing the results:
				.	.	.
			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, r1cond, rIcond, rmin, rmax, v
   REAL(8) :: da1norm, da2norm, daInorm

   INTEGER :: N, Np1, Nm1, LA, np, pos
   INTEGER :: i, j, rc, c0, c1, cycles_per_second, max_cycles
   
   ! 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..." 
   CALL SYSTEM_CLOCK(COUNT=c0, COUNT_RATE=cycles_per_second, &
                    COUNT_MAX=max_cycles)

   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
   CALL SYSTEM_CLOCK(COUNT=c1)
   WRITE(*,*) "Seconds for computing matrix norms: ",DBLE(c1-c0)/DBLE(cycles_per_second)
   
   WRITE(*,*) "Computing matrix norm with DLANGE ..." 
   CALL SYSTEM_CLOCK(COUNT=c0)
   da1norm = dlange("1",LA,LA,A,LA,WORK)
   da2norm = dlange("F",LA,LA,A,LA,WORK)
   daInorm = dlange("I",LA,LA,A,LA,WORK)
   CALL SYSTEM_CLOCK(COUNT=c1)
   WRITE(*,*) "Seconds for computing matrix norms with LAPACK: ", &
      & DBLE(c1-c0)/DBLE(cycles_per_second)

   
   WRITE(*,*) "Computing matrix factorization..."   
   CALL SYSTEM_CLOCK(COUNT=c0)
   Ipiv = 0.0D0
   !! Call DGETRF to factorize matrix A:
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGETRF INFO = ",RC
   END IF
   CALL SYSTEM_CLOCK(COUNT=c1)
   WRITE(*,*) "Seconds for matrix factorization: ", &
      & DBLE(c1-c0)/DBLE(cycles_per_second)

   
   IF ( LA <= 5 ) THEN
      print*," Factorized A = "
      DO i = 1, LA
         PRINT*,A(i,:)
      END DO
      print*,""
   END IF
   
   WRITE(*,*) "Computing Condition number:"  
   CALL SYSTEM_CLOCK(COUNT=c1)
   work = 0.0D0; iwork = 0.0D0
   !! call DGECON to compute r1cond with a1norm:
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGECON INFO = ",RC
   END IF
   work = 0.0D0; iwork = 0.0D0
   !! call DGECON to compute rIcond with aInorm:
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGECON INFO = ",RC
   END IF
   CALL SYSTEM_CLOCK(COUNT=c1)
   
   WRITE(*,*) " A1norm, Da1norm, Rcond = ", a1norm, da1norm, r1cond 
   WRITE(*,*) " A2norm, Da2norm        = ", a2norm, da2norm 
   WRITE(*,*) " AInorm, DaInorm, Rcond = ", aInorm, daInorm, rIcond  

   WRITE(*,*) "Seconds for calculating C.N.: ", &
      & DBLE(c1-c0)/DBLE(cycles_per_second)
   
   STOP
END PROGRAM CondNumber