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
	  
	  SUBROUTINE DGETRS(TRANS, N, NRHS, A, LDA, IPIV, B, LDB, INFO) 	
	    ! DGETRS solves a system of linear equations
        ! A * X = B  or  A' * X = B
        ! with a general N-by-N matrix A using the LU factorization computed
        ! by DGETRF.
		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 LinearEquation
   USE Prototypes
   IMPLICIT NONE
   !
   ! Let's consider the 2D-square 1 < X < 2, 1 < Y < 2
   ! Suppose we divide the square in N-1 x N-1 portions 
   ! Given a function F(x,y) suppose we know that in every point:
   !   d2F(x,x)/dx2 + d2F(x,y)/dy2 = x+y
   ! Let's define a system of equation so that we can compute the value
   ! of F(x,y) in every point of the 2D-square
   !
   REAL(8), ALLOCATABLE, DIMENSION(:,:) :: A, Points
   REAL(8), ALLOCATABLE, DIMENSION(:) :: S, B, work
   INTEGER, ALLOCATABLE, DIMENSION(:) :: Ipiv, iwork
   REAL(8) :: x, y, xl, xr, yl, yu, dx, dy, Eps
   REAL(8) :: anorm, rcond, rmin, rmax, v
!  Define 2D-square:
   REAL(8), PARAMETER :: x0 = 1.0D0, y0 = 1.0D0, x1 = 2.0D0, y1 = 2.0D0

   INTEGER :: N, Np1, Nm1, LA, np, pos
   INTEGER :: i, j, rc
   
   WRITE(*,*) " How many intervals?"
   READ(*,*) N
   IF ( N < 0 .OR. N > 1000 ) N = 1000
   
   Np1 = N+1  ! Grid points per side
   Nm1 = N-1  ! Grid points minus boundary
   LA = Nm1*Nm1 ! unknown points
   
   ALLOCATE(A(LA,LA),Points(LA,2),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
   
   ALLOCATE(B(LA),S(LA),STAT=rc)   ! Unknown and RHS vectors
   IF ( rc /= 0 ) THEN
   	   WRITE(*,*) "Error allocating B, X vectors"
   	   STOP
   END IF
   
   CALL GridDef(x0,x1,y0,y1,N,Points,rc)
   
   CALL EqsDef(x0,x1,y0,y1,N,LA,A,B,Points,rc)

   CALL PrintMatrix(LA,LA,A," A = ")
   
   WRITE(*,*) "Computing matrix norm..." 
   anorm = MatrixNorm("1",LA,A,LA,rc)
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After MatrixNorm('1') INFO = ",RC
   END IF
   
   WRITE(*,*) "Computing matrix factorization..."   
   Ipiv = 0.0D0
   ! call DGETRF to compute matrix factorization
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGETRF INFO = ",RC
   END IF
   
   CALL PrintMatrix(LA,LA,A," FactorizedA = ")
   
   WRITE(*,*) "Computing condition number..."  
   work = 0.0D0; iwork = 0.0D0
   ! call DGECON to compute condition number rcond using anorm
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGECON INFO = ",RC
   END IF
   WRITE(*,"(A,2(1X,F16.5))") " Anorm, Rcond = ", anorm, rcond  
   
   WRITE(*,*) "Computing the solution..." 
   ! call DGETRS to compute solution
		.	.	.
   IF ( rc /= 0 ) THEN
   	  WRITE(*,*) " After DGETRS INFO = ",RC
   END IF
   
   CALL PrintMatrix(1,LA,B," Solution = ")
   
   ! Check the solution
   rmin = 9999.0D0; rmax = 0.0D0
   DO np = 1, LA
      x = Points(np,1)
      y = Points(np,2)
      v = ABS(Solution(x,y) - B(np))
      IF ( rmin > v ) rmin = v 
      IF ( rmax < v ) rmax = v 
   END DO
   WRITE(*,"(A,2(1X,F16.5))") " Solution checking: rmin, rmax = ",rmin, rmax
   
   STOP
CONTAINS
   FUNCTION Solution(x,y) RESULT (f)
      IMPLICIT NONE
      REAL(8), INTENT(IN) :: x, y   ! Point coordinates
      REAL(8) :: f
      
      f = (x*x*x+y*y*y)/6.0D0
      
      RETURN
   END FUNCTION Solution
   	   
   SUBROUTINE GridDef(x0,x1,y0,y1,N,Pts,Info)
   		IMPLICIT NONE
   		REAL(8), INTENT(IN) :: x0,x1,y0,y1  ! Define grid
   		INTEGER, INTENT(IN) :: N
   		REAL(8), DIMENSION((N-1)*(N-1),2), INTENT(OUT) :: Pts   ! inner grid point Coordinates
   		INTEGER, INTENT(OUT) :: Info	! Info < 0 if errors occur 
!
		REAL(8) :: dx, dy
		INTEGER :: np, Nm1
		
		Nm1 = N-1
   		dx = (x1-x0)/DBLE(N); dy = (y1-y0)/DBLE(N)
   		np = 0
   		DO i = 1, Nm1
   			DO j = 1, Nm1
        		 np = np + 1
        		 IF ( np > Nm1*Nm1 ) THEN
        		 	 WRITE(*,*) " Error: NP = ",np," > N*N = ",Nm1*Nm1
        		 	 STOP
        		 END IF
        		 x = x0 + dx * DBLE(i)
        		 y = y0 + dy * DBLE(j)
        		 Pts(np,1) = x
        		 Pts(np,2) = y
        	END DO
        END DO

   		RETURN
   END SUBROUTINE GridDef
   
      
   SUBROUTINE EqsDef(x0,x1,y0,y1,N,LA,A,Rhs,Pts,Info)
        IMPLICIT NONE
   		REAL(8), INTENT(IN) :: x0,x1,y0,y1  ! Define grid
   		INTEGER, INTENT(IN) :: N, LA   ! Extent of grid side and number of grid points
   		REAL(8), DIMENSION(LA,2), INTENT(IN) :: Pts   ! inner grid point Coordinates
   		REAL(8), DIMENSION(LA), INTENT(OUT) :: Rhs   ! Linear equation RHS
   		REAL(8), DIMENSION(LA,LA), INTENT(OUT) :: A   ! Linear equation matrix
   		INTEGER, INTENT(OUT) :: Info	! Info < 0 if errors occur 
!
		REAL(8) :: dx, dy
		INTEGER :: np, Nm1
!
!  Define A matrix and RHS
!
   Nm1 = N-1
   dx = (x1-x0)/DBLE(N); dy = (y1-y0)/DBLE(N)

   A = 0.0D0; Rhs = 0.0D0
   DO np = 1, LA
      x = Points(np,1)
      y = Points(np,2)
      if ( LA < 6 ) WRITE(*,"(A,2(1X,F8.3))") "x, y = ",x,y
      A(np,np) = -4.0D0
      Rhs(np) = (x + y)*dx*dy
      ! define Eps function of grid dimensions 
      Eps = (dx+dy)/20.0D0
      ! where is P(x-dx,y) ? 
      IF ( ABS((x-dx)-x0) < Eps ) THEN
      	   Rhs(np) = Rhs(np) - Solution(x0,y)
      ELSE
      	   ! Find pos = position of P(x-dx,y)
      	   pos = np - Nm1
      	   IF ( ABS(Points(pos,1)-(x-dx)) > Eps ) THEN
      	   	   WRITE(*,*) " Error x-dx: ",ABS(Points(pos,1)-(x-dx))
      	   END IF
      	   A(np,pos) = 1.0D0
      END IF
      ! where is P(x+dx,y) ? 
      IF ( ABS((x+dx)-x1) < Eps ) THEN
      	   Rhs(np) = Rhs(np) - Solution(x1,y)
      ELSE
   	      ! Find pos = position of P(x+dx,y)
      	   pos = np + Nm1
      	   IF ( ABS(Points(pos,1)-(x+dx)) > Eps ) THEN
      	   	   WRITE(*,*) " Error x+dx: ",ABS(Points(pos,1)-(x+dx))
      	   END IF
      	   A(np,pos) = 1.0D0
      END IF
      ! where is P(x,y-dy) ? 
      IF ( ABS((y-dy)-y0) < Eps ) THEN
      	   Rhs(np) = Rhs(np) - Solution(x,y0)
      ELSE
      	   ! Find pos = position of P(x,y-dy)
           pos = np - 1
      	   IF ( ABS(Points(pos,2)-(y-dy)) > Eps ) THEN
      	   	   WRITE(*,*) " Error y-dy: ",ABS(Points(pos,2)-(y-dy))
      	   END IF
      	   A(np,pos) = 1.0D0
      END IF
      ! where is P(x,y+dy) ? 
      IF ( ABS((y+dy)-y1) < Eps ) THEN
      	   Rhs(np) = Rhs(np) - Solution(x,y1)
      ELSE
      	   ! Find pos = position of P(x,y+dy)
           pos = np + 1
      	   IF ( ABS(Points(pos,2)-(y+dy)) > Eps ) THEN
      	   	   WRITE(*,*) " Error y+dy: ",ABS(Points(pos,2)-(y+dy))
      	   END IF
      	   A(np,pos) = 1.0D0
      END IF
      IF ( LA <= 5 ) WRITE(*,"(A,I4,A,1X,F8.3)") " Rhs(",np,") = ",Rhs(np)
   END DO
   END SUBROUTINE EqsDef   
   
   SUBROUTINE PrintMatrix(LA,N,A,msg)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: LA, N
      REAL(8), DIMENSION(LA,N), INTENT(IN) :: A
      CHARACTER(*), INTENT(IN) :: msg
      INTEGER :: i, j

      IF ( LA*N <= 25 ) THEN
         WRITE(*,*) TRIM(msg)
         DO i = 1, LA
            DO j = 1, N
               WRITE(*,"(1X,F8.3)", ADVANCE="NO") A(i,j)
            END DO
            WRITE(*,"(/)")
         END DO
         WRITE(*,"(/)")
      END IF

      RETURN      
   END SUBROUTINE PrintMatrix
   
END PROGRAM LinearEquation