133 nCells += lduMatrices[i].size();
138 convert(lduMatrices);
146 convert(ldum, interfaceCoeffs, interfaces);
153 const label numRows = m();
154 const label numCols =
n();
156 Pout<<
"LUscalarMatrix : size:" << numRows <<
endl;
157 for (label rowi = 0; rowi < numRows; ++rowi)
159 const scalar* row = operator[](rowi);
161 Pout<<
"cell:" << rowi <<
" diagCoeff:" << row[rowi] <<
endl;
163 Pout<<
" connects to upper cells :";
164 for (label coli = rowi+1; coli < numCols; ++coli)
166 if (
mag(row[coli]) > SMALL)
168 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
172 Pout<<
" connects to lower cells :";
173 for (label coli = 0; coli < rowi; ++coli)
175 if (
mag(row[coli]) > SMALL)
177 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
185 pivotIndices_.setSize(m());
193 void Foam::LUscalarMatrix::convert
203 const scalar* __restrict__ diagPtr = ldum.
diag().begin();
204 const scalar* __restrict__ upperPtr = ldum.
upper().begin();
205 const scalar* __restrict__ lowerPtr = ldum.
lower().begin();
207 const label nCells = ldum.
diag().size();
208 const label nFaces = ldum.
upper().size();
215 for (label face=0; face<nFaces; face++)
217 label uCell = uPtr[face];
218 label lCell = lPtr[face];
220 operator[](uCell)[lCell] = lowerPtr[face];
221 operator[](lCell)[uCell] = upperPtr[face];
226 if (interfaces.
set(inti))
228 const lduInterface&
interface = interfaces[inti].
interface();
232 const label* __restrict__ lPtr =
interface.faceCells().begin();
234 const cyclicLduInterface& cycInterface =
235 refCast<const cyclicLduInterface>(
interface);
236 label nbrInt = cycInterface.neighbPatchID();
237 const label* __restrict__ uPtr =
238 interfaces[nbrInt].interface().faceCells().
begin();
240 const scalar* __restrict__ nbrUpperLowerPtr =
241 interfaceCoeffs[nbrInt].begin();
243 label inFaces =
interface.faceCells().size();
245 for (label face=0; face<inFaces; face++)
247 label uCell = lPtr[face];
248 label lCell = uPtr[face];
250 operator[](uCell)[lCell] -= nbrUpperLowerPtr[face];
257 void Foam::LUscalarMatrix::convert
259 const PtrList<procLduMatrix>& lduMatrices
262 procOffsets_.setSize(lduMatrices.size() + 1);
265 forAll(lduMatrices, ldumi)
267 procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
270 forAll(lduMatrices, ldumi)
272 const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
273 label offset = procOffsets_[ldumi];
275 const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
276 const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
278 const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
279 const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
280 const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
282 const label nCells = lduMatrixi.size();
283 const label nFaces = lduMatrixi.upper_.size();
285 for (label cell=0; cell<nCells; cell++)
287 label globalCell = cell + offset;
288 operator[](globalCell)[globalCell] = diagPtr[cell];
291 for (label face=0; face<nFaces; face++)
293 label uCell = uPtr[face] + offset;
294 label lCell = lPtr[face] + offset;
296 operator[](uCell)[lCell] = lowerPtr[face];
297 operator[](lCell)[uCell] = upperPtr[face];
300 const PtrList<procLduInterface>& interfaces =
301 lduMatrixi.interfaces_;
305 const procLduInterface&
interface = interfaces[inti];
309 const label* __restrict__ ulPtr =
interface.faceCells_.begin();
311 const scalar* __restrict__ upperLowerPtr =
314 label inFaces =
interface.faceCells_.size()/2;
316 for (label face=0; face<inFaces; face++)
318 label uCell = ulPtr[face] + offset;
319 label lCell = ulPtr[face + inFaces] + offset;
321 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
322 operator[](lCell)[uCell] -= upperLowerPtr[face];
332 const PtrList<procLduInterface>& neiInterfaces =
333 lduMatrices[
interface.neighbProcNo_].interfaces_;
335 label neiInterfacei = -1;
337 forAll(neiInterfaces, ninti)
342 neiInterfaces[ninti].neighbProcNo_
345 && (neiInterfaces[ninti].tag_ ==
interface.tag_)
348 neiInterfacei = ninti;
353 if (neiInterfacei == -1)
358 const procLduInterface& neiInterface =
359 neiInterfaces[neiInterfacei];
361 const label* __restrict__ uPtr =
interface.faceCells_.begin();
362 const label* __restrict__ lPtr =
363 neiInterface.faceCells_.begin();
365 const scalar* __restrict__ upperPtr =
interface.coeffs_.begin();
366 const scalar* __restrict__ lowerPtr =
367 neiInterface.coeffs_.begin();
369 label inFaces =
interface.faceCells_.size();
370 label neiOffset = procOffsets_[
interface.neighbProcNo_];
372 for (label face=0; face<inFaces; face++)
374 label uCell = uPtr[face] + offset;
375 label lCell = lPtr[face] + neiOffset;
377 operator[](uCell)[lCell] -= lowerPtr[face];
378 operator[](lCell)[uCell] -= upperPtr[face];
386 void Foam::LUscalarMatrix::printDiagonalDominance()
const
388 for (label i=0; i<m(); i++)
391 for (label j=0; j<m(); j++)
395 sum += operator[](i)[j];
406 pivotIndices_.setSize(m());
415 for (label j=0; j<m(); j++)
420 for (label i=0; i<m(); i++)