/* Vladimir Kotlyar Sepember 1995
 *
 * Matrix manipulation
 *
 */

#ifndef _a_matrix_h
#define _a_matrix_h

#include "alpha_matrix_private.h"

typedef xa_vec_t *	AlphaVector;
typedef xa_mx_t *	AlphaMatrix;


#define a_size_vec(v)		xa_size_vec(v)
#define a_length_vec(v)		xa_size_vec(v)
#define a_elt_vec(v,i)		xa_elt_vec(v,i)

#define a_rows_mx(A)		xa_rows_mx(A)
#define a_cols_mx(A)		xa_cols_mx(A)
#define a_col_mx(A,j)		xa_col_mx((A),(j))
#define a_elt_mx(A,i,j)		xa_elt_mx((A),(i),(j))

AlphaVector	a_make_vector (arena_info *arena, int n);
void		a_set_vec (AlphaVector v, int x);
void		a_neg_vec (AlphaVector x);
void		a_saxpy (int a, AlphaVector x,
			 int b, AlphaVector y,
			 AlphaVector z);	/* z := a x + b y */
void		a_copy_vec (AlphaVector dest, AlphaVector src);
void		a_swap_vec (AlphaVector *px, AlphaVector *py);
void		a_resize_vec (arena_info *arena, AlphaVector *px, int newn);
int		a_is_zero_vec (AlphaVector x);
int		a_gcd_vec (AlphaVector x);
void		a_extract_vec (AlphaVector x, int i1, int i2, AlphaVector y);

void		a_read_vector (FILE *fp, AlphaVector x);
void		a_dump_vector (FILE *fp, AlphaVector x);

AlphaMatrix	a_make_matrix (arena_info *arena, int m, int n);
void		a_copy_mx (AlphaMatrix dst, AlphaMatrix src);
void		a_swap_col (AlphaMatrix A, int i, int j);
void		a_swap_row (AlphaMatrix A, int i, int j);
void		a_set2ident (AlphaMatrix A);
int		a_is_zero_matrix (AlphaMatrix A);
int		a_is_ident_matrix (AlphaMatrix A);
void		a_matmult (arena_info *arena, AlphaMatrix A,  AlphaMatrix B, 
			   AlphaMatrix C); 				/* C = A * B */
void		a_mvm (arena_info *arena, AlphaMatrix A, AlphaVector x, 
		       AlphaVector y); /* y = A x */
void		a_transpose (AlphaMatrix A, AlphaMatrix At); /* At = A^T */

void		a_extract (AlphaMatrix A, 
			   int i1, int i2, int j1, int j2, 
			   AlphaMatrix B);  /* B = A[i1:i2, j1:j2] */
void		a_replace (AlphaMatrix A, 
			   int i1, int i2, int j1, int j2, 
			   AlphaMatrix B);  /* A[i1:i2, j1:j2]  = B*/
void		a_add_cols (arena_info *arena, AlphaMatrix *pA, int c);
void		a_trunc_cols (arena_info *arena, AlphaMatrix *pA, int c);

void		a_lower_echelon (arena_info *arena, AlphaMatrix A, AlphaMatrix U, int *prank);
void		a_snf (arena_info *arena, AlphaMatrix A, AlphaMatrix U, AlphaMatrix V, int *prank);
int		a_solve_linear (arena_info *arena, AlphaMatrix A, AlphaVector B, 
				AlphaVector X0, AlphaMatrix H, int *pnull_dim);
int		a_solve_echelon (AlphaMatrix L, AlphaVector B, AlphaVector X);

void            a_unit_backsolve(AlphaMatrix U, AlphaVector b, AlphaVector x);
                    /* x = U^{-1} * b, U is unit-upper-tri */

void            a_block_unit_backsolve(AlphaMatrix U, AlphaMatrix B, AlphaMatrix X);
                    /* X = U^{-1} * B, U is unit-upper-tri */


/* A matrix is read BY ROW */
void		a_read_matrix (FILE *fp, AlphaMatrix A);
void		a_dump_matrix (FILE *fp, AlphaMatrix A);
void		a_print_flat_matrix (FILE *fp, AlphaMatrix A);

#define A_FOR_COLS_MX(A,col,j)  XA_FOR_COLS_MX((A),col,j)
#define A_END_COLS_MX           XA_END_COLS_MX

#endif


