
/** @file  psSparse.h
 *
 * functions to manipulate sparse matrices equations
 *  
 */

# ifndef PS_SPARSE_H
# define PS_SPARSE_H

// constraints to limit the range of the matrix equation solution
typedef struct
{
    double paramDelta;
    double paramMin;
    double paramMax;
}
psSparseConstraint;

// a sparse matrix equation: A x = Bf
// Aij contains the populated elements of the matrix
// Bfj contains the elements of the vector Bf
// Qii contains the diagonal elements of Aij
// Si contains the i-index values of Aij
// Sj contains the j-index values of Aij
typedef struct
{
    psVector *Aij;
    psVector *Bfj;
    psVector *Qii;
    psVector *Si;
    psVector *Sj;
    int Nelem;
    int Nrows;
}
psSparse;

// allocate a sparse matrix structure
psSparse *psSparseAlloc (int Nrows, int Nelem);

// add a new matrix element
// user should only add elements above the diagonal
void psSparseMatrixElement (psSparse *sparse, int i, int j, float value);

// define a new sparse matrix equation vector element
void psSparseVectorElement (psSparse *sparse, int i, float value);

// perform the operation matrix * vector on a sparse matrix and a vector
psVector *psSparseMatrixTimesVector (psVector *output, psSparse *matrix, psVector *vector);

// re-sort a sparse matrix to have all elements in index order rather than insertion order
// call this before solving, but after populating matrix and vector
void psSparseResort (psSparse *sparse);

// solve the equation A x = Bf for the value of x
// a good starting guess is the vector Bf
psVector *psSparseSolve (psVector *guess, psSparseConstraint constraint, psSparse *sparse, int Niter);

// test of the sparse matrix solutions (move to test suite)
void psSparseMatrixTest ();

# endif /* PS_SPARSE_H */
