text:
function hb_to_mmm ( input_file, output_file )

%*****************************************************************************80
%
%% HB_TO_MM transfers sparse matrix information from an HB file to an MM file.
%
%  Discussion:
%
%    An HB file is a sparse matrix stored in the Harwell Boeing format.
%
%    An MM file is a sparse matrix stored in the Matrix Market format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE, the name of the HB file to be read.
%
%    Input, string OUTPUT_FILE, the name of the MM file to be created.
%
  a = hb_to_msm ( input_file );

  msm_to_mm_coordinate_real_general ( output_file, a );

  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 a = hb_to_msm ( input_file )

%*****************************************************************************80
%
%% HB_TO_MSM reads a Harwell-Boeing Sparse Matrix File, and returns it in MATLAB Sparse format.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE, the name of the file containing the
%    Harwell-Boeing matrix data.
%
%    Output, real A, a sparse matrix in MATLAB sparse matrix format.
%
  input_unit = fopen ( input_file );

  if ( input_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_TO_MSM - 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

  a = sparse ( rowind, colind, values, nrow, ncol );

  fclose ( input_unit );

  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 msm_to_mm_coordinate_real_general ( output_filename, a )

%*****************************************************************************80
%
%% MSM_TO_MM_COORDINATE_REAL_GENERAL writes a "matrix coordinate real general" Matrix Market file.
%
%  Modified:
%
%    02 November 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string OUTPUT_FILENAME, the name of the file to which the information
%    should be written.
%
%    Input, sparse matrix A, the NROW by NCOL matrix, stored in MATLAB sparse
%    matrix format, which is to be written to the file.
%
  [ nrow, ncol ] = size ( a );
  nnzeros = nnz ( a );

  fid = fopen ( output_filename, 'wt+' );

  if ( fid < 0 );
    fprintf ( 1, '\n' );
    fprintf ( 1, 'MSM_TO_MM_COORDINATE_REAL_GENERAL - Fatal error!\n' );
    fprintf ( 1, '  Cannot open the output file.\n' );
    error ( 'MSM_TO_MM_COORDINATE_REAL_GENERAL - Fatal error!');
  end;

  fprintf ( fid, '%%%%MatrixMarket matrix coordinate real general\n');
  fprintf ( fid, '%%%%  Created by MSM_TO_MM_COORDINATE_REAL_GENERAL.M\n' );
  fprintf ( fid, '  %d  %d  %d\n', nrow, ncol, nnzeros );

  for j = 1 : ncol
    [ rows, temp, vals ] = find ( a(:,j) );
    sz = size ( rows );
    for k2 = 1 : sz
      fprintf ( fid, '  %d  %d  %f\n', rows(k2), j, vals(k2) );
    end

  end

  fclose ( fid );

  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 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

author: johannesgerer

repo: jburkardt-m


text:
function f = cumsum2(f, dims)
%CUMSUM2   Double indefinite integral of a CHEBFUN3.
%   F = CUMSUM2(F) returns the double indefinite integral of a CHEBFUN3 F.
%   By default, that means cumsum in the first 2 variables, i.e., x and y:
%                   y  x
%                  /  /
%   CUMSUM2(F) =  |  |   F(x,y,z) dx dy
%                 /  /
%                c  a
%
%   where [a,b] x [c,d] x [e,g] is the domain of F.
%
%   DIMS is a vector containing two of the three indices 1,2,3 to show
%   which two of the dimensions are to be used.
%
% See also CHEBFUN3/CUMSUM, CHEBFUN3/CUMSUM3, CHEBFUN3/SUM, CHEBFUN3/SUM2
% and CHEBFUN3/SUM3.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% Check for empty:
if ( isempty(f) )
    f = [];
    return
end

% Default to cumsum in the first two variables x and y:
if ( nargin == 1 )
    dims = [1, 2];
end

if ( numel(dims) ~= 2 )
    error('CHEBFUN:CHEBFUN3:cumsum2:dims', 'Dims should have 2 entries.');
end

if ismember(1, dims)
    % cumsum along the 1st variable
    f.cols = cumsum(f.cols);
end
if ismember(2, dims)
    % cumsum along the 2nd variable
    f.rows = cumsum(f.rows);
end
if ismember(3, dims)
    % cumsum along the 3rd variable
    f.tubes = cumsum(f.tubes);
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
% Creates a small neural network to check the backpropagation gradients.
% It will output the analytical gradients
% produced by your backprop code and the numerical gradients (computed
% using computeNumericalGradient). These two gradient computations should
% result in very similar values.
function debug_nn_gradients(lambda)
    input_layer_size = 3;
    hidden_layer_size = 5;
    num_labels = 3;
    m = 5;

    layers = [
        input_layer_size,
        hidden_layer_size,
        num_labels
    ];

    % We generate some 'random' test data.
    Theta1 = debug_initialize_weights(hidden_layer_size, input_layer_size);
    Theta2 = debug_initialize_weights(num_labels, hidden_layer_size);

    % Reusing debug_initialize_weights to generate X.
    X  = debug_initialize_weights(m, input_layer_size - 1);
    y  = 1 + mod(1:m, num_labels)';

    % Unroll parameters
    nn_params_unrolled = [Theta1(:); Theta2(:)];

    % Short hand for cost function
    gradient_step = @(p) nn_gradient_step(p, layers, X, y, lambda);

    [cost, grad] = gradient_step(nn_params_unrolled);
    numgrad = debug_numerical_gradient(gradient_step, nn_params_unrolled);

    % Visually examine the two gradient computations.  The two columns
    % you get should be very similar.
    disp([numgrad grad]);
    fprintf(['The above two columns you get should be very similar.\n' ...
            '(Left - Your Numerical Gradient, Right - Analytical Gradient)\n\n']);

    % Evaluate the norm of the difference between two solutions.
    % If you have a correct implementation, and assuming you used EPSILON = 0.0001
    % in debug_numerical_gradient.m, then diff below should be less than 1e-9.
    diff = norm(numgrad - grad) / norm(numgrad + grad);

    fprintf(['If your backpropagation implementation is correct, then \n' ...
            'the relative difference will be small (less than 1e-9). \n' ...
            '\nRelative Difference: %g\n'], diff);
end

author: trekhleb

repo: machine-learning-octave


text:
% Copyright Mohammad SAFEEA, 17th-Aug-2017

% Simulation of controlling KUKA iiwa robot using game pad.
% In this script, the code is testsed, a graph showing the joints angles is
% showed to confirm the correctness of the code

clear all;
close all;
clc;
% Start the joystick, make sure that the joystick is connected
instrreset;
ID=1;
joy=vrjoystick(ID);
% Initial configuration
jPos={pi / 180 * 30, pi / 180 * 30, 0, -pi / 180 * 60, 0,...
                        pi / 180 * 90, 0};
% add the KST directory
cDir = pwd;
cDir=getTheKSTDirectory(cDir);
addpath(cDir);
% IK solver parameters
numberOfIterations=10;
lambda=0.1;
TefTool=eye(4);

qin=zeros(7,1);
for i=1:7
    qin(i)=jPos{i};
end

[Tt,j]=directKinematics(qin,TefTool); % EEF frame transformation amtrix

% Draw the interface
%[figureHandle,anglesTextHandlesCellArray,labelTextHandlesCellArray]=constructInterface();
% max angular velocity
w=5*pi/180; % rad/sec
% max linear velocity
v=0.05; % m/sec
% Joint space control
firstExecution=0;

% setBackGroundColor(labelTextHandlesCellArray, 1 );
figure();
bufferSize=100;

tempVec=[1:bufferSize];
qVec=zeros(7,bufferSize);

pause(0.1);
% frameUpdateCount=0;
% maxDisp=1*pi/180;
while true
    joyStatus=read(joy);
    % Remove the bias in the analog signal
    analogPrecission=1/20;
    for i=1:4
        if abs( joyStatus(i))<analogPrecission
            joyStatus(i)=0;
        end
    end

    % About the variable (joyStatus)
    % joyStatus: is a 4x1 vector, the first and the second elements of this
    % vector correspond to the analog values of the left analog-stick
    % psoition. the third and the fourth elements of this
    % vector correspond to the analog values of the right analog-stick
    % psoition.
    %
    % Only the seconds and the third elements of this vector are used for
    % controlling the joints motion of the robot one at a time.
    %
    % The 2d element is used to select the robot axes, this element corresponds
    % to the (up,down) directions of the left analog stick of the gamepad.
    %
    % The 3rd element is used to control the motion of the selcted axes,
    % this element corresponds to the (right, left) directions of the right
    % analog stick of the gamepad.

    % Calculate the elapsed time
    if(firstExecution==0)
        firstExecution=1;
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=0; % elapsed is zero at first excution
        time0=timeNow;
    else
        a=datevec(now);
        timeNow=a(6)+a(5)*60+a(4)*60*60; % calculate time at this instant
        dt=timeNow-time0; % calculate elapsed time interval
        time0=timeNow;
    end
 % Construct motion control command
     n=[joyStatus(3);
         joyStatus(4);
         0];
     k=[0;0;1];
     w_control=w*cross(n,k);
     v_control=v*[joyStatus(1);
         joyStatus(2);
         0];
% Calculate target transform using command, and current transform
     Tt=updateTransform(Tt,w_control,v_control,dt);

     qt=kukaDLSSolver( qin, Tt, TefTool,numberOfIterations,lambda );

     flag=jointLimitReached(qt);
     if flag==true;
         qt=qin;
     end

     % A condition to break the loop, when x button of the gamepad is
     % pressed
     b=button(joy,1);
     if(b==1)
         break;
     end

        qVec(:,1)=[];
        qVec=[qVec,columnVec(qt)];
        plot(tempVec,qVec);
        pause(0.001)
end

repo: KST-Kuka-Sunrise-Toolbox


text:
function ind = mohsst5_points_to_grid_index(C)

%data = mohsst5_loaddata();

lon = C(1,:);
lat = C(2,:);

N_lon = 72;
N_lat = 36;

% Random locations to 5x5-grid
ind_lon = floor(lon/5)+(N_lon/2+1);
ind_lat = N_lat - (floor(lat/5)+(N_lat/2+1)) + 1;

ind = (ind_lon-1)*N_lat + ind_lat;

author: jluttine

repo: matlab

save_path: github-repos/MATLAB/jluttine-matlab


text:
function [l,d,perm] = mchol(A,mu)
% Compute a modified LDL factorization of A
% (MEX ME!)

if nargin < 2
    mu = 1e-12;
end

n = size(A,1);
l = eye(n);
d = zeros(n,1);
perm = 1:n;

for i = 1:n
    c(i,i) = A(i,i);
end

% Compute modification parameters
gamma = max(abs(diag(A)));
xi = max(max(abs(setdiag(A,0))));
delta = mu*max(gamma+xi,1);
if n > 1
    beta = sqrt(max([gamma xi/sqrt(n^2-1) mu]));
else
    beta = sqrt(max([gamma mu]));
end

for j = 1:n

    % Find q that results in Best Permutation with j
    [maxVal maxPos] = max(abs(diag(c(j:end,j:end))));
    q = maxPos+j-1;

    % Permute d,c,l,a
    d([j q]) = d([q j]);
    perm([j q]) = perm([q j]);
    c([j q],:) = c([q j],:);
    c(:,[j q]) = c(:,[q j]);
    l([j q],:) = l([q j],:);
    l(:,[j q]) = l(:,[q j]);
    A([j q],:) = A([q j],:);
    A(:,[j q]) = A(:,[q j]);

    for s = 1:j-1
        l(j,s) = c(j,s)/d(s);
    end
    for i = j+1:n
        c(i,j) = A(i,j) - sum(l(j,1:j-1).*c(i,1:j-1));
    end
    theta = 0;
    if j < n && j > 1
        theta = max(abs(c(j+1:n,j)));
    end
    d(j) = max([abs(c(j,j)) (theta/beta)^2 delta]);
    if j < n
        for i = j+1:n
            c(i,i) = c(i,i) - (c(i,j)^2)/d(j);
        end
    end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [engine, clqtoroot] = push(engine, clq, pushdom)
%PUSH_POT push the variables in putshdom which is subset of clq to the clique toword the root and get new engine
%pushdom is pushed variables set
%clq is the index of the clique that pushdom belongs to

clqdom = engine.cliques{clq};
assert( mysubset(pushdom, clqdom));
clqtoroot = parents(engine.jtree, clq);
%sepdom = engine.separator{clq, clqtoroot};
sepdom = engine.separator{clqtoroot, clq};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the strong marginal of the union of pushdom and and the separatordomain and  %
% the corresponding complement                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%[margpot, comppot] = complement_pot(engine.clpot{clq}, pushdom);
newsepdom = myunion(pushdom,sepdom);
[margpot,comppot] = complement_pot(engine.clpot{clq}, newsepdom);
engine.clpot{clqtoroot} = direct_combine_pots(engine.clpot{clqtoroot}, margpot);
engine.clpot{clq} = comppot;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of the new separator and separatorpotential of the junction tree %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
engine.seppot{clqtoroot, clq} = direct_combine_pots(engine.seppot{clqtoroot, clq}, margpot);
engine.separator{clqtoroot, clq} = myunion(engine.separator{clqtoroot, clq}, pushdom);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Add pushdomain to the clique towards the root %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
engine.cliques{clqtoroot} = myunion(engine.cliques{clqtoroot}, pushdom);

num_cliques = length(engine.cliques);
B = sparse(num_cliques, 1);
for i=1:num_cliques
  B(i, engine.cliques{i}) = 1;
end
engine.cliques_bitv = B;

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function X = rocholTransMultiply(ch, Y);

% ROCHOLTRANSMULTIPLY Multiply by the transposed version of the rank one Cholesky.
% ROCHOL

%/~
% This would be the long way of doing it.
%L = rocholExtract(ch);
%X1 = L'*Y;
%~/
if size(Y, 1) ~= ch.n
  error('Inner matrix dimensions do not match');
end
X = zeros(size(Y));
t = zeros(1, size(Y, 2));
X(ch.n, :) = Y(ch.n, :)*ch.s(ch.n);
for i = ch.n-1:-1:1
  t = t + ch.v(i+1)*Y(i+1, :);
  X(i, :) = ch.s(i)*(Y(i, :)+ ch.u(i)*t);
end

%/~
%disp(max(max(X - X1)));
%~/

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
%% LSHAPESTOKES Poiseuille flow in a L-shaped domain
%
% Stokes equations on the square [-1,1]^2.
%
% Reference: page 237 in Finite Elements and Fast Iterative Solvers with
% Applications in Incompressible Fluid Dynamics. by Howard C. Elman, David
% J. Silvester, and Andrew J. Wathen.
%
% See also squareStokes, StokesP2P1
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

close all; clear all;
%% Parameters
maxIt = 4; N = zeros(maxIt,1);
erru = zeros(maxIt,1); errp = zeros(maxIt,1);

%% Generate an initial mesh
[node,elem] = squaremesh([-1 1 -1 1], 0.25);
bdFlag = setboundary(node,elem,'Dirichlet');

%% PDE and options
pde = Stokesdata1;

%% Finite Element Method
for k = 1:maxIt
    % solve the equation
    [u,p,edge,A] = StokesP2P1(node,elem,pde,bdFlag);
    N(k) = length(u)+length(p);
    if N(k) < 2e3 % show solution for small size
        figure(1);  showresult(node,elem,p);
    end
    % compute error
    uI = pde.exactu([node; (node(edge(:,1),:)+node(edge(:,2),:))/2]);
    erru(k) = sqrt((u-uI(:))'*A*(u-uI(:)));
    errp(k) = getL2error(node,elem,pde.exactp,p);
    % refine mesh
    [node,elem,bdFlag] = uniformrefine(node,elem,bdFlag);
end

%% Plot convergence rates
figure(2); clf;
showrate2(N,erru,2,'-*','||Du_I-Du_h||',N,errp,2,'k-+','|| p - p_h||');

repo: Aero-matlab


text:
function minVal=minOverDimIdx(C,selDim,selIdx)
%%MINOVERDIMIDX Given an S-dimensional matrix C, find the minimum value in
%       the matrix when the index in dimension selDim is fixed to selIdx.
%       For example, if C is 5 dimensional and selDim=4this function
%       evaluates the equivalent of min(vec(C(:,:,:,selIdx,:))). This
%       function can be useful when the dimensionality of the matrix in
%       question can vary.
%
%INPUTS: C An S-dimensional hypermatrix.
%   selDim The dimension of the hypermatrix that is selected.
%   selIdx The index in the selected dimension that is fixed.
%
%OUTPUTS: minVal The minimum value in the matrix when the specified
%                dimension is fixed to the specified index.
%
%For an arbitrary-dimensional matrix, the dimensions before the selected
%dimension can be collapsed into a single dimension and those after the
%selected dimension can also be collapsed into a single dimension. Thus,
%the problem reduces to the minimum of a 3D matrix where the middle index
%is fixed. In the even that one selects the first or last dimension, then
%the problem is the minimum of a 2D matrix where the first or last
%dimension is fixed.
%
%EXAMPLE:
%One will see that both of the value obtained by minOverDimIdx in this
%example is the same as the direct Matlab expression.
% C=randn(12,26,36,12,18);
% selDim=3;
% selIdx=8;
% minVal=minOverDimIdx(C,selDim,selIdx)
% min(vec(C(:,:,selIdx,:,:)))
%
%November 2020 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

nDims=size(C);
S=length(nDims);

if(S==1)
    %The special case of an array.
    minVal=C(selIdx);
elseif(selDim==1)
    startDim=nDims(1);
    endDim=prod(nDims(2:S));

    C=reshape(C,[startDim,endDim]);
    minVal=min(C(selIdx,:));
elseif(selDim==S)
    startDim=prod(nDims(1:(S-1)));
    endDim=nDims(S);

    C=reshape(C,[startDim,endDim]);
    minVal=min(C(:,selIdx));
else
    startDim=prod(nDims(1:(selDim-1)));
    endDim=prod(nDims((selDim+1):S));
    C=reshape(C,[startDim,nDims(selDim),endDim]);
    minVal=min(vec(C(:,selIdx,:)));
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function [X, funVal, ValueL]=LeastR(A, Y, z, opts)
%
%%
% Function LeastR
%      Least Squares Loss with the L1-norm Regularization
%
%% Problem
%
%  min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * ||x||_1
%
%  By default, rsL2=0.
%  When rsL2 is nonzero, this corresponds to the well-know elastic net.
%
%% Input parameters:
%
%  A-         Matrix of size m x n
%                A can be a dense matrix
%                         a sparse matrix
%                         or a DCT matrix
%  y -        Response vector (of size mx1)
%  z -        L_1 norm regularization parameter (z >=0)
%  opts-      Optional inputs (default value: opts=[])
%
%% Output parameters:
%  x-         Solution
%  funVal-    Function value during iterations
%
%% Copyright (C) 2009-2010 Jun Liu, and Jieping Ye
%
% You are suggested to first read the Manual.
%
% For any problem, please contact with Jun Liu via j.liu@asu.edu
%
% Last modified on February 18, 2010.
%
%% Related papers
%
% [1]  Jun Liu and Jieping Ye, Efficient Euclidean Projections
%      in Linear Time, ICML 2009.
%
% [2]  Jun Liu and Jieping Ye, Sparse Learning with Efficient Euclidean
%      Projections onto the L1 Ball, Technical Report ASU, 2008.
%
%% Related functions
%
%  sll_opts, initFactor, pathSolutionLeast
%  LeastC, nnLeastR, nnLeastC,
%
%%

%% Verify and initialize the parameters
%%

% Verify the number of input parameters
if (nargin <3)
    error('\n Inputs: A, y and z should be specified!\n');
elseif (nargin==3)
    opts=[];
end

% Get the size of the matrix A
[m,n]=size(A);% m -- sample size,n -- feature dimension%******************************
k=size(Y,2);
%sigma=diag(diag(rand(k)));
sigma=opts.sigma;
% Verify the length of y
if(size(Y,1) ~=m)%********************************
%if (length(y) ~=m)
    error('\n Check the length of y!\n');
end

% Verify the value of z
if (z<0)
    error('\n z should be nonnegative!\n');
end

% run sll_opts to set default values (flags)
opts=sll_opts(opts);

%% Detailed initialization
%% Normalization

% Please refer to sll_opts for the definitions of mu, nu and nFlag
%
% If .nFlag =1, the input matrix A is normalized to
%                     A= ( A- repmat(mu, m,1) ) * diag(nu)^{-1}
%
% If .nFlag =2, the input matrix A is normalized to
%                     A= diag(nu)^{-1} * ( A- repmat(mu, m,1) )
%
% Such normalization is done implicitly
%     This implicit normalization is suggested for the sparse matrix
%                                    but not for the dense matrix
%

if (opts.nFlag~=0)
    if (isfield(opts,'mu'))
        mu=opts.mu;
        if(size(mu,2)~=n)
            error('\n Check the input .mu');
        end
    else
        mu=mean(A,1);
    end

    if (opts.nFlag==1)
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=n)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,1)/m).^(0.5); nu=nu';
        end
    else % .nFlag=2
        if (isfield(opts,'nu'))
            nu=opts.nu;
            if(size(nu,1)~=m)
                error('\n Check the input .nu!');
            end
        else
            nu=(sum(A.^2,2)/n).^(0.5);
        end
    end

    ind_zero=find(abs(nu)<= 1e-10);    nu(ind_zero)=1;
    % If some values in nu is typically small, it might be that,
    % the entries in a given row or column in A are all close to zero.
    % For numerical stability, we set the corresponding value to 1.
end

if (~issparse(A)) && (opts.nFlag~=0)
    fprintf('\n -----------------------------------------------------');
    fprintf('\n The data is not sparse or not stored in sparse format');
    fprintf('\n The code still works.');
    fprintf('\n But we suggest you to normalize the data directly,');
    fprintf('\n for achieving better efficiency.');
    fprintf('\n -----------------------------------------------------');
end

%% Starting point initialization

% compute AT y
if (opts.nFlag==0)
    ATY=A'*Y; %A-->mxn;y-->mxk****************************
elseif (opts.nFlag==1)
    ATy=A'*y - sum(y) * mu';  ATy=ATy./nu;
else
    invNu=y./nu;              ATy=A'*invNu-sum(invNu)*mu';
end

% process the regularization parameter

% L2 norm regularization
if isfield(opts,'rsL2')
    rsL2=opts.rsL2;
    if (rsL2<0)
        error('\n opts.rsL2 should be nonnegative!');
    end
else
    rsL2=0;
end

% L1 norm regularization
if (opts.rFlag==0)
    lambda=z;
else % z here is the scaling factor lying in [0,1]
    if (z<0 || z>1)
        error('\n opts.rFlag=1, and z should be in [0,1]');
    end

    lambda_max=max(abs(ATy));
    lambda=z*lambda_max;

    rsL2=rsL2*lambda_max; % the input rsL2 is a ratio of lambda_max
end

% initialize a starting point
if opts.init==2
    X=zeros(n,k);%*************************************
    %x=zeros(n,1);
else
    if isfield(opts,'x0')
        x=opts.x0;
        if (length(x)~=n)
            error('\n Check the input .x0');
        end
    else
        X=ATY;  % if .x0 is not specified, we use ratio*ATy,
        % where ratio is a positive value
    end
end

% compute A x
if (opts.nFlag==0)
    AX=A* X;
elseif (opts.nFlag==1)
    invNu=x./nu; mu_invNu=mu * invNu;
    Ax=A*invNu -repmat(mu_invNu, m, 1);
else
    Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
end

if (opts.init==0) % If .init=0, we set x=ratio*x by "initFactor"
    % Please refer to the function initFactor for detail

    x_norm=sum(abs(X(:))); x_2norm=trace(X'*X);%*****************************corresponding to L1norm,Fnorm
    if x_norm>=1e-6
        ratio=initFactorqiao(x_norm, AX, Y, lambda,'LeastR', rsL2, x_2norm);
        X=ratio*X;    AX=ratio*AX;
    end
end

%% The main program

%% The Armijo Goldstein line search scheme + accelearted gradient descent
if (opts.mFlag==0 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    Xp=X; AXp=AX; XXp=zeros(n,k);

    % alphap and alpha are used for computing the weight in forming search point
    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    S=X + beta* XXp;

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        AS=AX + beta* (AX-AXp);

        % compute AT As
        if (opts.nFlag==0)
            ATAS=A'*AS;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        G=ATAS-ATY + rsL2 * S*inv(sigma);%************************** +rsL2*s*inv(sigma)

        % copy x and Ax to xp and Axp
        Xp=X;    AXp=AX;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the l1-norm regularized projection
            V=S-G/L;

            % L1-norm regularized projection
            X=sign(V).*max(abs(V)-lambda / L,0);

            V=X-S;  % the difference between the new approximate solution x
            % and the search point s

            % compute A x
            if (opts.nFlag==0)
                AX=A* X;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            AV=AX -AS;
            r_sum=trace(V'*V); l_sum=trace(AV'*AV);

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 <= (L - rsL2) * ||v||_2^2
            if(l_sum <= r_sum * (L-rsL2))
                break;
            else
                L=max(2*L, l_sum/r_sum + rsL2);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        XXp=X-Xp;   AXY=AX-Y;
        funVal(iterStep)=trace(AXY'* AXY)/2 + rsL2/2 * trace(X'*X) + sum(abs(X(:))) * lambda;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp);    norm_xxp=sqrt(xxp'*xxp);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%% Reformulated problem + Nemirovski's scheme

% .mFlag=1, and .lFlag=0
%  refomulate the problem as the constrained convex optimization
%  problem, and then apply Armijo Goldstein line search scheme

% Problem:
%    min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * t' * 1
%    s.t.   |x| <= t

if(opts.mFlag==1 && opts.lFlag==0)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

    % assign xp with x, and Axp with Ax
    xp=x; Axp=Ax; xxp=zeros(n,1);
    t=abs(x); tp=t;
    % t is the upper bound of absolute value of x

    % alphap and alpha are used for computing the weight in forming search point
    alphap=0; alpha=1;

    for iterStep=1:opts.maxIter
        % --------------------------- step 1 ---------------------------
        % compute search point s based on xp and x (with beta)
        beta=(alphap-1)/alpha;    s=x + beta* xxp; s_t= t + beta * (t -tp);

        % --------------------------- step 2 ---------------------------
        % line search for L and compute the new approximate solution x

        % compute the gradient (g) at s
        As=Ax + beta* (Ax-Axp);

        % compute AT As
        if (opts.nFlag==0)
            ATAs=A'*As;
        elseif (opts.nFlag==1)
            ATAs=A'*As - sum(As) * mu';  ATAs=ATAs./nu;
        else
            invNu=As./nu;                ATAs=A'*invNu-sum(invNu)*mu';
        end

        % obtain the gradient g
        g=ATAs-ATy + rsL2 * s;

        % copy x and Ax to xp and Axp
        xp=x;    Axp=Ax;
        tp=t;

        while (1)
            % let s walk in a step in the antigradient of s to get v
            % and then do the l1-norm regularized projection

            u=s-g/L;
            v= s_t - lambda / L;

            % projection
            [x, t]=ep1R(u, v, n);

            v=x-s;  % the difference between the new approximate solution x
                       % and the search point s

            v_t=t-s_t;

            % compute A x
            if (opts.nFlag==0)
                Ax=A* x;
            elseif (opts.nFlag==1)
                invNu=x./nu; mu_invNu=mu * invNu;
                Ax=A*invNu -repmat(mu_invNu, m, 1);
            else
                Ax=A*x-repmat(mu*x, m, 1);     Ax=Ax./nu;
            end

            Av=Ax -As;
            r_sum=v'*v + v_t'*v_t; l_sum=Av'*Av + v'*v * rsL2;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 + rsL2 * ||v||_2^2
            %                       <= L * (||v||_2^2 + ||v_t|| _2^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end

        ValueL(iterStep)=L;

        % --------------------------- step 3 ---------------------------
        % update alpha and alphap, and check whether converge
        alphap=alpha; alpha= (1+ sqrt(4*alpha*alpha +1))/2;

        xxp=x-xp;   Axy=Ax-y;
        funVal(iterStep)=Axy'* Axy/2 + rsL2/2 * x'*x + sum(t) * lambda;

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp + norm(t-tp)^2);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp + tp'*tp);    norm_xxp=sqrt(xxp'*xxp+ norm(t-tp)^2);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end

end

%% adaptive line search

% .mFlag=1, and .lFlag=1
%  refomulate the problem as the constrained convex optimization
%  problem, and then apply adaptive line search scheme

% Problem:
%    min  1/2 || A x - y||^2 + 1/2 rsL2 * ||x||_2^2 + z * t' * 1
%    s.t.   |x| <= t

if(opts.mFlag==1 && opts.lFlag==1)

    bFlag=0; % this flag tests whether the gradient step only changes a little

    L=1 + rsL2;
    % We assume that the maximum eigenvalue of A'A is over 1

    gamma=1;
    % we shall set the value of gamma = L,
    % where L is appropriate for the starting point

    xp=x; Axp=Ax;
    % store x and Ax
    xxp=zeros(n,1);
    % the difference of x and xp
    t=abs(x); tp=t;
    % t is the upper bound of absolute value of x

    % compute AT Ax
    if (opts.nFlag==0)
        ATAx=A'*Ax;
    elseif (opts.nFlag==1)
        ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
    else
        invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
    end

    % We begin the adaptive line search in the following
    %
    % Note that, in the line search, L and beta are changing

    for iterStep=1:opts.maxIter

        ATAxp=ATAx;
        % store ATAx to ATAxp

        if (iterStep~=1)
            % compute AT Ax
            if (opts.nFlag==0)
                ATAx=A'*Ax;
            elseif (opts.nFlag==1)
                ATAx=A'*Ax - sum(Ax) * mu';  ATAx=ATAx./nu;
            else
                invNu=Ax./nu;                ATAx=A'*invNu-sum(invNu)*mu';
            end
        end

        %--------- Line Search for L begins
        while (1)
            if (iterStep~=1)
                alpha= (-gamma+ sqrt(gamma*gamma + 4* L * gamma)) / (2*L);
                beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
                % beta is the coefficient for generating search point s

                s=x + beta* xxp;   s_t= t + beta * (t -tp);
                As=Ax + beta* (Ax-Axp);
                ATAs=ATAx + beta * (ATAx- ATAxp);
                % compute the search point s, A * s, and A' * A * s
            else
                alpha= (-1+ sqrt(5)) / 2;
                beta=0; s=x; s_t=t; As=Ax; ATAs=ATAx;
            end

            g=ATAs-ATy + rsL2 * s;
            % compute the gradient g

            % let s walk in a step in the antigradient of s
            u=s-g/L;
            v= s_t - lambda / L;

            % projection
            [xnew, tnew]=ep1R(u,v,n);

            v=xnew-s;  % the difference between the new approximate solution x
                            % and the search point s
            v_t=tnew-s_t;

            % compute A xnew
            if (opts.nFlag==0)
                Axnew=A* xnew;
            elseif (opts.nFlag==1)
                invNu=xnew./nu; mu_invNu=mu * invNu;
                Axnew=A*invNu -repmat(mu_invNu, m, 1);
            else
                Axnew=A*xnew-repmat(mu*xnew, m, 1);     Axnew=Axnew./nu;
            end

            Av=Axnew -As;
            r_sum=v'*v + v_t'*v_t; l_sum=Av'*Av + v'*v * rsL2;

            if (r_sum <=1e-20)
                bFlag=1; % this shows that, the gradient step makes little improvement
                break;
            end

            % the condition is ||Av||_2^2 + rsL2 * ||v||_2^2
            %                       <= L * (||v||_2^2 + ||v_t|| _2^2 )
            if(l_sum <= r_sum * L)
                break;
            else
                L=max(2*L, l_sum/r_sum);
                % fprintf('\n L=%5.6f',L);
            end
        end
        %--------- Line Search for L ends

        gamma=L* alpha* alpha;    alphap=alpha;
        % update gamma, and alphap

        ValueL(iterStep)=L;

        tao=L * r_sum / l_sum;
        if (tao >=5)
            L=L*0.8;
        end
        % decrease the value of L

        xp=x;    x=xnew; xxp=x-xp;
        Axp=Ax;  Ax=Axnew;
        % update x and Ax with xnew and Axnew
        tp=t; t=tnew;
        % update tp and t

        Axy=Ax-y;
        funVal(iterStep)=Axy' * Axy/2 + rsL2/2 * x'*x + lambda * sum(t);
        % compute function value

        if (bFlag)
            % fprintf('\n The program terminates as the gradient step changes the solution very small.');
            break;
        end

        switch(opts.tFlag)
            case 0
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <= opts.tol)
                        break;
                    end
                end
            case 1
                if iterStep>=2
                    if (abs( funVal(iterStep) - funVal(iterStep-1) ) <=...
                            opts.tol* funVal(iterStep-1))
                        break;
                    end
                end
            case 2
                if ( funVal(iterStep)<= opts.tol)
                    break;
                end
            case 3
                norm_xxp=sqrt(xxp'*xxp+ norm(t-tp)^2);
                if ( norm_xxp <=opts.tol)
                    break;
                end
            case 4
                norm_xp=sqrt(xp'*xp + tp'*tp);    norm_xxp=sqrt(xxp'*xxp+ norm(t-tp)^2);
                if ( norm_xxp <=opts.tol * max(norm_xp,1))
                    break;
                end
            case 5
                if iterStep>=opts.maxIter
                    break;
                end
        end
    end
end

%%
if(opts.mFlag==0 && opts.lFlag==1)
    error('\n The function does not support opts.mFlag=0 & opts.lFlag=1!');
end

author: zzstefan

repo: BrainNetClass


text:
classdef StiffnessMatrixComputer < handle

    properties (Access = private)
        LHS
        stiffnessMatrix
    end

    properties (Access = private)
        dim
        mesh
        youngModulus
        inertiaMoment
        freeNodes
    end

    methods (Access = public)

        function obj = StiffnessMatrixComputer(cParams)
            obj.init(cParams);
        end

        function compute(obj)
            obj.createStiffnessMatrix();
            obj.computeStiffnessMatrix()
        end

        function Kfree = provideFreeStiffnessMatrix(obj)
            free = obj.freeNodes;
            K = obj.stiffnessMatrix;
            Kfree  = K(free,free);
        end

    end

    methods (Access = private)

        function obj = init(obj,cParams)
            obj.mesh          = cParams.mesh;
            obj.dim           = cParams.dim;
            obj.freeNodes     = cParams.freeNodes;
        end

        function createStiffnessMatrix(obj)
            s.type = 'StiffnessMatrixColumn';
            s.dim = obj.dim;
            s.mesh = obj.mesh;
            s.globalConnec = obj.mesh.connec;
            obj.LHS = LHSintegrator.create(s);
        end

        function computeStiffnessMatrix(obj)
            obj.stiffnessMatrix = obj.LHS.compute();
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function val = get_params(CPD, name)
% GET_PARAMS Get the parameters (fields) for a softmax_CPD object
% val = get_params(CPD, name)
%
% The following fields can be accessed
%
% weights - W(X,Y,Q)
% offset  - b(Y,Q)
%
% e.g., W = get_params(CPD, 'weights')

[W, b] = extract_params(CPD);
switch name
 case 'weights',   val = W;
 case 'offset',    val = b;
 otherwise,
  error(['invalid argument name ' name]);
end


text:
function oeprint1(mu, oev)

% print six classical orbital elements
% and orbital period in minutes

% input

%  mu     = gravitational constant (km**3/sec**2)
%  oev(1) = semimajor axis (kilometers)
%  oev(2) = orbital eccentricity (non-dimensional)
%           (0 <= eccentricity < 1)
%  oev(3) = orbital inclination (radians)
%           (0 <= inclination <= pi)
%  oev(4) = argument of perigee (radians)
%           (0 <= argument of perigee <= 2 pi)
%  oev(5) = right ascension of ascending node (radians)
%           (0 <= raan <= 2 pi)
%  oev(6) = true anomaly (radians)
%           (0 <= true anomaly <= 2 pi)

% Orbital Mechanics with Matlab

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rtd = 180 / pi;

% unload orbital elements array

sma = oev(1);
ecc = oev(2);
inc = oev(3);
argper = oev(4);
raan = oev(5);
tanom = oev(6);

arglat = mod(tanom + argper, 2.0 * pi);

if (sma > 0.0)
    period = 2.0d0 * pi * sma * sqrt(sma / mu);
else
    period = 99999.9;
end

% print orbital elements

fprintf ('\n        sma (km)              eccentricity          inclination (deg)         argper (deg)');

fprintf ('\n %+16.14e  %+16.14e  %+16.14e  %+16.14e \n', sma, ecc, inc * rtd, argper * rtd);

fprintf ('\n       raan (deg)          true anomaly (deg)         arglat (deg)            period (min)');

fprintf ('\n %+16.14e  %+16.14e  %+16.14e  %+16.14e \n', raan * rtd, tanom * rtd, arglat * rtd, period / 60);

author: Sable

repo: mcbench-benchmarks


text:
function pred = knnpred(Xtest,X,class,K,dist_type,pret_type)

% prediction of new samples with calculated model
%
% pred = knnpred(Xtest,X,class,K,dist_type,pret_type)
%
% ------------ INPUT ---------------------------------------------------
% Xtest:        dataset to be predicted [n_test x p] n objects, p variables
% X:            training data matrix (n x p)
% class:        training class vector (n x 1)
% K:            number of neighbors
% dist_type:    'euclidean' Euclidean distance
%               'mahalanobis' Mahalanobis distance
%               'cityblock' City Block metric
%               'minkowski' Minkowski metric
%               'sm' Sokal-Michener
%               'jt' Jaccard Tanimoto
%               'gle' Gleason-Dice
%               'ct4' Consonni-Todeschini
%               'ac' Austin-Colwell
% pret_type:    'cent' cenering
%               'scal' variance scaling
%               'auto' for autoscaling (centering + variance scaling)
%               'rang' range scaling (0-1)
%               'fp'   fingerprints
%
% ------------ OUTPUT --------------------------------------------------
% pred is a structure conyaining
% class_pred    predicted class vector [n_test x 1]
% neighbors     list of k neighbors for each predicted sample [n_test x k]
%
% version 1.0 - september 2009
% Davide Ballabio
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% version 2.0 - February 2012
% Kamel Mansouri
% Milano Chemometrics and QSAR Research Group
% www.disat.unimib.it/chm

% data check
if length(class)~=size(X,1)
    disp('the class input should be for the training set')
    %class_tr=input('class tr');
    %class=evalin(WS,);
    %keyboard
end

[n,p] = size(Xtest);
[X_scal_train,param] = data_pretreatment(X,pret_type);
X_scal = test_pretreatment(Xtest,param);
Xd = [X_scal;X_scal_train];
% D = pdist(Xd,model.set.dist_type);
% D = squareform(D);
D = knn_calc_dist(X_scal_train,X_scal,dist_type,pret_type);
neighbors = zeros(n,K);
dc=zeros(n,K);
class_calc=zeros(n,1);
class_calc_weighted=zeros(n,1);
for i=1:n
    D_in = D(i,:);
    [d_tmp,n_tmp] = sort(D_in);
    neighbors(i,:) = n_tmp(1:K);
    d_neighbors = d_tmp(1:K);
    class_calc(i) = knnclass(class(neighbors(i,:)),d_neighbors,max(class),K);
    [yc(i),class_calc_weighted(i),w(i,:)] = nnrcalcy(class(neighbors(i,:)),d_neighbors,K);
    class_calc_weighted(i)=round(class_calc_weighted(i));
    dc(i,:)=d_neighbors;
end

pred.neighbors  = neighbors;
pred.class_pred = class_calc';
pred.class_pred_w = class_calc_weighted';
pred.w=w;
pred.D=D;
pred.dc=dc;

author: kmansouri

repo: OPERA

save_path: github-repos/MATLAB/kmansouri-OPERA


text:
function [out,options] = FCLSU_fast2(HIM,M,options)
% Fully Constrained Linear Spectral Unmixing
% Perform a Linear least squares with nonnegativity constraints.
% --------------------------------------------------------------------
% Input:   HIM : input data [nrows x nchannels]
%          M   : set of p endmembers [nchannels x p].
%
% Output:  out : fractions [nrows x p]
%
%
% Copyright (2007) GRNPS group @ University of Extremadura, Spain.
%
% *** Edited version ***

[ns,nb] = size(HIM);
[l,p] = size(M);

Delta = 1/1000; % should be an small value

N = zeros(l+1,p);
N(1:l,1:p) = Delta*M;
N(l+1,:) = ones(1,p);
s = zeros(l+1,1);

OutputImage = zeros(ns,p);

go=(nargin==3);

for i = 1:ns
    s(1:l) = Delta*HIM(i,:)';
    s(l+1) = 1;
    if go==0
        [Abundances,options] = lsqnonneg_fast(N,s);
        go=1;
    else
        Abundances = lsqnonneg_fast(N,s,options);
    end
    OutputImage(i,:) = Abundances;
end

out = OutputImage;

author: zhouyuanzxcv

repo: Hyperspectral


text:
%% SET
function obj = set(obj, varargin)
    %SET Set properties for EventRate object(s)
    %   obj = set(obj,'property_name', val, ['property_name2', val2])
    %   SET is one of the two gateway functions of an object, such as EventRate.
    %   Properties that are changed through SET are typechecked and otherwise
    %   scrutinized before being stored within the EventRate object.  This
    %   ensures that the other EventRate methods are all retrieving valid data,
    %   thereby increasing the reliability of the code.
    %
    %   Another strong advantage to using SET and GET to change and retrieve
    %   properties, rather than just assigning them to EventRate object directly,
    %   is that the underlying data structure can change and grow without
    %   harming the code that is written based on the EventRate object.
    %
    %   For a list of valid property names, type:
    %       properties(obj)
    %
    %   If user-defined fields were added to the EventRate object (ie, through
    %   addField), these fieldnames are also available through set.
    %
    %   Examples:
    %       (1) Change the description property
    %           obj = obj.set('description','hello world');
    %
    %       (2) Add new a field called CLOSEST_STATION with
    %           % the value 'MBLG'
    %           obj = obj.addfield('CLOSEST_STATION','MBLG');
    %
    %           % change the value of the CLOSEST_STATION field
    %           obj = obj.set('CLOSEST_STATION','MBWH');
    %
    %  See also EventRate/get, EventRate/addfield

    Vidx = 1 : numel(varargin);

    while numel(Vidx) >= 2
        prop_name = upper(varargin{Vidx(1)});
        val = varargin{Vidx(2)};
        mask = strcmp(upper(prop_name),upper(properties(obj)));
        if any(mask)
            mc = metaclass(obj);
            i = find(mask);
            prop_name = mc.PropertyList(i).Name;
            if isempty(mc.PropertyList(i).GetMethod)
                %eval(sprintf('obj.%s=val;',prop_name));
                obj.(prop_name) = val;
            else
                warning('Property %s is a derived property and cannot be set',prop_name);
            end
        else
            switch prop_name
                case obj.misc_fields
                    mask = strcmp(prop_name,obj.misc_fields);
                    obj.misc_values(mask) = {val};
                otherwise
                    error('%s:set:unknownProperty',...
                        'can''t understand property name : %s', mfilename,prop_name);
            end
        end
        Vidx(1:2) = []; %done with those parameters, move to the next ones...
    end
end

author: geoscience-community-codes

repo: GISMO


text:
% **********************Zip file 2: AF feature calculation
% AF feature
% calculation code package, which you can use to calculate the 14 AF
% features as the input feature metrix of SVM model. The AF 14 features
% inlcude the NFEn feature I sent you before.
%
% ***********************

clear all
close all
clc

data=load('RR_example.txt');
features = AF_features(data(7,1:53),250);

author: cliffordlab

repo: PhysioNet-Cardiovascular-Signal-Toolbox


text:
function volume = sphere_imp_volume_nd ( dim_num, r )

%*****************************************************************************80
%
%% SPHERE_IMP_VOLUME_ND computes the volume of an implicit sphere in ND.
%
%  Discussion:
%
%    An implicit sphere in ND satisfies the equation:
%
%      sum ( ( X(1:N) - CENTER(1:N) )^2 ) = R^2
%
%    where R is the radius and CENTER is the center.
%
%    Results for the first few values of N are:
%
%    DIM_NUM  Volume
%    -     -----------------------
%    2                PI   * R^2
%    3     (4/3)    * PI   * R^3
%    4     (1/2)    * PI^2 * R^4
%    5     (8/15)   * PI^2 * R^5
%    6     (1/6)    * PI^3 * R^6
%    7     (16/105) * PI^3 * R^7
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the space.
%
%    Input, real R, the radius of the sphere.
%
%    Output, real VOLUME, the volume of the sphere.
%
  volume = r^dim_num * sphere_unit_volume_nd ( dim_num );

  return
end

repo: jburkardt-m


text:
function w = demean (w)
   %DEMEAN remove offset voltage from signal
   %   waveform = demean(waveform)
   %   Removes the mean signal from the waveform object
   %
   %   Input Arguments
   %       WAVEFORM: a waveform object   N-DIMENSIONAL

   % AUTHOR: Celso Reyes, Geophysical Institute, Univ. of Alaska Fairbanks
   % $Date$
   % $Revision$

   m = mean(w);
   w = w - m;
end