44 #include "Ifpack_SORa.h"
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);
70 IsInitialized_(false),
77 HaveOAZBoundaries_(false),
78 UseInterprocDamping_(false),
79 UseGlobalDamping_(false),
82 PowerMethodIters_(10),
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(){
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;
148 int Nmax=A_->MaxNumEntries();
150 std::vector<double> values(Nmax);
153 std::vector<int> indices_local(Nmax);
154 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
158 for(
int j=0;j<length;j++)
159 indices_local[j]= ColMap->
GID(indices_local[j]);
165 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
167 std::vector<int_type> indices_global(Nmax);
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);
186 IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,
false,1,*Acrs_,
true,-1,Askew2));
190 IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,
false,1,*Acrs_,
true,1,Aherm2));
191 Aherm2->FillComplete();
199 IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
205 if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
206 for(
int i=0;i<nnz2;i++)
207 if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
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);
231 int_type* gids=
new int_type [maxentries];
232 double* newvals=
new double[maxentries];
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);
257 IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));
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";
325 if(!IsComputed_)
return -1;
326 Time_.ResetStartTime();
327 bool initial_guess_is_zero=
false;
328 const int lclNumRows = W_->NumMyRows();
329 const int NumVectors = X.NumVectors();
332 double omega=GetOmega();
335 Teuchos::RCP<const Epetra_MultiVector> Xcopy;
336 if (X.Pointers()[0] == Y.Pointers()[0]){
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);
363 IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));
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];
382 for (
int m=0 ; m<NumVectors; m++) {
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;
414 double Ifpack_SORa::Condest(
const Ifpack_CondestType ,
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;
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;
454 ZeroCols->PutValue(0);
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;
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;
491 for (
int j=0; j<numEntries; j++) vals[j]=0.0;
495 for (
int i=0; i < Matrix.
NumMyRows(); i++) {
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;
520 if(RngSeed) x.SetSeed(*RngSeed);
523 if (norm == 0.0) IFPACK_CHK_ERR(-1);
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;
544 if (norm == 0.0) IFPACK_CHK_ERR(-1);
545 x.Update(1.0 / norm, y, 0.0);
553 NumSweeps_=NumSweepsBackup;
554 UseGlobalDamping_=UseGlobalDampingBackup;
bool GlobalIndicesLongLong() const
bool SameAs(const Epetra_BlockMap &Map) const
bool GlobalIndicesInt() const
const Epetra_Map & ColMap() const
int FillComplete(bool OptimizeDataStorage=true)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
int ReplaceGlobalValues(int NumEntries, const double *Values, const int *Indices)
const Epetra_Map & RowMap() const
const Epetra_Import * Importer() const
int MaxNumEntries() const
virtual int NumMyRows() const =0
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
const Epetra_Map & DomainMap() const
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
int ExtractCrsDataPointers(int *&IndexOffset, int *&Indices, double *&Values_in) const
int NumMyNonzeros() const