/*
 $Id: sysmat.cc,v 1.2 1996/10/04 15:07:38 roitzsch Exp $
 (C)opyright 1996 by Konrad-Zuse-Center, Berlin
 All rights reserved.
 Part of the Kaskade distribution
*/

#include "sysmat.h"
#include "numerics.h"

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

SystemMatrix:: SystemMatrix(symmetryType symmetry0) : symmetry(symmetry0) { }

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

void SystemMatrix:: notImplemented(const char* s)  const
{
    cout << "\n*** Class System Matrix: called function " 
    					<< s << " not implemented\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

void SystemMatrix:: callError(const char* s) const
{
    cout << "\n*** Error in calling " << s 
        	<< ": operation may not be done in place !\n";
    cout.flush(); abort();
}
//-------------------------------------------------------------------------

		// ... lots of dummy(!) definitions:

SystemMatrix:: SystemMatrix(SystemMatrix* /*A*/,SystemMatrix* /*B*/,Num /*lambda*/)
    : symmetry(sym)  		      { notImplemented("SystemMatrix"); }

Num& SystemMatrix:: Diag(int /*i*/)	      
		{ static Num dummy; notImplemented("Diag"); return dummy; }
void SystemMatrix:: Mult(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					{ notImplemented("Mult");}
void SystemMatrix:: ATMult(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					{ notImplemented("ATMult");}
void SystemMatrix:: smoothNode(int /*no*/, Vector<Num>& /*e*/, const Vector<Num>& /*r*/)
					{ notImplemented("smoothNode");}
Num  SystemMatrix:: MultRow(int /*row*/, const Vector<Num>& /*e*/)
                                        { notImplemented("MultRow"); return 0; }
void SystemMatrix:: DiagMult(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/) 
					{ notImplemented("DiagMult"); }
void SystemMatrix:: DiagDiv (Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)
			    		 { notImplemented("DiagDiv"); }
void SystemMatrix:: LMult (Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					{ notImplemented("LMult");}
void SystemMatrix:: UMult (Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					{ notImplemented("UMult");}
void SystemMatrix:: F  (Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)
					{ notImplemented("F");}
void SystemMatrix:: FT (Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)	
					{ notImplemented("FT");}
void SystemMatrix:: Fm1(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)
						{ notImplemented("Fm1");}
void SystemMatrix:: FmT(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/, Real /*omega*/)	
						{ notImplemented("Fmt");}
void SystemMatrix:: removeLU() {  notImplemented("removeLU"); }
void SystemMatrix:: symmetricExpansion() {notImplemented("symmetricExpansion");}
void SystemMatrix:: removeSymmetricExpansion() 
			       {notImplemented("removeSymmetricExpansion");}
void SystemMatrix:: print()  const	{ notImplemented("print"); }
void SystemMatrix:: printMatLabFormat()  const  
					{ notImplemented("printMatLabFormat"); }
//-------------------------------------------------------------------------

void SystemMatrix:: Decompose(Bool /*removeOriginal*/)   
					      { notImplemented("Decompose"); }
void SystemMatrix:: FBSubst(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					      { notImplemented("FBSubst"); }
void SystemMatrix:: ILUDecompose(Real /*dropTol*/){ notImplemented("ILUDecompose"); }
void SystemMatrix:: ILUFBSubst(Vector<Num>& /*lhs*/, Vector<Num>& /*rhs*/)
					      { notImplemented("ILUFBSubst"); }
//-------------------------------------------------------------------------

Real SystemMatrix:: omegaOpt(int /*mode*/) const { return 1.0; }

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

void SystemMatrix:: checkDiagonal()
{
    int i;
    Real d;
    Real maxD   =  ::real(Diag(1));
    Real minD   =  ::real(Diag(1));
    Real minAbs =  Abs(Diag(1));

    for (i=2; i<=Dim(); ++i)
    {
	d = ::real(Diag(i));

	if (d < minD) minD = d;
	if (d > maxD) maxD = d;

	d = Abs(d);
	if (d < minAbs) minAbs = d;

	//if (d <= 0.0) {
	//    cout << "\n*** SparseMatrix: Diagonal element of A <= 0:  D[" 
	//    << i << "] = " << D[i] << "\n"; 
	//    return; }
    }
    cout << "\n    --- Matrix Diagonal:  min=" << minD << "  max=" << maxD 
				     << "  min(abs)=" << minAbs << "\n";
}
//-------------------------------------------------------------------------


void SystemMatrix:: store(const Matrix<Num>& AElem, 
			  const Vector<int>& globalNodes)
{
    int i, j, row, col, jMax;
    const int dim = globalNodes.high();

    for (i=1; i<=dim; ++i)
    {
	if ((row=globalNodes[i]))
	{
	    jMax = (symmetry==sym)? i : dim;
	
	    for (j=1; j<=jMax; ++j)
	      if ((col=globalNodes[j])) (*this)(row,col) += AElem(i,j);
	}
    }
}