137 nCells += lduMatrices[i].size();
142 convert(lduMatrices);
150 convert(ldum, interfaceCoeffs, interfaces);
157 const label numRows = m();
158 const label numCols =
n();
160 Pout<<
"LUscalarMatrix : size:" << numRows <<
endl;
161 for (label rowi = 0; rowi < numRows; ++rowi)
163 const scalar* row = operator[](rowi);
165 Pout<<
"cell:" << rowi <<
" diagCoeff:" << row[rowi] <<
endl;
167 Pout<<
" connects to upper cells :";
168 for (label coli = rowi+1; coli < numCols; ++coli)
170 if (
mag(row[coli]) > SMALL)
172 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
176 Pout<<
" connects to lower cells :";
177 for (label coli = 0; coli < rowi; ++coli)
179 if (
mag(row[coli]) > SMALL)
181 Pout<<
' ' << coli <<
" (coeff:" << row[coli] <<
')';
189 pivotIndices_.setSize(m());
197 void Foam::LUscalarMatrix::convert
207 const scalar* __restrict__ diagPtr = ldum.
diag().begin();
208 const scalar* __restrict__ upperPtr = ldum.
upper().begin();
209 const scalar* __restrict__ lowerPtr = ldum.
lower().begin();
211 const label nCells = ldum.
diag().size();
212 const label nFaces = ldum.
upper().size();
219 for (label face=0; face<nFaces; face++)
221 label uCell = uPtr[face];
222 label lCell = lPtr[face];
224 operator[](uCell)[lCell] = lowerPtr[face];
225 operator[](lCell)[uCell] = upperPtr[face];
230 if (interfaces.
set(inti))
232 const lduInterface&
interface = interfaces[inti].
interface();
236 const label* __restrict__ lPtr =
interface.faceCells().begin();
238 const cyclicLduInterface& cycInterface =
239 refCast<const cyclicLduInterface>(
interface);
240 label nbrInt = cycInterface.neighbPatchID();
241 const label* __restrict__ uPtr =
242 interfaces[nbrInt].interface().faceCells().
begin();
244 const scalar* __restrict__ nbrUpperLowerPtr =
245 interfaceCoeffs[nbrInt].begin();
247 label inFaces =
interface.faceCells().size();
249 for (label face=0; face<inFaces; face++)
251 label uCell = lPtr[face];
252 label lCell = uPtr[face];
254 operator[](uCell)[lCell] -= nbrUpperLowerPtr[face];
261 void Foam::LUscalarMatrix::convert
263 const PtrList<procLduMatrix>& lduMatrices
266 procOffsets_.setSize(lduMatrices.size() + 1);
269 forAll(lduMatrices, ldumi)
271 procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
274 forAll(lduMatrices, ldumi)
276 const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
277 label offset = procOffsets_[ldumi];
279 const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
280 const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
282 const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
283 const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
284 const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
286 const label nCells = lduMatrixi.size();
287 const label nFaces = lduMatrixi.upper_.size();
289 for (label cell=0; cell<nCells; cell++)
291 label globalCell = cell + offset;
292 operator[](globalCell)[globalCell] = diagPtr[cell];
295 for (label face=0; face<nFaces; face++)
297 label uCell = uPtr[face] + offset;
298 label lCell = lPtr[face] + offset;
300 operator[](uCell)[lCell] = lowerPtr[face];
301 operator[](lCell)[uCell] = upperPtr[face];
304 const PtrList<procLduInterface>& interfaces =
305 lduMatrixi.interfaces_;
309 const procLduInterface&
interface = interfaces[inti];
313 const label* __restrict__ ulPtr =
interface.faceCells_.begin();
315 const scalar* __restrict__ upperLowerPtr =
318 label inFaces =
interface.faceCells_.size()/2;
320 for (label face=0; face<inFaces; face++)
322 label uCell = ulPtr[face] + offset;
323 label lCell = ulPtr[face + inFaces] + offset;
325 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
326 operator[](lCell)[uCell] -= upperLowerPtr[face];
336 const PtrList<procLduInterface>& neiInterfaces =
337 lduMatrices[
interface.neighbProcNo_].interfaces_;
339 label neiInterfacei = -1;
341 forAll(neiInterfaces, ninti)
346 neiInterfaces[ninti].neighbProcNo_
349 && (neiInterfaces[ninti].tag_ ==
interface.tag_)
352 neiInterfacei = ninti;
357 if (neiInterfacei == -1)
362 const procLduInterface& neiInterface =
363 neiInterfaces[neiInterfacei];
365 const label* __restrict__ uPtr =
interface.faceCells_.begin();
366 const label* __restrict__ lPtr =
367 neiInterface.faceCells_.begin();
369 const scalar* __restrict__ upperPtr =
interface.coeffs_.begin();
370 const scalar* __restrict__ lowerPtr =
371 neiInterface.coeffs_.begin();
373 label inFaces =
interface.faceCells_.size();
374 label neiOffset = procOffsets_[
interface.neighbProcNo_];
376 for (label face=0; face<inFaces; face++)
378 label uCell = uPtr[face] + offset;
379 label lCell = lPtr[face] + neiOffset;
381 operator[](uCell)[lCell] -= lowerPtr[face];
382 operator[](lCell)[uCell] -= upperPtr[face];
390 void Foam::LUscalarMatrix::printDiagonalDominance()
const
392 for (label i=0; i<m(); i++)
395 for (label j=0; j<m(); j++)
399 sum += operator[](i)[j];
410 pivotIndices_.setSize(m());
419 for (label j=0; j<m(); j++)
424 for (label i=0; i<m(); i++)