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);
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),
83 Epetra_CrsMatrix *Acrs=
dynamic_cast<Epetra_CrsMatrix*
>(A);
85 Acrs_= Teuchos::rcp (Acrs,
false);
87 A_ = Teuchos::rcp (A,
false);
90 void Ifpack_SORa::Destroy(){
95 int Ifpack_SORa::Initialize(){
96 Alpha_ = List_.get(
"sora: alpha", Alpha_);
97 Gamma_ = List_.get(
"sora: gamma",Gamma_);
98 NumSweeps_ = List_.get(
"sora: sweeps",NumSweeps_);
99 HaveOAZBoundaries_= List_.get(
"sora: oaz boundaries", HaveOAZBoundaries_);
100 UseInterprocDamping_ = List_.get(
"sora: use interproc damping",UseInterprocDamping_);
101 UseGlobalDamping_ = List_.get(
"sora: use global damping",UseGlobalDamping_);
103 if (A_->Comm().NumProc() != 1) IsParallel_ =
true;
107 UseInterprocDamping_=
false;
116 int Ifpack_SORa::SetParameters(Teuchos::ParameterList& parameterlist){
122 template<
typename int_type>
123 int Ifpack_SORa::TCompute(){
127 Epetra_Map *RowMap=
const_cast<Epetra_Map*
>(&A_->RowMatrixRowMap());
128 Epetra_Vector Adiag(*RowMap);
129 Epetra_CrsMatrix *Askew2=0, *Aherm2=0,*W=0;
130 int *rowptr_s,*colind_s,*rowptr_h,*colind_h;
131 double *vals_s,*vals_h;
132 bool RowMatrixMode=(Acrs_==Teuchos::null);
135 sprintf(Label_,
"IFPACK SORa (alpha=%5.2f gamma=%5.2f)",Alpha_,Gamma_);
138 if(!A_->Comm().MyPID()) cout<<
"SORa: RowMatrix mode enabled"<<endl;
140 Epetra_RowMatrix *Arow=&*A_;
141 Epetra_Map *ColMap=
const_cast<Epetra_Map*
>(&A_->RowMatrixColMap());
143 int Nmax=A_->MaxNumEntries();
145 std::vector<double> values(Nmax);
146 Epetra_CrsMatrix *Acrs=
new Epetra_CrsMatrix(Copy,*RowMap,Nmax);
148 std::vector<int> indices_local(Nmax);
149 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 150 if(RowMap->GlobalIndicesInt()) {
151 for(
int i=0;i<Arow->NumMyRows();i++) {
152 Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
153 for(
int j=0;j<length;j++)
154 indices_local[j]= ColMap->GID(indices_local[j]);
155 Acrs->InsertGlobalValues(RowMap->GID(i),length,&values[0],&indices_local[0]);
160 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 161 if(RowMap->GlobalIndicesLongLong()) {
162 std::vector<int_type> indices_global(Nmax);
163 for(
int i=0;i<Arow->NumMyRows();i++) {
164 Arow->ExtractMyRowCopy(i,Nmax,length,&values[0],&indices_local[0]);
165 for(
int j=0;j<length;j++)
166 indices_global[j]= (int_type) ColMap->GID64(indices_local[j]);
167 Acrs->InsertGlobalValues((int_type) RowMap->GID64(i),length,&values[0],&indices_global[0]);
172 throw "Ifpack_SORa::TCompute: GlobalIndices type unknown";
174 Acrs->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
175 Acrs_ = Teuchos::rcp (Acrs,
true);
181 IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,
false,1,*Acrs_,
true,-1,Askew2));
182 Askew2->FillComplete();
185 IFPACK_CHK_ERR(EpetraExt::MatrixMatrix::Add(*Acrs_,
false,1,*Acrs_,
true,1,Aherm2));
186 Aherm2->FillComplete();
188 int nnz2=Askew2->NumMyNonzeros();
189 int N=Askew2->NumMyRows();
193 IFPACK_CHK_ERR(Askew2->ExtractCrsDataPointers(rowptr_s,colind_s,vals_s));
194 IFPACK_CHK_ERR(Aherm2->ExtractCrsDataPointers(rowptr_h,colind_h,vals_h));
200 if(rowptr_s[i]!=rowptr_h[i]) IFPACK_CHK_ERR(-2);
201 for(
int i=0;i<nnz2;i++)
202 if(colind_s[i]!=colind_h[i]) IFPACK_CHK_ERR(-3);
207 if(HaveOAZBoundaries_){
209 int* dirRows=FindLocalDiricheltRowsFromOnesAndZeros(*Acrs_,numBCRows);
210 Epetra_IntVector* dirCols=FindLocalDirichletColumnsFromRows(dirRows,numBCRows,*Aherm2);
211 Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Aherm2);
212 Apply_BCsToMatrixRowsAndColumns(dirRows,numBCRows,*dirCols,*Askew2);
218 A_->ExtractDiagonalCopy(Adiag);
221 Epetra_Vector *Wdiag =
new Epetra_Vector(*RowMap);
225 int maxentries=Askew2->MaxNumEntries();
226 int_type* gids=
new int_type [maxentries];
227 double* newvals=
new double[maxentries];
228 W=
new Epetra_CrsMatrix(Copy,*RowMap,0);
229 for(
int i=0;i<N;i++){
231 int_type rowgid = (int_type) Acrs_->GRID64(i);
236 for(
int j=rowptr_s[i];j<rowptr_s[i+1];j++){
237 int_type colgid = (int_type) Askew2->GCID64(colind_s[j]);
238 c_data+=fabs(vals_s[j]);
241 if(colind_s[j] < N) {
243 newvals[idx]=vals_h[j]/2 + Alpha_ * vals_s[j]/2;
247 ipdamp+=fabs(vals_h[j]/2 + Alpha_ * vals_s[j]/2);
252 IFPACK_CHK_ERR(W->InsertGlobalValues(rowgid,idx,newvals,gids));
255 double w_val= c_data*Alpha_*Gamma_/4 + Adiag[Acrs_->LRID(rowgid)];
256 if(UseInterprocDamping_) w_val+=ipdamp;
258 W->InsertGlobalValues(rowgid,1,&w_val,&rowgid);
259 IFPACK_CHK_ERR(Wdiag->ReplaceGlobalValues(1,&w_val,&rowgid));
261 W->FillComplete(A_->OperatorDomainMap(),A_->OperatorRangeMap());
263 Wdiag_= Teuchos::rcp(Wdiag);
269 if(UseGlobalDamping_) {
270 PowerMethod(10,LambdaMax_);
271 if(!A_->Comm().MyPID()) printf(
"SORa: Global damping parameter = %6.4e (lmax=%6.4e)\n",GetOmega(),LambdaMax_);
285 ComputeTime_ += Time_.ElapsedTime();
289 int Ifpack_SORa::Compute(){
290 if(!IsInitialized_) Initialize();
292 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 293 if(A_->RowMatrixRowMap().GlobalIndicesInt()) {
294 ret_val = TCompute<int>();
298 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 299 if(A_->RowMatrixRowMap().GlobalIndicesLongLong()) {
300 ret_val = TCompute<long long>();
304 throw "Ifpack_SORa::Compute: GlobalIndices type unknown";
310 int Ifpack_SORa::ApplyInverse(
const Epetra_MultiVector& X, Epetra_MultiVector& Y)
const{
311 if(!IsComputed_)
return -1;
312 Time_.ResetStartTime();
313 bool initial_guess_is_zero=
false;
314 const int lclNumRows = W_->NumMyRows();
315 const int NumVectors = X.NumVectors();
316 Epetra_MultiVector Temp(A_->RowMatrixRowMap(),NumVectors);
318 double omega=GetOmega();
321 Teuchos::RCP<const Epetra_MultiVector> Xcopy;
322 if (X.Pointers()[0] == Y.Pointers()[0]){
323 Xcopy = Teuchos::rcp(
new Epetra_MultiVector(X) );
326 initial_guess_is_zero=
true;
329 Xcopy = Teuchos::rcp( &X,
false );
331 Teuchos::RCP< Epetra_MultiVector > T2;
337 if (IsParallel_ && W_->Importer()) T2 = Teuchos::rcp(
new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors,
true));
338 else T2 = Teuchos::rcp(
new Epetra_MultiVector(A_->RowMatrixRowMap(),NumVectors,
true));
343 double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
344 T2->ExtractView(&t2_ptr);
345 Y.ExtractView(&y_ptr);
346 Temp.ExtractView(&t_ptr);
347 Xcopy->ExtractView(&x_ptr);
348 Wdiag_->ExtractView(&d_ptr);
349 IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));
352 for(
int i=0; i<NumSweeps_; i++){
354 if(!initial_guess_is_zero || i > 0) {
356 Temp.Update(1.0,*Xcopy,-1.0);
359 Temp.Update(1.0,*Xcopy,0.0);
366 for(
int j=0; j<lclNumRows; j++){
367 double diag=d_ptr[j];
368 for (
int m=0 ; m<NumVectors; m++) {
372 for(
int k=rowptr[j];k<rowptr[j+1];k++){
373 dtmp+= values[k]*t2_ptr[m][colind[k]];
376 t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;
377 y_ptr[m][j] += omega*t2_ptr[m][j];
384 ApplyInverseTime_ += Time_.ElapsedTime();
389 std::ostream& Ifpack_SORa::Print(std::ostream& os)
const{
392 os<<
"Ifpack_SORa"<<endl;
393 os<<
" alpha = "<<Alpha_<<endl;
394 os<<
" gamma = "<<Gamma_<<endl;
400 double Ifpack_SORa::Condest(
const Ifpack_CondestType CT,
403 Epetra_RowMatrix* Matrix_in){
412 inline int* FindLocalDiricheltRowsFromOnesAndZeros(
const Epetra_CrsMatrix & Matrix,
int &numBCRows){
413 int *dirichletRows =
new int[Matrix.NumMyRows()];
415 for (
int i=0; i<Matrix.NumMyRows(); i++) {
416 int numEntries, *cols;
418 int ierr = Matrix.ExtractMyRowView(i,numEntries,vals,cols);
421 for (
int j=0; j<numEntries; j++)
if (vals[j]!=0) nz++;
422 if (nz==1) dirichletRows[numBCRows++] = i;
424 if(nz==2) dirichletRows[numBCRows++] = i;
427 return dirichletRows;
433 inline Epetra_IntVector * FindLocalDirichletColumnsFromRows(
const int *dirichletRows,
int numBCRows,
const Epetra_CrsMatrix & Matrix){
437 Epetra_IntVector ZeroRows(Matrix.RowMap());
438 Epetra_IntVector *ZeroCols=
new Epetra_IntVector(Matrix.ColMap());
439 ZeroRows.PutValue(0);
440 ZeroCols->PutValue(0);
443 if(Matrix.RowMap().SameAs(Matrix.ColMap())){
444 for (
int i=0; i < numBCRows; i++)
445 (*ZeroCols)[dirichletRows[i]]=1;
450 for (
int i=0; i < numBCRows; i++)
451 ZeroRows[dirichletRows[i]]=1;
454 if(Matrix.RowMap().SameAs(Matrix.DomainMap())){
456 ZeroCols->Import(ZeroRows,*Matrix.Importer(),Insert);
460 Epetra_Import Importer(Matrix.ColMap(),Matrix.RowMap());
461 ZeroCols->Import(ZeroRows,Importer,Insert);
468 inline void Apply_BCsToMatrixRowsAndColumns(
const int *dirichletRows,
int numBCRows,
const Epetra_IntVector &dirichletColumns,
const Epetra_CrsMatrix & Matrix){
473 for(
int i=0;i<numBCRows;i++){
474 int numEntries, *cols;
476 Matrix.ExtractMyRowView(dirichletRows[i],numEntries,vals,cols);
477 for (
int j=0; j<numEntries; j++) vals[j]=0.0;
481 for (
int i=0; i < Matrix.NumMyRows(); i++) {
485 Matrix.ExtractMyRowView(i,numEntries,vals,cols);
486 for (
int j=0; j < numEntries; j++) {
487 if (dirichletColumns[ cols[j] ] > 0) vals[j] = 0.0;
497 PowerMethod(
const int MaximumIterations,
double& lambda_max)
501 double RQ_top, RQ_bottom, norm;
502 Epetra_Vector x(A_->OperatorDomainMap());
503 Epetra_Vector y(A_->OperatorRangeMap());
504 Epetra_Vector z(A_->OperatorRangeMap());
507 if (norm == 0.0) IFPACK_CHK_ERR(-1);
512 int NumSweepsBackup=NumSweeps_;
513 bool UseGlobalDampingBackup=UseGlobalDamping_;
514 NumSweeps_=1;UseGlobalDamping_=
false;
516 for (
int iter = 0; iter < MaximumIterations; ++iter)
521 y.Update(1.0,x,-1.0);
525 x.Dot(x, &RQ_bottom);
526 lambda_max = RQ_top / RQ_bottom;
528 if (norm == 0.0) IFPACK_CHK_ERR(-1);
529 x.Update(1.0 / norm, y, 0.0);
537 NumSweeps_=NumSweepsBackup;
538 UseGlobalDamping_=UseGlobalDampingBackup;