<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Here is an example I put together that uses the Gnu Scientific
Library to do complex matrix arithmetic using the GSL</p>
<p>#include <stdio.h><br>
#include <math.h><br>
#include <gsl/gsl_blas.h><br>
#include <gsl/gsl_matrix.h><br>
#include <gsl/gsl_complex.h><br>
#include <gsl/gsl_complex_math.h><br>
#include <gsl/gsl_linalg.h><br>
<br>
// Compile with gcc gslztest.c -lgsl -lm <br>
<br>
int main (void)<br>
{<br>
<br>
const int N = 10;<br>
int i, j, s;<br>
double pi, NCUBE;<br>
double alpha, beta;<br>
gsl_complex sumdiag;<br>
<br>
<br>
// Create pointers and allocate space for complex matrices<br>
gsl_matrix_complex *A = NULL;<br>
gsl_matrix_complex *B = NULL;<br>
gsl_matrix_complex *C = NULL;<br>
gsl_matrix_complex *INV = NULL;<br>
<br>
A = gsl_matrix_complex_alloc(N,N);<br>
B = gsl_matrix_complex_alloc(N,N);<br>
C = gsl_matrix_complex_alloc(N,N);<br>
INV = gsl_matrix_complex_alloc(N,N);<br>
<br>
// Also need a permutation matrix for PLU algorithm<br>
<br>
gsl_permutation *P = gsl_permutation_alloc(N);<br>
<br>
// Set the value of Pi to machine precision<br>
<br>
pi = acos(-1.0);<br>
<br>
// Set the value of N^3 since we use is a lot<br>
<br>
NCUBE = (double) N * (double) N * (double) N;<br>
<br>
// Now fill the matrices with the same info<br>
<br>
for (i=1; i<=N; i++)<br>
for (j=1; j<=N; j++) {<br>
alpha = sqrt( pow(cos ( (double) ( i*i*j) * pi /
NCUBE ),2) +<br>
pow(sin ( (double) (2*i*i*j) * pi /
NCUBE ),2));<br>
alpha = (double) 1.0 / alpha;<br>
gsl_matrix_complex_set( A, i-1, j-1,<br>
gsl_complex_rect(<br>
alpha * cos( (double) ( i*i*j) * pi /
NCUBE ),<br>
alpha * sin( (double) (2*i*i*j) * pi /
NCUBE )));<br>
gsl_matrix_complex_set( B, i-1, j-1,<br>
gsl_matrix_complex_get(A, i-1, j-1));<br>
}<br>
<br>
// Use the gsl routines to compute the inverse of B and place
it in INV<br>
<br>
gsl_linalg_complex_LU_decomp( B, P, &s);<br>
gsl_linalg_complex_LU_invert( B, P, INV);<br>
<br>
// Now multiply A and its inverse (stored in INV) together
using<br>
// blas routines <br>
<br>
gsl_blas_zgemm( CblasNoTrans, CblasNoTrans,<br>
GSL_COMPLEX_ONE, A, INV,<br>
GSL_COMPLEX_ZERO, C );<br>
<br>
sumdiag = GSL_COMPLEX_ZERO;<br>
for (i=0;i<N;i++) sumdiag = gsl_complex_add(sumdiag,<br>
gsl_matrix_complex_get(C, i, i));<br>
<br>
<br>
printf("Sum of diagonal is %g + i %g\n", GSL_REAL(sumdiag),
GSL_IMAG(sumdiag));<br>
<br>
}<br>
<br>
</p>
<div class="moz-signature">-- <br>
<b>Andrew J. Pounds, Ph.D.</b><br>
<i>Professor of Chemistry and Computer Science</i><br>
<i>Director of the Computational Science Program</i><br>
<i>Mercer University, Macon, GA, 31207 (478) 301-5627 </i></div>
</body>
</html>