8.4.13 Example programs for matrices

The program below shows how to allocate, initialize and read from a matrix using the functions gsl_matrix_alloc, gsl_matrix_set and gsl_matrix_get.

     #include <stdio.h>
     #include <gsl/gsl_matrix.h>
     
     int
     main (void)
     {
       int i, j; 
       gsl_matrix * m = gsl_matrix_alloc (10, 3);
       
       for (i = 0; i < 10; i++)
         for (j = 0; j < 3; j++)
           gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
       
       for (i = 0; i < 100; i++)  /* OUT OF RANGE ERROR */
         for (j = 0; j < 3; j++)
           printf ("m(%d,%d) = %g\n", i, j, 
                   gsl_matrix_get (m, i, j));
     
       gsl_matrix_free (m);
     
       return 0;
     }

Here is the output from the program. The final loop attempts to read outside the range of the matrix m, and the error is trapped by the range-checking code in gsl_matrix_get.

     $ ./a.out
     m(0,0) = 0.23
     m(0,1) = 1.23
     m(0,2) = 2.23
     m(1,0) = 100.23
     m(1,1) = 101.23
     m(1,2) = 102.23
     ...
     m(9,2) = 902.23
     gsl: matrix_source.c:13: ERROR: first index out of range
     Default GSL error handler invoked.
     Aborted (core dumped)

The next program shows how to write a matrix to a file.

     #include <stdio.h>
     #include <gsl/gsl_matrix.h>
     
     int
     main (void)
     {
       int i, j, k = 0; 
       gsl_matrix * m = gsl_matrix_alloc (100, 100);
       gsl_matrix * a = gsl_matrix_alloc (100, 100);
       
       for (i = 0; i < 100; i++)
         for (j = 0; j < 100; j++)
           gsl_matrix_set (m, i, j, 0.23 + i + j);
     
       {  
          FILE * f = fopen ("test.dat", "wb");
          gsl_matrix_fwrite (f, m);
          fclose (f);
       }
     
       {  
          FILE * f = fopen ("test.dat", "rb");
          gsl_matrix_fread (f, a);
          fclose (f);
       }
     
       for (i = 0; i < 100; i++)
         for (j = 0; j < 100; j++)
           {
             double mij = gsl_matrix_get (m, i, j);
             double aij = gsl_matrix_get (a, i, j);
             if (mij != aij) k++;
           }
     
       gsl_matrix_free (m);
       gsl_matrix_free (a);
     
       printf ("differences = %d (should be zero)\n", k);
       return (k > 0);
     }

After running this program the file test.dat should contain the elements of m, written in binary format. The matrix which is read back in using the function gsl_matrix_fread should be exactly equal to the original matrix.

The following program demonstrates the use of vector views. The program computes the column norms of a matrix.

     #include <math.h>
     #include <stdio.h>
     #include <gsl/gsl_matrix.h>
     #include <gsl/gsl_blas.h>
     
     int
     main (void)
     {
       size_t i,j;
     
       gsl_matrix *m = gsl_matrix_alloc (10, 10);
     
       for (i = 0; i < 10; i++)
         for (j = 0; j < 10; j++)
           gsl_matrix_set (m, i, j, sin (i) + cos (j));
     
       for (j = 0; j < 10; j++)
         {
           gsl_vector_view column = gsl_matrix_column (m, j);
           double d;
     
           d = gsl_blas_dnrm2 (&column.vector);
     
           printf ("matrix column %d, norm = %g\n", j, d);
         }
     
       gsl_matrix_free (m);
     
       return 0;
     }

Here is the output of the program,

     $ ./a.out
     matrix column 0, norm = 4.31461
     matrix column 1, norm = 3.1205
     matrix column 2, norm = 2.19316
     matrix column 3, norm = 3.26114
     matrix column 4, norm = 2.53416
     matrix column 5, norm = 2.57281
     matrix column 6, norm = 4.20469
     matrix column 7, norm = 3.65202
     matrix column 8, norm = 2.08524
     matrix column 9, norm = 3.07313

The results can be confirmed using gnu octave,

     $ octave
     GNU Octave, version 2.0.16.92
     octave> m = sin(0:9)' * ones(1,10)
                    + ones(10,1) * cos(0:9);
     octave> sqrt(sum(m.^2))
     ans =
       4.3146  3.1205  2.1932  3.2611  2.5342  2.5728
       4.2047  3.6520  2.0852  3.0731