#include <stdbool.h>
#include <stddef.h>

#ifndef CA_UWATERLOO_DWHARDER_LINEAR_ALGEBRA
#define CA_UWATERLOO_DWHARDER_LINEAR_ALGEBRA

typedef enum {
    ONE_NORM,
    MEAN_ERROR,
    TWO_NORM,
    TWO_NORM_SQUARED,
    ROOT_MEAN_SQUARED_ERROR,
    MEAN_SQUARED_ERROR,
    INFINITY_NORM
} norm_t;

typedef enum {
    MAPLE,
    MATLAB,
    C,
    MATHEMATICA
} output_t;

typedef enum {
    UNIFORM_POSITIVE,
    UNIFORM_SYMMETRIC,
    STANDARD_NORMAL
} distribution_t;

typedef enum {
    RANDOM_DIAGONAL,
    DENSE_DIAGONAL,
    ZERO_DIAGONAL
} diagonal_t;

typedef enum {
    SYMMETRIC,
    SKEW_SYMMETRIC,
    UPPER_TRIANGULAR,
    LOWER_TRIANGULAR,
    GENERIC
} matrix_shape_t;

/**********************************************************
 * Vector structure definitions and function declarations *
 **********************************************************/

typedef struct {
    size_t dimension;
    double *a_values;
} vector_t;

// Vector initialization and clean up
  void vector_init( vector_t *, size_t );
  void vector_init_scalar( vector_t *, size_t, double );
  void vector_init_vector( vector_t *, vector_t * );
  void vector_destroy( vector_t * );

// Vector set, get, assignment and queries
double vector_get( vector_t *, size_t );
void vector_set( vector_t *, size_t, double );

  void vector_vector_assign( vector_t *, vector_t * );
  void vector_scalar_assign( vector_t *, double );

size_t vector_dimension( vector_t * );

// Vector operations
  void vector_scalar_multiply( vector_t *, double );
  void vector_scalar_add( vector_t *, double );
  void vector_vector_add( vector_t *, vector_t * );
  void vector_vector_scalar_add( vector_t *, vector_t *, double );

// Vector functions
double inner_product( vector_t *, vector_t * );
  void vector_project( vector_t *, vector_t * );
  void vector_perpendicularize( vector_t *, vector_t * );
double vector_norm( vector_t *, norm_t );
  void vector_normalize( vector_t *, norm_t );
double vector_distance( vector_t *, vector_t *, norm_t );

  void gram_schmidt( vector_t **, size_t );

// Random vector generator

void random_vector( vector_t *, distribution_t );

// Vector printing
   int vector_print_language( char *, char *, char *, output_t );
   int vector_printf( vector_t *, char *, output_t );

/*****************************************************************
 * Sparse matrix structure definitions and function declarations *
 *****************************************************************/

typedef struct {
    size_t dimension;
    size_t max_off_diagonal_values;
    size_t *a_row_indices;
    size_t *a_column_indices;
    double *a_diagonal_values;
    double *a_off_diagonal_values;
} sparse_matrix_t;

// Matrix initialization and clean up
  void matrix_init( sparse_matrix_t *, size_t, size_t );
  void matrix_init_scalar( sparse_matrix_t *, size_t, size_t, double );
  void matrix_init_vector( sparse_matrix_t *, vector_t *, size_t );
  void matrix_destroy( sparse_matrix_t * );

// Matrix set, get, append, assignment and queries
double matrix_get( sparse_matrix_t *, size_t, size_t );
  void matrix_set( sparse_matrix_t *, size_t, size_t, double );
  void matrix_append( sparse_matrix_t *, size_t, size_t, double );

  void matrix_matrix_assign( sparse_matrix_t *, sparse_matrix_t * );
  void matrix_vector_assign( sparse_matrix_t *, vector_t * );
  void matrix_scalar_assign( sparse_matrix_t *, double );

size_t matrix_dimension( sparse_matrix_t * );
double matrix_density( sparse_matrix_t * );
  bool is_diagonal( sparse_matrix_t * );
  bool is_upper_triangular( sparse_matrix_t * );
  bool is_lower_triangular( sparse_matrix_t * );
  bool is_symmetric( sparse_matrix_t * );
  bool is_skew_symmetric( sparse_matrix_t * );
  bool is_valid( sparse_matrix_t * );

// Matrix operations
  void matrix_scalar_multiply( sparse_matrix_t *, double );
  void matrix_matrix_add( sparse_matrix_t *, sparse_matrix_t *, sparse_matrix_t * );
  void matrix_vector_add( sparse_matrix_t *, vector_t * );
  void matrix_scalar_add( sparse_matrix_t *, double );

// Matrix printing
   int matrix_data_structure_printf( sparse_matrix_t *, char * );
   int matrix_density_print( sparse_matrix_t * );
   int matrix_printf( sparse_matrix_t *, char *, output_t );

// Matrix functions and solvers
  void matrix_vector_multiply( vector_t *, sparse_matrix_t *, vector_t * );
  void matrix_transpose_vector_multiply( vector_t *, sparse_matrix_t *, vector_t * );

  void      diagonal_solve( vector_t *, sparse_matrix_t *, vector_t * );
  void backward_substitute( vector_t *, sparse_matrix_t *, vector_t * );
  void  forward_substitute( vector_t *, sparse_matrix_t *, vector_t * );

  size_t                         jacobi( vector_t *, sparse_matrix_t *, vector_t *, vector_t *, size_t, double );
  size_t                   gauss_seidel( vector_t *, sparse_matrix_t *, vector_t *, vector_t *, size_t, double );
  size_t               steepest_descent( vector_t *, sparse_matrix_t *, vector_t *,
                                         vector_t *, vector_t *, size_t, double );
  size_t               minimal_residual( vector_t *, sparse_matrix_t *, vector_t *,
                                         vector_t *, vector_t *, size_t, double );
  size_t residual_norm_steepest_descent( vector_t *, sparse_matrix_t *, vector_t *,
                                         vector_t *, vector_t *, vector_t *, size_t, double );

  void lu( sparse_matrix_t *, sparse_matrix_t *, sparse_matrix_t * );

// Random sparse matrix generator
void               random_matrix( sparse_matrix_t *, matrix_shape_t, distribution_t, diagonal_t, double );
void      random_diagonal_matrix( sparse_matrix_t *, distribution_t, diagonal_t, double );
void   random_tridiagonal_matrix( sparse_matrix_t *, matrix_shape_t, distribution_t );
void random_band_diagonal_matrix( sparse_matrix_t *, size_t, matrix_shape_t, distribution_t, diagonal_t, double );
#endif