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

#include "sysmatma28.h"
#include "utils.h"
#include "numerics.h"

#include "cmdpars.h"
extern CmdPars Cmd;

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

extern "C" 
{ 
    void F77NAME(ma28adint,MA28ADINT)(int* N, int* NZ, double* A, int* LICN, int* IRN, 
			    int* LIRN, int* ICN, double* U, int* IKEEP, 
			    int* IW, double* W, int* FLAG, int* nsearch, 
			    int* rncp, int* cncp, int* iMinIRN, int* iMinICN, 
			    double* TOL, int* NDROP);
    void F77NAME(ma28cdint,MA28CDINT)(int* N, double* A, int* LICN, int* ICN,
			    int* IKEEP, double* rhs, double* W, int* MTYPE);

    void F77NAME(madc28adint,MAC28ADINT)(int* N, int* NZ, Complex* A, int* LICN, int* IRN, 
                              int* LIRN, int* ICN, double* U, int* IKEEP,
			      int* IW, Complex* W, int* FLAG, int* nsearch, 
			      int* rncp, int* cncp, int* iMinIRN, int* iMinICN,
			      double* TOL, int* NDROP);
    void F77NAME(madc28cdint,MADC28CDINT)(int* N, Complex* A, int* LICN, int* ICN,
                              int* IKEEP, Complex* rhs, Complex* W, int* MTYPE);
};
//-------------------------------------------------------------------------
//-------------------------------------------------------------------------


MA28Matrix:: MA28Matrix(int spaceDim0)

	: spaceDim(spaceDim0),
	  A(0,0), W(0,0), IRN(0,0), ICN(0,0), IKEEP(0,0), IW(0,0)
{ 
    NZ = LICN = LIRN = 0;

    defDropTol = 1e-8;		Cmd.get("defDropTol", &defDropTol);

    UPivot = 0.1; 		Cmd.get("UPivot", &UPivot);

    ICNFactor = 1.0;		Cmd.get("ICNFactor", &ICNFactor);
    IRNFactor = 1.0;		Cmd.get("IRNFactor", &IRNFactor);
    NSearch   = 1;		Cmd.get("nsearch",   &NSearch);

    infoLinSystem = 0;       	Cmd.get("infoLinSystem",  &infoLinSystem);
}
//-------------------------------------------------------------------------


void MA28Matrix:: MA28RemoveLUFactorization()
{
    A.resize(0,0); 
    W.resize(0,0);
    IW.resize(0,0); 
    IRN.resize(0,0); 
    ICN.resize(0,0); 
    IKEEP.resize(0,0); 
}
//-------------------------------------------------------------------------


void MA28Matrix:: MA28Decompose(Real dropTol)
{
    int rncp, cncp, MINIRN, MINICN, NDROP, FLAG, extend, maxExtend = 3;


    countEntries(&N, &NZ);
    MA28ResizeArrays(dropTol);
    fillMA28Vectors();

    if (equal(dropTol, 0.0)) dropTol = defDropTol;
    Real absTol = dropTol * minDiag();    

    
    for (extend=0; extend<=maxExtend; ++extend)
    {
	#if (ComplexVersion)

	F77NAME(madc28adint,MADC28ADINT)(&N, &NZ, A.v, &LICN, IRN.v, &LIRN, ICN.v,
			     &UPivot, IKEEP.v, IW.v, W.v, &FLAG, &NSearch, 
			     &rncp, &cncp, &MINIRN, &MINICN, &absTol, &NDROP);
	#else

	F77NAME(ma28adint,MA28ADINT)(&N, &NZ, A.v, &LICN, IRN.v, &LIRN, ICN.v,
			   &UPivot, IKEEP.v, IW.v, W.v, &FLAG, &NSearch, 
			   &rncp, &cncp, &MINIRN, &MINICN, &absTol, &NDROP);
	#endif

	if (infoLinSystem)
	{
	    double space = double(MA28MemSpace())/1e6;
	    cout << Form("\n    MA28 Decomposition: FLAG = %1d", FLAG);
	    cout << Form(" N = %1d Entries = %1d (%1.6f MB)", N, NZ, space);
	}

	if (FLAG == 0) break;		      // successful decomposition

	else if (FLAG <= -3 && FLAG >= -6)    // space in IRN/ICN not sufficient
	{
	    MA28ExtendArrays(FLAG); 
	    fillMA28Vectors();
	}    
	else 
	{
	    cout << "\n*** MA28Decompose failed: FLAG = " << FLAG << "\n";
	    cout.flush(); abort();
	}
    }

    if (extend > maxExtend)
    {
	cout << "\n*** MA28Decompose failed. Increase ICNFact and/or IRNFact!\n";
	cout.flush(); abort();
    }

    IRN.resize(0,0); 

    if (infoLinSystem)
    {
	cout << Form("\n    U = %1.2f  NSearch = %1d", UPivot, NSearch);
	cout << Form("  irncp=%1d  icncp=%1d  fill-in: %1d (f=%1.2f)", 
		       rncp, cncp, MINICN-NZ, Real(MINICN-NZ)/Real(NZ));
	cout << Form("\n    MINICN=%1d (%1.2f LICN=%1d)", 
		       MINICN, Real(MINICN)/LICN, LICN);
	cout << Form("  MINIRN=%1d (%1.2f LIRN=%1d)\n", 
		       MINIRN, Real(MINIRN)/LIRN, LIRN);
	cout << Form("    dropTol = %5.3g (absTol = %5.3g)   NDROP = %1d\n",
		       dropTol, absTol, NDROP);
    }
}
//-------------------------------------------------------------------------


void MA28Matrix:: MA28ResizeArrays(Real dropTol)
{
    W.resize(0,N);
    IKEEP.resize(0,5*N);

    if (NSearch <= N) IW.resize(0,7*N);
    else 	      IW.resize(0,8*N);

    Real fact = 0.05;

    if (spaceDim > 1)
    {
        if      (dropTol > 1e-1) fact *= 0.1;
        else if (dropTol > 1e-2) fact *= 0.2;
        else if (dropTol > 1e-3) fact *= 0.3;
        else if (dropTol > 1e-4) fact *= 0.4;
        else if (dropTol > 1e-5) fact *= 0.5;
        else if (dropTol > 1e-6) fact *= 0.6;
        else if (dropTol > 1e-7) fact *= 0.7;
    }

    if (spaceDim == 1) 
    {
	LICN = int(ICNFactor * NZ * 1.5);
	LIRN = int(IRNFactor * NZ * 1.2);
    }
    else if (spaceDim == 2) 
    {
	LICN = int(ICNFactor * NZ * (4.0 + fact*int(pow(double(N),0.50))));
	LIRN = int(IRNFactor * NZ * (2.0 + fact*int(pow(double(N),0.25))));
    }
    else if (spaceDim == 3) 
    {
	LICN = int(ICNFactor * NZ * (6.0 + fact*int(pow(double(N),0.66))));
	LIRN = int(IRNFactor * NZ * (3.0 + fact*int(pow(double(N),0.50))));
    }


    A.resize  (0,LICN);
    ICN.resize(0,LICN);
    IRN.resize(0,LIRN);
}
//-------------------------------------------------------------------------


void MA28Matrix:: MA28ExtendArrays(int FLAG)
{
    if (FLAG == -3 || FLAG == -6) 
    {  
	LIRN = int(1.5*LIRN); 
	IRN.resize(0,LIRN); 
	if (infoLinSystem) cout << "\n\t * MA28 FLAG " << FLAG 
	    			<< ": IRN Working Space extended\n";
    }
    if (FLAG == -5 || FLAG == -6) 
    {
	LICN = int(1.5*LICN); 
	ICN.resize(0,LICN); 
	A.resize  (0,LICN);
	if (infoLinSystem) cout << "\n\t * MA28 FLAG " << FLAG 
	    			<< ": ICN Working Space extended\n";
    }
    if (FLAG == -4)
    {
	LICN = int(3*LICN); 
	ICN.resize(0,LICN); 
	A.resize  (0,LICN);
	if (infoLinSystem) cout << "\n\t * MA28 FLAG " << FLAG 
	    		    << ": ICN Working Space extended (far too small)\n";

    }
}
//-------------------------------------------------------------------------


void MA28Matrix:: MA28FBSubst(Num* rhs)
{
    int MTYPE = 1;

    #if (ComplexVersion)
        F77NAME(madc28cdint,MADC28CDINT)(&N, A.v, &LICN, ICN.v, IKEEP.v, rhs, W.v, &MTYPE);
    #else
        F77NAME(ma28cdint,MA28CDINT)  (&N, A.v, &LICN, ICN.v, IKEEP.v, rhs, W.v, &MTYPE);
    #endif
}
//-------------------------------------------------------------------------


Real MA28Matrix:: minDiag()
{
    int i;
    Real minD = Abs(A[0]);
    for (i=1; i<N; ++i) minD = Min(minD,Abs(A[i]));
    return minD;
}
//-------------------------------------------------------------------------


int MA28Matrix:: MA28MemSpace(int /*print*/) const
{
    int space = A.memSpace() + W.memSpace() + IRN.memSpace() + ICN.memSpace() + 
   	        IKEEP.memSpace() + IW.memSpace();
    return space;
}
//-------------------------------------------------------------------------

/*
   void MA28Matrix:: sortRows()
   {
   int i, n, k;
   
   for (i=2; i<=Dim(); ++i) 
   {
   n = end[i]-end[i-1];
   k = end[i-1];
   sortRow(n, &(col.v[k]), &(L->v[k]));
   }
   }
   //-------------------------------------------------------------------------
   
   // -- sort by straight insertion; first element in vector is j=1
   
   void MA28Matrix:: sortRow(int n, int* col, Num* L)  
   {
   int i,j;
   int a;
   Num u;
   
   for (j=2; j<=n; j++) 
   {
   a = col[j];
   u = L[j];
   
   i = j-1;
   
   while (i > 0 && col[i] > a) 
   {
   col[i+1] = col[i];
   L  [i+1] = L  [i];
   --i;
   }
   col[i+1] = a;
   L  [i+1] = u;
   }
}
*/
//-------------------------------------------------------------------------