//   ======================< VECTOR.CPP >=======================
//   * Class Vector                                            *
//   *       for operations between matrices and vectors       *
//   * Description: Chapter 12                                 *
//   * Scientific C++ Building Numerical Libraries             *
//   *                the Object-Oriented Way                  *
//   * by G. Buzzi-Ferraris                                    *
//   * Addison-Wesley (1993)                                   *
//   ===========================================================

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <math.h>
#include "utility.hpp"
#include "vector.hpp"
#include "matrix.hpp"

const char *const Vector::ERROR =
          "\n========>>> Vector error!!!!!!\n";
int Vector::count = 0;

//   ===========================================================
//   ================= Private functions =======================
//   ===========================================================

//   ********************< Initialize >*************************
//   * Purpose: For initialising various constructors          *
//   * Description: Internal use                               *
//   * Example: Initialize(nc);                                *
//   ***********************************************************
void Vector::Initialize(int size)
   {
   count++;
   whoAmI = count;
   if(size < 1)
      {
      dimensions = 0;
      vector = 0;
      }
   else
      {
      dimensions = size;
      vector = new float [size+1]; 
      if(!vector)
        Error("%s%s%sInitialize",
        Vector::ERROR,ERR_SPACE,ERR_FUNCTION);
      }
}

//   ****************< private constructor >********************
//   * Purpose: Sizing a Vector without initialisation         *
//   *         for internal use                                *
//   * Example: Vector x('*',3);                               *
//   ***********************************************************
Vector::Vector(char,int nc)
   {
   Initialize(nc);
   }

//   ===========================================================
//   ================== Public functions =======================
//   ===========================================================

//   ===========================================================
//   *******************   constructors   **********************
//   ===========================================================

//   **********************< default >**************************
//   * Purpose: Initialising a Vector which can then receive   *
//   *               an assignment                             *
//   * Example: Vector x;                                      *
//   *          x = y;                                         *
//   ***********************************************************
Vector::Vector(void)
   {
   Initialize(0);
   }

//   ******************< Copy-Initializer >*********************
//   * Purpose: Initialising from Vector                       *
//   *          Used in returns                                *
//   * Example: Vector x = b; // simple initialisation         *
//   *          Vector x = b+2.*y;  // with operations         *
//   *          return x; (used implicitly)                    *
//   ***********************************************************
Vector::Vector(const Vector &rval)
   {
   Initialize(rval.dimensions);
   if(dimensions != 0)
      memcpy(vector,rval.vector,
            (dimensions+1)*sizeof(float));
   }

//   ********************< sizing >*****************************
//   * Purpose: Sizing a Vector                                *
//   *          Initialises it to zero                         *
//   * Example: Vector x(3);                                   *
//   *          x(1) = 2.; x(2) = 4.; x(3) = 5.;               *
//   ***********************************************************
Vector::Vector(int nc)
   {
   Initialize(nc);
   if(dimensions != 0)
      memset(vector,0,(dimensions+1)*sizeof(float));
   }

//   **************< sizing and initialising >******************
//   * Purpose: Initialising a Vector                          *
//   * Example: Vector x(3,1.,2.,3.);                          *
//   ***********************************************************
Vector::Vector(int nc,double v1,...)
   {
   Initialize(nc);
   float *w = vector + 1;
   va_list puntaList;
   va_start(puntaList,v1);
   int i;
   *w = Control(v1);
   for(i = 2;i <= nc;i++)
      *++w = Control(va_arg(puntaList,double));
   va_end(puntaList);
   }

//   *********************< from array >************************
//   * Purpose: Initialising a Vector from an array            *
//   * Example: float w[] = {1.,2.,3.,4.,5.,6.};               *
//   *          Vector x(6,w);                                 *
//   ***********************************************************
Vector::Vector(int nc,float *initvalues)
   {
   Initialize(nc);
   if(dimensions != 0)
      {
      float *w = vector;
      w++;
      memcpy(w,initvalues,nc*sizeof(float));
      }
   }

//   ********************< ASCII FILE >*************************
//   * Purpose: Initialising a Vector from FILE                *
//   * Example: Vector x("VECT.DAT");                          *
//   ***********************************************************
Vector::Vector(char *filevector)
   {
   if((bzzFileIn = fopen(filevector,"r"))==NULL)
       Error("%s%s%s%sConstructor",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_FUNCTION);
   int nc;
   fscanf(bzzFileIn,"%d",&nc);
   Initialize(nc);
   float *w = vector;
   int i;
   if(dimensions != 0)
      for(i = 1;i <= nc;i++)
        fscanf(bzzFileIn,"%f",++w);
   fclose(bzzFileIn);
   }

//   ********************< from binary FILE >*******************
//   * Purpose: Initialising a Vector from unformatted FILE    *
//   * Description: Constructor for recovering a Vector        *
//   *              saved on FILE with Save('*',file name)     *
//   * Example: Vector v('*',"VET.BIN");                       *
//   ***********************************************************
Vector::Vector(char,char *filevector)
   {
   if((bzzFileIn=fopen(filevector,"rb"))==NULL)
       Error("%s%s%s%sFILE",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_CONSTRUCTOR);
   int nc;
   if(fread(&nc,sizeof(int),1,bzzFileIn) != 1)
      Error("when reading dimensions");
   Initialize(nc);
   if(dimensions != 0)
      if(fread(vector+1,sizeof(float)*nc,1,bzzFileIn) != 1)
         Error("%s%s%sConstructor",
         Vector::ERROR,ERR_READING_FILE,ERR_CONSTRUCTOR);
   fclose(bzzFileIn);
   }


//   ********************< subvector >**************************
//   * Purpose: Subvector with the same beginning and nc <= n  *
//   * Example: Vector vNew(5,vOld);                           *
//   ***********************************************************
Vector::Vector(int nc,const Vector &rval)
   {
   if(nc < 1 || nc > rval.dimensions)
      Error("%s%s%s subvector",
      Vector::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   Initialize(nc);
   memcpy(vector,rval.vector,(nc+1)*sizeof(float));
   } 

//   ********************< subvector >**************************
//   * Purpose: Subvector with beginning at ielem              *
//   *          and length nc<=n-ielem+1                       *
//   * Example: Vector vNew(3,2,vOld);                         *
//   ***********************************************************
Vector::Vector(int nc,int ielem,const Vector &rval)
   {
   if(nc < 1 || ielem < 1 || (nc+ielem-1) > rval.dimensions)
      Error("%s%s%sSubvector ",
      Vector::ERROR,ERR_CHECK_DIMENSION,ERR_CONSTRUCTOR);
   Initialize(nc);
   float *r = rval.vector+ielem-1;
   memcpy(vector,r,(dimensions+1)*sizeof(float));
   }

//   ********************< destructor >*************************
//   * Purpose: Deinitialising Vector                          *
//   ***********************************************************
Vector::~Vector(void)
   {
   delete vector;
   }

//   ===========================================================
//   ***********  Non-modifying access functions    ************
//   ===========================================================

//   *********************< GetValue >**************************
//   * Purpose: Receiving values with control                  *
//   * Description: Return is not a reference!                 *
//   *              The user cannot access the data.           *
//   * Example: x = y.GetValue(3);                             *
//   ***********************************************************
float Vector::GetValue(int i) const
   {
   if(i < 1 || i > dimensions )
      Error("%s%s%sValue",
      Vector::ERROR,ERR_RANGE,ERR_FUNCTION);
   return vector[i];
   }

//   ********************< GetVector >**************************
//   * Purpose: Creating a subvector                           *
//   *        (similar to constructor Vector(nc,ielem,Vector)  *
//   *        It is better to use constructor                  *
//   * Description: nc: dimensions of new vector               *
//   *              ielem: the element from where it starts    *
//   *              rval: Vector to be used                    *
//   * Example: x = GetVector(5,2,y);                          *
//   *          x has dimension 5                              *
//   *          and starts from element 2 of y                 *
//   ***********************************************************
Vector GetVector(int nc,int ielem,const Vector &rval)
   {
   if(nc < 1 || ielem < 1 || (nc+ielem-1) > rval.dimensions)
      Error("%s%s%sGetVector",
      Vector::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   Vector result('*',nc);
   float *r = rval.vector+ielem-1;
   memcpy(result.vector,r,(nc+1)*sizeof(float));
   return result;
   }

//   ===========================================================
//   *************    Modifying access functions ***************
//   ===========================================================

//   ********************< operator () >************************
//   * Purpose: Receiving and assigning values with control    *
//   * Example: x = v(5); v(7) = 3.;                           *
//   ***********************************************************
   float &Vector::operator () (int i)
       {
       if(i < 1 || i > dimensions )
          Error("%s%s%s()",
          Vector::ERROR,ERR_RANGE,ERR_OPERATOR);
       return vector[i];
       }

//   *********************< SetValue >**************************
//   * Purpose: Assigning values with control                  *
//   * Description: Not with a reference!                      *
//   *              The user cannot access the data.           *
//   * Example: x.SetValue(3,3.14);                            *
//   ***********************************************************
void Vector::SetValue(int i,float val)
   {
   if(i < 1 || i > dimensions)
      Error("%s%s%sSetValue",
      Vector::ERROR,ERR_RANGE,ERR_FUNCTION);
   vector[i] = val;
   }

//   ********************< SetVector >**************************
//   * Purpose: Assigning coefficients for a portion           *
//   * Description: v.SetVector(3,w);                          *
//   *         substitutes w starting from coefficient 3 of v  *
//   ***********************************************************
void Vector::SetVector(int ielem,const Vector &rval)
   {
   if(ielem < 1 || ielem > dimensions
                  ||
      rval.dimensions > dimensions-ielem+1)
        Error("%s%s%sSetVector",
        Vector::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
      float *r = vector+ielem-1;
      memcpy(r,rval.vector,
           (rval.dimensions+1)*sizeof(float));
   }

//   ===========================================================
//   ****************   assignment operator ********************
//   ===========================================================

//   *************************< = >*****************************
//   * Purpose: Assignment                                     *
//   * Example: x = b;    // simple                            *
//   *          x = y + z;    // with operations               *
//   *          x = b = c + d;  // multiple                    *
//   ***********************************************************
Vector &Vector::operator = (const Vector &rval)
   {
   int who = whoAmI;
   if(dimensions != rval.dimensions)
      {
      delete vector;
      Initialize(rval.dimensions);
      count--;
      }
   whoAmI = who;
   memcpy(vector,rval.vector,(dimensions+1)*sizeof(float));
   return *this;
   }

//   ===========================================================
//   ************   operator for creating vectors   ************
//   ===========================================================

//   ************************< && >*****************************
//   * Purpose: Constructing one Vector from two Vectors       *
//   *          Putting them in series                         *
//   * Description: I use && and not || because Vector         *
//   *              is a column vector                         *
//   * Example:  x = a&&b&&c;                                  *
//   ***********************************************************
Vector operator && (const Vector &lval,const Vector &rval)
   {
   // sizing without initialisation
   Vector result('*',lval.dimensions+rval.dimensions);
   float *w = lval.vector;
   float *z = rval.vector;
   float *r = result.vector;
   memcpy(r,w,(lval.dimensions+1)*sizeof(float));
   r = result.vector+lval.dimensions+1;
   z = z+1;
   memcpy(r,z,(rval.dimensions+1)*sizeof(float));
   return result;
   }

//   ===========================================================
//   ****************   operators for tests   ******************
//   ===========================================================

//   ************************< == >*****************************
//   * Purpose: Comparing two Vectors                          *
//   * Description: Uses the function memcmp                   *
//   * Example: if(w == v)                                     *
//   ***********************************************************
char operator == (const Vector &lval,const Vector &rval)
   {
   if(lval.whoAmI == rval.whoAmI) return 1;
   char ch = 1;
   if(lval.dimensions != rval.dimensions) ch = 0;
   else
      {
      if(memcmp(lval.vector+1,rval.vector+1,
              rval.dimensions*sizeof(float))==0)ch=1;
      else ch=0;
      }
   return ch;
   }

//   ************************< != >*****************************
//   * Purpose: Testing whether two Vectors are different      *
//   * Description: Uses memcmp                                *
//   * Example: if(v != w)                                     *
//   ***********************************************************
char operator != (const Vector &lval,const Vector &rval)
   {
   if(lval.whoAmI == rval.whoAmI) return 0;
   char ch = 0;
   if(lval.dimensions != rval.dimensions)ch = 1;
   else
      ch = memcmp(lval.vector+1,rval.vector+1,
              rval.dimensions*sizeof(float));
   return ch;
   }

//   ===========================================================
//   ====================   OPERATIONS   =======================
//   ===========================================================

//   ===========================================================
//   ***********************   Sum   ***************************
//   ===========================================================

//   ************************< Sum >****************************
//   * Purpose: Addition between vectors                       *
//   * Description: Uses Sum from utility.cpp                  *
//   * Example: Sum(a,b,&c); c = a + b;                        *
//   ***********************************************************
void Sum(const Vector &lval,const Vector &rval,
        Vector *result)
   {
   if(lval.dimensions != rval.dimensions)
       Error("%s%s%s+",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(result->whoAmI == lval.whoAmI)
      (*result) += rval;
   else if(result->whoAmI == rval.whoAmI)
      Sum(lval,result);
   else
      {
      ChangeDimensions(lval.dimensions,result,1);
      Sum(result->dimensions,lval.vector+1,
         rval.vector+1,result->vector+1); // UTILITY.CPP
      }
   }

//   *************************< + >*****************************
//   * Purpose: Addition of two Vectors                        *
//   * Description: Uses Sum(x,y,&z);                          *
//   * Example: y = a + x;                                     *
//   ***********************************************************
Vector operator + (const Vector &lval,const Vector &rval)
   {
   Vector result('*',lval.dimensions);
   Sum(lval,rval,&result);
   return result;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition between vectors                       *
//   * Description: The result goes in the left vector         *
//   * Example: Sum(&a,b); a += b;                             *
//   ***********************************************************
void Sum(Vector *lvalAndResult,const Vector &rval)
   {
   (*lvalAndResult) += rval;
   }

//   ************************< += >*****************************
//   * Purpose: Addition of two Vectors with +=                *
//   * Description: Better than a = a + x;                     *
//   * Example: a += x;                                        *
//   ***********************************************************
Vector &Vector::operator += (const Vector &rval)
   {
   if(dimensions != rval.dimensions)
       Error("%s%s%s+=",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(whoAmI == rval.whoAmI)
      Sum(this);
   else
      Sum(dimensions,vector+1,rval.vector+1); // UTILITY.CPP
   return *this;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition between vectors                       *
//   * Description: The result goes in the right vector        *
//   * Example: Sum(b,&a); a = b + a;                          *
//   ***********************************************************
void Sum(const Vector &lval,Vector *rvalAndResult)
   {
   (*rvalAndResult) += lval;
   }

//   ************************< Sum >****************************
//   * Purpose: Addition between two identical vectors         *
//   * Description: The result replaces the vector             *
//   * Example: Sum(&a); a = a + a;                            *
//   ***********************************************************
void Sum(Vector *lvalRvalAndResult)
   {
   Sum(lvalRvalAndResult->dimensions,
       lvalRvalAndResult->vector+1);
   }

//   ===========================================================
//   ********************   Difference   ***********************
//   ===========================================================

//   **********************< Difference >***********************
//   * Purpose: Difference between vectors                     *
//   * Description: Uses Difference from utility.cpp           *
//   * Example: Difference(a,b,&c); c = a - b;                 *
//   ***********************************************************
void Difference(const Vector &lval,const Vector &rval,
              Vector *result)
   {
   if(lval.dimensions!=rval.dimensions)
       Error("%s%s%s-",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(result->whoAmI == lval.whoAmI)
      (*result) -= rval;
   else if(result->whoAmI == rval.whoAmI)
      Difference(lval,result);
   else
      {
      ChangeDimensions(lval.dimensions,result,1);
      Difference(result->dimensions,lval.vector + 1,
         rval.vector + 1,result->vector + 1); //UTILITY.CPP
      }
   }

//   *************************< - >*****************************
//   * Purpose: Difference between two vectors                 *
//   * Example: y = z - x;                                     *
//   ***********************************************************
Vector operator - (const Vector &lval,const Vector &rval)
   {
   Vector result('*',lval.dimensions);
   Difference(lval,rval,&result);
   return result;
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between vectors                     *
//   * Description: The result goes in the left vector         *
//   * Example: Difference(&a,b); a = a - b;                   *
//   ***********************************************************
void Difference(Vector *lvalAndResult,const Vector &rval)
   {
   (*lvalAndResult) -= rval;
   }

//   ************************< -= >*****************************
//   * Purpose: Difference between two vectors with -=         *
//   * Description: Better than a = a - x;                     *
//   * Example: a -= x;                                        *
//   ***********************************************************
Vector &Vector::operator -= (const Vector &rval)
   {
   if(dimensions != rval.dimensions)
       Error("%s%s%s+=",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(whoAmI == rval.whoAmI)
      Difference(this);
   else
      Difference(dimensions,vector + 1,
             rval.vector +1); // in UTILITY.CPP
   return *this;
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between vectors                     *
//   * Description: The result goes in the right vector        *
//   * Example: Difference(a,&b); b = a - b;                   *
//   ***********************************************************
void Difference(const Vector &lval,Vector *rvalAndResult)
   {
   if(lval.dimensions != rvalAndResult->dimensions)
       Error("%s%s%sDiffernece",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   Difference(lval.dimensions,lval.vector+1,
             rvalAndResult->vector+1,1);
   }

//   **********************< Difference >***********************
//   * Purpose: Difference between identical vectors           *
//   * Description: null Vector                                *
//   * Example: Difference(&a); a = a - a;                     *
//   ***********************************************************
void Difference(Vector *lvalRvalAndResult)
   {
   float *w = lvalRvalAndResult->vector;
   for(int i = 1;i <= lvalRvalAndResult->dimensions;i++)
       *++w = 0.;
   }

//   ===========================================================
//   ***********************   Minus   *************************
//   ===========================================================

//   ************************< Minus >**************************
//   * Purpose: Unary minus           *
//   * Example: Minus(a,&b); b = -a;                           *
//   ***********************************************************
void Minus (const Vector &rval,Vector *result)
   {
   ChangeDimensions(rval.dimensions,result,1);
   float *w = rval.vector;
   float *r = result->vector;
   for(int i = 1;i <= rval.dimensions;i++)
       *++r = -(*++w);
   }

//   *************************< - >*****************************
//   * Purpose: Changing the sign of a Vector                  *
//   * Example: y = -x;                                        *
//   ***********************************************************
Vector Vector::operator -() // unary minus
   {
   Vector result('*',dimensions);
   Minus(*this,&result);
   return result;
   }

//   ************************< Minus >**************************
//   * Purpose: Unary minus           *
//   * Example: Minus(&a); a = -a;                             *
//   ***********************************************************
void Minus (Vector *rvalAndResult)
   {
   float *r = rvalAndResult->vector;
   for(int i = 1;i <= rvalAndResult->dimensions;i++)
       r[i] = -r[i];
   }    

//   ===========================================================
//   ***********************   Product   ***********************
//   ===========================================================

//   ***********************< Product >*************************
//   * Purpose: The product between a matrix and a vector      *
//   * Example: Product(A,x,&y); y = A*x;                      *
//   ***********************************************************
void Product(const Matrix &lval,const Vector &rval,
            Vector *result)
   {
   if(lval.numColumns!=rval.dimensions)
       Error("%s%s%s*",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   if(rval.whoAmI == result->whoAmI)
      {
      Vector aux;
      Product(lval,rval,&aux); // recursive call 
      Swap(&aux,result); // avoids copying
      }
   else
      {
      ChangeDimensions(lval.numRows,result,1);
      float *r = result->vector;
      for(int i = 1;i <= lval.numRows;i++)
        r[i] = Dot(rval.dimensions,lval.matrix[i] +1,
                rval.vector + 1);
      }
   }

//   *************************< * >*****************************
//   * Purpose: The product between a matrix and a vector      *
//   * Example: y = A*x;                                       *
//   ***********************************************************
Vector operator * (const Matrix &lval,const Vector &rval)
   {
    Vector result('*',lval.numRows);
    Product(lval,rval,&result);
    return result;
   }

//   ***********************< Product >*************************
//   * Purpose: The product between a matrix and a vector      *
//   * Example: Product(A,&x); x = A*x;                        *
//   ***********************************************************
void Product(const Matrix &lval,Vector *rvalAndResult)
   {
   Product(lval,*rvalAndResult,rvalAndResult);
   }

//   ***********************< Product >*************************
//   * Purpose: The product between a float and a vector       *
//   * Example: Product(3.,x,&y); y = 3.*x;                    *
//   ***********************************************************
void Product(float lval,const Vector &rval,Vector *result)
   {
   if(rval.whoAmI == result->whoAmI)
      {
   Product(result->dimensions,lval,
          result->vector + 1); // in UTILITY.CPP
      }
   else
      {
      ChangeDimensions(rval.dimensions,result,1);
      Product(rval.dimensions,lval,rval.vector + 1,
                   result->vector + 1); // in UTILITY.CPP
      }
   }

//   *************************< * >*****************************
//   * Purpose: The product between a float and a vector       *
//   * Example: y = 3.*x;                                      *
//   ***********************************************************
Vector operator * (float lval,const Vector &rval)
   {
   Vector result('*',rval.dimensions);
   Product(lval,rval,&result);
   return result;
   }

//   ***********************< Product >*************************
//   * Purpose: The product between a float and a vector       *
//   * Example: Product(3.,&x); x = 3.*x;                      *
//   ***********************************************************
void Product(float lval,Vector *rvalAndResult)
   {
   Product(rvalAndResult->dimensions,lval,
          rvalAndResult->vector + 1); // in UTILITY.CPP
   }

//   ************************< *= >*****************************
//   * Purpose: The product between a float and                *
//   *          a vector with *=                               *
//   * Description: To be used preferably!                     *
//   * Example: y *= 3.;                                       *
//   ***********************************************************
Vector &Vector::operator *= (float rval)
   {
   Product(dimensions,rval,vector +1); // in UTILITY.CPP
   return *this;
   }

//   ===========================================================
//   ******************   TProduct, Dot   **********************
//   ===========================================================

//   *********************< TProduct >**************************
//   * Purpose: Dot product between two vectors (xTy)          *
//   * Description: Uses Dot of UTILITY                        *
//   * Example: TProduct(x,y,&c); c = xTy;                     *
//   ***********************************************************
void TProduct(const Vector &lval,const Vector &rval,
             float *result)
   {
   if(lval.dimensions!=rval.dimensions)
       Error("%s%s%sDot",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_OPERATOR);
   *result = Dot(lval.dimensions,lval.vector + 1,
                 rval.vector +1); // in UTILITY.CPP
   }

//   ************************< Dot >****************************
//   * Purpose: Dot product between two vectors (xTy)          *
//   * Description: Uses TProduct                              *
//   * Example: c = Dot(x,y);                                  *
//   ***********************************************************
float Dot(const Vector &lval,const Vector &rval)
   {
   float result;
   TProduct(lval,rval,&result);
   return result;
   }

//   *************************< % >*****************************
//   * Purpose: Dot product between two vectors (xTy)          *
//   * Description: Uses TProduct                              *
//   * Example: w = x%y;                                       *
//   ***********************************************************
float operator % (const Vector &lval,const Vector &rval)
   {
   float result;
   TProduct(lval,rval,&result);
   return result;
   }

//   *********************< TProduct >**************************
//   * Purpose: ATx: Matrix A transposed and multiplied        *
//   *               by Vector x                               *
//   * Description: the matrix is scanned row by row           *
//   * Example: TProduct(A,x,&y);                              *
//   ***********************************************************
void TProduct(const Matrix &lval,
             const Vector &rval,Vector *result)
   {
   if(lval.numRows != rval.dimensions)
       Error("%s%s%s%TProduct",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   if(rval.whoAmI == result->whoAmI)
      {
      Vector aux;
      TProduct(lval,rval,&aux); // recursive call
      Swap(&aux,result); // avoids copying
      }
   else
      {
      ChangeDimensions(lval.numColumns,result,1);
      int i,j;
      float *r = result->vector;
      float *left = lval.matrix[1] + 1;
      float *right = rval.vector + 1;
      for(j = 1;j <= lval.numColumns;j++)
        r[j] = 0.;
      for(i = 1;i <= lval.numRows;i++)
        {
	for(j = 1;j<= lval.numColumns;j++)
           r[j] = Control(r[j] + (*left++)*(*right));
        right++;
        }
      }
   }

//   *************************< % >*****************************
//   * Purpose: ATx: Matrix A transposed and multiplied        *
//   *               by Vector x                               *
//   * Description: The matrix is scanned row by row           *
//   * Example: y = A%x;                                       *
//   ***********************************************************
Vector operator % (const Matrix &lval,const Vector &rval)
   {
   Vector result('*',lval.numColumns);
   TProduct(lval,rval,&result);
   return result;
   }

//   ***********************< TProduct >************************
//   * Purpose: The product between a transposed matrix        *
//   *          and a vector                                   *
//   * Example: TProduct(A,&x); // x = ATx;                    *
//   ***********************************************************
void TProduct(const Matrix &lval,Vector *rvalAndResult)
   {
   TProduct(lval,*rvalAndResult,rvalAndResult);
   }

//   ===========================================================
//   ********************   ProductT   *************************
//   ===========================================================

//   *********************< ProductT >**************************
//   * Purpose: The product between a column vector and        *
//   *         a row vector                                    *
//   * Example: ProductT(x,y,&A); A = xyT;                     *
//   ***********************************************************
void ProductT(const Vector &lval,const Vector &rval,
             Matrix *result)
   {
   if(lval.dimensions != rval.dimensions)
       Error("%s%s%sProductT",
       Vector::ERROR,ERR_CHECK_DIMENSION,ERR_FUNCTION);
   ChangeDimensions(lval.dimensions,
                  lval.dimensions,result,1);
   float *r = result->matrix[1] + 1;
   float *left = lval.vector + 1;
   float *right = rval.vector;
   for(int i = 1;i <= lval.dimensions;i++)
     {
     for(int j = 1;j <= rval.dimensions;j++)
       *r++ = Control(*left * right[j]);
     left++;
     }
   }

//   ************************< ->* >****************************
//   * Purpose: The product between a Vector                   *
//   *          and a transposed Vector                        *
//   * Example: Matrix A = y ->* x;                            *
//   ***********************************************************
Matrix operator ->* (const Vector &lval,const Vector &rval)
   {
   Matrix result('*',lval.dimensions,lval.dimensions);
   ProductT(lval,rval,&result);
   return result;
   }

//   ===========================================================
//   **********************   Division   ***********************
//   ===========================================================

//   ***********************< Division >************************
//   * Purpose: Dividing the coefficients of a vector          *
//   *          by a float                                     *
//   * Description: Uses Division of utility.cpp               *
//   * Example: Division(x,3.,&y); y = x/3.;                   *
//   ***********************************************************
void Division(const Vector &lval,float rval,Vector *result)
   {
   if(lval.whoAmI == result->whoAmI)
      {
      Division(result->dimensions,
           result->vector+1,rval); //in UTILITY.CPP
      }
   else
      {
      ChangeDimensions(lval.dimensions,result,1);
      Division(lval.dimensions,lval.vector + 1,rval,
                   result->vector + 1); // in UTILITY.CPP
      }
   }

//   *************************< / >*****************************
//   * Purpose: Dividing the coefficients of a vector          *
//   *          by a float                                     *
//   * Description: Uses Division of utility.cpp               *
//   * Example: y = x/3.;                                      *
//   ***********************************************************
Vector operator / (const Vector &lval,float rval)
   {
   Vector result('*',lval.dimensions);
   Division(lval.dimensions,lval.vector + 1,rval,
                   result.vector + 1); // in UTILITY.CPP
   return result;
   }

//   ***********************< Division >************************
//   * Purpose: Dividing the coefficients of a vector          *
//   *          by a float                                     *
//   * Description: Uses Division of utility.cpp               *
//   * Example: Division(&x,3.); x = x/3.;                     *
//   ***********************************************************
void Division(Vector *lvalAndResult,float rval)
   {
   Division(lvalAndResult->dimensions,
           lvalAndResult->vector+1,rval);
   }

//   ************************< /= >*****************************
//   * Purpose: Divides and modifies a Vector by a float       *
//   * Description: Uses Division of utility.cpp               *
//   * Example: x /= 2.; // x = x/2.;                          *
//   ***********************************************************
Vector &Vector::operator /= (float rval)
   {
   Division(dimensions,vector + 1,rval);//in UTILITY.CPP
   return *this;
   }

//   ===========================================================
//   *                 Non-modifying functions                 *
//   ===========================================================

//   ***********************< Print >***************************
//   * Purpose: Prints on bzzFileOut the elements of a Vector  *
//   * Description: Prints on bzzFileOut                       *
//   * Example: v.Print("Vector v");                           *
//   ***********************************************************
void Vector::Print(char *msg)
   {
   fprintf(bzzFileOut,"\nVector No.%d",whoAmI);
   if(*msg)fprintf(bzzFileOut,"  %s",msg);
   fprintf(bzzFileOut,"\nSize %d\n",dimensions);
   for(int i = 1;i <= dimensions;i++)
        fprintf(bzzFileOut," %12.4e",vector[i]);
   }

//   ************************< Save >***************************
//   * Purpose: Saving on file a Vector formatted in ASCII     *
//   * Description: Can be read with a constructor             *
//   *              or with Recover                            *
//   * Example: v.Save("VET.DAT");                             *
//   ***********************************************************
void Vector::Save(char *filevector) // formatted
   {
   FILE *bzz;
   if((bzz = fopen(filevector,"w")) == NULL)
       Error("%s%s%s%sSave",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_FUNCTION);

   fprintf(bzz," %d\n",dimensions);
   for(int i = 1;i <= dimensions;i++)
      fprintf(bzz," %15.7e",vector[i]);
   fclose(bzz);
   }

//   ************************< Save >***************************
//   * Purpose: Saving a Vector on an unformatted file         *
//   * Description: Can be read with a constructor             *
//   *              or with Recover                            *
//   * Example: v.Save('*',VET.BIN");                          *
//   ***********************************************************
void Vector::Save(char,char *filevector) // unformatted
   {
   FILE *bzz;
   if((bzz = fopen(filevector,"wb")) == NULL)
       Error("%s%s%s%sSave",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_FUNCTION);

   if(fwrite(&dimensions,sizeof(int),1,bzz) != 1)
      Error("I could not write ");

   if(fwrite(vector+1,sizeof(float)*dimensions,1,bzz) != 1)
      Error("I could not write ");

   fclose(bzz);
   }

//   ===========================================================
//   *                       Max Min                           *
//   ===========================================================

//   ************************< Max >****************************
//   * Purpose: The maximum of a Vector                        *
//   * Description: Uses the function defined  in UTILITY.CPP  *
//   * Example: float xmax = v.Max();                          *
//   *          float xmax = v.Max(&imax);                     *
//   ***********************************************************
float Vector::Max(int *imax)
   {
   float xmax = ::Max(dimensions,vector + 1,imax);
   if(imax != 0)(*imax)++;
   return xmax;
   }

//   **********************< MaxAbs >***************************
//   * Purpose: The absolute maximum of a Vector               *
//   * Description: Uses the function defined  in UTILITY.CPP  *
//   * Example: float xmax = v.MaxAbs();                       *
//   *          float xmax = v.MaxAbs(&imax);                  *
//   ***********************************************************
float Vector::MaxAbs(int *imax)
   {
   float xmax = ::MaxAbs(dimensions,vector + 1,imax);
   if(imax != 0) (*imax)++;
   return xmax;
   }

//   ************************< Min >****************************
//   * Purpose: The minimum of a Vector                        *
//   * Description: Uses the function defined  in UTILITY.CPP  *
//   * Example: float xmin = v.Min();                          *
//   *          float xmin = v.Min(&imin);                     *
//   ***********************************************************
float Vector::Min(int *imin)
   {
   float xmin = ::Min(dimensions,vector + 1,imin);
   if(imin != 0) (*imin)++;
   return xmin;
   }

//   *********************** MinAbs >***************************
//   * Purpose: The absolute minimum of a Vector               *
//   * Description: Uses the function described in UTILITY.CPP *
//   * Example: float xmin = v.MinAbs();                       *
//   *          float xmin = v.MinAbs(&imin);                  *
//   ***********************************************************
float Vector::MinAbs(int *imin)
   {
   float xmin = ::MinAbs(dimensions,vector + 1,imin);
   if(imin != 0) (*imin)++;
   return xmin;
   }

//   ===========================================================
//   *                        Norms                            *
//   ===========================================================

//   ***********************< Norm1 >***************************
//   * Purpose: Calculating the norm-1                         *
//   * Description: Addition of absolute values of v           *
//   * Example: float norm1 = v.Norm1();                       *
//   ***********************************************************
float Vector::Norm1(void)
{
float norm = 0.,*x = vector;
for(int i = 1;i <= dimensions;i++)
   norm+=Abs(*++x);
return norm;
}

//   ***********************< Norm2 >***************************
//   * Purpose: norm-2 or Euclidean norm                       *
//   * Description: square root of the sum of squares          *
//   *              Chapter 8                                  *
//   * Example: float norm2 = v.Norm2();                       *
//   ***********************************************************
float Vector::Norm2(void)
   {
   float *x = vector + 1;
   return SqrtSumSqr(dimensions,x); //in UTILITY.CPP
   }

//   ***********************< NormI >***************************
//   * Purpose: Norm-I                                         *
//   * Description: Absolute maximum                           *
//   * Example: float normi = v.NormI();                       *
//   ***********************************************************
float Vector::NormI(void)
{
return MaxAbs();
}

//   ===========================================================
//   *                 Modifying functions                     *
//   ===========================================================

//   **********************< Delete >***************************
//   * Purpose: Eliminate a vector of no use                   *
//   *          without leaving the scope                      *
//   * Example: Delete(&x);                                    *
//   ***********************************************************
void Delete(Vector *result)
   {
   delete result->vector;
   result->dimensions = 0;
   result->vector = 0;
   }

//   ***************< ChangeDimensions >************************
//   * Purpose: Change the dimensions of a Vector              *
//   * Example: ChangeDimensions(dim,&a);                      *
//   ***********************************************************
void ChangeDimensions(int dim,Vector *result,char zero)
   {
   int who = result->whoAmI;
   if(dim != result->dimensions)
        {
        delete result->vector;
        result->Initialize(dim);
        Vector::count--;
        }
   result->whoAmI = who;
   if(zero == 0 && result->dimensions != 0)
      memset(result->vector,0,
            (result->dimensions + 1)*sizeof(float));
   }

//   **********************< Recover >**************************
//   * Purpose: Recovering a vector stored with                *
//   *          formatted Save(file name)                      *
//   * Description: It is not necessary to have the same       *
//   *              dimensions as the saved vector             *
//   * Example: Recover(&x,nomefile);                          *
//   ***********************************************************
void Recover(Vector *result,char *filevector)
   {
   if((bzzFileIn = fopen(filevector,"r"))==NULL)
       Error("%s%s%s%sRecover",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_FUNCTION);
   int nc;
   fscanf(bzzFileIn,"%d",&nc);
   ChangeDimensions(nc,result,1);
   float *w = result->vector;
   int i;
   for(i = 1;i <= nc;i++)
       fscanf(bzzFileIn,"%f",++w);
   fclose(bzzFileIn);
   }

//   **********************< Recover >**************************
//   * Purpose: Recovering a vector stored by                  *
//   *          unformatted Save('*', file name)               *
//   * Description: It is not necessary to have the same       *
//   *              dimensions as the saved vector             *
//   * Example: Recover(&x,'*',file name);                     *
//   ***********************************************************
void Recover(Vector *result,char,char *filevector)
   {
   if((bzzFileIn=fopen(filevector,"rb"))==NULL)
       Error("%s%s%s%sRecover",
       Vector::ERROR,ERR_OPEN_FILE,filevector,ERR_FUNCTION);
   int nc;
   if(fread(&nc,sizeof(int),1,bzzFileIn) != 1)
      Error("reading in dimensions");
   ChangeDimensions(nc,result,1);
   if(nc < 1)
       {fclose(bzzFileIn);return;}
   float *v = result->vector;
   if(fread(v+1,sizeof(float)*nc,1,bzzFileIn) != 1)
       Error("%s%s%sRecover",
       Vector::ERROR,ERR_READING_FILE,ERR_FUNCTION);
   fclose(bzzFileIn);
   }

//   *********************< Normalize >*************************
//   * Purpose: Normalising a vector (xTx = 1)                 *
//   *          Returns the norm value                         *
//   * Example: float norm = Normalize(&x);                    *
//   ***********************************************************
float Normalize(Vector *result)
   {
   float norm = result->Norm2();
   if(norm==0.)
      return norm;
   float *w = result->vector;
   for(int i = 1;i <= result->dimensions;i++)
      *++w /= norm;
   return norm;
   }

//   **********************< Reverse >**************************
//   * Purpose: Inverting the elements of a Vector             *
//   * Description: Swaps extreme elements                     *
//   * Example: Reverse(&x);                                   *
//   ***********************************************************
void Reverse(Vector *result)
   {
   int m = (result->dimensions+1)/2;
   int n = result->dimensions;
   int i,j;
   if(n < 2)return;
   float *v = result->vector;
   if(n == 2) Swap(v+1,v+2);
   else
      for(i = 1,j = n;i <= m;i++,j--)
        Swap(v+i,v+j); // in UTILITY.CPP
   }

//   ***********************< Sort >****************************
//   * Purpose: Ordering a Vector                              *
//   * Description: Uses heap sort of UTILITY.CPP (Chapter 8)  *
//   * Example: Vector x(150,w); Sort(&x);                     *
//   ***********************************************************
void Sort(Vector *result)
   {
   Sort(result->dimensions,result->vector + 1); //UTILITY.cpp
   }

//   ************************< Swap >***************************
//   * Purpose: Swapping the elements of any two Vectors       *
//   * Description: Provides an efficient method of copying    *
//   *              When a Vector remains in the purpose and   *
//   *              one leaves, they swap                      *
//   * Example: Swap(&x,&y);                                   *
//   ***********************************************************
void Swap(Vector *lval,Vector *rval)
   {
   Swap(&lval->vector,&rval->vector);
   Swap(&lval->dimensions,&rval->dimensions);
   }