6.4 Linear Systems 181
#include"lib.h"
using namespacestd;
/function declarations/
voidinverse(double*,int);
/
This program sets up a simple 3x3 symmetric matrix
and finds its determinant and inverse
*/
intmain()
{
int i, j, k, result, n = 3;
double matr, sum,
a[3][3] ={ {1.0, 3.0, 4.0},
{3.0, 4.0, 6.0},
{4.0, 6.0, 8.0}};
// memory for inverse matrix
matr = (double) matrix(n, n,sizeof(double));
// various print statements in the original code are omitted
inverse(matr, n);// calculate and return inverse matrix
....
return0;
}// End: function main()
In order to use the program library you need to include thelib.hfile using the#include"lib.h"
statement. This function utilizes the library functionmatrixandfree_matrixto allocate and
free memory during execution. The matrixa[ 3 ][ 3 ]is set at compilation time. Alternatively, you
could have used either Blitz++ or Armadillo.
The corresponding Fortran program for the inverse of a matrix reads
http://folk.uio.no/compphys/programs/FYS3150/f90library/f90lib.f90
!
! Routines to do mtx inversion, from Numerical
! Recipes, Teukolsky et al. Routines included
! below are MATINV, LUDCMP and LUBKSB. See chap 2
! of Numerical Recipes for further details
!
SUBROUTINEmatinv(a,n, indx, d)
IMPLICITNONE
INTEGER,INTENT(IN) :: n
INTEGER:: i, j
REAL(DP),DIMENSION(n,n),INTENT(INOUT) :: a
REAL(DP),ALLOCATABLE:: y(:,:)
REAL(DP) :: d
INTEGER, ,INTENT(INOUT) :: indx(n)
ALLOCATE(y( n, n))
y=0.
! setup identity matrix
DOi=1,n
y(i,i)=1.
ENDDO
! LU decompose the matrix just once
CALLlu_decompose(a,n,indx,d)
! Find inverse by columns