/* 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; }
Interested in particle astrophysics and plasma astrophysics. This blog is my research/private notebook.
Wednesday, October 31, 2018
Strengthen the matrix in Gun Scientific Library (GSL)
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment