407 long long rows = (*Matrix_).NumGlobalRows64();
408 long long cols = (*Matrix_).NumGlobalCols64();
409 int num_edges = ((*Matrix_).NumMyNonzeros() - (*Matrix_).NumMyDiagonals())/2;
410 std::cout <<
"global num rows " << rows << std::endl;
415 if(rows > std::numeric_limits<int>::max())
417 std::cerr <<
"Ifpack_SupportGraph<T>::FindSupport: global num rows won't fit an int. " << rows << std::endl;
422 int num_verts = (int) rows;
425 E *edge_array =
new E[num_edges];
426 double *weights =
new double[num_edges];
429 int max_num_entries = (*Matrix_).MaxNumEntries();
430 double *values =
new double[max_num_entries];
431 int *indices =
new int[max_num_entries];
433 double * diagonal =
new double[num_verts];
436 for(
int i = 0; i < max_num_entries; i++)
444 for(
int i = 0; i < num_verts; i++)
446 (*Matrix_).ExtractMyRowCopy(i,max_num_entries,num_entries,values,indices);
448 for(
int j = 0; j < num_entries; j++)
453 diagonal[i] = values[j];
456 diagonal[i] *= DiagPertRel_;
458 diagonal[i] += DiagPertAbs_;
463 edge_array[k] =
E(i,indices[j]);
464 weights[k] = values[j];
468 weights[k] *= (1.0 + 1e-8 * drand48());
477 Graph g(edge_array, edge_array + num_edges, weights, num_verts);
480 property_map < Graph, edge_weight_t >::type weight = get(edge_weight, g);
482 std::vector < Edge > spanning_tree;
485 kruskal_minimum_spanning_tree(g, std::back_inserter(spanning_tree));
488 std::vector<int> NumNz(num_verts,1);
491 for (std::vector < Edge >::iterator ei = spanning_tree.begin();
492 ei != spanning_tree.end(); ++ei)
494 NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
495 NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
500 std::vector< std::vector< int > > Indices(num_verts);
505 std::vector< std::vector< double > > Values(num_verts);
507 for(
int i = 0; i < num_verts; i++)
509 std::vector<int> temp(NumNz[i],0);
510 std::vector<double> temp2(NumNz[i],0);
515 int *l =
new int[num_verts];
516 for(
int i = 0; i < num_verts; i++)
524 for(
int i = 0; i < NumForests_; i++)
528 spanning_tree.clear();
529 kruskal_minimum_spanning_tree(g,std::back_inserter(spanning_tree));
530 for(std::vector < Edge >::iterator ei = spanning_tree.begin();
531 ei != spanning_tree.end(); ++ei)
533 NumNz[source(*ei,g)] = NumNz[source(*ei,g)] + 1;
534 NumNz[target(*ei,g)] = NumNz[target(*ei,g)] + 1;
536 for(
int i = 0; i < num_verts; i++)
538 Indices[i].resize(NumNz[i]);
539 Values[i].resize(NumNz[i]);
543 for (std::vector < Edge >::iterator ei = spanning_tree.begin();
544 ei != spanning_tree.end(); ++ei)
548 Indices[source(*ei,g)][0] = source(*ei,g);
549 Values[source(*ei,g)][0] = Values[source(*ei,g)][0] - weight[*ei];
550 Indices[target(*ei,g)][0] = target(*ei,g);
551 Values[target(*ei,g)][0] = Values[target(*ei,g)][0] - weight[*ei];
553 Indices[source(*ei,g)][l[source(*ei,g)]] = target(*ei,g);
554 Values[source(*ei,g)][l[source(*ei,g)]] = weight[*ei];
555 l[source(*ei,g)] = l[source(*ei,g)] + 1;
557 Indices[target(*ei,g)][l[target(*ei,g)]] = source(*ei,g);
558 Values[target(*ei,g)][l[target(*ei,g)]] = weight[*ei];
559 l[target(*ei,g)] = l[target(*ei,g)] + 1;
576 Matrix_->Multiply(
false, ones, surplus);
578 for(
int i = 0; i < num_verts; i++)
580 Values[i][0] += surplus[i];
581 Values[i][0] = KeepDiag_*diagonal[i] +
582 (1.-KeepDiag_) * Values[i][0];
588#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
589 if((*Matrix_).RowMatrixRowMap().GlobalIndicesLongLong())
592 for(
int i = 0; i < num_verts; i++)
594 std::vector<long long> IndicesLL(l[i]);
595 for(
int k = 0; k < l[i]; ++k)
596 IndicesLL[k] = Indices[i][k];
598 (*Support_).InsertGlobalValues(i,l[i],&Values[i][0],&IndicesLL[0]);
603#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
604 if((*Matrix_).RowMatrixRowMap().GlobalIndicesInt())
607 for(
int i = 0; i < num_verts; i++)
609 (*Support_).InsertGlobalValues(i,l[i],&Values[i][0],&Indices[i][0]);
614 throw "Ifpack_SupportGraph::FindSupport: GlobalIndices unknown.";;
616 (*Support_).FillComplete();
739Print(std::ostream& os)
const
741 os <<
"================================================================================" << std::endl;
742 os <<
"Ifpack_SupportGraph: " << Label () << std::endl << std::endl;
743 os <<
"Condition number estimate = " << Condest() << std::endl;
744 os <<
"Global number of rows = " << Matrix_->NumGlobalRows64() << std::endl;
745 os <<
"Number of edges in support graph = " << (Support_->NumGlobalNonzeros64()-Support_->NumGlobalDiagonals64())/2 << std::endl;
746 os <<
"Fraction of off diagonals of support graph/off diagonals of original = "
747 << ((double)Support_->NumGlobalNonzeros64()-Support_->NumGlobalDiagonals64())/(Matrix_->NumGlobalNonzeros64()-Matrix_->NumGlobalDiagonals64());
749 os <<
"Phase # calls Total Time (s) Total MFlops MFlops/s" << std::endl;
750 os <<
"----- ------- -------------- ------------ --------" << std::endl;
751 os <<
"Initialize() " << std::setw(10) << NumInitialize_
752 <<
" " << std::setw(15) << InitializeTime_
753 <<
" 0.0 0.0" << std::endl;
754 os <<
"Compute() " << std::setw(10) << NumCompute_
755 <<
" " << std::setw(22) << ComputeTime_
756 <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_;
757 if (ComputeTime_ != 0.0)
758 os <<
" " << std::setw(15) << 1.0e-6 * ComputeFlops_ / ComputeTime_ << std::endl;
760 os <<
" " << std::setw(15) << 0.0 << std::endl;
761 os <<
"ApplyInverse() " << std::setw(10) << NumApplyInverse_
762 <<
" " << std::setw(22) << ApplyInverseTime_
763 <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_;
764 if (ApplyInverseTime_ != 0.0)
765 os <<
" " << std::setw(15) << 1.0e-6 * ApplyInverseFlops_ / ApplyInverseTime_ << std::endl;
767 os <<
" " << std::setw(15) << 0.0 << std::endl;
769 os << std::endl << std::endl;
770 os <<
"Now calling the underlying preconditioner's print()" << std::endl;