#ifndef COVMATRIX_H
#define COVMATRIX_H

// Definitions for generalized CovMatrix class
// $Id: CovMatrix.h,v 1.15 2002/09/19 16:45:18 sachs Exp $
// David Sachs - Fermilab - April 2002

// A CovMatrix is a positive definite symmetric matrix,
// which is stored with the redundant below diagnonal
// elements omitted.  The code for this class is optimized
// for matrices of fairly small order.
//
// There is no specific Vector class associated with the
// CovMatrix class, because the word "Vector" is too common
// to use in this manner, and the standard vector<double>
// class works properly for this purpose.

#include "CovMatrices/BiVector.h"
// BiVector.h includes needed standard headers

namespace CovMatrices
{
  // ********** Error Handler **********
  // The following routine may be replaced by the user
  // The default version prints the error message and exits with code 2
  // This function is called for LOGIC errors only.
  void CovMatrixThrow(const std::string& errormessage); // should not return

  // Covariant vMatrix class
  class CovMatrix
  {
    private:
    int n;	// Order of Matrix
    double* M;	// Elements of Matrix in triagular order

    //  *** Constructors ***
    public:

    // Construct Empty CovMatrix (Default Constructor
    CovMatrix() : n(0), M(0) {}	

    // Construct Uninitialized Matrix of size n
    explicit CovMatrix(int n);

    // Copy Constructor
    CovMatrix(const CovMatrix&);

    // Construct from vector<double>
    // The size of the vector must be of the form (n*(n+1))/2
    explicit CovMatrix(const std::vector<double>&);

    // Construct from BiVector
    // The result is an order 2 Matrix representing the product
    // of the BiVector with itself
    explicit CovMatrix(const BiVector&);

    // Destructor
    ~CovMatrix() { delete[] M; }

    //  *** Assignment operators *** May resize CovMatrix ***
    CovMatrix& operator= (const CovMatrix&);
    CovMatrix& operator= (const std::vector<double>&); 
    CovMatrix& operator= (const BiVector&);

    // Get order of CovMatrix
    inline int size() const {return n;}

    // Resize without initializing elements
    //
    void resize(int);

    // *** Element Access ***
    // These are implemented as functions because subscripting
    // allows only 1 argument
    // These functions abort if a subscript is out of range
    double& operator() (int i, int j);
    const double& operator() (int i, int j) const;

    // Dump elements to vector<double>
    void dump(std::vector<double>&) const;

    // Matrix Inversion
    // Return is false if inversion algorithm failed
    bool invert();	
    bool inverse(CovMatrix& C) const { C = *this; return (C.invert()); }

    // Matrix Determinant - Returns 0. if CovMatrix is not positive definite
    double determinant() const;

    // Eigenvalues and Eigenvectors
    // Returns false if Eigenvalue algorithm failed to converge
    bool diagonalize(std::vector<double>&) const;	// Eignevalues only
    bool diagonalize(std::vector<double>&, std::vector<std::vector<double> >&) const;

    // Print CovMatrix
    void put( std::ostream&) const;

    // Static function to display a vector<double>
    // This function is proveded for BiVector and CovMatrix classes
    // because vector<double> does not have output operator defined
    static void printv( std::ostream&, const std::vector<double>& );

    // Matrix arithmetic
    CovMatrix operator+(const CovMatrix&) const;
    CovMatrix& operator+=(const CovMatrix&);
    CovMatrix operator-(const CovMatrix&) const; // Possibly dangerous
    CovMatrix& operator-=(const CovMatrix&); // Possibly dangerous

    // friend functions
    friend BiVector operator* (const CovMatrix&, const BiVector&);
    friend BiVector operator* (const BiVector&, const CovMatrix&);
    friend std::vector<double> operator* (const std::vector<double>&, const CovMatrix&);
    friend std::vector<double> operator* (const CovMatrix&, const std::vector<double>&);

  };

  // Declarations of non-member functions
  BiVector operator* (const CovMatrix&, const BiVector&);
  BiVector operator* (const BiVector&, const CovMatrix&);
  std::vector<double> operator* (const std::vector<double>&, const CovMatrix&);
  std::vector<double> operator* (const CovMatrix&, const std::vector<double>&);
  inline std::ostream& operator<<(std::ostream& o, const CovMatrix& C)
    { C.put(o); return o; }
}  // namespace CovMatrices
#endif