Geografic-Information-System/include/la.h

/**
 * \file la.h
 *
 * \brief Wrapper headers for BLAS/LAPACK.
 *
 * This program is free software under the GNU General Public License
 * (>=v2). Read the file COPYING that comes with GRASS for details.
 *
 * \author David D. Gray <ddgray AT armacde demon co uk>
 * \author GRASS GIS Development Team
 *
 * \date 2000-2007
 */

#ifndef LA_H_
#define LA_H_

/* QUESTION: On some systems there appears to be no default link
   to this. Do we create a symlink to
   /usr/lib/gcc-lib/<platform>/<vers_num>/include/g2c.h

   or link to any old f2c.h that happens to hanging around?

   A job for autoconf

   [Also a consideration for -lg2c]
*/

#include <grass/config.h>
#include <stdio.h>

#ifdef HAVE_G2C_H
#include <g2c.h>
#else /* for gcc4+ */
typedef int                       integer;
typedef unsigned int              uinteger;
typedef char *                    address;
typedef short                     shortint;
typedef float                     real;
typedef double                    doublereal;
typedef struct {real r, i;}       complex;
typedef struct {doublereal r, i;} doublecomplex;
typedef int                       logical;
typedef short                     shortlogical;
typedef char                      logical1;
typedef char                      integer1;
typedef long                      longint;
typedef unsigned long             ulongint;

/* IO stuff */
typedef int ftnlen;

/* procedure parameter types for -A */
typedef int          (*U_fp)();
typedef shortint     (*J_fp)();
typedef integer      (*I_fp)();
typedef real         (*R_fp)();
typedef doublereal   (*D_fp)(), (*E_fp)();
typedef void         (*C_fp)();
typedef void         (*Z_fp)();
typedef logical      (*L_fp)();
typedef shortlogical (*K_fp)();
typedef void         (*H_fp)();
typedef int          (*S_fp)();

/* E_fp is for real functions when -R is not specified */
typedef void       C_f;   /* complex function                    */
typedef void       H_f;   /* character function                  */
typedef void       Z_f;   /* double complex function             */
typedef doublereal E_f;   /* real function with -R not specified */
#endif /* HAVE_G2C_H */

/* The following may have to be selectively installed according
   to platform, at least partly
*/

#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
#include <grass/blas.h>
#include <grass/lapack.h>
#endif


/* Useful defines */

#define MAX_POS          1   /* Indicates maximum value         */
#define MAX_NEG         -1   /* Indicates minimum value         */
#define MAX_ABS          0   /* Indicates absolute value        */

#define DO_COMPACT       0   /* Elliminate unnecessary rows (cols) in matrix  */
#define NO_COMPACT       1   /* ... or not                                    */


/* define macros for fortran symbols (called directly). Needed because 
   of platform invariance on fortran->C symbol translations
*/

#if defined(HAVE_LIBBLAS) && defined(HAVE_LIBLAPACK)
#define f77_dgesv                   dgesv_
#define f77_dgemm                   dgemm_
#define f77_dnrm2                   dnrm2_
#endif

/* Operations should know type of coefficient matrix, so that
   they can call the right driver
*/

typedef enum { NONSYM, SYM, HERMITIAN } mat_type;
typedef enum { MATRIX_, ROWVEC_, COLVEC_ } mat_spec;
typedef enum { RVEC, CVEC }             vtype;


/************************************************************
 *                                                          *
 * A general matrix wrapper for use with BLAS / LAPACK      *
 *  routines, and perhaps elsewhere                         *
 *                                                          *
 ************************************************************/

typedef struct matrix_ {

  mat_spec type;    /* matrix, row vector or column vector? */
  int v_indx;       /* If a vector, which row(column) is active?
                     * If a matrix this is ignored. If value is < 0,
                     * the first row(column) is assumed, ie. index 0.  */
  int rows, cols;   /* Rows and columns of matrix */
  int ldim;         /* Lead dimension of matrix. How many `rows' are
                     * alloc'ed? May exceed real number of rows `rows' */
  doublereal *vals; /* The values (should be dimensioned to lda * cols */
  int is_init;      /* Is  matrix initialised: values array
                     * is allocated and parameters set ?               */
}  mat_struct;

typedef mat_struct vec_struct;


/* Prototypes */

/* Matrix routines corresponding to BLAS Level III */

mat_struct *G_matrix_init(int, int, int);
int G_matrix_zero(mat_struct *);
int G_matrix_set(mat_struct *, int, int, int);
mat_struct *G_matrix_copy(const mat_struct *);
mat_struct *G_matrix_add(mat_struct *, mat_struct *);
mat_struct *G_matrix_subtract(mat_struct *, mat_struct *);
mat_struct *G_matrix_scale(mat_struct *, const double);
mat_struct *G__matrix_add(mat_struct *, mat_struct *, const double, const double);
mat_struct *G_matrix_product(mat_struct *, mat_struct *);
mat_struct *G_matrix_transpose(mat_struct *);
int G_matrix_LU_solve(const mat_struct *, mat_struct **, const mat_struct *, mat_type);
mat_struct *G_matrix_inverse(mat_struct *);
void G_matrix_free(mat_struct *);
void G_matrix_print(mat_struct *);
int G_matrix_set_element(mat_struct *, int, int, double);
double G_matrix_get_element(mat_struct *, int, int);


/* Matrix-vector routines corresponding to BLAS Level II */

vec_struct *G_matvect_get_column(mat_struct *, int);
vec_struct *G_matvect_get_row(mat_struct *, int);
int G_matvect_extract_vector( mat_struct *, vtype, int);
int G_matvect_retrieve_matrix(vec_struct *);


/* Vector routines corresponding to BLAS Level I */

vec_struct *G_vector_init(int, int, vtype);
int G_vector_set(vec_struct *, int, int, vtype, int);
double G_vector_norm_euclid(vec_struct *);
double G_vector_norm_maxval(vec_struct *, int);
vec_struct *G_vector_copy(const vec_struct *, int);

/* Matrix and vector routines corresponding to ?? */

void G_vector_free(vec_struct *);
vec_struct *G_vector_sub(vec_struct *, vec_struct *, vec_struct *);
double G_vector_norm1(vec_struct *);
int G_matrix_read(FILE *, mat_struct *);
int G_matrix_stdin(mat_struct *);
int G_matrix_eigen_sort(vec_struct *, mat_struct *);

#endif /* LA_H_ */