47 #include "Teuchos_ParameterList.hpp" 48 #include "Teuchos_RCP.hpp" 49 #include "Epetra_Import.h" 50 #include "Epetra_Export.h" 51 #include "Epetra_IntVector.h" 53 #ifdef HAVE_IFPACK_EPETRAEXT 54 #include "EpetraExt_MatrixMatrix.h" 55 #include "EpetraExt_RowMatrixOut.h" 56 #include "EpetraExt_MultiVectorOut.h" 57 #include "EpetraExt_Transpose_RowMatrix.h" 60 #define ABS(x) ((x)>=0?(x):-(x)) 61 #define MIN(x,y) ((x)<(y)?(x):(y)) 62 #define MAX(x,y) ((x)>(y)?(x):(y)) 65 int* FindLocalDiricheltRowsFromOnesAndZeros(
const Epetra_CrsMatrix & Matrix,
int &numBCRows);
66 void Apply_BCsToMatrixRowsAndColumns(
const int *dirichletRows,
int numBCRows,
const Epetra_IntVector &dirichletColumns,
const Epetra_CrsMatrix & Matrix);
67 Epetra_IntVector * FindLocalDirichletColumnsFromRows(
const int *dirichletRows,
int numBCRows,
const Epetra_CrsMatrix & Matrix);
69 Ifpack_SORa::Ifpack_SORa(Epetra_RowMatrix* A):
70 IsInitialized_(
false),
77 HaveOAZBoundaries_(
false),
78 UseInterprocDamping_(
false),
79 UseGlobalDamping_(
false),
82 PowerMethodIters_(10),
85 Epetra_CrsMatrix *Acrs=
dynamic_cast<Epetra_CrsMatrix*
>(A);
87 Acrs_= Teuchos::rcp (Acrs,
false);
89 A_ = Teuchos::rcp (A,
false);
92 void Ifpack_SORa::Destroy(){
97 int Ifpack_SORa::Initialize(){
98 Alpha_ = List_.get(
"sora: alpha", Alpha_);
99 Gamma_ = List_.get(
"sora: gamma",Gamma_);
100 NumSweeps_ = List_.get(
"sora: sweeps",NumSweeps_);
101 HaveOAZBoundaries_ = List_.get(
"sora: oaz boundaries", HaveOAZBoundaries_);
102 UseInterprocDamping_= List_.get(
"sora: use interproc damping",UseInterprocDamping_);
103 UseGlobalDamping_ = List_.get(
"sora: use global damping",UseGlobalDamping_);
104 LambdaMax_ = List_.get(
"sora: max eigenvalue",LambdaMax_);
105 LambdaMaxBoost_ = List_.get(
"sora: eigen-analysis: boost",LambdaMaxBoost_);
106 PowerMethodIters_ = List_.get(
"sora: eigen-analysis: max iters",PowerMethodIters_);
108 if (A_->Comm().NumProc() != 1) IsParallel_ =
true;
112 UseInterprocDamping_=
false;
121 int Ifpack_SORa::SetParameters(Teuchos::ParameterList& parameterlist){
127 template<
typename int_type>
128 int Ifpack_SORa::TCompute(){
132 Epetra_Map *RowMap=
const_cast<Epetra_Map*
>(&A_->RowMatrixRowMap());
133 Epetra_Vector Adiag(*RowMap);
134 Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
135 int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
136 double *vals_s,*vals_h;
137 bool RowMatrixMode=(Acrs_==Teuchos::null);
140 sprintf(Label_,
"IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_);
143 if(!A_->Comm().MyPID()) cout<<
"SORa: RowMatrix mode enabled"<<endl;
145 Epetra_RowMatrix *Arow=&*A_;
146 Epetra_Map *ColMap=
const_cast<Epetra_Map*
>(&A_->RowMatrixColMap());
148 int Nmax=A_->MaxNumEntries();
150 std::vector<double> values(Nmax);
151 Epetra_CrsMatrix *Acrs=
new Epetra_CrsMatrix(Copy,*RowMap,Nmax);
153 std::vector<int> indices_local(Nmax);
154 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 155 if(RowMap->GlobalIndicesInt()) {
156 for(
int i=0;i<Arow->NumMyRows();i++) {
157 Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
158 for(
int j=0;j<length;j++)
159 indices_local[j]= ColMap->GID(indices_local[j]);
160 Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices_local[0]);
165 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 166 if(RowMap->GlobalIndicesLongLong()) {
167 std::vector<int_type> indices_global(Nmax);
168 for(
int i=0;i<Arow->NumMyRows();i++) {
169 Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
170 for(
int j=0;j<length;j++)
171 indices_global[j]= (int_type) ColMap->GID64(indices_local[j]);
172 Acrs->InsertGlobalValues((int_type) RowMap->GID64(i),length,&values[0],&indices_global[0]);
177 throw "Ifpack_SORa::TCompute: GlobalIndices type unknown";
179 Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
180 Acrs_ = Teuchos::rcp (Acrs,
true);
187 Askew2->FillComplete();
191 Aherm2->FillComplete();
193 int nnz2=Askew2->NumMyNonzeros();
194 int N=Askew2->NumMyRows();
198 IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
199 IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
206 for(
int i=0;i<nnz2;i++)
212 if(HaveOAZBoundaries_){
214 int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
215 Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
216 Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
217 Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
223 A_->ExtractDiagonalCopy(Adiag);
226 Epetra_Vector *Wdiag =
new Epetra_Vector(*RowMap);
230 int maxentries=Askew2->MaxNumEntries();
231 int_type* gids=
new int_type [maxentries];
232 double* newvals=
new double[maxentries];
233 W=
new Epetra_CrsMatrix(Copy,*RowMap,0);
234 for(
int i=0;i<
N;i++){
236 int_type rowgid = (int_type) Acrs_->GRID64(i);
241 for(
int j=rowptr_s[i];j<rowptr_s[i+1];j++){
242 int_type colgid = (int_type) Askew2->GCID64(colind_s[j]);
243 c_data+=fabs(vals_s[j]);
246 if(colind_s[j] <
N) {
248 newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
252 ipdamp+=fabs(vals_h[j]/2 + Alpha_ * vals_s[j]/2);
260 double w_val= c_data*Alpha_*Gamma_/4 + Adiag[Acrs_->LRID(rowgid)];
261 if(UseInterprocDamping_) w_val+=ipdamp;
263 W->InsertGlobalValues(rowgid,1,&w_val,&rowgid);
266 W->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
268 Wdiag_= Teuchos::rcp(Wdiag);
274 if(UseGlobalDamping_ && LambdaMax_ == -1) {
275 if(List_.isParameter(
"sora: eigen-analysis: random seed")) {
277 seed = List_.get(
"sora: eigen-analysis: random seed",seed);
278 PowerMethod(PowerMethodIters_,LambdaMax_,&seed);
281 PowerMethod(PowerMethodIters_,LambdaMax_);
284 LambdaMax_*=LambdaMaxBoost_;
285 if(!A_->Comm().MyPID()) printf(
"SORa: Global damping parameter = %6.4e (lmax=%6.4e)\n",GetOmega(),LambdaMax_);
299 ComputeTime_ += Time_.ElapsedTime();
303 int Ifpack_SORa::Compute(){
304 if(!IsInitialized_) Initialize();
306 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 307 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
308 ret_val = TCompute<int>();
312 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 313 if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
314 ret_val = TCompute<long long>();
318 throw "Ifpack_SORa::Compute: GlobalIndices type unknown";
324 int Ifpack_SORa::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
325 if(!IsComputed_)
return -1;
326 Time_.ResetStartTime();
327 bool initial_guess_is_zero=
false;
328 const int lclNumRows = W_->NumMyRows();
330 Epetra_MultiVector Temp(A_->RowMatrixRowMap(),
NumVectors);
332 double omega=GetOmega();
335 Teuchos::RCP<const Epetra_MultiVector> Xcopy;
336 if (X.Pointers()[0] == Y.Pointers()[0]){
337 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
340 initial_guess_is_zero=
true;
343 Xcopy = Teuchos::rcp( &X,
false );
345 Teuchos::RCP< Epetra_MultiVector > T2;
351 if (IsParallel_ && W_->Importer()) T2 = Teuchos::rcp(
new Epetra_MultiVector(W_->Importer()->TargetMap(),
NumVectors,
true));
352 else T2 = Teuchos::rcp(
new Epetra_MultiVector(A_->RowMatrixRowMap(),
NumVectors,
true));
357 double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
358 T2->ExtractView(&t2_ptr);
359 Y.ExtractView(&y_ptr);
360 Temp.ExtractView(&t_ptr);
361 Xcopy->ExtractView(&x_ptr);
362 Wdiag_->ExtractView(&d_ptr);
366 for(
int i=0; i<NumSweeps_; i++){
368 if(!initial_guess_is_zero || i > 0) {
370 Temp.Update(1.0,*Xcopy,-1.0);
373 Temp.Update(1.0,*Xcopy,0.0);
380 for(
int j=0; j<lclNumRows; j++){
381 double diag=d_ptr[j];
386 for(
int k=rowptr[j];k<rowptr[j+1];k++){
387 dtmp+= values[k]*t2_ptr[m][colind[k]];
390 t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;
391 y_ptr[m][j] += omega*t2_ptr[m][j];
398 ApplyInverseTime_ += Time_.ElapsedTime();
403 std::ostream& Ifpack_SORa::Print(std::ostream& os)
const{
406 os<<
"Ifpack_SORa"<<endl;
407 os<<
" alpha = "<<Alpha_<<endl;
408 os<<
" gamma = "<<Gamma_<<endl;
417 Epetra_RowMatrix* Matrix_in){
426 inline int* FindLocalDiricheltRowsFromOnesAndZeros(
const Epetra_CrsMatrix & Matrix,
int &numBCRows){
427 int *dirichletRows =
new int[Matrix.NumMyRows()];
429 for (
int i=0; i<Matrix.NumMyRows(); i++) {
430 int numEntries, *cols;
432 int ierr = Matrix.ExtractMyRowView(i,numEntries,vals,cols);
435 for (
int j=0; j<numEntries; j++)
if (vals[j]!=0) nz++;
436 if (nz==1) dirichletRows[numBCRows++] = i;
438 if(nz==2) dirichletRows[numBCRows++] = i;
441 return dirichletRows;
447 inline Epetra_IntVector * FindLocalDirichletColumnsFromRows(
const int *dirichletRows,
int numBCRows,
const Epetra_CrsMatrix & Matrix){
451 Epetra_IntVector ZeroRows(Matrix.RowMap());
452 Epetra_IntVector *ZeroCols=
new Epetra_IntVector(Matrix.ColMap());
453 ZeroRows.PutValue(0);
454 ZeroCols->PutValue(0);
457 if(Matrix.RowMap().SameAs(Matrix.ColMap())){
458 for (
int i=0; i < numBCRows; i++)
459 (*ZeroCols)[dirichletRows[i]]=1;
464 for (
int i=0; i < numBCRows; i++)
465 ZeroRows[dirichletRows[i]]=1;
468 if(Matrix.RowMap().SameAs(Matrix.DomainMap())){
470 ZeroCols->Import(ZeroRows,*Matrix.Importer(),Insert);
474 Epetra_Import Importer(Matrix.ColMap(),Matrix.RowMap());
475 ZeroCols->Import(ZeroRows,Importer,Insert);
482 inline void Apply_BCsToMatrixRowsAndColumns(
const int *dirichletRows,
int numBCRows,
const Epetra_IntVector &dirichletColumns,
const Epetra_CrsMatrix & Matrix){
487 for(
int i=0;i<numBCRows;i++){
488 int numEntries, *cols;
490 Matrix.ExtractMyRowView(dirichletRows[i],numEntries,vals,cols);
491 for (
int j=0; j<numEntries; j++) vals[j]=0.0;
495 for (
int i=0; i < Matrix.NumMyRows(); i++) {
499 Matrix.ExtractMyRowView(i,numEntries,vals,cols);
500 for (
int j=0; j < numEntries; j++) {
501 if (dirichletColumns[ cols[j] ] > 0) vals[j] = 0.0;
511 PowerMethod(
const int MaximumIterations,
double& lambda_max,
const unsigned int * RngSeed)
515 double RQ_top, RQ_bottom, norm;
516 Epetra_Vector x(A_->OperatorDomainMap());
517 Epetra_Vector y(A_->OperatorRangeMap());
518 Epetra_Vector z(A_->OperatorRangeMap());
520 if(RngSeed) x.SetSeed(*RngSeed);
528 int NumSweepsBackup=NumSweeps_;
529 bool UseGlobalDampingBackup=UseGlobalDamping_;
530 NumSweeps_=1;UseGlobalDamping_=
false;
532 for (
int iter = 0; iter < MaximumIterations; ++iter)
537 y.Update(1.0,x,-1.0);
541 x.Dot(x, &RQ_bottom);
542 lambda_max = RQ_top / RQ_bottom;
545 x.Update(1.0 / norm, y, 0.0);
553 NumSweeps_=NumSweepsBackup;
554 UseGlobalDamping_=UseGlobalDampingBackup;
Ifpack_CondestType
Ifpack_CondestType: enum to define the type of condition number estimate.
static int Add(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
Teuchos::RefCountPtr< Epetra_RowMatrix > A_
Pointer to the matrix to be preconditioned.
#define IFPACK_CHK_ERR(ifpack_err)