/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**			     Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/


/*
	Header file for ``matrix2.a'' library file
*/


#ifndef MATRIX2H
#define MATRIX2H

#include "matrix1.h"

/* Unless otherwise specified, factorisation routines overwrite the
   matrix that is being factorised */

#ifndef ANSI_C

extern	MAT	*BKPfactor(), *CHfactor(), *LUfactor(), *QRfactor(),
		*QRCPfactor(), *LDLfactor(), *Hfactor(), *MCHfactor(),
		*m_inverse();
extern	double	LUcondest(), QRcondest();
extern	MAT	*makeQ(), *makeR(), *makeHQ(), *makeH();
extern	MAT	*LDLupdate(), *QRupdate();

extern	VEC	*BKPsolve(), *CHsolve(), *LUsolve(), *_Qsolve(), *QRsolve(),
		*LDLsolve(), *Usolve(), *Lsolve(), *Dsolve(), *LTsolve(),
		*UTsolve(), *LUTsolve(), *QRCPsolve();

extern  BAND    *bdLUfactor(), *bdLDLfactor();
extern  VEC     *bdLUsolve(), *bdLDLsolve();

extern	VEC	*hhvec();
extern	VEC	*hhtrvec();
extern	MAT	*hhtrrows();
extern	MAT	*hhtrcols();

extern	void	givens();
extern	VEC	*rot_vec();	/* in situ */
extern	MAT	*rot_rows();	/* in situ */
extern	MAT	*rot_cols();	/* in situ */


/* eigenvalue routines */
extern	VEC	*trieig(), *symmeig();
extern	MAT	*schur();
extern	void	schur_evals();
extern	MAT	*schur_vecs();

/* singular value decomposition */
extern	VEC	*bisvd(), *svd();

/* matrix powers and exponent */
MAT  *_m_pow();
MAT  *m_pow();
MAT  *m_exp(), *_m_exp();
MAT  *m_poly();

/* FFT */
void fft();
void ifft();


#else

                 /* forms Bunch-Kaufman-Parlett factorisation for
                        symmetric indefinite matrices */
extern	MAT	*BKPfactor(MAT *A,PERM *pivot,PERM *blocks),
                 /* Cholesky factorisation of A
                        (symmetric, positive definite) */
		*CHfactor(MAT *A),
                /* LU factorisation of A (with partial pivoting) */ 
                *LUfactor(MAT *A,PERM *pivot),
                /* QR factorisation of A; need dim(diag) >= # rows of A */
		*QRfactor(MAT *A,VEC *diag),
                /* QR factorisation of A with column pivoting */
		*QRCPfactor(MAT *A,VEC *diag,PERM *pivot),
                /* L.D.L^T factorisation of A */
		*LDLfactor(MAT *A), 
                /* Hessenberg factorisation of A -- for schur() */
                *Hfactor(MAT *A,VEC *diag1,VEC *diag2),
                /* modified Cholesky factorisation of A;
                        actually factors A+D, D diagonal with no
                        diagonal entry in the factor < sqrt(tol) */
                *MCHfactor(MAT *A,double tol),
		*m_inverse(MAT *A,MAT *out);

                /* returns condition estimate for A after LUfactor() */
extern	double	LUcondest(MAT *A,PERM *pivot),
                /* returns condition estimate for Q after QRfactor() */
                QRcondest(MAT *A);

/* Note: The make..() and ..update() routines assume that the factorisation
        has already been carried out */

     /* Qout is the "Q" (orthongonal) matrix from QR factorisation */
extern	MAT	*makeQ(MAT *A,VEC *diag,MAT *Qout),
                /* Rout is the "R" (upper triangular) matrix
                        from QR factorisation */
		*makeR(MAT *A,MAT *Rout),
                /* Qout is orthogonal matrix in Hessenberg factorisation */
		*makeHQ(MAT *A,VEC *diag1,VEC *diag2,MAT *Qout),
                /* Hout is the Hessenberg matrix in Hessenberg factorisation */
		*makeH(MAT *A,MAT *Hout);

                /* updates L.D.L^T factorisation for A <- A + alpha.u.u^T */
extern	MAT	*LDLupdate(MAT *A,VEC *u,double alpha),
                /* updates QR factorisation for QR <- Q.(R+u.v^T)
		   Note: we need explicit Q & R matrices,
                        from makeQ() and makeR() */
		*QRupdate(MAT *Q,MAT *R,VEC *u,VEC *v);

/* Solve routines assume that the corresponding factorisation routine
        has already been applied to the matrix along with auxiliary
        objects (such as pivot permutations)

        These solve the system A.x = b,
        except for LUTsolve and QRTsolve which solve the transposed system
                                A^T.x. = b.
        If x is NULL on entry, then it is created.
*/

extern	VEC	*BKPsolve(MAT *A,PERM *pivot,PERM *blocks,VEC *b,VEC *x),
		*CHsolve(MAT *A,VEC *b,VEC *x),
		*LDLsolve(MAT *A,VEC *b,VEC *x),
		*LUsolve(MAT *A,PERM *pivot,VEC *b,VEC *x),
		*_Qsolve(MAT *A,VEC *,VEC *,VEC *, VEC *),
		*QRsolve(MAT *A,VEC *,VEC *b,VEC *x),
    		*QRTsolve(MAT *A,VEC *,VEC *b,VEC *x),


     /* Triangular equations solve routines;
        U for upper triangular, L for lower traingular, D for diagonal
        if diag_val == 0.0 use that values in the matrix */

		*Usolve(MAT *A,VEC *b,VEC *x,double diag_val),
		*Lsolve(MAT *A,VEC *b,VEC *x,double diag_val),
		*Dsolve(MAT *A,VEC *b,VEC *x),
		*LTsolve(MAT *A,VEC *b,VEC *x,double diag_val),
		*UTsolve(MAT *A,VEC *b,VEC *x,double diag_val),
                *LUTsolve(MAT *A,PERM *,VEC *,VEC *),
                *QRCPsolve(MAT *QR,VEC *diag,PERM *pivot,VEC *b,VEC *x);

extern  BAND    *bdLUfactor(BAND *A,PERM *pivot),
                *bdLDLfactor(BAND *A);
extern  VEC     *bdLUsolve(BAND *A,PERM *pivot,VEC *b,VEC *x),
                *bdLDLsolve(BAND *A,VEC *b,VEC *x);



extern	VEC	*hhvec(VEC *,u_int,Real *,VEC *,Real *);
extern	VEC	*hhtrvec(VEC *,double,u_int,VEC *,VEC *);
extern	MAT	*hhtrrows(MAT *,u_int,u_int,VEC *,double);
extern	MAT	*hhtrcols(MAT *,u_int,u_int,VEC *,double);

extern	void	givens(double,double,Real *,Real *);
extern	VEC	*rot_vec(VEC *,u_int,u_int,double,double,VEC *); /* in situ */
extern	MAT	*rot_rows(MAT *,u_int,u_int,double,double,MAT *); /* in situ */
extern	MAT	*rot_cols(MAT *,u_int,u_int,double,double,MAT *); /* in situ */


/* eigenvalue routines */

               /* compute eigenvalues of tridiagonal matrix
                  with diagonal entries a[i], super & sub diagonal entries
                  b[i]; eigenvectors stored in Q (if not NULL) */
extern	VEC	*trieig(VEC *a,VEC *b,MAT *Q),
                 /* sets out to be vector of eigenvectors; eigenvectors
                   stored in Q (if not NULL). A is unchanged */
		*symmeig(MAT *A,MAT *Q,VEC *out);

               /* computes real Schur form = Q^T.A.Q */
extern	MAT	*schur(MAT *A,MAT *Q);
         /* computes real and imaginary parts of the eigenvalues
                        of A after schur() */
extern	void	schur_evals(MAT *A,VEC *re_part,VEC *im_part);
          /* computes real and imaginary parts of the eigenvectors
                        of A after schur() */
extern	MAT	*schur_vecs(MAT *T,MAT *Q,MAT *X_re,MAT *X_im);


/* singular value decomposition */

        /* computes singular values of bi-diagonal matrix with
                   diagonal entries a[i] and superdiagonal entries b[i];
                   singular vectors stored in U and V (if not NULL) */
VEC	*bisvd(VEC *a,VEC *b,MAT *U,MAT *V),
               /* sets out to be vector of singular values;
                   singular vectors stored in U and V */
	*svd(MAT *A,MAT *U,MAT *V,VEC *out);

/* matrix powers and exponent */
MAT  *_m_pow(MAT *,int,MAT *,MAT *);
MAT  *m_pow(MAT *,int, MAT *);
MAT  *m_exp(MAT *,double,MAT *);
MAT  *_m_exp(MAT *,double,MAT *,int *,int *);
MAT  *m_poly(MAT *,VEC *,MAT *);

/* FFT */
void fft(VEC *,VEC *);
void ifft(VEC *,VEC *);

#endif


#endif