[CSC 335] Using C and the GSL to answer problem 2 on the homework

[Andrew J. Pounds]

Here is an example I put together that uses the Gnu Scientific Library 
to do complex matrix arithmetic using the GSL

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_linalg.h>

// Compile with gcc gslztest.c -lgsl -lm

int main (void)
{

     const int N = 10;
     int i, j, s;
     double pi, NCUBE;
     double alpha, beta;
     gsl_complex sumdiag;


     // Create pointers and allocate space for complex matrices
     gsl_matrix_complex *A   = NULL;
     gsl_matrix_complex *B   = NULL;
     gsl_matrix_complex *C   = NULL;
     gsl_matrix_complex *INV = NULL;

     A   = gsl_matrix_complex_alloc(N,N);
     B   = gsl_matrix_complex_alloc(N,N);
     C   = gsl_matrix_complex_alloc(N,N);
     INV = gsl_matrix_complex_alloc(N,N);

     // Also need a permutation matrix for PLU algorithm

     gsl_permutation *P = gsl_permutation_alloc(N);

     // Set the value of Pi to machine precision

     pi = acos(-1.0);

     // Set the value of N^3 since we use is a lot

     NCUBE = (double) N * (double) N * (double) N;

     // Now fill the matrices with the same info

     for (i=1; i<=N; i++)
         for (j=1; j<=N; j++) {
             alpha =  sqrt( pow(cos ( (double) (  i*i*j) * pi / NCUBE ),2) +
                            pow(sin ( (double) (2*i*i*j) * pi / NCUBE ),2));
             alpha = (double) 1.0 / alpha;
             gsl_matrix_complex_set( A, i-1, j-1,
                     gsl_complex_rect(
                         alpha * cos( (double) (  i*i*j) * pi / NCUBE ),
                         alpha * sin( (double) (2*i*i*j) * pi / NCUBE )));
             gsl_matrix_complex_set( B, i-1, j-1,
                     gsl_matrix_complex_get(A, i-1, j-1));
         }

     // Use the gsl routines to compute the inverse of B and place it in INV

     gsl_linalg_complex_LU_decomp( B, P, &s);
     gsl_linalg_complex_LU_invert( B, P, INV);

     // Now multiply A and its inverse (stored in INV) together using
     // blas routines

     gsl_blas_zgemm( CblasNoTrans, CblasNoTrans,
             GSL_COMPLEX_ONE, A, INV,
             GSL_COMPLEX_ZERO, C );

     sumdiag = GSL_COMPLEX_ZERO;
     for (i=0;i<N;i++) sumdiag = gsl_complex_add(sumdiag,
             gsl_matrix_complex_get(C, i, i));


     printf("Sum of diagonal is %g + i %g\n", GSL_REAL(sumdiag), 
GSL_IMAG(sumdiag));

}

-- 
*Andrew J. Pounds, Ph.D.*
/Professor of Chemistry and Computer Science/
/Director of the Computational Science Program/
/Mercer University, Macon, GA, 31207 (478) 301-5627 /
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://theochem.mercer.edu/pipermail/csc335/attachments/20231114/820d8054/attachment.html>