// Copyright (C) 2003-2009 Marc Duruflé
// Copyright (C) 2001-2009 Vivien Mallet
//
// This file is part of the linear-algebra library Seldon,
// http://seldon.sourceforge.net/.
//
// Seldon is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 2.1 of the License, or (at your option)
// any later version.
//
// Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
// more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Seldon. If not, see http://www.gnu.org/licenses/.


#ifndef SELDON_FILE_FUNCTIONS_ARRAYS_CXX

/*
  Functions defined in this file
  
  QuickSort(m, n, X);
  QuickSort(m, n, X, Y);
  QuickSort(m, n, X, Y, Z);

  MergeSort(m, n, X);
  MergeSort(m, n, X, Y);
  MergeSort(m, n, X, Y, Z);

  Sort(m, n, X);
  Sort(m, n, X, Y);
  Sort(m, n, X, Y, Z);
  Sort(m, X);
  Sort(m, X, Y);
  Sort(m, X, Y, Z);
  Sort(X);
  Sort(X, Y);
  Sort(X, Y, Z);

  Assemble(m, X);
  Assemble(m, X, Y);
  
  RemoveDuplicate(m, X);
  RemoveDuplicate(m, X, Y);
  RemoveDuplicate(X);
  RemoveDuplicate(X, Y);
*/

namespace Seldon
{
  
  
  ////////////
  //  SORT  //
  
  
  //! Intermediary function used for quick sort algorithm.
  template<class T, class Storage, class Allocator>
  int PartitionQuickSort(int m, int n,
			 Vector<T, Storage, Allocator>& t)
  {
    T temp, v;
    v = t(m);
    int i = m - 1;
    int j = n + 1;
    
    while (true)
      {
	do
	  {
	    j--;
	  }
	while (t(j) > v);
	
	do
	  {
	    i++;
	  }
	while (t(i) < v);
	
	if (i < j)
	  {
	    temp = t(i);
	    t(i) = t(j);
	    t(j) = temp;
	  }
	else
	  {
	    return j;
	  }
      }
  }
  
  
  //! Vector \a t is sorted by using QuickSort algorithm.
  /*!
    Sorts array \a t between position m and n.
  */
  template<class T, class Storage, class Allocator>
  void QuickSort(int m, int n,
		 Vector<T, Storage, Allocator>& t)
  {
    if (m < n)
      {
	int p = PartitionQuickSort(m, n, t);
	QuickSort(m, p, t);
	QuickSort(p+1, n, t);
      }
  }
  
  
  //! Intermediary function used for quick sort algorithm.
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  int PartitionQuickSort(int m, int n,
			 Vector<T1, Storage1, Allocator1>& t1,
			 Vector<T2, Storage2, Allocator2>& t2)
  {
    T1 temp1, v;
    T2 temp2;
    v = t1(m);
    int i = m - 1;
    int j = n + 1;
    
    while (true)
      {
	do
	  {
	    j--;
	  }
	while (t1(j) > v);
	do
	  {
	    i++;
	  }
	while (t1(i) < v);
	
	if (i < j)
	  {
	    temp1 = t1(i);
	    t1(i) = t1(j);
	    t1(j) = temp1;
	    temp2 = t2(i);
	    t2(i) = t2(j);
	    t2(j) = temp2;
	  }
	else
	  {
	    return j;
	  }
      }
  }
  
  
  //! Vector \a t1 is sorted by using QuickSort algorithm.
  /*! Sorts array \a t2 between position m and n, the sorting operation
    affects vector \a t2.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void QuickSort(int m, int n,
		 Vector<T1, Storage1, Allocator1>& t1,
		 Vector<T2, Storage2, Allocator2>& t2)
  {
    if (m < n)
      {
	int p = PartitionQuickSort(m, n, t1, t2);
	QuickSort(m, p, t1, t2);
	QuickSort(p+1, n, t1, t2);
      }
  }
  
  
  //! Intermediary function used for quick sort algorithm.
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  int PartitionQuickSort(int m, int n,
			 Vector<T1, Storage1, Allocator1>& t1,
			 Vector<T2, Storage2, Allocator2>& t2,
			 Vector<T3, Storage3, Allocator3>& t3)
  {
    T1 temp1, v;
    T2 temp2;
    T3 temp3;
    v = t1(m);
    int i = m - 1;
    int j = n + 1;
    
    while (true)
      {
	do
	  {
	    j--;
	  }
	while (t1(j) > v);
	
	do
	  {
	    i++;
	  }
	while (t1(i) < v);
	
	if (i < j)
	  {
	    temp1 = t1(i);
	    t1(i) = t1(j);
	    t1(j) = temp1;
	    temp2 = t2(i);
	    t2(i) = t2(j);
	    t2(j) = temp2;
	    temp3 = t3(i);
	    t3(i) = t3(j);
	    t3(j) = temp3;
	  }
	else
	  {
	    return j;
	  }
      }
  }
  
  
  //! Vector \a t1 is sorted by using QuickSort algorithm.
  /*! Sorts array \a t1 between position m and n, the sorting operation
    affects vectors \a t2 and \a t3.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  void QuickSort(int m, int n,
		 Vector<T1, Storage1, Allocator1>& t1,
		 Vector<T2, Storage2, Allocator2>& t2,
		 Vector<T3, Storage3, Allocator3>& t3)
  {
    if (m < n)
      {
	int p = PartitionQuickSort(m, n, t1, t2, t3);
	QuickSort(m, p, t1, t2, t3);
	QuickSort(p+1, n, t1, t2, t3);
      }
  }
  
  
  //! Vector \a tab1 is sorted by using MergeSort algorithm.
  /*!
    Sorts array \a tab1 between position m and n.
  */
  template<class T, class Storage, class Allocator>
  void MergeSort(int m, int n, Vector<T, Storage, Allocator>& tab1)
  {
    if (m <= n)
      return;
    
    int inc = 1, ind = 0, current, i, j, sup;
    Vector<T, Storage, Allocator> tab1t(m-n+1);
    // Performs a merge sort with a recurrence.
    // inc = 1, 2, 4, 8, ...
    while (inc < n)
      {
	// 'i' is the first index of the sub array of size 2*inc.
	// A loop is performed on these sub-arrays.
	// Each sub array is divided in two sub-arrays of size 'inc'.
	// These two sub-arrays are then merged.
	for (i = m; i < n - inc; i += 2 * inc)
	  {
	    ind = i;
	    current = i + inc; // Index of the second sub-array.
	    sup = i + 2 * inc; // End of the merged array.
	    if (sup >= n)
	      sup = n;
	    j = i;
	    // Loop on values of the first sub-array.
	    while (j < i + inc)
	      {
		// If the second sub-array has still unsorted elements.
		if (current < sup)
		  {
		    // Insert the elements of the second sub-array in the
		    // merged array until tab1(j) < tab1(current).
		    while (current < sup && tab1(j) > tab1(current))
		      {
			tab1t(ind-m) = tab1(current);
			current++;
			ind++;
		      }
		    // Inserts the element of the first sub-array now.
		    tab1t(ind - m) = tab1(j);
		    ind++;
		    j++;
		    // If the first sub-array is sorted, all remaining
		    // elements of the second sub-array are inserted.
		    if (j == i + inc)
		      {
			for (j = current; j < sup; j++)
			  {
			    tab1t(ind - m) = tab1(j);
			    ind++;
			  }
		      }
		  }
		else
		  {
		    // If the second sub-array is sorted, all remaining
		    // elements of the first sub-array are inserted.
		    for (current = j; current < i + inc; current++)
		      {
			tab1t(ind - m) = tab1(current);
			ind++;
		      }
		    j = current + 1;
		  }
	      }
	  }
	
	for (i = m; i < ind; i++)
	  tab1(i) = tab1t(i - m);
	
	inc = 2 * inc;
      }
  }
  
  
  //! Vector \a tab1 is sorted by using MergeSort algorithm.
  /*! Sorts array \a tab1 between position m and n. The sort operation affects
    \a tab2.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void MergeSort(int m, int n, Vector<T1, Storage1, Allocator1>& tab1,
		 Vector<T2, Storage2, Allocator2>& tab2)
  {
    if (m <= n)
      return;
    
    int inc = 1, ind = 0, current, i, j, sup;
    Vector<T1, Storage1, Allocator1> tab1t(m - n + 1);
    Vector<T2, Storage2, Allocator2> tab2t(m - n + 1);
    
    while (inc < n)
      {
	for (i = 0; i < n - inc; i += 2 * inc)
	  {
	    ind = i;
	    current = i + inc;
	    sup = i + 2 * inc;
	    if (sup >= n)
	      sup = n;
	    j = i;
	    while (j < i + inc)
	      {
		if (current < sup)
		  {
		    while (current < sup && tab1(j) > tab1(current))
		      {
			tab1t(ind - m) = tab1(current);
			tab2t(ind - m) = tab2(current);
			current++;
			ind++;
		      }
		    tab1t(ind - m) = tab1(j);
		    tab2t(ind - m) = tab2(j);
		    ind++;
		    j++;
		    if (j == i + inc)
		      {
			for (j = current; j < sup; j++)
			  {
			    tab1t(ind - m) = tab1(j);
			    tab2t(ind - m) = tab2(j);
			    ind++;
			  }
		      }
		  }
		else
		  {
		    for (current = j; current < i + inc; current++)
		      {
			tab1t(ind - m) = tab1(current);
			tab2t(ind - m) = tab2(current);
			ind++;
		      }
		    j = current + 1;
		  }
	      }
	  }
	for (i = 0; i < ind; i++)
	  {
	    tab1(i) = tab1t(i - m);
	    tab2(i) = tab2t(i - m);
	  }
	inc = 2 * inc;
      }
  }
  
  
  //! Vector \a tab1 is sorted by using MergeSort algorithm.
  /*! Sorts array \a tab1 between position m and n. The sort operation affects
    \a tab2 and \a tab3.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  void MergeSort(int m, int n, Vector<T1, Storage1, Allocator1>& tab1,
		 Vector<T2, Storage2, Allocator2>& tab2,
		 Vector<T3, Storage3, Allocator3>& tab3)
  {
    if (m <= n)
      return;
    
    int inc = 1, ind = 0, current, i, j, sup;
    Vector<T1, Storage1, Allocator1> tab1t(n - m + 1);
    Vector<T2, Storage2, Allocator2> tab2t(n - m + 1);
    Vector<T3, Storage3, Allocator3> tab3t(n - m + 1);
    
    while (inc < n)
      {
	for (i = 0; i < n - inc; i += 2 * inc)
	  {
	    ind = i;
	    current = i + inc;
	    sup = i + 2 * inc;
	    if (sup >= n)
	      sup = n;
	    j = i;
	    while (j < i + inc)
	      {
		if (current < sup)
		  {
		    while (current < sup && tab1(j) > tab1(current))
		      {
			tab1t(ind - m) = tab1(current);
			tab2t(ind - m) = tab2(current);
			tab3t(ind - m) = tab3(current);
			current++;
			ind++;
		      }
		    tab1t(ind - m) = tab1(j);
		    tab2t(ind - m) = tab2(j);
		    tab3t(ind - m) = tab3(j);
		    ind++;
		    j++;
		    if (j == i + inc)
		      {
			for (j = current; j < sup; j++)
			  {
			    tab1t(ind - m) = tab1(j);
			    tab2t(ind - m) = tab2(j);
			    tab3t(ind - m) = tab3(j);
			    ind++;
			  }
		      }
		  }
		else
		  {
		    for (current = j; current < i + inc; current++)
		      {
			tab1t(ind - m) = tab1(current);
			tab2t(ind - m) = tab2(current);
			tab3t(ind - m) = tab3(current);
			ind++;
		      }
		    j = current+1;
		  }
	      }
	  }
	for (i = 0; i < ind; i++)
	  {
	    tab1(i) = tab1t(i - m);
	    tab2(i) = tab2t(i - m);
	    tab3(i) = tab3t(i - m);
	  }
	inc = 2 * inc;
      }
  }
  
  
  //! Assembles a sparse vector.
  /*! The function sorts the indices in \a Node and adds the corresponding
    values of \a Vect corresponding.

    For example, if Node = [3 2 2 0], Vect = [1.0 0.4 0.4 -0.3], the function
    will return: n = 3, Node = [0 2 3], Vect = [-0.3 0.8 1.0].

    \param[in,out] n on entry, the number of elements to assemble; on exit,
    the number of values after assembling.
    \param[in,out] Node positions in the vector.
    \param[in,out] Vect values.
    \warning Vectors \a Node and \a Vect are not resized.
  */
  template<class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2 >
  void Assemble(int& n, Vector<int, Storage1, Allocator1>& Node,
		Vector<T2, Storage2, Allocator2>& Vect)
  {
    if (n <= 1)
      return;
    
    Sort(n, Node, Vect);
    int prec = Node(0);
    int nb = 0;
    for (int i = 1; i < n; i++)
      if (Node(i) == prec)
	{
	  Vect(nb) += Vect(i);
	}
      else
	{
	  nb++;
	  Node(nb) = Node(i);
	  Vect(nb) = Vect(i);
	  prec = Node(nb);
	}
    n = nb + 1;
  }
  
  
  //! Sorts and removes duplicate entries of a vector.
  /*!
    \param[in,out] n on entry, the number of elements to assemble; on exit,
    the number of values after assembling.
    \param[in,out] Node vector to be assembled.
    \warning The vector \a Node is not resized.
  */
  template<class T, class Storage1, class Allocator1>
  void Assemble(int& n, Vector<T, Storage1, Allocator1>& Node)
  {
    if (n <= 1)
      return;
    
    Sort(n, Node);
    T prec = Node(0);
    int nb = 1;
    for (int i = 1; i < n; i++)
      if (Node(i) != prec)
	{
	  Node(nb) = Node(i);
	  prec = Node(nb);
	  nb++;
	}
    n = nb;
  }
  

  //! Sorts and removes duplicate entries of a vector.
  template<class T, class Storage1, class Allocator1>
  void Assemble(Vector<T, Storage1, Allocator1>& Node)
  {
    int nb = Node.GetM();
    Assemble(nb, Node);
    Node.Resize(nb);
  }
  
  
  //! Sorts and removes duplicate entries of a vector.
  /*!
    Sorting operations on \a Node also affect \a Node2.
  */
  template<class T, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void RemoveDuplicate(int& n, Vector<T, Storage1, Allocator1>& Node,
		       Vector<T2, Storage2, Allocator2>& Node2)
  {
    if (n <= 1)
      return;
    
    Sort(n, Node, Node2);
    T prec = Node(0);
    int nb = 1;
    for (int i = 1; i < n; i++)
      if (Node(i) != prec)
	{
	  Node(nb) = Node(i);
	  Node2(nb) = Node2(i);
	  prec = Node(nb);
	  nb++;
	}
    n = nb;
  }
  
  
  //! Sorts and removes duplicate entries of a vector.
  template<class T, class Storage1, class Allocator1>
  void RemoveDuplicate(int& n, Vector<T, Storage1, Allocator1>& Node)
  {
    Assemble(n, Node);
  }
  
  
  //! Sorts and removes duplicate entries of a vector.
  /*!
    Sorting operations of \a Node also affect \a Node2.
  */
  template<class T, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void RemoveDuplicate(Vector<T, Storage1, Allocator1>& Node,
		       Vector<T2, Storage2, Allocator2>& Node2)
  {
    int n = Node.GetM();
    if (n <= 1)
      return;
    
    RemoveDuplicate(n, Node, Node2);
    Node.Resize(n);
    Node2.Resize(n);
  }

  
  //! Sorts and removes duplicate entries of a vector.
  template<class T, class Storage1, class Allocator1>
  void RemoveDuplicate(Vector<T, Storage1, Allocator1>& Node)
  {
    int n = Node.GetM();
    if (n <= 1)
      return;
    
    Assemble(n, Node);
    Node.Resize(n);
  }
  
  
  //! Sorts vector \a V between a start position and an end position.
  template<class T, class Storage, class Allocator>
  void Sort(int m, int n, Vector<T, Storage, Allocator>& V)
  {
    QuickSort(m, n, V);
  }
  
  
  //! Sorts vector \a V between a start position and an end position.
  /*!
    The sorting operation of \a V also affects \a V2.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void Sort(int m, int n, Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2)
  {
    QuickSort(m, n, V, V2);
  }
  
  
  //! Sorts vector \a V between a start position and an end position.
  /*!
    The sorting operation of \a V also affects \a V2 and \a V3.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  void Sort(int m, int n, Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2,
	    Vector<T3, Storage3, Allocator3>& V3)
  {
    QuickSort(m, n, V, V2, V3);
  }
  
  
  //! Sorts the \a n first elements of \a V.
  template<class T, class Storage, class Allocator>
  void Sort(int n, Vector<T, Storage, Allocator>& V)
  {
    Sort(0, n - 1, V);
  }
  
  
  //! Sorts the \a n first elements of \a V.
  /*!
    The sorting operation of \a V also affects \a V2.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void Sort(int n, Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2)
  {
    Sort(0, n - 1, V, V2);
  }
  
  
  //! Sorts the \a n first elements of \a V.
  /*!
    The sorting operation of \a V also affects \a V2 and \a V3.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  void Sort(int n, Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2,
	    Vector<T3, Storage3, Allocator3>& V3)
  {
    Sort(0, n - 1, V, V2, V3);
  }
  
  
  //! Sorts vector \a V.
  template<class T, class Storage, class Allocator>
  void Sort(Vector<T, Storage, Allocator>& V)
  {
    Sort(0, V.GetM() - 1, V);
  }
  
  
  //! Sorts vector \a V.
  /*!
    The sorting operation of \a V also affects \a V2.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2>
  void Sort(Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2)
  {
    Sort(0, V.GetM() - 1, V, V2);
  }
  
  
  //! Sorts vector \a V.
  /*!
    The sorting operation of \a V also affects \a V2 and \a V3.
  */
  template<class T1, class Storage1, class Allocator1,
	   class T2, class Storage2, class Allocator2,
	   class T3, class Storage3, class Allocator3>
  void Sort(Vector<T1, Storage1, Allocator1>& V,
	    Vector<T2, Storage2, Allocator2>& V2,
	    Vector<T3, Storage3, Allocator3>& V3)
  {
    Sort(0, V.GetM() - 1, V, V2, V3);
  }
  
  
  //  SORT  //
  ////////////
  
  
} // namespace Seldon

#define SELDON_FILE_FUNCTIONS_ARRAYS_CXX
#endif