/*********************************************************************/
/* File:   evcoarse.cc                                               */
/* Author: Joachim Schoeberl                                         */
/* Date:   5. Nov. 2000                                             */
/*********************************************************************/

/* 
   Coarse Grid Preconditioner based on a few eigenvectors
*/

#include <multigrid.hpp>

namespace ngmg_impl
{

  EVCoarseGrid :: EVCoarseGrid (const BaseMatrix & mat,
				const Smoother & sm, 
				const ARRAY<BaseVector*> & coarsevecs,
				int maxv, double maxits)
    : smoother(sm)
  {
    /*
      int i, j, k;
      double lam;
      BaseVector & r = *mat.CreateVector();
      BaseVector & hv = *mat.CreateVector();
      ARRAY<BaseVector*> avecs;
      
      
      for (i = 1; i <= coarsevecs.Size(); i++)
      {
	vecs.Append (coarsevecs.Get(i)->Copy());
	avecs.Append (coarsevecs.Get(i)->Copy());
	
	vecs.Elem(i) -> Set (1, *coarsevecs.Get(i));
	mat.Mult (*vecs.Get(i), *avecs.Elem(i));

	double scal = sqrt (*vecs.Get(i) * *avecs.Get(i));
	*vecs.Elem(i) *= 1/scal;
	*avecs.Elem(i) *= 1/scal;
      }

    for (i = vecs.Size()+1; i <= maxv; i++)
      {
	BaseVector & w = *mat.CreateVector();
	BaseVector & aw = *mat.CreateVector();

//	if (i < coarsevecs.Size())
//	  w.Set (1, *coarsevecs.Get(i));
//	else

	w.SetRandom();

	for (j = 1; j <= 1000; j++)
	  {
	    mat.Mult (w, aw);
	    hv.SetScalar (0);
	    sm.PreSmooth (1, hv, aw, 1);
	    sm.PostSmooth (1, hv, aw, 1);

	    double waw = w * aw;
	    double wacaw = hv * aw;
	    lam = wacaw / waw;

	    cout << "i = " << i 
		 << ", j = " << j
		 << ", lam = " << wacaw / waw << endl;
	    w.Add (-1, hv);
	    
	    w *= 1/sqrt (waw);

	    //	    mat.Mult (w, aw);
	    for (k = 1; k <= vecs.Size(); k++)
	      w.Add (- (w * *avecs.Get(k)), *vecs.Get(k));
	  }

	if (lam > 1 / maxcond)
	  {
	    delete &w;
	    delete &aw;
	    break;
	  }
	mat.Mult (w, aw);
	w *= 1 / sqrt (w * aw);

	mat.Mult (w, aw);
	vecs.Append (&w);
	avecs.Append (&aw);
      }

    DenseMatrix locmat;
    locmat.SetSize(vecs.Size());
    for (i = 1; i <= vecs.Size(); i++)
      for (j = 1; j <= vecs.Size(); j++)
	locmat.Elem(i, j) = *avecs.Get(i) * *vecs.Get(j);

    if (vecs.Size())
      invlocmat = locmat.InverseMatrix();
    else
      invlocmat = NULL;

    for (i = 1; i <= avecs.Size(); i++)
      delete avecs.Elem(i);

    delete &r;
    delete &hv;
    */

    
    
    int i, j, k;
    double lam;
    
    cout.precision(8);

    maxv++;

    BaseVector & r = *mat.CreateVector();
    BaseVector & hv = *mat.CreateVector();

    BaseVector & av1 = *mat.CreateVector();
    BaseVector & av2 = *mat.CreateVector();
    BaseVector & dav1 = *mat.CreateVector();
    BaseVector & dav2 = *mat.CreateVector();

    ARRAY<BaseVector*> avecs;
    ARRAY<BaseVector*> davecs;
    DenseMatrix vadav(maxv);
    DenseMatrix vadav2(2);

    Vector lami(maxv);
    DenseMatrix evecs(maxv);

    /*
    const EdgeSmoother & esm = dynamic_cast<const EdgeSmoother&> (sm);
    BaseMatrix & pre = *esm.jac.Last();
    */
    SmoothingPreconditioner pre (sm, 1);

    EigenSystem eigen(mat, pre);
    eigen.Calc();
    double lammax = eigen.EigenValue (eigen.NumEigenValues());
    
    cout << "lammax = " << lammax << endl;

    lammax *= 1.5;  // safety factor

    ChebyshevIteration cheby(mat, pre, 20);
    cheby.SetBounds (1-1*lammax, 0.98);


    for (i = 1; i <= maxv; i++)
      {
	vecs.Append (mat.CreateVector());
	avecs.Append (mat.CreateVector());
	davecs.Append (mat.CreateVector());

	if (i <= coarsevecs.Size())
	  vecs.Elem(i)->Set (1, *coarsevecs.Get(i));
	else
	  vecs.Elem(i)->SetRandom();
      }

    // orthogonalize
    for (j = 1; j <= vecs.Size(); j++)
      {
	mat.Mult (*vecs.Elem(j), hv);
	double s = hv * (*vecs.Elem(j));
	
	(*vecs.Elem(j)) *= 1 / sqrt (s);
	hv *= 1 / sqrt(s);
	
	for (k = j+1; k <= vecs.Size(); k++)
	  {
	    double s = hv * (*vecs.Elem(k));
	    vecs.Elem(k) -> Add (-s, *vecs.Elem(j));
	  }
      }

	


    for (j = 1; j <= vecs.Size(); j++)
      mat.Mult (*vecs.Elem(j), *avecs.Elem(j));

    
    for (i = 1; i <= maxits; i++)
      {    
	cout << "i = " << i << " ";
	(*testout) << "i = " << i << " ";

	// vecs are A - orthogonal

	cout << "li = ";
	(*testout) << "li = ";
	for (j = 1; j <= maxv; j++)
	  {
	    mat.Mult (*vecs.Elem(j), hv);
	    pre.Mult (hv, r);
	    cout << (hv * r) << " ";
	    (*testout) << (hv * r) << " ";
	  }
	cout << endl;
	(*testout) << endl;

	// update upper bound
	mat.Mult (*vecs.Last(), hv);
	pre.Mult (hv, r);
	double lnext = 1 - hv * r;
	//	cout << "upper bound = " << lnext << endl;
	
	if (i >= 2)
	  cheby.SetBounds (1-1*lammax, lnext);


	for (j = 1; j <= vecs.Size(); j++)
	  {
	    /*
	    if (i % 2 == 0)
	      {
		pre.Mult (*avecs.Elem(j), hv);
		vecs.Elem(j) -> Add (-1/lammax, hv);
	      }
	    else
	    */
	      {
		cheby.Mult (*avecs.Elem(j), hv);
		vecs.Elem(j) -> Add (-1, hv);
	      }
	  }

	// V = p(I-D^{-1} A) V 
	// check V^t (A - A D^{-1} A) V

	
	for (j = 1; j <= vecs.Size(); j++)
	  mat.Mult (*vecs.Elem(j), *avecs.Elem(j));


	for (j = 1; j <= vecs.Size(); j++)
	  for (k = 1; k <= vecs.Size(); k++)
	    vadav.Elem(j,k) = - (*vecs.Elem(j) * *avecs.Elem(k));
	for (j = 1; j <= vecs.Size(); j++)
	  vadav.Elem(j,j) += 1; 

	(*testout) << "vadav = " << endl << vadav << endl;

	/*
	// do A-orthogonalization

	for (j = vecs.Size(); j >= 1; j--)
	  {
	    mat.Mult (*vecs.Elem(j), *avecs.Elem(j));
	    double s = sqrt (*vecs.Elem(j) * *avecs.Elem(j));
	    *vecs.Elem(j) /= s;
	    *avecs.Elem(j) /= s;

	    for (k = 1; k < j; k++)
	      {
		s = *avecs.Elem(j) * *vecs.Elem(k);
		vecs.Elem(k) -> Add (-1, *vecs.Elem(j));
	      }
	  }

	for (j = 1; j <= vecs.Size(); j++)
	  {
	    sm.Precond (1, *avecs.Elem(j), hv);
	    vecs.Elem(j) -> Add (-1/lammax, hv);
	  }

	for (j = 1; j <= vecs.Size(); j++)
	  for (k = 1; k <= vecs.Size(); k++)
	    vadav.Elem(j,k) = - (*vecs.Elem(j) * *avecs.Elem(k));
	for (j = 1; j <= vecs.Size(); j++)
	  vadav.Elem(j,j) += 1;
	*/
	
	vadav.EigenSystem (lami, evecs);

	(*testout) << "lami = " << lami << endl;
	(*testout) << "evecs = " << evecs << endl;

	//	if (i % 2 == 0)
	//	  cout << "lami = " << lami << endl;
	//	cout << "evecs = " << evecs << endl;

	for (j = 1; j <= vecs.Size(); j++)
	  {
	    avecs.Elem(j) -> SetScalar (0);
	    for (k = 1; k <= vecs.Size(); k++)
	      avecs.Elem(j) -> Add (evecs.Elem(j,k), *vecs.Elem(k));
	  }
	for (j = 1; j <= vecs.Size(); j++)
	  {
	    mat.Mult (*avecs.Elem(j), hv);
	    double s = sqrt (*avecs.Elem(j) * hv);
	    vecs.Elem(j) -> Set (1/s, *avecs.Elem(j));
	    avecs.Elem(j) -> Set (1/s, hv);
	  }
      }


    delete vecs.Last();
    vecs.SetSize (vecs.Size()-1);
    

    for (j = 1; j <= vecs.Size(); j++)
      mat.Mult (*vecs.Elem(j), *avecs.Elem(j));

    DenseMatrix locmat;
    locmat.SetSize(vecs.Size());
    for (i = 1; i <= vecs.Size(); i++)
      for (j = 1; j <= vecs.Size(); j++)
	locmat.Elem(i, j) = *avecs.Get(i) * *vecs.Get(j);

    if (vecs.Size())
      invlocmat = locmat.InverseMatrix();
    else
      invlocmat = NULL;

    for (i = 1; i <= avecs.Size(); i++)
      delete avecs.Elem(i);

  }
  
  
  
  void EVCoarseGrid :: Mult (const BaseVector & x, BaseVector & y) const
  {
    BaseVector & res = *smoother.CreateVector(1);
    y.SetScalar (0);
    smoother.PreSmooth (1, y, x, 1);
    smoother.Residuum (1, y, x, res);

    Vector lx(vecs.Size());
    Vector ly(vecs.Size());

    int k;
    for (k = 1; k <= vecs.Size(); k++)
      lx.Elem(k) = (res * *vecs.Get(k));

    if (invlocmat)
      invlocmat->Mult (lx, ly);
    
    for (k = 1; k <= vecs.Size(); k++)
      y.Add (ly.Get(k), *vecs.Get(k));

    smoother.PostSmooth (1, y, x, 1);

    delete &res;


    /*
    y.SetScalar (0);
    smoother.PreSmooth (1, y, x, 1);
    smoother.PostSmooth (1, y, x, 1);

    Vector lx(vecs.Size());
    Vector ly(vecs.Size());

    int k;
    for (k = 1; k <= vecs.Size(); k++)
      lx.Elem(k) = (x * *vecs.Get(k));

    if (invlocmat)
      invlocmat->Mult (lx, ly);
    
    for (k = 1; k <= vecs.Size(); k++)
      y.Add (ly.Get(k), *vecs.Get(k));
    */
  }
}