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