eg_dmatrix.ex.c

/* EGlib "Efficient General Library" provides some basic structures and
 * algorithms commons in many optimization algorithms.
 *
 * Copyright (C) 2005 Daniel Espinoza and Marcos Goycoolea.
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation; either version 2.1 of the License, or (at your
 * option) any later version.
 *
 * This library is distributed in the hope that it will be useful, but 
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 
 * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public 
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA 
 * */
/** @file 
 * @ingroup EGdMatrix
 * @brief This file shows some examples of the use of #EGdMatrix_t structure 
 * and functions. as a test-set we use the Hilber Matrix wich is a matrix
 * with coefficients \f$H^n_{i,j} = \frac1{i+j-1}\f$. It is well known that the
 * Hilbert matrix is extremmaly ill conditioned, and is very hard to invert,
 * here are some values for it's determinant for different dimmensions \f$n\f$.
 * - \f$\mathrm{det}(H_1) = 1 \f$. 
 * - \f$\mathrm{det}(H_2) = \frac1{12} \f$. 
 * - \f$\mathrm{det}(H_3) = \frac1{2160} \f$. 
 * - \f$\mathrm{det}(H_4) = \frac1{604800} \f$. 
 * - \f$\mathrm{det}(H_5) = \frac1{266716800000} \f$. 
 * Since a lot is known about the Hilber Matrix, some accuracy tests are
 * possible. For more details see <A HREF=http://mathworld.wolfram.com/HilbertMatrix.html TARGET=_top>MathWorld</A>
 * This test program takes as input the number of columns and rows for the
 * hilbert matrix, and so some operations on it and display it to the standard
 * output.
 * */
/** @addtogroup EGdMatrix */
/** @{ */
/* ========================================================================= */
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include "eg_dmatrix.h"
#include "eg_dbasis_red.h"

/* ========================================================================= */
/** @brief size of the table of eigenvalues of the hilbert matrix */
#define HILBERT_TABLE_SIZE 9U

/* ========================================================================= */
/** @brief table containing 1/eigenvalues of the hilbert matrix. */
static unsigned dmatrix_hilbert_eigenvalues[HILBERT_TABLE_SIZE] = {
  1,
  12,
  180,
  2800,
  44100,
  698544,
  11099088,
  176679360,
  2815827300
};

/* ========================================================================= */
/** @brief display the usage message for this program */
void dmatrix_usage (char *program)
{
  fprintf (stdout, "Usage: %s [options]\n", program);
  fprintf (stdout, "Options:\n");
  fprintf (stdout, "     -n n   number of rows in the hilbert matrix.\n");
  fprintf (stdout, "     -m n   number of columns in the hilbert matrixt.\n");

}

/* ========================================================================= */
/** @brief parse the input argumenbts for the program */
int dmatrix_parseargs (int argc,
                       char **argv,
                       size_t * n,
                       size_t * m)
{
  int c;
  while ((c = getopt (argc, argv, "n:m:")) != EOF)
  {
    switch (c)
    {
    case 'n':
      *n = atoi (optarg);
      break;
    case 'm':
      *m = atoi (optarg);
      break;
    default:
      dmatrix_usage (argv[0]);
      return 1;
    }
  }
  /* reporting the options */
  if (!(*m) || !(*n))
  {
    dmatrix_usage (argv[0]);
    return 1;
  }
  fprintf (stdout, "Parsed Options:\n");
  fprintf (stdout, "n : %u\n", *n);
  fprintf (stdout, "m : %u\n", *m);
  return 0;
}

/* ========================================================================= */
/** @brief main function */
int main (int argc,
          char **argv)
{

  int rval = 0,
    status;
  size_t n = 1,
    m = 1;
  register unsigned i = 0,
    j = 0;
  unsigned rank;
  EGdMatrix_t dmatrix;
  EGdBsRed_t bsred;
  EGlpNum_t error;
  /* start */
  EGlpNumStart ();
  EGlpNumInitVar (error);
  EGdMatrixInit (&dmatrix);
  EGdBsRedInit (&bsred);
  /* parse args */
  rval = dmatrix_parseargs (argc, argv, &n, &m);
  CHECKRVALG (rval, CLEANUP);
  /* set the matrix size */
  EGdMatrixSetDimension (&dmatrix, n, m);
  /* set the coefficients of the hilbert matrix */
  for (i = n; i--;)
    for (j = m; j--;)
    {
      EGlpNumOne (dmatrix.matrow[i][j]);
      EGlpNumDivUiTo (dmatrix.matrow[i][j], i + j + 1);
    }
  /* Display the hilbert matrix */
  fprintf (stdout, "Hilber ");
  EGdMatrixDisplay (&dmatrix, 1, stdout);
  fflush (stdout);
  /* now we call gaussian elimination */
  rval =
    EGdMatrixGaussianElimination (&dmatrix, 0, 0, &rank, epsLpNum, &status);
  CHECKRVALG (rval, CLEANUP);
  fprintf (stdout, "Gaussian Elimination ended with status %d, rank %u and ",
           status, rank);
  EGdMatrixDisplay (&dmatrix, 0, stdout);
  fflush (stdout);
  /* now we compute the worst error in the diagonal up to
   * EGmin(HILBERT_TABLE_SIZE,rank) */
  fprintf (stdout, "Relative errors in eigenvalues:\n");
  for (i = EGmin (HILBERT_TABLE_SIZE, rank); i--;)
  {
    EGlpNumCopy (error, dmatrix.matrow[dmatrix.row_ord[i]][dmatrix.col_ord[i]]);
    fprintf (stdout, "Lambda_%u = %lg", i + 1, EGlpNumToLf (error));
    EGlpNumMultUiTo (error, dmatrix_hilbert_eigenvalues[i]);
    EGlpNumSubTo (error, oneLpNum);
    fprintf (stdout, " error = %lg\n", EGlpNumToLf (error));
  }
  /* set the coefficients of the hilbert matrix */
  for (i = n; i--;)
    for (j = m; j--;)
    {
      EGlpNumOne (dmatrix.matrow[i][j]);
      EGlpNumDivUiTo (dmatrix.matrow[i][j], i + j + 1);
    }
  fprintf (stdout, "Using Basis Reduction on ");
  EGdMatrixDisplay (&dmatrix, 1, stdout);
  /* now we do basis reduction over the rows of the matrix */
  /* ending */
  EGdBsRedSetLength (&bsred, m);
  for (i = 0; i < rank; i++)
    EGdBsRedAddElement (&bsred, dmatrix.matrow[dmatrix.row_ord[i]]);
  /* now call the Basis Reduction Algorithm and print the reduced basis */
  rval = EGdBsRed (&bsred, &rank, epsLpNum, &status);
  CHECKRVALG (rval, CLEANUP);
  fprintf (stdout, "Basis Reduction ended with status %d, rank %u and ",
           status, rank);
  EGdMatrixDisplay (&dmatrix, 1, stdout);
  fflush (stdout);
CLEANUP:
  EGdBsRedProfile (stderr);
  EGlpNumClearVar (error);
  EGdBsRedClear (&bsred);
  EGdMatrixClear (&dmatrix);
  EGlpNumExit ();
  return rval;
}

/* ========================================================================= */
/** @} */

---