GSLUtility.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include <string>
#include <vector>
#include <numeric>
#include <iostream>
#include <sstream>
 
#include <gsl/gsl_errno.h>
 
#include "GSLUtility.h"
#include "IException.h"
 
using namespace std;
 
namespace Isis {
  namespace GSL {
 
    GSLUtility *GSLUtility::_instance = 0;
 
    GSLUtility::GSLUtility() {
      gsl_set_error_handler(handler);
    }
 
 
    GSLUtility *GSLUtility::getInstance() {
      if(_instance == 0) {
        _instance = new GSLUtility();
      }
      return (_instance);
    }
 
 
    gsl_vector *GSLUtility::vector(size_t n, bool zero) const {
      if(zero) {
        return (gsl_vector_calloc(n));
      }
      else {
        return (gsl_vector_alloc(n));
      }
    }
 
    gsl_matrix *GSLUtility::matrix(size_t n1, size_t n2, bool zero) const {
      if(zero) {
        return (gsl_matrix_calloc(n1, n2));
      }
      else {
        return (gsl_matrix_alloc(n1, n2));
      }
    }
 
    gsl_matrix *GSLUtility::identity(size_t n1, size_t n2) const {
      gsl_matrix *i = gsl_matrix_alloc(n1, n2);
      gsl_matrix_set_identity(i);
      return (i);
    }
 
    void GSLUtility::setIdentity(gsl_matrix *m) const {
      gsl_matrix_set_identity(m);
      return;
    }
 
    void GSLUtility::free(gsl_vector *v) const {
      gsl_vector_free(v);
      return;
    }
 
    void GSLUtility::free(gsl_matrix *m) const {
      gsl_matrix_free(m);
      return;
    }
 
    GSLUtility::GSLVector GSLUtility::gslToGSL(const gsl_vector *v) const {
      size_t n = size(v);;
      GSLVector Nv(n);
      for(size_t i = 0 ; i < n ; i++) {
        Nv[i] = gsl_vector_get(v, i);
      }
      return (Nv);
    }
 
    GSLUtility::GSLMatrix GSLUtility::gslToGSL(const gsl_matrix *m) const {
      size_t nrows = Rows(m);
      size_t ncols = Columns(m);
      GSLMatrix Nm(nrows, ncols);
      for(size_t i = 0 ; i < nrows ; i++) {
        for(size_t j = 0 ; j < ncols ; j++) {
          Nm[i][j] = gsl_matrix_get(m, i, j);
        }
      }
      return (Nm);
    }
 
 
    gsl_vector *GSLUtility::GSLTogsl(const GSLUtility::GSLVector &v,
                                     gsl_vector *gv) const {
      if(gv == 0) {
        gv = gsl_vector_alloc(v.dim());
      }
      else if(size(gv) != (size_t) v.dim()) {
        ostringstream mess;
        mess << "Size of NL vector (" << v.dim() << ") not same as GSL vector ("
             << gv->size << ")";
        throw IException(IException::Programmer,
                         mess.str().c_str(),
                         _FILEINFO_);
      }
 
      for(int i = 0 ; i < v.dim() ; i++) {
        gsl_vector_set(gv, i, v[i]);
      }
      return (gv);
    }
 
    gsl_matrix *GSLUtility::GSLTogsl(const GSLUtility::GSLMatrix &m,
                                     gsl_matrix *gm) const {
      if(gm == 0) {
        gm = gsl_matrix_alloc(m.dim1(), m.dim2());
      }
      else if((Rows(gm) != (size_t) m.dim1()) &&
              (Columns(gm) != (size_t) m.dim2())) {
        ostringstream mess;
        mess << "Size of NL matrix (" << m.dim1() << "," << m.dim2()
             << ") not same as GSL matrix (" << Rows(gm) << "," << Columns(gm)
             << ")";
        throw IException(IException::Programmer,
                         mess.str().c_str(),
                         _FILEINFO_);
      }
 
      for(int i = 0 ; i < m.dim1() ; i++) {
        for(int j = 0 ; j < m.dim2() ; j++) {
          gsl_matrix_set(gm, i, j, m[i][j]);
        }
      }
      return (gm);
    }
 
    size_t GSLUtility::Rows(const gsl_matrix *m) const {
      return (m->size1);
    }
 
    size_t GSLUtility::Columns(const gsl_matrix *m) const {
      return (m->size2);
    }
 
    size_t GSLUtility::Columns(const GSLMatrix &m) const {
      return (m.dim1());
    }
 
    size_t GSLUtility::Rows(const GSLMatrix &m) const {
      return (m.dim2());
    }
 
    size_t GSLUtility::size(const gsl_vector *v) const {
      return (v->size);
    }
 
    size_t GSLUtility::size(const gsl_matrix *m) const {
      return (Rows(m) * Columns(m));
    }
 
    void GSLUtility::check(int gsl_status, const char *src, int line) const {
      if(gsl_status != GSL_SUCCESS) {
        string msg = "GSL error occured: " + string(gsl_strerror(gsl_status));
        throw IException(IException::Programmer, msg.c_str(), src, line);
      }
      return;
    }
 
    void GSLUtility::handler(const char *reason, const char *file, int line,
                             int gsl_errno) {
      ostringstream mess;
      mess << "GSLError (" << gsl_errno << ") -> " << reason;
      throw IException(IException::Programmer, mess.str().c_str(),
                       file, line);
    }
 
  }
} // namespace ISIS::GSL