/** @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
typedef struct
{
    psVector *Aij;                      // Aij contains the populated elements of the matrix
    psVector *Bfj;                      // Bfj contains the elements of the vector Bf
    psVector *Qii;                      // Qii contains the diagonal elements of Aij
    psVector *Si;                       // Si contains the i-index values of Aij
    psVector *Sj;                       // Sj contains the j-index values of Aij
    int Nelem;                         // Number of elements
    int Nrows;                         // Number of rows
}
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
bool psSparseMatrixElement(psSparse *sparse, // Matrix to which to add
                           int i, int j, // Matrix indices at which to add
                           float value  // Value to add
                          );

// define a new sparse matrix equation vector element
void psSparseVectorElement(psSparse *sparse, // Matrix to which to add
                           int i,      // Index to add
                           float value  // Value to add
                          );

// perform the operation matrix * vector on a sparse matrix and a vector
psVector *psSparseMatrixTimesVector(psVector *output, // Output vector, or NULL
                                    const psSparse *matrix, // Sparse matrix
                                    const psVector *vector // Corresponding 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
bool psSparseResort(psSparse *sparse    // Matrix to re-sort
                   );

// solve the equation A x = Bf for the value of x
// a good starting guess is the vector Bf
psVector *psSparseSolve(psVector *output,// The output vector, or NULL
                        psSparseConstraint constraint, // Constraint to limit the range of the solution
                        const psSparse *sparse, // Sparse matrix
                        int Niter       // Number of iterations
                       );

#endif /* PS_SPARSE_H */
