# include <math.h>
# include <stdbool.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

# include "condition.h"
# include "r8lib.h"

/******************************************************************************/

double *combin ( double alpha, double beta, int n )

/******************************************************************************/
/*
  Purpose:

    COMBIN returns the COMBIN matrix.

  Discussion:

    This matrix is known as the combinatorial matrix.

  Formula:

    If ( I = J ) then
      A(I,J) = ALPHA + BETA
    else
      A(I,J) = BETA

  Example:

    N = 5, ALPHA = 2, BETA = 3

    5 3 3 3 3
    3 5 3 3 3
    3 3 5 3 3
    3 3 3 5 3
    3 3 3 3 5

  Properties:

    A is symmetric: A' = A.

    Because A is symmetric, it is normal.

    Because A is normal, it is diagonalizable.

    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).

    A is a circulant matrix: each row is shifted once to get the next row.

    det ( A ) = ALPHA^(N-1) * ( ALPHA + N * BETA ).

    A has constant row sums.

    Because A has constant row sums,
    it has an eigenvalue with this value,
    and a (right) eigenvector of ( 1, 1, 1, ..., 1 ).

    A has constant column sums.

    Because A has constant column sums,
    it has an eigenvalue with this value,
    and a (left) eigenvector of ( 1, 1, 1, ..., 1 ).

    LAMBDA(1:N-1) = ALPHA,
    LAMBDA(N) = ALPHA + N * BETA.

    The eigenvector associated with LAMBDA(N) is (1,1,1,...,1)/sqrt(N).

    The other N-1 eigenvectors are simply any (orthonormal) basis
    for the space perpendicular to (1,1,1,...,1).

    A is nonsingular if ALPHA /= 0 and ALPHA + N * BETA /= 0.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 June 2008

  Author:

    John Burkardt

  Reference:

    Robert Gregory, David Karney,
    A Collection of Matrices for Testing Computational Algorithms,
    Wiley, 1969,
    ISBN: 0882756494,
    LC: QA263.68

    Donald Knuth,
    The Art of Computer Programming,
    Volume 1, Fundamental Algorithms, Second Edition,
    Addison-Wesley, Reading, Massachusetts, 1973, page 36.

  Parameters:

    Input, double ALPHA, BETA, scalars that define A.

    Input, int N, the order of the matrix.

    Output, double COMBIN[N*N], the matrix.
*/
{
  double *a;
  int i;
  int j;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == j )
      {
        a[i+j*n] = alpha + beta;
      }
      else
      {
        a[i+j*n] = beta;
      }
    }
  }
  return a;
}
/******************************************************************************/

double *combin_inverse ( double alpha, double beta, int n )

/******************************************************************************/
/*
  Purpose:

    COMBIN_INVERSE returns the inverse of the COMBIN matrix.

  Formula:

    if ( I = J )
      A(I,J) = (ALPHA+(N-1)*BETA) / (ALPHA*(ALPHA+N*BETA))
    else
      A(I,J) =             - BETA / (ALPHA*(ALPHA+N*BETA))

  Example:

    N = 5, ALPHA = 2, BETA = 3

           14 -3 -3 -3 -3
           -3 14 -3 -3 -3
   1/34 *  -3 -3 14 -3 -3
           -3 -3 -3 14 -3
           -3 -3 -3 -3 14

  Properties:

    A is symmetric: A' = A.

    Because A is symmetric, it is normal.

    Because A is normal, it is diagonalizable.

    A is persymmetric: A(I,J) = A(N+1-J,N+1-I).

    A is a circulant matrix: each row is shifted once to get the next row.

    A is Toeplitz: constant along diagonals.

    det ( A ) = 1 / (ALPHA^(N-1) * (ALPHA+N*BETA)).

    A is well defined if ALPHA /= 0 and ALPHA+N*BETA /= 0.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 June 2008

  Author:

    John Burkardt

  Reference:

    Donald Knuth,
    The Art of Computer Programming,
    Volume 1, Fundamental Algorithms, Second Edition,
    Addison-Wesley, Reading, Massachusetts, 1973, page 36.

  Parameters:

    Input, double ALPHA, BETA, scalars that define the matrix.

    Input, int N, the order of the matrix.

    Output, double COMBIN_INVERSE[N*N], the matrix.
*/
{
  double *a;
  double bot;
  int i;
  int j;

  if ( alpha == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "COMBIN_INVERSE - Fatal error!\n" );
    fprintf ( stderr, "  The entries of the matrix are undefined\n" );
    fprintf ( stderr, "  because ALPHA = 0.\n" );
    exit ( 1 );
  }
  else if ( alpha + n * beta == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "COMBIN_INVERSE - Fatal error!\n" );
    fprintf ( stderr, "  The entries of the matrix are undefined\n" );
    fprintf ( stderr, "  because ALPHA+N*BETA is zero.\n" );
    exit ( 1 );
  }

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  bot = alpha * ( alpha + ( double ) ( n ) * beta );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == j )
      {
        a[i+j*n] = ( alpha + ( double ) ( n - 1 ) * beta ) / bot;
      }
      else
      {
        a[i+j*n] = - beta / bot;
      }
    }
  }
  return a;
}
/******************************************************************************/

double condition_hager ( int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    CONDITION_HAGER estimates the L1 condition number of a matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    12 April 2012

  Author:

    John Burkardt

  Reference:

    William Hager,
    Condition Estimates,
    SIAM Journal on Scientific and Statistical Computing,
    Volume 5, Number 2, June 1984, pages 311-316.

  Parameters:

    Input, int N, the dimension of the matrix.

    Input, double A[N*N], the matrix.

    Output, double CONDITION_HAGER, the estimated L1 condition.
*/
{
  double *a_lu;
  double *b;
  double c1;
  double c2;
  double cond;
  int i;
  int i1;
  int i2;
  int job;
  int *pivot;

  i1 = -1;
  c1 = 0.0;
/*
  Factor the matrix.
*/
  a_lu = r8mat_copy_new ( n, n, a );

  pivot = ( int * ) malloc ( n * sizeof ( int ) );

  r8ge_fa ( n, a_lu, pivot );

  b = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    b[i] = 1.0 / ( double ) ( n );
  }

  while ( 1 )
  {
    job = 0;
    r8ge_sl ( n, a_lu, pivot, b, job );

    c2 = r8vec_norm_l1 ( n, b );

    for ( i = 0; i < n; i++ )
    {
      b[i] = r8_sign ( b[i] );
    }

    job = 1;
    r8ge_sl ( n, a_lu, pivot, b, job );

    i2 = r8vec_max_abs_index ( n, b );

    if ( 0 <= i1 )
    {
      if ( i1 == i2 || c2 <= c1 )
      {
        break;
      }
    }
    i1 = i2;
    c1 = c2;

    for ( i = 0; i < n; i++ )
    {
      b[i] = 0.0;
    }
    b[i1] = 1.0;
  }
  cond = c2 * r8mat_norm_l1 ( n, n, a );

  free ( a_lu );
  free ( b );
  free ( pivot );

  return cond;
}
/******************************************************************************/

double condition_linpack ( int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    CONDITION_LINPACK estimates the L1 condition number.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

    For the system A * X = B, relative perturbations in A and B
    of size EPSILON may cause relative perturbations in X of size
    EPSILON * COND.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    09 April 2012

  Author:

    Original FORTRAN77 version by Dongarra, Bunch, Moler, Stewart.
    C++ version by John Burkardt.

  Reference:

    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
    LINPACK User's Guide,
    SIAM, 1979,
    ISBN13: 978-0-898711-72-1,
    LC: QA214.L56.

  Parameters:

    Input, int N, the order of the matrix A.

    Input/output, double A[N*N].  On input, a matrix to be factored.
    On output, the LU factorization of the matrix.

    Output, double CONDITION_LINPACK, the estimated L1 condition.
*/
{
  double anorm;
  double cond;
  double ek;
  int i;
  int info;
  int j;
  int k;
  int l;
  int *pivot;
  double s;
  double sm;
  double t;
  double wk;
  double wkm;
  double ynorm;
  double *z;
/*
  Compute the L1 norm of A.
*/
  anorm = 0.0;
  for ( j = 0; j < n; j++ )
  {
    s = 0.0;
    for ( i = 0; i < n; i++ )
    {
      s = s + fabs ( a[i+j*n] );
    }
    anorm = fmax ( anorm, s );
  }
/*
  Compute the LU factorization.
*/
  pivot = ( int * ) malloc ( n * sizeof ( int ) );

  info = r8ge_fa ( n, a, pivot );

  if ( info != 0 ) 
  {
    free ( pivot );
    cond = 0.0;
    return cond;
  }
/*
  COND = norm(A) * (estimate of norm(inverse(A)))

  estimate of norm(inverse(A)) = norm(Z) / norm(Y)

  where
    A * Z = Y
  and
    A' * Y = E

  The components of E are chosen to cause maximum local growth in the
  elements of W, where U'*W = E.  The vectors are frequently rescaled
  to avoid overflow.

  Solve U' * W = E.
*/
  ek = 1.0;
  z = ( double * ) malloc ( n * sizeof ( double ) );
  for ( i = 0; i < n; i++ )
  {
    z[i] = 0.0;
  }

  for ( k = 0; k < n; k++ )
  {
    if ( z[k] != 0.0 ) 
    {
      ek = - r8_sign2 ( ek, z[k] );
    }

    if ( fabs ( a[k+k*n] ) < fabs ( ek - z[k] ) )
    {
      s = fabs ( a[k+k*n] ) / fabs ( ek - z[k] );
      for ( i = 0; i < n; i++ )
      {
        z[i] = s * z[i];
      }
      ek = s * ek;
    }

    wk = ek - z[k];
    wkm = -ek - z[k];
    s = fabs ( wk );
    sm = fabs ( wkm );

    if ( a[k+k*n] != 0.0 )
    {
      wk = wk / a[k+k*n];
      wkm = wkm / a[k+k*n];
    }
    else
    {
      wk = 1.0;
      wkm = 1.0;
    }

    if ( k + 2 <= n )
    {
      for ( j = k+1; j < n; j++ )
      {
        sm = sm + fabs ( z[j] + wkm * a[k+j*n] );
        z[j] = z[j] + wk * a[k+j*n];
        s = s + fabs ( z[j] );
      }

      if ( s < sm )
      {
        t = wkm - wk;
        wk = wkm;
        for ( j = k+1; j < n; j++ )
        {
          z[j] = z[j] + t * a[k+j*n];
        }
      }
    }
    z[k] = wk;
  }

  s = 0.0;
  for ( i = 0; i < n; i++ )
  {
    s = s + fabs ( z[i] );
  }

  for ( i = 0; i < n; i++ )
  {
    z[i] = z[i] / s;
  }
/*
  Solve L' * Y = W
*/
  for ( k = n-1; 0 <= k; k-- )
  {
    for ( i = k+1; i < n; i++ )
    {
      z[k] = z[k] + z[i] * a[i+k*n];
    }
    t = fabs ( z[k] );
    if ( 1.0 < t )
    {
      for ( i = 0; i < n; i++ )
      {
        z[i] = z[i] / t;
      }
    }

    l = pivot[k] - 1;

    t    = z[l];
    z[l] = z[k];
    z[k] = t;
  }

  t = 0.0;
  for ( i = 0; i < n; i++ )
  {
    t = t + fabs ( z[i] );
  }
  for ( i = 0; i < n; i++ )
  {
    z[i] = z[i] / t;
  }

  ynorm = 1.0;
/*
  Solve L * V = Y.
*/
  for ( k = 0; k < n; k++ )
  {
    l = pivot[k] - 1;

    t    = z[l];
    z[l] = z[k];
    z[k] = t;

    for ( i = k+1; i < n; i++ )
    {
      z[i] = z[i] + t * a[i+k*n];
    }

    t = fabs ( z[k] );

    if ( 1.0 < t )
    {
      ynorm = ynorm / t;
      for ( i = 0; i < n; i++ )
      {
        z[i] = z[i] / t;
      }
    }
  }
  s = 0.0;
  for ( i = 0; i < n; i++ )
  {
    s = s + fabs ( z[i] );
  }

  for ( i = 0; i < n; i++ )
  {
    z[i] = z[i] / s;
  }
  ynorm = ynorm / s;
/*
  Solve U * Z = V.
*/
  for ( k = n-1; 0 <= k; k-- )
  {
    if ( fabs ( a[k+k*n] ) < fabs ( z[k] ) )
    {
      s = fabs ( a[k+k*n] ) / fabs ( z[k] );
      for ( i = 0; i < n; i++ )
      {
        z[i] = s * z[i];
      }
      ynorm = s * ynorm;
    }

    if ( a[k+k*n] != 0.0 )
    {
      z[k] = z[k] / a[k+k*n];
    }
    else
    {
      z[k] = 1.0;
    }

    for ( i = 0; i < k; i++ )
    {
      z[i] = z[i] - a[i+k*n] * z[k];
    }
  }
/*
  Normalize Z in the L1 norm.
*/
  s = 0.0;
  for ( i = 0; i < n; i++ )
  {
    s = s + fabs ( z[i] );
  }
  s = 1.0 / s;

  for ( i = 0; i < n; i++ )
  {
    z[i] = s * z[i];
  }
  ynorm = s * ynorm;

  cond = anorm / ynorm;

  free ( pivot );
  free ( z );

  return cond;
}
/******************************************************************************/

double condition_sample1 ( int n, double a[], int m )

/******************************************************************************/
/*
  Purpose:

    CONDITION_SAMPLE1 estimates the L1 condition number of a matrix.

  Discussion:

    A naive sampling method is used.

    Only "forward" sampling is used, that is, we only look at results
    of the form y=A*x.

    Presumably, solving systems A*y=x would give us a better idea of 
    the inverse matrix.

    Moreover, a power sequence y1 = A*x, y2 = A*y1, ... and the same for
    the inverse might work better too.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    09 April 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the dimension of the matrix.

    Input, double A[N*N], the matrix.

    Input, int M, the number of samples to use.

    Output, double CONDITION_SAMPLE1, the estimated L1 condition.
*/
{
  double a_norm;
  double ainv_norm;
  double *ax;
  double ax_norm;
  double cond;
  int i;
  int seed;
  double *x;
  double x_norm;

  a_norm = 0.0;
  ainv_norm = 0.0;
  seed = 123456789;

  for ( i = 1; i <= m; i++ )
  {
    x = r8vec_uniform_unit_new ( n, &seed );
    x_norm = r8vec_norm_l1 ( n, x );
    ax = r8mat_mv_new ( n, n, a, x );
    ax_norm = r8vec_norm_l1 ( n, ax );

    if ( ax_norm == 0.0 )
    {
      cond = 0.0;
      return cond;
    }

    a_norm    = fmax ( a_norm,    ax_norm / x_norm  );
    ainv_norm = fmax ( ainv_norm, x_norm  / ax_norm );

    free ( ax );
    free ( x );
  }

  cond = a_norm * ainv_norm;

  return cond;
}
/******************************************************************************/

double *conex1 ( double alpha )

/******************************************************************************/
/*
  Purpose:

    CONEX1 returns the CONEX1 matrix.

  Discussion:

    The CONEX1 matrix is a counterexample to the LINPACK condition
    number estimator RCOND available in the LINPACK routine DGECO.

  Formula:

    1  -1 -2*ALPHA   0
    0   1    ALPHA    -ALPHA
    0   1  1+ALPHA  -1-ALPHA
    0   0  0           ALPHA

  Example:

    ALPHA = 100

    1  -1  -200     0
    0   1   100  -100
    0   1   101  -101
    0   0     0   100

  Properties:

    A is generally not symmetric: A' /= A.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 August 2008

  Author:

    John Burkardt

  Reference:

    Alan Cline, RK Rew,
    A set of counterexamples to three condition number estimators,
    SIAM Journal on Scientific and Statistical Computing,
    Volume 4, 1983, pages 602-611.

  Parameters:

    Input, double ALPHA, the scalar defining A.  
    A common value is 100.0.

    Output, double CONEX1[4*4], the matrix.
*/
{
  double *a;
  int n = 4;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  a[0+0*n] = 1.0;
  a[1+0*n] = 0.0;
  a[2+0*n] = 0.0;
  a[3+0*n] = 0.0;

  a[0+1*n] = -1.0;
  a[1+1*n] = 1.0;
  a[2+1*n] = 1.0;
  a[3+1*n] = 0.0;

  a[0+2*n] = -2.0 * alpha;
  a[1+2*n] = alpha;
  a[2+2*n] = 1.0 + alpha;
  a[3+2*n] = 0.0;

  a[0+3*n] = 0.0;
  a[1+3*n] = -alpha;
  a[2+3*n] = -1.0 - alpha;
  a[3+3*n] = alpha;

  return a;
}
/******************************************************************************/

double *conex1_inverse ( double alpha )

/******************************************************************************/
/*
  Purpose:

    CONEX1_INVERSE returns the inverse of the CONEX1 matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 September 2008

  Author:

    John Burkardt

  Parameters:

    Input, double ALPHA, the scalar defining A.  

    Output, double CONEX1_INVERSE[4*4], the matrix.
*/
{
  double *a;
  int n = 4;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  a[0+0*n] =  1.0;
  a[1+0*n] =  0.0;
  a[2+0*n] =  0.0;
  a[3+0*n] =  0.0;

  a[0+1*n] =  1.0 - alpha;
  a[1+1*n] =  1.0 + alpha;
  a[2+1*n] = -1.0;
  a[3+1*n] =  0.0;

  a[0+2*n] =        alpha;
  a[1+2*n] =      - alpha;
  a[2+2*n] =  1.0;
  a[3+2*n] =  0.0;

  a[0+3*n] =  2.0;
  a[1+3*n] =  0.0;
  a[2+3*n] =  1.0 / alpha;
  a[3+3*n] =  1.0 / alpha;

  return a;
}
/******************************************************************************/

double *conex2 ( double alpha )

/******************************************************************************/
/*
  Purpose:

    CONEX2 returns the CONEX2 matrix.

  Formula:

    1   1-1/ALPHA^2  -2
    0   1/ALPHA      -1/ALPHA
    0   0             1

  Example:

    ALPHA = 100

    1  0.9999  -2
    0  0.01    -0.01
    0  0        1

  Properties:

    A is generally not symmetric: A' /= A.

    A is upper triangular.

    det ( A ) = 1 / ALPHA.

    LAMBDA = ( 1, 1/ALPHA, 1 )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 September 2008

  Author:

    John Burkardt

  Reference:

    Alan Cline, RK Rew,
    A set of counterexamples to three condition number estimators,
    SIAM Journal on Scientific and Statistical Computing,
    Volume 4, 1983, pages 602-611.

  Parameters:

    Input, double ALPHA, the scalar defining A.  
    A common value is 100.0.  ALPHA must not be zero.

    Output, double CONEX2[3*3], the matrix.
*/
{
  double *a;
  int n = 3;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  if ( alpha == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "CONEX2 - Fatal error!\n" );
    fprintf ( stderr, "  The input value of ALPHA was zero.\n" );
    exit ( 1 );
  }

  a[0+0*n] =  1.0;
  a[1+0*n] =  0.0;
  a[2+0*n] =  0.0;

  a[0+1*n] = ( alpha - 1.0 ) * ( alpha + 1.0 ) / alpha / alpha;
  a[1+1*n] =  1.0 / alpha;
  a[2+1*n] =  0.0;

  a[0+2*n] = -2.0;
  a[1+2*n] = -1.0 / alpha;
  a[2+2*n] =  1.0;

  return a;
}
/******************************************************************************/

double *conex2_inverse ( double alpha )

/******************************************************************************/
/*
  Purpose:

    CONEX2_INVERSE returns the inverse of the CONEX2 matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 September 2008

  Author:

    John Burkardt

  Parameters:

    Input, double ALPHA, the scalar defining A.  
    A common value is 100.0.  ALPHA must not be zero.

    Output, double CONEX2_INVERSE[3*3], the matrix.
*/
{
  double *a;
  int n = 3;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  if ( alpha == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "CONEX2_INVERSE - Fatal error!\n" );
    fprintf ( stderr, "  The input value of ALPHA was zero.\n" );
    exit ( 1 );
  }

  a[0+0*n] = 1.0;
  a[1+0*n] = 0.0;
  a[2+0*n] = 0.0;

  a[0+1*n] = ( 1.0 - alpha ) * ( 1.0 + alpha ) / alpha;
  a[1+1*n] = alpha;
  a[2+1*n] = 0.0;

  a[0+2*n] = ( 1.0 + alpha * alpha ) / alpha / alpha;
  a[1+2*n] = 1.0;
  a[2+2*n] = 1.0;

  return a;
}
/******************************************************************************/

double *conex3 ( int n )

/******************************************************************************/
/*
  Purpose:

    CONEX3 returns the CONEX3 matrix.

  Formula:

    if ( I = J and I < N )
      A(I,J) =  1.0 for 1<=I<N
    else if ( I = J = N )
      A(I,J) = -1.0
    else if ( J < I )
      A(I,J) = -1.0
    else
      A(I,J) =  0.0

  Example:

    N = 5

     1  0  0  0  0
    -1  1  0  0  0
    -1 -1  1  0  0
    -1 -1 -1  1  0
    -1 -1 -1 -1 -1

  Properties:

    A is generally not symmetric: A' /= A.

    A is integral, therefore det ( A ) is integral, and 
    det ( A ) * inverse ( A ) is integral.

    A is lower triangular.

    det ( A ) = -1.

    A is unimodular.

    LAMBDA = ( 1, 1, 1, 1, -1 )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 July 2008

  Author:

    John Burkardt

  Reference:

    Alan Cline, RK Rew,
    A set of counterexamples to three condition number estimators,
    SIAM Journal on Scientific and Statistical Computing,
    Volume 4, 1983, pages 602-611.

  Parameters:

    Input, int N, the order of the matrix.

    Output, double CONEX3[N*N], the matrix.
*/
{
  double *a;
  int i;
  int j;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( j < i )
      {
        a[i+j*n] = - 1.0;
      }
      else if ( j == i && i != n - 1 )
      {
        a[i+j*n] = 1.0;
      }
      else if ( j == i && i == n - 1 )
      {
        a[i+j*n] = - 1.0;
      }
      else
      {
        a[i+j*n] = 0.0;
      }
    }
  }
  return a;
}
/******************************************************************************/

double *conex3_inverse ( int n )

/******************************************************************************/
/*
  Purpose:

    CONEX3_INVERSE returns the inverse of the CONEX3 matrix.

  Example:

    N = 5

     1  0  0  0  0
     1  1  0  0  0
     2  1  1  0  0
     4  2  1  1  0
    -8 -4 -2 -1 -1

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 July 2008

  Author:

    John Burkardt

  Reference:

    Alan Cline, RK Rew,
    A set of counterexamples to three condition number estimators,
    SIAM Journal on Scientific and Statistical Computing,
    Volume 4, 1983, pages 602-611.

  Parameters:

    Input, int N, the order of the matrix.

    Output, double CONEX3_INVERSE[N*N], the matrix.
*/
{
  double *a;
  int i;
  int j;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i < n - 1 )
      {
        if ( j < i )
        {
          a[i+j*n] = pow ( 2.0, i - j - 1 );
        }
        else if ( i == j )
        {
          a[i+j*n] = 1.0;
        }
        else
        {
          a[i+j*n] = 0.0;
        }
	  }
      else if ( i == n - 1 )
      {
        if ( j < i )
        {
          a[i+j*n] = - pow ( 2.0, i - j - 1 );
        }
        else
        {
          a[i+j*n] = - 1.0;
        }
      }
    }
  }
  return a;
}
/******************************************************************************/

double *conex4 ( )

/******************************************************************************/
/*
  Purpose:

    CONEX4 returns the CONEX4 matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 April 2012

  Author:

    John Burkardt

  Parameters:

    Output, double CONEX4[4*4], the matrix.
*/
{
  double *a;
  static double a_save[4*4] = {
     7.0,  6.0,  5.0,  5.0,
    10.0,  8.0,  7.0,  7.0,
     8.0, 10.0,  9.0,  6.0,
     7.0,  9.0, 10.0,  5.0 };

  a = r8mat_copy_new ( 4, 4, a_save );

  return a;
}
/******************************************************************************/

double *conex4_inverse ( )

/******************************************************************************/
/*
  Purpose:

    CONEX4_INVERSE returns the inverse of the CONEX4 matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 April 2012

  Author:

    John Burkardt

  Parameters:

    Output, double CONEX4_INVERSE[4*4], the matrix.
*/
{
  double *a;
  static double a_save[4*4] = {
   -41.0,  25.0,  10.0, -6.0,
   -17.0,  10.0,   5.0, -3.0,
    10.0,  -6.0,  -3.0,  2.0,
    68.0, -41.0, -17.0, 10.0 };

  a = r8mat_copy_new ( 4, 4, a_save );

  return a;
}
/******************************************************************************/

double *kahan ( double alpha, int m, int n )

/******************************************************************************/
/*
  Purpose:

    KAHAN returns the KAHAN matrix.

  Formula:

    if ( I = J )
      A(I,I) =  sin(ALPHA)^(I)
    else if ( I < J )
      A(I,J) = - sin(ALPHA)^(I) * cos(ALPHA)
    else
      A(I,J) = 0

  Example:

    ALPHA = 0.25, N = 4

    S  -C*S    -C*S      -C*S
    0     S^2  -C*S^2    -C*S^2
    0     0       S^3    -C*S^3
    0     0       0         S^4

    where

      S = sin(ALPHA), C=COS(ALPHA)

  Properties:

    A is upper triangular.

    A = B * C, where B is a diagonal matrix and C is unit upper triangular.
    For instance, for the case M = 3, N = 4:

    A = | S 0    0    |  * | 1 -C -C  -C |
        | 0 S^2  0    |    | 0  1 -C  -C |
        | 0 0    S^3  |    | 0  0  1  -C |

    A is generally not symmetric: A' /= A.

    A has some interesting properties regarding estimation of
    condition and rank.

    det ( A ) = sin(ALPHA^(N*(N+1)/2).

    LAMBDA(I) = sin ( ALPHA )^I

    A is nonsingular if and only if sin ( ALPHA ) =/= 0.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    27 June 2011

  Author:

    John Burkardt

  Reference:

    Nicholas Higham,
    A survey of condition number estimation for triangular matrices,
    SIAM Review,
    Volume 9, 1987, pages 575-596.

    W Kahan,
    Numerical Linear Algebra,
    Canadian Mathematical Bulletin,
    Volume 9, 1966, pages 757-801.

  Parameters:

    Input, double ALPHA, the scalar that defines A.  A typical
    value is 1.2.  The "interesting" range of ALPHA is 0 < ALPHA < PI.

    Input, int M, N, the order of the matrix.

    Output, double KAHAN[M*N], the matrix.
*/
{
  double *a;
  double csi;
  int i;
  int j;
  double si;

  a = ( double * ) malloc ( m * n * sizeof ( double ) );

  for ( i = 0; i < m; i++ )
  {
    si = pow ( sin ( alpha ), i + 1 );
    csi = - cos ( alpha ) * si;
    for ( j = 0; j < n; j++ )
    {
      if ( j < i )
      {
        a[i+j*m] = 0.0;
      }
      else if ( j == i )
      {
        a[i+j*m] = si;
      }
      else
      {
        a[i+j*m] = csi;
      }
    }
  }
  return a;
}
/******************************************************************************/

double *kahan_inverse ( double alpha, int n )

/******************************************************************************/
/*
  Purpose:

    KAHAN_INVERSE returns the inverse of the KAHAN matrix.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    29 June 2008

  Author:

    John Burkardt

  Parameters:

    Input, double ALPHA, the scalar that defines A.  A typical 
    value is 1.2.  The "interesting" range of ALPHA is 0 < ALPHA < PI.

    Input, int N, the order of the matrix.

    Output, double KAHAN_INVERSE[N*N], the matrix.
*/
{
  double *a;
  double ci;
  int i;
  int j;
  double si;

  a = ( double * ) malloc ( n * n * sizeof ( double ) );

  ci = cos ( alpha );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == j )
      {
        a[i+j*n] = 1.0;
      }
      else if ( i == j - 1 )
      {
        a[i+j*n] = ci;
      }
      else if ( i < j )
      {
        a[i+j*n] = ci * pow ( 1.0 + ci, j - i - 1 );
      }
      else
      {
        a[i+j*n] = 0.0;
      }
    }
  }
/*
  Scale the columns.
*/
  for ( j = 0; j < n; j++ )
  {
    si = pow ( sin ( alpha ), j + 1 );
    for ( i = 0; i < n; i++ )
    {
      a[i+j*n] = a[i+j*n] / si;
    }
  }
  return a;
}
/******************************************************************************/

int r8ge_fa ( int n, double a[], int pivot[] )

/******************************************************************************/
/*
  Purpose:

    R8GE_FA performs a LINPACK-style PLU factorization of a R8GE matrix.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

    R8GE_FA is a simplified version of the LINPACK routine SGEFA.

    The two dimensional array is stored by columns in a one dimensional
    array.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    11 September 2003

  Author:

    John Burkardt

  Reference:

    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
    LINPACK User's Guide,
    SIAM, 1979,
    ISBN13: 978-0-898711-72-1,
    LC: QA214.L56.

  Parameters:

    Input, int N, the order of the matrix.
    N must be positive.

    Input/output, double A[N*N], the matrix to be factored.
    On output, A contains an upper triangular matrix and the multipliers
    which were used to obtain it.  The factorization can be written
    A = L * U, where L is a product of permutation and unit lower
    triangular matrices and U is upper triangular.

    Output, int PIVOT[N], a vector of pivot indices.

    Output, int R8GE_FA, singularity flag.
    0, no singularity detected.
    nonzero, the factorization failed on the INFO-th step.
*/
{
  int i;
  int j;
  int k;
  int l;
  double t;

  for ( k = 1; k <= n-1; k++ )
  {
/*
  Find L, the index of the pivot row.
*/
    l = k;

    for ( i = k+1; i <= n; i++ )
    {
      if ( fabs ( a[l-1+(k-1)*n] ) < fabs ( a[i-1+(k-1)*n] ) )
      {
        l = i;
      }
    }

    pivot[k-1] = l;
/*
  If the pivot index is zero, the algorithm has failed.
*/
    if ( a[l-1+(k-1)*n] == 0.0 )
    {
      printf ( "\n" );
      printf ( "R8GE_FA - Fatal error!\n" );
      printf ( "  Zero pivot on step %d\n", k );
      exit ( 1 );
    }
/*
  Interchange rows L and K if necessary.
*/
    if ( l != k )
    {
      t              = a[l-1+(k-1)*n];
      a[l-1+(k-1)*n] = a[k-1+(k-1)*n];
      a[k-1+(k-1)*n] = t;
    }
/*
  Normalize the values that lie below the pivot entry A(K,K).
*/
    for ( i = k+1; i <= n; i++ )
    {
      a[i-1+(k-1)*n] = -a[i-1+(k-1)*n] / a[k-1+(k-1)*n];
    }
/*
  Row elimination with column indexing.
*/
    for ( j = k+1; j <= n; j++ )
    {
      if ( l != k )
      {
        t              = a[l-1+(j-1)*n];
        a[l-1+(j-1)*n] = a[k-1+(j-1)*n];
        a[k-1+(j-1)*n] = t;
      }

      for ( i = k+1; i <= n; i++ )
      {
        a[i-1+(j-1)*n] = a[i-1+(j-1)*n] + a[i-1+(k-1)*n] * a[k-1+(j-1)*n];
      }

    }

  }

  pivot[n-1] = n;

  if ( a[n-1+(n-1)*n] == 0.0 )
  {
    printf ( "\n" );
    printf ( "R8GE_FA - Fatal error!\n" );
    printf ( "  Zero pivot on step %d\n", n );
    exit ( 1 );
  }

  return 0;
}
/******************************************************************************/

double *r8ge_inverse ( int n, double a[], int pivot[] )

/******************************************************************************/
/*
  Purpose:

    R8GE_INVERSE computes the inverse of a R8GE matrix factored by R8GE_FA.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

    R8GE_INVERSE is a simplified standalone version of the LINPACK routine
    SGEDI.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    15 September 2003

  Author:

    John Burkardt

  Parameters:

    Input, int N, the order of the matrix A.

    Input, double A[N*N], the factor information computed by R8GE_FA.

    Input, int PIVOT(N), the pivot vector from R8GE_FA.

    Output, double R8GE_INVERSE[N*N], the inverse matrix.
*/
{
  double *b;
  int i;
  int j;
  int k;
  double temp;

  b = ( double * ) malloc ( n * n * sizeof ( double ) );
/*
  Compute Inverse(U).
*/
  for ( k = 1; k <= n; k++ )
  {
    for ( i = 1; i <= k-1; i++ )
    {
      b[i-1+(k-1)*n] = -b[i-1+(k-1)*n] / a[k-1+(k-1)*n];
    }
    b[k-1+(k-1)*n] = 1.0 / a[k-1+(k-1)*n];

    for ( j = k+1; j <= n; j++ )
    {
      b[k-1+(j-1)*n] = 0.0;
      for ( i = 1; i <= k; i++ )
      {
        b[i-1+(j-1)*n] = b[i-1+(j-1)*n] + b[i-1+(k-1)*n] * a[k-1+(j-1)*n];
      }
    }
  }
/*
  Multiply Inverse(U) by Inverse(L).
*/
  for ( k = n-1; 1 <= k; k-- )
  {
    for ( i = k+1; i <= n; i++ )
    {
      b[i-1+(k-1)*n] = 0.0;
    }

    for ( j = k+1; j <= n; j++ )
    {
      for ( i = 1; i <= n; i++ )
      {
        b[i-1+(k-1)*n] = b[i-1+(k-1)*n] + b[i-1+(j-1)*n] * a[j-1+(k-1)*n];
      }
    }

    if ( pivot[k-1] != k )
    {
      for ( i = 1; i <= n; i++ )
      {
        temp = b[i-1+(k-1)*n];
        b[i-1+(k-1)*n] = b[i-1+(pivot[k-1]-1)*n];
        b[i-1+(pivot[k-1]-1)*n] = temp;
      }

    }

  }

  return b;
}
/******************************************************************************/

void r8ge_sl ( int n, double a_lu[], int pivot[], double x[], int job )

/******************************************************************************/
/*
  Purpose:

    R8GE_SL solves a R8GE system factored by R8GE_FA.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

    R8GE_SL is a simplified version of the LINPACK routine SGESL.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    09 April 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the order of the matrix.
    N must be positive.

    Input, double A_LU[N*N], the LU factors from R8GE_FA.

    Input, int PIVOT[N], the pivot vector from R8GE_FA.

    Input/output, double X[N], on input, the right hand side vector.
    On output, the solution vector.

    Input, int JOB, specifies the operation.
    0, solve A * x = b.
    nonzero, solve A' * x = b.
*/
{
  int i;
  int k;
  int l;
  double t;
/*
  Solve A * x = b.
*/
  if ( job == 0 )
  {
/*
  Solve PL * Y = B.
*/
    for ( k = 1; k <= n-1; k++ )
    {
      l = pivot[k-1];

      if ( l != k )
      {
        t      = x[l-1];
        x[l-1] = x[k-1];
        x[k-1] = t;
      }
      for ( i = k+1; i <= n; i++ )
      {
        x[i-1] = x[i-1] + a_lu[i-1+(k-1)*n] * x[k-1];
      }
    }
/*
  Solve U * X = Y.
*/
    for ( k = n; 1 <= k; k-- )
    {
      x[k-1] = x[k-1] / a_lu[k-1+(k-1)*n];
      for ( i = 1; i <= k-1; i++ )
      {
        x[i-1] = x[i-1] - a_lu[i-1+(k-1)*n] * x[k-1];
      }
    }
  }
/*
  Solve A' * X = B.
*/
  else
  {
/*
  Solve U' * Y = B.
*/
    for ( k = 1; k <= n; k++ )
    {
      t = 0.0;
      for ( i = 1; i <= k-1; i++ )
      {
        t = t + x[i-1] * a_lu[i-1+(k-1)*n];
      }
      x[k-1] = ( x[k-1] - t ) / a_lu[k-1+(k-1)*n];
    }
/*
  Solve ( PL )' * X = Y.
*/
    for ( k = n-1; 1 <= k; k-- )
    {
      t = 0.0;
      for ( i = k+1; i <= n; i++ )
      {
        t = t + x[i-1] * a_lu[i-1+(k-1)*n];
      }
      x[k-1] = x[k-1] + t;

      l = pivot[k-1];

      if ( l != k )
      {
        t      = x[l-1];
        x[l-1] = x[k-1];
        x[k-1] = t;
      }
    }
  }
  return;
}