LAP_Vector.hpp

Go to the documentation of this file.

#ifndef LAP_Vector_HPP
#define LAP_Vector_HPP


using namespace std;



template<class T>
const tFlag LAP_Vector<T>::ASCENT=1;
template<class T>
const tFlag LAP_Vector<T>::DESCENT=2;

template<class T>
T LAP_Vector<T>::maxValue(tLVectorIndex& i) const {
  tLVectorIndex n=getSize();
  tLVectorIncrement incx =getIncrement();
  const T* v=&(*this)(0);
  tLVectorIndex k=0;
  if (n==0) return 0;
  T maxValue=*v;
  i=0;
  v++;
  for (k=1;k<n;k++) {
    if (*v>maxValue) {
      maxValue=*v;
      i=k;
    }   
    v+=incx;
  }
  return maxValue;
}

template<class T>
T LAP_Vector<T>::minValue(tLVectorIndex& i) const {
  tLVectorIndex n=getSize();
  tLVectorIncrement incx =getIncrement();
  const T* v=&(*this)(0);
  tLVectorIndex k=0;
  if (n==0) return 0;
  T minValue=*v;
  i=0;
  v++;
  for (k=1;k<n;k++) {
    if (*v<minValue) {
      minValue=*v;
      i=k;
    }   
    v+=incx;
  }
  return minValue;
}


template<class T>
tLVectorIndex LAP_Vector<T>::indexMin() const {
  tLVectorIndex n = getSize();
  tLVectorIndex sx=getStart();
  tLVectorIncrement incx = getIncrement();
  tLVectorIndex idamin=0;
  tLVectorIndex i=0,index=sx;
  T v,vmin;
  if ((n<=1) || (incx==0)) return idamin;
  if (incx==1) {
    vmin=fabs((*this)[index++]);
    for (i=1;i<n;i++) {
      v=fabs((*this)[index]);
      if (v<=vmin) {
        idamin=i;
        vmin=v;
      }
      index++;
      
    }
  } else {
    vmin=fabs((*this)[index]);
    index+=incx;
    for (i=1;i<n;i++) {
      v=fabs((*this)[index]);
      if (v<=vmin) {
        idamin=i;
        vmin=v;
      }
      index+=incx;
      
    }
  }
  return idamin;
}


template<class T>
void LAP_Vector<T>::merge(const T &alpha,const LAP_Vector<T>& x,
                          const T &beta,const LAP_Vector<T>& y) {
    tLVectorIndex i,n=x.getSize();
    tLVectorIndex p=y.getSize();
    setSize(n+p);
    
    T *vs=&(*this)(0);
    const T *xs=&x(0);
    const T *ys=&y(0);
    
    tLVectorIndex inc=getIncrement();
    tLVectorIncrement incx=x.getIncrement();
    tLVectorIncrement incy=y.getIncrement();
    
    for (i=0;i<n-1;i++) {
        *vs=*xs*alpha;
        vs+=inc;
        xs+=incx;
    }
    //i=n
    *vs=*xs*alpha;
    if (p>0) {
        vs+=inc;
        for (i=0;i<p-1;i++) {
            *vs=*ys*beta;
            vs+=inc;
            ys+=incy;
        }
        *vs=*ys*beta;
    }
       
}


template<class T>
tReal LAP_Vector<T>::distance2(const LAP_Vector<T>&  a,const LAP_Vector<T>& b) {
    if (a.getSize()!=b.getSize())
      throw LAP_Exception("math/linalg/core",
                          "LAP_Vector::distance2",
                          "incompatible dimension of vectors");
    tReal  scale,ssq,distance=0;
    tLVectorIndex n=a.getSize();
    tLVectorIndex sx=a.getStart();
    tLVectorIndex sy=b.getStart();
    tLVectorIncrement incx=a.getIncrement();
    tLVectorIncrement incy=b.getIncrement();
    tLVectorIndex i,ix,iy;
    double d,absd;
    if ((incx==0) || incy==0) distance=0;
    else if (n==1) distance=fabs(a[sx]-b[sy]);
    else {
      scale=0;
      ssq=1;
      ix=sx;
      iy=sy;
      for (i=0;i<n;i++) {
        d=a[ix]-b[iy];
        if (d!=0) {
          absd=fabs(d);
          if (scale <absd) {
            ssq=1+ssq*(scale*scale/(absd*absd));
          } else {
            ssq+=(absd*absd)/(scale*scale);
          }
        }
        ix+=incx;
        iy+=incy;
      }
      distance=scale*sqrt(ssq);
    }
    return distance;
}

template<class T>
T LAP_Vector<T>::sum() const {
  
  tLVectorIndex n=getSize();
  if (n==0) return 0;
  tLVectorIndex sx=getStart();
  tLVectorIncrement incx=getIncrement();
  T s=0;
  const T *iter=&(*this)[sx];
  const T *end=&(*this)[sx+n-1];
  if (incx==1) {
    while( iter<=end) {
      s+=*iter;
      iter++;
    }
  } else {
    while( iter<=end) {
      s+=*iter;
      iter+=incx;
    }
  }
  return s;
}


template<class T>
void LAP_Vector<T>::sort(const tFlag& order) {

  tLVectorIndex h=1;
  tLVectorIndex i,j;
  T v;
  tLVectorIndex N=getSize();
  T *values=&(*this)(0);
  tLVectorIncrement inc=getIncrement();
  
  do {
      h=3*h+1; 
  } while (h<=N);
  
  do {
    h/=3;
    for (i=h;i<N;i++) {
      if (compare(values[i*inc],
                  values[(i-h)*inc],
                  order)) {
        v=values[i*inc];
        j=i;
        do {
          j=j*inc;
          values[j*inc]=values[(j-h)*inc];
          j=j-h;
        } while ((j>=h) &&
                 (!compare(values[(j-h)*inc],
                           v,
                           order)));
        (*this)[j]=v;
      }
    }
  }  while (h>1);
}

template<class T>
void LAP_Vector<T>::quickSort() {
  int n=getSize();
  if (n<=1) return;
  quickSort(&(*this)(0),getIncrement(),0,n-1);
}

template<class T>
void LAP_Vector<T>::quickSort(T *v,
                              const tLVectorIncrement& incv,
                              const tLVectorIndex& from,
                              const tLVectorIndex& to) {
  // no need to sort a vector of zero or one elements
  if (from>=to) return;
  
  
  // select the pivot value
  tLVectorIndex pivotLLVectorIndex = (from+to) / 2;
  T pivotValue=v[pivotLLVectorIndex*incv];
  // partition the vector
  tLVectorIndex left=from;
  tLVectorIndex right=to;
  while (left<=right) {
    while ((left<=to) && (v[left*incv]>pivotValue)) left++;
    while ((right>from) && (v[right*incv]<pivotValue)) right--;
    if (left<=right) {
      std::swap(v[left*incv],v[right*incv]);
      left++;
      if (right>0) right--;
    }
  }
  
  // sort the two sub vectors
  quickSort(v, incv,from, right);
  quickSort(v, incv,left, to);
}

template<class T>
void LAP_Vector<T>::reverse() {

    tLVectorIndex n =getSize();
    tLVectorIndex sx= getStart();
    tLVectorIncrement incx = getIncrement();
    tLVectorIndex i=sx,k=sx+(n-1)*incx;
    T temp;
    //i=sx+i*incx
    while (i<k) {
        temp=(*this)[k];
        (*this)[k]=(*this)[i];
        (*this)[i]=temp;
        i+=incx;
        k-=incx;
    }
}

#endif