// Copyright: Jon Wilkening, 1999.

#ifndef __SYM_MAT_H
#define __SYM_MAT_H

#ifdef __cplusplus
extern "C" {
#endif

#ifndef HIDE_FROM_DOXYGEN

typedef struct sym_mat_elt_struct {
    int row;
    double val;
} sym_mat_elt;

typedef struct sym_mat_list_elt_struct {
    int next;
    int row;
    double val;
} sym_mat_list_elt;

typedef struct sym_mat_hit_list_elt_struct {
    int next;
    int prev;
    double val;
} sym_mat_hit_list_elt;

typedef struct sym_mat_elt_list_struct {
    int n;
    int size, maxsize;
    sym_mat_list_elt *space;
} sym_mat_elt_list;

typedef struct sym_mat_row_elt_struct {
    int next;
    int idx;
} sym_mat_row_elt;

#endif // HIDE_FROM_DOXYGEN

typedef struct sym_mat_struct {
    int n;
    double flops;    // flops to factor A (not counting equil, solve, etc.)
    int nnzA;        // diagonal elements counted even if they are zero
    int nnzL;        // (but they will be nonzero if A is SPD)
    double scond;    // smallest eq / largest eq
    
    int no_more_entries;
    int already_factored;
    int b_is_set_up;

    sym_mat_elt_list v;
    sym_mat_elt *col;              // size n
    sym_mat_hit_list_elt *hits;    // size n
    sym_mat_row_elt *row;          // size n
    int *pa;                       // size n+1
    int *a;                        // size 1 or nnz_A
    int *pb;                       // size n+1
    sym_mat_elt *b;                // size 1 or nnz_A
    int *pc;                       // size n+1
    sym_mat_elt *c;                // size 1 or nnz_L
    int *perm1;                    // size n+1
    int *perm2;                    // size n
    double *eq;                    // size n
} sym_mat;

// allocate space, initial set-up
void sym_mat_init(sym_mat *A, int n);

// free all space allocated by a sym_mat routine
void sym_mat_free(sym_mat *A);

// add val to entry A(i,j).  indices are 0 based.
// Only lower triangle is accessed, so nothing is done if i<j
void sym_mat_add_entry(sym_mat *A, int i, int j, double val);

// return A(i,j) (original indices)
double sym_mat_get_entry(sym_mat *A, int i, int j);

// equilibrate if equil[0] != 'n' or 'N"
// re-order rows and columns if re-order != 'n' or 'N'
// do sparse cholesky decomposition of matrix
void sym_mat_factor(sym_mat *A, const char *equil, const char *re_order);

// replace b with A*x. (in-situ OK)
void sym_mat_apply(sym_mat *A, double *x, double *b);

// dumps lower triangle of B=PAP' (equilibration is reversed if necessary)
// dumps permutation P to get A back (in matlab: A = B(P,P))
// adds one to all indices so matlab can read the files
void sym_mat_dump_A(sym_mat *A,
		    const char *filename_A, const char *filename_perm,
		    const char *filename_matlab, int precision);

// replace x by solution to Ax = b. (in-situ OK)
void sym_mat_solve(sym_mat *A, double *b, double *x);
void sym_mat_apply_invL(sym_mat *A, double *b, double *x);

// compute LAPACK condition number estimate of (equilibrated) A
// (based on LAPACK's  dpocon.f  and  dlansy.f)
double sym_mat_rcond(sym_mat *A);

// do iterative refinement and compute forward and backward error estimates.
// (based on LAPACK's  dporfs.f)
void sym_mat_refine(sym_mat *A, double *b, double *x,
		    double *ferr, double *berr);

// compute largest and smallest eigenvalue
// (mode = 1: largest eigenvalue of A, mode = 2: largest eigenvalue of A^{-1}
void sym_mat_arpack(sym_mat *A, int mode);

#ifdef __cplusplus
} // end extern "C" block
#endif

#endif