/*
 $Id: sysmatsp.cc,v 1.4 1997/07/11 09:46:14 bzferdma Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "sysmatsp.h"

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

#include "cmdpars.h"
extern CmdPars Cmd;


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


SparseMatrix:: SparseMatrix(ConnectionPattern& cPattern, symmetryType symmetry0,
			    int spaceDim0)

      : SystemMatrix(symmetry0),  MA28Matrix(spaceDim0),
	dimension(0), decomposed(False), ILUDecomposed(False),
	D(1), L(0), U(0), col(1), end(1),
	dirichlet(1), symU(0), symCol(1), symEnd(1)
{
    
    allocate(cPattern);
    dimension = D.high();
}
//-------------------------------------------------------------------------


// --	construct A for eigenvalue problems (A = A0 + lambda*B0)
//      or simply make a copy (B=0)


SparseMatrix:: SparseMatrix(int spaceDim0, SystemMatrix* A0, SystemMatrix* B0, 
			    Num lambda)

	  : SystemMatrix(A0->symmetry), MA28Matrix(spaceDim0),
	    D(1), L(0), U(0), col(1), end(1), dirichlet(1), 
	    symU(0), symCol(1), symEnd(1)
{
    int i, k;
    SparseMatrix* Asp, *B;

    Asp = A0->castToSparseMatrix();

    if (B0) B = B0->castToSparseMatrix();
    else    B = 0;

    const int nlow  = (Asp->D).l;
    const int nhigh = (Asp->D).h;

    D.resize(nlow, nhigh);

    if (B) FORALL(D,i) D[i] = (Asp->D)[i] - lambda * (B->D)[i];
    else   FORALL(D,i) D[i] = (Asp->D)[i];

    int noOfEntries = (Asp->L)->size();

    if (noOfEntries == 0)  			// no off-diagonal terms
    {
	L = U = 0;
	return;
    }

    end.resize(nlow-1,nhigh);
    col.resize(noOfEntries);
    L = new Vector<Num>(noOfEntries);
    if (symmetry == asym) U = new Vector<Num>(noOfEntries);
    else 		  U = L;
    
    FORALL(Asp->end,i)  end[i] = (Asp->end)[i];

    if (B)
    {
	FORALL(D,i) 
	    for (k=end[i-1]+1; k<=end[i]; ++k)
	    {
		col[k] = Asp->col[k];
		(*L)[k] = (*Asp->L)[k] - lambda * (*B->L)[k];
		if (symmetry == asym) (*U)[k] = (*Asp->U)[k] - lambda * (*B->U)[k];
	    }
    }
    else
    {
	FORALL(D,i) 
	    for (k=end[i-1]+1; k<=end[i]; ++k)
	    {
		col[k] = Asp->col[k];
		(*L)[k] = (*Asp->L)[k];
		(*U)[k] = (*Asp->U)[k];
	    }
    }

    dimension = D.high() - D.low() + 1;
}
//-------------------------------------------------------------------------


SparseMatrix:: ~SparseMatrix()
{ 
    removeSymmetricExpansion();

    if (L == U) delete L; 
    else 
    { 
	delete L; 
	delete U; 
    }
}
//-------------------------------------------------------------------------


void SparseMatrix:: allocate(ConnectionPattern& cPattern)
{
    int i, n;

    const int nlow  = cPattern.l();
    const int nhigh = cPattern.h();

    D.resize(nlow, nhigh);

    int noOfEntries = 0;
    for (i=nlow; i<=nhigh; ++i) noOfEntries += cPattern.noOfElements(i);

    if (noOfEntries == 0)		// no off-diagonal terms
    { 
	L = U = 0;
	return;
    }

    end.resize(nlow-1,nhigh);
    col.resize(noOfEntries);

    L = new Vector<Num>(noOfEntries);
    if (symmetry == asym) U = new Vector<Num>(noOfEntries);
    else 		  U = L;

    NeighbourNode* nnp;

    int index = 0;
    end[nlow-1] = index;		// care for end[i-1] in case i=nlow !

    for (i=nlow; i<=nhigh; ++i) 	// all rows
    {
	for (nnp=cPattern.first(i); nnp; nnp=nnp->next) 
	{
	    n = nnp->neighbour;
	    if (n >= i)
	    {
		cout << "\n*** Error in sparse node pattern: row = " << i
		     << "  column = " << n << "\n";
		cout.flush(); abort();
	    }
	    col[++index] = n;
	}
	end[i] = index;
    }
}
//-------------------------------------------------------------------------

void SparseMatrix:: removeLU()
{
    if (U != L) delete U; 
    delete L;
    L = U = 0;
    end.resize(1);
    col.resize(1);
}
//-------------------------------------------------------------------------

void SparseMatrix:: printMatLabFormat() const
{
    char s[10];
    printf ("\n---  print A ?  <No> ");  fgets(s,10,stdin);
    strToLower(s);
    if (!strchr(s,'y'))  return;

    FILE* f;
    f = fopen("sparsematrix","w");

    int i, j;
    FORALL(D,i) 
    {
	for (j=end[i-1]+1; j<=end[i]; ++j)  
	if (L) fprintf(f,"a(%1d,%1d)=%g\n", i, col[j], ::real((*L)[j]));
	fprintf(f,"a(%1d,%1d)=%g\n", i, i, ::real(D[i]));
    }
    fclose(f);
    printf("--- file 'sparsematrix' created \n");
}
//-------------------------------------------------------------------------

void SparseMatrix:: print() const
{
    int i, j;

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

    cout << "\n--- Sparse Matrix diagonal: \n\n";  D.print();
    
    cout << "\nLower Triangle:";

    if (L)
      FORALL(D,i) 
      {
	  cout << "\nrow " << i << " : ";
	  for (j=end[i-1]+1; j<=end[i]; ++j)  
	    cout  << " " << (*L)[j] << "[" << col[j] << "]";
      }

    if (U != L)
    {
	cout << "\n\nUpper Triangle:";
	FORALL(D,i) 
	{
	    cout << "\ncol " << i << " : ";
	    for (j=end[i-1]+1; j<=end[i]; ++j)  
	      cout  << " " << (*U)[j] << "[" << col[j] << "]";
	}
    }
    
    if (symU)
    {
	cout << "\n\nSymmetric Upper Triangle:";
	FORALL(D,i) 
	{
	    cout << "\nrow " << i << " : ";
	    for (j=symEnd[i-1]+1; j<=symEnd[i]; ++j)  
	      cout  << " " << (*symU)[j] << "[" << symCol[j] << "]";
	}
    }

    cout << "\n";
}
//-------------------------------------------------------------------------


int SparseMatrix:: memSpace(int print) const
{
    int space = D.memSpace() + col.memSpace() + end.memSpace();
    if (L) space += L->memSpace();
    if (symmetry == asym && U) space += U->memSpace();

    space += dirichlet.memSpace();
    space += MA28MemSpace();

    if (print)
      cout << "\n    Memory allocated for sparse matrix (dim=" << Dim()
           << "): " << Form("%1.6f", float(space)/1e6) << " MB\n";

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


void SparseMatrix:: 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 SparseMatrix:: setASymDirichletBCs(DirichletBCs& dirichletBCs,
					Vector<Num>& b, Vector<Num>& bSave, 
					Num& Fct0, Num& E0)
{
    int i, j, 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 (j=end[i-1]+1; j<=end[i]; ++j) 		// along row
	    {
		icol = col[j];

		if (dirichletBCs.isSet(icol)) 
			     Fct0 += 0.5 * (*L)[j] * value *
			    				dirichletBCs.value(icol);
		else  bSave[icol] += 0.5 * (*L)[j] * value;
		(*L)[j] = 0.0;
	    }
	    Fct0 -= b[i]*value;			// - u0*b0
	    Fct0 += D[i]*value*value;  		// u0*Add*u0 

	    E0 += b[i]*value;
	    E0 -= D[i]*value*value;

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

	for (j=end[i-1]+1; j<=end[i]; ++j) 		// column up
	{
	    icol = col[j];			// here icol is the row index

	    if (dirichletBCs.isSet(icol))
	    {
		if (dirichletBCs.isSet(i))
		          Fct0 += 0.5* (*U)[j] * dirichletBCs.value(i)
			                          * dirichletBCs.value(icol);
		else  bSave[i] += 0.5* (*U)[j] * dirichletBCs.value(icol);

		(*U)[j] = 0.0;
	    }
	}
    }
}
//-------------------------------------------------------------------------


void SparseMatrix:: setSymDirichletBCs(DirichletBCs& dirichletBCs,
				       Vector<Num>& b, Vector<Num>& /*bSave*/, 
				       Num& Fct0, Num& /*E0*/)
{
    int i, j, icol;
    Num value;

    Fct0 = 0.0;

    FORALL(D,i)
    {
 	if (dirichletBCs.isSet(i))			// row with Dir. BC
 	{ 
	    value = dirichletBCs.value(i);
	    
	    for (j=end[i-1]+1; j<=end[i]; ++j) 
	    {
		icol = col[j];

		if (dirichletBCs.isSet(icol)) 
			 Fct0 += (*L)[j] * value * dirichletBCs.value(icol);
		else  b[icol] -= (*L)[j] * value;

		(*L)[j] = 0.0;
	    }
	    Fct0 -= b[i]*value;		// - u0*b0
	    Fct0 += D[i]*value*value;  	// u0*Add*u0    (0.5*... + 0.5*...)
	    b[i]  = D[i]*value;
    	}
    	else
	{
	    for (j=end[i-1]+1; j<=end[i]; ++j)
	    {
		icol = col[j];
		if (dirichletBCs.isSet(icol))  
		{
		    b[i] -= (*L)[j] * dirichletBCs.value(icol);	   // ADF*u0
		    (*L)[j] = 0.0;
		}
	    }    
	}
    }
}
//-------------------------------------------------------------------------


Num& SparseMatrix:: operator() (int n1, int n2)
{
    int i;

    //if (n1==0 || n2==0) 	// !!!
    //	{ cout << "\n*** SparseMatrix::operator(): node number 0\n"; cout.flush(); abort(); }

    if (n1 == n2) return D[n1];
    if (n1 > n2)			// lower triangle
    {
	i = end[n1];
	while (col[i] != n2) --i;
	return (*L)[i];
    }
    else 
    {
	i = end[n2];                
	while (col[i] != n1) --i;   
	return (*U)[i];            
    }
}
//-------------------------------------------------------------------------

void SparseMatrix:: reset()
{
    int i, k;
    FORALL(D,i) D[i] = 0.0;
    FORALL(D,i) for (k=end[i-1]+1; k<=end[i]; ++k)  (*L)[k] = (*U)[k] = 0.0;
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void SparseMatrix:: Mult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    int i, k;
    if (lhs.v == rhs.v) callError("Mult");

    FORALL(D,i) lhs[i] = D[i] * rhs[i];
    
    FORALL(D,i) 
        for (k=end[i-1]+1; k<=end[i]; ++k)
	{
	    lhs[i]      += (*L)[k] * rhs[col[k]];
	    lhs[col[k]] += (*U)[k] * rhs[i];
	}
}
//-------------------------------------------------------------------------

void SparseMatrix:: ATMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    int i, k;
    if (lhs.v == rhs.v) callError("ATMult");

    FORALL(D,i) lhs[i] = conj(D[i]) * rhs[i];
    
    FORALL(D,i) 
        for (k=end[i-1]+1; k<=end[i]; ++k)
	{
	    lhs[i]      += conj((*U)[k]) * rhs[col[k]];
	    lhs[col[k]] += conj((*L)[k]) * rhs[i];
	}
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


void SparseMatrix:: smoothNode(int no, Vector<Num>& e, const Vector<Num>& r)
{
    int k;
    Num sum = 0.0;

    for (k=end[no-1]+1;    k<=end[no];    ++k)  sum += (*L)[k]    * e[col[k]];

    if (symU)
    for (k=symEnd[no-1]+1; k<=symEnd[no]; ++k)  sum += (*symU)[k] * e[symCol[k]];

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


Num SparseMatrix:: MultRow(int row, const Vector<Num>& e)
{
    int k;
    Num sum = D[row] * e[row];

    for (k=end[row-1]+1;    k<=end[row];    ++k) sum += (*L)[k]   * e[col[k]];

    if (symU)
    for (k=symEnd[row-1]+1; k<=symEnd[row]; ++k) sum += (*symU)[k]* e[symCol[k]];

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

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

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


			//  lhs = L * rhs: 

void SparseMatrix:: LMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    int i, k;

    if (lhs.v == rhs.v) callError("LMult");
    FORALL(D,i) 
    {
	lhs[i] = 0.0;
        for (k=end[i-1]+1; k<=end[i]; ++k) lhs[i] += (*L)[k] * rhs[col[k]];
    }
}
//-------------------------------------------------------------------------

			//  lhs = UMult * rhs

void SparseMatrix:: UMult(Vector<Num>& lhs, Vector<Num>& rhs)
{
    int i, k;

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

    FORALL(D,i) lhs[i] = 0.0;
    FORALL(D,i)
        for (k=end[i-1]+1; k<=end[i]; ++k)  lhs[col[k]] += (*U)[k] * rhs[i];
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------

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


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

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

    if (equal(omega,1.0))
    {
	FORALL(D,i) {
	    sum = 0.0;
	    for (k=end[i-1]+1; k<=end[i]; ++k)  sum += (*L)[k] *rhs[col[k]];
	    lhs[i] += sum;
	}
    }
    else
    {
	FORALL(D,i) {
	    sum = 0.0;
	    for (k=end[i-1]+1; k<=end[i]; ++k)  sum += (*L)[k] *rhs[col[k]];
	    lhs[i] += omega*sum;
	}
    }
}
//-------------------------------------------------------------------------

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


void SparseMatrix:: FT(Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    int i, k;
    if (lhs.v == rhs.v) callError("FT");

    FORALL(D,i) lhs[i] = 0.0;
    FORALL(D,i)
        for (k=end[i-1]+1; k<=end[i]; ++k)  lhs[col[k]] += (*U)[k] * 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 SparseMatrix:: Fm1 (Vector<Num>& lhs, Vector<Num>& rhs, Real omega)
{
    int i, k;
    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 (k=end[i-1]+1; k<=end[i]; ++k)  sum += (*L)[k] * lhs[col[k]];
	    lhs[i] = (lhs[i] - sum) / D[i];
	}
    }
    else
    {
	FORALL(D,i) {
	    sum = 0.0;
	    for (k=end[i-1]+1; k<=end[i]; ++k)  sum += (*L)[k] * lhs[col[k]];
	    lhs[i] = (lhs[i] - omega*sum) / D[i];
	}
    }
}
//-------------------------------------------------------------------------

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

void SparseMatrix:: FmT (Vector<Num>& lhs, Vector<Num>& rhs, Real omega) 
{
    int i, k;
    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 (k=end[i-1]+1; k<=end[i]; ++k)  lhs[col[k]] -= sum * (*U)[k];
	}
    }
    else
    {
	FORALL_DOWNWARD(D,i)
	{
	    lhs[i] /= D[i];
	    sum = omega*lhs[i];
	    for (k=end[i-1]+1; k<=end[i]; ++k)  lhs[col[k]] -= sum * (*U)[k];
	}
    }
}
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


Real SparseMatrix:: omegaOpt(int mode)  const
{
    int i, k;
    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.h; i++)
        for (k=end[i-1]+1; k<=end[i]; ++k)
        {
	    sym  = 0.5*((*L)[k] + (*U)[k]);
	    skew = 0.5*((*L)[k] - (*U)[k]);
	    sum = Abs(sym) + Abs(skew);
	    z[i]      += sum;
	    z[col[k]] += 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 SparseMatrix:: resetDirichletEntries(Vector<Num>& r) const
{
    int i;
    FORALL(dirichlet,i) { if (dirichlet[i]) r[i] = 0.0; }
}
//-------------------------------------------------------------------------


// create the upper triangle of the matrix (a formal symmetric expansion of U)


void SparseMatrix:: symmetricExpansion()
{
    int i, k, pos;

    delete symU;
    symU = new Vector<Num>(U->low(),U->high());
    symCol.resize(col.low(),col.high());
    symEnd.resize(end.low(),end.high());

    Vector<int> noOfEntries(Dim());

    FORALL(noOfEntries,i) noOfEntries[i] = 0;

    FORALL(D,i) 		// count entries in each row in upper triangle
    {
        for (k=end[i-1]+1; k<=end[i]; ++k)  ++noOfEntries[col[k]];
    }
  
    symEnd[symEnd.l] = 0;				// for row '-1'
    FORALL(noOfEntries,i) symEnd[i] = symEnd[i-1] + noOfEntries[i];


    FORALL(noOfEntries,i) noOfEntries[i] = 0;

    FORALL(D,i) 		// fill the upper triangle
    {
        for (k=end[i-1]+1; k<=end[i]; ++k)  
	{
	    pos = symEnd[col[k]-1] + 1 + noOfEntries[col[k]];
	    ++noOfEntries[col[k]];

	    symCol[pos]  = i;
	    (*symU)[pos] = (*U)[k];
	}
    }
}
//-------------------------------------------------------------------------

void SparseMatrix:: removeSymmetricExpansion()
{
    if (symU==0) return;

    delete symU;  symU = 0;
    symCol.resize(1);
    symEnd.resize(1);
}
//-------------------------------------------------------------------------

/*   //             Sparse Incomplete LU-Decomposition: A = L*U = L*D*LT

   void SparseMatrix:: ILUDecompose(Real weight)
   {
   int  i, k, m, n, r;
   Bool found;
   Num  e;
   
   symmetricExpansion();
   
   if (Cmd.isTrue("printILU"))  
   { cout << "\nMatrix before Factorization:\n";  print(); }
   
   for (r=1; r < Dim(); ++r)
   {
   for (i=r+1; i <= Dim(); ++i)	               // all subsequent rows i
   {
   for (k=end[i-1]+1; k<=end[i]; ++k) 	       // search A(i,r), the
   { 						         // 'row factor'
   if (col[k]==r && !equal((*L)[k],0.0))	    // A(i,r) != 0
   {
   e = (*L)[k]/D[r];       	      // e = l(i,r)=A(i,r)/A(r,r)
   (*L)[k] = e;
   
   // --  apply e to the rest of this row 
   // 					in the upper triangle:
   
   for (m=symEnd[r-1]+1; m<=symEnd[r]; ++m)
   {
   if (!equal((*symU)[m],0.0))   	// (*symU)[m]=A(r,j)
   {
   const j = symCol[m];
   
   // --  search corresponding element of column j:
   
   found = False;
   
   if (i==j)  			// is it the diagonal ?
   {
   D[i] -= e * (*symU)[m];
   found = True;
   }
   if (!found)			// search upper triangle
   {
   for (n=symEnd[i-1]+1; n<=symEnd[i]; ++n)  
   {		 
   if (symCol[n]==j)
   {
   (*symU)[n] -= e * (*symU)[m];
   found = True;
   }
   }
   }
   if (!found)			// ... and lower triangle
   {
   for (n=end[i-1]+1; n<=end[i]; ++n)  
   {		 
   if (col[n]==j)
   {
   (*L)[n] -= e * (*symU)[m];
   found = True;
   }
   }
   }
   if (!found)  D[i] -= weight * e * (*symU)[m];
   }
   } 
   break;
   }
   }
   }
   }
   
   if (Cmd.isTrue("printILU"))  
   { cout << "\nMatrix after Factorization:\n";  print(); }
   
   if (symmetry != sym) 
   {
   //  store symU to U !!!
   cout << "\n*** SparseMatrix:: ILUDecompose not implemented for"
   << " non-symmetric matrices\n";
   cout.flush(); abort();
   }    
   
   removeSymmetricExpansion();
   }
   //-------------------------------------------------------------------------
   
   //   A**(-1) = L**(-T) * D**(-1) * L**(-1)
   
   void SparseMatrix:: ILUFBSubst(Vector<Num>& lhs, Vector<Num>& rhs) const 
   {
   int i, k;
   Num sum;
   
   if (lhs.v != rhs.v)  FORALL(D,i) lhs[i] = rhs[i];
   
   FORALL(D,i) 
   {
   sum = 0.0;
   for (k=end[i-1]+1; k<=end[i]; ++k)  sum += (*L)[k] * lhs[col[k]];
   lhs[i] = lhs[i] - sum;
   }
   
   FORALL(D,i)  lhs[i] /= D[i];
   
   FORALL_DOWNWARD(D,i)
   {
   sum = lhs[i];
   for (k=end[i-1]+1; k<=end[i]; ++k)  lhs[col[k]] -= sum * (*U)[k];
   }
   }
   */
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------
   
   
void SparseMatrix:: countEntries(int* dim, int* nEntries)
{
    int i, k;

    int n = dimension;
    FORALL(D,i) { for (k=end[i-1]+1; k<=end[i]; ++k)  n += 2; }	

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


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

    MA28Decompose();

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


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

    MA28Decompose(dropTol);

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


void SparseMatrix:: FBSubst(Vector<Num>& lhs, Vector<Num>& rhs)
{
    if (!decomposed)
    {
	cout << "\n*** SparseMatrix::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 SparseMatrix:: 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 SparseMatrix:: fillMA28Vectors()
{
    int i,k,n;

    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 (k=end[i-1]+1; k<=end[i]; ++k)
	{
	    ++n;
	    A  [n] = (*L)[k];
	    IRN[n] = i;
	    ICN[n] = col[k];
	    ++n;
	    A  [n] = (*U)[k];
	    IRN[n] = col[k];
	    ICN[n] = i;
	}
    }    
}