#include "machdefs.h"
#include <string.h>
#include "cvector.h"
#include "generalfit.h"
#include "tvector.h"

//++
// Class	OVector
// Lib	Outils++
// include	cvector.h
//
//    Vecteur colonne, en fait une matrice avec une seule colonne, pour faire
//    des opérations matricielles.      
//--

//++
// Links	Parents
// OMatrix
//--

//++
// Titre	Constructeurs
//--

//++
OVector::OVector(int n)
//
//	Constructeur, n = nombre d'éléments.
//--
: OMatrix(n, 1)
{END_CONSTRUCTOR}


//++
OVector::OVector(int n, double* values)
//
//	Constructeur, à partir des valeurs. Pas d'allocation.
//--
: OMatrix(n, 1, values)
{END_CONSTRUCTOR}


//++
OVector::OVector(const OVector& v)
//
//	Constructeur par copie.
//--
: OMatrix(v)
{END_CONSTRUCTOR}


//++
OVector::OVector(const OMatrix& a)
//
//	Constructeur par "copie" à partir d'une matrice, qui doit avoir une seule colonne.
//
//	*Exceptions :*
//	sizeMismatchErr    si la matrice a plus d'une colonne
//--
: OMatrix(a.nr,1)      // OMatrix(a) ne marche pas, pourquoi ???? bug gcc ???
{
  if (a.nc != 1) THROW(sizeMismatchErr);
  memcpy(data, a.data, a.nr * sizeof(double));
  END_CONSTRUCTOR
}

//++
OVector::OVector(const TVector<r_8>& v)
//
//	Constructeur par "copie" a partir d'un TVector<r_8>.
//	Attention, les donnees sont partagees.
//--
: OMatrix(v)
{
}

//++
// Titre	Opérateurs
//--

//++
// double& operator()(int n)
// double const& operator()(int n) const
//	Accès aux éléments.
//--


//++
OVector& OVector::operator = (double x)
//
//	Affectation à partir d'un scalaire. Tous les éléments prennent cette valeur.
//--
{
  for (int i=0; i<nr; i++)
    (*this)(i) = x ;
  return *this;
}

//++
double operator* (const OVector& v1, const OVector& v2)
//
//	Produit scalaire entre deux vecteurs.
//--
{
  if (v1.nr != v2.nr) THROW(sizeMismatchErr);
  double *p    = v1.data;
  double *pEnd = v1.data + v1.ndata;
  double *q    = v2.data;

  double r = 0;

  while (p<pEnd)
    r += *p++ * *q++;

  return r;
}

//++
// OVector operator* (const OVector& v, double b)
//	multiplication par un scalaire.
//--

//++
// OVector operator* (const OMatrix& a, const OVector& b)
//	produit matrice*vecteur.
//--

//++
// Titres	Méthodes
//--

//++
// int NElts() const
//
//	Nombre de lignes du vecteur.
//--

//++
// void Realloc(int n, bool force=false)
//
//	Change la taille du vecteur. Réallocation physique seulement si
//	pas assez de place, ou forcée si force=true.
//--

//++
double LinSolveInPlace(OMatrix& a, OVector& b)
//
//	Résolution du système A*C = B
//--
{
  if (a.NCol() != b.NRows()) THROW(sizeMismatchErr);
  if (a.NCol() != a.NRows()) THROW(sizeMismatchErr);

  return OMatrix::GausPiv(a,b);
}


//++
double LinSolve(const OMatrix& a, const OVector& b, OVector& c)
//
//	Résolution du système A*C = B, avec C retourné dans B
//--
{ 
  if (a.NCol() != b.NRows()) THROW(sizeMismatchErr);
  if (a.NCol() != a.NRows()) THROW(sizeMismatchErr);
  c = b;
  OMatrix a1(a);
  return OMatrix::GausPiv(a1,c);
}

//////////////////////////////////////////////////////////
//++
OVector OVector::FitResidus(GeneralFit& gfit,double xorg,double dx)
//
//	Retourne une classe contenant les residus du fit ``gfit''.
//	La coordonnee de l'element (i) est -> x = xorg + i*dx
//--
{
if(NElts()<=0)
  throw(SzMismatchError("OVector::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("OVector::FitResidus: NULL pointer\n"));
OVector par = gfit.GetParm();
OVector v(*this);
for(int i=0;i<NElts();i++) {
  double x = xorg+i*dx;
  v(i) -= f->Value(&x,par.Data());
}
return v;
}

//++
OVector OVector::FitFunction(GeneralFit& gfit,double xorg,double dx)
//
//	Retourne une classe contenant la fonction du fit ``gfit''.
//	La coordonnee de l'element (i) est -> x = xorg + i*dx
//--
{
if(NElts()<=0)
  throw(SzMismatchError("OVector::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("OVector::FitResidus: NULL pointer\n"));
OVector par = gfit.GetParm();
OVector v(*this);
for(int i=0;i<NElts();i++) {
  double x = xorg+i*dx;
  v(i) = f->Value(&x,par.Data());
}
return v;
}