Labels

Wednesday, October 31, 2018

Strengthen the matrix in Gun Scientific Library (GSL)

/* MYMATRIX.h */

 #include <stdio.h>
 #include <string.h>
 #include <vector>
 #include <assert.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
 /* arbitrary maximum size */
 #define MAXSIZE 5000*10000 

 void print_matrix( const gsl_matrix *,const char *);
 void fprint_matrix( const char *filename, const gsl_matrix *m, const char *comment ="");
 gsl_matrix * fread_matrix( const char *filen, int *rows, int *clos, int transpose);
 gsl_matrix * vector2matrix(const std::vector<double> x,const std::vector<double>y);



/* MYMATRIX.cc */
#include "MYMATRIX.h"

void print_matrix(const gsl_matrix *m, const char *name) 
{
  int r,c;
  int rows = m->size1;
  int cols = m->size2;
  printf ("Matrix %s (%d x %d):\n", name, rows, cols );
  for( r = 0; r < rows; r++ ) { 
    for( c = 0; c < cols; c++ ) { 
      printf( "%8e", gsl_matrix_get (m, r, c) );
    }
    printf ("\n");
  }
}


/*
 * fprint_matrix - print matrix to a filename (free format, 8 decimals)
 */
void fprint_matrix( const char *filename, const gsl_matrix *m, 
 const char *comment ) {
  int r,c;
  FILE * fh;
  int rows = m->size1;
  int cols = m->size2;
  fh = fopen( filename, "w" );

  if( fh == NULL ) { 
    fprintf( stderr, "Cannot open \"%s\"!\n", filename );
    exit(-1);
  }

  if (strlen(comment) !=0){
    fprintf ( fh, "#%s %d %d\n",comment, rows, cols );} // For reading to matrix next time
  else{
    fprintf ( fh, "# %d %d\n", rows, cols );
  }
  for( r = 0; r < rows; r++ ) { 
    for( c = 0; c < cols; c++ ) { 
      fprintf( fh, "%e ", gsl_matrix_get (m, r, c) );
    }
    fprintf ( fh, "\n");
  }
  fclose(fh);

}


gsl_matrix * fread_matrix(const char *filen, int *rows, int *cols, int transpose ) {
  int r,c,rc;
  int n;
  double temp; /* temporary read variable */
  gsl_matrix *m;

  /* open the file */
  FILE *fi;
  fi = fopen(filen, "r");
  if( fi == NULL ) { 
    printf( "Cannot open \"%s\"\n", filen );
    return(NULL);
  }

  /* read numbers of rows and columns */
  /* required the first line in input file is :"comment rows cols " format!!! */
  char *comment;
  fscanf( fi, "#%s %d %d",comment, rows, cols ); 

  rc = (*rows) * (*cols);

  /* check for valid; max size is pretty arbitrary */
  if( *rows <= 0 || *cols <= 0 || rc > MAXSIZE ) { 
//    printf("Read invalid row and/or column from first line from %c %s %c (max 
//size=%d).\n",'\"', filen,'\"', MAXSIZE );
   printf("Read invalid row and/or column from first line from %c %s %c (max size=%d).\n",'\"', filen,'\"', MAXSIZE );
    return(NULL);
  }

  /* request memory for this matrix */
  m = gsl_matrix_alloc ( *rows, *cols );
  if( m == NULL ) { 
    printf( "Cannot allocate a matrix with %d elements (%d rows and %d cols).", 
rc, *rows, *cols );
    return(NULL);
  }

  /* read in the matrix */
  /* what errors will break this and should they be trapped? */
  n = 0;
  r = c = 0;
  while ( fscanf(fi, "%lf", &temp) != EOF ) {
    n++;
    if( transpose ) {
      gsl_matrix_set( m, c, r, temp ); /* TRANSPOSE */
    } else {
      gsl_matrix_set( m, r, c, temp );
    }
    c++;
    if ( c >= *cols ) { 
      c = 0;
      r++;
      if( r >= *rows ) 
        break; /* should this be an error condition? */
    }
  }

  /* check that we read the right amount of data */
  if( n != rc ) { 
    printf( "Error: Expected %d elements for %d x %d matrix, but read %d!\n", 
rc, *rows, *cols, n );
    return NULL;
  }

  /* close the file */
  fclose(fi);

  return( m );
}

gsl_matrix * vector2matrix(std::vector<double> x,std::vector<double>y)
{
    assert(("two vector must be equal length",x.size()==y.size()));
    gsl_matrix * M = gsl_matrix_alloc(2,x.size());
    for(std::size_t i=0;i<x.size();i++)
    {
        gsl_matrix_set(M,0,x[i]);
        gsl_matrix_set(M,1,y[i]);
    }
    return M;
}

No comments:

Post a Comment