IFPACK  Development
 All Classes Namespaces Files Functions Variables Enumerations Friends Pages
Ifpack_BlockRelaxation.h
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack: Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2002) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #ifndef IFPACK_BLOCKPRECONDITIONER_H
45 #define IFPACK_BLOCKPRECONDITIONER_H
46 
47 #include "Ifpack_ConfigDefs.h"
48 #include "Ifpack_Preconditioner.h"
49 #include "Ifpack_Partitioner.h"
50 #include "Ifpack_LinePartitioner.h"
51 #include "Ifpack_LinearPartitioner.h"
52 #include "Ifpack_GreedyPartitioner.h"
53 #include "Ifpack_METISPartitioner.h"
54 #include "Ifpack_EquationPartitioner.h"
55 #include "Ifpack_UserPartitioner.h"
56 #include "Ifpack_Graph_Epetra_RowMatrix.h"
57 #include "Ifpack_DenseContainer.h"
58 #include "Ifpack_Utils.h"
59 #include "Teuchos_ParameterList.hpp"
60 #include "Teuchos_RefCountPtr.hpp"
61 
62 #include "Epetra_Map.h"
63 #include "Epetra_RowMatrix.h"
64 #include "Epetra_MultiVector.h"
65 #include "Epetra_Vector.h"
66 #include "Epetra_Time.h"
67 #include "Epetra_Import.h"
68 
69 static const int IFPACK_JACOBI = 0;
70 static const int IFPACK_GS = 1;
71 static const int IFPACK_SGS = 2;
72 
73 
75 
138 template<typename T>
140 
141 public:
142 
144 
151 
152  virtual ~Ifpack_BlockRelaxation();
153 
155 
157 
159 
167  virtual int Apply(const Epetra_MultiVector& X,
168  Epetra_MultiVector& Y) const;
169 
171 
180  virtual int ApplyInverse(const Epetra_MultiVector& X,
181  Epetra_MultiVector& Y) const;
182 
184  virtual double NormInf() const
185  {
186  return(-1.0);
187  }
189 
191 
192  virtual int SetUseTranspose(bool UseTranspose_in)
193  {
194  if (UseTranspose_in)
195  IFPACK_CHK_ERR(-98); // FIXME: can I work with the transpose?
196  return(0);
197  }
198 
199  virtual const char* Label() const;
200 
202  virtual bool UseTranspose() const
203  {
204  return(false);
205  }
206 
208  virtual bool HasNormInf() const
209  {
210  return(false);
211  }
212 
214  virtual const Epetra_Comm & Comm() const;
215 
217  virtual const Epetra_Map & OperatorDomainMap() const;
218 
220  virtual const Epetra_Map & OperatorRangeMap() const;
222 
224  int NumLocalBlocks() const
225  {
226  return(NumLocalBlocks_);
227  }
228 
230  virtual bool IsInitialized() const
231  {
232  return(IsInitialized_);
233  }
234 
236  virtual bool IsComputed() const
237  {
238  return(IsComputed_);
239  }
240 
242  virtual int SetParameters(Teuchos::ParameterList& List);
243 
245  virtual int Initialize();
246 
248  virtual int Compute();
249 
250  virtual const Epetra_RowMatrix& Matrix() const
251  {
252  return(*Matrix_);
253  }
254 
255  virtual double Condest(const Ifpack_CondestType /* CT */ = Ifpack_Cheap,
256  const int /* MaxIters */ = 1550,
257  const double /* Tol */ = 1e-9,
258  Epetra_RowMatrix* /* Matrix_in */ = 0)
259  {
260  return(-1.0);
261  }
262 
263  virtual double Condest() const
264  {
265  return(-1.0);
266  }
267 
268  std::ostream& Print(std::ostream& os) const;
269 
271  virtual int NumInitialize() const
272  {
273  return(NumInitialize_);
274  }
275 
277  virtual int NumCompute() const
278  {
279  return(NumCompute_);
280  }
281 
283  virtual int NumApplyInverse() const
284  {
285  return(NumApplyInverse_);
286  }
287 
289  virtual double InitializeTime() const
290  {
291  return(InitializeTime_);
292  }
293 
295  virtual double ComputeTime() const
296  {
297  return(ComputeTime_);
298  }
299 
301  virtual double ApplyInverseTime() const
302  {
303  return(ApplyInverseTime_);
304  }
305 
307  virtual double InitializeFlops() const
308  {
309 #ifdef IFPACK_FLOPCOUNTERS
310  if (Containers_.size() == 0)
311  return(0.0);
312 
313  // the total number of flops is computed each time InitializeFlops() is
314  // called. This is becase I also have to add the contribution from each
315  // container.
316  double total = InitializeFlops_;
317  for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
318  if(Containers_[i]) total += Containers_[i]->InitializeFlops();
319  return(total);
320 #else
321  return(0.0);
322 #endif
323  }
324 
325  virtual double ComputeFlops() const
326  {
327 #ifdef IFPACK_FLOPCOUNTERS
328  if (Containers_.size() == 0)
329  return(0.0);
330 
331  double total = ComputeFlops_;
332  for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
333  if(Containers_[i]) total += Containers_[i]->ComputeFlops();
334  return(total);
335 #else
336  return(0.0);
337 #endif
338  }
339 
340  virtual double ApplyInverseFlops() const
341  {
342 #ifdef IFPACK_FLOPCOUNTERS
343  if (Containers_.size() == 0)
344  return(0.0);
345 
346  double total = ApplyInverseFlops_;
347  for (unsigned int i = 0 ; i < Containers_.size() ; ++i) {
348  if(Containers_[i]) total += Containers_[i]->ApplyInverseFlops();
349  }
350  return(total);
351 #else
352  return(0.0);
353 #endif
354  }
355 
356 private:
357 
360 
362  Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation& rhs)
363  {
364  return(*this);
365  }
366 
367  virtual int ApplyInverseJacobi(const Epetra_MultiVector& X,
368  Epetra_MultiVector& Y) const;
369 
370  virtual int DoJacobi(const Epetra_MultiVector& X,
371  Epetra_MultiVector& Y) const;
372 
373  virtual int ApplyInverseGS(const Epetra_MultiVector& X,
374  Epetra_MultiVector& Y) const;
375 
376  virtual int DoGaussSeidel(Epetra_MultiVector& X,
377  Epetra_MultiVector& Y) const;
378 
379  virtual int ApplyInverseSGS(const Epetra_MultiVector& X,
380  Epetra_MultiVector& Y) const;
381 
382  virtual int DoSGS(const Epetra_MultiVector& X,
383  Epetra_MultiVector& Xtmp,
384  Epetra_MultiVector& Y) const;
385 
386  int ExtractSubmatrices();
387 
388  // @{ Initializations, timing and flops
389 
391  bool IsInitialized_;
393  bool IsComputed_;
395  int NumInitialize_;
397  int NumCompute_;
399  mutable int NumApplyInverse_;
401  double InitializeTime_;
403  double ComputeTime_;
405  mutable double ApplyInverseTime_;
407  double InitializeFlops_;
409  double ComputeFlops_;
411  mutable double ApplyInverseFlops_;
412  // @}
413 
414  // @{ Settings
416  int NumSweeps_;
418  double DampingFactor_;
420  int NumLocalBlocks_;
422  Teuchos::ParameterList List_;
423  // @}
424 
425  // @{ Other data
428  Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_;
429  mutable std::vector<Teuchos::RefCountPtr<T> > Containers_;
430  Epetra_Vector Diagonal_ ;
431 
433  Teuchos::RefCountPtr<Ifpack_Partitioner> Partitioner_;
434  std::string PartitionerType_;
435  int PrecType_;
437  std::string Label_;
439  bool ZeroStartingSolution_;
440  Teuchos::RefCountPtr<Ifpack_Graph> Graph_;
442  Teuchos::RefCountPtr<Epetra_Vector> W_;
443  // Level of overlap among blocks (for Jacobi only).
444  int OverlapLevel_;
445  mutable Epetra_Time Time_;
446  bool IsParallel_;
447  Teuchos::RefCountPtr<Epetra_Import> Importer_;
448  // @}
449 
450 }; // class Ifpack_BlockRelaxation
451 
452 //==============================================================================
453 template<typename T>
456  IsInitialized_(false),
457  IsComputed_(false),
458  NumInitialize_(0),
459  NumCompute_(0),
460  NumApplyInverse_(0),
461  InitializeTime_(0.0),
462  ComputeTime_(0.0),
463  ApplyInverseTime_(0.0),
464  InitializeFlops_(0.0),
465  ComputeFlops_(0.0),
466  ApplyInverseFlops_(0.0),
467  NumSweeps_(1),
468  DampingFactor_(1.0),
469  NumLocalBlocks_(1),
470  Matrix_(Teuchos::rcp(Matrix_in,false)),
471  Diagonal_( Matrix_in->Map()),
472  PartitionerType_("greedy"),
473  PrecType_(IFPACK_JACOBI),
474  ZeroStartingSolution_(true),
475  OverlapLevel_(0),
476  Time_(Comm()),
477  IsParallel_(false)
478 {
479  if (Matrix_in->Comm().NumProc() != 1)
480  IsParallel_ = true;
481 }
482 
483 //==============================================================================
484 template<typename T>
486 {
487 }
488 
489 //==============================================================================
490 template<typename T>
491 const char* Ifpack_BlockRelaxation<T>::Label() const
492 {
493  return(Label_.c_str());
494 }
495 
496 //==============================================================================
497 template<typename T>
500 {
501  int ierr = Matrix().Apply(X,Y);
502  IFPACK_RETURN(ierr);
503 }
504 
505 //==============================================================================
506 template<typename T>
508 Comm() const
509 {
510  return(Matrix().Comm());
511 }
512 
513 //==============================================================================
514 template<typename T>
517 {
518  return(Matrix().OperatorDomainMap());
519 }
520 
521 //==============================================================================
522 template<typename T>
525 {
526  return(Matrix().OperatorRangeMap());
527 }
528 
529 //==============================================================================
530 template<typename T>
532 {
533 
534  if (Partitioner_ == Teuchos::null)
535  IFPACK_CHK_ERR(-3);
536 
537  NumLocalBlocks_ = Partitioner_->NumLocalParts();
538 
539  Containers_.resize(NumLocalBlocks());
540 
541  Diagonal_ = Epetra_Vector(Matrix_->Map());
542  Matrix_->ExtractDiagonalCopy(Diagonal_);
543 
544  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
545 
546  int rows = Partitioner_->NumRowsInPart(i);
547  // if rows == 1, then this is a singleton block, and should not be
548  // created. For now, allow creation, and just force the compute step below.
549 
550  if( rows != 1 ) {
551  Containers_[i] = Teuchos::rcp( new T(rows) );
552 
553  IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
554  IFPACK_CHK_ERR(Containers_[i]->Initialize());
555  // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
556 
557  // set "global" ID of each partitioner row
558  for (int j = 0 ; j < rows ; ++j) {
559  int LRID = (*Partitioner_)(i,j);
560  Containers_[i]->ID(j) = LRID;
561  }
562 
563  IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
564  }
565  // otherwise leave Containers_[i] as null
566  }
567 
568 #ifdef SINGLETON_DEBUG
569  int issing = 0;
570 
571  for (int i = 0 ; i < NumLocalBlocks() ; ++i)
572  issing += (int) ( Partitioner_->NumRowsInPart(i) == 1);
573  printf( " %d of %d containers are singleton \n",issing,NumLocalBlocks());
574 #endif
575 
576  return(0);
577 }
578 
579 //==============================================================================
580 template<typename T>
582 {
583 
584  if (!IsInitialized())
585  IFPACK_CHK_ERR(Initialize());
586 
587  Time_.ResetStartTime();
588 
589  IsComputed_ = false;
590 
591  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
592  IFPACK_CHK_ERR(-2); // only square matrices
593 
594  IFPACK_CHK_ERR(ExtractSubmatrices());
595 
596  if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
597  // not needed by Jacobi (done by matvec)
598  Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
599  Matrix().RowMatrixRowMap()) );
600 
601  if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
602  }
603  IsComputed_ = true;
604  ComputeTime_ += Time_.ElapsedTime();
605  ++NumCompute_;
606 
607  return(0);
608 
609 }
610 
611 //==============================================================================
612 template<typename T>
615 {
616  if (!IsComputed())
617  IFPACK_CHK_ERR(-3);
618 
619  if (X.NumVectors() != Y.NumVectors())
620  IFPACK_CHK_ERR(-2);
621 
622  Time_.ResetStartTime();
623 
624  // AztecOO gives X and Y pointing to the same memory location,
625  // need to create an auxiliary vector, Xcopy
626  Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
627  if (X.Pointers()[0] == Y.Pointers()[0])
628  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
629  else
630  Xcopy = Teuchos::rcp( &X, false );
631 
632  switch (PrecType_) {
633  case IFPACK_JACOBI:
634  IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
635  break;
636  case IFPACK_GS:
637  IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
638  break;
639  case IFPACK_SGS:
640  IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
641  break;
642  }
643 
644  ApplyInverseTime_ += Time_.ElapsedTime();
645  ++NumApplyInverse_;
646 
647  return(0);
648 }
649 
650 //==============================================================================
651 // This method in general will not work with AztecOO if used
652 // outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
653 //
654 template<typename T>
657  Epetra_MultiVector& Y) const
658 {
659 
660  if (ZeroStartingSolution_)
661  Y.PutScalar(0.0);
662 
663  // do not compute the residual in this case
664  if (NumSweeps_ == 1 && ZeroStartingSolution_) {
665  int ierr = DoJacobi(X,Y);
666  IFPACK_RETURN(ierr);
667  }
668 
669  Epetra_MultiVector AX(Y);
670 
671  for (int j = 0; j < NumSweeps_ ; j++) {
672  IFPACK_CHK_ERR(Apply(Y,AX));
673  ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros64();
674  IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
675  ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows64();
676  IFPACK_CHK_ERR(DoJacobi(AX,Y));
677  // flops counted in DoJacobi()
678  }
679 
680 
681  return(0);
682 }
683 
684 //==============================================================================
685 template<typename T>
688 {
689  int NumVectors = X.NumVectors();
690 
691  if (OverlapLevel_ == 0) {
692 
693  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
694 
695  int rows = Partitioner_->NumRowsInPart(i);
696  // may happen that a partition is empty
697  if (rows == 0)
698  continue;
699 
700  if(rows != 1) {
701  int LID;
702 
703  // extract RHS from X
704  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
705  LID = Containers_[i]->ID(j);
706  for (int k = 0 ; k < NumVectors ; ++k) {
707  Containers_[i]->RHS(j,k) = X[k][LID];
708  }
709  }
710 
711  // apply the inverse of each block. NOTE: flops occurred
712  // in ApplyInverse() of each block are summed up in method
713  // ApplyInverseFlops().
714 
715  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
716 
717  // copy back into solution vector Y
718  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
719  LID = Containers_[i]->ID(j);
720  for (int k = 0 ; k < NumVectors ; ++k) {
721  Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
722  }
723  }
724  } //end if(rows !=1)
725  else {
726  // rows == 1, this is a singleton. compute directly.
727  int LRID = (*Partitioner_)(i,0);
728  double b = X[0][LRID];
729  double a = Diagonal_[LRID];
730  Y[0][LRID] += DampingFactor_* b/a;
731  }
732  // NOTE: flops for ApplyInverse() of each block are summed up
733  // in method ApplyInverseFlops()
734 #ifdef IFPACK_FLOPCOUNTERS
735  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
736 #endif
737 
738  }
739  }
740  else { // overlap test
741 
742  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
743 
744  int rows = Partitioner_->NumRowsInPart(i);
745 
746  // may happen that a partition is empty
747  if (rows == 0)
748  continue;
749  if(rows != 1) {
750  int LID;
751 
752  // extract RHS from X
753  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
754  LID = Containers_[i]->ID(j);
755  for (int k = 0 ; k < NumVectors ; ++k) {
756  Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
757  }
758  }
759 
760  // apply the inverse of each block
761  // if(rows != 1)
762  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
763 
764  // copy back into solution vector Y
765  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
766  LID = Containers_[i]->ID(j);
767  for (int k = 0 ; k < NumVectors ; ++k) {
768  Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
769  }
770  }
771  } // end if(rows != 1)
772  else { // rows == 1, this is a singleton. compute directly.
773  int LRID = (*Partitioner_)(i,0);
774  double w = (*W_)[LRID];
775  double b = w * X[0][LRID];
776  double a = Diagonal_[LRID];
777 
778  Y[0][LRID] += DampingFactor_ * w * b / a;
779  }
780  }
781  // NOTE: flops for ApplyInverse() of each block are summed up
782  // in method ApplyInverseFlops()
783  // NOTE: do not count for simplicity the flops due to overlapping rows
784 #ifdef IFPACK_FLOPCOUNTERS
785  ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
786 #endif
787  }
788 
789  return(0);
790 }
791 
792 //==============================================================================
793 template<typename T>
796  Epetra_MultiVector& Y) const
797 {
798 
799  if (ZeroStartingSolution_)
800  Y.PutScalar(0.0);
801 
802  Epetra_MultiVector Xcopy(X);
803  for (int j = 0; j < NumSweeps_ ; j++) {
804  IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
805  if (j != NumSweeps_ - 1)
806  Xcopy = X;
807  }
808 
809  return(0);
810 
811 }
812 
813 //==============================================================================
814 template<typename T>
817 {
818 
819  // cycle over all local subdomains
820 
821  int Length = Matrix().MaxNumEntries();
822  std::vector<int> Indices(Length);
823  std::vector<double> Values(Length);
824 
825  int NumMyRows = Matrix().NumMyRows();
826  int NumVectors = X.NumVectors();
827 
828  // an additonal vector is needed by parallel computations
829  // (note that applications through Ifpack_AdditiveSchwarz
830  // are always seen are serial)
831  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
832  if (IsParallel_)
833  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
834  else
835  Y2 = Teuchos::rcp( &Y, false );
836 
837  double** y_ptr;
838  double** y2_ptr;
839  Y.ExtractView(&y_ptr);
840  Y2->ExtractView(&y2_ptr);
841 
842  // data exchange is here, once per sweep
843  if (IsParallel_)
844  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
845 
846  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
847  int rows = Partitioner_->NumRowsInPart(i);
848 
849  // may happen that a partition is empty, but if rows == 1, the container is null
850  if (rows!=1 && Containers_[i]->NumRows() == 0)
851  continue;
852 
853  int LID;
854 
855  // update from previous block
856 
857  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
858  LID = (*Partitioner_)(i,j);
859 
860  int NumEntries;
861  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
862  &Values[0], &Indices[0]));
863 
864  for (int k = 0 ; k < NumEntries ; ++k) {
865  int col = Indices[k];
866 
867  for (int kk = 0 ; kk < NumVectors ; ++kk) {
868  X[kk][LID] -= Values[k] * y2_ptr[kk][col];
869  }
870  }
871  }
872 
873  if(rows != 1) {
874  // solve with this block
875 
876  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
877  LID = Containers_[i]->ID(j);
878  for (int k = 0 ; k < NumVectors ; ++k) {
879  Containers_[i]->RHS(j,k) = X[k][LID];
880  }
881  }
882 
883  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
884 #ifdef IFPACK_FLOPCOUNTERS
885  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
886 #endif
887 
888  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
889  LID = Containers_[i]->ID(j);
890  for (int k = 0 ; k < NumVectors ; ++k) {
891  double temp = DampingFactor_ * Containers_[i]->LHS(j,k);
892  y2_ptr[k][LID] += temp;
893  }
894  }
895  } // end if(rows != 1)
896  else {
897  int LRID = (*Partitioner_)(i,0);
898  double b = X[0][LRID];
899  double a = Diagonal_[LRID];
900  y2_ptr[0][LRID]+= DampingFactor_* b/a;
901  }
902  }
903  // operations for all getrow()'s
904  // NOTE: flops for ApplyInverse() of each block are summed up
905  // in method ApplyInverseFlops()
906 #ifdef IFPACK_FLOPCOUNTERS
907  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
908  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
909 #endif
910 
911  // Attention: this is delicate... Not all combinations
912  // of Y2 and Y will always work (tough for ML it should be ok)
913  if (IsParallel_)
914  for (int m = 0 ; m < NumVectors ; ++m)
915  for (int i = 0 ; i < NumMyRows ; ++i)
916  y_ptr[m][i] = y2_ptr[m][i];
917 
918  return(0);
919  }
920 
921 //==============================================================================
922 template<typename T>
925 {
926 
927  if (ZeroStartingSolution_)
928  Y.PutScalar(0.0);
929 
930  Epetra_MultiVector Xcopy(X);
931  for (int j = 0; j < NumSweeps_ ; j++) {
932  IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
933  if (j != NumSweeps_ - 1)
934  Xcopy = X;
935  }
936  return(0);
937 }
938 
939 //==============================================================================
940 template<typename T>
943  Epetra_MultiVector& Y) const
944 {
945 
946  int NumMyRows = Matrix().NumMyRows();
947  int NumVectors = X.NumVectors();
948 
949  int Length = Matrix().MaxNumEntries();
950  std::vector<int> Indices;
951  std::vector<double> Values;
952  Indices.resize(Length);
953  Values.resize(Length);
954 
955  // an additonal vector is needed by parallel computations
956  // (note that applications through Ifpack_AdditiveSchwarz
957  // are always seen are serial)
958  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
959  if (IsParallel_)
960  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
961  else
962  Y2 = Teuchos::rcp( &Y, false );
963 
964  double** y_ptr;
965  double** y2_ptr;
966  Y.ExtractView(&y_ptr);
967  Y2->ExtractView(&y2_ptr);
968 
969  // data exchange is here, once per sweep
970  if (IsParallel_)
971  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
972 
973  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
974  int rows = Partitioner_->NumRowsInPart(i);
975  // may happen that a partition is empty
976  if (rows !=1 && Containers_[i]->NumRows() == 0)
977  continue;
978 
979  int LID;
980 
981  // update from previous block
982 
983  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
984  LID = (*Partitioner_)(i,j);
985  int NumEntries;
986  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
987  &Values[0], &Indices[0]));
988 
989  for (int k = 0 ; k < NumEntries ; ++k) {
990  int col = Indices[k];
991 
992  for (int kk = 0 ; kk < NumVectors ; ++kk) {
993  Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
994  }
995  }
996  }
997 
998  // solve with this block
999 
1000  if(rows != 1) {
1001  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1002  LID = Containers_[i]->ID(j);
1003  for (int k = 0 ; k < NumVectors ; ++k) {
1004  Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1005  }
1006  }
1007 
1008  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
1009 #ifdef IFPACK_FLOPCOUNTERS
1010  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1011 #endif
1012 
1013  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1014  LID = Containers_[i]->ID(j);
1015  for (int k = 0 ; k < NumVectors ; ++k) {
1016  y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1017  }
1018  }
1019  }
1020  else {
1021  int LRID = (*Partitioner_)(i,0);
1022  double b = Xcopy[0][LRID];
1023  double a = Diagonal_[LRID];
1024  y2_ptr[0][LRID]+= DampingFactor_* b/a;
1025  }
1026  }
1027 
1028 #ifdef IFPACK_FLOPCOUNTERS
1029  // operations for all getrow()'s
1030  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1031  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
1032 #endif
1033 
1034  Xcopy = X;
1035 
1036  for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
1037  int rows = Partitioner_->NumRowsInPart(i);
1038  if (rows != 1 &&Containers_[i]->NumRows() == 0)
1039  continue;
1040 
1041  int LID;
1042 
1043  // update from previous block
1044 
1045  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1046  LID = (*Partitioner_)(i,j);
1047 
1048  int NumEntries;
1049  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
1050  &Values[0], &Indices[0]));
1051 
1052  for (int k = 0 ; k < NumEntries ; ++k) {
1053  int col = Indices[k];
1054 
1055  for (int kk = 0 ; kk < NumVectors ; ++kk) {
1056  Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
1057  }
1058  }
1059  }
1060 
1061  // solve with this block
1062  if(rows != 1) {
1063  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1064  LID = Containers_[i]->ID(j);
1065  for (int k = 0 ; k < NumVectors ; ++k) {
1066  Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1067  }
1068  }
1069 
1070  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
1071 #ifdef IFPACK_FLOPCOUNTERS
1072  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1073 #endif
1074 
1075  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1076  LID = Containers_[i]->ID(j);
1077  for (int k = 0 ; k < NumVectors ; ++k) {
1078  y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1079  }
1080  }
1081  }
1082  else {
1083  int LRID = (*Partitioner_)(i,0);
1084  double b = Xcopy[0][LRID];
1085  double a = Diagonal_[LRID];
1086  y2_ptr[0][LRID]+= DampingFactor_* b/a;
1087  }
1088  }
1089 
1090 #ifdef IFPACK_FLOPCOUNTERS
1091  // operations for all getrow()'s
1092  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1093  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
1094 #endif
1095 
1096  // Attention: this is delicate... Not all combinations
1097  // of Y2 and Y will always work (tough for ML it should be ok)
1098  if (IsParallel_)
1099  for (int m = 0 ; m < NumVectors ; ++m)
1100  for (int i = 0 ; i < NumMyRows ; ++i)
1101  y_ptr[m][i] = y2_ptr[m][i];
1102 
1103  return(0);
1104 }
1105 
1106 //==============================================================================
1107 template<typename T>
1108 std::ostream& Ifpack_BlockRelaxation<T>::Print(std::ostream & os) const
1109 {
1110  using std::endl;
1111 
1112  std::string PT;
1113  if (PrecType_ == IFPACK_JACOBI)
1114  PT = "Jacobi";
1115  else if (PrecType_ == IFPACK_GS)
1116  PT = "Gauss-Seidel";
1117  else if (PrecType_ == IFPACK_SGS)
1118  PT = "symmetric Gauss-Seidel";
1119 
1120  if (!Comm().MyPID()) {
1121  os << endl;
1122  os << "================================================================================" << endl;
1123  os << "Ifpack_BlockRelaxation, " << PT << endl;
1124  os << "Sweeps = " << NumSweeps_ << endl;
1125  os << "Damping factor = " << DampingFactor_;
1126  if (ZeroStartingSolution_)
1127  os << ", using zero starting solution" << endl;
1128  else
1129  os << ", using input starting solution" << endl;
1130  os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
1131  //os << "Condition number estimate = " << Condest_ << endl;
1132  os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1133  os << endl;
1134  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1135  os << "----- ------- -------------- ------------ --------" << endl;
1136  os << "Initialize() " << std::setw(5) << NumInitialize()
1137  << " " << std::setw(15) << InitializeTime()
1138  << " " << std::setw(15) << 1.0e-6 * InitializeFlops();
1139  if (InitializeTime() != 0.0)
1140  os << " " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
1141  else
1142  os << " " << std::setw(15) << 0.0 << endl;
1143  os << "Compute() " << std::setw(5) << NumCompute()
1144  << " " << std::setw(15) << ComputeTime()
1145  << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
1146  if (ComputeTime() != 0.0)
1147  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
1148  else
1149  os << " " << std::setw(15) << 0.0 << endl;
1150  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
1151  << " " << std::setw(15) << ApplyInverseTime()
1152  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
1153  if (ApplyInverseTime() != 0.0)
1154  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
1155  else
1156  os << " " << std::setw(15) << 0.0 << endl;
1157  os << "================================================================================" << endl;
1158  os << endl;
1159  }
1160 
1161  return(os);
1162 }
1163 
1164 //==============================================================================
1165 template<typename T>
1166 int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
1167 {
1168  using std::cerr;
1169  using std::endl;
1170 
1171  std::string PT;
1172  if (PrecType_ == IFPACK_JACOBI)
1173  PT = "Jacobi";
1174  else if (PrecType_ == IFPACK_GS)
1175  PT = "Gauss-Seidel";
1176  else if (PrecType_ == IFPACK_SGS)
1177  PT = "symmetric Gauss-Seidel";
1178 
1179  PT = List.get("relaxation: type", PT);
1180 
1181  if (PT == "Jacobi") {
1182  PrecType_ = IFPACK_JACOBI;
1183  }
1184  else if (PT == "Gauss-Seidel") {
1185  PrecType_ = IFPACK_GS;
1186  }
1187  else if (PT == "symmetric Gauss-Seidel") {
1188  PrecType_ = IFPACK_SGS;
1189  } else {
1190  cerr << "Option `relaxation: type' has an incorrect value ("
1191  << PT << ")'" << endl;
1192  cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
1193  exit(EXIT_FAILURE);
1194  }
1195 
1196  NumSweeps_ = List.get("relaxation: sweeps", NumSweeps_);
1197  DampingFactor_ = List.get("relaxation: damping factor",
1198  DampingFactor_);
1199  ZeroStartingSolution_ = List.get("relaxation: zero starting solution",
1200  ZeroStartingSolution_);
1201  PartitionerType_ = List.get("partitioner: type",
1202  PartitionerType_);
1203  NumLocalBlocks_ = List.get("partitioner: local parts",
1204  NumLocalBlocks_);
1205  // only Jacobi can work with overlap among local domains,
1206  OverlapLevel_ = List.get("partitioner: overlap",
1207  OverlapLevel_);
1208 
1209  // check parameters
1210  if (PrecType_ != IFPACK_JACOBI)
1211  OverlapLevel_ = 0;
1212  if (NumLocalBlocks_ < 0)
1213  NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
1214  // other checks are performed in Partitioner_
1215 
1216  // copy the list as each subblock's constructor will
1217  // require it later
1218  List_ = List;
1219 
1220  // set the label
1221  std::string PT2;
1222  if (PrecType_ == IFPACK_JACOBI)
1223  PT2 = "BJ";
1224  else if (PrecType_ == IFPACK_GS)
1225  PT2 = "BGS";
1226  else if (PrecType_ == IFPACK_SGS)
1227  PT2 = "BSGS";
1228  Label_ = "IFPACK (" + PT2 + ", sweeps="
1229  + Ifpack_toString(NumSweeps_) + ", damping="
1230  + Ifpack_toString(DampingFactor_) + ", blocks="
1231  + Ifpack_toString(NumLocalBlocks()) + ")";
1232 
1233  return(0);
1234 }
1235 
1236 //==============================================================================
1237 template<typename T>
1239 {
1240  IsInitialized_ = false;
1241  Time_.ResetStartTime();
1242 
1243  Graph_ = Teuchos::rcp( new Ifpack_Graph_Epetra_RowMatrix(Teuchos::rcp(&Matrix(),false)) );
1244  if (Graph_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1245 
1246  if (PartitionerType_ == "linear")
1247  Partitioner_ = Teuchos::rcp( new Ifpack_LinearPartitioner(&*Graph_) );
1248  else if (PartitionerType_ == "greedy")
1249  Partitioner_ = Teuchos::rcp( new Ifpack_GreedyPartitioner(&*Graph_) );
1250  else if (PartitionerType_ == "metis")
1251  Partitioner_ = Teuchos::rcp( new Ifpack_METISPartitioner(&*Graph_) );
1252  else if (PartitionerType_ == "equation")
1253  Partitioner_ = Teuchos::rcp( new Ifpack_EquationPartitioner(&*Graph_) );
1254  else if (PartitionerType_ == "user")
1255  Partitioner_ = Teuchos::rcp( new Ifpack_UserPartitioner(&*Graph_) );
1256  else if (PartitionerType_ == "line")
1257  Partitioner_ = Teuchos::rcp( new Ifpack_LinePartitioner(&Matrix()) );
1258  else
1259  IFPACK_CHK_ERR(-2);
1260 
1261  if (Partitioner_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1262 
1263  // need to partition the graph of A
1264  IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
1265  IFPACK_CHK_ERR(Partitioner_->Compute());
1266 
1267  // get actual number of partitions
1268  NumLocalBlocks_ = Partitioner_->NumLocalParts();
1269 
1270  // weight of each vertex
1271  W_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
1272  W_->PutScalar(0.0);
1273 
1274  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
1275 
1276  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1277  int LID = (*Partitioner_)(i,j);
1278  (*W_)[LID]++;
1279  }
1280  }
1281  W_->Reciprocal(*W_);
1282 
1283  // Update Label_ if line smoothing
1284  if (PartitionerType_ == "line") {
1285  // set the label
1286  std::string PT2;
1287  if (PrecType_ == IFPACK_JACOBI)
1288  PT2 = "BJ";
1289  else if (PrecType_ == IFPACK_GS)
1290  PT2 = "BGS";
1291  else if (PrecType_ == IFPACK_SGS)
1292  PT2 = "BSGS";
1293  Label_ = "IFPACK (" + PT2 + ", auto-line, sweeps="
1294  + Ifpack_toString(NumSweeps_) + ", damping="
1295  + Ifpack_toString(DampingFactor_) + ", blocks="
1296  + Ifpack_toString(NumLocalBlocks()) + ")";
1297  }
1298 
1299  InitializeTime_ += Time_.ElapsedTime();
1300  IsInitialized_ = true;
1301  ++NumInitialize_;
1302 
1303  return(0);
1304 }
1305 
1306 //==============================================================================
1307 #endif // IFPACK_BLOCKPRECONDITIONER_H
virtual const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this operator.
Ifpack_BlockRelaxation: a class to define block relaxation preconditioners of Epetra_RowMatrix&#39;s.
Ifpack_UserPartitioner: A class to define linear partitions.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this operator.
virtual double ApplyInverseFlops() const
Returns the number of flops in the application of the preconditioner.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual int NumInitialize() const
Returns the number of calls to Initialize().
virtual double ApplyInverseTime() const
Returns the time spent in ApplyInverse().
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the block Jacobi preconditioner to X, returns the result in Y.
int NumLocalBlocks() const
Returns the number local blocks.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
virtual double ComputeTime() const
Returns the time spent in Compute().
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual bool IsInitialized() const
Returns true if the preconditioner has been successfully computed.
Ifpack_EquationPartitioner: A class to decompose an Ifpack_Graph so that each block will contain all ...
virtual double ComputeFlops() const
Returns the number of flops in the computation phase.
virtual const Epetra_RowMatrix & Matrix() const
Returns a pointer to the matrix to be preconditioned.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Applies the matrix to an Epetra_MultiVector.
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual int MaxNumEntries() const =0
virtual const Epetra_Comm & Comm() const =0
virtual int Compute()
Computes the preconditioner.
Ifpack_METISPartitioner: A class to decompose Ifpack_Graph&#39;s using METIS.
Ifpack_Preconditioner: basic class for preconditioning in Ifpack.
Ifpack_BlockRelaxation(const Epetra_RowMatrix *Matrix)
Ifpack_BlockRelaxation constructor with given Epetra_RowMatrix.
virtual int SetParameters(Teuchos::ParameterList &List)
Sets all the parameters for the preconditioner.
virtual int NumMyRows() const =0
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator&lt;&lt;.
virtual int NumCompute() const
Returns the number of calls to Compute().
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this operator.
virtual double InitializeTime() const
Returns the time spent in Initialize().
virtual int NumProc() const =0
virtual double Condest() const
Returns the computed condition number estimate, or -1.0 if not computed.
virtual int Initialize()
Initializes the preconditioner.
Ifpack_LinearPartitioner: A class to define linear partitions.
virtual double InitializeFlops() const
Returns the number of flops in the initialization phase.
virtual int NumApplyInverse() const
Returns the number of calls to ApplyInverse().
virtual double Condest(const Ifpack_CondestType=Ifpack_Cheap, const int=1550, const double=1e-9, Epetra_RowMatrix *=0)
Computes the condition number estimate, returns its value.
virtual double NormInf() const
Returns the infinity norm of the global matrix (not implemented)
std::string Ifpack_toString(const int &x)
Converts an integer to std::string.
Ifpack_GreedyPartitioner: A class to decompose Ifpack_Graph&#39;s using a simple greedy algorithm...