! test_inverse.f90  tough matrix to invert accurately

program test_inverse
  implicit none
  integer, parameter :: sz = 1024
  double precision, dimension(sz,sz) :: AA
  double precision, dimension(sz,sz) :: AI
  double precision, dimension(sz,sz) :: II
  double precision :: err, sum, val
  integer :: i, j, k, n

  interface
    subroutine inverse(n, sz, A, AI)
      integer, intent(in) :: n  ! number of equations
      integer, intent(in) :: sz ! dimension of arrays
      double precision, dimension(sz,sz), intent(in) :: A
      double precision, dimension(sz,sz), intent(in) :: AI
    end subroutine inverse 
  end interface

  print *, 'test_inverse.f90 on tough to invert matrix'
  n = 1
  do while(n .lt. 1024)
    n=n*2
    print *, 'initializing big matrix, n=',n,', n*n=', n*n
    do i=1,n
      do j=1,n
        AA(i,j) = 1.0d0/dble(i+j)
      end do ! j
    end do ! i
    if(n<=4) then
      print *, 'data'
      do i=1,n
        do j=1,n
          print *, 'AA(',i,',',j,')=',AA(i,j)
        end do ! j
      end do ! i
    end if 
    call inverse( n, sz, AA , AI )
    if(n<=4) then
      print *, 'inverse'
      do i=1,n
        do j=1,n
          print *, 'AI(',i,',',j,')=',AI(i,j)
        end do ! j
      end do ! i
    end if 
    err = 0.0D0
    do i=1,n
      do j=1,n
        II(i,j) = 0.0D0
        do k=1,n
          II(i,j) = II(i,j) + AA(i,k)*AI(k,j)
        end do ! k
        if(i .eq. j) then
          err = err + abs(II(i,j)-1.0D0)
        else
          err = err + abs(II(i,j))
        end if
      end do ! j
    end do ! i
    print *, 'sum of error=',err,', avg error=', err/dble(n*n)
  end do ! while
  print *, 'finished test_inverse.f90'
end program test_inverse