100 nCells += lduMatrices[i].
size();
105 convert(lduMatrices);
113 convert(ldum, interfaceCoeffs, interfaces);
118 pivotIndices_.setSize(n());
126 void Foam::LUscalarMatrix::convert
136 const scalar* __restrict__ diagPtr = ldum.
diag().
begin();
137 const scalar* __restrict__ upperPtr = ldum.
upper().
begin();
138 const scalar* __restrict__ lowerPtr = ldum.
lower().
begin();
140 register const label nCells = ldum.
diag().
size();
141 register const label nFaces = ldum.
upper().
size();
148 for (
register label face=0; face<nFaces; face++)
150 label uCell = uPtr[face];
151 label lCell = lPtr[face];
153 operator[](uCell)[lCell] = lowerPtr[face];
154 operator[](lCell)[uCell] = upperPtr[face];
159 if (interfaces.
set(inti))
161 const lduInterface&
interface = interfaces[inti].
interface();
163 const label* __restrict__ ulPtr =
interface.faceCells().begin();
164 const scalar* __restrict__ upperLowerPtr =
165 interfaceCoeffs[inti].
begin();
167 register label inFaces =
interface.faceCells().size()/2;
169 for (
register label face=0; face<inFaces; face++)
171 label uCell = ulPtr[face];
172 label lCell = ulPtr[face + inFaces];
174 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
175 operator[](lCell)[uCell] -= upperLowerPtr[face];
184 void Foam::LUscalarMatrix::convert
186 const PtrList<procLduMatrix>& lduMatrices
189 procOffsets_.setSize(lduMatrices.size() + 1);
192 forAll(lduMatrices, ldumi)
194 procOffsets_[ldumi+1] = procOffsets_[ldumi] + lduMatrices[ldumi].size();
197 forAll(lduMatrices, ldumi)
199 const procLduMatrix& lduMatrixi = lduMatrices[ldumi];
200 label offset = procOffsets_[ldumi];
202 const label* __restrict__ uPtr = lduMatrixi.upperAddr_.begin();
203 const label* __restrict__ lPtr = lduMatrixi.lowerAddr_.begin();
205 const scalar* __restrict__ diagPtr = lduMatrixi.diag_.begin();
206 const scalar* __restrict__ upperPtr = lduMatrixi.upper_.begin();
207 const scalar* __restrict__ lowerPtr = lduMatrixi.lower_.begin();
209 register const label nCells = lduMatrixi.
size();
210 register const label nFaces = lduMatrixi.upper_.size();
212 for (
register label cell=0; cell<nCells; cell++)
214 label globalCell = cell + offset;
215 operator[](globalCell)[globalCell] = diagPtr[cell];
218 for (
register label face=0; face<nFaces; face++)
220 label uCell = uPtr[face] + offset;
221 label lCell = lPtr[face] + offset;
223 operator[](uCell)[lCell] = lowerPtr[face];
224 operator[](lCell)[uCell] = upperPtr[face];
227 const PtrList<procLduInterface>& interfaces =
228 lduMatrixi.interfaces_;
232 const procLduInterface&
interface = interfaces[inti];
236 const label* __restrict__ ulPtr =
interface.faceCells_.begin();
238 const scalar* __restrict__ upperLowerPtr =
241 register label inFaces =
interface.faceCells_.size()/2;
243 for (
register label face=0; face<inFaces; face++)
245 label uCell = ulPtr[face] + offset;
246 label lCell = ulPtr[face + inFaces] + offset;
248 operator[](uCell)[lCell] -= upperLowerPtr[face + inFaces];
249 operator[](lCell)[uCell] -= upperLowerPtr[face];
254 const PtrList<procLduInterface>& neiInterfaces =
255 lduMatrices[
interface.neighbProcNo_].interfaces_;
257 label neiInterfacei = -1;
259 forAll(neiInterfaces, ninti)
263 neiInterfaces[ninti].neighbProcNo_
267 neiInterfacei = ninti;
272 if (neiInterfacei == -1)
277 const procLduInterface& neiInterface =
278 neiInterfaces[neiInterfacei];
280 const label* __restrict__ uPtr =
interface.faceCells_.begin();
281 const label* __restrict__ lPtr = neiInterface.faceCells_.begin();
283 const scalar* __restrict__ upperPtr =
interface.coeffs_.begin();
284 const scalar* __restrict__ lowerPtr = neiInterface.coeffs_.begin();
286 register label inFaces =
interface.faceCells_.size();
287 label neiOffset = procOffsets_[
interface.neighbProcNo_];
289 for (
register label face=0; face<inFaces; face++)
291 label uCell = uPtr[face] + offset;
292 label lCell = lPtr[face] + neiOffset;
294 operator[](uCell)[lCell] -= lowerPtr[face];
295 operator[](lCell)[uCell] -= upperPtr[face];
305 void Foam::LUscalarMatrix::printDiagonalDominance()
const
307 for (label i=0; i<n(); i++)
310 for (label j=0; j<n(); j++)
314 sum += operator[](i)[j];