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

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

//++
// Links	Parents
// Matrix
//--

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

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


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


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


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

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

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

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


//++
Vector& Vector::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 Vector& v1, const Vector& 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;
}

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

//++
// Vector operator* (const Matrix& a, const Vector& 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(Matrix& a, Vector& b)
//
//	Résolution du système A*C = B
//--
{
  if (a.NCol() != b.NRows()) THROW(sizeMismatchErr);
  if (a.NCol() != a.NRows()) THROW(sizeMismatchErr);

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


//++
double LinSolve(const Matrix& a, const Vector& b, Vector& 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;
  Matrix a1(a);
  return Matrix::GausPiv(a1,c);
}

//////////////////////////////////////////////////////////
//++
Vector Vector::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("Vector::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Vector::FitResidus: NULL pointer\n"));
Vector par = gfit.GetParm();
Vector v(*this);
for(int i=0;i<NElts();i++) {
  double x = xorg+i*dx;
  v(i) -= f->Value(&x,par.Data());
}
return v;
}

//++
Vector Vector::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("Vector::FitResidus: size mismatch\n"));
GeneralFunction* f = gfit.GetFunction();
if(f==NULL)
  throw(NullPtrError("Vector::FitResidus: NULL pointer\n"));
Vector par = gfit.GetParm();
Vector v(*this);
for(int i=0;i<NElts();i++) {
  double x = xorg+i*dx;
  v(i) = f->Value(&x,par.Data());
}
return v;
}