/* This is -*- C -*- */
/* vim: set sw=2: */
/* $Id: guppi-matrix.h,v 1.1 2001/04/14 01:19:21 trow Exp $ */

/*
 * guppi-matrix.h
 *
 * Copyright (C) 2001 The Free Software Foundation, Inc.
 *
 * Developed by Jon Trowbridge <trow@gnu.org>
 */

/*
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation; either version 2 of the
 * License, or (at your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 */

/* 
   WARNING!   WARNING!   WARNING!   WARNING!   WARNING!   WARNING!

   This implementation is not suitable for any sort of heavy-duty
   (or even medium-duty) linear algebra computations.  It is only
   a convenience to be used by Guppi internally, and shouldn't be
   directly exposed to the user in any way.

   If you want to do real linear algebra calculations, use a real
   linear algebra package --- not this hacked-up mess.

   WARNING!   WARNING!   WARNING!   WARNING!   WARNING!   WARNING!
*/


#ifndef __GUPPI_MATRIX_H__
#define __GUPPI_MATRIX_H__

#include <glib.h>
#include "guppi-vector.h"

typedef struct _GuppiMatrix GuppiMatrix;
struct _GuppiMatrix {
  gint r, c;
  double *data;
  double epsilon;

  GuppiMatrix *LU;
  gint *perm;
};

#define guppi_matrix_rows(m) ((m)->r)
#define guppi_matrix_cols(m) ((m)->c)
#define guppi_matrix_is_square(m) ((m)->c == (m)->r)

GuppiMatrix *guppi_matrix_new  (gint r, gint c);
GuppiMatrix *guppi_matrix_copy (GuppiMatrix *);
void         guppi_matrix_free (GuppiMatrix *);

/* If you've directly manipulated the GuppiMatrix entries, you
   need to call guppi_matrix_touch () to make sure that 
   cached data (like the LU decomposition) is reset. */
void         guppi_matrix_touch (GuppiMatrix *);

void         guppi_matrix_set_constant (GuppiMatrix *, double);

#define guppi_matrix_ptr(m, i, j)       ((m)->data + ((i)*((m)->c)) + (j))
#define guppi_matrix_ptr_col_incr(m, p) ((p)+1)
#define guppi_matrix_ptr_row_incr(m, p) ((p)+((m)->c))
#define guppi_matrix_entry(m, i, j)       (*(guppi_matrix_ptr (m,i,j)))

GuppiVector *guppi_matrix_get_row (GuppiMatrix *, gint r);
GuppiVector *guppi_matrix_get_col (GuppiMatrix *, gint c);

void     guppi_matrix_normalize_row  (GuppiMatrix *, gint r);

gboolean guppi_matrix_row_is_nonzero (GuppiMatrix *, gint r);
gboolean guppi_matrix_column_is_nonzero (GuppiMatrix *, gint r);

double   guppi_matrix_row_dot (GuppiMatrix *, gint r1, gint r2);

void guppi_matrix_subtract_scaled_row_from_row (GuppiMatrix *m, double scale,
						gint r, gint r_sub_from);

GuppiVector *guppi_matrix_solve               (GuppiMatrix *, GuppiVector *vec);
GuppiVector *guppi_matrix_solve_with_fallback (GuppiMatrix *, GuppiVector *vec,
					       gboolean (*fallback)(GuppiMatrix *, GuppiVector *vec, gint i, gpointer),
					       gpointer);

GuppiMatrix *guppi_matrix_invert (GuppiMatrix *);
GuppiVector *guppi_matrix_apply  (GuppiMatrix *, GuppiVector *);

void guppi_matrix_spew (GuppiMatrix *);

#endif /* __GUPPI_MATRIX_H__ */