/*
 * -----------------------------------------------------------------
 * $Revision: 1.7.2.1 $
 * $Date: 2005/01/24 21:50:35 $
 * -----------------------------------------------------------------
 * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and
 *                Radu Serban @ LLNL
 * -----------------------------------------------------------------
 * 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.  The routines listed in this file all use type
 * DenseMat, defined below, for matrices.  These routines in turn
 * call routines in the smalldense.h/smalldense.c module, which
 * use the type realtype** for matrices.  This separation allows
 * for possible modifications in which matrices of type DenseMat
 * may not be stored contiguously, while small matrices can still
 * be treated with the routines in smalldense.
 *
 * Routines that work with the type DenseMat begin with "Dense".
 * The DenseAllocMat function allocates a dense matrix for use in
 * the other DenseMat routines listed in this file. Matrix
 * storage details are given in the documentation for the type
 * DenseMat. The DenseAllocPiv function allocates memory for
 * pivot information. The storage allocated by DenseAllocMat and
 * DenseAllocPiv is deallocated by the routines DenseFreeMat and
 * DenseFreePiv, respectively. The DenseFactor and DenseBacksolve
 * routines perform the actual solution of a dense linear system.
 *
 * Routines that work with realtype** 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 in smalldense.h
 * -----------------------------------------------------------------
 */

#ifndef _DENSE_H
#define _DENSE_H

#ifdef __cplusplus  /* wrapper to enable C++ usage */
extern "C" {
#endif

#include "sundialstypes.h"
#include "smalldense.h"

/*
 * -----------------------------------------------------------------
 * Type : DenseMat
 * -----------------------------------------------------------------
 * The type DenseMat is defined to be a pointer to a structure
 * with a size and a data field. The size field indicates the
 * number of columns (== number of rows) of a dense matrix, while
 * the data field is a two dimensional array used for component
 * storage. The elements of a dense matrix are stored columnwise
 * (i.e columns are stored one on top of the other in memory). If
 * A is of type DenseMat, then the (i,j)th element of A (with
 * 0 <= i,j <= size-1) is given by the expression (A->data)[j][i]
 * or by the expression (A->data)[0][j*n+i]. The macros below
 * allow a user to access efficiently individual matrix
 * elements without writing out explicit data structure
 * references and without knowing too much about the underlying
 * element storage. The only storage assumption needed is that
 * elements are stored columnwise and that a pointer to the jth
 * column of elements can be obtained via the DENSE_COL macro.
 * Users should use these macros whenever possible.
 * -----------------------------------------------------------------
 */

typedef struct _DenseMat {
  long int size;
  realtype **data;
} *DenseMat;

/* DenseMat accessor macros */

/*
 * -----------------------------------------------------------------
 * Macro : DENSE_ELEM
 * -----------------------------------------------------------------
 * Usage : DENSE_ELEM(A,i,j) = a_ij;  OR
 *         a_ij = DENSE_ELEM(A,i,j);
 * -----------------------------------------------------------------
 * DENSE_ELEM(A,i,j) references the (i,j)th element of the N by N
 * DenseMat A, 0 <= i,j <= N-1.
 * -----------------------------------------------------------------
 */

#define DENSE_ELEM(A,i,j) ((A->data)[j][i])

/*
 * -----------------------------------------------------------------
 * Macro : DENSE_COL
 * -----------------------------------------------------------------
 * Usage : col_j = DENSE_COL(A,j);
 * -----------------------------------------------------------------
 * DENSE_COL(A,j) references the jth column of the N by N
 * DenseMat A, 0 <= j <= N-1. The type of the expression
 * DENSE_COL(A,j) is realtype *. After the assignment in the usage
 * above, col_j may be treated as an array indexed from 0 to N-1.
 * The (i,j)th element of A is referenced by col_j[i].
 * -----------------------------------------------------------------
 */

#define DENSE_COL(A,j) ((A->data)[j])

/* Functions that use the DenseMat representation for a dense matrix */

/*
 * -----------------------------------------------------------------
 * Function : DenseAllocMat
 * -----------------------------------------------------------------
 * Usage : A = DenseAllocMat(N);
 *         if (A == NULL) ... memory request failed
 * -----------------------------------------------------------------
 * DenseAllocMat allocates memory for an N by N dense matrix and
 * returns the storage allocated (type DenseMat). DenseAllocMat
 * returns NULL if the request for matrix storage cannot be
 * satisfied. See the above documentation for the type DenseMat
 * for matrix storage details.
 * -----------------------------------------------------------------
 */

DenseMat DenseAllocMat(long int N);

/*
 * -----------------------------------------------------------------
 * Function : DenseAllocPiv
 * -----------------------------------------------------------------
 * Usage : p = DenseAllocPiv(N);
 *         if (p == NULL) ... memory request failed
 * -----------------------------------------------------------------
 * DenseAllocPiv allocates memory for pivot information to be
 * filled in by the DenseFactor routine during the factorization
 * of an N by N dense matrix. The underlying type for pivot
 * information is an array of N integers and this routine returns
 * the pointer to the memory it allocates. If the request for
 * pivot storage cannot be satisfied, DenseAllocPiv returns NULL.
 * -----------------------------------------------------------------
 */

long int *DenseAllocPiv(long int N);

/*
 * -----------------------------------------------------------------
 * Function : DenseFactor
 * -----------------------------------------------------------------
 * Usage : ier = DenseFactor(A, p);
 *         if (ier != 0) ... A is singular
 * -----------------------------------------------------------------
 * DenseFactor performs the LU factorization of the N by N dense
 * matrix A. This is done using standard Gaussian elimination
 * with partial pivoting.
 *
 * 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.
 *
 * DenseFactor 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.
 * -----------------------------------------------------------------
 */

long int DenseFactor(DenseMat A, long int *p);

/*
 * -----------------------------------------------------------------
 * Function : DenseBacksolve
 * -----------------------------------------------------------------
 * Usage : DenseBacksolve(A, p, b);
 * -----------------------------------------------------------------
 * DenseBacksolve solves the N-dimensional system A x = b using
 * the LU factorization in A and the pivot information in p
 * computed in DenseFactor. The solution x is returned in b. This
 * routine cannot fail if the corresponding call to DenseFactor
 * did not fail.
 * -----------------------------------------------------------------
 */

void DenseBacksolve(DenseMat A, long int *p, realtype *b);

/*
 * -----------------------------------------------------------------
 * Function : DenseZero
 * -----------------------------------------------------------------
 * Usage : DenseZero(A);
 * -----------------------------------------------------------------
 * DenseZero sets all the elements of the N by N matrix A to 0.0.
 * -----------------------------------------------------------------
 */

void DenseZero(DenseMat A);

/*
 * -----------------------------------------------------------------
 * Function : DenseCopy
 * -----------------------------------------------------------------
 * Usage : DenseCopy(A, B);
 * -----------------------------------------------------------------
 * DenseCopy copies the contents of the N by N matrix A into the
 * N by N matrix B.
 * -----------------------------------------------------------------
 */

void DenseCopy(DenseMat A, DenseMat B);

/*
 * -----------------------------------------------------------------
 * Function: DenseScale
 * -----------------------------------------------------------------
 * Usage : DenseScale(c, A);
 * -----------------------------------------------------------------
 * DenseScale scales the elements of the N by N matrix A by the
 * constant c and stores the result back in A.
 * -----------------------------------------------------------------
 */

void DenseScale(realtype c, DenseMat A);

/*
 * -----------------------------------------------------------------
 * Function : DenseAddI
 * -----------------------------------------------------------------
 * Usage : DenseAddI(A);
 * -----------------------------------------------------------------
 * DenseAddI adds the identity matrix to A and stores the result
 * back in A.
 * -----------------------------------------------------------------
 */

void DenseAddI(DenseMat A);

/*
 * -----------------------------------------------------------------
 * Function : DenseFreeMat
 * -----------------------------------------------------------------
 * Usage : DenseFreeMat(A);
 * -----------------------------------------------------------------
 * DenseFreeMat frees the memory allocated by DenseAllocMat for
 * the N by N matrix A.
 * -----------------------------------------------------------------
 */

void DenseFreeMat(DenseMat A);

/*
 * -----------------------------------------------------------------
 * Function : DenseFreePiv
 * -----------------------------------------------------------------
 * Usage : DenseFreePiv(p);
 * -----------------------------------------------------------------
 * DenseFreePiv frees the memory allocated by DenseAllocPiv for
 * the pivot information array p.
 * -----------------------------------------------------------------
 */

void DenseFreePiv(long int *p);

/*
 * -----------------------------------------------------------------
 * Function : DensePrint
 * -----------------------------------------------------------------
 * Usage : DensePrint(A);
 * -----------------------------------------------------------------
 * This routine 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/lg/Lg option. A blank line is printed
 * before and after the matrix.
 * -----------------------------------------------------------------
 */

void DensePrint(DenseMat A);

#ifdef __cplusplus
}
#endif

#endif
