function nsasm_interface ( p_file, t_file, e_file, np0, nu )

%% NSASM_INTERFACE is a sample program for NSASM.
%
%  Usage:
%
%    nsasm_test ( p_file, t_file, e_file )
%
%    where
%
%      P_FILE contains 2 rows and NP columns of (X,Y) node coordinates.
%
%      T_FILE contains 6 rows and NT columns, with each column containing
%      the (1-based) indices of nodes forming the triangles, in a particular
%      order.
%
%      E_FILE contains 3 rows and NE columns, defining constraints on
%      the data, including Dirichlet boundary values in particular.
%      Item #1 is a node, #2 is a variable index (0 = horizontal velocity,
%      1 = vertical velocity, 2 = pressure) and #3 is an associated value.
%
%      NP0 is the number of pressure nodes.
%
%      NU is the fluid viscosity.
%
%  Modified:
%
%    05 September 2006
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Per-Olof Persson,
%    Implementation of Finite Element-Based Navier-Stokes Solver,
%    April 2002
%

%  Local Parameters:
%
%    Local, integer NE, the number of constraints.
%
%    Local, integer NP, the number of nodes.
%
%    Local, integer NP0, the number of pressure nodes (midside nodes
%    of triangles.)
%
%    Local, integer NT, the number of triangles.
%
%    Local, real NU, the viscosity.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'NSASM_INTERFACE:\n' );
  fprintf ( 1, '  A simple MATLAB interface to the C routine NSASM\n' );
  fprintf ( 1, '  after it has been compiled by MATLAB''s MEX compiler.\n' );

% np0 = 9;
% nu = 100.0;


  fprintf ( 1, '\n' );
  fprintf ( 1, '  Loading user node data from "%s".\n', p_file );

% p = load ( 'p.txt' );

  p = load ( p_file );
  [ dim_num, np ] = size ( p );
  fprintf ( 1, '  Number of nodes = %d\n', np );
  fprintf ( 1, '  Number of pressure nodes NP0 = %d\n', np0 );

  fprintf ( 1, '  Loading user triangle data from "%s".\n', t_file );

  t = load ( t_file );
% t = load ( 't.txt' );
  [ triangle_order, nt ] = size ( t );
  fprintf ( 1, '  Number of triangles NT = %d\n', nt );

  fprintf ( 1, '  Loading user constraint data from "%s".\n', e_file );

  e = load ( e_file );
% e = load ( 'e.txt' );
  [ e_order, ne ] = size ( e );
  fprintf ( 1, '  Number of constraints NE = %d\n', ne );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Viscosity NU = %f\n', nu );
%
%  U is the set of finite element coefficients.
%  We set this to zero.
%  In typical usage, NSASM is used inside a Newton iteration,
%  and we are trying to improve the estimate of U.
%
  ndof = 2 * np + np0 + ne;
  fprintf ( 1, '  Degrees of freedom NDOF =  %d\n', ndof );

  u = zeros ( 1, ndof );
  
  [ K, L ] = nsasm ( p, t, np0, e, u, nu );

  l_norm = norm ( L, inf );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L-Infinity norm of L = %f\n', l_norm );
%
%  Print (some of) L.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Beginning and end of L vector:\n' );
  fprintf ( 1, '\n' );

  for i = 1 : ndof

    if ( i < 10 | ndof - 10 < i )
      fprintf ( 1, '  %8d  %14f\n', i, L(i) );
    end

    if ( i == 10 )
      fprintf ( 1, '..(skipping some entries)...\n' );
    end

  end
%
%  Get the number of nonzero entries in K.
%
  nz_num = nnz ( K );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  nz_num = nnz ( K )\n' );
  fprintf ( 1, '  Matrix nonzeros NZ_NUM = %d\n', nz_num );
%
%  Get the sparse triplet representation of K as a set of row indices,
%  column indices, and values.  Notice how the data is sorted.
%
  [ row, col, val ] = find ( K );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  [ row, col, val ] = find ( K )\n' );
  fprintf ( 1, '  Matrix sparse triplet representation:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '       ROW       COL     VAL\n' );
  fprintf ( 1, '\n' );
%
%  Print (some of) the sparse triplet data.
%
  for i = 1 : nz_num

    if ( i < 10 | nz_num - 10 < i )
      fprintf ( 1, '  %8d  %8d  %14f\n', row(i), col(i), val(i) );
    end

    if ( i == 10 )
      fprintf ( 1, '..(skipping some entries)...\n' );
    end

  end
%
%  Print (some of) the matrix.
%
  dsp_print_some ( ndof, ndof, nz_num, col, row, val, 1, 1, 10, ...
    10, '  Initial part of K as a Matrix' );