program multiplication

      integer, parameter :: size= 20
      external:: dgesv        ! subroutines from BLAS
      integer, dimension(:), allocatable :: ipiv
      real, dimension(:,:), allocatable :: A,AA, b
      real, dimension(:), allocatable :: x, r, rr
      real :: t1, t2, t3, tt, mflops, trace, diff
      integer :: n, loops, i,j,l, case, m

      real, dimension(:,:), allocatable :: D_AF, xx
      real, dimension(:), allocatable :: D_R, D_C, WORK , FERR, BERR
      integer, dimension(:), allocatable ::  IWORK
      integer :: LWORK
      real ::  RCOND
      character*1 :: EQUED

      m = 1
      
      allocate( A(size, size), AA(size, size) )
      allocate( x(1:size), b(1:size, 1:m),  r(1:size), rr(1:size), ipiv (1:size) )


      n = size

      do k=1, n
         ipiv(k) = k
      enddo


      
      loops = 1

      A = 0.
      do i = 1, n
         !do j = 1, n
         do j = max(i-1, 1), min(i+1, n)
            A(i,j) = i+j
         end do
      end do

      do i = 1, n
         rr(i) = 1.
      enddo
      
      b(1:n, 1) = rr(1:n)
      
        
      AA(:,:) = A(:,:)


      call DGESV(n, m, AA(1:n, 1:n), n, ipiv, b(1:n, 1:m), n, info)
      if(INFO /= 0) print*,'Warning in lap_sub.f90 LAPACK (A): DGESV,  INFO = ',info


      x(1:n) = b(1:n, 1)

      ! residual
      r(1:n) = rr(1:n) - matmul( A(1:n, 1:n), x(1:n) )

      do i=1, n
         write(*,'(i5,a2, 2es12.4)') i, ': ', x(i), r(i)
         write(11,'(i5,a2, 2es12.4)') i, ': ', x(i), r(i)
      enddo
      

      ! VERSION with error estimates

      LWORK = 4 * n
      allocate(D_AF(1:n, 1:n), D_R(1:n), D_C(1:n), xx(1:n, 1:m) )
      allocate(FERR(1:m), BERR(1:m), WORK(LWORK), IWORK(n) )

      AA(:,:) = A(:,:)
      b(1:n, 1) = rr(1:n)
      
      call DGESVX( 'E', 'N', n, m, AA(1:n, 1:n), n, D_AF(1:n, 1:n), n, ipiv(1:n),  &
           EQUED, D_R(1:n), D_C(1:n), b(1:n, 1:m), n, xx(1:n, 1:m), n, RCOND, FERR, BERR, &
           WORK, IWORK, INFO )
      if(INFO /= 0) print*,'Warning in lap_sub.f90 LAPACK (A): DGESVX,  INFO = ',info

      ! residual
      r(1:n) = rr(1:n) - matmul( A(1:n, 1:n), xx(1:n,1) )

      do i=1, n
         write(*,'(i5,a2, 8es12.4)') i, ': ', xx(i,1), r(i), ferr(1), berr(1), rcond, 1./rcond
         write(12,'(i5,a2, 2es12.4)') i, ': ', xx(i,1), r(i)
      enddo

      
      
      deallocate(A, AA, b, r, x)

    end program Multiplication