Module Precision_Module
   integer,parameter :: $rp=8
   integer,parameter :: $ip=8
   Type RowType
      integer($ip) :: n
      integer($ip) :: CenterDiag
      integer($ip),allocatable :: ID(:)
      real($rp),allocatable  :: a(:)
   End Type RowType

   integer($ip) :: NoOfSuperDiag
   integer($ip) :: NoOfSubDiag
   integer($ip) :: NoOfBands

   integer($ip),allocatable :: BandsID(:)

End Module Precision_Module


Module Vectors
      use Precision_Module
      implicit none
Contains
   Function MatTransposeMultVect(a,b)
      implicit none
      real($rp),intent(in), dimension(:,:) :: a
      real($rp),intent(in), dimension(:) :: b
      real($rp) :: MatTransposeMultVect(size(a,1))

      real($rp) :: sum

      integer(4) :: n
      integer(4) :: m
      integer(4) :: p
      integer(4) :: i
      integer(4) :: j
      integer(4) :: k

      m=size(a,1)
      n=size(a,2)
      p=size(b,1)

      MatTransposeMultVect=0.0_$rp

      if(m .ne. p) return

      do j=1,n
         sum=0.0_$rp
         do i=1,m
               sum=sum+a(i,j)*b(i)
         end do
         MatTransposeMultVect(j)=sum
      end do
      return
   End Function MatTransposeMultVect


   Function MatSparseMultVect(row,b)
      implicit none
      type(RowType),intent(in), dimension(:) :: row
      real($rp),intent(in), dimension(:) :: b
      real($rp) :: MatSparseMultVect(size(b))

      real($rp) :: sum

      integer(4) :: n
      integer(4) :: m
      integer(4) :: p
      integer(4) :: i
      integer(4) :: j
      integer(4) :: k

      m=size(b,1)

      MatSparseMultVect=0.0_$rp

      do i=1,m
         sum=0.0_$rp
         n=row(i)%n
         do j=1,n
               sum=sum+row(i)%a(j)*b(row(i)%ID(j))
         end do
         MatSparseMultVect(i)=sum
      end do
      return
   End Function MatSparseMultVect
End Module Vectors

Module ErrorModule
Contains
   Subroutine Error_BMFS(n,x,row,b,errPk,errL2)

      use Precision_Module
      use Vectors

      implicit none

      integer($ip)   :: n
      real($rp)      :: x(:)
      type(RowType)  :: row(:)
      real($rp)      :: b(:)

      real($rp)      :: errpk
      real($rp)      :: errL2

      real($rp),allocatable :: err(:)

      allocate (err(n))

      errpk = 0.0_$rp
      err = 0.0_$rp

      err = b -  MatSparseMultVect(row, x) 
      errpk = maxval(abs(err))
      errL2 = sqrt ( sum ( err(:n)*err(:n) ))   ! L2 norm

      deallocate(err)

      return
   End Subroutine Error_BMFS

   Subroutine Error_Pardiso(n,      &
                            nelem,  &
                            ia,     &
                            ja,     &
                            a,      &
                            b,      &
                            x,      &
                            errPk,  &
                           errL2)
      use Precision_Module
      IMPLICIT NONE
      include 'mkl.fi'

      integer(8) :: n
      integer(8) :: nelem
      integer(8) :: ia(n)
      integer(8) :: ja(nelem+1)
      real(8) :: a(nelem)
      real(8) :: b(n)
      real(8) :: x(n)

      real(8),allocatable :: y(:)
      real(8),allocatable :: err(:)

      real(8) :: errpk
      real(8) :: errL2

      allocate (y(n))
      allocate (err(n))

      errpk = 0.0_$rp
      err = 0.0_$rp

      call mkl_dcsrgemv('n', n, a, ia, ja, x, y)

      err = b - y
      errpk = maxval(abs(err))
      errL2 = sqrt ( sum ( err(:n)*err(:n) ))   ! L2 norm

      deallocate(y)
      deallocate(err)

      return
   End Subroutine Error_Pardiso
End Module ErrorModule


!*******************************************************************************
!   Copyright(C) 2004-2015 Intel Corporation. All Rights Reserved.
!
!   The source code, information  and  material ("Material") contained herein is
!   owned  by Intel Corporation or its suppliers or licensors, and title to such
!   Material remains  with Intel Corporation  or its suppliers or licensors. The
!   Material  contains proprietary information  of  Intel or  its  suppliers and
!   licensors. The  Material is protected by worldwide copyright laws and treaty
!   provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
!   modified, published, uploaded, posted, transmitted, distributed or disclosed
!   in any way  without Intel's  prior  express written  permission. No  license
!   under  any patent, copyright  or  other intellectual property rights  in the
!   Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
!   implication, inducement,  estoppel or  otherwise.  Any  license  under  such
!   intellectual  property  rights must  be express  and  approved  by  Intel in
!   writing.
!
!   *Third Party trademarks are the property of their respective owners.
!
!   Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
!   this  notice or  any other notice embedded  in Materials by Intel or Intel's
!   suppliers or licensors in any way.
!
!*******************************************************************************
!   Content : MKL PARDISO Fortran-90 header file
!
!           Contains PARDISO routine definition.
!           For CDECL use only.
!
!*******************************************************************************
!DEC$ IF .NOT. DEFINED( __MKL_PARDISO_F90 )

!DEC$ DEFINE __MKL_PARDISO_F90

      MODULE MKL_PARDISO_PRIVATE

        TYPE MKL_PARDISO_HANDLE; INTEGER(KIND=8) DUMMY; END TYPE

        INTEGER, PARAMETER :: PARDISO_OOC_FILE_NAME = 1

      END MODULE MKL_PARDISO_PRIVATE

      MODULE MKL_PARDISO
        USE MKL_PARDISO_PRIVATE

!
! Subroutine prototype for PARDISO
!

      INTERFACE PARDISO
        SUBROUTINE PARDISO_S( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR

          REAL(KIND=4),     INTENT(IN)    :: A(*)
          REAL(KIND=4),     INTENT(INOUT) :: B(*)
          REAL(KIND=4),     INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_S

        SUBROUTINE PARDISO_D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR

          REAL(KIND=8),     INTENT(IN)    :: A(*)
          REAL(KIND=8),     INTENT(INOUT) :: B(*)
          REAL(KIND=8),     INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_D

        SUBROUTINE PARDISO_SC( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR

          COMPLEX(KIND=4),  INTENT(IN)    :: A(*)
          COMPLEX(KIND=4),  INTENT(INOUT) :: B(*)
          COMPLEX(KIND=4),  INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_SC

        SUBROUTINE PARDISO_DC( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=8),  INTENT(IN)    :: A(*)
          COMPLEX(KIND=8),  INTENT(INOUT) :: B(*)
          COMPLEX(KIND=8),  INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_DC


        SUBROUTINE PARDISO_S_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR
          REAL(KIND=4),     INTENT(IN)    :: A(*)
          REAL(KIND=4),     INTENT(INOUT) :: B(N,*)
          REAL(KIND=4),     INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_S_2D

        SUBROUTINE PARDISO_D_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR
          REAL(KIND=8),     INTENT(IN)    :: A(*)
          REAL(KIND=8),     INTENT(INOUT) :: B(N,*)
          REAL(KIND=8),     INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_D_2D

        SUBROUTINE PARDISO_SC_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=4),  INTENT(IN)    :: A(*)
          COMPLEX(KIND=4),  INTENT(INOUT) :: B(N,*)
          COMPLEX(KIND=4),  INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_SC_2D

        SUBROUTINE PARDISO_DC_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MAXFCT
          INTEGER,          INTENT(IN)    :: MNUM
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(IN)    :: PHASE
          INTEGER,          INTENT(IN)    :: N
          INTEGER,          INTENT(IN)    :: IA(*)
          INTEGER,          INTENT(IN)    :: JA(*)
          INTEGER,          INTENT(INOUT) :: PERM(*)
          INTEGER,          INTENT(IN)    :: NRHS
          INTEGER,          INTENT(INOUT) :: IPARM(*)
          INTEGER,          INTENT(IN)    :: MSGLVL
          INTEGER,          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=8),  INTENT(IN)    :: A(*)
          COMPLEX(KIND=8),  INTENT(INOUT) :: B(N,*)
          COMPLEX(KIND=8),  INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_DC_2D
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_STORE( PT, STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_STORE
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_STORE_64( PT, STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_STORE_64
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_RESTORE( PT, STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_RESTORE
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_RESTORE_64( PT, STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_RESTORE_64
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_DELETE( STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_DELETE
      END INTERFACE

      INTERFACE
        SUBROUTINE PARDISO_HANDLE_DELETE_64( STORE_NAME, ERROR )
          USE MKL_PARDISO_PRIVATE
          CHARACTER, INTENT(IN)  :: STORE_NAME(*)
          INTEGER,   INTENT(OUT) :: ERROR
        END SUBROUTINE PARDISO_HANDLE_DELETE_64
      END INTERFACE

!
! Subroutine prototype for PARDISO_64
!
! Note: The pardiso_64 interface is not supported on IA-32 architecture.
!       If called on IA-32, error = -12 is returned.
!
      INTERFACE PARDISO_64
        SUBROUTINE PARDISO_S_64( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          REAL(KIND=4),             INTENT(IN)    :: A(*)
          REAL(KIND=4),             INTENT(INOUT) :: B(*)
          REAL(KIND=4),             INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_S_64

        SUBROUTINE PARDISO_D_64( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          REAL(KIND=8),             INTENT(IN)    :: A(*)
          REAL(KIND=8),             INTENT(INOUT) :: B(*)
          REAL(KIND=8),             INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_D_64

        SUBROUTINE PARDISO_SC_64( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=4),          INTENT(IN)    :: A(*)
          COMPLEX(KIND=4),          INTENT(INOUT) :: B(*)
          COMPLEX(KIND=4),          INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_SC_64

        SUBROUTINE PARDISO_DC_64( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=8),          INTENT(IN)    :: A(*)
          COMPLEX(KIND=8),          INTENT(INOUT) :: B(*)
          COMPLEX(KIND=8),          INTENT(OUT)   :: X(*)
        END SUBROUTINE PARDISO_DC_64

        SUBROUTINE PARDISO_S_64_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          REAL(KIND=4),             INTENT(IN)    :: A(*)
          REAL(KIND=4),             INTENT(INOUT) :: B(N,*)
          REAL(KIND=4),             INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_S_64_2D

        SUBROUTINE PARDISO_D_64_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          REAL(KIND=8),             INTENT(IN)    :: A(*)
          REAL(KIND=8),             INTENT(INOUT) :: B(N,*)
          REAL(KIND=8),             INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_D_64_2D

        SUBROUTINE PARDISO_SC_64_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=4),          INTENT(IN)    :: A(*)
          COMPLEX(KIND=4),          INTENT(INOUT) :: B(N,*)
          COMPLEX(KIND=4),          INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_SC_64_2D

        SUBROUTINE PARDISO_DC_64_2D( PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, PERM, NRHS, IPARM, MSGLVL, B, X, ERROR )
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MAXFCT
          INTEGER(KIND=8),          INTENT(IN)    :: MNUM
          INTEGER(KIND=8),          INTENT(IN)    :: MTYPE
          INTEGER(KIND=8),          INTENT(IN)    :: PHASE
          INTEGER(KIND=8),          INTENT(IN)    :: N
          INTEGER(KIND=8),          INTENT(IN)    :: IA(*)
          INTEGER(KIND=8),          INTENT(IN)    :: JA(*)
          INTEGER(KIND=8),          INTENT(INOUT) :: PERM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: NRHS
          INTEGER(KIND=8),          INTENT(INOUT) :: IPARM(*)
          INTEGER(KIND=8),          INTENT(IN)    :: MSGLVL
          INTEGER(KIND=8),          INTENT(OUT)   :: ERROR
          COMPLEX(KIND=8),          INTENT(IN)    :: A(*)
          COMPLEX(KIND=8),          INTENT(INOUT) :: B(N,*)
          COMPLEX(KIND=8),          INTENT(OUT)   :: X(N,*)
        END SUBROUTINE PARDISO_DC_64_2D

      END INTERFACE

      INTERFACE

       SUBROUTINE PARDISOINIT(PT, MTYPE, IPARM)
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(OUT) :: PT(*)
          INTEGER,          INTENT(IN)    :: MTYPE
          INTEGER,          INTENT(OUT) :: IPARM(*)
       END SUBROUTINE PARDISOINIT

      END INTERFACE

      INTERFACE PARDISO_GET

       FUNCTION PARDISO_GETENV(PT, OptName, StrVal)
          USE MKL_PARDISO_PRIVATE
          INTEGER PARDISO_GETENV
          TYPE(MKL_PARDISO_HANDLE), INTENT(IN)  :: PT(*)
          INTEGER,                  INTENT(IN)  :: OptName
          CHARACTER(*),             INTENT(OUT)  :: StrVal
       END FUNCTION PARDISO_GETENV

      END INTERFACE

      INTERFACE PARDISO_SET

       FUNCTION PARDISO_SETENV(PT, OptName, StrVal)
          USE MKL_PARDISO_PRIVATE
          INTEGER PARDISO_SETENV
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          INTEGER,                  INTENT(IN)    :: OptName
          INTEGER,                  INTENT(IN)    :: StrVal(*)
       END FUNCTION PARDISO_SETENV

      END INTERFACE

      INTERFACE PARDISO_PIV
      
      FUNCTION MKL_PARDISO_PIVOT( AII, BII, EPS)
         REAL(KIND=8)   :: AII, BII, EPS
         INTEGER   MKL_PARDISO_PIVOT
         END
      END INTERFACE PARDISO_PIV

      INTERFACE PARDISO_GETDIAG

         SUBROUTINE PARDISO_GETDIAG_D(PT, DIAG_FACT, DIAG_A, MNUM, EPS)
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          REAL(KIND=8), INTENT(INOUT)             :: DIAG_FACT, DIAG_A
          INTEGER, INTENT(IN)            :: MNUM
          INTEGER, INTENT(INOUT)         :: EPS
         END

         SUBROUTINE PARDISO_GETDIAG_Z(PT, DIAG_FACT, DIAG_A, MNUM, EPS)
          USE MKL_PARDISO_PRIVATE
          TYPE(MKL_PARDISO_HANDLE), INTENT(INOUT) :: PT(*)
          COMPLEX(KIND=8), INTENT(INOUT)          :: DIAG_FACT, DIAG_A
          INTEGER, INTENT(IN)            :: MNUM
          INTEGER, INTENT(INOUT)         :: EPS
         END

      END INTERFACE PARDISO_GETDIAG




      END MODULE MKL_PARDISO

!DEC$ ENDIF
      Subroutine pardiso_unsym(n,      &
                               nelem,  &
                               ia,     &
                               ja,     &
                               a,      &
                               b,      &
                               x,      &
                               timer)
         use Precision_Module
         use  MKL_PARDISO
         use omp_lib

         IMPLICIT NONE

         !integer(4) :: n8
         !integer(4) :: nelem8
         !integer(4) :: ia8(n8)
         !integer(4) :: ja8(nelem8+1)

         integer(8) :: n
         integer(8) :: nelem
         integer(8) :: ia(n)
         integer(8) :: ja(nelem+1)
         real(8) :: a(nelem)
         real(8) :: b(n)
         real(8) :: x(n)
         real(8) :: timer

         real(8) :: tic
         real(8) :: toc

!C.. Internal solver memory pointer for 64-bit architectures
!C.. INTEGER*8 pt(64)
!C.. Internal solver memory pointer for 32-bit architectures
!C.. INTEGER*4 pt(64)
!C.. This is OK in both cases
         TYPE(MKL_PARDISO_HANDLE) :: pt(64)
!C.. All other variables
         INTEGER(8) :: maxfct
         INTEGER(8) :: mnum
         INTEGER(8) :: mtype
         INTEGER(8) :: phase
         INTEGER(8) :: nrhs
         INTEGER(8) :: error
         INTEGER(8) :: msglvl
         INTEGER(8) :: iparm(64)


         INTEGER(8) :: i
         INTEGER(8) :: idum(1)
         REAL(8) :: ddum(1)
!C.. Fill all arrays containing matrix data.
         DATA nrhs /1/, maxfct /1/, mnum /1/

!C..
!C.. Setup PARDISO control parameter
!C..
         iparm    = 0
         iparm(1) = 1 ! no solver default
         iparm(2) = 2 ! fill-in reordering from METIS
         iparm(3) = 8 ! numbers of processors
         iparm(4) = 0 ! no iterative-direct algorithm
         iparm(5) = 0 ! no user fill-in reducing permutation
         iparm(6) = 0 ! =0 solution on the first n components of x
         iparm(7) = 0 ! not in use
         iparm(8) = 9 ! numbers of iterative refinement steps
         iparm(9) = 0 ! not in use
         iparm(10) = 13 ! perturb the pivot elements with 1E-13
         iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
         iparm(12) = 0 ! not in use
         iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
         iparm(14) = 0 ! Output: number of perturbed pivots
         iparm(15) = 0 ! not in use
         iparm(16) = 0 ! not in use
         iparm(17) = 0 ! not in use
         iparm(18) = -1 ! Output: number of nonzeros in the factor LU
         iparm(19) = -1 ! Output: Mflops for LU factorization
         iparm(20) = 0 ! Output: Numbers of CG Iterations
         error = 0 ! initialize error flag
         msglvl = 0!1 ! print statistical information
         mtype = 11 ! real unsymmetric

         !n     = n8
         !nelem = nelem8
         !ia    = ia8
         !ja    = ja8

         tic = omp_get_wtime()

!C.. Initialize the internal solver memory pointer. This is only
!C necessary for the FIRST call of the PARDISO solver.


         pt = MKL_PARDISO_HANDLE(0)

!C.. Reordering and Symbolic Factorization, This step also allocates
!C   all memory that is necessary for the factorization

         phase = 11 ! only reordering and symbolic factorization
         CALL pardiso_64 (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)

         !WRITE(*,*) 'Reordering completed ... '
         IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
         END IF

         !WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
         !WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)

!C.. Factorization.
         phase = 22 ! only factorization
         CALL pardiso_64 (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)

         !WRITE(*,*) 'Factorization completed ... '
         IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
         END IF


!C.. Back substitution and iterative refinement
         iparm(8) = 2 ! max numbers of iterative refinement steps
         phase = 33 ! only factorization
         DO i = 1, n
            b(i) = 1.0D0
         END DO

         CALL pardiso_64 (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, b, x, error)
         !WRITE(*,*) 'Solve completed ... '
         !WRITE(*,*) 'The solution of the system is '
         !DO i = 1, n
         !   WRITE(*,*) ' x(',i,') = ', x(i)
         !END DO


!C.. Termination and release of memory
         phase = -1 ! release internal memory
         CALL pardiso_64 (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, idum, nrhs, iparm, msglvl, ddum, ddum, error)

         toc = omp_get_wtime()
         timer = toc - tic

         return
      END Subroutine pardiso_unsym
Module BFMS_Transpose1
Contains

   Subroutine BMFS_Algorithm(x,row,b,p,q,n)

      use Precision_Module
      use lin_sol_gen_int
      implicit none
    
      integer($ip) :: p
      integer($ip) :: q
      integer($ip) :: n
      integer($ip) :: i
      integer($ip) :: j
      integer($ip) :: k

      type(RowType)  :: row(:)
      real($rp)      :: b(:)
      real($rp)      :: x(:)
    
      real($rp),allocatable :: alpha(:,:)
      real($rp),allocatable :: u(:,:)
      real($rp),allocatable :: zeta(:)         ! basis vector
      real($rp),allocatable :: RM(:,:)           ! space for constraint system RHS
      real($rp),allocatable :: RMx(:,:)           ! space for constraint system RHS
      real($rp),allocatable :: PM(:,:)         ! space for constraint system LHS

      real($rp) :: tsum = 0.0_$rp

!$OMP SECTIONS 
      allocate(alpha(q,q+1))
      allocate(u(n,q+1))
      allocate(zeta(q))         ! basis vector
      allocate(RM(q,1))           ! space for constraint system RHS
      allocate(RMx(q,1))           ! space for constraint system RHS
      allocate(PM(q,q))         ! space for constraint system LHS
!$OMP END SECTIONS

!$OMP SECTIONS 
      zeta           = 1.0_$rp
      u(1:q,1:q+1)   = 1.0_$rp
      u(q+1:n,1:q+1) = 0.0_$rp
      alpha          = 0.0_$rp
      !
      !write(6,*)'calling MFS with b vector'
!$OMP END SECTIONS
      k=1
      call MFS(u,alpha,zeta,row,b,p,q,n,k)

!$OMP PARALLEL NUM_THREADS(8) SHARED(u,alpha,zeta,row,p,q,n)
   !$OMP DO SCHEDULE(STATIC) PRIVATE(k)

      do k=2,q+1  !concurrent (k=2:q+1)
         !write(6,*)'calling MFS with k=',k
         call FindU(u,alpha,zeta,row,p,q,n,k)
      end do

  !$OMP END DO NOWAIT 
!$OMP END PARALLEL 

      RM(:,1) = alpha(:,1)


!$OMP PARALLEL NUM_THREADS(8) SHARED(PM,zeta,alpha,q)
   !$OMP DO SCHEDULE(STATIC) PRIVATE(i,j)

      do i=1,q !concurrent (i = 1:q)
         do j=1,q !concurrent (j = 1:q)
            PM(i,j) =         - zeta(j)*alpha(i,j+1)
         end do
      end do

      !do i=1,q !concurrent (i = 1:q)
      !   do j=1,q !concurrent (j = 1:q)
      !      if (i .eq. j) then
      !         PM(i,j) = 1.0_$rp - zeta(j)*alpha(i,j+1)
      !      else
      !         PM(i,j) =         - zeta(j)*alpha(i,j+1)
      !      end if
      !   end do
      !end do

  !$OMP END DO NOWAIT 
!$OMP END PARALLEL 

!$OMP PARALLEL NUM_THREADS(8) SHARED(PM,zeta,alpha,q)
   !$OMP DO SCHEDULE(STATIC) PRIVATE(i,j)

      do i=1,q !concurrent (i = 1:q)
         PM(i,i) = 1.0_$rp - zeta(i)*alpha(i,i+1)
      end do

  !$OMP END DO NOWAIT 
!$OMP END PARALLEL 

      CALL LIN_SOL_GEN(PM,RM,RMx)


!$OMP PARALLEL NUM_THREADS(8) SHARED(x,u,zeta,RMx,q,n)
   !$OMP DO SCHEDULE(STATIC) PRIVATE(i,j,tsum)

      do i = 1,n
         tsum = 0.0_$rp
         do j = 1,q
            tsum = tsum + zeta(j)*RMx(j,1)*u(i,j+1)
         end do
         x(i) = u(i,1) + tsum
      end do

  !$OMP END DO NOWAIT 
!$OMP END PARALLEL 

      deallocate(alpha)
      deallocate(u)
      deallocate(zeta)        ! basis vector
      deallocate(RM)          ! space for constraint system RHS
      deallocate(RMx)          ! space for constraint system RHS
      deallocate(PM)          ! space for constraint system LHS

      return
   End Subroutine BMFS_Algorithm

   Pure Subroutine FindU(u,alpha,zeta,row,p,q,n,k)
      use Precision_Module
      implicit none
      real($rp),intent(inout) :: u(:,:)
      real($rp),intent(out)   :: alpha(:,:)
      real($rp),intent(in)    :: zeta(:)
      type(RowType),intent(in):: row(:)
      integer($ip),intent(in) :: p
      integer($ip),intent(in) :: q
      integer($ip),intent(in) :: n
      integer($ip),intent(in) :: k

      integer($ip) :: i,j

      real($rp),allocatable :: e(:)       ! temporary basis vector
      allocate(e(n))                      ! temporary basis vector

      if(q .ne. 1)u(k-1,k)=0.0_$rp
      e = 0.0_$rp
      e(n-q+k-1) = 1.0_$rp
      call MFS(u,alpha,zeta,row,e,p,q,n,k)
      deallocate(e)
      return
   End Subroutine FindU

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   Pure Subroutine MFS(u,alpha,zeta,row,b,p,q,n,k)
      use Precision_Module
      implicit none

      real($rp),intent(inout) :: u(:,:)
      real($rp),intent(out)   :: alpha(:,:)
      real($rp),intent(in)    :: zeta(:)
      type(RowType),intent(in):: row(:)
      real($rp),intent(in)    :: b(:)

      integer($ip),intent(in) :: p
      integer($ip),intent(in) :: q
      integer($ip),intent(in) :: n
      integer($ip),intent(in) :: k
      integer($ip)            :: irow

      real($rp) :: tsum
      integer($ip) :: i, j, m, nBands

      do m=1,n-q
         tsum = 0.0_$rp
         nBands=row(m)%n
         do j=1,nBands-1
            tsum = tsum + row(m)%a(j) * u(row(m)%ID(j),k)
         end do
         u(row(m)%ID(nBands),k) = (b(m) - tsum) / row(m)%a(nBands)
      end do

!$OMP PARALLEL NUM_THREADS(8) SHARED(q,n,row,u,k,b,zeta,alpha)
   !$OMP DO SCHEDULE(STATIC) PRIVATE(i,j,tsum,m)

      do i=1,q
         m=n-q+i
         tsum = 0.0_$rp
         do j=1,row(m)%n
            tsum=tsum+row(m)%a(j)*u(row(m)%ID(j),k)
         end do
         alpha(i,k) = (b(m) - tsum)/zeta(i)
      end do

  !$OMP END DO NOWAIT 
!$OMP END PARALLEL 
      return
   End Subroutine  MFS


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

End Module BFMS_Transpose1

Module BFMS_Interface
Contains

   Subroutine call_BMFS_Algorithm(x,row,b,p,q,n,timer)

      use Precision_Module
      use BFMS_Transpose1
      use ifport
      use omp_lib
      use Vectors
      implicit none

      real($rp)      :: x(:)
      type(RowType)  :: row(:)
      real($rp)      :: b(:)
      real($rp)      :: timer

      integer($ip) :: n
      integer($ip) :: p
      integer($ip) :: q

      real($rp)      :: tic
      real($rp)      :: toc

      integer($ip) :: i
      integer($ip) :: j


      tic = omp_get_wtime()
      call BMFS_Algorithm(x,row,b,p,q,n)
      toc = omp_get_wtime()
      timer=toc - tic

      return
   End Subroutine call_BMFS_Algorithm
End Module BFMS_Interface


Program RunBenchmark
   use Precision_Module
   use BFMS_Interface
   use ErrorModule
   implicit none

   integer($ip) :: m
   integer($ip) :: i
   integer($ip) :: j
   integer($ip) :: k
   integer($ip) :: icase
   integer($ip) :: jcase
   integer($ip) :: BandID

   real($rp) :: timer_BMFS
   real($rp) :: errPK_BMFS
   real($rp) :: errL2_BMFS

   real($rp) :: timer_Pardiso
   real($rp) :: errPK_Pardiso
   real($rp) :: errL2_Pardiso

   real($rp) :: acoeff
   real($rp) :: bcoeff

   integer($ip) :: n
   integer($ip) :: p
   integer($ip) :: q

   integer($ip) :: L
   integer($ip) :: U

   integer($ip) :: NDiag

   integer($ip) :: TotalNumberOfCoeff

   real($rp),allocatable :: b(:)
   real($rp),allocatable :: x_BMFS(:)
   type(RowType),allocatable :: row(:)


   integer($ip) :: TotalNumberOfElements
   integer($ip) :: iElement
   integer($ip),allocatable :: RowIndex(:)
   integer($ip),allocatable :: ColIndex(:)
   real($rp),allocatable :: a(:)
   real($rp),allocatable :: x_Pardiso(:)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   open(3,file='fort.3',status='unknown',ACCESS='append')
   NoOfSuperDiag=1
   NoOfSubDiag=1

   NoOfbands = NoOfSubDiag + NoOfSuperDiag + 1
   allocate(BandsID(-NoOfSubDiag:NoOfSuperDiag))

   acoeff = -1
   bcoeff = 2
   do jcase=1,8
      select case(jcase)
         case(1)
            n=1
         case(2)
            n=10
         case(3)
            n=20
         case(4)
            n=50
         case(5)
            n=100
         case(6)
            n=500
         case(7)
            n=1000
         case(8)
            n=5000
      end select
      n=n*1000
      do icase=1,12
         select case(icase)
            case(1)
               BandId=1
               BandsID=[-BandId,0,BandId]
            case(2)
               BandId=100
               BandsID=[-BandId,0,BandId]
            case(3)
               BandId=1000
               BandsID=[-BandId,0,BandId]
            case(4)
               BandId=2000
               BandsID=[-BandId,0,BandId]
            case(5)
               BandId=3000
               BandsID=[-BandId,0,BandId]
            case(6)
               BandId=4000
               BandsID=[-BandId,0,BandId]
            case(7)
               BandId=5000
               BandsID=[-BandId,0,BandId]
            case(8)
               BandId=6000
               BandsID=[-BandId,0,BandId]
            case(9)
               BandId=7000
               BandsID=[-BandId,0,BandId]
            case(10)
               BandId=8000
               BandsID=[-BandId,0,BandId]
            case(11)
               BandId=9000
               BandsID=[-BandId,0,BandId]
            case(12)
               BandId=10000
               BandsID=[-BandId,0,BandId]
         end select

         if(BandId .gt. n)exit

         p = minval(BandsID)  !NoOfSubDiag   !
         q = maxval(BandsID)  !NoOfSuperDiag !

         allocate (row(n))

         Ndiag=NoOfSuperDiag+1
         L=-1
         U=NoOfSuperDiag
         TotalNumberOfElements=0

         !write(6,*)' hello here we go'

         do i = 1,n
            if(i .eq. abs(BandsID(L))+1)then
               Ndiag=Ndiag+1
               L=max(L-1,-NoOfSubDiag)
            end if
            if(i .eq. n-BandsID(U)+1)then
               Ndiag=Ndiag-1
               U=max(U-1,1)
            end if
            row(i)%n = Ndiag

            allocate (row(i)%a(row(i)%n))
            allocate (row(i)%ID(row(i)%n))

            TotalNumberOfElements=TotalNumberOfElements+row(i)%n

         end do

         !write(6,*)' done finding number of diagonals per row'

         allocate (b(n))
         allocate (a(TotalNumberOfElements))

         allocate (x_BMFS(n))
         allocate (x_Pardiso(TotalNumberOfElements))

         allocate (RowIndex(n+1))
         allocate (ColIndex(TotalNumberOfElements))

         x_BMFS      = 0.0_$rp
         x_Pardiso   = 0.0_$rp
         b           = 0.0_$rp

         iElement=1
         RowIndex(1)=1
         do i = 1,n
            b(i) = 1.0_$rp
            m=1
            do j = 1,n
               do k=-NoOfSubDiag,NoOfSuperDiag
                  if(j .eq. BandsID(k)+i)then
                     if(k .eq. 0)then
                        row(i)%a(m) = bcoeff
                        row(i)%CenterDiag = m
                     else
                        row(i)%a(m) = acoeff
                     end if

                     row(i)%ID(m)=j
                     a(iElement)=row(i)%a(m)
                     ColIndex(iElement)=j
                     iElement=iElement+1
                     RowIndex(i+1)=iElement
                     m=m+1

                  end if
               end do
            end do
         end do

         write(6,*)'bands assembled   -   entering BMFS_Algorithm'

         call call_BMFS_Algorithm(x_BMFS,row,b,p,q,n,timer_BMFS)

         call Error_BMFS(n,x_BMFS,row,b,errPk_BMFS,errL2_BMFS)

         call pardiso_unsym(n,iElement-1,RowIndex,ColIndex,a,b,x_Pardiso,timer_Pardiso)
         
         call Error_Pardiso(n,iElement-1,RowIndex,ColIndex,a,b,x_Pardiso,errPk_Pardiso,errL2_Pardiso)

         write(3,1000)n,q,float(q)/float(n),timer_BMFS,timer_Pardiso,ErrPk_BMFS,ErrPk_Pardiso,ErrL2_BMFS,ErrL2_Pardiso
         write(6,1000)n,q,float(q)/float(n),timer_BMFS,timer_Pardiso,ErrPk_BMFS,ErrPk_Pardiso,ErrL2_BMFS,ErrL2_Pardiso

         !write(3,*) 'x='
         !write(3,'(1x,i3,f13.3,1x,f13.3)')(i,x(i),xPardiso(i),i=1,n)

         deallocate (row)
         deallocate (RowIndex)
         deallocate (ColIndex)
         deallocate (a)

         deallocate (b)
         deallocate (x_BMFS)
         deallocate (x_Pardiso)
      end do
   end do
1000 format(1x,i9,1x,i5,7(1x,G13.5))
   stop
End Program RunBenchmark

