/*
  Module name: lmat.h
  Created by: Ljubomir Buturovic
  Created: 04/19/2002
  Purpose: declarations for lmat.c.
*/

/*
  Copyright 2004 Ljubomir J. Buturovic

  Permission is hereby granted, free of charge, to any person
  obtaining a copy of this software and associated documentation files
  (the "Software"), to deal in the Software without restriction,
  including without limitation the rights to use, copy, modify, merge,
  publish, distribute, sublicense, and/or sell copies of the Software,
  and to permit persons to whom the Software is furnished to do so,
  subject to the following conditions:

  The above copyright notice and this permission notice shall be
  included in all copies or substantial portions of the Software.

  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/

/*
  Estimate mean given 'nvec' 'd'-dimensional input vectors in 'x'. If
  'x' is (float **) 0, or 'nvec' < 1, or 'd' < 1, or in case of
  malloc() error, returns (float **) 0, and sets errno in the latter
  case.
*/
float *mest(float **x, int d, int nvec);

/*
  Constants for cest().
*/

#define COVARIANCE_MATRIX  0
#define CORRELATION_MATRIX 1

/*
  Estimate covariance ('imx' == COVARIANCE_MATRIX) or correlation
  matrix ('imx' == CORRELATION_MATRIX) for 'nvec' 'd'-dimensional
  input vectors in 'x'. If 'x' is (float **) 0, or 'nvec' <= 0, or 'd'
  <= 0, or in case of malloc() error, returns (float **) 0, and sets
  errno.
*/
float **cest(float **x, int d, int nvec, int imx);

/*
  Convert covariance matrix 'sigma' to matrix of linear correlation
  coefficients.
*/
int correst(float **sigma, int d);

/*
  Convert symmetric 'matrix' to Symmetric Storage Mode (SSM). Return
  (float *) 0 in case of malloc() error and set errno.
*/
float *fmx_ssm(float **matrix, int nmx);

/*
  Compute scalar product of vectors 'avec' and 'bvec' of length
  'n'. If 'errc' is not NULL, set 'errc' to -1 in case of illegal
  arguments, otherwise set it to 0.
*/
float fvec_dot(float *avec, float *bvec, int n, int *errc);

/*
  Compute scalar product of vectors 'avec' and 'bvec' of length
  'n'. If 'errc' is not NULL, set 'errc' to -1 in case of illegal
  arguments, otherwise set it to 0.
*/
double dvec_dot(double *avec, double *bvec, int n, int *errc);

/*
  Set elements of 'vec' to 'value'.
*/
void dvec_set(double *vec, int len, double value);

/*
  Copy double vector 'src' into 'dest'.
*/
void fvec_dcopy(float *dest, double *src, int len);

/*
  Copy float vector 'src' into 'dest'.
*/
void dvec_fcopy(double *dest, float *src, int len);

/*
  Compute difference between vectors 'avec' and 'bvec' of length 'n'
  and store it to 'outvec'. If 'errc' is not NULL, set 'errc' to -1 in
  case of illegal arguments, otherwise set it to 0.
*/
void fvec_diff(float *avec, float *bvec, float *outvec, int n, int *errc);

/*
  Multiply vector 'vector' of length 'd' and symmetric 'matrix'
  d*d. The matrix is stored in IMSL Symmetric Storage Mode, and thus
  the array 'matrix' has d(d+1)/2 elements.  Store result in vector
  'mx' (must be preallocated).
*/
void fvec_smx(float *vector, float *matrix, int d, float *mx);

/*
  Multiply vector 'vector' of length 'd' and symmetric 'matrix'
  d*d. The matrix is stored in IMSL Symmetric Storage Mode, and thus
  the array 'matrix' has d(d+1)/2 elements. Return the result (which
  is a vector of length 'd').
  
  In case of error, return NULL and set errno. The errors are EINVAL
  and malloc(() errors.
*/
float *fvec_ssm(float *vector, float *matrix, int d);

/*
  Return norm of 'vector' of length 'len'.
*/
float fvec_norm(float *vector, int len);

/*
  Return norm of double 'vector' of length 'len'.
*/
double dvec_norm(double *vector, int len);

/*
  Return normalized 'vector' of length 'len' (i.e., the returned
  vector has norm 1). Return (float *) 0 in case of malloc() error and
  set errno.
*/
float *fvec_normalize(float *vector, int len);

/*
  Return transpose of 'matrix'.
*/
float **fmx_transpose(float **matrix, int rows, int columns);

/*
  LAPACK declarations.
*/
void ssyev_(char *jobz, char *uplo, const int *n, const float *a,
	    const int *lda, const float *w, const float *work,
	    const int *lwork, int *info);

void ssysv_(char *uplo, const int *n, int *nrhs, const float *a, const int *lda,
	    const int *ipiv, const float *b, const int *ldb, const float *work,
	    const int *lwork, int *info);

void ssytrf_(char *uplo, const int *n, const float *a, const int *lda,
	     const int *ipiv, const float *work, const int *lwork, int *info);

void ssytri_(char *uplo, const int *n, const float *a, const int *lda,
	     const int *ipiv, const float *work, int *info);

void sgebal_(char *job, const int *n, const float *a, const int *lda, 
	     const int *ilo, const int *ihi, const float *scale, int *info);

void sgehrd_(const int *n, const int *ilo, const int *ihi, const float *a, 
	     const int *lda, const float *tau, const float *work, const int *lwork,
	     int *info);

void sorghr_(const int *n, const int *ilo, const int *ihi, const float *a, 
	     const int *lda, const float *tau, const float *work, const int *lwork,
	     int *info);

void shseqr_(char *job, char *compz, const int *n, const int *ilo, const int *ihi,
	     const float *h, const int *ldh, const float *wr, const float *wi, 
	     const float *z, const int *ldz, const float *work, const int *lwork, int *info);

void strevc_(char *side, char *howmny, char *select, const int *n, const float *t, 
	     const int *ldt, const float *vl, const int *ldvl, const float *vr, 
	     const int *ldvr, const int *mm, int *m, const float *work, int *info);

void sgeevx_(char *balanc, char *jobvl, char *jobvr, char *sense, const int *n, 
	     const float *a, const int *lda, const float *wr, const float *wi, 
	     const float *vl, const int *ldvl, const float *vr, const int *ldvr,
	     const int *ilo, const int *ihi, const float *scale, const float *abnrm,
	     const float *rconde, const float *rcondv, const float *work, 
	     const int *lwork, const int *iwork, int *info);

/*
  Calculate eigenvalues and eigenvectors of an input 'imx' by 'imx'
  real symmetric matrix 'amx' using LAPACK function ssyev_().  The
  eigenvalues are returned in 'evl' in descending order, and the
  corresponding eigenvectors in columns of 'evx'.  Column (**evx)[i]
  corresponds to eigenvalue (*evl)[i]. In other words, upon return,
  eigenvector corresponding to largest eigenvalue will be in
  evx[0][0], evx[1][0], ... evx[imx-1][0]. Eigenvector corresponding
  to next largest eigenvalue will be in evx[0][1], evx[1][1],
  ... evx[imx-1][1], etc.

  Return 0 for success, -1 for malloc() failure, and set 'errc'. The
  errors are malloc() errors, LERR_ITMAX or EINVAL. LERR_ITMAX means
  the LAPACK function ssyev_() failed to converge.  EINVAL means bad
  input arguments: for example, 'imx' <= 0, or aml/evl NULL.

  Note: eigenvectors returned as columns for compatibility. Use
  eigsn() to get eigenvectors in rows.
*/
int eigsn_column(float **amx, int imx, float **evl, float ***evx, int *errc);

/*
  Compute eigenvectors and eigenvalues of a real symmetric
  matrix. This is a row-based version of eigsn_column().

  TBD: this should actually be the real computing function;
  eigsn_column() should transpose the result. Right now, it is the
  other way round. Later, not critical. ljb, 02/13/2003.
*/
int eigsn(float **amx, int imx, float **evl, float ***evx, int *errc);

/*
  Calculate eigenvalues and right eigenvectors of an input 'imx' by
  'imx' real matrix 'fmx' using LAPACK function sgeevx_().  The
  eigenvalues are returned in 'evl' in descending order, and the
  corresponding eigenvectors in rows of 'evx'.

  Return 0 for success, -1 for malloc() failure, and set 'errc'. The
  errors are malloc() errors, EINVAL, or LERR_ITMAX if the LAPACK
  function fails to converge.
*/
int eigen(float **fmx, int imx, float **evl, float ***evx, int *errc);

/*
  Invert real symmetric matrix 'fmx' using LAPACK functions ssytrf_()
  and ssytri_(). If 'dsign' and 'dexp' are not (float *) 0, return
  also the determinant: 'dsign' is the sign of the determinant
  (+/-1.0), and 'dexp' is natural logarithm of the absolute value, so
  the final determinant is dsign*exp(dexp).  The input matrix 'fmx' is
  not modified (the function allocates space for the inverse).

  Return NULL for failure, and set 'errc'. The errors are malloc()
  errors, EINVAL for bad arguments and LERR_SINGULAR for singular
  matrix.
*/
float **fmx_inv(float **fmx, int imx, float *dsign, float *dexp, int *errc);

/*
  Calculate determinant of a real symmetric matrix 'fmx' with 'imx'
  rows using LAPACK functions...  'dsign' is sign of the determinant
  (+/-1.0), and 'dexp' is natural logarithm of the absolute value, so
  the final determinant is dsign*exp(dexp).

  Return 0 for success, -1 for failure, and set 'errc'. The errors are
  malloc() errors and EINVAL for bad arguments.
*/
int fmx_det(float **fmx, int imx, float *dsign, float *dexp, int *errc);

/*
  Multiply 'amatrix' with 'bmatrix' and transpose the
  result. 

  'amatrix' is 'arows' by 'acolumns'. 'bmatrix' must be 'acolumns' by
  'bdim'. The product is 'bdim' by 'arows'.

  If 'tflag' is 1, multiply 'amatrix' with transpose of 'bmatrix'. In
  this case 'bmatrix' must be 'bdim' by 'acolumns' (meaning that the
  transpose is 'acolumns' by 'bdim'). The result is again 'bdim' by
  'arows'.

  The function allocates space for the product.

  In case of error, return (float **) 0 and set errno. The only
  possible errors are memory allocation errors.
*/
float **fmx_tmult(float **amatrix, int arows, int acolumns, float **bmatrix, 
		  int bdim, int tflag);

/*
  Multiply 'amatrix' with 'bmatrix', transpose the result and store it
  in 'fmx'. 'amatrix' is 'arows' by 'acolumns'. 'bmatrix' must be
  'acolumns' by 'bdim'. The product is 'bdim' by 'arows'.

  If 'tflag' is 1, multiply 'amatrix' with transpose of 'bmatrix'. In
  this case 'bmatrix' must be 'bdim' by 'acolumns' (meaning that the
  transpose is 'acolumns' by 'bdim'). The result is again 'bdim' by
  'arows'.

  The function assumes that the caller has allocated sufficient space
  in 'bdim' by 'arows' matrix 'fmx', for example by:

  fmx = fmx_alloc(bdim, arows);

  In case of success, return 0. If 'amatrix', 'bmatrix' or 'fmx' are
  NULL, or dimensions of 'fmx' are incompatible with dimensions of
  'amatrix', 'bmatrix', the function returns -1.
*/
int fmx_tmulta(float **fmx, float **amatrix, int arows, int acolumns, 
	       float **bmatrix, int bdim, int tflag);

/*
  Multiply 'amatrix' with 'bmatrix'. 'amatrix' is 'arows' by
  'acolumns'. 'bmatrix' must be 'acolumns' by 'bdim'. The product is
  'arows' by 'bdim'.

  If 'tflag' is 1, multiply 'amatrix' with transpose of 'bmatrix'. In
  this case 'bmatrix' must be 'bdim' by 'acolumns' (meaning that the
  transpose is 'acolumns' by 'bdim'). The result is again 'arows' by
  'bdim'.

  The function allocates space for the product.

  In case of error, return (float **) 0 and set errno. The only
  possible errors are memory allocation errors.
*/
float **fmx_mult(float **amatrix, int arows, int acolumns, float **bmatrix, 
		 int bdim, int tflag);

/*
  Multiply 'amatrix' with 'bmatrix' and store result in
  'fmx'. 'amatrix' is 'arows' by 'acolumns'. 'bmatrix' must be
  'acolumns' by 'bdim'. The product is 'arows' by 'bdim'.

  If 'tflag' is 1, multiply 'amatrix' with transpose of 'bmatrix'. In
  this case 'bmatrix' must be 'bdim' by 'acolumns' (meaning that the
  transpose is 'acolumns' by 'bdim'). The result is again 'arows' by
  'bdim'.

  The function assumes that the caller has allocated sufficient space
  in 'arows' by 'bdim' matrix 'fmx', for example by:

  fmx = fmx_alloc(arows, bdim);

  In case of success, return 0. If 'amatrix', 'bmatrix' or 'fmx' are
  NULL, or dimensions of 'fmx' are incompatible with dimensions of
  'amatrix', 'bmatrix', the function returns -1.
*/
int fmx_multa(float **fmx, float **amatrix, int arows, int acolumns, 
	      float **bmatrix, int bdim, int tflag);

/*
  Return mean value of column 'col' in 'fmx'.
*/
float fmx_col_mean(float **fmx, int rows, int col);

/*
  Convert vector 'fvec' of length 'len' to 'rows' by 'columns' matrix.

  Return the matrix in case of success, otherwise return NULL and set
  errno. The errors are bad arguments and memory allocation errors.
*/
float **fmx_fvec(float *fvec, int len, int rows, int columns);

/*
  Return index of min. element in 'vec'. Return -1 in case of bad
  arguments.
*/
int fvec_argmin(float *fvec, int len);

/*
  Return index of max. element in 'vec'. Return -1 in case of bad
  arguments.
*/
int fvec_argmax(float *fvec, int len);

/*
  Return index, and optionally 'value', of max. element in 'vec'.
*/
int fvec_valmax(float *vec, int len, float *value);

/*
  Find max. element in 'matrix'. The indexes are 0-based.
*/
void fmx_max(float **matrix, int rows, int columns, int *max_row, int *max_column,
	     float *value);