182 int Length =
A_.MaxNumEntries();
183 std::vector<double> RowValuesL(Length);
184 std::vector<int> RowIndicesU(Length);
185 std::vector<int_type> RowIndicesU_LL;
186 RowIndicesU_LL.resize(Length);
187 std::vector<double> RowValuesU(Length);
194 if ((
L_.get() == 0) || (
U_.get() == 0))
199 &RowValuesU[0],&RowIndicesU[0]));
201 bool distributed = (
Comm().NumProc() > 1)?
true:
false;
206 for (
int i = 0 ;i < RowNnzU ; ++i)
209 RowIndicesU[count] = RowIndicesU[i];
210 RowValuesU[count] = RowValuesU[i];
220 for (
int i = 0 ;i < RowNnzU ; ++i) {
221 if (RowIndicesU[i] == 0) {
222 double& v = RowValuesU[i];
228 std::copy(&(RowIndicesU[0]), &(RowIndicesU[0]) + RowNnzU, RowIndicesU_LL.begin());
230 &(RowIndicesU_LL[0])));
234 int_type RowIndicesU_0_templ = RowIndicesU[0];
236 &RowIndicesU_0_templ));
238 int max_keys = (int) 10 *
A_.MaxNumEntries() *
LevelOfFill() ;
242 int hash_size = SingleRowU.getRecommendedHashSize(max_keys) ;
243 std::vector<int> keys; keys.reserve(hash_size * 10);
244 std::vector<double> values; values.reserve(hash_size * 10);
245 std::vector<double> AbsRow; AbsRow.reserve(hash_size * 10);
251#ifdef IFPACK_FLOPCOUNTERS
252 double this_proc_flops = 0.0;
255 for (
int row_i = 1 ; row_i <
NumMyRows_ ; ++row_i)
259 &RowValuesU[0],&RowIndicesU[0]));
264 for (
int i = 0 ;i < RowNnzU ; ++i)
267 RowIndicesU[count] = RowIndicesU[i];
268 RowValuesU[count] = RowValuesU[i];
280 for (
int i = 0 ;i < RowNnzU ; ++i) {
281 if (RowIndicesU[i] < row_i)
283 else if (RowIndicesU[i] == row_i) {
286 double& v = RowValuesU[i];
295 if (FillL == 0) FillL = 1;
296 if (FillU == 0) FillU = 1;
301 for (
int i = 0 ; i < RowNnzU ; ++i) {
302 SingleRowU.set(RowIndicesU[i], RowValuesU[i]);
309 for (
int i = 0 ; i < RowNnzU ; ++i)
310 start_col =
EPETRA_MIN(start_col, RowIndicesU[i]);
312#ifdef IFPACK_FLOPCOUNTERS
316 for (
int col_k = start_col ; col_k < row_i ; ++col_k) {
318 int_type* ColIndicesK;
322 double xxx = SingleRowU.get(col_k);
330 double DiagonalValueK = 0.0;
331 for (
int i = 0 ; i < ColNnzK ; ++i) {
332 if (ColIndicesK[i] == col_k) {
333 DiagonalValueK = ColValuesK[i];
339 if (DiagonalValueK == 0.0)
342 double yyy = xxx / DiagonalValueK ;
343 SingleRowL.set(col_k, yyy);
344#ifdef IFPACK_FLOPCOUNTERS
348 for (
int j = 0 ; yyy != 0.0 && j < ColNnzK ; ++j)
350 int_type col_j = ColIndicesK[j];
352 if (col_j < col_k)
continue;
354 SingleRowU.set(col_j, -yyy * ColValuesK[j],
true);
355#ifdef IFPACK_FLOPCOUNTERS
362#ifdef IFPACK_FLOPCOUNTERS
363 this_proc_flops += 1.0 * flops;
367 double DiscardedElements = 0.0;
372 int sizeL = SingleRowL.getNumEntries();
374 values.resize(sizeL);
376 AbsRow.resize(sizeL);
379 keys.size() ? &keys[0] : 0,
380 values.size() ? &values[0] : 0
382 for (
int i = 0; i < sizeL; ++i)
388 nth_element(AbsRow.begin(), AbsRow.begin() + FillL, AbsRow.begin() + count,
389 std::greater<double>());
390 cutoff = AbsRow[FillL];
393 for (
int i = 0; i < sizeL; ++i) {
395 int_type key_templ = keys[i];
399 DiscardedElements += values[i];
405 int_type row_i_templ = row_i;
410 int sizeU = SingleRowU.getNumEntries();
411 AbsRow.resize(sizeU + 1);
412 keys.resize(sizeU + 1);
413 values.resize(sizeU + 1);
415 SingleRowU.arrayify(&keys[0], &values[0]);
417 for (
int i = 0; i < sizeU; ++i)
424 nth_element(AbsRow.begin(), AbsRow.begin() + FillU, AbsRow.begin() + count,
425 std::greater<double>());
426 cutoff = AbsRow[FillU];
430 for (
int i = 0; i < sizeU; ++i)
433 double val = values[i];
436 if (
IFPACK_ABS(val) >= cutoff || row_i == col) {
437 int_type col_templ = col;
441 DiscardedElements += val;
453#ifdef IFPACK_FLOPCOUNTERS
455 Comm().SumAll(&this_proc_flops, &tf, 1);
470 cout <<
"A = " <<
A_.NumGlobalNonzeros() << endl;
471 cout <<
"L = " <<
L_->NumGlobalNonzeros() << endl;
472 cout <<
"U = " <<
U_->NumGlobalNonzeros() << endl;
474 A_.Multiply(
false,LHS,RHS1);
475 U_->Multiply(
false,LHS,RHS2);
476 L_->Multiply(
false,RHS2,RHS3);
478 RHS1.Update(-1.0, RHS3, 1.0);
483 long long MyNonzeros =
L_->NumGlobalNonzeros64() +
U_->NumGlobalNonzeros64();