#include <comp.hpp>

#include <multigrid.hpp>
using namespace ngmg;

#include <solve.hpp>


namespace ngcomp
{
  using namespace ngcomp;

  Preconditioner :: Preconditioner ()
  {
    test = 0;
    timing = 0;
    print = 0;
  }
  
  Preconditioner :: Preconditioner (const Flags & flags)
  {
    test = flags.GetDefineFlag ("test");
    timing = flags.GetDefineFlag ("timing");
    print = flags.GetDefineFlag ("print");
  }
  
  Preconditioner :: ~Preconditioner ()
  {
    ;
  }

  void Preconditioner :: Test () const
  {
    cout << "Compute eigenvalues" << endl;
    const BaseMatrix & amat = GetAMatrix();
    const BaseMatrix & pre = GetMatrix();

    EigenSystem eigen (amat, pre);
    eigen.SetPrecision(1e-30);
    eigen.SetMaxSteps(1000); 
    eigen.Calc();
    eigen.PrintEigenValues (cout);
    (*testout) << " Min Eigenvalue : "  << eigen.EigenValue(1) << endl; 
    (*testout) << " Max Eigenvalue : " << eigen.MaxEigenValue() << endl; 
    (*testout) << " Condition   " << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl; 
  }



  void Preconditioner :: Timing () const
  {
    cout << "Timing Preconditioner ... " << flush;
    const BaseMatrix & amat = GetAMatrix();
    const BaseMatrix & pre = GetMatrix();

    clock_t starttime;
    double time;
    starttime = clock();

    BaseVector & vecf = *pre.CreateVector();
    BaseVector & vecu = *pre.CreateVector();

    vecf = 1;
    int steps = 0;
    do
      {
	vecu = pre * vecf;
	steps++;
	time = double(clock() - starttime) / CLOCKS_PER_SEC;
      }
    while (time < 2.0);

    cout << " 1 step takes " 
	 << time / steps
	 << " seconds" << endl;



    starttime = clock();
    steps = 0;
    do
      {
	vecu = amat * vecf;
	steps++;
	time = double(clock() - starttime) / CLOCKS_PER_SEC;
      }
    while (time < 2.0);

    cout << ", 1 matrix takes " 
	 << time / steps
	 << " seconds" << endl;

  }





  MGPreconditioner :: MGPreconditioner (PDE * pde, Flags & flags)
    : Preconditioner (flags)
  {
    const MeshAccess & ma = pde->GetMeshAccess();
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));

    const LinearForm * lfconstraint = 
      pde->GetLinearForm (flags.GetStringFlag ("constraint", ""),1);
    const FESpace & fes = bfa->GetFESpace();


    const BilinearForm * lo_bfa = bfa;
    const FESpace * lo_fes = &fes;

    if (&bfa->GetLowOrderBilinearForm())
      {
	lo_bfa = &bfa->GetLowOrderBilinearForm();
	lo_fes = &fes.LowOrderFESpace();
      }
    
    Smoother * sm = NULL;
    //    const char * smoother = flags.GetStringFlag ("smoother", "point");
    smoothertype = flags.GetStringFlag ("smoother", "point");

    if (smoothertype == "point")
      {
	sm = new GSSmoother (ma, *lo_bfa);
      }
    else if (smoothertype == "line")
      {
	sm = new AnisotropicSmoother (ma, *lo_bfa);
      }
    else if (smoothertype == "block") 
      {
	if (!lfconstraint)
	  sm = new BlockSmoother (ma, *lo_bfa);
	else
	  sm = new BlockSmoother (ma, *lo_bfa, *lfconstraint);
      }
    else if (smoothertype == "potential")
      {
	sm = new PotentialSmoother (ma, *lo_bfa);
      }

    else
      cerr << "Unknown Smoother " << smoothertype << endl;


    const Prolongation * prol = lo_fes->GetProlongation();

    mgp = new MultigridPreconditioner (ma, *lo_fes, *lo_bfa, sm, 
				       const_cast<Prolongation*> (prol));

    mgp->SetOwnProlongation (0);
    mgp->SetSmoothingSteps (int(flags.GetNumFlag ("smoothingsteps", 1)));
    mgp->SetCycle (int(flags.GetNumFlag ("cycle", 1)));
    mgp->SetIncreaseSmoothingSteps (int(flags.GetNumFlag ("increasesmoothingsteps", 1)));
    mgp->SetCoarseSmoothingSteps (int(flags.GetNumFlag ("coarsesmoothingsteps", 1)));

    MultigridPreconditioner::COARSETYPE ct = MultigridPreconditioner::EXACT_COARSE;
    const char * coarse = flags.GetStringFlag ("coarsetype", "direct");
    if (strcmp (coarse, "smoothing") == 0)
      ct = MultigridPreconditioner::SMOOTHING_COARSE;
    else if (strcmp (coarse, "cg") == 0)
      ct = MultigridPreconditioner::CG_COARSE;
    
    mgp->SetCoarseType (ct);
    

    coarse_pre = 
      pde->GetPreconditioner (flags.GetStringFlag ("coarseprecond", ""), 1);
    if (coarse_pre)
      mgp->SetCoarseType (MultigridPreconditioner::USER_COARSE);

    finesmoothingsteps = int (flags.GetNumFlag ("finesmoothingsteps", 1));

    tlp = 0;
  }


 
  MGPreconditioner :: ~MGPreconditioner()
  {
    delete mgp;
    delete tlp;
  }
  

  void MGPreconditioner :: Update ()
  {
    mgp->Update();

    if (coarse_pre)
      mgp->SetCoarseGridPreconditioner (&coarse_pre->GetMatrix());

    if (&bfa->GetLowOrderBilinearForm())
      {
	delete tlp;

	cout << "now prepare 2-level" << endl;

	Smoother * fine_smoother;
	if (smoothertype == "potential")
	  fine_smoother = new PotentialSmoother (bfa->GetMeshAccess(), *bfa);
	else
	  fine_smoother = new BlockSmoother (bfa->GetMeshAccess(), *bfa);
	  // fine_smoother = new GSSmoother (bfa->GetMeshAccess(), *bfa);

	tlp = new TwoLevelMatrix (&bfa->GetMatrix(),
				  mgp,
				  fine_smoother,
				  bfa->GetMeshAccess().GetNLevels()-1);
	
	tlp -> SetSmoothingSteps (finesmoothingsteps);
	tlp -> Update();
      }
    else
      tlp = 0;

    if (timing) Timing();
    if (test) Test();
  }
  
  const ngla::BaseMatrix & MGPreconditioner :: GetMatrix() const
  {
    if (tlp)
      return *tlp;
    else
      return *mgp; 
  }

  void MGPreconditioner::PrintReport (ostream & ost)
  {
    ost << "Multigrid preconditioner" << endl
	<< "bilinear-form = " << bfa->GetName() << endl
	<< "smoothertype = " << smoothertype << endl;
  }


  void MGPreconditioner::MemoryUsage (ARRAY<MemoryUsageStruct*> & mu) const
  {
    int olds = mu.Size();

    if (&GetMatrix()) GetMatrix().MemoryUsage (mu);;

    for (int i = olds; i < mu.Size(); i++)
      mu[i]->AddName (string(" mgpre ")); // +GetName());
  }
    








  DirectPreconditioner :: DirectPreconditioner (PDE * pde, Flags & flags)
    : Preconditioner(flags)
  {
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    inverse = NULL;
  }

  DirectPreconditioner :: ~DirectPreconditioner()
  {
    delete inverse;
  }

  void DirectPreconditioner :: Update ()
  {
    delete inverse;
    try
      {
	inverse = dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix())
	  .InverseMatrix();
	if (print)
	  (*testout) << "inverse = " << endl << (*inverse) << endl;
      }
    catch (exception & e)
      {
	throw Exception ("DirectPreconditioner: needs a sparse matrix");
      }
    // (*testout) << "mat = " << bfa->GetMatrix() << endl;
    // (*testout) << "inv = " << (*inverse) << endl;

    //    if (test) Test();
  }

  const ngla::BaseMatrix & DirectPreconditioner :: GetMatrix() const
  {
    return *inverse;
  }







  DNDDPreconditioner :: DNDDPreconditioner (PDE * apde, Flags & flags)
  {
    pde = apde;
    bfa = apde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    inverse = NULL;
  }

  DNDDPreconditioner :: ~DNDDPreconditioner()
  {
    delete inverse;
  }

  void DNDDPreconditioner :: Update ()
  {
    cout << "Dirichlet-Neumann DD Preconditioner" << endl;

    int i, j;

    const MeshAccess & ma = pde->GetMeshAccess();

    int np = bfa->GetFESpace().GetNDof();
    int ne = ma.GetNE();
    ARRAY<int> domain(np), dnums;
    const FESpace & space = bfa->GetFESpace();

    for (i = 0; i < np; i++)
      domain[i] = 1;

    for (i = 0; i < ne; i++)
      {
	space.GetDofNrs (i, dnums);
	int elind = ma.GetElIndex(i);
	if (elind == 2 || elind == 4)
	  for (j = 0; j < dnums.Size(); j++)
	    domain[dnums[j]] = 2;
      }

    delete inverse;
    //  inverse = bfa->GetMatrix().InverseMatrix();
    // inverse = new SparseSystemLDL<SysMatrixC1d,SysVectorC1d> (dynamic_cast<const SparseSystemMatrix<SysMatrixC1d,SysVectorC1d>&> (bfa->GetMatrix()), NULL, &domain);
  }

  const ngla::BaseMatrix & DNDDPreconditioner :: GetMatrix() const
  {
    return *inverse;
  }







  LocalPreconditioner :: LocalPreconditioner (PDE * pde, Flags & flags)
    : Preconditioner (flags)
  {
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    block = flags.GetDefineFlag ("block");
    jacobi = NULL;
  }

  LocalPreconditioner :: ~LocalPreconditioner()
  {
    delete jacobi;
  }

  void LocalPreconditioner :: Update ()
  {
    cout << "Update Local Preconditioner" << flush;
    //  delete jacobi;

    if (block)
      {

	Table<int> * blocks = bfa->GetFESpace().CreateSmoothingBlocks();

	jacobi = dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix())
	  .CreateBlockJacobiPrecond(*blocks);

      }
    else
      {
	jacobi = dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix())
	  .CreateJacobiPrecond();
      }

    if (test) Test();
    //    cout << " done" << endl;
  }

  const ngla::BaseMatrix & LocalPreconditioner :: GetMatrix() const
  {
    return *jacobi;
  }


#ifdef OLD
  VEFC_Preconditioner :: VEFC_Preconditioner (PDE * pde, Flags & flags)
    : Preconditioner (flags)
  {
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    jacobi = NULL;
  }

  VEFC_Preconditioner :: ~VEFC_Preconditioner()
  {
    delete jacobi;
  }

  void VEFC_Preconditioner :: Update ()
  {
    cout << "Update VEFC_ Preconditioner" << flush;
    //  delete jacobi;

    
    jacobi = new VEFC_Matrix (bfa->GetMatrix(),
			      bfa->GetFESpace());
    /*
    Table<int> * blocks = bfa->GetFESpace().CreateSmoothingBlocks();

    jacobi = dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix())
      .CreateBlockJacobiPrecond(*blocks);
    */

    if (test) Test();
  }
#endif

  const ngla::BaseMatrix & VEFC_Preconditioner :: GetMatrix() const
  {
    return *jacobi;
  }


  TwoLevelPreconditioner :: TwoLevelPreconditioner (PDE * apde, Flags & flags)
  {
    pde = apde;
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", NULL));
    cpre = pde->GetPreconditioner (flags.GetStringFlag ("coarsepreconditioner", NULL));
    smoothingsteps = int (flags.GetNumFlag ("smoothingsteps", 1));
    premat = NULL;
  }

  TwoLevelPreconditioner :: ~TwoLevelPreconditioner()
  {
    delete premat;
  }

  void TwoLevelPreconditioner :: Update ()
  {
    /*
      delete premat;

      Smoother * smoother = new GSSmoother (pde->GetMeshAccess(), *bfa);
      //  Smoother * smoother = new EBESmoother (pde->GetMeshAccess(), *bfa);

      premat = new TwoLevelMatrix (&bfa->GetMatrix(),
      &cpre->GetMatrix(),
      smoother,
      pde->GetMeshAccess().GetNLevels());

      //			       bfa->GetMatrix().CreateJacobiPrecond());
      premat -> SetSmoothingSteps (smoothingsteps);
      premat -> Update();
    */

    /*
      cout << "2-Level Preconditioner" << endl;
      EigenSystem eigen (bfa->GetMatrix(), *premat);
      eigen.Calc();
      eigen.PrintEigenValues (cout);
    */
  }



  ComplexPreconditioner :: ComplexPreconditioner (PDE * apde, Flags & flags)
  {
    dim = int (flags.GetNumFlag ("dim", 1));
    cm = 0;
    creal = apde->GetPreconditioner (flags.GetStringFlag ("realpreconditioner",""));
  }

  ComplexPreconditioner :: ~ComplexPreconditioner()
  { 
    ;
  }

  void ComplexPreconditioner :: Update ()
  {
    delete cm;
    switch (dim)
      {
      case 1:
	cm = new Real2ComplexMatrix<double,Complex> (&creal->GetMatrix());
	break;
      case 2:
	cm = new Real2ComplexMatrix<Vec<2,double>,Vec<2,Complex> > (&creal->GetMatrix());
	break;
      case 3:
	cm = new Real2ComplexMatrix<Vec<3,double>,Vec<3,Complex> > (&creal->GetMatrix());
	break;
      case 4:
	cm = new Real2ComplexMatrix<Vec<4,double>,Vec<4,Complex> > (&creal->GetMatrix());
	break;
      default:
	cout << "Error: dimension " << dim << " for complex preconditioner not supported!" << endl;
      }
  }



  ChebychevPreconditioner :: ChebychevPreconditioner (PDE * apde, Flags & flags)
  {
    steps = int (flags.GetNumFlag ("steps",10.));
    cm = 0;
    csimple = apde->GetPreconditioner (flags.GetStringFlag ("csimple",""));
    bfa = apde->GetBilinearForm (flags.GetStringFlag ("bilinearform",""));

  }

  ChebychevPreconditioner :: ~ChebychevPreconditioner()
  { 
    ;
  }

  void ChebychevPreconditioner :: Update ()
  {
    delete cm;

    cout << "Compute eigenvalues csimple" << endl;
    const BaseMatrix & amat = bfa->GetMatrix();
    const BaseMatrix & pre = csimple->GetMatrix();

    EigenSystem eigen (amat, pre);
    eigen.SetPrecision(1e-30);
    eigen.SetMaxSteps(1000); 
    eigen.Calc();

    double lmin = eigen.EigenValue(1); 
    double lmax = eigen.MaxEigenValue(); 
    (*testout) << " Min Eigenvalue cs: "  << eigen.EigenValue(1) << endl; 
    (*testout) << " Max Eigenvalue cs : " << eigen.MaxEigenValue() << endl; 
    (cout) << " Min Eigenvalue cs: "  << eigen.EigenValue(1) << endl; 
    (cout)<< " Max Eigenvalue cs: " << eigen.MaxEigenValue() << endl; 
    (*testout) << " Condition csimple  " << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl; 
    (cout) << " Condition csimple" << eigen.MaxEigenValue()/eigen.EigenValue(1) << endl; 
  
    cm = new ChebyshevIteration(amat, pre, steps); 
    cm->SetBounds(1-lmax,1-lmin); 
    // Test(); 

  }





  ////////////////////////////////////////////////////////////////////////////////
  // added 08/19/2003, FB

  NonsymmetricPreconditioner :: NonsymmetricPreconditioner (PDE * apde, Flags & flags)
  {
    dim = int (flags.GetNumFlag ("dim", 1));
    cm = 0;
    cbase = apde->GetPreconditioner (flags.GetStringFlag ("basepreconditioner",""));
  }

  NonsymmetricPreconditioner :: ~NonsymmetricPreconditioner()
  { 
    ;
  }

  void NonsymmetricPreconditioner :: Update ()
  {
    delete cm;
    switch (dim)
      {
      case 2:
	cm = new Small2BigNonSymMatrix<double, Vec<2,double> > (&cbase->GetMatrix());
	break;
      case 4:
	cm = new Small2BigNonSymMatrix<Vec<2,double>, Vec<4,double> > (&cbase->GetMatrix());
	break;
      case 6:
	cm = new Small2BigNonSymMatrix<Vec<3,double>, Vec<6,double> > (&cbase->GetMatrix());
	break;
       case 8:
 	cm = new Small2BigNonSymMatrix<Vec<4,double>, Vec<8,double> > (&cbase->GetMatrix());
 	break;
//       case 2:
// 	cm = new Sym2NonSymMatrix<Vec<2,double> > (&cbase->GetMatrix());
// 	break;
//       case 4:
// 	cm = new Sym2NonSymMatrix<Vec<4,double> > (&cbase->GetMatrix());
// 	break;
//       case 6:
// 	cm = new Sym2NonSymMatrix<Vec<6,double> > (&cbase->GetMatrix());
// 	break;
//        case 8:
//  	cm = new Sym2NonSymMatrix<Vec<8,double> > (&cbase->GetMatrix());
//  	break;
      default:
	cout << "Error: dimension " << dim;
	cout << " for nonsymmetric preconditioner not supported!" << endl;
      }
  }





  ////////////////////////////////////////////////////////////////////////////////


  CommutingAMGPreconditioner :: CommutingAMGPreconditioner (PDE * apde, Flags & flags)
    : Preconditioner (flags), pde(apde)
  {
    bfa = pde->GetBilinearForm (flags.GetStringFlag ("bilinearform", ""));
    while (&bfa->GetLowOrderBilinearForm())
      bfa = &bfa->GetLowOrderBilinearForm();

    coefse = pde->GetCoefficientFunction (flags.GetStringFlag ("coefse", ""),1);    
    coefe = pde->GetCoefficientFunction (flags.GetStringFlag ("coefe", ""),1);    
    coeff = pde->GetCoefficientFunction (flags.GetStringFlag ("coeff", ""),1);    

    hcurl = dynamic_cast<const NedelecFESpace*> (&bfa->GetFESpace()) != 0;
    levels = int (flags.GetNumFlag ("levels", -1));
    coarsegrid = flags.GetDefineFlag ("coarsegrid");

    amg = 0;
  }

  CommutingAMGPreconditioner :: ~CommutingAMGPreconditioner ()
  {
    delete amg;
  }

  void CommutingAMGPreconditioner :: Update ()
  {
    cout << "Update amg" << endl;

    int i, j;
    const MeshAccess & ma = pde->GetMeshAccess();

    int nedge = ma.GetNEdges(); 
    int nface = ma.GetNFaces(); 
    int nel = ma.GetNE();

    if (coarsegrid && ma.GetNLevels() > 1)
      return;


    delete amg;

    ARRAY<INT<2> > e2v (nedge);
    for (i = 0; i < nedge; i++)
      {
	if (ma.GetClusterRepEdge (i) >= 0)
	  ma.GetEdgePNums (i, e2v[i][0], e2v[i][1]);
	else
	  e2v[i][0] = e2v[i][1] = -1;
      }


    ARRAY<INT<4> > f2v (nface);
    ARRAY<int> fpnums;
    for (i = 0; i < nface; i++)
      {
	if (ma.GetClusterRepFace (i) >= 0)
	  {
	    ma.GetFacePNums (i, fpnums);
	    
	    f2v[i][3] = -1;
	    for (j = 0; j < fpnums.Size(); j++)
	      f2v[i][j] = fpnums[j];
	  }
	else
	  {
	    cout << "coarse face" << endl;
	    for (j = 0; j < 4; j++)
	      f2v[i][j] = -1;
	  }
      }


    // compute edge and face weights:

    ARRAY<double> weighte(nedge), weightf(nface);
    ARRAY<double> hi(nedge), ai(nface);

    for (i = 0; i < nedge; i++)
      {
	Vec<3> p1, p2, v;
	ma.GetPoint (e2v[i][0], p1);
	ma.GetPoint (e2v[i][1], p2);
	v = p2 - p1;
	hi[i] = L2Norm (v);
      }

    if (hcurl)
    for (i = 0; i < nface; i++)
      {
	Vec<3> p1, p2, p3, v1, v2, vn;
	ma.GetPoint (f2v[i][0], p1);
	ma.GetPoint (f2v[i][1], p2);
	ma.GetPoint (f2v[i][2], p3);
	v1 = p2 - p1;
	v2 = p3 - p1;
	vn = Cross (v1, v2);
	ai[i] = L2Norm (vn);
      }

    ARRAY<int> ednums(12), edorient(12);
    ARRAY<int> fanums(12), faorient(12);
    LocalHeap lh (10000);
    ElementTransformation eltrans;
    IntegrationPoint ip(0, 0, 0, 0);


    weighte = 0;
    weightf = 0;
    for (i = 0; i < nel; i++)
      {
	lh.CleanUp();
	ma.GetElEdges (i, ednums, edorient);
	ma.GetElFaces (i, fanums, faorient);

	ma.GetElementTransformation (i, eltrans, lh);
	SpecificIntegrationPoint<3,3> sip(ip, eltrans, lh);

	double vol = ma.ElementVolume (i);
	double vale = Evaluate (*coefe, sip);

	for (j = 0; j < ednums.Size(); j++)
	  weighte[ednums[j]] += vale * vol / sqr (hi[ednums[j]]);

	if (hcurl)
	  {
	    double valf = Evaluate (*coeff, sip);
	    for (j = 0; j < fanums.Size(); j++)
	      weightf[fanums[j]] += valf * vol / sqr (ai[fanums[j]]);
	  }
      }

    int nsel = ma.GetNSE();
    if (coefse)
      for (i = 0; i < nsel; i++)
	{
	  lh.CleanUp();
	  ma.GetSElEdges (i, ednums, edorient);

	  ma.GetSurfaceElementTransformation (i, eltrans, lh);
	  SpecificIntegrationPoint<2,3> sip(ip, eltrans, lh);

	  double vol = ma.SurfaceElementVolume (i);
	  double vale = Evaluate (*coefse, sip);

	  for (j = 0; j < ednums.Size(); j++)
	    weighte[ednums[j]] += vale * vol / sqr (hi[ednums[j]]);
	}



    clock_t starttime, endtime, soltime;
    starttime = clock();

    CommutingAMG * amgmat;
    if (hcurl)
      amgmat = new AMG_HCurl (e2v, f2v, weighte, weightf, levels);
    else
      amgmat = new AMG_H1 (e2v, weighte, levels);

    endtime = clock();
    cout << "AMG coarsening time = " << double(endtime - starttime)/CLOCKS_PER_SEC << " sec" << endl;

    amgmat->ComputeMatrices (dynamic_cast<const BaseSparseMatrix&> (bfa->GetMatrix()));

    endtime = clock();
    cout << "AMG projection time = " << double(endtime - starttime)/CLOCKS_PER_SEC << " sec" << endl;

    cout << "Total NZE = " << amgmat->NZE() << endl;

    amg = amgmat;

    /*
    ChebyshevIteration * cheby =
      new ChebyshevIteration (bfa->GetMatrix(), *amgmat, 10);
    cheby -> SetBounds (0, 0.98);
    amg = cheby;
    */

    cout << "matrices done" << endl;

    if (timing) Timing();
    if (test) Test();
  }







  PreconditionerClasses::PreconditionerInfo::
  PreconditionerInfo (const string & aname,
	       Preconditioner* (*acreator)(const PDE & pde, const Flags & flags))
    : name(aname), creator(acreator)
  {
    ;
  }
  
  PreconditionerClasses :: PreconditionerClasses ()
  {
    ;
  }

  PreconditionerClasses :: ~PreconditionerClasses()
  {
    ;
  }
  
  void PreconditionerClasses :: 
  AddPreconditioner (const string & aname,
	      Preconditioner* (*acreator)(const PDE & pde, const Flags & flags))
  {
    prea.Append (new PreconditionerInfo(aname, acreator));
  }

  const PreconditionerClasses::PreconditionerInfo * 
  PreconditionerClasses::GetPreconditioner(const string & name)
  {
    for (int i = 0; i < prea.Size(); i++)
      {
	if (name == prea[i]->name)
	  return prea[i];
      }
    return 0;
  }

  void PreconditionerClasses :: Print (ostream & ost) const
  {
    ost << endl << "Preconditioners:" << endl;
    ost <<         "---------" << endl;
    ost << setw(20) << "Name" << endl;
    for (int i = 0; i < prea.Size(); i++)
      ost << setw(20) << prea[i]->name << endl;
  }

 
  PreconditionerClasses & GetPreconditionerClasses ()
  {
    static PreconditionerClasses fecl;
    return fecl;
  }




}