/*
 *  Copyright 2002 by Soos, Antal 
 *  All rights reserved. Property of Soos, Antal. 
 *  Restricted rights to use, duplicate or disclose this code are 
 *  granted through contract.   
 */ 
/***************************************************************************/ 
/*                                                                         */ 
/*     MatrixCal . hpp                                                     */ 
/*                                                                         */ 
/*     Basic C++ standard matrix calculations.                             */ 
/*                                                                         */ 
/*                                                                         */ 
/***************************************************************************/ 
 
#ifndef _MATRIX_CAL_HPP_ 
#define _MATRIX_CAL_HPP_    
//----------------------------------------------------------------------------- 
//----------------------------------------------------------------------------- 
// 
//		The used include files: 
//-----------------------------------------------------------------------------
#include <unistd.h>

//#include <iostream.h>
//#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
#include <ctype.h> 
#include <cstdlib>
//#include <conio.h>   // For getch() function 
#include <math.h> 
#include <cmath> // For the sqrt() function  
 
#include <new> 
//#include <alloc.h>   // <- in older compilers !!!!
//----------------------------------------------------------------------------- 
// 
//		All the locally used definitions 
//----------------------------------------------------------------------------- 
 
#define DEBUG_PRINT_ON 1  //Prints the matrices if 1, for print_matrix() 
 
#define NR_END 1 
#define FREE_ARG char* 
 
#define DIMENSION_MISMACH 10 
#define SMALL_MATRIX_ALOCATION 11 
 
 
#define POSITIV_ZERO (float)1e-15 
#define NEGATIV_ZERO -POSITIV_ZERO 

#define TINY_NUM 1.0e-20;    // A small number.
#define HUGE_NUM 1.0e+20;  // A big number

//----------------------------------------------------------------------------- 
// 
//		The nrutil.h definitions from the Numerical Recepies 
//----------------------------------------------------------------------------- 
 
// ----------------------------------------------------------------------------------------
// nrutil.h::

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)

static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ? (dmaxarg1) : (dmaxarg2))

static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ? (dminarg1) : (dminarg2))

static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2))

static float minarg1,minarg2;
#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ? (minarg1) : (minarg2))

static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ? (lmaxarg1) : (lmaxarg2))

static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ? (lminarg1) : (lminarg2))

static int imaxarg1,imaxarg2;
#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ? (imaxarg1) : (imaxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ? (iminarg1) : (iminarg2))

#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))


//----------------------------------------------------------------------------- 
// 
//		The used Matrix structure definition 
//----------------------------------------------------------------------------- 
 
typedef struct {  // definies the Matrix structure: 
					int	 row;	// number of rowes {1...row} 
					int	 col;	//number of coloms {1...col} 
					float  **mtx;   // the alocated matrix [row x col] 
					int   al_row;   // allocated rows 
					int   al_col;   // allocated coloms 
				} Matrix;  
                  		 
//----------------------------------------------------------------------------- 
//----------------------------------------------------------------------------- 
// 
//		The used Tensor structure definition 
//----------------------------------------------------------------------------- 
 
typedef struct {  // definies the Tensor structure: 
					int	  row;	// number of rowes {1...row} 
					int	  col;	//number of coloms {1...col}
					int	  dep;	//number of coloms {1...set} 
					float  ***tns;  // the alocated tensor [row x col x ] 
					int    al_row;  // allocated rows 
					int    al_col;  // allocated coloms
                                        int    al_dep;  // allocated coloms
				} Tensor;  
                  		 
//-----------------------------------------------------------------------------   
//----------------------------------------------------------------------------- 
// 
//		The function defined in MatrixCal.cpp 
//----------------------------------------------------------------------------- 
 
void matrix_alloc(Matrix *x, int row, int col);	/*Allocate float zero atrix*/ 
void eye_matrix_alloc(Matrix *x, int dim);	/*Allocate float unity Matrix*/  
float **matrixa(int nrl, int nrh, int ncl, int nch); /*Alocate a float matrix*/ 
void matrix_delete(Matrix *x);				/*Remouve the alocated matrix*/  
void free_matrixa(float **m, int nrl, int nrh, int ncl, int nch); 

void print_matrix(Matrix *x,char s[80]); /*  usage : printmat(&x,"x"); */ 
void printlog(char s[80]);	/* Mesage printing::usage: printlog("mesage");*/ 
void sg_error(int x);		     /* Error handling function */ 
void sg_warning(int x);		/* Warning handling function */ 
void addmat(Matrix *r,Matrix *x,Matrix *y); /* Addition of two matrices */ 
void submat(Matrix *r,Matrix *x,Matrix *y); /*Substraction of two matrices*/ 
void mulmat(Matrix *r,Matrix *x,Matrix *y); /*Multiplixation of two matrices*/ 
void mulmat_t1(Matrix *r, Matrix *x, Matrix *y);/*Multiplixation two matrices*/ 
void mulmat_t2(Matrix *r, Matrix *x, Matrix *y);/*Multiplixation two matrices*/ 
void cmulmat(Matrix *r,float a, Matrix *x); /* Multiplixation with constant*/ 
int  minv(Matrix *r,Matrix *x); /* Matrix inverse calculation r = 1/x  */
int  ludcmp(float *a, int n, int *indx, float *d); /* for minv() LU decomposition */
void lubksb(float *a, int n, int *indx, float b[]); /* for minv() */
void transpose(Matrix *r, Matrix *x );	/* Transpose  of matrice */ 
void squaremat(Matrix *r, Matrix *x); /* Power of 2::  r = x'*I*x */ 
void sqrtmat(Matrix *r, Matrix *x);  /* sqrt of a matrix */ 
void zeromat(Matrix *r);	    /* Zeroes the matrix */ 
void eyemat(Matrix *r);              /* Unity matrix */
void diagn(Matrix *r,int n,float x); /* Diagonal  matrix nxn <- x*/ 
void maxmatel(Matrix *r, Matrix *x); /* Maximum matrix element selection */ 
void minmatel(Matrix *r, Matrix *x); /* Minimum matrix element selection */ 
void copymat(Matrix *r, Matrix *x); /* Copy matrix x to r:: r=x; */ 
void concat_row(Matrix *r,Matrix *x,Matrix *y);/*Concatenation of a row matrice*/ 
void concat_col(Matrix *r,Matrix *x,Matrix *y);/*Concatenation of a colum matrice*/ 
void extract_row(Matrix *r,int rs,int re,Matrix *x);/* extract from rs to re rews*/ 
void extract_col(Matrix *r,int cs,int ce,Matrix *x);/* extract from rs to re lolums*/ 
float frobenius_norm(Matrix *x);		/* NormF of a matrice or a vector */ 
float covar(Matrix *x);					/* Covariance for vector X{[1]*[N]} */ 
float random_1(long int *idum);		/* random number generation < (+/- 1) */  
//void jacobi(Matrix *rd,Matrix *rv,Matrix *x); /*matrix eigenvalue/eigenvectors calculation */
int  ptest(Matrix *x);  /*Positivity test for a square matrix */
float sign(float x); /* return the signum of x:(-1,0,1) */ 
float *vectora(int nl, int nh);	/* alocate a float vector v[nl,..nh] */

void nrerror(char error_text[]);
float *vector(long nl, long nh);
float **matrix(long nrl, long nrh, long ncl, long nch);
void free_vector(float *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
 
//***************************************************************************** 
 
 #endif   // _MATRIX_CAL_HPP_