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
|
