gsl_vector_matrix_ops.h

Go to the documentation of this file.

//

//--------------------------------------------------------------------------
//--------------------------------------------------------------------------
//  
//  libMoM Library - a library for stochastic simulations in engineering.
// 
//  Copyright (C) 2010 Rochan R. Upadhyay
// 
//  This library is free software; you can redistribute it and/or
//  modify it under the terms of the Version 2.1 GNU Lesser General
//  Public License as published by the Free Software Foundation.
// 
//  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 Street, Fifth Floor, 
//  Boston, MA  02110-1301  USA
// 
//--------------------------------------------------------------------------

// Some functions are taken from the gslwrap project whose Copyright notice is below.

//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C)  ULP-IPB Strasbourg

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.

//  This program 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.


#ifndef _GSL_VECTOR_MATRIX_OPS_H
#define _GSL_VECTOR_MATRIX_OPS_H

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_eigen.h>
#include "gsl_vector_double.h"
#include "gsl_vector_integer.h"
#include "gsl_matrix_double.h"
#include "gsl_matrix_integer.h"
#include "gsl_permutation.h"
#include <glpk.h>

//
//

inline 
vector operator*(const matrix& m, const vector& v)
{
        vector y(m.get_rows());
        gsl_blas_dgemv(CblasNoTrans, 1.0, m.gslobj(), v.gslobj(), 0.0, y.gslobj());
        return y;
}

//
//

inline
vector LU_solve (const matrix &m, const vector &v)
{
   permutation p;
   vector y(m.get_cols());
   matrix m_lud(m.get_rows(), m.get_cols());
   m_lud = m.LU_decomp(&p);
   gsl_linalg_LU_solve(m_lud.gslobj(), p.gslperm_obj(), v.gslobj(), y.gslobj()); 
   return y;
}

//
//

inline
vector singular_values (const matrix &m)
{
   vector s(m.get_cols());
   vector work(m.get_cols());
   matrix v(m.get_cols(),m.get_cols());
   matrix a(m.get_rows(),m.get_cols());
   a=m;
   gsl_linalg_SV_decomp(a.gslobj(), v.gslobj(), s.gslobj(), work.gslobj());
   return s;
}

//
//

inline
double condition_number (const matrix &m)
{
   vector s(m.get_cols());
   double smax,smin,cond_number;
   s = singular_values(m);
   smax = s.max();
   smin = s.min();
   if (smin > 0 || smin < 0) {
   cond_number = smax/smin;
   return cond_number;}
   else return 1.e300; 
}

//
//

inline
matrix_int addrow(const size_t &irow, const size_t &itot, const vector_int &v)
{
   matrix_int y(itot, v.size());  
   for(size_t j = 0; j < v.size(); j++ ) {
      y.set_element( irow, j, v.get(j) );}
   return y;
}

//
//

inline
matrix addrow(const size_t &irow, const size_t &itot, const vector &v)
{
   matrix y(itot, v.size());  
   for(size_t j = 0; j < v.size(); j++ ) {
      y.set_element( irow, j, v.get(j) );}
   return y;
}

//
//

inline
matrix extract_cols(const matrix &m, const vector_int &v)
{
  size_t num_rows = m.get_rows();
  size_t num_cols = v.size();
  matrix y(num_rows,num_cols);
  for(size_t i = 0; i < num_rows; i++){
    for(size_t j = 0; j < num_cols; j++){
      y.set_element(i,j,m.get_element(i,v.get(j)));}}
  return y;
}  

//
//

inline
matrix tridiag_eigen_solve (const vector &v1, const vector &v2)
{
  size_t msize = v1.size();
  matrix m( msize, msize );
  vector eval( msize );
  matrix evec( msize, msize );
  matrix eigvalvec( msize+1, msize);
  vector colvec;
  m.set_tridiag(v1, v2, v2);
  gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc ( msize );
  gsl_eigen_symmv ( m.gslobj(), eval.gslobj(), evec.gslobj(), w );   
  gsl_eigen_symmv_free ( w );
  gsl_eigen_symmv_sort ( eval.gslobj(), evec.gslobj(), GSL_EIGEN_SORT_ABS_ASC );
  eigvalvec = addrow(0, msize+1, eval);
  for(size_t i=0; i < msize; i++){
    colvec = evec.row(i);
    eigvalvec += addrow(i+1, msize+1, colvec);}
  return ( eigvalvec );
}

//
//

inline
vector subvector_elements (const vector &v, const vector_int &index)
{
  vector y(index.size()); 
  for (size_t i = 0; i < index.size(); i++){
    y.set(i,v.get(index.get(i)));}
  return y;  
} 

//
//

inline
vector simplex_method (const matrix &A, const vector &b, const vector &c)
{
  size_t m = A.get_rows();
  size_t n = A.get_cols();
  vector xsimp(n);
  int* ia = NULL;
  int *ja = NULL;
  double *ar = NULL;
  ia = new int[1+m*n];
  ja = new int[1+m*n];
  ar = new double[1+m*n];
  for (int ir=0; ir<m*n+1; ir++) {
    ia[ir] = 0;
    ja[ir] = 0;
    ar[ir] = 0.0;}
  glp_prob *lp;
  glp_smcp parm;
  lp = glp_create_prob();
  glp_init_smcp(&parm);
  parm.msg_lev = GLP_MSG_ERR;
  parm.it_lim = 50000;
  parm.meth = GLP_DUALP;
  parm.tol_bnd = 1.0e-7;
  parm.tol_dj = 1.0e-7;
  glp_set_obj_dir(lp, GLP_MIN);
  glp_add_rows(lp, m);
  glp_add_cols(lp, n);
  for (size_t ir = 1; ir<m+1; ir++){
    glp_set_row_bnds(lp, ir, GLP_FX, b.get(ir-1), b.get(ir-1));}
  for (size_t ic = 1; ic<n+1; ic++){
     glp_set_col_bnds(lp, ic, GLP_LO, 0.0, 0.0);
     glp_set_obj_coef(lp, ic, c.get(ic-1));}
  size_t ircnt = 1;
  for (size_t ira = 1; ira<m+1; ira++){
    for (size_t jca = 1; jca<n+1; jca++){
       ia[ircnt]=ira;
       ja[ircnt]=jca;
       ar[ircnt]=A.get_element(ira-1,jca-1);
       ircnt++;}}
  glp_load_matrix(lp, m*n, ia, ja, ar);
  glp_simplex(lp, &parm);
  for (size_t ic = 1; ic<n+1; ic++){
    xsimp.set(ic-1,glp_get_col_prim(lp, ic));}
  glp_delete_prob(lp);
  delete [] ia;
  ia = NULL;
  delete [] ja;
  ja = NULL;
  delete [] ar;
  ar = NULL;
  return xsimp;
}

//
//

inline
matrix_int glex_dtuples( size_t dim, size_t total)
{
  matrix_int g(total,dim);
  vector_int v(dim);
  for(size_t j = 0; j < dim; j++){
    v.set(j,0);}
  for(size_t i = 0; i < total; i++){
    for(size_t j = 0; j < dim; j++){
      g.set_element(i,j,v.get(j));}
    if (dim==1) v.set(0,i+1);  
    if (dim>1) v=v.glex();}
  return g;
}

#endif // GSL_VECTOR_MATRIX_OPS_H