text:
function hb_to_st ( input_file, output_file )

%*****************************************************************************80
%
%% HB_TO_ST reads a Harwell-Boeing Sparse Matrix File, and writes an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE, the name of the file containing the
%    Harwell-Boeing matrix data.
%
%    Input, string OUTPUT_FILE, the name of the ST file to be created.
%
  input_unit = fopen ( input_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_TO_ST - Fatal error!\n' );
    fprintf ( 1, '  Error opening the file.\n' );
    error ( 'HB_TO_MSM - Fatal error!' );
  end

  [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
    nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
    nrhs, nrhsix, colptr, rowind, values, rhsval, rhsptr, rhsind, guess, ...
    exact ] = hb_file_read ( input_unit );

  for col = 1 : ncol
    for k = colptr(col) : colptr(col+1)-1
      colind(k) = col;
    end
  end

  fclose ( input_unit );

  base0 = 0;
  base1 = 1;
  rowind = st_rebase ( base1, base0, nnzero, rowind );
  colind = st_rebase ( base1, base0, nnzero, colind );

  st_header_print ( nrow, ncol, nnzero );

  st_header_write ( output_file, nrow, ncol, nnzero );

  st_data_print ( nrow, ncol, nnzero, rowind, colind, values, ...
    '  Data read from HB file' );

  st_data_write ( output_file, nrow, ncol, nnzero, rowind, colind, values );

  return
end
function c = ch_cap ( c )

%*****************************************************************************80
%
%% CH_CAP capitalizes a single character.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to capitalize.
%
%    Output, character C, the capitalized character.
%
  if ( 'a' <= c & c <= 'z' )
    c = c + 'A' - 'a';
  end

  return
end
function truefalse = ch_eqi ( c1, c2 )

%*****************************************************************************80
%
%% CH_EQI is a case insensitive comparison of two characters for equality.
%
%  Example:
%
%    CH_EQI ( 'A', 'a' ) is TRUE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 July 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C1, C2, the characters to compare.
%
%    Output, logical TRUEFALSE, is TRUE (1) if the characters are equal.
%
  FALSE = 0;
  TRUE = 1;

  if ( ch_cap ( c1 ) == ch_cap ( c2 ) )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function truefalse = ch_is_digit ( c )

%*****************************************************************************80
%
% CH_IS_DIGIT returns TRUE if the character C is a digit.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 December 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, a character.
%
%    Output, integer TRUEFALSE, is TRUE (1) if C is a digit, FALSE (0) otherwise.
%
  TRUE = 1;
  FALSE = 0;

  if ( '0' <= c & c <= '9' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function truefalse = ch_is_format_code ( c )

%*****************************************************************************80
%
%% CH_IS_FORMAT_CODE returns TRUE if a character is a FORTRAN format code.
%
%  Discussion:
%
%    The format codes accepted here are not the only legal format
%    codes in FORTRAN90.  However, they are more than sufficient
%    for my needs!
%
%  Table:
%
%    A  Character
%    B  Binary digits
%    D  Real number, exponential representation
%    E  Real number, exponential representation
%    F  Real number, fixed point
%    G  General format
%    I  Integer
%    L  Logical variable
%    O  Octal digits
%    Z  Hexadecimal digits
%    *  Free format
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the character to be analyzed.
%
%    Output, logical TRUEFALSE, is TRUE (1) if C is a FORTRAN format code
%     and FALSE (0) otherwise.
%
  FALSE = 0;
  TRUE = 1;

      if ( ch_eqi ( c, 'A' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'B' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'D' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'E' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'F' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'G' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'I' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'L' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'O' ) )
    truefalse = TRUE;
  elseif ( ch_eqi ( c, 'Z' ) )
    truefalse = TRUE;
  elseif ( c == '*' )
    truefalse = TRUE;
  else
    truefalse = FALSE;
  end

  return
end
function digit = ch_to_digit ( c )

%*****************************************************************************80
%
%% CH_TO_DIGIT returns the integer value of a base 10 digit.
%
%  Example:
%
%     C   DIGIT
%    ---  -----
%    '0'    0
%    '1'    1
%    ...  ...
%    '9'    9
%    ' '    0
%    'X'   -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, character C, the decimal digit, '0' through '9' or blank
%    are legal.
%
%    Output, integer DIGIT, the corresponding integer value.  If C was
%    'illegal', then DIGIT is -1.
%
  if ( '0' <= c & c <= '9' )

    digit = c - '0';

  elseif ( c == ' ' )

    digit = 0;

  else

    digit = -1;

  end

  return
end
function exact = hb_exact_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp )

%*****************************************************************************80
%
%% HB_EXACT_READ reads the exact solution vectors in an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer RHSCRD, the number of lines in the file for
%    right hand sides.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Output, real EXACT(NROW,NRHS), the exact solution vectors.
%
  exact = [];

  if ( 0 < rhscrd )

    if ( rhstyp(3) == 'X' )

      [ p, code, w, m ] = s_to_format ( rhsfmt );

      line_num = 1 + floor ( ( nrow - 1 ) / p );

      for irhs = 1 : nrhs

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrow );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            exact(j,irhs) = s_to_r8 ( s );
          end
        end
      end

    end

  end

  return
end
function [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
  nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
  nrhs, nrhsix, colptr, rowind, values, rhsval, rhsptr, rhsind, guess, ...
  exact ] = hb_file_read ( input_unit )

%*****************************************************************************80
%
%% HB_FILE_READ reads an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from data is read.
%
%    Output, character ( len = 72 ) TITLE, a title for the matrix.
%
%    Output, character ( len = 8 ) KEY, an identifier for the matrix.
%
%    Output, integer TOTCRD, the total number of lines of data.
%
%    Output, integer PTRCRD, the number of input lines for pointers.
%
%    Output, integer INDCRD, the number of input lines for row indices.
%
%    Output, integer VALCRD, the number of input lines for numerical values.
%
%    Output, integer RHSCRD, the number of input lines for right hand sides.
%
%    Output, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Output, integer NROW, the number of rows or variables.
%
%    Output, integer NCOL, the number of columns or elements.
%
%    Output, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Output, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Output, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Output, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Output, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Output, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Output, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Output, integer NRHS, the number of right hand sides.
%
%    Output, integer NRHSIX, the number of row indices (set to 0
%    in the case of unassembled matrices.)  Ignored if NRHS = 0.
%
%    Output, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Output, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%
%    Output, real VALUES(NNZERO) or VALUES(NELTVL), the nonzero values
%    of the matrix.
%
%    If RHSTYP(1:1) == 'F':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NROW,NRHS), contains NRHS dense right hand
%      side vectors.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'A':
%
%      Output, integer RHSPTR(NRHS+1), RHSPTR(I) points to the location of
%      the first entry of right hand side I in the sparse right hand
%      side vector.
%
%      Output, integer RHSIND(NRHSIX), indicates, for each entry of
%      RHSVAL, the corresponding row index.
%
%      Output, real RHSVAL(NRHSIX), contains the value of the right hand
%      side entries.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'E':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NNZERO,NRHS), contains NRHS unassembled
%      finite element vector right hand sides.
%
%    Output, real GUESS(NROW,NRHS), the starting guess vectors.
%
%    Output, real EXACT(NROW,NRHS), the exact solution vectors.
%

%
%  Read the header block.
%
  [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
    nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
    nrhs, nrhsix ] = hb_header_read ( input_unit );
%
%  Read the matrix structure.
%
  [ colptr, rowind ] = hb_structure_read ( input_unit, ncol, mxtype, ...
    nnzero, neltvl, ptrcrd, ptrfmt, indcrd, indfmt );
%
%  Read the matrix values.
%
  values = hb_values_read ( input_unit, valcrd, mxtype, nnzero, neltvl, ...
    valfmt );
%
%  Read the right hand sides.
%
  [ rhsval, rhsptr, rhsind ] = hb_rhs_read ( input_unit, nrow, nnzero, ...
    nrhs, nrhsix, rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp );
%
%  Read the starting guesses.
%
  guess = hb_guess_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp );
%
%  Read the exact solutions.
%
  exact = hb_exact_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp );

  return
end
function guess = hb_guess_read ( input_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp )

%*****************************************************************************80
%
%% HB_GUESS_READ reads the starting guess vectors in an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer RHSCRD, the number of lines in the file for
%    right hand sides.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    Output, real GUESS(NROW,NRHS), the starting guess vectors.
%
  guess = [];

  if ( 0 < rhscrd )

    if ( rhstyp(2) == 'G' )

      [ p, code, w, m ] = s_to_format ( rhsfmt );

      line_num = 1 + floor ( ( nrow - 1 ) / p );

      for irhs = 1 : nrhs

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrow );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            guess(j,irhs) = s_to_r8 ( s );
          end
        end
      end

    end

  end

  return
end
function [ title, key, totcrd, ptrcrd, indcrd, valcrd, rhscrd, mxtype, ...
  nrow, ncol, nnzero, neltvl, ptrfmt, indfmt, valfmt, rhsfmt, rhstyp, ...
  nrhs, nrhsix ] = hb_header_read ( input_unit )

%*****************************************************************************80
%
%% HB_HEADER_READ reads the header of an HB file.
%
%  Discussion:
%
%    The user should already have opened the file, and positioned it
%    to the first record.
%
%    Thanks to Shahadat Hossain for pointing out an error in the determination
%    of INDFMT, on 24 June 2004.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 June 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Output, character ( len = 72 ) TITLE, a title for the matrix.
%
%    Output, character ( len = 8 ) KEY, an identifier for the matrix.
%
%    Output, integer TOTCRD, the total number of lines of data.
%
%    Output, integer PTRCRD, the number of input lines for pointers.
%
%    Output, integer INDCRD, the number of input lines for row indices.
%
%    Output, integer VALCRD, the number of input lines for numerical values.
%
%    Output, integer RHSCRD, the number of input lines for right hand sides.
%
%    Output, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Output, integer NROW, the number of rows or variables.
%
%    Output, integer NCOL, the number of columns or elements.
%
%    Output, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Output, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Output, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Output, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Output, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Output, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Output, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%
%    Output, integer NRHS, the number of right hand sides.
%
%    Output, integer NRHSIX, the number of entries of storage for right
%    hand side values, in the case where RHSTYP(1:1) = 'M' and
%    MXTYPE(3:3) = 'A'.
%

%
%  Read the header block.
%  Use FGETL rather that FGETS, because we don't want the line terminator character!
%  If fewer than 80 characters were read, you need to pad the line out.
%
  line = fgetl ( input_unit );

  if ( line == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  I/O error reading header line 1.\n' );
    return;
  end

  line_num = length ( line );
  for i = line_num+1 : 80
    line(i) = ' ';
  end

  title = line(1:72);
  title = title(1:s_len_trim(title));

  key = line(73:80);
  key = key(1:s_len_trim(key));

  line = fgetl ( input_unit );

  if ( line == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  I/O error reading header line 2.\n' );
    return;
  end

  line_num = length ( line );
  for i = line_num+1 : 80
    line(i) = ' ';
  end

  totcrd = s_to_i4 ( line( 1:14) );
  ptrcrd = s_to_i4 ( line(15:28) );
  indcrd = s_to_i4 ( line(29:42) );
  valcrd = s_to_i4 ( line(43:56) );
  rhscrd = s_to_i4 ( line(57:70) );

  line = fgetl ( input_unit );

  if ( line == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  I/O error reading header line 3.\n' );
    return;
  end

  line_num = length ( line );
  for i = line_num+1 : 80
    line(i) = ' ';
  end

  mxtype =           line( 1: 3);
  nrow   = s_to_i4 ( line(15:28) );
  ncol   = s_to_i4 ( line(29:42) );
  nnzero = s_to_i4 ( line(43:56) );
  neltvl = s_to_i4 ( line(57:70) );

  line = fgetl ( input_unit );

  if ( line == -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_HEADER_READ - Fatal error!\n' );
    fprintf ( 1, '  I/O error reading header line 4.\n' );
    return;
  end

  line_num = length ( line );
  for i = line_num+1 : 80
    line(i) = ' ';
  end

  ptrfmt = line( 1:16);
  ptrfmt = ptrfmt(1:s_len_trim(ptrfmt));
  indfmt = line(17:32);
  indfmt = indfmt(1:s_len_trim(indfmt));
  valfmt = line(33:52);
  valfmt = valfmt(1:s_len_trim(valfmt));
  rhsfmt = line(53:72);
  rhsfmt = rhsfmt(1:s_len_trim(rhsfmt));

  if ( 0 < rhscrd )

    line = fgetl ( input_unit );

    if ( line == -1 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'HB_HEADER_READ - Fatal error!\n' );
      fprintf ( 1, '  I/O error reading header line 5.\n' );
      return;
    end

    line_num = length ( line );
    for i = line_num+1 : 80
      line(i) = ' ';
    end

    rhstyp =           line( 1: 3);
    nrhs   = s_to_i4 ( line(15:28) );
    nrhsix = s_to_i4 ( line(29:42) );

  else

    rhstyp = ' ';
    nrhs = 0;
    nrhsix = 0;

  end

  return
end
function [ rhsval, rhsptr, rhsind ] = hb_rhs_read ( input_unit, nrow, ...
  nnzero, nrhs, nrhsix, rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp )

%*****************************************************************************80
%
%% HB_RHS_READ reads the right hand side information in an HB file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NROW, the number of rows or variables.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NRHS, the number of right hand sides.
%
%    Input, integer NRHSIX, the number of entries of storage for right
%    hand side values, in the case where RHSTYP(1:1) = 'M' and
%    MXTYPE(3:3) = 'A'.
%
%    Input, integer RHSCRD, the number of lines in the file for
%    right hand sides.
%
%    Input, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Input, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Input, character ( len = 20 ) RHSFMT, the format for reading values
%    of the right hand side.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, character ( len = 3 ) RHSTYP, the right hand side type.
%    First character is F for full storage or M for same as matrix.
%    Second character is G if starting "guess" vectors are supplied.
%    Third character is X if exact solution vectors are supplied.
%    Ignored if NRHS = 0.
%
%    If RHSTYP(1:1) == 'F':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NROW,NRHS), contains NRHS dense right hand
%      side vectors.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'A':
%
%      Output, integer RHSPTR(NRHS+1), RHSPTR(I) points to the location of
%      the first entry of right hand side I in the sparse right hand
%      side vector.
%
%      Output, integer RHSIND(NRHSIX), indicates, for each entry of
%      RHSVAL, the corresponding row index.
%
%      Output, real RHSVAL(NRHSIX), contains the value of the right hand
%      side entries.
%
%    If RHSTYP(1:1) = 'M' and MXTYPE(3:3) = 'E':
%
%      Output, integer RHSPTR(*), is not used.
%
%      Output, integer RHSIND(*), is not used.
%
%      Output, real RHSVAL(NNZERO,NRHS), contains NRHS unassembled
%      finite element vector right hand sides.
%
  rhsptr = [];
  rhsind = [];
  rhsval = [];
%
%  Read the right hand sides.
%    case F                             = "full" or "dense";
%    case not F + matrix storage is "A" = sparse pointer RHS
%    case not F + matrix storage is "E" = finite element RHS
%
  if ( 0 < rhscrd )
%
%  Dense right hand sides:
%
    if ( rhstyp(1) == 'F' )

      [ p, code, w, m ] = s_to_format ( rhsfmt );

      line_num = 1 + floor ( ( nrow - 1 ) / p );

      for irhs = 1 : nrhs

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrow );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            rhsval(j,irhs) = s_to_r8 ( s );
          end
        end
      end
%
%  Sparse right-hand sides stored like the matrix.
%  Read pointer array, indices, and values.
%
    elseif ( rhstyp(1) == 'M' )

      if ( mxtype(3) == 'A' )

        [ p, code, w, m ] = s_to_format ( ptrfmt );

        line_num = 1 + floor ( ( nrhs + 1 - 1 ) / p );

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrhs + 1 );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            rhsptr(j) = s_to_i4 ( s );
          end
        end

        [ p, code, w, m ] = s_to_format ( indfmt );

        line_num = 1 + floor ( ( nrhsix - 1 ) / p );

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrhsix );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            rhsind(j) = s_to_i4 ( s );
          end
        end

        [ p, code, w, m ] = s_to_format ( rhsfmt );

        line_num = 1 + floor ( ( nrhsix - 1 ) / p );

        jhi = 0;

        for i = 1 : line_num

          line = fgetl ( input_unit );
          jlo = jhi + 1;
          jhi = min ( jlo + p - 1, nrhsix );

          khi = 0;

          for j = jlo : jhi
            klo = khi + 1;
            khi = min ( klo + w - 1, length ( line ) );
            s = line(klo:khi);
            rhsval(j) = s_to_r8 ( s );
          end
        end
%
%  Sparse right hand sides in finite element format.
%
      elseif ( mxtype(3) == 'E' )

        [ p, code, w, m ] = s_to_format ( rhsfmt );

        line_num = 1 + floor ( ( nnzero - 1 ) / p );

        for irhs = 1 : nrhs

          jhi = 0;

          for i = 1 : line_num

            line = fgetl ( input_unit );
            jlo = jhi + 1;
            jhi = min ( jlo + p - 1, nnzero );

            khi = 0;

            for j = jlo : jhi
              klo = khi + 1;
              khi = min ( klo + w - 1, length ( line ) );
              s = line(klo:khi);
              rhsval(j,irhs) = s_to_r8 ( s );
            end
          end
        end

      else

        fprintf ( 1, '\n' );
        fprintf ( 1, 'HB_RHS_READ - Fatal error!\n' );
        fprintf ( 1, '  Illegal value of MXTYPE(3)!\n' );
        error ( 'HB_RHS_READ - Fatal error!' );

      end
%
%  0 < RHSCARD, but RHSTYP not recognized.
%
    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'HB_RHS_READ - Fatal error!\n' );
      fprintf ( 1, '  Illegal value of RHSTYP(1)!\n' );
      error ( 'HB_RHS_READ - Fatal error!' );

    end

  end

  return
end
function [ colptr, rowind ] = hb_structure_read ( input_unit, ncol, mxtype, ...
  nnzero, neltvl, ptrcrd, ptrfmt, indcrd, indfmt )

%*****************************************************************************80
%
%% HB_STRUCTURE_READ reads the structure of an HB matrix.
%
%  Discussion:
%
%    The user should already have opened the file, and positioned it
%    to just after the header records.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 April 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer NCOL, the number of columns.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, integer PTRCRD, the number of pointer records.
%
%    Input, character ( len = 16 ) PTRFMT, the format for reading pointers.
%
%    Input, integer INDCRD, the number of index records.
%
%    Input, character ( len = 16 ) INDFMT, the format for reading indices.
%
%    Output, integer COLPTR(NCOL+1), COLPTR(I) points to the location of
%    the first entry of column I in the sparse matrix structure.
%
%    Output, integer ROWIND(NNZERO) or ROWIND(NELTVL), the row index of
%    each item.
%
  [ p, code, w, m ] = s_to_format ( ptrfmt );

  if ( mxtype(3) == 'A' )
    line_num = 1 + floor ( ( ( ncol + 1 ) - 1 ) / p );
  else
    line_num = 1 + floor ( ( ( ncol     ) - 1 ) / p );
  end

  jhi = 0;

  for i = 1 : line_num
    line = fgetl ( input_unit );
    jlo = jhi + 1;
    if ( mxtype(3) == 'A' )
      jhi = min ( jlo + p - 1, ncol + 1 );
    else
      jhi = min (  jlo + p - 1, ncol );
    end

    khi = 0;

    for j = jlo : jhi
      klo = khi + 1;
      khi = min ( klo + w - 1, length ( line ) );
      s = line(klo:khi);
      colptr(j) = s_to_i4 ( s );
    end

  end

  if ( mxtype(3) == 'A' )

    [ p, code, w, m ] = s_to_format ( indfmt );

    line_num = 1 + floor ( ( nnzero - 1 ) / p );

    jhi = 0;

    for i = 1 : line_num
      line = fgetl ( input_unit );
      jlo = jhi + 1;
      jhi = min ( jlo + p - 1, nnzero );

      khi = 0;

      for j = jlo : jhi
        klo = khi + 1;
        khi = min ( klo + w - 1, length ( line ) );
        s = line(klo:khi);
        rowind(j) = s_to_i4 ( s );
      end

    end

  elseif ( mxtype(3) == 'E' )

    [ p, code, w, m ] = s_to_format ( indfmt );
    number = colptr(ncol) - colptr(1);
    line_num = 1 + floor ( ( number - 1 ) / p );

    jhi = 0;

    for i = 1 : line_num
      line = fgetl ( input_unit );
      jlo = jhi + 1;
      jhi = min ( jlo + p - 1, number );

      khi = 0;

      for j = jlo : jhi
        klo = khi + 1;
        khi = min ( klo + w - 1, length ( line ) );
        s = line(klo:khi);
        rowind(j) = s_to_i4 ( s );
      end

    end

  else

    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_STRUCTURE_READ - Fatal error!\n' );
    fprintf ( 1, '  Illegal value of MXTYPE(3).\n' );
    return

  end

  return
end
function values = hb_values_read ( input_unit, valcrd, mxtype, nnzero, ...
  neltvl, valfmt )

%*****************************************************************************80
%
%% HB_VALUES_READ reads the values of an HB matrix.
%
%  Discussion:
%
%    The user should already have opened the file, and positioned it
%    to just after the header and structure records.
%
%    Values are contained in an HB file if the VALCRD parameter
%    is nonzero.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Iain Duff, Roger Grimes, John Lewis,
%    User's Guide for the Harwell-Boeing Sparse Matrix Collection,
%    October 1992.
%
%  Parameters:
%
%    Input, integer INPUT_UNIT, the unit from which data is read.
%
%    Input, integer VALCRD, the number of input lines for numerical values.
%
%    Input, character ( len = 3 ) MXTYPE, the matrix type.
%    First character is R for Real, C for complex, P for pattern only.
%    Second character is S for symmetric, U for unsymmetric, H for
%      Hermitian, Z for skew symmetric, R for rectangular.
%    Third character is A for assembled and E for unassembled
%      finite element matrices.
%
%    Input, integer NNZERO.  In the case of assembled sparse matrices,
%    this is the number of nonzeroes.  In the case of unassembled finite
%    element matrices, in which the right hand side vectors are also
%    stored as unassembled finite element vectors, this is the total
%    number of entries in a single unassembled right hand side vector.
%
%    Input, integer NELTVL, the number of finite element matrix entries,
%    set to 0 in the case of assembled matrices.
%
%    Input, character ( len = 20 ) VALFMT, the format for reading values.
%
%    Output, real VALUES(NNZERO) or VALUES(NELTVL), the nonzero values
%    of the matrix.
%
  [ p, code, w, m ] = s_to_format ( valfmt );
%
%  Read the matrix values.
%    case "A" = assembled;
%    case "E" = unassembled finite element matrices.
%
  if ( 0 < valcrd )

    if ( mxtype(3:3) == 'A' )

      line_num = 1 + floor ( ( nnzero - 1 ) / p );

      jhi = 0;

      for i = 1 : line_num

        line = fgetl ( input_unit );
        jlo = jhi + 1;
        jhi = min ( jlo + p - 1, nnzero );

        khi = 0;

        for j = jlo : jhi
          klo = khi + 1;
          khi = min ( klo + w - 1, length ( line ) );
          s = line(klo:khi);
          values(j) = s_to_r8 ( s );
        end

      end

    elseif ( mxtype(3) == 'E' )

      line_num = 1 + floor ( ( neltvl - 1 ) / p );

      jhi = 0;

      for i = 1 : line_num

        line = fgetl ( input_unit );
        jlo = jhi + 1;
        jhi = min ( jlo + p - 1, neltvl );

        khi = 0;

        for j = jlo : jhi
          klo = khi + 1;
          khi = min ( klo + w - 1, length ( line ) );
          s = line(klo:khi);
          values(j) = s_to_r8 ( s );
        end

      end

    else

      fprintf ( 1, '\n' );
      fprintf ( 1, 'HB_VALUES_READ - Fatal error!\n' );
      fprintf ( 1, '  Illegal value of MXTYPE(3).\n' );
      error ( 'HB_VALUES_READ - Fatal error!' );

    end

  end

  return
end
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a character string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end
function [ r, code, w, m ] = s_to_format ( s )

%*****************************************************************************80
%
%% S_TO_FORMAT reads a FORTRAN format from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the format.  This routine is limited in its ability to
%    recognize FORTRAN formats.  In particular, we are only expecting
%    a single format specification, and cannot handle extra features
%    such as 'ES' and 'EN' codes, '5X' spacing, and so on.
%
%    Legal input is:
%
%       0 nothing
%       1 blanks
%       2 optional '('
%       3 blanks
%       4 optional repeat factor R
%       5 blanks
%       6 CODE ( 'A', 'B', 'E', 'F', 'G', 'I', 'L', 'O', 'Z', '*' )
%       7 blanks
%       8 width W
%       9 optional decimal point
%      10 optional mantissa M
%      11 blanks
%      12 optional ')'
%      13 blanks
%
%  Example:
%
%    S                 R   CODE   W    M
%
%    'I12              1   I      12   0
%    'E8.0'            1   E       8   0
%    'F10.5'           1   F      10   5
%    '2G14.6'          2   G      14   6
%    '*'               1   *      -1  -1
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read.
%
%    Output, integer R, the repetition factor, which defaults to 1.
%
%    Output, character CODE, the format code.
%
%    Output, integer W, the field width.
%
%    Output, integer M, the mantissa width.
%
  LEFT = 1;
  RIGHT = -1;

  state = 0;
  paren_sum = 0;
  pos = 0;
  s_length = s_len_trim ( s );

  r = 0;
  w = 0;
  code = '?';
  m = 0;

  while ( pos < s_length )

    pos = pos + 1;
    c = s(pos);
%
%  BLANK character:
%
    if ( c == ' ' )

      if ( state == 4 )
        state = 5;
      elseif ( state == 6 )
        state = 7;
      elseif ( state == 10 )
        state = 11;
      elseif ( state == 12 )
        state = 13;
      end
%
%  LEFT PAREN
%
    elseif ( c == '(' )

      if ( state < 2 )
        paren_sum = paren_sum + LEFT;
      else
        state = -1;
        break;
      end
%
%  DIGIT (R, F, or W)
%
    elseif ( ch_is_digit ( c ) )

      if ( state <= 3 )
        state = 4;
        r = ch_to_digit ( c );
      elseif ( state == 4 )
        d = ch_to_digit ( c );
        r = 10 * r + d;
      elseif ( state == 6 | state == 7 )
        if ( code == '*' )
          state = -1;
          break;
        end
        state = 8;
        w = ch_to_digit ( c );
      elseif ( state == 8 )
        d = ch_to_digit ( c );
        w = 10 * w + d;
      elseif ( state == 9 )
        state = 10;
        m = ch_to_digit ( c );
      elseif ( state == 10 )
        d = ch_to_digit ( c );
        m = 10 * m + d;
      else
        state = -1;
        break;
      end
%
%  DECIMAL POINT
%
    elseif ( c == '.' )

      if ( state == 8 )
        state = 9;
      else
        state = -1;
        break;
      end
%
%  RIGHT PAREN
%
    elseif ( c == ')' )

      paren_sum = paren_sum + RIGHT;

      if ( paren_sum ~= 0 )
        state = -1;
        break;
      end

      if ( state == 6 & code == '*' )
        state = 12;
      elseif ( 6 <= state )
        state = 12;
      else
        state = -1;
        break;
      end
%
%  Code
%
    elseif ( ch_is_format_code ( c ) )

      if ( state < 6 )
        state = 6;
        code = c;
      else
        state = -1;
        break;
      end
%
%  Unexpected character
%
    else

      state = -1;
      break;

    end

  end

  if ( paren_sum ~= 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parentheses mismatch.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( state < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_FORMAT - Fatal error!\n' );
    fprintf ( 1, '  Parsing error.\n' );
    error ( 'S_TO_FORMAT - Fatal error!' );
  end

  if ( r == 0 )
    r = 1;
  end

  return
end
function ival = s_to_i4 ( s )

%*****************************************************************************80
%
%% S_TO_I4 reads an integer value from a string.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, a string to be examined.
%
%    Output, integer IVAL, the integer value read from the string.
%
  sgn = 1;
  state = 0;
  ival = 0;

  i = 0;

  while ( i < s_len_trim ( s ) )

    i = i + 1;
    c = s(i);

    if ( state == 0 )

      if ( c == ' ' )

      elseif ( c == '-' )
        state = 1;
        sgn = -1;
      elseif ( c == '+' )
        state = 1;
        sgn = +1;
      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return;
      end
%
%  Have read the sign, now expecting the first digit.
%
    elseif ( state == 1 )

      if ( c == ' ' )

      elseif ( '0' <= c & c <= '9' )
        state = 2;
        ival = c - '0';
      else
        fprintf ( '\n' );
        fprintf ( 'S_TO_I4 - Fatal error!\n' );
        fprintf ( '  Illegal character %c while in state %d.\n', c, state );
        return
      end
%
%  Have read at least one digit, expecting more.
%
    elseif ( state == 2 )

      if ( '0' <= c & c <= '9' )
        ival = 10 * ival + c - '0';
      else
        ival = sgn * ival;
        return;
      end

    end

  end
%
%  If we read all the characters in the string, see if we're OK.
%
  if ( state ~= 2 )
    fprintf ( '\n' );
    fprintf ( 'S_TO_I4 - Fatal error!\n' );
    fprintf ( '  Did not read enough information to define an integer!\n' );
    return;
  end

  ival = sgn * ival;

  return
end
function [ r, lchar, ierror ] = s_to_r8 ( s )

%*****************************************************************************80
%
%% S_TO_R8 reads an R8 from a string.
%
%  Discussion:
%
%    This routine will read as many characters as possible until it reaches
%    the end of the string, or encounters a character which cannot be
%    part of the real number.
%
%    Legal input is:
%
%       1 blanks,
%       2 '+' or '-' sign,
%       2.5 spaces
%       3 integer part,
%       4 decimal point,
%       5 fraction part,
%       6 'E' or 'e' or 'D' or 'd', exponent marker,
%       7 exponent sign,
%       8 exponent integer part,
%       9 exponent decimal point,
%      10 exponent fraction part,
%      11 blanks,
%      12 final comma or semicolon.
%
%    with most quantities optional.
%
%  Example:
%
%    S                 R
%
%    '1'               1.0
%    '     1   '       1.0
%    '1A'              1.0
%    '12,34,56'        12.0
%    '  34 7'          34.0
%    '-1E2ABCD'        -100.0
%    '-1X2ABCD'        -1.0
%    ' 2E-1'           0.2
%    '23.45'           23.45
%    '-4.2E+2'         -420.0
%    '17d2'            1700.0
%    '-14e-2'         -0.14
%    'e2'              100.0
%    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    22 November 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string containing the
%    data to be read.  Reading will begin at position 1 and
%    terminate at the end of the string, or when no more
%    characters can be read to form a legal real.  Blanks,
%    commas, or other nonnumeric data will, in particular,
%    cause the conversion to halt.
%
%    Output, real R, the value that was read from the string.
%
%    Output, integer LCHAR, the number of characters of S that were used to form R.
%
%    Output, integer IERROR, is 0 if no error occurred.
%
  s_length = s_len_trim ( s );
  ierror = 0;
  r = 0.0;
  lchar = -1;
  isgn = 1;
  rtop = 0.0;
  rbot = 1.0;
  jsgn = 1;
  jtop = 0;
  jbot = 1;
  ihave = 1;
  iterm = 0;

  while ( 1 )

    lchar = lchar + 1;
    c = s(lchar+1);
%
%  Blank character.
%
    if ( c == ' ' )

      if ( ihave == 2 )

      elseif ( ihave == 6 | ihave == 7 )
        iterm = 1;
      elseif ( 1 < ihave )
        ihave = 11;
      end
%
%  Comma.
%
    elseif ( c == ',' | c == ';' )

      if ( ihave ~= 1 )
        iterm = 1;
        ihave = 12;
        lchar = lchar + 1;
      end
%
%  Minus sign.
%
    elseif ( c == '-' )

      if ( ihave == 1 );
        ihave = 2;
        isgn = -1;
      elseif ( ihave == 6 )
        ihave = 7;
        jsgn = -1;
      else
        iterm = 1;
      end
%
%  Plus sign.
%
    elseif ( c == '+' )

      if ( ihave == 1 )
        ihave = 2;
      elseif ( ihave == 6 )
        ihave = 7;
      else
        iterm = 1;
      end
%
%  Decimal point.
%
    elseif ( c == '.' )

      if ( ihave < 4 )
        ihave = 4;
      elseif ( 6 <= ihave & ihave <= 8 )
        ihave = 9;
      else
        iterm = 1;
      end
%
%  Exponent marker.
%
    elseif ( ch_eqi ( c, 'E' ) | ch_eqi ( c, 'D' ) )

      if ( ihave < 6 )
        ihave = 6;
      else
        iterm = 1;
      end
%
%  Digit.
%
    elseif ( ihave < 11 & ch_is_digit ( c ) )

      if ( ihave <= 2 )
        ihave = 3;
      elseif ( ihave == 4 )
        ihave = 5;
      elseif ( ihave == 6 | ihave == 7 )
        ihave = 8;
      elseif ( ihave == 9 )
        ihave = 10;
      end

      d = ch_to_digit ( c );

      if ( ihave == 3 )
        rtop = 10.0 * rtop + d;
      elseif ( ihave == 5 )
        rtop = 10.0 * rtop + d;
        rbot = 10.0 * rbot;
      elseif ( ihave == 8 )
        jtop = 10 * jtop + d;
      elseif ( ihave == 10 )
        jtop = 10 * jtop + d;
        jbot = 10 * jbot;
      end
%
%  Anything else is regarded as a terminator.
%
    else
      iterm = 1;
    end
%
%  If we haven't seen a terminator, and we haven't examined the
%  entire string, go get the next character.
%
    if ( iterm == 1 | s_length <= lchar + 1 )
      break;
    end

  end
%
%  If we haven't seen a terminator, and we have examined the
%  entire string, then we're done, and LCHAR is equal to S_LENGTH.
%
  if ( iterm ~= 1 & lchar + 1 == s_length )
    lchar = s_length;
  end
%
%  Number seems to have terminated.  Have we got a legal number?
%  Not if we terminated in states 1, 2, 6 or 7!
%
  if ( ihave == 1 | ihave == 2 | ihave == 6 | ihave == 7 )
    ierror = ihave;
    fprintf ( 1, '\n' );
    fprintf ( 1, 'S_TO_R8 - Fatal error!\n' );
    fprintf ( 1, '  IHAVE = %d\n', ihave );
    error ( 'S_TO_R8 - Fatal error!' );
  end
%
%  Number seems OK.  Form it.
%
  if ( jtop == 0 )
    rexp = 1.0;
  else

    if ( jbot == 1 )
      rexp = 10.0^( jsgn * jtop );
    else
      rexp = jsgn * jtop;
      rexp = rexp / jbot;
      rexp = 10.0^rexp;
    end

  end

  r = isgn * rexp * rtop / rbot;

  return
end
function st_data_print ( nrow, ncol, nnzero, row, col, a, title )

%*****************************************************************************80
%
%% ST_DATA_PRINT prints the data in an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NROW, the assumed number of rows in the matrix.
%
%    Input, integer NCOL, the assumed number of columns in the matrix.
%
%    Input, integer NNZERO, the assumed number of nonzeros in the matrix.
%
%    Input, integer ROW(NNZERO), COL(NNZERO), the 0-based row and column
%    index of a nonzero matrix entry.
%
%    Input, real A(NNZERO), the value of a nonzero matrix entry.
%
%    Input, string TITLE, a title.
%
  st_data_print_some ( 0, nrow-1, 0, ncol-1, nnzero, row, col, a, title );

  return
end
function st_data_print_some ( row1, row2, col1, col2, nnzero, row, col, a, ...
  title )

%*****************************************************************************80
%
%% ST_DATA_PRINT_SOME prints some of the data in an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer ROW1, ROW2, the first and last rows to print.
%
%    Input, integer COL1, COL2, the first and last columns to print.
%
%    Input, integer NNZERO, the assumed number of nonzeros in the matrix.
%
%    Input, integer ROW(NNZERO), COL(NNZERO), the 0-based row and column
%    index of a nonzero matrix entry.
%
%    Input, real A(NNZERO), the value of a nonzero matrix entry.
%
%    Input, string TITLE, a title.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '%s\n', title );
  fprintf ( 1, '\n' );

  for k = 1 : nnzero
    if ( row1 <= row(k) & row(k) <= row2 & col1 <= col(k) & col(k) <= col2 )
      fprintf ( 1, '  %8d  %8d  %8d  %16.8f\n', k, row(k), col(k), a(k) );
    end
  end

  return
end
function st_data_write ( output_filename, nrow, ncol, nnzero, row, col, a )

%*****************************************************************************80
%
%% ST_DATA_WRITE writes the data of an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the ST file.
%
%    Input, integer NROW, the assumed number of rows in the matrix.
%
%    Input, integer NCOL, the assumed number of columns in the matrix.
%
%    Input, integer NNZERO, the assumed number of nonzeros in the matrix.
%
%    Input, integer ROW(NNZERO), COL(NNZERO), the 0-based row and column
%    index of a nonzero matrix entry.
%
%    Input, real A(NNZERO), the value of a nonzero matrix entry.
%
  output_unit = fopen ( output_filename, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ST_DATA_WRITE - Error!\n' );
    fprintf ( 1, '  Could not open the output file.\n' );
    error ( 'ST_DATA_WRITE - Error!' );
    return;
  end

  for k = 1 : nnzero

    fprintf ( output_unit, '  %8d  %8d  %16.8f\n', row(k), col(k), a(k) );

  end

  fclose ( output_unit );

  return
end
function st_header_print ( nrow, ncol, nnzero )

%*****************************************************************************80
%
%% ST_HEADER_PRINT prints the header of an ST file.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NROW, the assumed number of rows in the matrix.
%
%    Input, integer NCOL, the assumed number of columns in the matrix.
%
%    Input, integer NNZERO, the assumed number of nonzeros in the matrix.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  ST header information:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of rows =     %d\n', nrow );
  fprintf ( 1, '  Number of columns =  %d\n', ncol );
  fprintf ( 1, '  Number of nonzeros = %d\n', nnzero );

  return
end
function st_header_write ( output_filename, nrow, ncol, nnzero )

%*****************************************************************************80
%
%% ST_HEADER_WRITE writes the header of an ST file.
%
%  Discussion:
%
%    There is no header for an ST file.  But for completeness, we supply
%    this dummy routine.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the ST file.
%
%    Input, integer NROW, the assumed number of rows in the matrix.
%
%    Input, integer NCOL, the assumed number of columns in the matrix.
%
%    Input, integer NNZERO, the assumed number of nonzeros in the matrix.
%
  return
end
function indx = st_rebase ( base1, base2, nnzero, indx )

%*****************************************************************************80
%
%% ST_REBASE changes the base of an index array.
%
%  Discussion:
%
%    Both the ROW and COL arrays are ordinarily 0-based in the ST format.
%    FORTRAN and MATLAB expect 1-based indices.
%
%    To convert ROW and COL from 0-based to 1-based form, call this routine
%    with BASE1=0, BASE2=1.
%
%    If ROW and COL from FORTRAN or MATLAB are to be converted to ST format,
%    call this routine with BASE1=1 and BASE2=0.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    01 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer BASE1, BASE2, the old and new bases.
%
%    Input, integer NNZERO, the number of nonzeros in the matrix.
%
%    Input, integer INDX(NNZERO), the index vector to be rebased.
%
%    Output, integer INDX(NNZERO), the index vector
%    to be rebased.
%
  indx(1:nnzero) = indx(1:nnzero) - base1 + base2;

  return
end
function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end


text:
function varargout = nii_stc(nii, out, timing)
% Perform slice timing correction to the input NIfTI data.
%
% Syntax:
%  nii_stc(filename_in, filename_out);
%  nii_out = NII_STC(filename_in);
%  nii_out = NII_STC(nii_in);
%  nii_out = NII_STC(nii_in, [], FSL_slice_timing);
%
% The mandatory first input is a NIfTI file name or a struct returned by
% nii_tool('load').
%
% If the second optional input, a NIfTI file name, is provided, the corrected
% data will be saved into the file (overwrite if exists).
%
% The third optional input, slice timing, is needed only if there is no timing
% information in nii data. It should be in FSL custom timing style, one number
% per slice in unit of TR, with a positive number meaning to shift forwards in
% time. If this input is provided, it will override the parameters in nii
% header.
%
% The optional output, if requested, returns the nii struct after corrected.
%
% The dicm2nii converter stores slice timing information into nii hdr and ext.
% The numbers range from -0.5 to 0.5 in unit of TR. This means the reference
% slice is the slice acquired at the middle of TR. If nifti from dicm2nii is
% used, there is no need to provide slice order/timing information for NII_STC.
%
% See also NII_VIEWER, NII_TOOL, DICM2NII

% Xiangrui Li (xiangrui.li@gmail.com)
% 160517 Wrote it

toSave = nargin>1 && ~isempty(out);
if toSave && ~ischar(out)
    error('Second input must be nii file name to save data.');
end
if ~toSave && nargout<1
    [out, pth] = uiputfile('*.nii;*.nii.gz', 'Input file name to save result');
    out = fullfile(pth, out);
    toSave = true;
end

if ischar(nii), nii = nii_tool('load', nii); end % file name
if ~isstruct(nii) || ~all(isfield(nii, {'hdr' 'img'}))
    error('Input must be nii struct or existing nii file name.');
end

nSL = nii.hdr.dim(4);
if nargin>2 % error check for timing
    if ~isnumeric(timing) || numel(timing)~=nSL || ~isvector(timing)
        help(mfilename);
        error('timing should have one number per slice.');
    elseif min(timing)<-1 || max(timing)>1 || max(timing)-min(timing)>1
        error('timing out of range: must be in unit of TR.');
    end
end

try
    t = nii.ext.edata_decoded.SliceTiming; % FSL style timing by dicm2nii.m
catch % get it fro hdr, but no valid code for multiband sequence
    TR = nii.hdr.pixdim(5);
    dur = nii.hdr.slice_duration;
    if dur<=0 || dur>TR/nSL, dur = TR / nSL; end % best guess
    t = 0.5 - (0:nSL-1)' * dur / TR;

    switch nii.hdr.slice_code
        case 1                                % ascending
        case 2,	t(nSL:-1:1) = t;              % descending
        case 3,	t([1:2:nSL 2:2:nSL]) = t;     % interleaved ascending, start 1
        case 4, t([nSL:-2:1 nSL-1:-2:1]) = t; % interleaved descending, start n
        case 5,	t([2:2:nSL 1:2:nSL]) = t;     % interleaved ascending, start 2
        case 6,	t([nSL-1:-2:1 nSL:-2:1]) = t; % interleaved descending, start n-1
        otherwise, t = [];
    end
end

if isempty(t)
    if nargin<3
        error('There is no slice order/timing info in nii header.');
    else
        t = timing;
    end
elseif nargin>2 && any(abs(diff(t(:)-timing(:)))>0.01) % timing provided
    warning(['Provided timing is inconsistent with header timing.' char(10 )...
        'The provided timing will be used for ' nii.hdr.file_name]);
    t = timing; % warn but use the provided timing
end

nVol = single(nii.hdr.dim(5));
nFFT = 2^ceil(log2(nVol+40)); % arbituary extra padding to reduce edge problem
F = (0:nFFT/2) / nFFT; % freq in unit of points: last is Nyquist
F = [F -F(end-1:-1:2)] * 2*pi; % symmetric around Nyquist, except 1st 0
F = permute(exp(F * 1i), [1 3 4 2]); % [1 1 1 nFFT]
ramp = single(permute(linspace(1, 0, nFFT-nVol), [1 3 4 2])); % [1 1 1 nRamp]
t = single(t);

nii.img = single(nii.img); % avoid integer type
for i = 1:nSL % slice by slice fft/ifft seems faster than all together
    y = nii.img(:,:,i,:); % [r c 1 nVol]
    pad = bsxfun(@times, y(:,:,1,nVol)-y(:,:,1,1), ramp); % [r c 1 nRamp]
    pad = bsxfun(@plus, pad, y(:,:,1,1)); % linear ramp from last to 1st point
    y = cat(4, y, pad); %  [r c 1 nFFT]
    y = fft(y, [], 4);
    y = bsxfun(@times, y, F.^t(i)); % shift phase
    y = real(ifft(y, [], 4));
    nii.img(:,:,i,:) = y(:,:,1,1:nVol); % update img while drop padding
end

if toSave, nii_tool('save', nii, out); end
if nargout, varargout{1} = nii_tool('update', nii); end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function vw = viewSetTimeSeries(vw,param,val,varargin)
%Organize methods for setting view parameters.
%
% This function is wrapped by viewSet. It should not be called by anything
% else other than viewSet.
%
% This function retrieves information from the view that relates to a
% specific component of the application.
%
% We assume that input comes to us already fixed and does not need to be
% formatted again.

if notDefined('vw'),  error('No view defined.'); end
if notDefined('param'), error('No parameter defined'); end
if notDefined('val'),   val = []; end

mrGlobals;

switch param

    case 'tseries'
        vw.tSeries = val;
    case 'tseriesslice'
        vw.tSeriesSlice = val;
    case 'tseriesscan'
        vw.tSeriesScan = val;

    otherwise
        error('Unknown view parameter %s.', param);

end %switch

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% GCMEX An efficient graph-cut based energy minimization
%
% [LABELS ENERGY ENERGYAFTER] =
%   GCMEX(CLASS, UNARY, PAIRWISE, LABELCOST,EXPANSION)
%
%   Runs a minimization starting with the labels for each node defined
%   by CLASS, with unary potentials UNARY and the structure of the
%   graph and pairwise potentials defined by PAIRWISE. LABELCOST
%   determines data costs in terms of the labels of adjacent nodes.
%
% Parameters:
%   CLASS:: A 1xN vector which specifies the initial labels of each
%     of the N nodes in the graph
%   UNARY:: A CxN matrix specifying the potentials (data term) for
%     each of the C possible classes at each of the N nodes.
%   PAIRWISE:: An NxN sparse matrix specifying the graph structure and
%     cost for each link between nodes in the graph.
%   LABELCOST:: A CxC matrix specifying the fixed label cost for the
%     labels of each adjacent node in the graph.
%   EXPANSION:: A 0-1 flag which determines if the swap or expansion
%     method is used to solve the minimization. 0 == swap,
%     1 == expansion. If ommitted, defaults to swap.
%
% Outputs:
%   LABELS:: A 1xN vector of the final labels.
%   ENERGY:: The energy of the initial labeling contained in CLASS
%   ENERGYAFTER:: The energy of the final labels LABELS
%
% How do I know if I should use swap or expansion? From GC_README.txt:
%   The expansion algorithm for energy minimization can be used
%   whenever for any 3 labels a,b,c V(a,a) + V(b,c) <= V(a,c)+V(b,a).
%   In other words, expansion algorithm can be used if the binary
%   energy for the expansion algorithm step is regular, using V.
%   Kolmogorov's terminology.
%
%   The swap algorithm for energy minimization can be used whenever
%   for any 2 labels a,b V(a,a) + V(b,b) <= V(a,b)+V(b,a). In other
%   words, swap algorithm can be used if the binary energy for the
%   swap algorithm step is regular, using V. Kolmogorov's terminology.
%
% GCMex Version 2.3.0

author: Cloud-CV

repo: object-proposals


text:
function Vdata = voltageCurve(Tdata,x,Res,ThVal,ThBeta)
%% Calculate Voltage Curve given Indices (x)
% Copyright (c) 2012, MathWorks, Inc.
%% Index into Resistor and Thermistor arrays to extract component values
y(1) = Res(x(1));
y(2) = Res(x(2));
y(3) = Res(x(3));
y(4) = Res(x(4));
y(5) = ThVal(x(5));
y(6) = ThBeta(x(5));
y(7) = ThVal(x(6));
y(8) = ThBeta(x(6));

%% Calculate temperature curve for a particular set of components
Vdata = tempCompCurve(y, Tdata);

author: Sable

repo: mcbench-benchmarks


text:
%MVTBDEMO 	Machine Vision toolbox demonstrations
%
% Displays popup menu of toolbox demonstration scripts that illustrate:
%   - image processing
%   - feature extraction
%   - visual servoing
%
% Notes::
% - The scripts require the user to periodically hit <Enter> in order to move
%   through the explanation.
% - Set PAUSE OFF if you want the scripts to run completely automatically.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
function mvtbdemo(timeout)
    echo off
    close all
    echo off
    clear all
    delete( get(0, 'Children') );

    % find the path to the demos
    if exist('mvtbdemo', 'file') == 2
        tbpath = fileparts(which('mvtbdemo'));
        demopath = fullfile(tbpath, 'demos');
    end

    opts = {'begin', 'path', demopath};

    % if a timeout interval is given, add this to the options
    if nargin > 0
        opts = {opts, 'delay', timeout};
    end

    fprintf('------------------------------------------------------------\n');
    fprintf('Many of these demos print tutorial text and MATLAB commmands in the console window.\n');
    fprintf('Read the text and press <enter> to move on to the next command\n');
    fprintf('At the end of the tutorial/demo you can choose the next one from the graphical menu.\n');
    fprintf('------------------------------------------------------------\n');

    demos = {
        'Basics', 'files';
        'Image/Movie file', 'moviegrab';
        'Image/Webcam', 'webcam';
        'Image/Google Earth', 'earthgrab';
        'Image/Live video', 'livevideo';
        'Image/Track object', 'trackblue';
        'Application/Mosaicing', 'mosaic';
        'Application/ICP', 'icp';              % broken
        'Application/Line Fit', 'linefit';     % broken
        'Application/Hough transform', 'hough';
        'Application/Chroma key', 'chromakey';
        'Application/RANSAC', 'ransac';         % broken
        'Application/Background estimation', 'bgest';
        'Application/Where''s Wally', 'wally';
        'Multiview/Stereo', 'stereo';
        % blobs
        % camera
        % convdemo
        % fly
        % harris
        % homdemo
        % imagegrab
        % kernel
        % lensanim
        % opticalflow
        % poseest
        % ransacline
        % seg
        % skel
        % track
        % track2
        % vservo_points
        % xx
        'VisualServo/Mosaicing', '';
        'VisualServo/Mosaicing', 'mosaic';
        'VisualServo/Mosaicing', 'mosaic';
        'Exit', '';
        };

    while true
        selection = menu('Machine Vision Toolbox demonstrations', demos{:,1});
        if strcmp(demos{selection,1}, 'Exit')
            % quit now
            delete( get(0, 'Children') );
            break;
        else
            % run the appropriate script
            script = demos{selection,2}
            runscript(script, opts{:})
        end
    end
end

author: petercorke

repo: machinevision-toolbox-matlab


text:
%rhsta simulates the examples of Sections 3.2.3 and 3.3.4
%    ____________________________________________________________________
%  /Program:  rhsta.m													   \
% / Description:  This program runs the simulation of the Rhors example     \
%|  	found in sections 3.2.3 and 3.3.4 of the text.  The program allows	 |
%|  	the user to choose between the example of feedforward around 		 |
%|	plant or feedforward in both the plant and model.  The program 	     |
%|	prompts the user to enter values for the variables that are adjusted |
%|	within the text.													 |
%|	This program requires the simulink diagram rhsim.m in order to run   |
% \ 		/
%  \_______________________________________________________________________/

%Copyright

% Programmed Summer 1996 by Kenyon Thayer, RPI, Troy, NY

a=exist('D');
if a==1
	D_old=D;
	tau_old=tau;
	Ti_old=Ti;
	Tp_old=Tp;
else
	D_old=0.1;
	tau_old=0.2;
	Ti_old=10;
	Tp_old=10;
end

disp('  **This programs simulates the text book examples 3.2.3 and 3.3.4.**');
disp('  ');
disp('Note: example 3.2.3 is feedforward in the plant only');
disp('      example 3.3.4 is feedforward in the plant and model');

str=input('Is this feedforward in both plant and model? Y/N: ','s');
if str=='Y' | str=='y'
	ff=1;
else
	ff=0;
end
% State space form of a Model
% **************************
num1=[3];den1=[1 3];
[Am,Bm,Cm,Dm]=tf2ss(num1,den1);

% Initial Conditions
% *************************************
Key0=0;
Kx0=0;
Ku0=0;
Ki0=[Key0 Kx0 Ku0];
xm0=0;

% State space form of a FF Compensator
% *************************************
D=input(sprintf('Enter a value for D (default value is %5.3f): ',D_old));
if isempty(D)
	D=D_old;
end

tau=input(sprintf('Enter a value for tau (default value is %5.3f): ',...
	tau_old));
if isempty(tau)
	tau=tau_old;
end

Ti=input(sprintf('Enter a value for Ti (default value is %d): ',...
	Ti_old));
if isempty(Ti)
	Ti=Ti_old;
end

Tp=input(sprintf('Enter a value for Tp (default value is %d): ',...
	Tp_old));
if isempty(Tp)
	Tp=Tp_old;
end

numc=D;
denc=[tau 1];
[Af,Bf,Cf,Df]=tf2ss(numc,denc);
xf0=[0];

% State space form of a Plant
%*****************************
a=30.0;
nump=[458];denp=[1 (a+1) (a+229) 229];
[Ap,Bp,Cp,Dp]=tf2ss(nump,denp);
xp0=[0 0 0]';

%************ Ti and Tp ************
Tie=Ti;Tix=Ti;Tiu=Ti;
Tpe=Tp;Tpx=Tp;Tpu=Tp;

%*********************************************
%[t,x,y]=gear('rhsim',40);
[t,x,y]=rk45('rhsim',40,[],[1e-6,1e-2,1]);
%[t,x,y]=linsim('rhsim',10,[],[1e-5,1e-2,1-2]);

plot(t,y(:,1),'--',t,y(:,2))
axis([0 40 -.4 .4])
text(.41, 1.05, 'Ym', 'color', [1 0 1], 'FontSize', 12, 'Units', 'normal');
text(.51, 1.05, 'Yp', 'color', [1 1 0], 'FontSize', 12, 'Units', 'normal');


text:
function freq_grid = fftshiftfreqgrid(N,Fs)
%fftshiftfreqgrid Generate frequency grid

freq_res = Fs/N;
freq_grid = (0:N-1).'*freq_res;
Nyq = Fs/2;
half_res = freq_res/2;
if rem(N,2) % odd
    idx = 1:(N-1)/2;
    halfpts = (N+1)/2;
    freq_grid(halfpts) = Nyq-half_res;
    freq_grid(halfpts+1) = Nyq+half_res;
else
    idx = 1:N/2;
    hafpts = N/2+1;
    freq_grid(hafpts) = Nyq;
end

freq_grid = fftshift(freq_grid);
freq_grid(idx) = freq_grid(idx)-Fs;

end

author: Xiangyu-Gao


text:
function res = afni_isdigit(s)
% res = afni_isdigit(s);
% res(i) = 1 is s(1) is a digit between 0 and 9 (inclusive)
%
digit_msk=s >= '0' & s <= '9';
res = zeros(size(s));
res(digit_msk) = 1;

return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function mmse(filename,outfile,SPU)

%
%  Implements the MMSE algorithm [1].
%
%  Usage:  mmse(noisyFile, outputFile, SPU)
%
%         infile - noisy speech file in .wav format
%         outputFile - enhanced output file in .wav format
%         SPU  - if 1, includes speech-presence uncertainty
%                if 0, doesnt include speech-presence uncertainty
%
%
%  Example call:  mmse('sp04_babble_sn10.wav','out_mmse.wav',1);
%
%  References:
%   [1] Ephraim, Y. and Malah, D. (1985). Speech enhancement using a minimum
%       mean-square error log-spectral amplitude estimator. IEEE Trans. Acoust.,
%       Speech, Signal Process., ASSP-23(2), 443-445.
%
% Authors: Philipos C. Loizou
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $  $Date: 10/09/2006 $
%-------------------------------------------------------------------------

if nargin<3
    fprintf('Usage: mmse(infile.wav,outfile.wav,SPU) \n');
    fprintf('where SPU=1 - includes speech presence uncertainty\n');
    fprintf('      SPU=0 - does not includes speech presence uncertainty\n\n');
    return;
end;

if SPU~=1 & SPU~=0
    error('ERROR: SPU needs to be either 1 or 0.');
end

[x, Srate, bits]= wavread( filename);

% =============== Initialize variables ===============

len=floor(20*Srate/1000); % Frame size in samples
if rem(len,2)==1, len=len+1; end;
PERC=50; % window overlap in percent of frame size
len1=floor(len*PERC/100);
len2=len-len1;

win=hanning(len);  % define window
win = win*len2/sum(win);  % normalize window for equal level output

% Noise magnitude calculations - assuming that the first 6 frames is noise/silence
%
nFFT=2*len;
j=1;
noise_mean=zeros(nFFT,1);
for k=1:6
    noise_mean=noise_mean+abs(fft(win.*x(j:j+len-1),nFFT));
    j=j+len;
end
noise_mu=noise_mean/6;
noise_mu2=noise_mu.^2;

%--- allocate memory and initialize various variables

k=1;
img=sqrt(-1);
x_old=zeros(len1,1);
Nframes=floor(length(x)/len2)-1;
xfinal=zeros(Nframes*len2,1);

% --------------- Initialize parameters ------------
%
k=1;
aa=0.98;
eta= 0.15;
mu=0.98;
c=sqrt(pi)/2;
qk=0.3;
qkr=(1-qk)/qk;
ksi_min=10^(-25/10); % note that in Chap. 7, ref. [17], ksi_min (dB)=-15 dB is recommended

%===============================  Start Processing =======================================================
%
for n=1:Nframes

    insign=win.*x(k:k+len-1);

    %--- Take fourier transform of  frame
    %
    spec=fft(insign,nFFT);
    sig=abs(spec); % compute the magnitude
    sig2=sig.^2;

    gammak=min(sig2./noise_mu2,40);  % posteriori SNR
    if n==1
        ksi=aa+(1-aa)*max(gammak-1,0);
    else
        ksi=aa*Xk_prev./noise_mu2 + (1-aa)*max(gammak-1,0);
        % decision-direct estimate of a priori SNR
        ksi=max(ksi_min,ksi);  % limit ksi to -25 dB
    end

    log_sigma_k= gammak.* ksi./ (1+ ksi)- log(1+ ksi);
    vad_decision= sum( log_sigma_k)/nFFT;
    if (vad_decision< eta) % noise only frame found
        noise_mu2= mu* noise_mu2+ (1- mu)* sig2;
    end
    % ===end of vad===

    vk=ksi.*gammak./(1+ksi);
    [j0,err]=besseli(0,vk/2);
    [j1,err2]=besseli(1,vk/2);
    if any(err) | any(err2)
        fprintf('ERROR! Overflow in Bessel calculation in frame: %d \n',n);
    else
        C=exp(-0.5*vk);
        A=((c*(vk.^0.5)).*C)./gammak;
        B=(1+vk).*j0+vk.*j1;
        hw=A.*B;
    end

    % --- estimate speech presence probability
    %
    if SPU==1
        evk=exp(vk);
        Lambda=qkr*evk./(1+ksi);
        pSAP=Lambda./(1+Lambda);
       sig=sig.*hw.*pSAP;
    else
        sig=sig.*hw;
    end

    Xk_prev=sig.^2;  % save for estimation of a priori SNR in next frame

    xi_w= ifft( sig .* exp(img*angle(spec)),nFFT);

    xi_w= real( xi_w);

    xfinal(k:k+ len2-1)= x_old+ xi_w(1:len1);
    x_old= xi_w(len1+ 1: len);

    k=k+len2;

end
%========================================================================================

wavwrite(xfinal,Srate,16,outfile);

author: jtkim-kaist

repo: Speech-enhancement


text:
function T2 = sumv(T1, sum_over)
%
% Like the built in sum, but will sum over several dimensions and then squeeze the result.

T2 = T1;
for i=1:length(sum_over)
  if sum_over(i) <= ndims(T2) % prevent summing over non-existent dimensions
    T2=sum(T2, sum_over(i));
  end
end
T2 = squeeze(T2);


text:
function M = imfiltermtx(F, sz, bndry);
%IMFILTERMTX   N-D image filtering matrix
%   M = IMFILTERMTX(F, SZ[, BNDRY]) returns the image filtering matrix
%   for the matrix F.  SZ gives the size of the array that the filtering
%   (convolution) should be applied to.  The returned matrix M is sparse.
%   If X is of size SZ, then reshape(M * X(:), SZ + size(F) - 1) is the
%   same as imfilter(X, F, BNDRY, 'conv').
%
%   The optional parameter BNDRY controls the boundary handling:
%     * 'replicate':  Input array values outside the bounds of the array
%                     are assumed to equal the nearest array border
%                     value.
%     * 'circular':   Input array values outside the bounds of the array
%                     wrap around to the opposite bound.
%     * '0':          (default) Input array values outside of the bounds
%                     of the array are assumed to be zero.
%
%   See also IMFILTER, MAKE_IMFILTER_MAT.
%
%   Author:  Stefan Roth, Department of Computer Science, TU Darmstadt
%   Contact: sroth@cs.tu-darmstadt.de
%   $Date: 2007-03-27 14:09:11 -0400 (Tue, 27 Mar 2007) $
%   $Revision: 252 $

% Copyright 2004-2007, Brown University, Providence, RI. USA
% Copyright 2007-2010 TU Darmstadt, Darmstadt, Germany.
%
%                          All Rights Reserved
%
% All commercial use of this software, whether direct or indirect, is
% strictly prohibited including, without limitation, incorporation into in
% a commercial product, use in a commercial service, or production of other
% artifacts for commercial purposes.
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for research purposes is hereby granted without fee,
% provided that the above copyright notice appears in all copies and that
% both that copyright notice and this permission notice appear in
% supporting documentation, and that the name of the author and Brown
% University not be used in advertising or publicity pertaining to
% distribution of the software without specific, written prior permission.
%
% For commercial uses contact the Technology Venture Office of Brown University
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO
% THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
% FITNESS FOR ANY PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR
% BROWN UNIVERSITY BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
% DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
% PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
% ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
% THIS SOFTWARE.

  % Size of padding for boundary handling
  pad_size = 0.5 * (size(F) - 1);

  % Corners of the original image in the padded image
  corner_lo = 1 + pad_size;
  corner_hi = size(F) + sz - 1 - pad_size;

  % Pixel index of the output image M
  M_idx = zeros(sz + size(F) - 1);
  M_idx(:) = 1:numel(M_idx);

  % Pixel index of the "input" image X
  X_idx = zeros(sz);
  X_idx(:) = 1:numel(X_idx);
  % Pad it to the size of M
  X_idx = [zeros(pad_size(1), sz(2) + size(F, 2) - 1); ...
           zeros(sz(1), pad_size(2)), X_idx, zeros(sz(1), pad_size(2)); ...
           zeros(pad_size(1), sz(2) + size(F, 2) - 1)];

  switch(bndry)
    case 'symmetric'
      error('Not implemented yet!')

    case 'replicate'
      % Replicate pixel values from the edges & corners

      % Fill all 4 corners
      for c = 0:3

        % Position or corner and range of block next to corner
        pos = cell(1, 2);
        rng = cell(1, 2);
        for d = 1:2
          if bitand(c, 2^(d-1))
            pos{d} = corner_hi(d);
            rng{d} = corner_hi(d)+1:size(M_idx, d);
          else
            pos{d} = corner_lo(d);
            rng{d} = 1:corner_lo(d)-1;
          end
        end

        % Get index of corner pixel
        idx = sub2ind(size(M_idx), pos{:});

        % Set entire block to corner pixel index
        M_idx(rng{:}) = idx;

      end

      % Fill all 4 edges
      for d = 1:2
        for c = 0:3

          % I don't quite understand this code anymore, but it appears to
          % work fine :)
          pos = cell(1, 2);
          rng = cell(1, 2);
          for e = 1:2
            if (d == e)
              pos{d} = corner_lo(e):corner_hi(e);
              rng{d} = corner_lo(e):corner_hi(e);
            else
              if bitand(c, 2^(e-1))
                pos{e} = repmat(corner_hi(e), 1, sz(d));
                rng{e} = corner_hi(e)+1:size(M_idx, e);
              else
                pos{e} = repmat(corner_lo(e), 1, sz(d));
                rng{e} = 1:corner_lo(e)-1;
              end
            end
          end

          % Get indices of edge pixels
          idx = sub2ind(size(M_idx), pos{:});

          % Generate fill indices from the nearest edge pixel
          tmp_sz = pad_size;
          tmp_sz(d) = 1;
          tmp2_sz = ones(1, 2);
          tmp2_sz(d) = sz(d);
          M_idx(rng{:}) = repmat(reshape(idx, tmp2_sz), tmp_sz);
        end
      end

      % Convert indices from the output image to the input image
      M_idx = X_idx(M_idx);

    case 'circular'
      % Create padded indices of rows and columns
      r = -pad_size(1)+1:sz(1)+pad_size(1);
      c = -pad_size(2)+1:sz(2)+pad_size(2);

      % Take the modolus to make sure they are in the valid bounds
      r = 1 + mod(r - 1, sz(1));
      c = 1 + mod(c - 1, sz(2));

      M_idx = sub2ind(sz, repmat(r(:), 1, length(c)), ...
                      repmat(c(:)', length(r), 1));

    otherwise
      % Assume all zeros outside of the boundaries
      M = convmtxn(F, sz);
      return;

  end

  % The pixel indices give a mapping of the columns of the padded image
  % to the output image
  col_map = M_idx(:);

  % Get convolution matrix and split the resulting sparse array
  M = make_convn_mat(F, sz + size(F) - 1, 'same');
  [rows, cols, vals] = find(M);

  % Re-assemble the sparse array, but use the column map to accumulate
  % rows appropriately
  M = sparse(rows, col_map(cols), vals, prod(size(F) + sz - 1), prod(sz));

author: kristinbranson

repo: JAABA


text:
function spectralWidth = getSpectralWidth(MRSIStruct)
    spectralWidth = MRSIStruct.spectralWidth;
end

author: CIC-methods

repo: FID-A

save_path: github-repos/MATLAB/CIC-methods-FID-A


text:
function pf = mtimes(arg1,arg2)
% scaling of PoleFigures, implements pf1 * b and a * pf2
%
% overload the * operator, i.e. one can now write x * pf in order to
% scale the @PoleFigure pf by the factor x
%
% See also
% PoleFigure/PoleFigure PoleFigure/plus PoleFigure/minus

pf = arg1 .* arg2;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
classdef dbn
% Deep Belief Network object.
%-----------------------------------------------------------------------------
% Initialize, train, and fine-tune a deep belief network.
% Note this function calls the rbm.m object to form layers.
%-----------------------------------------------------------------------------
% Dustin E. Stansbury
% stan_s_bury@berkeley.edu

properties
	class = 'dbn';			% OBJECT CLASS
	arch;					% ARCHITECTURE/OPTIONS
	inputType='binary';		% DEFAULT INPUT TYPE
	classifier=false;		% CLASSIFIER LAYER AT TOP
	nLayers;				% TOTAL # OF LAYERS (INCLUDING INPUT)
	nRBMLayers;				% # OF RBM LAYERS
	rbmLayers;				% RBM LAYERS
	verbose = 1;			% DISPLAY OUTPUT
	saveDir
end

methods
	function self = dbn(arch)
	% d = dbn(arch)
	%--------------------------------------------------------------------------
	% DBN constructor
	%--------------------------------------------------------------------------
	% INPUT:
	%  <arch>:  - a set of arguments defining the DBN architecture.
	%
	% OUTPUT:
	%     <d>:  - an DBN model object.
	%--------------------------------------------------------------------------
		self = self.init(arch);
	end

	function print(self)
	% print()
	%--------------------------------------------------------------------------
	% Display the properties and methods for the DBN object class.
	%--------------------------------------------------------------------------
		properties(self)
		methods(self)
	end

	function self = train(self,X,targets);
	% self = train(X,[targets])
	%--------------------------------------------------------------------------
	% Perform greedy layer-wise training of of a DBN
	%--------------------------------------------------------------------------
	% INPUT:
	%        <X>:  - to-be modeled data |X| = [#Obs x #Vis]
	%  <targets>:  - category targets. can be either [#Obs x #Class] as 1 of
	%                K representation or [#Obs x 1], where each entry is a
	%                numerical category label. (optional)
	% OUTPUT:
	%     <self>:  - trained DBN object.
	%--------------------------------------------------------------------------
		for iL = 1:self.nRBMLayers
			if self.verbose,
				self.printProgress('layerTrain',iL);
			end

			% TRAIN AN RBM LAYER
			if iL == self.nRBMLayers & self.classifier
				self.rbmLayers{iL} = self.rbmLayers{iL}.train(X,targets);
			else

				self.rbmLayers{iL} = self.rbmLayers{iL}.train(X,[]);
			end

			% PREPROCESS INPUT FOR NEXT LAYER
			X = self.rbmLayers{iL}.hidExpect(X);
		end
		if ~isempty(self.saveDir),
			d = self;
			save(fullfile(self.saveDir,'dbn.mat'),'d','-v7.3');
		end
	end

	function net = fineTune(self,X,targets,ftArch);
	%  net = fineTune(X,targets,ftArch);
	%--------------------------------------------------------------------------
	% Initialize a multi-layer neural network with trained DBN features and
	% fine tune for a task using backprop (stochastic gradient descent)
	%--------------------------------------------------------------------------
	% INPUT:
	%        <X>:  - to-be modeled data |X| = [#Obs x #Vis]
	%  <targets>:  - target variables.
	%   <ftArch>:  - architecture for the corresponding neural network (see
	%                mlnn.m)
	%
	% OUTPUT:
	%      <net>:  - a trained multi-layer neural network
	%--------------------------------------------------------------------------
		[nObs,nOut] = size(targets);
		if self.verbose
			self.printProgress('fineTune')
		end

		ftArch.size = [self.arch.size,nOut];

		net = mlnn(ftArch);
		% INITIALIZE MLNN WITH DBN WEIGHTS
		for iL = 1:numel(self.arch.size)-1
			net.layers{iL}.W = self.rbmLayers{iL}.W';
			net.layers{iL}.b = self.rbmLayers{iL}.c';
		end
		% FINE TUNE
		net = net.train(X,targets);
	end

	function [outAct,outProb,layerOut] = fProp(self,X,targets);
	% [outAct,outProb] = fProp(X,[targets]);
	%--------------------------------------------------------------------------
	% Forward propagate signal up through the layers of the dbn.
	%--------------------------------------------------------------------------
	% OUTPUT:
	%  <outAct>:  - activation of the top layer/output hidden units
	% <outProb>:  - probabilities of top layer hidden units
	%--------------------------------------------------------------------------
		if notDefined('targets'),
			targets = 0;
		end

		for iL = 1:self.nRBMLayers-1
			X = self.rbmLayers{iL}.hidExpect(X);
		end
		layerOut = self.rbmLayers{end}.hidGivVis(X,targets,1);
		outAct = layerOut.aHid;
		outProb = layerOut.pHid;
	end

	function [pred,error,misClass] = classify(self,X,targets)
	%  [pred,error,misClass] = classify(X,[targets])
	%--------------------------------------------------------------------------
	% Forward propagate data <X> through the layers of the DBN and classify.
	% If <targets> provided, also returns the classification error <error>,
	% and the indices of the misclassified inputs <misClass>.
	%--------------------------------------------------------------------------

		if notDefined('targets')
			targets = [];
		else
			if isvector(targets);
				targets = self.oneOfK(targets);
			end
		end

		% PROPAGATE DATA SIGNALS UP
		for iL = 1:self.nRBMLayers-1
			X = self.rbmLayers{iL}.hidExpect(X);
		end

		% CALCULATE TOP LAYER ACTIVATIONS
		nClasses = self.rbmLayers{end}.nClasses;
		tmp = self.rbmLayers{end}.hidGivVis(X,zeros(size(X,1),nClasses),1);

		% CALCULATE CLASS PROBABILITY
		pClass = tmp.softMax(bsxfun(@plus,tmp.aHid*tmp.classW',tmp.d));

		% WINNER-TAKE-ALL CLASSIFICATION
		[~, pred] = max(pClass,[],2);

		if ~notDefined('targets') && nargout > 1
			[~,targets] = max(targets,[],2);

			% CALCULATE MISSCLASSIFICATION RATE
			misClass = find(pred~=targets);
			error = numel(misClass)/size(X,1);
		end
	end

	function self = init(self,arch)
	% dbn = init(arch)
	%--------------------------------------------------------------------------
	% Initialize deep belief network, based on architecture structure <arch>.
	%--------------------------------------------------------------------------
	% INPUT:
	% <arch>:  - Is an archtiecture struct with the fields:
	%            .shape    --> [#Vis #Hid1, ..., #HidN, #out]
	%            .lRate    --> the learning rate for each layer
	%            .sparsity --> the sparsity target for hidden units
	%            .dropout  --> the hidden unit dropout rate for each rbm layer
	%--------------------------------------------------------------------------
		arch = self.ensureArchetecture(arch);

		% GLOBAL OPTIONS
		if isfield(arch,'opts');
			opts = arch.opts;
			fn = fieldnames(self);
			for iA = 1:2:numel(opts)
				if ~isstr(opts{iA})
					error('<opts> must be a cell array string-value pairs.');
				elseif sum(strcmp(fn,opts{iA}))
					self.(opts{iA})=opts{iA+1};
				end
			end
		end

		self.arch = arch;
		self.nLayers = numel(arch.size);
		self.nRBMLayers = self.nLayers - 1;

		% INITIALIZE RBM LAYERS
		for iL = 2:self.nLayers
			layerArch.size = [arch.size(iL-1),arch.size(iL)];
			if iL == 2
				layerArch.inputType = self.inputType;
			else
				layerArch.inputType = 'binary';
			end

			self.rbmLayers{iL-1} = rbm(layerArch);
			self.rbmLayers{iL-1}.lRate = arch.lRate(iL-1);
			self.rbmLayers{iL-1}.nEpoch = arch.nEpoch(iL-1);
			self.rbmLayers{iL-1}.sparsity = arch.sparsity(iL-1);
			self.rbmLayers{iL-1}.dropout = arch.dropout(iL-1);

			if self.classifier & (iL == self.nLayers)
				self.rbmLayers{end}.classifier = true;
			end
		end
	end

	function arch = ensureArchetecture(self,arch)
	% arch = ensureArchitecture(arch)
	%--------------------------------------------------------------------------
	% Preprocess the provided architecture structure <arch>.
	%--------------------------------------------------------------------------
	% INPUT:
	% <arch>:  - is either a [2 x 1] vector giving the [#Vis x #Hid], in which
	%            case we use the default model parameters, or it is a strucure
	%            with any of the fields:
	%             .size                --> Network size; [#Vis x # Hid]
	%             .nEpoch(optional)    --> # of training epochs for each layer
	%             .lRate (optional)    --> learing rate for each layer
	%             .sparsity (optional) --> sparsity rate for each layer
	%             .dropout (optional)  --> dropout rate for each layer
	%--------------------------------------------------------------------------
		if ~isstruct(arch),arch.size = arch; end

		nLayers = numel(arch.size);

		% CHECK # EPOCHS
		if ~isfield(arch,'nEpoch')
			arch.lRate = repmat(100,1,nLayers-1);
		elseif numel(arch.nEpoch) == 1
			arch.nEpoch = repmat(arch.nEpoch,1,nLayers-1);
		end

		% CHECK LEARNING RATES
		if ~isfield(arch,'lRate')
			arch.lRate = repmat(.01,1,nLayers-1);
		elseif numel(arch.lRate) == 1
			arch.lRate = repmat(arch.lRate,1,nLayers-1);
		end

		% CHECK SPARSITY TARGET RATES
		if ~isfield(arch,'sparsity')
			arch.sparsity = repmat(0,1,nLayers-1);
		elseif numel(arch.sparsity) == 1
			arch.sparsity = repmat(arch.sparsity,1,nLayers-1);
		end

		% CHECK DROPOUT
		if ~isfield(arch,'dropout')
			arch.dropout = repmat(0,1,nLayers-1);
		elseif numel(arch.dropout) == 1
			arch.dropout = repmat(arch.dropout,1,nLayers-1);
		end
	end

	function printProgress(self,type,aux)
	% VERBOSE
		switch type
			case 'layerTrain'
				fprintf('\n\nTraining RBM Layer: %i / %i\n',aux(1),self.nLayers-1);
			case 'fineTune'
				fprintf('\n\nFine-tuning using backprop:\n')
			case 'save'
				fprintf('\nSaving...\n\n');
		end
	end
end

end

author: dustinstansbury

repo: medal


text:
function [F10ycp, F01ycp] = Q1001ycpowp(F10, F01, p, c)
%------------------------------------------------------------------------------
%
% This function computes the gridfunction {F10ycp U F01ycp} of size {F10 U F01}
% with values (x-c)^p at its gridpoints. Piecewise constant approximation is
% assumed.
%
%   Orientation
%
%         x
%     o---->
%     |
%   y |
%     v
%
%
% Design and implementation by:
% Dr. Paul M. de Zeeuw <Paul.de.Zeeuw@cwi.nl>  http://homepages.cwi.nl/~pauldz/
% Last Revision: December 12, 2000.
% (c) 1999-2000 Stichting CWI, Amsterdam
%------------------------------------------------------------------------------
if isempty(F10) || isempty(F01)
  error(' Q1001ycpowp - at least one colour not present (empty) ')
else
  [n10, m10] = size(F10);
  [n01, m01] = size(F01);
  if p == 0
    F10ycp = ones(n10, m10);
    F01ycp = ones(n01, m01);
  else
    [hx, hy] = Q1001gridfdims(F10, F01);
%
    yfirst =  0 + hy/2 - c;
    ylast  =  yfirst + (m01 - 1) * 2 * hy;
    lineofycp = linspace(yfirst, ylast, n01)';
    F01ycp = lineofycp(:,ones(m01,1));
    if p ~= 1
      F01ycp = F01ycp.^p;
    end
%
    yfirst =  0 + 3*hy/2 - c;
    ylast  =  yfirst + (m10 - 1) * 2 * hy;
    lineofycp = linspace(yfirst, ylast, n10)';
    F10ycp = lineofycp(:,ones(m10,1));
    if p ~= 1
      F10ycp = F10ycp.^p;
    end
  end
end
%------------------------------------------------------------------------------


text:
% Packs an object of any complexity (structure, cell array ect.)
% into a serialized form.
%
% mex-function, written by Dima Mihaylov (dmih@in-solve.ru)
%
% See also DESERIALIZE


text:
%% CUBEAFEMQUADCURL quad curl equations on the unit cube
%
%   CUBEAFEMQUADCURL computes ND0-(CR-P0)-ND0 nonconforming approximations
%   of the quad curl equations in the unit cube. The mesh is refined
%   adaptively guided by a residual-based estimator.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear;

%% Set up
maxN = 5e5;
maxIt = 20;
theta = 0.3;
N = zeros(maxIt,1);
h = zeros(maxIt,1);
erru = zeros(maxIt,1);
errw = zeros(maxIt,1);
errphi = zeros(maxIt,1);
etaTotal = zeros(maxIt,1);
etaw = zeros(maxIt,1);
etaphi = zeros(maxIt,1);
etau = zeros(maxIt,1);

%% Generate initial mesh
H = pi;
[node,elem,HB] = cubemesh([0,H,0,H,0,H],H/2);
bdFlag = setboundary3(node,elem,'Dirichlet');

%% PDE and options
pde = quadCurlDataSmooth1;

option.solver = 'direct';
option.printlevel = 0;

%% Finite Element Method
for k = 1:maxIt

    %% SOLVE
    [soln,eqn] = quadcurl3NC(node,elem,bdFlag,pde,option);
    N(k) = 2*length(soln.u)+3*length(soln.phi);
    %% ESTIMATE
    [erru(k), errK] = getHcurlerror3ND(node,elem,pde.curlu,soln.u);
    [etaK, est] = estimatequadcurl(node,elem,soln,pde,option);
    etaTotal(k) = sqrt(sum(etaK.^2));
    etaw(k) = sqrt(sum((est.eta1).^2));
    etaphi(k) = sqrt(sum((est.eta2).^2));
    etau(k) = sqrt(sum((est.eta3).^2));
    fprintf("Iter %2d, # Dof= %6d; \n||curl(u-u_h)|| = %8.6g,eta = %8.6g \n \n",...
           k, N(k), erru(k), etaTotal(k))
    %% Visualize
    figure(1);
    set(gcf, 'Position', [100 600 1000 300])
    subplot(1,3,1)
    showboundary3(node,elem,'z<pi/2','facealpha',0.5);

    subplot(1,3,2)
    errorE2V = accumarray(elem(:),repmat(errK,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,errorE2V,'z<pi/2 ','EdgeColor','k');

    subplot(1,3,3)
    etaE2V = accumarray(elem(:),repmat(etaK,[4,1]),[size(node,1), 1]);
    showsolution3(node,elem,etaE2V,'z<pi/2 ','EdgeColor','k');

    drawnow;
    %% MARK and REFINE
    if N(k) > maxN
        break;
    end

    markedElem = mark(elem,etaK,theta);
    if k < maxIt
        [node,elem,bdFlag,HB] = bisect3(node,elem,markedElem,bdFlag,HB);
        [elem,bdFlag] = sortelem3(elem,bdFlag);
    end

end

%%
close all;
figure(1);
r1 = showrate(N,erru,10,'r-o');
r2 = showrate(N,etaTotal,10,'b-*');
r3 = showrate(N,etaw,10,'color',[0.2, 0.8,0.7],'marker','*');
r4 = showrate(N,etaphi,10,'color',[0.2, 0.8,0.7],'marker','*');
r5 = showrate(N,etau,10,'color',[0.2, 0.8,0.7],'marker','*');
set(gca,'xscale','log','yscale','log');
T = title('Convergence');
XL = xlabel('$\#$ DoFs');
YL = ylabel('Error Magnitudes');
L = legend('$\Vert  \nabla\times(u-u_h)\Vert$', ...
    ['$N^{' num2str(r1) '}$'],...
    '$\eta(w_h,\phi_h,u_h)$', ...
    ['$N^{' num2str(r2) '}$'],...
    '$\eta_1(w_h)$', ['$N^{' num2str(r3) '}$'],...
    '$\eta_2(\phi_h)$', ['$N^{' num2str(r4) '}$'],...
    '$\eta_3(u_h)$', ['$N^{' num2str(r5) '}$'],...
    'LOCATION','Best');
set([XL,YL,L],'Interpreter','latex','FontSize', 12);
set(T,'Interpreter','latex','FontSize',16);
grid on;
set(gcf,'color','w','Position', [100 200 350 500])
drawnow;

%% sanity check
if N<5e3
    figure(2);
    subplot(1,2,1);
    option.scale = 2;
    option.plot = 'quiver3';
    center = (node(elem(:,1),:) + node(elem(:,2),:) + ...
        node(elem(:,3),:) + node(elem(:,4),:))/4;
    showvector3(center,pde.curlu(center),option);
    subplot(1,2,2);
    curluh = curlu3(node,elem,soln.u);
    showvector3(center,curluh,option);
    % showvector3(center,curlwh_comp2,option);

    set(gcf,'color','w','Position', [500 600 750 300])
    drawnow;
end

repo: ifem