// This may look like C code, but it is really -*- C++ -*-
//                         C.Magneville          05/99
#ifndef TVECTOR_SEEN
#define TVECTOR_SEEN

#include "tmatrix.h"

namespace PlanckDPC {

class GeneralFit;

template <class T>
class TVector : public TMatrix<T> {
public:

  // Creation / destruction 
  TVector(uint_4 n=1);
  TVector(uint_4 n, T* values,Bridge* br=NULL);
  TVector(const TVector<T>& v);
  TVector(const TVector<T>& v,bool share);
  TVector(const TMatrix<T>& a);

  // Gestion taille/Remplissage
  inline void ReSize(uint_4 n) {TMatrix<T>::ReSize(n,1);} // Reallocation de place

  // Informations pointeur/data
  inline uint_4 NElts() const {return NRows();}
  
  // Acces aux elements
  inline T& operator()(uint_4 n) {return (*this)[n];}
  inline T const& operator()(uint_4 n) const {return (*this)[n];}

  // Operateur d'affectation
  inline TVector& operator = (const TVector& v)
                  {TMatrix<T>::operator =(v); return *this;}
  inline TVector& operator = (T x)
                  {for(uint_4 i=0;i<NRows();i++) (*this)(i)=x; return *this;}

  // Residus et fonction fittees.
  TVector<T> FitResidus(GeneralFit& gfit,double xorg=0.,double dx=1.);
  TVector<T> FitFunction(GeneralFit& gfit,double xorg=0.,double dx=1.);

};

// produit scalaire, matrice*vecteur
template <class T> inline T operator* (const TVector<T>& v1, const TVector<T>& v2)
  {if(v1.NRows() != v2.NRows())
     throw(SzMismatchError("TVector::operator*(TVector& v1,TVector v2) size mismatch"));
   T *p = const_cast<T *>(v1.Data()), *pEnd = p+v1.NElts(),
     *q = const_cast<T *>(v2.Data()), r = 0;
   while (p<pEnd) r += *p++ * *q++;
   return r;}

template <class T> inline TVector<T> operator* (const TMatrix<T>& a, const TVector<T>& b)
                   {return TVector<T>(a * ((TMatrix<T> const&)(b)));}

// Resolution du systeme A*C = B
inline r_8 LinSolveInPlace(TMatrix<r_8>& a, TVector<r_8>& b)
{
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolveInPlace(TMatrix<r_8>,TVector<r_8>) size mismatch"));
return TMatrix<r_8>::GausPiv(a,b);
}

// Resolution du systeme A*C = B, avec C retourne dans B
inline r_8 LinSolve(const TMatrix<r_8>& a, const TVector<r_8>& b, TVector<r_8>& c)
{ 
if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
  throw(SzMismatchError("LinSolve(TMatrix<r_8>,TVector<r_8>) size mismatch"));
c = b;
TMatrix<r_8> a1(a);
return TMatrix<r_8>::GausPiv(a1,c);
}

/////////////////////////////////////////////////////////////////////////
// Classe pour la gestion de persistance
template <class T>
class FIO_TVector : public  PPersist  {
public:
  FIO_TVector();
  FIO_TVector(string const & filename); 
  FIO_TVector(const TVector<T> & obj);
  FIO_TVector(TVector<T> * obj);
  virtual ~FIO_TVector();
  virtual AnyDataObj* DataObj();
  inline operator TVector<T>() { return(*dobj); }
protected :
  virtual void ReadSelf(PInPersist&);           
  virtual void WriteSelf(POutPersist&) const;  
  TVector<T> * dobj;
  bool ownobj;
};

} // Fin du namespace

#endif