/*
$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;
}
}
*/
//-------------------------------------------------------------------------
syntax highlighted by Code2HTML, v. 0.9.1