### 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++)
        for (j = 0; j < 3; j++)
          printf("m(%d,%d) = %g\n", i, j,
                 gsl_matrix_get (m, i, j));

      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`.

    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
    IOT trap/Abort (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", "w");
         gsl_matrix_fwrite (f, m);
         fclose (f);
      }

      {
         FILE * f = fopen("test.dat", "r");
         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++;
          }

      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);
    }

Here is the output of the program, which can be confirmed using GNU OCTAVE,

    $ ./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

    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