eg_dmatrix.h

Go to the documentation of this file.

/* 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 
 * */
#ifndef __EG_DMATRIX_H__
#define __EG_DMATRIX_H__
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <limits.h>
#include <math.h>
#include <float.h>
#include "eg_macros.h"
#include "eg_mem.h"
#include "eg_lpnum.h"
#include "eg_numutil.h"

/* ========================================================================= */
/** @defgroup EGdMatrix Dense Matrices
 * Here we define a common interface for dense matrices (i.e. a structure), and
 * some common operations over dense matrices. The definition uses EGlpNum as
 * reference number type, this allow for template initializations.
 * 
 * @par History:
 * Revision 0.0.2
 *  - 2005-10-27
 *            - First implementation.
 * */
/** @{*/
/** @file
 * @brief This file provide the user interface and function definitions for
 * Dense Matrices.
 * */
/** @example eg_dmatrix.ex.c */
/* ========================================================================= */
/** @brief structure to hold a dense matrix, we choose a row representation
 * of the matrix, and we allow row and column permutations. All actual values 
 * in the matrix are stored in #EGdMatrix_t::matval, and the rows in
 * #EGdMatrix_t::matrow. */
typedef struct EGdMatrix_t
{
  size_t col_sz;    /**< @brief Number of columns in the matrix. */
  size_t row_sz;    /**< @brief Number of rows in the matrix */
  EGlpNum_t **matrow;
                    /**< @brief Array of size #EGdMatrix_t::row_sz containing 
                         all rows of the matrix */
  EGlpNum_t *matval;/**< @brief Values for all entries */
  int *col_ord;     /**< @brief Array of size at least #EGdMatrix_t::col_sz 
                         containing the order ammong all columns i.e. it is a 
                         permutation of {0,....,col_sz-1} which is how the 
                         matrix is treated internally */
  int *row_ord;     /**< @brief Array of size at least #EGdMatrix_t::row_sz 
                         containing the order ammong all rows, i.e. it is a 
                         permutation of {0,...,row_sz-1} which is how the 
                         matrix is treated internally */
}
EGdMatrix_t;

/* ========================================================================= */
/** @brief Initialize (as a dense matrix of dimension 0x0) an #EGdMatrix_t
 * structure.
 * @param dmatrix dense matrix structure pointer.
 * */
#define EGdMatrixInit(dmatrix) memset(dmatrix,0,sizeof(EGdMatrix_t))

/* ========================================================================= */
/** @brief Clear a dense matrix structure, i.e. free all internally allocated
 * data of the structure. Note that no further use of the structure can be made
 * unless it is re-initialized and set to a suitable size.
 * @param dmatrix dense matrix structure pointer.
 * */
#define EGdMatrixClear(dmatrix) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  EGlpNumFreeArray(__EGdm->matval);\
  EGfree(__EGdm->matrow);\
  int_EGlpNumFreeArray(__EGdm->col_ord);\
  int_EGlpNumFreeArray(__EGdm->row_ord);} while(0)

/* ========================================================================= */
/** @brief Set new dimensions for a dense matrix structure.
 * @param dmatrix dense matrix structure pointer.
 * @param new_nrows number of rows in the matrix.
 * @param new_ncols number of columns in the matrix.
 * @note Take care that the values stored in the matrix are not initialized to
 * any particular number. Also the ordering (for both column and row) is reset
 * to the standard ordering 0,....,n.
 * */
#define EGdMatrixSetDimension(dmatrix,new_nrows,new_ncols) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  register unsigned __EGdmi;\
  __EGdm->col_sz = (new_ncols);\
  __EGdm->row_sz = (new_nrows);\
  EGlpNumReallocArray(&(__EGdm->matval),__EGdm->col_sz * __EGdm->row_sz);\
  EGrealloc(__EGdm->matrow,__EGdm->row_sz * sizeof(EGlpNum_t*));\
  int_EGlpNumReallocArray(&(__EGdm->col_ord),__EGdm->col_sz);\
  int_EGlpNumReallocArray(&(__EGdm->row_ord),__EGdm->row_sz);\
  __EGdmi = __EGdm->col_sz;\
  while(__EGdmi--) __EGdm->col_ord[__EGdmi] = __EGdmi;\
  __EGdmi = __EGdm->row_sz;\
  while(__EGdmi--) \
    __EGdm->matrow[__EGdmi] = __EGdm->matval + (__EGdmi * __EGdm->col_sz);\
  __EGdmi = __EGdm->row_sz;\
  while(__EGdmi--) __EGdm->row_ord[__EGdmi] = __EGdmi;} while(0)

/* ========================================================================= */
/** @brief Display a given #EGdMatrix_t structure contents.
 * @param dmatrix dense matrix structure pointer.
 * @param nat_order if set to one, display the matrix using the natural 
 * internal order, i.e. we discard the order of columns and rows as defined in
 * #EGdMatrix_t::col_ord and #EGdMatrix_t::row_ord. Otherwise, use such orders.
 * @param out_file pointer to a FILE structure where we want the output to be
 * printed.
 * */
#define EGdMatrixDisplay(dmatrix,nat_order,out_file) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  char* __EGdmstr = 0;\
  size_t __EGdmi, __EGdmj;\
  fprintf(out_file,"Matrix %p\nDimensions: %zd rows, %zd columns\n", (void*)__EGdm, __EGdm->row_sz, __EGdm->col_sz);\
  if(nat_order){\
    for(__EGdmi = 0 ; __EGdmi < __EGdm->row_sz ; __EGdmi++){\
      for(__EGdmj = 0 ; __EGdmj < __EGdm->col_sz ; __EGdmj++){\
        __EGdmstr = EGlpNumGetStr(__EGdm->matrow[__EGdmi][__EGdmj]);\
        fprintf(out_file,"%10s ", __EGdmstr);\
        EGfree(__EGdmstr);\
      }\
      fprintf(out_file,"\n");}\
  } else {\
    for(__EGdmi = 0 ; __EGdmi < __EGdm->row_sz ; __EGdmi++){\
      for(__EGdmj = 0 ; __EGdmj < __EGdm->col_sz ; __EGdmj++){\
        __EGdmstr = EGlpNumGetStr(__EGdm->matrow[__EGdm->row_ord[__EGdmi]][__EGdm->col_ord[__EGdmj]]);\
        fprintf(out_file,"%10s ", __EGdmstr);\
        EGfree(__EGdmstr);\
      }\
      fprintf(out_file,"\n");}\
  }} while(0)

/* ========================================================================= */
/** @brief Given a number 'num' and a two rows 'ori', 'dest', set rows 
 * 'dest' to 'dest' + 'ori' * 'num'. Note that the number MUST_NOT be stored 
 * in row 'dest', and note that rows 'ori' and 'dest' should be different.
 * This is needed because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param 'ori' index of the row whose multiple will be added to the 'dest'
 * row.
 * @param 'dest' row to be replaced by 'dest' + 'ori' * 'num'.
 * @param 'num' constant to be multiply to the 'ori' and be added to the 
 * 'dest' row.
 * @note The index of the row are taken as internal index, i.e. if we give row
 * 'k' we will use the row stored in #EGdMatrix_t::matrow[k], wich does not
 * mean that we will access the k-th row in the matrix (wich would need to use
 * as index the value #EGdMatrix_t::row_ord[k] instead). Note that we don't
 * test wether the given multiple is zero or not. we always perform the
 * operation.
 * */
#define EGdMatrixAddRowMultiple(dmatrix,dest,ori,num) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdest = (dest);\
  const size_t __EGori = (ori);\
  size_t __EGdmj = __EGdm->col_sz;\
  while(__EGdmj--) \
    EGlpNumAddInnProdTo(__EGdm->matrow[__EGdest][__EGdmj],\
                        __EGdm->matrow[__EGori][__EGdmj],num);\
  } while(0)

/* ========================================================================= */
/** @brief Given a number 'num' and a two rows 'ori', 'dest', set rows 
 * 'dest' to 'dest' - 'ori' * 'num'. Note that the number MUST_NOT be stored 
 * in row 'dest', and note that rows 'ori' and 'dest' should be different.
 * This is needed because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param 'ori' index of the row whose multiple will be added to the 'dest'
 * row.
 * @param 'dest' row to be replaced by 'dest' - 'ori' * 'num'.
 * @param 'num' constant to be multiply to the 'ori' and be added to the 
 * 'dest' row.
 * @note The index of the row are taken as internal index, i.e. if we give row
 * 'k' we will use the row stored in #EGdMatrix_t::matrow[k], wich does not
 * mean that we will access the k-th row in the matrix (wich would need to use
 * as index the value #EGdMatrix_t::row_ord[k] instead). Note that we don't
 * test wether the given multiple is zero or not. we always perform the
 * operation.
 * */
#define EGdMatrixSubRowMultiple(dmatrix,dest,ori,num) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdest = (dest);\
  const size_t __EGori = (ori);\
  size_t __EGdmj = __EGdm->col_sz;\
  while(__EGdmj--) \
    EGlpNumSubInnProdTo(__EGdm->matrow[__EGdest][__EGdmj],\
                        __EGdm->matrow[__EGori][__EGdmj],num);\
  } while(0)
/* ========================================================================= */
/** @brief Given a number and a row, multiply the complete row by the given
 * number. Note that the number MUST_NOT be stored in the row being multiplied,
 * this is because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param row_ind index of the row being multiplied, note that we will multiply
 * the row stored in #EGdMatrix_t::matrow[row_ind], wich is different to say
 * that we multiply the row in the row_ind-th position in the row ordering (to
 * do that, then row_ind should be #EGdMatrix_t::row_ord[k]).
 * @param multiple constant to be multiply to the row.
 * */
#define EGdMatrixMultiplyRow(dmatrix,row_ind,multiple) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdmi = (row_ind);\
  size_t __EGdmj = __EGdm->col_sz;\
  while(__EGdmj--) EGlpNumMultTo(__EGdm->matrow[__EGdmi][__EGdmj],multiple);\
  } while(0)

/* ========================================================================= */
/** @brief Given a number 'num' and a two rows 'ori', 'dest', set columns 
 * 'dest' to 'dest' + 'ori' * 'num'. Note that the number MUST_NOT be stored 
 * in column 'dest', and note that columns 'ori' and 'dest' should be 
 * different. This is needed because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param 'ori' index of the column whose multiple will be added to the 'dest'
 * column.
 * @param 'dest' column to be replaced by 'dest' + 'ori' * 'num'.
 * @param 'num' constant to be multiply to the 'ori' and be added to the 
 * 'dest' column.
 * @note The index of the column are taken as internal index, i.e. if we give 
 * column 'k' we will use the column stored in #EGdMatrix_t::matrow[*][k], 
 * wich does not mean that we will access the k-th column in the matrix (wich
 * would need to use as index the value #EGdMatrix_t::row_ord[k] instead). 
 * Note that we don't test wether the given multiple is zero or not. we 
 * always perform the operation.
 * */
#define EGdMatrixAddColMultiple(dmatrix,dest,ori,num) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdest = (dest);\
  const size_t __EGori = (ori);\
  size_t __EGdmj = __EGdm->row_sz;\
  while(__EGdmj--) \
    EGlpNumAddInnProdTo(__EGdm->matrow[__EGdmj][__EGdest],\
                        __EGdm->matrow[__EGdmj][__EGori],num);\
  } while(0)

/* ========================================================================= */
/** @brief Given a number 'num' and a two rows 'ori', 'dest', set columns 
 * 'dest' to 'dest' - 'ori' * 'num'. Note that the number MUST_NOT be stored 
 * in column 'dest', and note that columns 'ori' and 'dest' should be 
 * different. This is needed because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param 'ori' index of the column whose multiple will be added to the 'dest'
 * column.
 * @param 'dest' column to be replaced by 'dest' - 'ori' * 'num'.
 * @param 'num' constant to be multiply to the 'ori' and be added to the 
 * 'dest' column.
 * @note The index of the column are taken as internal index, i.e. if we give 
 * column 'k' we will use the column stored in #EGdMatrix_t::matrow[*][k], 
 * wich does not mean that we will access the k-th column in the matrix (wich 
 * would need to use as index the value #EGdMatrix_t::col_ord[k] instead). 
 * Note that we don't test wether the given multiple is zero or not. we 
 * always perform the operation.
 * */
#define EGdMatrixSubColMultiple(dmatrix,dest,ori,num) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdest = (dest);\
  const size_t __EGori = (ori);\
  size_t __EGdmj = __EGdm->row_sz;\
  while(__EGdmj--) \
    EGlpNumSubInnProdTo(__EGdm->matrow[__EGdmj][__EGdest],\
                        __EGdm->matrow[__EGdmj][__EGori],num);\
  } while(0)
/* ========================================================================= */
/** @brief Given a number and a column, multiply the complete column by the 
 * given number. Note that the number MUST_NOT be stored in the column being 
 * multiplied, this is because of the way GNU_MP interface works.
 * @param dmatrix dense matrix structure pointer.
 * @param col_ind index of the column being multiplied, note that we will 
 * multiply the column stored in #EGdMatrix_t::matrow[*][col_ind], wich is 
 * different to say that we multiply the column in the col_ind-th position in
 * the column ordering (to do that, then col_ind should be 
 * #EGdMatrix_t::col_ord[k]).
 * @param multiple constant to be multiply to the column.
 * */
#define EGdMatrixMultiplyCol(dmatrix,col_ind,multiple) do{\
  EGdMatrix_t*const __EGdm = (dmatrix);\
  const size_t __EGdmi = (col_ind);\
  size_t __EGdmj = __EGdm->row_sz;\
  while(__EGdmj--) EGlpNumMultTo(__EGdm->matrow[__EGdmj][__EGdmi],multiple);\
  } while(0)


/* ========================================================================= */
/** @brief This function performs gaussian elimination to the given matrix,
 * depending on the given options it may do row/columns permutations allong the
 * way to improve numerical stabillity.
 * @param dmatrix dense matrix structure pointer.
 * @param do_col_perm if set to one, the try columns permutation to improve
 * numericall stabillity, otherwise, not do column permutations at all.
 * @param do_row_perm if set to one, try row permutations to improve numericall
 * stabillity, otherwise, not do row permutations at all.
 * @param status pointer to where return an status, if the procedure finish all
 * the way (i.e. the matrix is full rank), then we return #EG_ALGSTAT_SUCCESS,
 * if the matrix is found to be partial rank, the status is
 * #EG_ALGSTAT_PARTIAL, otherwise, we return #EG_ALGSTAT_NUMERROR, wich means
 * that we stoped because a zero pivot was found (after checking for allowed
 * row/collumns permmutations).
 * @param rank where to return the (proven) rank of the matrix. This number is
 * accurate if the status is #EG_ALGSTAT_SUCCESS, or #EG_ALGSTAT_PARTIAL, but
 * is just a lower bound if the status is #EG_ALGSTAT_NUMERROR
 * @param zero_tol What is the threshold for a value to be considered zero.
 * @return if no error happen, we return zero, otherwise a non-zero valued is
 * returned. Note that the algorithm status is independent of the return value,
 * non zero values araise only if an error happen during execution, wich is
 * different to say that the algorithm didn't finish correctly. */
int EGdMatrixGaussianElimination (EGdMatrix_t * const dmatrix,
                                  const unsigned do_col_perm,
                                  const unsigned do_row_perm,
                                  unsigned *const rank,
                                  EGlpNum_t zero_tol,
                                  int *const status);

/* ========================================================================= */
/** @}*/
#endif

---