Re: [Help-gsl] Reading data from a file to a vector

[Alan Mead]

Steven Vancoillie wrote:
> It's not because you use GSL that you're forbidden to use just C :)
>
> With x and y are pointers to a gsl_vector, one can use:
>
>  FILE * input = fopen("filename", "r");
>  int i = 0; while (fscanf (input, "%lf%lf", x->data+i, y->data+i) != EOF) ++i;
>
> I assume you know the length of the vectors, in which case you could
> use a for loop of course.
>   

Being new to C and GSL, I've written the attached; however, this reading 
GSL structures doesn't seem trivial to me.

In your example, one cannot allocate the vectors x and y unless one 
knows how many elements but one doesn't know how many elements until you 
open the file.  In my attached code, I assume that the rows and columns 
are the first line of the file (if I didn't, I'd have to seek through 
the entire file to find its size, right?).  Unless you always allocate 
some default huge structure (that might be less an issue for vectors).  
How do non-noob C programmers handle this chicken-and-egg problem?

Also, since I only find the size of the matrix inside my read_matrix(), 
I must allocate memory there.  I'm given to understand that good form to 
allocate at the same level you will free, and that my code is "bad" but 
I don't see any alternative. 

Related to this, how would I return two pointers to an x & y vector?

-Alan

--
Alan D. Mead, Ph.D.
Assistant Professor of Industrial and Organization Psychology
Scientific Advisor to the Center for Research and Service
Illinois Institute of Technology
3101 South Dearborn, 2nd floor
Chicago IL 60616
+312.567.5933 (Campus)
+815.588.3846 (Home Office)
+312.567.3493 (Fax)

http://www.iit.edu/~mead
http://www.center.iit.edu
http://www.alanmead.org

Be not afraid of life. 
Believe that life is worth living, 
and your belief will help create the fact.
-- Henry James

#include "read_matrix.h"

void print_matrix( const gsl_matrix *m, char *name, int rows, int cols ) {
  int r,c;

  printf ("Matrix %s (%d x %d):\n", name, rows, cols );
  for( r = 0; r < rows; r++ ) { 
    for( c = 0; c < cols; c++ ) { 
      printf( "%8.3f", gsl_matrix_get (m, r, c) );
    }
    printf ("\n");
  }

}

/*
 * fprint_matrix - print matrix to a filename (free format, 8 decimals)
 */
void fprint_matrix( char *filename, const gsl_matrix *m, char *name, int rows, 
int cols ) {
  int r,c;
  FILE * fh;

  fh = fopen( filename, "w" );

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

  fprintf ( fh, "%d %d\n", rows, cols );
  for( r = 0; r < rows; r++ ) { 
    for( c = 0; c < cols; c++ ) { 
      fprintf( fh, "%.8f ", gsl_matrix_get (m, r, c) );
    }
    fprintf ( fh, "\n");
  }
  fclose(fh);

}


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

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

  /* read numbers of rows and columns */
  fscanf( fi, "%d %d", 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 \"%s\" (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 );
}




#include <stdio.h>
#include <string.h>
#include <gsl/gsl_linalg.h>

/* arbitrary maximum size */
#define MAXSIZE 5000*10000 

void print_matrix( const gsl_matrix *, char *, int, int );
void fprint_matrix( char *, const gsl_matrix *, char *, int, int );
gsl_matrix * fread_matrix( char *, int *, int *, int );