/***************************************************************************
* Matrix.h                                                                 *
*                                                                          *
* General Vector and Matrix classes, with all the associated methods.      *
*                                                                          *
*   HISTORY                                                                *
*      Name    Date        Description                                     *
*                                                                          *
*      arvo    08/16/2000    Revamped for CIT tools.                       *
*      arvo    10/31/1994    Combined RowVec & ColVec into Vector.         *
*      arvo    06/30/1993    Added singular value decomposition class.     *
*      arvo    06/25/1993    Major revisions.                              *
*      arvo    09/08/1991    Initial implementation.                       *
*                                                                          *
*--------------------------------------------------------------------------*
* Copyright (C) 2000, James Arvo                                           *
*                                                                          *
* 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.  See http://www.fsf.org/copyleft/gpl.html      *
*                                                                          *
* This program is distributed in the hope that it will be useful, but      *
* WITHOUT EXPRESS OR IMPLIED WARRANTY of merchantability or fitness for    *
* any particular purpose.  See the GNU General Public License for more     *
* details.                                                                 *
*                                                                          *
***************************************************************************/
#ifndef __MATRIX_INCLUDED__
#define __MATRIX_INCLUDED__

#include <stream.h>
#include <util/Vector.h>

class Matrix {
    public:
        Matrix( const Matrix & );
        Matrix( int num_rows = 0, int num_cols = 0, float value = 0.0 );
       ~Matrix();
        Matrix &operator=( const Matrix &M );
        Matrix &operator=( float s );
        Vector  GetCol( int col ) const;
        Vector  GetRow( int row ) const;
        void    SetCol( int col, const Vector & );
        void    SetRow( int row, const Vector & );
        Matrix  GetBlock( int imin, int imax, int jmin, int jmax ) const;
        void    SetBlock( int imin, int imax, int jmin, int jmax, const Matrix & );
        void    SetBlock( int imin, int imax, int jmin, int jmax, const Vector & );
        Matrix &SwapRows( int i1, int i2 );
        Matrix &SwapCols( int j1, int j2 );
        void    SetSize( int rows, int cols = 0 );
        double  Cofactor( int i, int j ) const;
        static  const Matrix Null;

    public: // Inlined functions.
        inline float  operator()( int i, int j ) const { return elem[ i * cols + j ]; }
        inline float &operator()( int i, int j )       { return elem[ i * cols + j ]; }
        inline int    Rows  () const { return rows; }
        inline int    Cols  () const { return cols; }
        inline float *Array () const { return elem; }

    private:
        int    rows; // Number of rows in the matrix.
        int    cols; // Number of columns in the matrix.
        float *elem; // Pointer to the actual data.
    };


extern Vector  operator *  ( const Matrix &, const Vector & );
extern Vector  operator *  ( const Vector &, const Matrix & );
extern Vector& operator *= (       Vector &, const Matrix & );
extern Matrix  Outer       ( const Vector &, const Vector & );  // Outer product.
extern Matrix  operator +  ( const Matrix &, const Matrix & );
extern Matrix  operator -  ( const Matrix &                 );
extern Matrix  operator -  ( const Matrix &, const Matrix & );
extern Matrix  operator *  ( const Matrix &, const Matrix & );
extern Matrix  operator *  ( const Matrix &,       float    );
extern Matrix  operator *  (       float  ,  const Matrix & );
extern Matrix  operator /  ( const Matrix &,       float    );
extern Matrix& operator += (       Matrix &, const Matrix & );
extern Matrix& operator *= (       Matrix &,       float    );
extern Matrix& operator *= (       Matrix &, const Matrix & );
extern Matrix& operator /= (       Matrix &,       float    );
extern Matrix  Ident       (       int    n );
extern Matrix  Householder ( const Vector & );
extern Matrix  Rotation    ( const Vector &Axis, float angle );
extern Matrix  Rotation    ( const Vector &Axis, const Vector &Origin, float angle );
extern Matrix  Rotation    (       float, float, float ); // For random 3D rotations.
extern Matrix  Xrotation   (       float    );
extern Matrix  Yrotation   (       float    );
extern Matrix  Zrotation   (       float    );
extern Matrix  Diag        ( const Vector & );
extern Vector  Diag        ( const Matrix & );
extern Matrix  Diag        ( float, float, float );
extern Matrix  Adjoint     ( const Matrix & );
extern Matrix  Adjoint     ( const Matrix &, double &det );
extern Matrix  AATransp    ( const Matrix & );
extern Matrix  ATranspA    ( const Matrix & );
extern double  OneNorm     ( const Matrix & );
extern double  SupNorm     ( const Matrix & );
extern double  Determinant ( const Matrix & );
extern double  Trace       ( const Matrix & );
extern Matrix  Transp      ( const Matrix & );
extern Matrix  Inverse     ( const Matrix & );
extern int     Null        ( const Matrix & );
extern int     Square      ( const Matrix & );
extern Vector  ToVector    ( const Matrix & ); // Only for vector-shaped matrices.

enum pivot_type {
    pivot_off,
    pivot_partial,
    pivot_total
    };

extern int GaussElimination( 
    const Matrix &A, 
    const Vector &b, // This is the right-hand side.
    Vector       &x, // This is the matrix we are solving for.
    pivot_type = pivot_off
    );

extern int LeastSquares( 
    const Matrix &A, 
    const Vector &b, 
    Vector       &x
    );

extern int WeightedLeastSquares( 
    const Matrix &A, 
    const Vector &b, 
    const Vector &w, 
    Vector       &x 
    );

ostream &operator<<( 
    ostream &out, 
    const Matrix &
    );

#endif