program lumatrix_test

  use lumatrix
  integer,parameter  :: idim=2, jdim=2, ndim=2

  real(8) :: alu3(jdim,ndim,ndim)
  real(8) :: b3(idim,jdim,ndim)
  integer :: kp2(jdim,ndim)
  real(8) :: b2(jdim,ndim)

  real(8) :: alu2(ndim,ndim)
  integer :: kp1(ndim)
  real(8) :: bb2(idim,ndim)
  real(8) :: bb1(ndim)

  alu3(1,1,:) = (/1,1/)
  alu3(1,2,:) = (/1,0/)
  alu3(2,1,:) = (/0,1/)
  alu3(2,2,:) = (/1,-1/)

  b3(1,1,:)     = (/2,0/)
  b3(2,1,:)     = (/4,0/)
  b3(1,2,:)     = (/2,0/)
  b3(2,2,:)     = (/4,0/)

  b2(1,:)     = (/2,0/)
  b2(2,:)     = (/2,0/)

  write(6,*) 'Test 323'
    call ludecomp(alu3,kp2)
    b3=lusolve(alu3,kp2,b3)

    write(6,*) b3(1,1,:)
    write(6,*) b3(2,1,:)
    write(6,*) 
    write(6,*) b3(1,2,:)
    write(6,*) b3(2,2,:)

  write(6,*) 
  write(6,*) 'Test 322'
    b2=lusolve(alu3,kp2,b2)
    write(6,*) b2(1,:)
    write(6,*) b2(2,:)

  alu2(1,:) = (/1,1/)
  alu2(2,:) = (/1,0/)

  bb2(1,:)     = (/2,0/)
  bb2(2,:)     = (/4,0/)

  bb1(:)     = (/2,0/)

  write(6,*) 
  write(6,*) 'Test 212'
    call ludecomp(alu2,kp1)
    bb2=lusolve(alu2,kp1,bb2)

    write(6,*) bb2(1,:)
    write(6,*) bb2(2,:)

  write(6,*) 
  write(6,*) 'Test 211'
    bb1=lusolve(alu2,kp1,bb1)

    write(6,*) bb1(:)

end program lumatrix_test

