/*
 $Id: sysmatml.cc,v 1.3 1996/10/11 09:10:11 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "sysmatml.h"

#include "utils.h"
#include "numerics.h"
#include "dirichlet.h"
#include "family.h"

//-------------------------------------------------------------------------

StaticAllocator<NodeNeighbour>     NodeNeighbour::     alloc(ElementsInBlock);
StaticAllocator<AsymNodeNeighbour> AsymNodeNeighbour:: alloc(ElementsInBlock);

//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MLSparseMatrix:: MLSparseMatrix(symmetryType symmetry0, int spaceDim0, int dim) 

	: MLMatrix(symmetry0), MA28Matrix(spaceDim0),
	  dimension(dim), decomposed(False), ILUDecomposed(False),
	  D(dim), L(dim), symU(0), mark(dim), dirichlet(dim), smooth(dim)
{ 
    int i;

    FORALL(mark,i)   mark[i]   = False;
    FORALL(smooth,i) smooth[i] = True;
    FORALL(dirichlet,i) dirichlet[i] = False;

    FORALL(L,i) L[i] = 0;
    reset();
}
//-------------------------------------------------------------------------

MLSparseMatrix:: ~MLSparseMatrix() 
{
    removeLU();
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: extend(int noOfNodes)
{
    int i;

    if (noOfNodes != dimension)		// matrix will be really extended
    {
	MA28RemoveLUFactorization();
	decomposed    = False;
	ILUDecomposed = False;
    }

    dimension = noOfNodes;

    removeLU();

    D.resize(noOfNodes);	
    L.resize(noOfNodes);	
    FORALL(L,i) L[i] = 0;

    dirichlet.resize(noOfNodes);
    mark.resize(noOfNodes);
    smooth.resize(noOfNodes);
    FORALL(smooth,i) smooth[i] = True;

    FORALL(mark,i)   mark[i]   = False;
    FORALL(smooth,i) smooth[i] = True;
    FORALL(dirichlet,i) dirichlet[i] = False;

    reset();
}
//-------------------------------------------------------------------------


Num& MLSparseMatrix:: operator() (int row, int col)
{
    NodeNeighbour *np;
    int r,c;

    if (row == col) return D[row];
    
    if (row > col) { r = row; c = col; }
    else	   { r = col; c = row; }

    for (np=L[r]; np; np=np->next)
    {
	if (np->node == c) 
	{
	    if (row > col) return np->Aij();	
	    else	   return np->Aji();
	}
    }    
    
    // -- 	neighbour entry not found: insert a new entry
    
    NodeNeighbour* newNeigh;
    if (symmetry == sym) newNeigh = new NodeNeighbour;
    else	     	 newNeigh = new AsymNodeNeighbour;
    
    newNeigh->node = c;
    newNeigh->next = L[r];
    L[r] = newNeigh;
    
    if (row > col) return newNeigh->Aij();	
    else	   return newNeigh->Aji();
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: reset()
{
    int i;
    FORALL(D,i) resetRow(i);
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: resetRow(int node)
{
    NodeNeighbour* np;

    D[node] = 0.0;
    for (np=L[node]; np; np=np->next)  np->Aij() = np->Aji() = 0.0;
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: removeRow(int node)
{
    NodeNeighbour* np, *toRem;

    if (L[node] == 0) return;

    np = L[node];  
    L[node] = 0;

    while (np)
    {
	toRem = np;
	np = np->next;
	delete toRem;
    }
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: removeLU()
{
    int i;
    removeSymmetricExpansion();
    FORALL(L,i) removeRow(i);
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: markNodes()   { int i; FORALL(mark,i) mark[i] = 1; }
void MLSparseMatrix:: unMarkNodes() { int i; FORALL(mark,i) mark[i] = 0; }

//-------------------------------------------------------------------------

/*
// --	reset coupling terms between all marked (==changed) nodes

void MLSparseMatrix:: resetMarkedEntries()	
{
    NodeNeighbour* np;
    int i;

    FORALL(D,i) 
	if (mark[i])
	{
	    D[i] = 0.0;
	    for (np=L[i]; np; np=np->next) np->Aij() = np->Aji() = 0.0;
	}
	else
	{
	    for (np=L[i]; np; np=np->next)  
	      if (mark[np->node])  np->Aij() = np->Aji() = 0.0;
	}
}
*/
//-------------------------------------------------------------------------


void MLSparseMatrix:: print() const
{
    NodeNeighbour* np;
    int i;

    if (decomposed || ILUDecomposed)
    {
	cout << "\n* Matrix decomposed: print() not available\n";
	return;
    }

    cout << "\n--- MLSparse Matrix diagonal: \n\n";  
    FORALL(D,i) cout << D[i] << " ";
    
    cout << "\nLower Triangle:";
    FORALL(D,i) 
    {
	cout << "\nrow " << i << " : ";
	for (np=L[i]; np; np=np->next) 
	  cout << " " << np->Aij() << "[" << np->node << "]";
    }

    if (symmetry == asym)
    {
	cout << "\nUpper Triangle:";
	FORALL(D,i) 
	{
	    cout << "\ncol " << i << " : ";
	    for (np=L[i]; np; np=np->next) 
	      cout << " " << np->Aji() << "[" << np->node << "]";
	}
    }

    if (symU)
    {
	cout << "\nSymmetric Expansion: \n";
	FORALL(*symU,i) 
	{
	    cout << "\nrow " << i << " : ";
	    for (np=(*symU)[i]; np; np=np->next) 
	        cout << " " << np->Aij() << "[" << np->node << "]";
	}
    }
    cout << "\n";
}
//-------------------------------------------------------------------------


int MLSparseMatrix:: memSpace(Bool print) const
{
    int space = 0;

    space += L.memSpace();
    space += D.memSpace();
    space += mark.memSpace();
    space += smooth.memSpace();
    space += dirichlet.memSpace();
    if (symU) space += symU->memSpace();

    space += MA28MemSpace();

    if (print)
      cout << "\n    Memory allocated for MLSparseMatrix without LU: " 
           << Form("%1.6f", float(space)/1e6) << " MB\n";

    return space;
}
//-------------------------------------------------------------------------


int MLSparseMatrix:: memSpaceLU(Bool print) const
{
    int i, space = 0;

    FORALL(L,i) if (L[i]) { space += L[i]->alloc.MemSpace(); break; }

    if (print)
      cout << "\n    Memory allocated for LU in MLSparseMatrices: " 
           << Form("%1.6f", float(space)/1e6) << " MB\n";

    return space;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: smoothNode(int no, Vector<Num>& e, const Vector<Num>& r)
{
    NodeNeighbour* np;
    Num sum = 0.0;

    if (dirichlet[no]) return;

    for (np=L[no]; 	 np; np=np->next)  sum += np->Aij() * e[np->node];
    for (np=(*symU)[no]; np; np=np->next)  sum += np->Aij() * e[np->node];

    e[no] = (r[no] - sum) / D[no];
}
//-------------------------------------------------------------------------


Num MLSparseMatrix:: MultRow(int no, const Vector<Num>& e) 
{
    NodeNeighbour* np;

    Num sum = D[no] * e[no];

    for (np=L[no]; 	 np; np=np->next)  sum += np->Aij() * e[np->node];
    for (np=(*symU)[no]; np; np=np->next)  sum += np->Aij() * e[np->node];

    return sum;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: DiagDiv(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    int i;
    FORALL(D,i) lhs[i] = rhs[i] / D[i];
    if (!equal(omega,1.0)) FORALL(D,i) lhs[i] *= omega;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: DiagMult(Vector<Num>& lhs, Vector<Num>& rhs,Real omega)
{
    int i;
    FORALL(D,i) lhs[i] = rhs[i] * D[i];
    if (!equal(omega,1.0)) FORALL(D,i) lhs[i] *= omega;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: Mult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    NodeNeighbour* np;
    int i;

    if (lhs.v == rhs.v) callError("Mult");

    FORALL(D,i) lhs[i] = D[i] * rhs[i];
    FORALL(D,i) 
    {
	for (np=L[i]; np; np=np->next)
	{
	    lhs[i] 	  += np->Aij() * rhs[np->node];
	    lhs[np->node] += np->Aji() * rhs[i];	
	}
    }
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: ATMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    NodeNeighbour* np;
    int i;

    if (lhs.v == rhs.v) callError("ATMult");

    FORALL(D,i) lhs[i] = D[i] * rhs[i];

    FORALL(D,i) 
    {
	for (np=L[i]; np; np=np->next)
	{
	    lhs[i] 	  += conj(np->Aji()) * rhs[np->node];
	    lhs[np->node] += conj(np->Aij()) * rhs[i];	
	}
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: LMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    NodeNeighbour* np;
    int i;

    if (lhs.v == rhs.v) callError("LMult");
    FORALL(D,i) 
    {
	lhs[i] = 0.0;
	for (np=L[i]; np; np=np->next) lhs[i] += np->Aij() * rhs[np->node];
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: UMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    NodeNeighbour* np;
    int i;

    if (lhs.v == rhs.v) callError("UMult");

    FORALL(D,i) lhs[i] = 0.0;
    FORALL(D,i) 
    {
	for (np=L[i]; np; np=np->next) lhs[np->node] += np->Aji() * rhs[i];
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

	//  F = D + omega*L  of Matrix A;  lhs = F*rhs;  


void MLSparseMatrix:: F(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    NodeNeighbour* np;
    int i;
    Num sum;
    if (lhs.v == rhs.v) callError("F"); 

    FORALL(D,i) lhs[i] = D[i] * rhs[i];

    if (equal(omega,1.0))
    {
	FORALL_DOWNWARD(D,i) 
	{
	    sum = 0.0;
	    for (np=L[i]; np; np=np->next)  sum += np->Aij() * rhs[np->node];
	    lhs[i] += sum;
	}
    }
    else
    {
	FORALL_DOWNWARD(D,i) 
	{
	    sum = 0.0;
	    for (np=L[i]; np; np=np->next)  sum += np->Aij() * rhs[np->node];
	    lhs[i] += omega*sum;
	}
    }
}
//-------------------------------------------------------------------------

	//  FT = D + omega*U  of Matrix A;  lhs = FT*rhs 


void MLSparseMatrix:: FT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    NodeNeighbour* np;
    int i;
    if (lhs.v == rhs.v) callError("FT");
    
    FORALL(D,i) lhs[i] = 0.0;
    FORALL(D,i)
    {
        for (np=L[i]; np; np=np->next) 	lhs[np->node] += np->Aji() * rhs[i];
    }

    if (equal(omega,1.0)) FORALL(D,i) lhs[i] = lhs[i] + D[i] * rhs[i]; 
    else  		  FORALL(D,i) lhs[i] = omega*lhs[i] + D[i] * rhs[i]; 
}
//-------------------------------------------------------------------------

	//  Fm1:  F = D + omega*L  of Matrix A;  lhs = F**(-1) * rhs;


void MLSparseMatrix:: Fm1(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    NodeNeighbour* np;
    int i;
    Num sum;

    if (lhs.v != rhs.v)  FORALL(D,i) lhs[i] = rhs[i];

    if (equal(omega,1.0))
    {
	FORALL(D,i) 
	{
	    sum = 0.0;
	    for (np=L[i]; np; np=np->next)  sum += np->Aij() * lhs[np->node];
	    lhs[i] = (lhs[i] - sum) / D[i];
	}
    }
    else
    {
	FORALL(D,i) 
	{
	    sum = 0.0;
	    for (np=L[i]; np; np=np->next)  sum += np->Aij() * lhs[np->node];
	    lhs[i] = (lhs[i] - omega*sum) / D[i];
	}
    }
}
//-------------------------------------------------------------------------

	//  FmT:  FT = D + omega*U  of Matrix A;  lhs = F**(-t) * rhs;

void MLSparseMatrix:: FmT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    NodeNeighbour* np;
    int i;
    Num sum;

    if (lhs.v != rhs.v)  FORALL(D,i) lhs[i] = rhs[i];

    if (equal(omega,1.0))
    {
	FORALL_DOWNWARD(D,i)
	{
	    lhs[i] /= D[i];
	    sum = lhs[i];
	    for (np=L[i]; np; np=np->next)  lhs[np->node] -= sum * np->Aji();
	}
    }
    else
    {
	FORALL_DOWNWARD(D,i)
	{
	    lhs[i] /= D[i];
	    sum = omega*lhs[i];
	    for (np=L[i]; np; np=np->next) lhs[np->node] -= sum * np->Aji();
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Real MLSparseMatrix:: omegaOpt(int mode)  const
{
    NodeNeighbour* np;
    int i;
    Vector<Real> z(Dim());
    Num  skew, sym;
    Real omega, sum;
    
    FORALL(D,i) z[i] = 0.0;
 
    z[1] = 1.0;
    for (i=2; i<=D.high(); ++i)
        for (np=L[i]; np; np=np->next)
        {
	    sym  = 0.5*(np->Aij() + np->Aji());
	    skew = 0.5*(np->Aij() - np->Aji());
	    sum = Abs(sym) + Abs(skew);
	    z[i]        += sum;
	    z[np->node] += sum;
        }

    omega = 1.4;
    FORALL(D,i) if (z[i] > 0.0)  omega = Min(omega, Abs(D[i])/z[i]); 

    if (mode==multiGrid || mode==multiLevel)  return omega;
    else return 1./(omega - 0.5*(omega-1));   	// reduce and invert value
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

/*
   void MLSparseMatrix:: restoreDirichletEntries()
   {
   NodeNeighbour* np, *nextDir;
   int i;
   
   FORALL(D,i)
   {
   if (L[i] == 0)  			// a 'dirichlet row'
   L[i] = nodes[i]->firstDirichlet; 	
   
   else 						// move entry by entry
   {
   nextDir = nodes[i]->firstDirichlet;
   while (nextDir)
   {
   np = nextDir;
   nextDir = nextDir->next;
   
   np->next = nodes[i]->first;
   nodes[i]->first = np;
   }
   }
   nodes[i]->firstDirichlet = 0;
   }
   }
   //-------------------------------------------------------------------------
   
   void MLSparseMatrix:: removeDirichletEntries()
   {
   NodeNeighbour* np, *prev=0, *toRem;
   int i;
   
   FORALL(D,i) 
   {
   if (nodes[i]->dirichlet)		// remove all entries of row
   {
   np = nodes[i]->first;  
   nodes[i]->first = 0;
   
   while (np)
   {
   toRem = np;
   np = np->next;
   delete toRem;
   }
   }
   else		    // remove entries coupled to a dirichlet neighbour
   {
   np = nodes[i]->first;
   while (np)
   {
   if (nodes[np->node]->dirichlet)
   {
   if (np == nodes[i]->first) nodes[i]->first = np->next;
   else  prev->next = np->next;
   
   toRem = np;
   np = np->next;
   delete toRem;
   }
   else
   {
   prev = np;
   np = np->next;
   }
   }
   }
   }
   }
*/
//-------------------------------------------------------------------------


void MLSparseMatrix:: setDirichletBCs(DirichletBCs& dirichletBCs, Vector<Num>& b,
				      Vector<Num>& bSave, Num& Fct0, Num& E0)
{
    int i;

    if (symmetry == sym)  setSymDirichletBCs(dirichletBCs, b, bSave, Fct0, E0);
    else		 setASymDirichletBCs(dirichletBCs, b, bSave, Fct0, E0);

    dirichlet.resize(Dim());
    FORALL(dirichlet,i)
    {
	if (dirichletBCs.isSet(i)) dirichlet[i] = True;
	else 			   dirichlet[i] = False;
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: setSymDirichletBCs(DirichletBCs& dirichletBCs, 
					 Vector<Num>& b, Vector<Num>& /*bS*/, 
					 Num& Fct0, Num& /*E0*/)
{
    NodeNeighbour* np, *next, *prev=0;
    int i, icol;
    Num value;
    
    Fct0 = 0.0;
    
    FORALL(D,i)
    {
 	if (dirichletBCs.isSet(i))			// row with Dir. BC
 	{ 
	    value = dirichletBCs.value(i);
	    
	    for (np=L[i]; np; np=np->next)
	    {
		icol = np->node;

		if (dirichletBCs.isSet(icol)) 
			 Fct0 += np->Aij() * value * dirichletBCs.value(icol);
		else  b[icol] -= np->Aij() * value;

	    }

	    Fct0 -= b[i] * value;		// - u0*b0
	    Fct0 += D[i] * value * value;  
	    b[i]  = D[i] * value;

	    removeRow(i);
    	}
    	else
	{
	    for (np=L[i]; np; np=next)
	    {
		next = np->next;
		icol = np->node;

		if (dirichletBCs.isSet(icol))  
		{
		    b[i] -= np->Aij() * dirichletBCs.value(icol);      // ADF*u0

		    // -- 	remove np:

		    if (np == L[i]) L[i] = np->next;
		    else  	    prev->next = np->next;
		    delete np;
		}
		else  prev = np;
	    }    
	}
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: setASymDirichletBCs(DirichletBCs& dirichletBCs, 
					  Vector<Num>& b, Vector<Num>& bSave, 
					  Num& Fct0, Num& E0)
{
    NodeNeighbour* np;
    int i, icol;
    Num value;
    
    FORALL(bSave,i) bSave[i] = 0.0;  
    Fct0 = 0.0;
    E0   = 0.0;
    
    FORALL(D,i)
    {
 	if (dirichletBCs.isSet(i))			// row with Dir. BC
 	{ 
	    value = dirichletBCs.value(i);

	    for (np=L[i]; np; np=np->next)
	    {
		icol = np->node;

		if (dirichletBCs.isSet(icol)) 
			 Fct0 += 0.5 * np->Aij() * value * 
			 				dirichletBCs.value(icol);
		else  bSave[icol] += 0.5 * np->Aij() * value;
		np->Aij() = 0.0;
	    }

	    Fct0 -= b[i] * value;		// - u0*b0
	    Fct0 += D[i] * value * value; 
 
	    E0 += b[i] * value;
	    E0 -= D[i] * value * value;

	    b[i] = D[i] * value;
    	}

	for (np=L[i]; np; np=np->next) 		// column up
	{
	    icol = np->node;

	    if (dirichletBCs.isSet(icol))  
	    {
		if (dirichletBCs.isSet(i))
		      Fct0 += 0.5 * np->Aji() * dirichletBCs.value(i)
			                          * dirichletBCs.value(icol);
		else  bSave[i] += 0.5 * np->Aji() * dirichletBCs.value(icol);

		np->Aji() = 0.0;
	    }
	}
    }
}
//-------------------------------------------------------------------------

void MLSparseMatrix:: resetDirichletEntries(Vector<Num>& r) const
{
    int i;
    FORALL(dirichlet,i) { if (dirichlet[i])  r[i] = 0.0; }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: GalerkinRestriction(MLMatrix& Ah0, Generation& generation,
					  const Vector<SBool>* critical)
{
    NodeNeighbour* np;
    Son* son, *son2;
    int  i, k, m, nFathers, nFathers2, father, father2;
    Real weight, weight2, fact;

    if (generation.castToMCGeneration()) 
    {
	cout << "\n*** MLSparseMatrix : generation is MCGeneration\n";
	cout.flush(); abort();
    }

    MLSparseMatrix* Ah = Ah0.castToMLSparseMatrix();
    removeLU();
    reset();

    FORALL(Ah->D,i) 
    {
	if (critical) if ((*critical)[i]) continue;

	son = generation.getSon(i);
	nFathers = son->NoOfFathers();


	// --		restrict diagonal term:

	if (nFathers == 1) 
	{
	    father = son->getFather();
	    if (mark[father]) D[father] += Ah->D[i];
	}
	else
	{
	    for (k=1; k<=nFathers; ++k)
	    for (m=1; m<=k; ++m)
	    {
		father  = son->getFather(k);
		father2 = son->getFather(m);
		
		if (dirichlet[father] || dirichlet[father2]) continue;

		if (mark[father] || mark[father2])
		{
		    weight  = son->getWeight(k);
		    weight2 = son->getWeight(m);
		    fact = weight * weight2;

		    (*this)(father,father2) += fact * Ah->D[i];

		    if (symmetry==asym && father!=father2) 
		      (*this)(father2,father) += fact * Ah->D[i];
		}
	    }
	}
    
    
	// --		off-diagonal terms:
	
	for (np=Ah->L[i]; np; np=np->next) 
	{
	    if (critical) if ((*critical)[np->node]) continue;

	    son2 = generation.getSon(np->node);
	    nFathers2 = son2->NoOfFathers();
	    
	    for (k=1; k<=nFathers;  ++k)
	    for (m=1; m<=nFathers2; ++m)
	    {	
		father  = son ->getFather(k);
		father2 = son2->getFather(m);

		if (dirichlet[father] || dirichlet[father2]) continue;

		if (mark[father] || mark[father2]) 
		{	       
		    weight  = son ->getWeight(k);
		    weight2 = son2->getWeight(m);

		    fact = weight*weight2;

		    if (symmetry == sym)
		    {
			if (father == father2) fact *= 2.0; 
			(*this)(father,father2) += fact * np->Aij();
		    }
		    else
		    {
			(*this)(father,father2) += fact * np->Aij();
			(*this)(father2,father) += fact * np->Aji();
		    }
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: symmetricExpansion()
{
    NodeNeighbour* np, *newN;
    int i;

    removeSymmetricExpansion();

    symU = new Vector<NodeNeighbour*>(Dim());
    FORALL(*symU,i) (*symU)[i] = 0;

    FORALL(D,i) 
    {
	for (np=L[i]; np; np=np->next)
	{
	    newN = new NodeNeighbour;
	    
	    newN->Aij() = np->Aji();
	    newN->node = i;			// exchange column and row
	    
	    newN->next = (*symU)[np->node];
	    (*symU)[np->node] = newN;
	}
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: removeSymmetricExpansion()
{
    int i;
    NodeNeighbour* np, *toRem;

    if (symU == 0) return;

    FORALL(*symU,i)
    {
	np = (*symU)[i];
	while(np)
	{
	    toRem = np;
	    np = np->next;
	    delete toRem;
	}
    }

    delete symU;
    symU = 0;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void LocalMLSparseMatrix:: smoothNode(int node, Vector<Num>& e, 
				      const Vector<Num>& r)
{
    NodeNeighbour* np;
    Num sum = 0.0;

    if (!smooth[node])   return;
    if (dirichlet[node]) return;

    for (np=L[node];       np; np=np->next)  sum += np->Aij() * e[np->node];
    for (np=(*symU)[node]; np; np=np->next)  sum += np->Aji() * e[np->node];

    e[node] = (r[node] - sum) / D[node];
}
//-------------------------------------------------------------------------


void LocalMLSparseMatrix:: DiagDiv(Vector<Num>& lhs, Vector<Num>& rhs, 
				   Real omega)
{
    int i;
    if (equal(omega,1.0)) 
      FORALL(D,i) 
      {
	  if (smooth[i]) lhs[i] = rhs[i] / D[i];
	  else 		 lhs[i] = 0.0;
      }
    else
      FORALL(D,i) 
      {
	  if (smooth[i]) lhs[i] = omega * rhs[i] / D[i];
	  else 		 lhs[i] = 0.0;
      }
}
//-------------------------------------------------------------------------


void LocalMLSparseMatrix:: DiagMult(Vector<Num>& lhs, Vector<Num>& rhs,
				    Real omega)
{
    int i;
    if (equal(omega,1.0)) 
      FORALL(D,i) 
      {
	  if (smooth[i]) lhs[i] = rhs[i] * D[i];
	  //else lhs[i] = 0.0;
      }
    else	          
      FORALL(D,i) 
      {
	  if (smooth[i]) lhs[i] = omega * rhs[i] * D[i];
	  //else lhs[i] = 0.0;
      }
}
//-------------------------------------------------------------------------

	//  F = D + omega*L  of Matrix A;  lhs = F*rhs;  

void LocalMLSparseMatrix:: F(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)
{
    notImplemented("LocalMLSparseMatrix:: F");
}
//-------------------------------------------------------------------------

	//  FT = D + omega*U  of Matrix A;  lhs = FT*rhs 

void LocalMLSparseMatrix:: FT(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)
{
    notImplemented("LocalMLSparseMatrix:: F");
}
//-------------------------------------------------------------------------


	//  Fm1:  F = D + omega*L  of Matrix A;  lhs = F**(-1) * rhs;


void LocalMLSparseMatrix:: Fm1(Vector<Num>& lhs, Vector<Num>& rhs, Real omega) 
{
    NodeNeighbour* np;
    int i;
    Num sum;

    FORALL(D,i) { if (smooth[i]) lhs[i] = rhs[i]; }

    if (equal(omega,1.0))
    {
	FORALL(D,i) 
	    if (smooth[i])
	    {
		sum = 0.0;
		for (np=L[i]; np; np=np->next)  sum += np->Aij() * lhs[np->node];
		lhs[i] = (lhs[i] - sum) / D[i];
	    }
    }
    else
    {
	FORALL(D,i) 
	    if (smooth[i])
	    {
		sum = 0.0;
		for (np=L[i]; np; np=np->next)  sum += np->Aij() * lhs[np->node];
		lhs[i] = (lhs[i] - omega*sum) / D[i];
	    }
    }
}
//-------------------------------------------------------------------------

	//  FmT:  FT = D + omega*U  of Matrix A;  lhs = F**(-T) * rhs;

void LocalMLSparseMatrix:: FmT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    NodeNeighbour* np;
    int i;
    Num sum;

    FORALL(D,i) { if (smooth[i]) lhs[i] = rhs[i]; }
      
    if (equal(omega,1.0))
    {
	FORALL_DOWNWARD(D,i)
	{
	    if (smooth[i]) lhs[i] /= D[i];
	    sum = lhs[i];
	    for (np=L[i]; np; np=np->next)  
	        if (smooth[np->node]) lhs[np->node] -= sum * np->Aji();
	}
    }
    else
    {
	FORALL_DOWNWARD(D,i)
	{
	    if (smooth[i]) lhs[i] /= D[i];
	    sum = omega*lhs[i];
	    for (np=L[i]; np; np=np->next)  
	        if (smooth[np->node]) lhs[np->node] -= sum * np->Aji();
	}
    }
}
//-------------------------------------------------------------------------


void LocalMLSparseMatrix:: removeNonSmoothedEntries()
{
    NodeNeighbour* np, *prev=0, *toRem;
    int i;

    FORALL(D,i) 
    {
	if (!smooth[i])	
	{
	    D[i] = 0.0;			// may never be used
	    np = L[i];

	    while (np)
	    {
		if (!smooth[np->node])
		{
		    if (np == L[i]) L[i] = np->next;
		    else  	    prev->next = np->next;

		    toRem = np;
		    np = np->next;
		    delete toRem;
		}
		else
		{
		    prev = np;
		    np = np->next;
		}
	    }
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void MLSparseMatrix:: countEntries(int* dim, int* nEntries)
{
    int i;
    NodeNeighbour* np;

    int n = dimension;
    FORALL(D,i) { for (np=L[i]; np; np=np->next)  n += 2; }

    *nEntries = n;
    *dim = dimension;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: Decompose(Bool removeOriginal)
{
    if (decomposed) return;

    MA28Decompose();

    if (removeOriginal) { removeLU(); D.resize(1); }
    decomposed = True;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: ILUDecompose(Real dropTol)
{
    if (ILUDecomposed || decomposed) return;

    MA28Decompose(dropTol);

    ILUDecomposed = True;
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: FBSubst(Vector<Num>& lhs, Vector<Num>& rhs)
{
    if (!decomposed)
    {
	cout << "\n*** MLSparseMatrix::FBSubst(...) : Matrix not decomposed\n";
	cout.flush(); abort();
    }

    int i;
    FORALL(lhs,i) lhs[i] = rhs[i];
    Num* b = lhs.v + lhs.low();

    MA28FBSubst(b);
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: ILUFBSubst(Vector<Num>& lhs, Vector<Num>& rhs)
{
    int i;

    if (!ILUDecomposed) FBSubst(lhs, rhs);
    else
    {
	FORALL(lhs,i) lhs[i] = rhs[i];
	Num* b = lhs.v + lhs.low();
	
	MA28FBSubst(b);
    }
}
//-------------------------------------------------------------------------


void MLSparseMatrix:: fillMA28Vectors()
{
    int i, n;
    NodeNeighbour* np;

    n = -1;			// start with 0 in ma28-arrays

    FORALL(D,i)  		// array A must start with diagonal 
    {
	++n;
	A  [n] = D[i];
	IRN[n] = i;
	ICN[n] = i;
    }

    FORALL(D,i) 
    {
        for (np=L[i]; np; np=np->next)
	{
	    ++n;
	    A  [n] = np->Aij();
	    IRN[n] = i;
	    ICN[n] = np->node;
	    ++n;
	    A  [n] = np->Aji();
	    IRN[n] = np->node;
	    ICN[n] = i;
	}
    }    
}