File: Example-programs-for-matrices.html

package info (click to toggle)
gsl-ref-html 1.16-1
links: PTS
area: non-free
in suites: jessie, jessie-kfreebsd, stretch
size: 5,816 kB
ctags: 4,130
sloc: makefile: 35
file content (240 lines) | stat: -rw-r--r-- 7,704 bytes
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 The GSL Team.

Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with no
Invariant Sections and no cover texts.  A copy of the license is
included in the section entitled "GNU Free Documentation License". -->
<!-- Created by GNU Texinfo 5.1, http://www.gnu.org/software/texinfo/ -->
<head>
<title>GNU Scientific Library &ndash; Reference Manual: Example programs for matrices</title>

<meta name="description" content="GNU Scientific Library &ndash; Reference Manual: Example programs for matrices">
<meta name="keywords" content="GNU Scientific Library &ndash; Reference Manual: Example programs for matrices">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Function-Index.html#Function-Index" rel="index" title="Function Index">
<link href="Matrices.html#Matrices" rel="up" title="Matrices">
<link href="Vector-and-Matrix-References-and-Further-Reading.html#Vector-and-Matrix-References-and-Further-Reading" rel="next" title="Vector and Matrix References and Further Reading">
<link href="Matrix-properties.html#Matrix-properties" rel="previous" title="Matrix properties">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.indentedblock {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smallindentedblock {margin-left: 3.2em; font-size: smaller}
div.smalllisp {margin-left: 3.2em}
kbd {font-style:oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nocodebreak {white-space:nowrap}
span.nolinebreak {white-space:nowrap}
span.roman {font-family:serif; font-weight:normal}
span.sansserif {font-family:sans-serif; font-weight:normal}
ul.no-bullet {list-style: none}
-->
</style>


</head>

<body lang="en" bgcolor="#FFFFFF" text="#000000" link="#0000FF" vlink="#800080" alink="#FF0000">
<a name="Example-programs-for-matrices"></a>
<div class="header">
<p>
Previous: <a href="Matrix-properties.html#Matrix-properties" accesskey="p" rel="previous">Matrix properties</a>, Up: <a href="Matrices.html#Matrices" accesskey="u" rel="up">Matrices</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Example-programs-for-matrices-1"></a>
<h4 class="subsection">8.4.13 Example programs for matrices</h4>

<p>The program below shows how to allocate, initialize and read from a matrix
using the functions <code>gsl_matrix_alloc</code>, <code>gsl_matrix_set</code> and
<code>gsl_matrix_get</code>.
</p>
<div class="example">
<pre class="verbatim">#include &lt;stdio.h&gt;
#include &lt;gsl/gsl_matrix.h&gt;

int
main (void)
{
  int i, j; 
  gsl_matrix * m = gsl_matrix_alloc (10, 3);
  
  for (i = 0; i &lt; 10; i++)
    for (j = 0; j &lt; 3; j++)
      gsl_matrix_set (m, i, j, 0.23 + 100*i + j);
  
  for (i = 0; i &lt; 100; i++)  /* OUT OF RANGE ERROR */
    for (j = 0; j &lt; 3; j++)
      printf (&quot;m(%d,%d) = %g\n&quot;, i, j, 
              gsl_matrix_get (m, i, j));

  gsl_matrix_free (m);

  return 0;
}
</pre></div>

<p>Here is the output from the program.  The final loop attempts to read
outside the range of the matrix <code>m</code>, and the error is trapped by
the range-checking code in <code>gsl_matrix_get</code>.
</p>
<div class="example">
<pre class="example">$ ./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)
</pre></div>

<p>The next program shows how to write a matrix to a file.
</p>
<div class="example">
<pre class="verbatim">#include &lt;stdio.h&gt;
#include &lt;gsl/gsl_matrix.h&gt;

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 &lt; 100; i++)
    for (j = 0; j &lt; 100; j++)
      gsl_matrix_set (m, i, j, 0.23 + i + j);

  {  
     FILE * f = fopen (&quot;test.dat&quot;, &quot;wb&quot;);
     gsl_matrix_fwrite (f, m);
     fclose (f);
  }

  {  
     FILE * f = fopen (&quot;test.dat&quot;, &quot;rb&quot;);
     gsl_matrix_fread (f, a);
     fclose (f);
  }

  for (i = 0; i &lt; 100; i++)
    for (j = 0; j &lt; 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 (&quot;differences = %d (should be zero)\n&quot;, k);
  return (k &gt; 0);
}
</pre></div>

<p>After running this program the file <samp>test.dat</samp> should contain the
elements of <code>m</code>, written in binary format.  The matrix which is read
back in using the function <code>gsl_matrix_fread</code> should be exactly
equal to the original matrix.
</p>
<p>The following program demonstrates the use of vector views.  The program
computes the column norms of a matrix.
</p>
<div class="example">
<pre class="verbatim">#include &lt;math.h&gt;
#include &lt;stdio.h&gt;
#include &lt;gsl/gsl_matrix.h&gt;
#include &lt;gsl/gsl_blas.h&gt;

int
main (void)
{
  size_t i,j;

  gsl_matrix *m = gsl_matrix_alloc (10, 10);

  for (i = 0; i &lt; 10; i++)
    for (j = 0; j &lt; 10; j++)
      gsl_matrix_set (m, i, j, sin (i) + cos (j));

  for (j = 0; j &lt; 10; j++)
    {
      gsl_vector_view column = gsl_matrix_column (m, j);
      double d;

      d = gsl_blas_dnrm2 (&amp;column.vector);

      printf (&quot;matrix column %d, norm = %g\n&quot;, j, d);
    }

  gsl_matrix_free (m);

  return 0;
}
</pre></div>

<p>Here is the output of the program, 
</p>
<div class="example">
<pre class="example">$ ./a.out
</pre><pre class="verbatim">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
</pre></div>

<p>The results can be confirmed using <small>GNU OCTAVE</small>,
</p>
<div class="example">
<pre class="example">$ octave
GNU Octave, version 2.0.16.92
octave&gt; m = sin(0:9)' * ones(1,10) 
               + ones(10,1) * cos(0:9); 
octave&gt; 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
</pre></div>


<hr>
<div class="header">
<p>
Previous: <a href="Matrix-properties.html#Matrix-properties" accesskey="p" rel="previous">Matrix properties</a>, Up: <a href="Matrices.html#Matrices" accesskey="u" rel="up">Matrices</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>



</body>
</html>