/* *************************************************************************/
/* File:   order.cc                                                        */
/* Author: Joachim Schoeberl                                               */
/* Date:   25. Mar. 2000                                                   */
/* *************************************************************************/



#include <la.hpp>
namespace ngla
{
  using namespace ngla;


  /* 

  Finds ordering for sparse Cholesky Factorization
  Algorithm: Minimum Degree
              
  See:       The evolution of the 
  minimum degree ordering algorithm

  A. George and J.W.H. Liu
  SIAM Review, Vol31, 1989, pp1-19

  */


  MDOVertex :: MDOVertex(int ma)
  {
    Init (ma);
  }

  void MDOVertex :: Init (int ma)
  {
    master = ma;
    nextslave = -1;
    flag = 0;
    eliminated = 0;
    used = 0;
  }


  MDOVertex :: ~MDOVertex()
  {
    ;
  }


  ngstd::BlockAllocator CliqueEl :: ball(sizeof (CliqueEl));

  void *  CliqueEl::operator new(size_t)
  {
    return ball.Alloc();
  }

  void CliqueEl :: operator delete (void * p)
  {
    ball.Free (p);
  }


  /*  
  CliqueEl :: CliqueEl ()
  {
    next = NULL;
    nextcl = NULL;
    eliminate = 0;
    flag = 0;
  }
  */




  MinimumDegreeOrdering :: MinimumDegreeOrdering (int an)
    :  priqueue(an, an+1)
  {
    int i;

    n = an;

    order.SetSize(n);
    blocknr.SetSize(n);
    cliques.SetSize (n);
    vertices.SetSize (n);

    for (i = 0; i < n; i++)
      cliques[i] = NULL;
    for (i = 0; i < n; i++)
      blocknr[i] = 0;
    for (i = 0; i < n; i++)
      order[i] = 0;
    for (i = 0; i < n; i++)
      vertices[i].Init(i);
  }



  void MinimumDegreeOrdering :: AddEdge (int v1, int v2)
  {
    if (v1 == v2) return;

    CliqueEl *p1, *p2;
    p1 = new CliqueEl;
    p2 = new CliqueEl;
  
    p1->vnr = v1;
    p2->vnr = v2;

    p1->eliminate = 0;
    p2->eliminate = 0;

    p1->next = p2;
    p2->next = p1;

    p1->flag = 0;
    p2->flag = 0;

    p1->nextcl = cliques[v1];
    cliques[v1] = p1;
    p2->nextcl = cliques[v2];
    cliques[v2] = p2;
  }


  void MinimumDegreeOrdering :: PrintCliques () 
  {
    CliqueEl *p1, *p2;

    for (int i = 0; i < n; i++)
      if (!vertices[i].Eliminated())
	{
	  (cout) << "Vertex " << i << ", degree = " 
		 << CalcDegree (i) 
	    //	       << " precaled: " << vertices.Get(i).Degree() 
		 << endl;
	  p1 = cliques[i];
	  while (p1)
	    {
	      p2 = p1;
	      (cout) << "( ";
	      do
		{
		  if (!vertices[p2->GetVertexNr()].Eliminated())
		    (cout) << p2->GetVertexNr() << " ";
		  p2 = p2->next;
		}
	      while (p2 != p1);
	      (cout) << ")";
	      p1 = p1->nextcl;
	    }
	  (cout) << endl;
	}
  }





  int MinimumDegreeOrdering :: CalcDegree (int v1)
  {
    CliqueEl *p1, *p2;
    int deg = 0;

    // count members of all cliques of vertex


    // clear flags to count vertices just once:
    p1 = cliques[v1];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    vertices[p2->GetVertexNr()].SetUsed(0);   
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }

    p1 = cliques[v1];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    if (!vertices[p2->GetVertexNr()].Used())
	      {
		deg++;
		vertices[p2->GetVertexNr()].SetUsed(1);
	      }
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }
  
    return deg;
  }
 


  void MinimumDegreeOrdering :: EliminateVertex (int v, bool master)
  {
    //  (cout) << "eliminate vertex " << v << endl;
    //  PrintCliques();

    CliqueEl *p1, *p2, *p3, *p4, *newp;

  
    // clear used-field
    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    vertices[p2->GetVertexNr()].SetUsed (0);
	    p2->eliminate = 1;
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }

  
    //  new clique is union of cliques of vertex v
    newp = NULL;
    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    if (!vertices[p2->GetVertexNr()].Used() && p2->GetVertexNr() != v)
	      {
		p3 = new CliqueEl;
		p3 -> next = newp;
		p3 -> vnr = p2->GetVertexNr();
		p3 -> eliminate = 0;
	      
		p3 -> nextcl = NULL;
		p3 -> flag = 0;

		newp = p3;

		vertices[p2->GetVertexNr()].SetUsed(1);
	      }
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }

    // close new clique
    if (newp)
      {
	p3 = newp;
	while (p3->next)
	  p3 = p3->next;
	p3->next = newp;
      }


    // find dominated cliques
    if (newp && master)
      {
	// set flag for new clique
	p3 = newp;
	do
	  {
	    vertices[p3->GetVertexNr()].SetFlag (1);
	    p3 = p3->next;
	  }
	while (p3 != newp);
      

	// check all cliques of all members of new clique
	p3 = newp;
	do
	  {
	    p1 = cliques[p3->GetVertexNr()];
	    while (p1)
	      {
		if (!p1->eliminate)
		  {
		    int dominated = 1;
		  
		    p2 = p1;
		    do
		      {
			if (!vertices[p2->GetVertexNr()].Flag())
			  {
			    dominated = 0;
			    break;
			  }
			p2 = p2->next;
		      }
		    while (p1 != p2);

		    // if clique is dominated, set eliminate flag
		    if (dominated)
		      {
			p2 = p1;
			do
			  {
			    p2->eliminate = 1;
			    p2 = p2->next;
			  }
			while (p1 != p2);
		      }
		  } 
		p1 = p1->nextcl;
	      }
	    p3 = p3->next;
	  }
	while (p3 != newp);
     
	// clear flag
	p3 = newp;
	do
	  {
	    vertices[p3->GetVertexNr()].SetFlag (0);
	    p3 = p3->next;
	  }
	while (p3 != newp);
      }
  
    // delete old cliques
    if (newp)
      {
	p3 = newp;

	do
	  {
	    // p3->next is first remaining clique:

	    p3->nextcl = cliques[p3->GetVertexNr()];
	    while (p3->nextcl && p3->nextcl->eliminate)
	      p3->nextcl = p3->nextcl->nextcl;


	    CliqueEl hcel;
	    hcel.nextcl = cliques[p3->GetVertexNr()];
	    p1 = &hcel;
	  
	    while (p1)
	      {
		while (p1->nextcl && p1->nextcl->eliminate)
		  {
		    CliqueEl * hp = p1->nextcl;
		    p1->nextcl = p1->nextcl->nextcl;
		    delete hp;
		  }
		p1 = p1->nextcl;
	      }

	    cliques[p3->GetVertexNr()] = p3;

	    p3 = p3->next;
	  }
	while (p3 != newp);
      }
  
    p1 = cliques[v];
    while (p1)
      {
	p2 = p1->GetNextClique();
	delete p1;
	p1 = p2;
      }



    // find connections
    if (newp && master)
      {
	p3 = newp;

	do
	  {
	    int nclp3;

	    if (vertices[p3->GetVertexNr()].Master() == p3->GetVertexNr())
	      {
		SetFlagCliques (p3->GetVertexNr());
		nclp3 = NumCliques (p3->GetVertexNr());
	      
		p4 = p3->next;
		while (p4 != newp)
		  {
		    // are p3 and p4 connected ?
		  
		    if (vertices[p4->GetVertexNr()].Master() == p4->GetVertexNr())
		      {
			bool samecl = 1;
		      
			if (NumCliques (p4->GetVertexNr()) == nclp3)
			  {
			    p1 = p4;
			    while (p1)
			      {
				if (!p1->flag)
				  samecl = 0;
				p1 = p1->nextcl;
			      }
			  }
			else
			  samecl = 0;
		      
			if (samecl)
			  {
			    // append p4 and its slave to slave list of p3

			    int hv = p3->GetVertexNr();
			    while (vertices[hv].NextSlave() != -1)
			      hv = vertices[hv].NextSlave();
			    
			    vertices[hv].SetNextSlave (p4->GetVertexNr());
			    while (hv != -1)
			      {
				vertices[hv].SetMaster (p3->GetVertexNr());
				hv = vertices[hv].NextSlave();
			      }
			  }
		      }

		  
		    p4 = p4->next;
		  }
		ClearFlagCliques (p3->GetVertexNr());
	      }

	    p3 = p3->next;
	  }
	while (p3 != newp);
      }

    if (newp)
      {
	if (master)
	  { // eliminated node was master
	  
	    // calc master degrees in new clique
	    p3 = newp;
	    do
	      {
		if (vertices[p3->GetVertexNr()].Master() == p3->GetVertexNr())
		  {
		    int deg = CalcDegree (p3->GetVertexNr());
		    priqueue.SetDegree (p3->GetVertexNr(), deg);
		  }
		p3 = p3->next;
	      }
	    while (p3 != newp);      
	  
	    // copy slave degrees
	    p3 = newp;
	    do
	      {
		int deg = priqueue.GetDegree (vertices[p3->GetVertexNr()].Master());
		priqueue.SetDegree (p3->GetVertexNr(), deg);
		p3 = p3->next;
	      }
	    while (p3 != newp);      
	  }
	else
	  {
	    // eliminated node was slave

	    p3 = newp;
	    do
	      {
		int deg = priqueue.GetDegree (p3->GetVertexNr())-1;
		priqueue.SetDegree (p3->GetVertexNr(), deg);
		p3 = p3->next;
	      }
	    while (p3 != newp); 
	  }
      }

    int cnt = 0;

    if (newp)
      {
	p3 = newp;
	do 
	  {
	    cnt++;
	    p3 = p3->next;
	  }
	while (p3 != newp);
      
	vertices[v].nconnected = cnt;
	vertices[v].connected = new int[cnt];
	cnt = 0;
      
	do 
	  {
	    vertices[v].connected[cnt] = p3->vnr;
	    cnt++;
	    p3 = p3->next;
	  }
	while (p3 != newp);
      }
    else
      {
	vertices[v].nconnected = 0;
	vertices[v].connected = NULL;
      }

    //  cliques.Elem(v) = savep;
    cliques[v] = NULL;
  }



















  void MinimumDegreeOrdering :: Order()
  {
    bool master;

    for (int j = 0; j < n; j++)
      priqueue.SetDegree(j, CalcDegree(j));


    int minj = -1;
    int lastel = -1;

    if (n > 5000)
      cout << "order " << flush;

    for (int i = 0; i < n; i++)
      {
	if (n > 5000 && i % 1000 == 999)
	  {
	    if (i % 10000 == 9999)
	      cout << "+" << flush;
	    else
	      cout << "." << flush;
	  }
      
	if (lastel != -1 && vertices[lastel].NextSlave() != -1)
	  {
	    master = 0;
	    minj = vertices[lastel].NextSlave();

	    if (vertices[minj].Eliminated())
	      cerr << "alread eliminated !!!" << endl;
	    priqueue.Invalidate(minj);
	  }
	else
	  {
	    master = 1;
	    do
	      {
		minj = priqueue.MinDegree();
		priqueue.Invalidate(minj); 
		if (vertices[minj].Master() != minj)
		  priqueue.SetDegree (minj, n);
	      }
	    while (vertices[minj].Master() != minj);
	  }

	order[i] = minj;

	if (master)
	  blocknr[i] = i;
	else
	  blocknr[i] = blocknr[i-1];

	//      (*testout) << "d(" << i << ") = " << blocknr.Get(i) << endl;

	EliminateVertex (minj, master);
	vertices[minj].SetEliminated (1);
	lastel = minj;
      }
    // PrintCliques();
  }



  MinimumDegreeOrdering:: ~MinimumDegreeOrdering ()
  {
    for (int i = 0; i < vertices.Size(); i++)
      delete [] vertices[i].connected;
  }




  int MinimumDegreeOrdering :: NumCliques (int v) const
  {
    CliqueEl *p1;
    int cnt = 0;

    p1 = cliques[v];
    while (p1)
      {
	p1 = p1->nextcl;
	cnt++;
      }
    return cnt;
  }


  void MinimumDegreeOrdering :: SetFlagNodes (int v)
  {
    CliqueEl *p1, *p2;

    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    vertices[p2->GetVertexNr()].SetFlag (1);
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }
  }

  void MinimumDegreeOrdering :: ClearFlagNodes (int v)
  {
    CliqueEl *p1, *p2;

    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    vertices[p2->GetVertexNr()].SetFlag (0);
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }
  }

  void MinimumDegreeOrdering :: SetFlagCliques (int v)
  {
    CliqueEl *p1, *p2;

    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    p2->flag = 1;
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }
  }

  void MinimumDegreeOrdering :: ClearFlagCliques (int v)
  {
    CliqueEl *p1, *p2;

    p1 = cliques[v];
    while (p1)
      {
	p2 = p1;
	do
	  {
	    p2->flag = 0;
	    p2 = p2->next;
	  }
	while (p2 != p1);
	p1 = p1->nextcl;
      }
  }


  int MinimumDegreeOrdering :: GetNZE() const
  {
    int cnt = 0;
    for (int i = 0; i < n; i++)
      cnt += vertices[i].nconnected;
    return cnt;
  }
    



  /*
    struct entry { 
    int degree, prev, next;
    }
    ARRAY<entry> list;
    ARRAY<int> first_in_class;
  */

  MDOPriorityQueue :: MDOPriorityQueue (int size, int maxdeg)
  {
    int i;
    list.SetSize (size);
    for (i = 0; i < size; i++)
      list[i].degree = 0;
    first_in_class.SetSize (maxdeg);
    for (i = 0; i < maxdeg; i++)
      first_in_class[i] = -1;
    //  (*testout) << "priqueue, size = " << size << endl;
  }

  MDOPriorityQueue :: ~MDOPriorityQueue ()
  {
    ;
  }

  int MDOPriorityQueue :: MinDegree () const
  {
    for (int i = 0; i < first_in_class.Size(); i++)
      {
	if (first_in_class[i] != -1)
	  return first_in_class[i];
      }
    return 0;
  }

  void MDOPriorityQueue :: SetDegree (int nr, int deg)
  {
    //    cout << "set deg ( " << nr << " ) : " << deg << endl;
    if (deg == 0)
      deg++;

    if (list[nr].degree > 0)
      Invalidate (nr);
    //  (*testout) << "set nr = " << nr << ", deg = " << deg << endl;
    if (first_in_class[deg] != -1)
      {
	int next = first_in_class[deg];
	int prev = list[next].prev;

	list[nr].next = next;
	list[nr].prev = prev;
	list[next].prev = nr;
	list[prev].next = nr;
      }
    else
      {
	list[nr].next = nr;
	list[nr].prev = nr;
	first_in_class[deg] = nr;
      }
    list[nr].degree = deg;
  }

  void MDOPriorityQueue :: Invalidate (int nr)
  {
    if (!list[nr].degree)
      cerr << "already eliminated" << endl;

    if (list[nr].next == nr)
      { // just on element in class
	first_in_class[list[nr].degree] = -1;
	//      (*testout) << "delete el last " << nr << endl;
      }
    else
      {
	int next = list[nr].next;
	int prev = list[nr].prev;
	/*
	  (*testout) << "delete el alt, nr = " << nr
	  << ", deg = " << list.Elem(nr).degree 
	  << ", first_in_cl = " << first_in_class.Elem(list.Elem(nr).degree)							      
	  << endl;
	*/
      
	list[prev].next = next;
	list[next].prev = prev;
	first_in_class[list[nr].degree] = next;
      }
    list[nr].degree = 0;
  }




}