test program for dgesvx

dgesvx_tester.f90

! Test program for dgesvx
! compile command used: 
!ifort dgesvx_tester.f90 -L/opt/intel/composer2018/mkl/lib/intel64 -lmkl_core -lmkl_intel_lp64 -lmkl_sequential -lpthread
program main
  implicit none
  external :: dgemv
  external :: dgesv
  external :: dgesvx
  integer,  parameter   :: dp = selected_real_kind(15,300)

  real(dp), allocatable :: m1(:,:)
  real(dp), allocatable :: v1(:), v2(:)
  real(dp) :: rnorm 
  integer  :: sz, ii

  ! linear solver vars
  real(dp), allocatable :: af(:,:)
  integer,  allocatable :: ipiv(:)
  real(dp), allocatable :: rs(:)
  real(dp), allocatable :: cs(:)
  real(dp), allocatable :: ferr(:)
  real(dp), allocatable :: berr(:)
  real(dp), allocatable :: rwork(:)
  integer,  allocatable :: iwork(:)
  real(dp) :: rcond
  integer  :: info
  character(1) :: equed

  sz = 5

  ! matrix and vector vars    
  allocate( v1(sz) )
  allocate( v2(sz) )
  allocate( m1(sz,sz) )

  ! linear solve vars
  allocate( af(sz,sz) )
  allocate( ipiv(sz) )
  allocate( rs(sz) )
  allocate( cs(sz) )
  allocate( ferr(sz) )
  allocate( berr(sz) )
  allocate( rwork(4*sz) )
  allocate( iwork(sz) )


  ! run dgesv
  call setup(m1, v2, sz) ! init the matrix and vector
  call dgesv(sz, 1, m1, sz, ipiv, v2, sz, info)

  if( info /= 0 )then 
      write(*,*) 'ERR in dgesv'
  else
      ! manually check solution
      call setup(m1, v1, sz) ! reset matrix and vector
      call dgemv('N', sz, sz, 1.0_dp, m1, sz, v2, 1, -1.0_dp, v1, 1)
      rnorm = sum(v1*v1)

      write(*,'(1x,"dgesv solution,  ||A.x-b||:",g12.5)') rnorm
      write(*,'(3x,"   x: ")',advance='no')
      write(*,'(100(2x,f12.7))') (v2(ii), ii=1,sz)

  end if

  ! run dgesvx
  call setup(m1, v1, sz) ! reset matrix and vector
  equed = 'N'
  call dgesvx('N','N', sz, 1, m1, sz, af, sz, ipiv, &
              equed, rs, cs, v1, sz, v2, sz, rcond, &
              ferr, berr, rwork, iwork, info)

  if( info /= 0 )then 
      write(*,*) 'ERR in dgesvx'
  else
      ! manually check solution
      call setup(m1, v1, sz) ! reset matrix and vector
      call dgemv('N', sz, sz, 1.0_dp, m1, sz, v2, 1, -1.0_dp, v1, 1)
      rnorm = sum(v1*v1)

      write(*,'(1x,"dgesvx solution, ||A.x-b||:",g12.5," rcond: ",g12.4)') &
              rnorm, rcond
      write(*,'(3x,"   x: ")',advance='no')
      write(*,'(100(2x,f12.7))') (v2(ii), ii=1,sz)
      write(*,'(3x,"ferr: ")',advance='no')
      write(*,'(100(2x,f12.7))') (ferr(ii), ii=1,sz)
      write(*,'(3x,"berr: ")',advance='no')
      write(*,'(100(2x,f12.7))') (berr(ii), ii=1,sz)
  end if


  deallocate( v1 )
  deallocate( v2 )
  deallocate( m1 )

  deallocate( af )
  deallocate( ipiv )
  deallocate( rs )
  deallocate( cs )
  deallocate( ferr )
  deallocate( berr )
  deallocate( rwork )
  deallocate( iwork ) 

  stop
  contains
  subroutine setup(mio, vio, isz)
    implicit none
    real(dp), intent(inout) :: mio(:,:)
    real(dp), intent(inout) :: vio(:)
    integer,  intent(in)    :: isz
    integer :: jj
    real(dp) :: rt
    
    mio(:,:) = 0.0_dp
    do jj = 1, isz
       vio(jj) = real(jj,dp)
       mio(jj,jj) = 1.0_dp
    end do
    rt = 1.0_dp / sqrt(2.0_dp)
    mio(1:2,1) = rt
    mio(1,2) = rt
    mio(2,2) = -rt

  end subroutine setup
end program main