[next] [prev] [prev-tail] [tail] [up]

4 Linear Algebra

4.1 Vectors

4.1.1 Overview

The structures and functions related to vectors are declared in pnl/pnl_vector.h.

Vectors are declared for several basic types : double, int, and dcomplex. In the following declarations, BASE must be replaced by one the previous types and the corresponding vector structures are respectively named PnlVect, PnlVectInt, PnlVectComplex 

typedef struct _PnlVect { 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlVect pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int size; /*!< size of the vector */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  double *array; /*!< pointer to store the data */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlVect; 
 
typedef struct _PnlVectInt { 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlVectInt pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int size; /*!< size of the vector */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  int *array; /*!< pointer to store the data */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlVectInt; 
 
typedef struct _PnlVectComplex { 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlVectComplex pointer to be cast 
   * to a PnlObject 
   */ 
  PnlObject object; 
  int size; /*!< size of the vector */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  dcomplex *array; /*!< pointer to store the data */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlVectComplex;

size is the size of the vector, array is a pointer containing the data and owner is an integer to know if the vector owns its array pointer (owner=1) or shares it with another structure (owner=0). mem_size is the number of elements the vector can hold at most.

4.1.2 Functions

General functions These functions exist for all types of vector no matter what the basic type is. The following conventions are used to name functions operating on vectors. Here is the table of prefixes used for the different basic types.

type prefix BASE
double pnl_vect double
int pnl_vect_int int
dcomplexpnl_vect_complexdcomplex

In this paragraph, we present the functions operating on PnlVect which exist for all types. To deduce the prototypes of these functions for other basic types, one must replace pnl_vect and double according the above table.

Constructors and destructors There are no special functions to access the size of a vector, instead the field size should be accessed directly.

Resizing vectors

Accessing elements If it is supported by the compiler, the following functions are declared inline. To speed up these functions, you can define the macro PNL_RANGE_CHECK_OFF, see Section 1.3.2 for an explanation.

Accessing elements of a vector is faster using the following macros

Printing vector

Applying external operation to vectors

Element wise operations

Scalar products and norms

Comparison functions

Ordering functions The following functions are not defined for PnlVectComplex because there is no total ordering on Complex numbers

Misc

Complex vector functions

There exist functions to directly access the real or imaginary parts of an element of a complex vector. These functions also have inlined versions that are used if the variable HAVE_INLINE was declared at compilation time.

Equivalently to these functions, there exist macros. When the compiler is able to handle inline code, there is no gain in using macros instead of inlined functions at least in principle.

4.2 Compact Vectors

4.2.1 Short description
typedef struct PnlVectCompact { 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlVectCompact pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int size; /* size of the vector */ 
  double val; /* single value */ 
  double *array; /* Pointer to double values */ 
  char convert; /* ’a’, ’d’ : array, double */ 
} PnlVectCompact;

4.2.2 Functions

4.3 Matrices

4.3.1 Overview

The structures and functions related to matrices are declared in pnl/pnl_matrix.h. 

typedef struct _PnlMat{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< nb rows */ 
  int n; /*!< nb columns */ 
  int mn; /*!< product m*n */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  double *array; /*!< pointer to store the data row-wise */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlMat; 
 
typedef struct _PnlMatInt{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlMatInt pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< nb rows */ 
  int n; /*!< nb columns */ 
  int mn; /*!< product m*n */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  int *array; /*!< pointer to store the data row-wise */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlMatInt; 
 
typedef struct _PnlMatComplex{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlMatComplex pointer to be cast 
   * to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< nb rows */ 
  int n; /*!< nb columns */ 
  int mn; /*!< product m*n */ 
  int mem_size; /*!< size of the memory block allocated for array */ 
  dcomplex *array; /*!< pointer to store the data row-wise */ 
  int owner; /*!< 1 if the object owns its array member, 0 otherwise */ 
} PnlMatComplex;

m is the number of rows, n is the number of columns. array is a pointer containing the data of the matrix stored line wise, The element (i, j) of the matrix is array[i*m+j]. owner is an integer to know if the matrix owns its array pointer (owner=1) or shares it with another structure (owner=0). mem_size is the number of elements the matrix can hold at most.

The following operations are implemented on matrices and vectors. alpha and beta are numbers, A and B are matrices and x and y are vectors.

pnl_mat_axpy B := alpha * A + B
pnl_mat_scalar_prod x’ A y
pnl_mat_dgemm C := alpha * op (A) * op (B) + beta * C
pnl_mat_mult_vect_transpose_inplacey = A’ * x
pnl_mat_mult_vect_inplace y = A * x
pnl_mat_lAxpby y := lambda * A * x + beta * y
pnl_mat_dgemv y := alpha * op (A) * x + beta * y
pnl_mat_dger A := alpha x * y’ + A

4.3.2 Generic Functions

These functions exist for all types of matrices no matter what the basic type is. The following conventions are used to name functions operating on matrices. Here is the table of prefixes used for the different basic types.

type prefix BASE
double pnl_mat double
int pnl_mat_int int
dcomplexpnl_mat_complexdcomplex

In this paragraph we present the functions operating on PnlMat which exist for all types. To deduce the prototypes of these functions for other basic types, one must replace pnl_mat and double according the above table.

Constructors and destructors There are no special functions to access the sizes of a matrix, instead the fields m, n and mn give direct access to the number of rows, columns and the size of the matrix.

Accessing elements. If it is supported by the compiler, the following functions are declared inline. To speed up these functions, you can define the macro PNL_RANGE_CHECK_OFF, see Section 1.3.2 for an explanation.

Accessing elements of a matrix is faster using the following macros

Printing Matrices

Applying external operations

Element wise operations

Comparison functions

Ordering operations

Standard matrix operations

Linear systems and matrix decompositions The following functions are designed to solve linear system of the from A x = b where A is a matrix and b is a vector except in the functions pnl_mat_syslin_mat, pnl_mat_lu_syslin_mat and pnl_mat_chol_syslin_mat which expect the right hand side member to be a matrix too. Whenever the vector b is not needed once the system is solved, you should consider using “inplace” functions.

All the functions described in this paragraph return OK if the computations have been carried out successfully and FAIL otherwise.

The following functions are designed to invert matrices. The authors provide these functions although they cannot find good reasons to use them. Note that to solve a linear system, one must used the syslin functions and not invert the system matrix because it is much longer.

4.3.3 Functions specific to base type double

Linear systems and matrix decompositions The following functions are designed to solve linear system of the from A x = b where A is a matrix and b is a vector except in the functions pnl_mat_syslin_mat, pnl_mat_lu_syslin_mat and pnl_mat_chol_syslin_mat which expect the right hand side member to be a matrix too. Whenever the vector b is not needed once the system is solved, you should consider using “inplace” functions.

All the functions described in this paragraph return OK if the computations have been carried out successfully and FAIL otherwise.

4.3.4 Functions specific to base type dcomplex

4.3.5 Permutations
typedef PnlVectInt PnlPermutation;

The PnlPermutation type is actually nothing else than a vector of integers, i.e. a PnlVectInt. It is used to store the partial pivoting with row interchanges transformation needed in the LU decomposition. We use the Blas convention for storing permutations. Consider a PnlPermutation p generated by a LU decomposition of a matrix A : to compute the decomposition, row i of A was interchanged with row p(i).

4.4 Tridiagonal Matrices

4.4.1 Overview

The structures and functions related to tridiagonal matrices are declared in pnl/pnl_tridiag_matrix.h.

We only store the three main diagonals as three vectors. 

typedef struct PnlTridiagMat{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlTridiagMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int size; /*!< number of rows, the matrix must be square */ 
  double *D; /*!< diagonal elements */ 
  double *DU; /*!< upper diagonal elements */ 
  double *DL; /*!< lower diagonal elements */ 
} PnlTridiagMat;

size is the size of the matrix, D is an array of size size containing the diagonal terms. DU, DL are two arrays of size size-1 containing respectively the upper diagonal (Mi,i+1) and the lower diagonal (Mi-1,i). 

typedef struct PnlTridiagMatLU{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlTridiagMatLU pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int size; /*!< number of rows, the matrix must be square */ 
  double *D; /*!< diagonal elements */ 
  double *DU; /*!< upper diagonal elements */ 
  double *DU2; /*!< second upper diagonal elements */ 
  double *DL; /*!< lower diagonal elements */ 
  int *ipiv; /*!< Permutation: row i has been interchanged with row ipiv(i) */ 
};

This type is used to store the LU decomposition of a tridiagonal matrix.

4.4.2 Functions

Constructors and destructors

Accessing elements. If it is supported by the compiler, the following functions are declared inline. To speed up these functions, you can use the macro constant PNL_RANGE_CHECK_OFF, see Section 1.3.2 for an explanation.

Printing Matrix

Algebra operations

Element-wise operations

Standard matrix operations & Linear systems

4.5 Band Matrices

4.5.1 Overview
typedef struct 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlBandMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< nb rows */ 
  int n; /*!< nb columns */ 
  int nu; /*!< nb of upperdiagonals */ 
  int nl; /*!< nb of lowerdiagonals */ 
  int m_band; /*!< nb rows of the band storage */ 
  int n_band; /*!< nb columns of the band storage */ 
  double *array;  /*!< a block to store the bands */ 
} PnlBandMat;

The structures and functions related to band matrices are declared in pnl/pnl_band_matrix.h.

4.5.2 Functions

Constructors and destructors

Accessing elements. If it is supported by the compiler, the following functions are declared inline. To speed up these functions, you can use the macro constant PNL_RANGE_CHECK_OFF, see Section 1.3.2 for an explanation.

Element wise operations

Standard matrix operations & Linear system

4.6 Sparse Matrices

4.6.1 Short description

The structures and functions related to matrices are declared in pnl/pnl_sp_matrix.h. 

typedef struct _PnlSpMat 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows a PnlSpMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< number of rows */ 
  int n; /*!< number of columns */ 
  int nz; /*!< number of non-zero elements */ 
  int *J; /*!< column indices, vector of size nzmax */ 
  int *I; /*!< row offset integer vector, 
            array[I[i]] is the first element of row i. 
            Vector of size (m+1) */ 
  double *array; /*!< pointer to store the data of size nzmax*/ 
  int nzmax; /*!< size of the memory block allocated for array */ 
} PnlSpMat; 
 
typedef struct _PnlSpMatInt 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows a PnlSpMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< number of rows */ 
  int n; /*!< number of columns */ 
  int nz; /*!< number of non-zero elements */ 
  int *J; /*!< column indices, vector of size nzmax */ 
  int *I; /*!< row offset integer vector, 
            array[I[i]] is the first element of row i. 
            Vector of size (m+1) */ 
  int *array; /*!< pointer to store the data of size nzmax */ 
  int nzmax; /*!< size of the memory block allocated for array */ 
} PnlSpMatInt; 
 
typedef struct _PnlSpMatComplex 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows a PnlSpMat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int m; /*!< number of rows */ 
  int n; /*!< number of columns */ 
  int nz; /*!< number of non-zero elements */ 
  int *J; /*!< column indices, vector of size nzmax */ 
  int *I; /*!< row offset integer vector, 
            array[I[i]] is the first element of row i. 
            Vector of size (m+1) */ 
  dcomplex *array; /*!< pointer to store the data of size nzmax */ 
  int nzmax; /*!< size of the memory block allocated for array */ 
} PnlSpMatComplex;

The non zero elements of row i are stored in array between the indices I[i] and I[i+1]-1. The array J contains the column indices of every element of array.

Sparse matrices are defined using the internal template approach and can be used for integer, float or complex base data according to the following table

base typeprefix type
double pnl_sp_mat PnlSpMat
int pnl_sp_mat_int PnlSpMatInt
dcomplex pnl_sp_mat_complexPnlSpMatComplex

4.6.2 Functions

Constructors and destructors

Accessing elements

Applying external operations

Standard matrix operations

Comparison functions

4.7 Hyper Matrices

4.7.1 Short description

The Hyper matrix types and related functions are defined in the header pnl/pnl_matrix.h. 

typedef struct PnlHmat{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlHmat pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int ndim; /*!< nb dimensions */ 
  int *dims; /*!< pointer to store the values of the ndim dimensions */ 
  int mn; /*!< product dim_1 *...*dim_ndim */ 
  int *pdims; /*!< array of size ndim, s.t. pdims[i] = dims[ndim-1] x ... dims[i+1] 
                with pdims[ndim - 1] = 1 */ 
  double *array; /*!< pointer to store */ 
} PnlHmat; 
 
typedef struct PnlHmatInt{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlHmatInt pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int ndim; /*!< nb dimensions */ 
  int *dims; /*!< pointer to store the value of the ndim dimensions */ 
  int mn; /*!< product dim_1 *...*dim_ndim */ 
  int *pdims; /*!< array of size ndim, s.t. pdims[i] = dims[ndim-1] x ... dims[i+1] 
                with pdims[ndim - 1] = 1 */ 
  int *array; /*!< pointer to store */ 
} PnlHmatInt; 
 
typedef struct PnlHmatComplex{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlHmatComplex pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int ndim; /*!< nb dimensions */ 
  int *dims; /*!< pointer to store the value of the ndim dimensions */ 
  int mn; /*!< product dim_1 *...*dim_ndim */ 
  int *pdims; /*!< array of size ndim, s.t. pdims[i] = dims[ndim-1] x ... dims[i+1] 
                with pdims[ndim - 1] = 1 */ 
  dcomplex *array; /*!< pointer to store */ 
} PnlHmatComplex;

ndim is the number of dimensions, dim is an array to store the size of each dimension and nm contains the product of the sizes of each dimension. array is an array of size mn containing the data. The integer array pdims is used to create the one–to–one map between the natural indexing and the linear indexing used in array.

4.7.2 Functions

These functions exist for all types of hypermatrices no matter what the basic type is. The following conventions are used to name functions operating on hypermatrices. Here is the table of prefixes used for the different basic types.

base typeprefix type
double pnl_hmat PnlHmat
int pnl_hmat_int PnlHmatInt
dcomplex pnl_hmat_complexPnlHmatComplex

In this paragraph, we present the functions operating on PnlHmat which exist for all types. To deduce the prototypes of these functions for other basic types, one must replace pnl_hmat and double according the above table.

Constructors and destructors

Accessing elements

Printing hypermatrices

Term by term operations

4.8 Iterative Solvers

4.8.1 Overview

The structures and functions related to solvers are declared in pnl/pnl_linalgsolver.h. 

typedef struct _PnlIterationBase PnlIterationBase; 
typedef struct _PnlCgSolver PnlCgSolver; 
typedef struct _PnlBicgSolver PnlBicgSolver; 
typedef struct _PnlGmresSolver PnlGmresSolver; 
 
struct _PnlIterationBase 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlVectXXX pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int iteration; 
  int max_iter; 
  double normb; 
  double tol_; 
  double resid; 
  int error; 
  /* char *  err_msg; */ 
}; 
 
/* When you repeatedly use iterative solvers, do not malloc each time */ 
struct _PnlCgSolver 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlCgSolver  pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  PnlVect * r; 
  PnlVect * z; 
  PnlVect * p; 
  PnlVect * q; 
  double rho; 
  double oldrho; 
  double beta; 
  double alpha; 
  PnlIterationBase * iter; 
} ; 
 
struct _PnlBicgSolver 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlBicgSolver pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  double rho_1, rho_2, alpha, beta, omega; 
  PnlVect * p; 
  PnlVect * phat; 
  PnlVect * s; 
  PnlVect * shat; 
  PnlVect * t; 
  PnlVect * v; 
  PnlVect * r; 
  PnlVect *  rtilde; 
  PnlIterationBase * iter; 
} ; 
 
struct _PnlGmresSolver 
{ 
  /** 
   * Must be the first element in order for the object mechanism to work 
   * properly. This allows any PnlGmresSolver pointer to be cast to a PnlObject 
   */ 
  PnlObject object; 
  int restart; 
  double beta; 
  PnlVect * s; 
  PnlVect * cs; 
  PnlVect * sn; 
  PnlVect * w; 
  PnlVect * r; 
  PnlMat * H; 
  PnlVect * v[MAX_RESTART]; 
  PnlIterationBase *iter; 
  PnlIterationBase *iter_inner; 
} ;

A Left preconditioner solves the problem :

PM x = P b,

and whereas right preconditioner solves

M Py = b, P y = x.

More information is given in Saad, Yousef (2003). Iterative methods for sparse linear systems (2nd ed. ed.). SIAM. ISBN 0898715342. OCLC 51266114. The reader will find in this book some discussion about right or/and left preconditioner and a description of the following algorithms.

These algorithms, we implemented with a left preconditioner. Right preconditioner can be easily computed changing matrix vector multiplication operator from M x to M PR x and solving PRy = x at the end of algorithm.

4.8.2 Functions

Three methods are implemented : Conjugate Gradient, BICGstab and GMRES with restart. For each of them a structure is created to store temporary vectors used in the algorithm. In some cases, we have to apply iterative methods more than once : for example to solve at each time step a discrete form of an elliptic problem come from parabolic problem. In the cases, do not call the constructor and destructor at each time, but instead use the initialization and solve procedures.

Formally we have, 

Create iterative method 
For each time step 
  Initialisation of iterative method 
  Solve linear system link to elliptic problem 
end for 
free iterative method

In these functions, we don’t use any particular matrix structure. We give the matrix vector multiplication as a parameter of the solver.

Conjugate Gradient method Only available for symmetric and positive matrices.

BICG stab

GMRES with restart See Saad, Yousef (2003) for a discussion about the restart parameter. For GMRES we need to store at the p-th iteration p vectors of the same size of the right and side. It could be very expensive in term of memory allocation. So GMRES with restart algorithm stop if p = restart and restarts the algorithm with the previously computed solution as initial guess.

Note that if restart equals m, we have a classical GMRES algorithm.

In the next paragraph, we write all the solvers for PnlMat . This will be done as follows: construct an application matrix vector. 

static void pnl_mat_mult_vect_applied(const void *mat, const PnlVect *vec, 
                                      const double a , const double b, 
                                      PnlVect *lhs) 
{pnl_mat_lAxpby(a, (PnlMat*)mat, vec, b, lhs);}

and give it as the parameter of the iterative method 

int pnl_mat_cg_solver_solve(const PnlMat * Matrix, const PnlMat * PC, 
                            PnlVect * x, const PnlVect *b, PnlCgSolver * Solver) 
{ return pnl_cg_solver_solve(pnl_mat_mult_vect_applied, 
                             Matrix, pnl_mat_mult_vect_applied, 
                             PC, x, b, Solver);}

In practice, we cannot define all iterative methods for all structures. With this implementation, the user can easily :

Iterative algorithms for PnlMat

[next] [prev] [prev-tail] [front] [up]