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