/*******************************************************************
* *
* File : smalldense.h *
* Programmers : Scott D. Cohen and Alan C. Hindmarsh @ LLNL *
* Version of : 26 June 2002 *
*-----------------------------------------------------------------*
* Copyright (c) 2002, The Regents of the University of California *
* Produced at the Lawrence Livermore National Laboratory *
* All rights reserved *
* For details, see sundials/shared/LICENSE *
*-----------------------------------------------------------------*
* This is the header file for a generic DENSE linear solver *
* package, intended for small dense matrices. These routines *
* use the type realtype** for dense matrix arguments. *
* *
* These routines begin with "den" (except for the factor and *
* solve routines which are called gefa and gesl, respectively). *
* The underlying matrix storage is described in the *
* documentation for denalloc. *
* *
*******************************************************************/
#ifdef __cplusplus /* wrapper to enable C++ usage */
extern "C" {
#endif
#ifndef _smalldense_h
#define _smalldense_h
#include "sundialstypes.h"
/******************************************************************
* *
* Function : denalloc *
* Usage : realtype **a; *
* a = denalloc(n); *
* if (a == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denalloc(n) allocates storage for an n by n dense matrix. It *
* returns a pointer to the newly allocated storage if *
* successful. If the memory request cannot be satisfied, then *
* denalloc returns NULL. The underlying type of the dense matrix *
* returned is realtype **. If we allocate a dense matrix *
* realtype **a by a = denalloc(n), then a[j][i] references the *
* (i,j)th element of the matrix a, 0 <= i,j <= n-1, and a[j] is *
* a pointer to the first element in the jth column of a. *
* The location a[0] contains a pointer to n^2 contiguous *
* locations which contain the elements of a. *
* *
******************************************************************/
realtype **denalloc(integertype n);
/******************************************************************
* *
* Function : denallocpiv *
* Usage : integertype *pivot; *
* pivot = denallocpiv(n); *
* if (pivot == NULL) ... memory request failed *
*----------------------------------------------------------------*
* denallocpiv(n) allocates an array of n integertype. It returns *
* a pointer to the first element in the array if successful. *
* It returns NULL if the memory request could not be satisfied. *
* *
******************************************************************/
integertype *denallocpiv(integertype n);
/******************************************************************
* *
* Function : gefa *
* Usage : integertype ier; *
* ier = gefa(a,n,p); *
* if (ier > 0) ... zero element encountered during *
* the factorization *
*----------------------------------------------------------------*
* gefa(a,n,p) factors the n by n dense matrix a. It overwrites *
* the elements of a with its LU factors and keeps track of the *
* pivot rows chosen in the pivot array p. *
* *
* A successful LU factorization leaves the matrix a and the *
* pivot array p with the following information: *
* *
* (1) p[k] contains the row number of the pivot element chosen *
* at the beginning of elimination step k, k=0, 1, ..., n-1. *
* *
* (2) If the unique LU factorization of a is given by Pa = LU, *
* where P is a permutation matrix, L is a lower triangular *
* matrix with all 1's on the diagonal, and U is an upper *
* triangular matrix, then the upper triangular part of a *
* (including its diagonal) contains U and the strictly lower *
* triangular part of a contains the multipliers, I-L. *
* *
* gefa returns 0 if successful. Otherwise it encountered a zero *
* diagonal element during the factorization. In this case it *
* returns the column index (numbered from one) at which it *
* encountered the zero. *
* *
******************************************************************/
integertype gefa(realtype **a, integertype n, integertype *p);
/******************************************************************
* *
* Function : gesl *
* Usage : realtype *b; *
* ier = gefa(a,n,p); *
* if (ier == 0) gesl(a,n,p,b); *
*----------------------------------------------------------------*
* gesl(a,n,p,b) solves the n by n linear system ax = b. It *
* assumes that a has been LU factored and the pivot array p has *
* been set by a successful call to gefa(a,n,p). The solution x *
* is written into the b array. *
* *
******************************************************************/
void gesl(realtype **a, integertype n, integertype *p, realtype *b);
/******************************************************************
* *
* Function : denzero *
* Usage : denzero(a,n); *
*----------------------------------------------------------------*
* denzero(a,n) sets all the elements of the n by n dense matrix *
* a to be 0.0. *
* *
******************************************************************/
void denzero(realtype **a, integertype n);
/******************************************************************
* *
* Function : dencopy *
* Usage : dencopy(a,b,n); *
*----------------------------------------------------------------*
* dencopy(a,b,n) copies the n by n dense matrix a into the *
* n by n dense matrix b. *
* *
******************************************************************/
void dencopy(realtype **a, realtype **b, integertype n);
/******************************************************************
* *
* Function : denscale *
* Usage : denscale(c,a,n); *
*----------------------------------------------------------------*
* denscale(c,a,n) scales every element in the n by n dense *
* matrix a by c. *
* *
******************************************************************/
void denscale(realtype c, realtype **a, integertype n);
/******************************************************************
* *
* Function : denaddI *
* Usage : denaddI(a,n); *
*----------------------------------------------------------------*
* denaddI(a,n) increments the n by n dense matrix a by the *
* identity matrix. *
* *
******************************************************************/
void denaddI(realtype **a, integertype n);
/******************************************************************
* *
* Function : denfreepiv *
* Usage : denfreepiv(p); *
*----------------------------------------------------------------*
* denfreepiv(p) frees the pivot array p allocated by *
* denallocpiv. *
* *
******************************************************************/
void denfreepiv(integertype *p);
/******************************************************************
* *
* Function : denfree *
* Usage : denfree(a); *
*----------------------------------------------------------------*
* denfree(a) frees the dense matrix a allocated by denalloc. *
* *
******************************************************************/
void denfree(realtype **a);
/******************************************************************
* *
* Function : denprint *
* Usage : denprint(a,n); *
*----------------------------------------------------------------*
* denprint(a,n) prints the n by n dense matrix a to standard *
* output as it would normally appear on paper. It is intended as *
* a debugging tool with small values of n. The elements are *
* printed using the %g option. A blank line is printed before *
* and after the matrix. *
* *
******************************************************************/
void denprint(realtype **a, integertype n);
#endif
#ifdef __cplusplus
}
#endif
