      program tumtrx
*
      integer idim, jdim, ndim
        parameter ( idim = 1 )
        parameter ( jdim = 1 )
        parameter ( ndim = 3 )
*
      real alu ( jdim, ndim, ndim )
      real xv  ( idim, jdim, ndim )
      integer kp( jdim, ndim )
*
      integer i, j, n, m
*
      do 1000 j = 1, jdim
         do 1000 n = 1, ndim
            alu ( j, n, n   ) = 0
 1000 continue 
*
      alu( 1, 1, 1 )=1
      alu( 1, 1, 2 )=1
      alu( 1, 1, 3 )=1
*
      alu( 1, 2, 1 )=1
      alu( 1, 2, 2 )=0
      alu( 1, 2, 3 )=-1
*
      alu( 1, 3, 1 )=0
      alu( 1, 3, 2 )=1
      alu( 1, 3, 3 )=-3
*
      write (6,*) 'coefficient matrix'
      do 1050 j = 1, jdim
         write ( 6,* ) 'J = ' , j
         do 1050 n = 1, ndim
         do 1050 m = 1, ndim
            write ( 6, * ) alu( j, n, m )
 1050 continue 
*
      CALL LUMAKE
     M         ( ALU   ,
     O           KP    ,
     D           JDIM  , NDIM   )
*
      write (6,*) 'LU  matrix'
      do 1100 j = 1, jdim
         write ( 6,* ) 'J = ' , j
         do 1100 n = 1, ndim
         do 1100 m = 1, ndim
            write ( 6, * ) alu( j, n, m )
 1100 continue 
*
      call reset( xv , idim*jdim*ndim )
      xv ( 1, 1, 1 ) =  1.
      xv ( 1, 1, 2 ) =  0.
      xv ( 1, 1, 3 ) =  6.
*
      CALL LUSOLV
     M         ( XV    ,
     I           ALU   , KP    ,
     D           IDIM  , JDIM  , NDIM   )
*
      write ( 6, * ) ' solution'
      do 3000 i = 1, idim
         do 3000 j = 1, jdim
            do 3000 n = 1, ndim
               write ( 6, * ) xv ( i, j, n )
 3000 continue 
*
      stop
      end
