program main

c*********************************************************************72
c
cc umfpack_2.0_test() tests umfpack_2.0().
c
c  Discussion:
c
c    factor and solve a 5-by-5 system, ax=b, using default parameters,
c    except with complete printing of all arguments on input and output,
c    where
c        [ 2  3  0  0  0 ]      [  8 ]                  [ 1 ]
c        [ 3  0  4  0  6 ]      [ 45 ]                  [ 2 ]
c    a = [ 0 -1 -3  2  0 ], b = [ -3 ]. solution is x = [ 3 ].
c        [ 0  0  1  0  0 ]      [  3 ]                  [ 4 ]
c        [ 0  4  2  0  1 ]      [ 19 ]                  [ 5 ]
c
c    solve a'x=b, with solution:
c      x = [  1.8158  1.4561 1.5000 -24.8509 10.2632 ]'
c    using the factors of a.
c
c    modify one entry (a (5,2) = 1.0e-6) and refactorize.  
c
c    solve ax=b both without and with iterative refinement,
c    with true solution (rounded to 4 digits past the decimal point):
c       x = [-15.0000 12.6667 3.0000   9.3333 13.0000 ]'
c
c  Modified:
c
c    03 November 2013
c
      integer nmax, nemax, lvalue, lindex
      parameter (nmax=20, nemax=100, lvalue=300, lindex=300)
      integer keep (20), index (lindex), info (40),
     &  i, icntl (20), n, ne, ai (2*nemax)
      real b (nmax), x (nmax), w (4*nmax), value (lvalue),
     & cntl (10), rinfo (20), ax (nemax)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'umfpack_2.0_test():'
      write ( *, '(a)' ) '  Test umfpack_2.0().'
c
c  read input matrix and right-hand side.  
c
c  keep a copy of the triplet form in ai and ax.
c
      n = 5
      ne = 12
      ai(1:2*ne) = (/
     &  1,1,2,2,2,3,3,3,4,5,5,5,
     &  1,2,1,3,5,2,3,4,3,2,3,5 /)
      ax(1:ne) = (/ 
     &  2.0, 3.0, 3.0, 4.0, 6.0, -1.0, -3.0,
     &  2.0, 1.0, 4.0, 2.0, 1.0 /)
      b(1:n) = (/ 8.0, 45.0, -3.0, 3.0, 19.0 /)

      do i = 1, ne
        index (i) = ai (i)
        index (ne+i) = ai (ne+i)
        value (i) = ax (i)
      end do
c
c initialize controls, and change default printing control.  note that
c this change from the default should only be used for test cases.  it
c can generate a lot of output for large matrices.
c
      call ums2in (icntl, cntl, keep)
      icntl (3) = 4
c
c  Factorize a, and print the factors.  input matrix is not preserved.
c
      call ums2fa (n, ne, 0, .false., lvalue, lindex, value, index,
     &               keep, cntl, icntl, info, rinfo)
      if (info (1) .lt. 0) then
        stop
      end if
c
c  Reset default printing control (ums2in could be called instead)
c
      icntl (3) = 2
c
c  Solve ax = b and print solution.
c
      call ums2so (n, 0, .false., lvalue, lindex, value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution of A * x = b:'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)

      if (info (1) .lt. 0) then
        stop
      end if
c
c  Solve a'x = b and print solution.
c
      call ums2so (n, 0, .true., lvalue, lindex,  value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)


      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution of A'' * x = b:'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)

      if (info (1) .lt. 0) then
        stop
      end if
c
c  modify one entry of a, and refactorize using ums2rf.
c
      do i = 1, ne
        index (i) = ai (i)
        index (ne+i) = ai (ne+i)
        value (i) = ax (i)
      end do
c
c  A(5,2) happens to be (paq)_22, the second pivot entry:
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Modify one entry of A, call it A#:'
      write ( *, '(a)' ) ' '

      value (10) = 1.0e-6
 
      call ums2rf (n, ne, 1, .false., lvalue, lindex, value, index,
     &               keep, cntl, icntl, info, rinfo)

      if (info (1) .lt. 0) then
        stop
      end if
c
c  Solve ax = b without iterative refinement, and print solution.
c  This will be very inaccurate due to the tiny second pivot entry.
 
      call ums2so (n, 0, .false., lvalue, lindex,  value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution of A# * x = b:'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)

      if (info (1) .lt. 0) then
        stop
      end if
c
c  Solve ax = b with iterative refinement, and print solution.
c  This is much more accurate.
 
      icntl (8) = 10
      call ums2so (n, 0, .false., lvalue, lindex,  value, index,
     &               keep, b, x, w, cntl, icntl, info, rinfo)
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Solution of A# * x = b'
      write ( *, '(a)' ) '  via iterative refinement.'
      write ( *, '(a)' ) ' '
      write (*, 30) (x (i), i = 1, n)

      if (info (1) .lt. 0) then
        stop
      end if

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'umfpack_2.0_test():'
      write ( *, '(a)' ) '  Normal end of execution'

      stop
30    format ('solution: ', 5(/,f20.16))
      end