REAL(4) H(4,4),A(4,4),A_(4,4),E(4,4), beta,cond,work(4),R(4,4)
      integer ndim/4/,n/4/,ipvt(4), i, j, k
      DATA H/2,1,0,0,
     *       4,2,1,0,
     *       8,4,2,1,
     *       16,8,4,2/
      DATA E/1,0,0,0,
     *       0,1,0,0,
     *       0,0,1,0,
     *       0,0,0,1/
      DATA A_/1,0,0,0,
     *       0,1,0,0,
     *       0,0,1,0,
     *       0,0,0,1/

      data beta/0.1/

      H = transpose(H)
c      A_ = transpose(E)
c      A = transpose(A)
      E = transpose(E)
      print *, "#################"
      print 102, beta
      print *, "#################"

      do i=1,4
          do j=1,4
              A(i,j) = beta * E(i,j) + H(i,j)
          end do
      end do
      print *, "matrix A = beta * E + H"
      print 101,((a(i,j),j=1,n), i=1,n)

      call decomp(ndim,n,a,cond,ipvt,work)

      do i=1,n
          call solve(ndim,n,a,A_(i,1:n),ipvt)
      end do

      do i=1,4
          do j=1,4
              A(i,j) = beta * E(i,j) + H(i,j)
          end do
      end do

      a_ = transpose(a_)
      print *, "inverse matrix A"
      print 101,((a_(i,j),j=1,n), i=1,n)
      r = 0.0d0
c     do i=1, n
c         do j=1, n
c             do k=1, n
c                 R(i,j) =r(i,j) + a_(j,k) * a(k,i)
c             end do
c         end do
c     end do
      r = matmul(a, a_)
c     r = transpose(r)
      do i=1, n
          do j =1, n
              R(i,j) = R(i,j) - E(i,j)
          end do
      end do

      print *, "matrix R = A*A_ - E"
      print 101,((r(i,j),j=1,n), i=1,n)
      print *, "||R|| = ", norm2(R)
      stop

  101 format(4X, 4(/1X, 4F15.4))
  102 format(1x, "beta = ", f6.4)
      end