IFPACK  Development
 All Classes 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_LinearPartitioner.h"
51 #include "Ifpack_GreedyPartitioner.h"
52 #include "Ifpack_METISPartitioner.h"
53 #include "Ifpack_EquationPartitioner.h"
54 #include "Ifpack_UserPartitioner.h"
55 #include "Ifpack_Graph_Epetra_RowMatrix.h"
56 #include "Ifpack_DenseContainer.h"
57 #include "Ifpack_Utils.h"
58 #include "Teuchos_ParameterList.hpp"
59 #include "Teuchos_RefCountPtr.hpp"
60 #include "Epetra_RowMatrix.h"
61 #include "Epetra_MultiVector.h"
62 #include "Epetra_Vector.h"
63 #include "Epetra_Time.h"
64 #include "Epetra_Import.h"
65 
66 static const int IFPACK_JACOBI = 0;
67 static const int IFPACK_GS = 1;
68 static const int IFPACK_SGS = 2;
69 
70 
72 
135 template<typename T>
137 
138 public:
139 
141 
147  Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix);
148 
149  virtual ~Ifpack_BlockRelaxation();
150 
152 
154 
156 
164  virtual int Apply(const Epetra_MultiVector& X,
165  Epetra_MultiVector& Y) const;
166 
168 
177  virtual int ApplyInverse(const Epetra_MultiVector& X,
178  Epetra_MultiVector& Y) const;
179 
181  virtual double NormInf() const
182  {
183  return(-1.0);
184  }
186 
188 
189  virtual int SetUseTranspose(bool UseTranspose_in)
190  {
191  if (UseTranspose_in)
192  IFPACK_CHK_ERR(-98); // FIXME: can I work with the transpose?
193  return(0);
194  }
195 
196  virtual const char* Label() const;
197 
199  virtual bool UseTranspose() const
200  {
201  return(false);
202  }
203 
205  virtual bool HasNormInf() const
206  {
207  return(false);
208  }
209 
211  virtual const Epetra_Comm & Comm() const;
212 
214  virtual const Epetra_Map & OperatorDomainMap() const;
215 
217  virtual const Epetra_Map & OperatorRangeMap() const;
219 
221  int NumLocalBlocks() const
222  {
223  return(NumLocalBlocks_);
224  }
225 
227  virtual bool IsInitialized() const
228  {
229  return(IsInitialized_);
230  }
231 
233  virtual bool IsComputed() const
234  {
235  return(IsComputed_);
236  }
237 
239  virtual int SetParameters(Teuchos::ParameterList& List);
240 
242  virtual int Initialize();
243 
245  virtual int Compute();
246 
247  virtual const Epetra_RowMatrix& Matrix() const
248  {
249  return(*Matrix_);
250  }
251 
252  virtual double Condest(const Ifpack_CondestType CT = Ifpack_Cheap,
253  const int MaxIters = 1550,
254  const double Tol = 1e-9,
255  Epetra_RowMatrix* Matrix_in = 0)
256  {
257  return(-1.0);
258  }
259 
260  virtual double Condest() const
261  {
262  return(-1.0);
263  }
264 
265  std::ostream& Print(std::ostream& os) const;
266 
268  virtual int NumInitialize() const
269  {
270  return(NumInitialize_);
271  }
272 
274  virtual int NumCompute() const
275  {
276  return(NumCompute_);
277  }
278 
280  virtual int NumApplyInverse() const
281  {
282  return(NumApplyInverse_);
283  }
284 
286  virtual double InitializeTime() const
287  {
288  return(InitializeTime_);
289  }
290 
292  virtual double ComputeTime() const
293  {
294  return(ComputeTime_);
295  }
296 
298  virtual double ApplyInverseTime() const
299  {
300  return(ApplyInverseTime_);
301  }
302 
304  virtual double InitializeFlops() const
305  {
306 #ifdef IFPACK_FLOPCOUNTERS
307  if (Containers_.size() == 0)
308  return(0.0);
309 
310  // the total number of flops is computed each time InitializeFlops() is
311  // called. This is becase I also have to add the contribution from each
312  // container.
313  double total = InitializeFlops_;
314  for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
315  total += Containers_[i]->InitializeFlops();
316  return(total);
317 #else
318  return(0.0);
319 #endif
320  }
321 
322  virtual double ComputeFlops() const
323  {
324 #ifdef IFPACK_FLOPCOUNTERS
325  if (Containers_.size() == 0)
326  return(0.0);
327 
328  double total = ComputeFlops_;
329  for (unsigned int i = 0 ; i < Containers_.size() ; ++i)
330  total += Containers_[i]->ComputeFlops();
331  return(total);
332 #else
333  return(0.0);
334 #endif
335  }
336 
337  virtual double ApplyInverseFlops() const
338  {
339 #ifdef IFPACK_FLOPCOUNTERS
340  if (Containers_.size() == 0)
341  return(0.0);
342 
343  double total = ApplyInverseFlops_;
344  for (unsigned int i = 0 ; i < Containers_.size() ; ++i) {
345  total += Containers_[i]->ApplyInverseFlops();
346  }
347  return(total);
348 #else
349  return(0.0);
350 #endif
351  }
352 
353 private:
354 
357 
359  Ifpack_BlockRelaxation & operator=(const Ifpack_BlockRelaxation& rhs)
360  {
361  return(*this);
362  }
363 
364  virtual int ApplyInverseJacobi(const Epetra_MultiVector& X,
365  Epetra_MultiVector& Y) const;
366 
367  virtual int DoJacobi(const Epetra_MultiVector& X,
368  Epetra_MultiVector& Y) const;
369 
370  virtual int ApplyInverseGS(const Epetra_MultiVector& X,
371  Epetra_MultiVector& Y) const;
372 
373  virtual int DoGaussSeidel(Epetra_MultiVector& X,
374  Epetra_MultiVector& Y) const;
375 
376  virtual int ApplyInverseSGS(const Epetra_MultiVector& X,
377  Epetra_MultiVector& Y) const;
378 
379  virtual int DoSGS(const Epetra_MultiVector& X,
380  Epetra_MultiVector& Xtmp,
381  Epetra_MultiVector& Y) const;
382 
383  int ExtractSubmatrices();
384 
385  // @{ Initializations, timing and flops
386 
388  bool IsInitialized_;
390  bool IsComputed_;
392  int NumInitialize_;
394  int NumCompute_;
396  mutable int NumApplyInverse_;
398  double InitializeTime_;
400  double ComputeTime_;
402  mutable double ApplyInverseTime_;
404  double InitializeFlops_;
406  double ComputeFlops_;
408  mutable double ApplyInverseFlops_;
409  // @}
410 
411  // @{ Settings
413  int NumSweeps_;
415  double DampingFactor_;
417  int NumLocalBlocks_;
419  Teuchos::ParameterList List_;
420  // @}
421 
422  // @{ Other data
425  Teuchos::RefCountPtr< const Epetra_RowMatrix > Matrix_;
426  mutable std::vector<Teuchos::RefCountPtr<T> > Containers_;
428  Teuchos::RefCountPtr<Ifpack_Partitioner> Partitioner_;
429  string PartitionerType_;
430  int PrecType_;
432  string Label_;
434  bool ZeroStartingSolution_;
435  Teuchos::RefCountPtr<Ifpack_Graph> Graph_;
437  Teuchos::RefCountPtr<Epetra_Vector> W_;
438  // Level of overlap among blocks (for Jacobi only).
439  int OverlapLevel_;
440  mutable Epetra_Time Time_;
441  bool IsParallel_;
442  Teuchos::RefCountPtr<Epetra_Import> Importer_;
443  // @}
444 
445 }; // class Ifpack_BlockRelaxation
446 
447 //==============================================================================
448 template<typename T>
450 Ifpack_BlockRelaxation(const Epetra_RowMatrix* Matrix_in) :
451  IsInitialized_(false),
452  IsComputed_(false),
453  NumInitialize_(0),
454  NumCompute_(0),
455  NumApplyInverse_(0),
456  InitializeTime_(0.0),
457  ComputeTime_(0.0),
458  ApplyInverseTime_(0.0),
459  InitializeFlops_(0.0),
460  ComputeFlops_(0.0),
461  ApplyInverseFlops_(0.0),
462  NumSweeps_(1),
463  DampingFactor_(1.0),
464  NumLocalBlocks_(1),
465  Matrix_(Teuchos::rcp(Matrix_in,false)),
466  PartitionerType_("greedy"),
467  PrecType_(IFPACK_JACOBI),
468  ZeroStartingSolution_(true),
469  OverlapLevel_(0),
470  Time_(Comm()),
471  IsParallel_(false)
472 {
473  if (Matrix_in->Comm().NumProc() != 1)
474  IsParallel_ = true;
475 }
476 
477 //==============================================================================
478 template<typename T>
480 {
481 }
482 
483 //==============================================================================
484 template<typename T>
485 const char* Ifpack_BlockRelaxation<T>::Label() const
486 {
487  return(Label_.c_str());
488 }
489 
490 //==============================================================================
491 template<typename T>
493 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
494 {
495  IFPACK_RETURN(Matrix().Apply(X,Y));
496 }
497 
498 //==============================================================================
499 template<typename T>
500 const Epetra_Comm& Ifpack_BlockRelaxation<T>::
501 Comm() const
502 {
503  return(Matrix().Comm());
504 }
505 
506 //==============================================================================
507 template<typename T>
508 const Epetra_Map& Ifpack_BlockRelaxation<T>::
510 {
511  return(Matrix().OperatorDomainMap());
512 }
513 
514 //==============================================================================
515 template<typename T>
516 const Epetra_Map& Ifpack_BlockRelaxation<T>::
518 {
519  return(Matrix().OperatorRangeMap());
520 }
521 
522 //==============================================================================
523 template<typename T>
525 {
526 
527  if (Partitioner_ == Teuchos::null)
528  IFPACK_CHK_ERR(-3);
529 
530  NumLocalBlocks_ = Partitioner_->NumLocalParts();
531 
532  Containers_.resize(NumLocalBlocks());
533 
534  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
535 
536  int rows = Partitioner_->NumRowsInPart(i);
537  Containers_[i] = Teuchos::rcp( new T(rows) );
538 
539  //Ifpack_DenseContainer* DC = 0;
540  //DC = dynamic_cast<Ifpack_DenseContainer*>(Containers_[i]);
541 
542  if (Containers_[i] == Teuchos::null)
543  IFPACK_CHK_ERR(-5);
544 
545  IFPACK_CHK_ERR(Containers_[i]->SetParameters(List_));
546  IFPACK_CHK_ERR(Containers_[i]->Initialize());
547  // flops in Initialize() will be computed on-the-fly in method InitializeFlops().
548 
549  // set "global" ID of each partitioner row
550  for (int j = 0 ; j < rows ; ++j) {
551  int LRID = (*Partitioner_)(i,j);
552  Containers_[i]->ID(j) = LRID;
553  }
554 
555  IFPACK_CHK_ERR(Containers_[i]->Compute(*Matrix_));
556  // flops in Compute() will be computed on-the-fly in method ComputeFlops().
557 
558  }
559 
560  return(0);
561 }
562 
563 //==============================================================================
564 template<typename T>
566 {
567 
568  if (!IsInitialized())
569  IFPACK_CHK_ERR(Initialize());
570 
571  Time_.ResetStartTime();
572 
573  IsComputed_ = false;
574 
575  if (Matrix().NumGlobalRows64() != Matrix().NumGlobalCols64())
576  IFPACK_CHK_ERR(-2); // only square matrices
577 
578  IFPACK_CHK_ERR(ExtractSubmatrices());
579 
580  if (IsParallel_ && PrecType_ != IFPACK_JACOBI) {
581  // not needed by Jacobi (done by matvec)
582  Importer_ = Teuchos::rcp( new Epetra_Import(Matrix().RowMatrixColMap(),
583  Matrix().RowMatrixRowMap()) );
584 
585  if (Importer_ == Teuchos::null) IFPACK_CHK_ERR(-5);
586  }
587  IsComputed_ = true;
588  ComputeTime_ += Time_.ElapsedTime();
589  ++NumCompute_;
590 
591  return(0);
592 
593 }
594 
595 //==============================================================================
596 template<typename T>
598 ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
599 {
600  if (!IsComputed())
601  IFPACK_CHK_ERR(-3);
602 
603  if (X.NumVectors() != Y.NumVectors())
604  IFPACK_CHK_ERR(-2);
605 
606  Time_.ResetStartTime();
607 
608  // AztecOO gives X and Y pointing to the same memory location,
609  // need to create an auxiliary vector, Xcopy
610  Teuchos::RefCountPtr< const Epetra_MultiVector > Xcopy;
611  if (X.Pointers()[0] == Y.Pointers()[0])
612  Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
613  else
614  Xcopy = Teuchos::rcp( &X, false );
615 
616  switch (PrecType_) {
617  case IFPACK_JACOBI:
618  IFPACK_CHK_ERR(ApplyInverseJacobi(*Xcopy,Y));
619  break;
620  case IFPACK_GS:
621  IFPACK_CHK_ERR(ApplyInverseGS(*Xcopy,Y));
622  break;
623  case IFPACK_SGS:
624  IFPACK_CHK_ERR(ApplyInverseSGS(*Xcopy,Y));
625  break;
626  }
627 
628  ApplyInverseTime_ += Time_.ElapsedTime();
629  ++NumApplyInverse_;
630 
631  return(0);
632 }
633 
634 //==============================================================================
635 // This method in general will not work with AztecOO if used
636 // outside Ifpack_AdditiveSchwarz and OverlapLevel_ != 0
637 //
638 template<typename T>
640 ApplyInverseJacobi(const Epetra_MultiVector& X,
641  Epetra_MultiVector& Y) const
642 {
643 
644  if (ZeroStartingSolution_)
645  Y.PutScalar(0.0);
646 
647  // do not compute the residual in this case
648  if (NumSweeps_ == 1 && ZeroStartingSolution_) {
649  IFPACK_RETURN(DoJacobi(X,Y));
650  }
651 
652  Epetra_MultiVector AX(Y);
653 
654  for (int j = 0; j < NumSweeps_ ; j++) {
655  IFPACK_CHK_ERR(Apply(Y,AX));
656  ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalNonzeros64();
657  IFPACK_CHK_ERR(AX.Update(1.0,X,-1.0));
658  ApplyInverseFlops_ += X.NumVectors() * 2 * Matrix_->NumGlobalRows64();
659  IFPACK_CHK_ERR(DoJacobi(AX,Y));
660  // flops counted in DoJacobi()
661  }
662 
663 
664  return(0);
665 }
666 
667 //==============================================================================
668 template<typename T>
670 DoJacobi(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
671 {
672  int NumVectors = X.NumVectors();
673 
674  if (OverlapLevel_ == 0) {
675 
676  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
677 
678  // may happen that a partition is empty
679  if (Containers_[i]->NumRows() == 0)
680  continue;
681 
682  int LID;
683 
684  // extract RHS from X
685  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
686  LID = Containers_[i]->ID(j);
687  for (int k = 0 ; k < NumVectors ; ++k) {
688  Containers_[i]->RHS(j,k) = X[k][LID];
689  }
690  }
691 
692  // apply the inverse of each block. NOTE: flops occurred
693  // in ApplyInverse() of each block are summed up in method
694  // ApplyInverseFlops().
695  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
696 
697  // copy back into solution vector Y
698  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
699  LID = Containers_[i]->ID(j);
700  for (int k = 0 ; k < NumVectors ; ++k) {
701  Y[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
702  }
703  }
704 
705  }
706  // NOTE: flops for ApplyInverse() of each block are summed up
707  // in method ApplyInverseFlops()
708 #ifdef IFPACK_FLOPCOUNTERS
709  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
710 #endif
711 
712  }
713  else {
714 
715  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
716 
717  // may happen that a partition is empty
718  if (Containers_[i]->NumRows() == 0)
719  continue;
720 
721  int LID;
722 
723  // extract RHS from X
724  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
725  LID = Containers_[i]->ID(j);
726  for (int k = 0 ; k < NumVectors ; ++k) {
727  Containers_[i]->RHS(j,k) = (*W_)[LID] * X[k][LID];
728  }
729  }
730 
731  // apply the inverse of each block
732  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
733 
734  // copy back into solution vector Y
735  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
736  LID = Containers_[i]->ID(j);
737  for (int k = 0 ; k < NumVectors ; ++k) {
738  Y[k][LID] += DampingFactor_ * (*W_)[LID] * Containers_[i]->LHS(j,k);
739  }
740  }
741 
742  }
743  // NOTE: flops for ApplyInverse() of each block are summed up
744  // in method ApplyInverseFlops()
745  // NOTE: do not count for simplicity the flops due to overlapping rows
746 #ifdef IFPACK_FLOPCOUNTERS
747  ApplyInverseFlops_ += NumVectors * 4 * Matrix_->NumGlobalRows();
748 #endif
749  }
750 
751  return(0);
752 }
753 
754 //==============================================================================
755 template<typename T>
757 ApplyInverseGS(const Epetra_MultiVector& X,
758  Epetra_MultiVector& Y) const
759 {
760 
761  if (ZeroStartingSolution_)
762  Y.PutScalar(0.0);
763 
764  Epetra_MultiVector Xcopy(X);
765  for (int j = 0; j < NumSweeps_ ; j++) {
766  IFPACK_CHK_ERR(DoGaussSeidel(Xcopy,Y));
767  if (j != NumSweeps_ - 1)
768  Xcopy = X;
769  }
770 
771  return(0);
772 
773 }
774 
775 //==============================================================================
776 template<typename T>
778 DoGaussSeidel(Epetra_MultiVector& X, Epetra_MultiVector& Y) const
779 {
780 
781  // cycle over all local subdomains
782 
783  int Length = Matrix().MaxNumEntries();
784  std::vector<int> Indices(Length);
785  std::vector<double> Values(Length);
786 
787  int NumMyRows = Matrix().NumMyRows();
788  int NumVectors = X.NumVectors();
789 
790  // an additonal vector is needed by parallel computations
791  // (note that applications through Ifpack_AdditiveSchwarz
792  // are always seen are serial)
793  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
794  if (IsParallel_)
795  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
796  else
797  Y2 = Teuchos::rcp( &Y, false );
798 
799  double** y_ptr;
800  double** y2_ptr;
801  Y.ExtractView(&y_ptr);
802  Y2->ExtractView(&y2_ptr);
803 
804  // data exchange is here, once per sweep
805  if (IsParallel_)
806  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
807 
808  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
809 
810  // may happen that a partition is empty
811  if (Containers_[i]->NumRows() == 0)
812  continue;
813 
814  int LID;
815 
816  // update from previous block
817 
818  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
819  LID = Containers_[i]->ID(j);
820 
821  int NumEntries;
822  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
823  &Values[0], &Indices[0]));
824 
825  for (int k = 0 ; k < NumEntries ; ++k) {
826  int col = Indices[k];
827 
828  for (int kk = 0 ; kk < NumVectors ; ++kk) {
829  X[kk][LID] -= Values[k] * y2_ptr[kk][col];
830  }
831  }
832  }
833 
834  // solve with this block
835 
836  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
837  LID = Containers_[i]->ID(j);
838  for (int k = 0 ; k < NumVectors ; ++k) {
839  Containers_[i]->RHS(j,k) = X[k][LID];
840  }
841  }
842 
843  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
844 #ifdef IFPACK_FLOPCOUNTERS
845  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
846 #endif
847 
848  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
849  LID = Containers_[i]->ID(j);
850  for (int k = 0 ; k < NumVectors ; ++k) {
851  y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
852  }
853  }
854 
855  }
856 
857  // operations for all getrow()'s
858  // NOTE: flops for ApplyInverse() of each block are summed up
859  // in method ApplyInverseFlops()
860 #ifdef IFPACK_FLOPCOUNTERS
861  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
862  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
863 #endif
864 
865  // Attention: this is delicate... Not all combinations
866  // of Y2 and Y will always work (tough for ML it should be ok)
867  if (IsParallel_)
868  for (int m = 0 ; m < NumVectors ; ++m)
869  for (int i = 0 ; i < NumMyRows ; ++i)
870  y_ptr[m][i] = y2_ptr[m][i];
871 
872  return(0);
873 }
874 
875 //==============================================================================
876 template<typename T>
878 ApplyInverseSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
879 {
880 
881  if (ZeroStartingSolution_)
882  Y.PutScalar(0.0);
883 
884  Epetra_MultiVector Xcopy(X);
885  for (int j = 0; j < NumSweeps_ ; j++) {
886  IFPACK_CHK_ERR(DoSGS(X,Xcopy,Y));
887  if (j != NumSweeps_ - 1)
888  Xcopy = X;
889  }
890  return(0);
891 }
892 
893 //==============================================================================
894 template<typename T>
896 DoSGS(const Epetra_MultiVector& X, Epetra_MultiVector& Xcopy,
897  Epetra_MultiVector& Y) const
898 {
899 
900  int NumMyRows = Matrix().NumMyRows();
901  int NumVectors = X.NumVectors();
902 
903  int Length = Matrix().MaxNumEntries();
904  std::vector<int> Indices;
905  std::vector<double> Values;
906  Indices.resize(Length);
907  Values.resize(Length);
908 
909  // an additonal vector is needed by parallel computations
910  // (note that applications through Ifpack_AdditiveSchwarz
911  // are always seen are serial)
912  Teuchos::RefCountPtr< Epetra_MultiVector > Y2;
913  if (IsParallel_)
914  Y2 = Teuchos::rcp( new Epetra_MultiVector(Importer_->TargetMap(), NumVectors) );
915  else
916  Y2 = Teuchos::rcp( &Y, false );
917 
918  double** y_ptr;
919  double** y2_ptr;
920  Y.ExtractView(&y_ptr);
921  Y2->ExtractView(&y2_ptr);
922 
923  // data exchange is here, once per sweep
924  if (IsParallel_)
925  IFPACK_CHK_ERR(Y2->Import(Y,*Importer_,Insert));
926 
927  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
928 
929  // may happen that a partition is empty
930  if (Containers_[i]->NumRows() == 0)
931  continue;
932 
933  int LID;
934 
935  // update from previous block
936 
937  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
938  LID = Containers_[i]->ID(j);
939 
940  int NumEntries;
941  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
942  &Values[0], &Indices[0]));
943 
944  for (int k = 0 ; k < NumEntries ; ++k) {
945  int col = Indices[k];
946 
947  for (int kk = 0 ; kk < NumVectors ; ++kk) {
948  Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
949  }
950  }
951  }
952 
953  // solve with this block
954 
955  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
956  LID = Containers_[i]->ID(j);
957  for (int k = 0 ; k < NumVectors ; ++k) {
958  Containers_[i]->RHS(j,k) = Xcopy[k][LID];
959  }
960  }
961 
962  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
963 #ifdef IFPACK_FLOPCOUNTERS
964  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
965 #endif
966 
967  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
968  LID = Containers_[i]->ID(j);
969  for (int k = 0 ; k < NumVectors ; ++k) {
970  y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
971  }
972  }
973  }
974 
975 #ifdef IFPACK_FLOPCOUNTERS
976  // operations for all getrow()'s
977  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
978  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
979 #endif
980 
981  Xcopy = X;
982 
983  for (int i = NumLocalBlocks() - 1; i >=0 ; --i) {
984 
985  if (Containers_[i]->NumRows() == 0)
986  continue;
987 
988  int LID;
989 
990  // update from previous block
991 
992  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
993  LID = Containers_[i]->ID(j);
994 
995  int NumEntries;
996  IFPACK_CHK_ERR(Matrix().ExtractMyRowCopy(LID, Length,NumEntries,
997  &Values[0], &Indices[0]));
998 
999  for (int k = 0 ; k < NumEntries ; ++k) {
1000  int col = Indices[k];
1001 
1002  for (int kk = 0 ; kk < NumVectors ; ++kk) {
1003  Xcopy[kk][LID] -= Values[k] * y2_ptr[kk][col];
1004  }
1005  }
1006  }
1007 
1008  // solve with this block
1009 
1010  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1011  LID = Containers_[i]->ID(j);
1012  for (int k = 0 ; k < NumVectors ; ++k) {
1013  Containers_[i]->RHS(j,k) = Xcopy[k][LID];
1014  }
1015  }
1016 
1017  IFPACK_CHK_ERR(Containers_[i]->ApplyInverse());
1018 #ifdef IFPACK_FLOPCOUNTERS
1019  ApplyInverseFlops_ += Containers_[i]->ApplyInverseFlops();
1020 #endif
1021 
1022  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1023  LID = Containers_[i]->ID(j);
1024  for (int k = 0 ; k < NumVectors ; ++k) {
1025  y2_ptr[k][LID] += DampingFactor_ * Containers_[i]->LHS(j,k);
1026  }
1027  }
1028  }
1029 
1030 #ifdef IFPACK_FLOPCOUNTERS
1031  // operations for all getrow()'s
1032  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalNonzeros();
1033  ApplyInverseFlops_ += NumVectors * 2 * Matrix_->NumGlobalRows();
1034 #endif
1035 
1036  // Attention: this is delicate... Not all combinations
1037  // of Y2 and Y will always work (tough for ML it should be ok)
1038  if (IsParallel_)
1039  for (int m = 0 ; m < NumVectors ; ++m)
1040  for (int i = 0 ; i < NumMyRows ; ++i)
1041  y_ptr[m][i] = y2_ptr[m][i];
1042 
1043  return(0);
1044 }
1045 
1046 //==============================================================================
1047 template<typename T>
1048 ostream& Ifpack_BlockRelaxation<T>::Print(ostream & os) const
1049 {
1050 
1051  string PT;
1052  if (PrecType_ == IFPACK_JACOBI)
1053  PT = "Jacobi";
1054  else if (PrecType_ == IFPACK_GS)
1055  PT = "Gauss-Seidel";
1056  else if (PrecType_ == IFPACK_SGS)
1057  PT = "symmetric Gauss-Seidel";
1058 
1059  if (!Comm().MyPID()) {
1060  os << endl;
1061  os << "================================================================================" << endl;
1062  os << "Ifpack_BlockRelaxation, " << PT << endl;
1063  os << "Sweeps = " << NumSweeps_ << endl;
1064  os << "Damping factor = " << DampingFactor_;
1065  if (ZeroStartingSolution_)
1066  os << ", using zero starting solution" << endl;
1067  else
1068  os << ", using input starting solution" << endl;
1069  os << "Number of local blocks = " << Partitioner_->NumLocalParts() << endl;
1070  //os << "Condition number estimate = " << Condest_ << endl;
1071  os << "Global number of rows = " << Matrix_->NumGlobalRows64() << endl;
1072  os << endl;
1073  os << "Phase # calls Total Time (s) Total MFlops MFlops/s" << endl;
1074  os << "----- ------- -------------- ------------ --------" << endl;
1075  os << "Initialize() " << std::setw(5) << NumInitialize()
1076  << " " << std::setw(15) << InitializeTime()
1077  << " " << std::setw(15) << 1.0e-6 * InitializeFlops();
1078  if (InitializeTime() != 0.0)
1079  os << " " << std::setw(15) << 1.0e-6 * InitializeFlops() / InitializeTime() << endl;
1080  else
1081  os << " " << std::setw(15) << 0.0 << endl;
1082  os << "Compute() " << std::setw(5) << NumCompute()
1083  << " " << std::setw(15) << ComputeTime()
1084  << " " << std::setw(15) << 1.0e-6 * ComputeFlops();
1085  if (ComputeTime() != 0.0)
1086  os << " " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
1087  else
1088  os << " " << std::setw(15) << 0.0 << endl;
1089  os << "ApplyInverse() " << std::setw(5) << NumApplyInverse()
1090  << " " << std::setw(15) << ApplyInverseTime()
1091  << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
1092  if (ApplyInverseTime() != 0.0)
1093  os << " " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
1094  else
1095  os << " " << std::setw(15) << 0.0 << endl;
1096  os << "================================================================================" << endl;
1097  os << endl;
1098  }
1099 
1100  return(os);
1101 }
1102 
1103 //==============================================================================
1104 template<typename T>
1105 int Ifpack_BlockRelaxation<T>::SetParameters(Teuchos::ParameterList& List)
1106 {
1107 
1108  string PT;
1109  if (PrecType_ == IFPACK_JACOBI)
1110  PT = "Jacobi";
1111  else if (PrecType_ == IFPACK_GS)
1112  PT = "Gauss-Seidel";
1113  else if (PrecType_ == IFPACK_SGS)
1114  PT = "symmetric Gauss-Seidel";
1115 
1116  PT = List.get("relaxation: type", PT);
1117 
1118  if (PT == "Jacobi") {
1119  PrecType_ = IFPACK_JACOBI;
1120  }
1121  else if (PT == "Gauss-Seidel") {
1122  PrecType_ = IFPACK_GS;
1123  }
1124  else if (PT == "symmetric Gauss-Seidel") {
1125  PrecType_ = IFPACK_SGS;
1126  } else {
1127  cerr << "Option `relaxation: type' has an incorrect value ("
1128  << PT << ")'" << endl;
1129  cerr << "(file " << __FILE__ << ", line " << __LINE__ << ")" << endl;
1130  exit(EXIT_FAILURE);
1131  }
1132 
1133  NumSweeps_ = List.get("relaxation: sweeps", NumSweeps_);
1134  DampingFactor_ = List.get("relaxation: damping factor",
1135  DampingFactor_);
1136  ZeroStartingSolution_ = List.get("relaxation: zero starting solution",
1137  ZeroStartingSolution_);
1138  PartitionerType_ = List.get("partitioner: type",
1139  PartitionerType_);
1140  NumLocalBlocks_ = List.get("partitioner: local parts",
1141  NumLocalBlocks_);
1142  // only Jacobi can work with overlap among local domains,
1143  OverlapLevel_ = List.get("partitioner: overlap",
1144  OverlapLevel_);
1145 
1146  // check parameters
1147  if (PrecType_ != IFPACK_JACOBI)
1148  OverlapLevel_ = 0;
1149  if (NumLocalBlocks_ < 0)
1150  NumLocalBlocks_ = Matrix().NumMyRows() / (-NumLocalBlocks_);
1151  // other checks are performed in Partitioner_
1152 
1153  // copy the list as each subblock's constructor will
1154  // require it later
1155  List_ = List;
1156 
1157  // set the label
1158  string PT2;
1159  if (PrecType_ == IFPACK_JACOBI)
1160  PT2 = "BJ";
1161  else if (PrecType_ == IFPACK_GS)
1162  PT2 = "BGS";
1163  else if (PrecType_ == IFPACK_SGS)
1164  PT2 = "BSGS";
1165  Label_ = "IFPACK (" + PT2 + ", sweeps="
1166  + Ifpack_toString(NumSweeps_) + ", damping="
1167  + Ifpack_toString(DampingFactor_) + ", blocks="
1168  + Ifpack_toString(NumLocalBlocks()) + ")";
1169 
1170  return(0);
1171 }
1172 
1173 //==============================================================================
1174 template<typename T>
1176 {
1177  IsInitialized_ = false;
1178  Time_.ResetStartTime();
1179 
1180  Graph_ = Teuchos::rcp( new Ifpack_Graph_Epetra_RowMatrix(Teuchos::rcp(&Matrix(),false)) );
1181  if (Graph_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1182 
1183  if (PartitionerType_ == "linear")
1184  Partitioner_ = Teuchos::rcp( new Ifpack_LinearPartitioner(&*Graph_) );
1185  else if (PartitionerType_ == "greedy")
1186  Partitioner_ = Teuchos::rcp( new Ifpack_GreedyPartitioner(&*Graph_) );
1187  else if (PartitionerType_ == "metis")
1188  Partitioner_ = Teuchos::rcp( new Ifpack_METISPartitioner(&*Graph_) );
1189  else if (PartitionerType_ == "equation")
1190  Partitioner_ = Teuchos::rcp( new Ifpack_EquationPartitioner(&*Graph_) );
1191  else if (PartitionerType_ == "user")
1192  Partitioner_ = Teuchos::rcp( new Ifpack_UserPartitioner(&*Graph_) );
1193  else
1194  IFPACK_CHK_ERR(-2);
1195 
1196  if (Partitioner_ == Teuchos::null) IFPACK_CHK_ERR(-5);
1197 
1198  // need to partition the graph of A
1199  IFPACK_CHK_ERR(Partitioner_->SetParameters(List_));
1200  IFPACK_CHK_ERR(Partitioner_->Compute());
1201 
1202  // get actual number of partitions
1203  NumLocalBlocks_ = Partitioner_->NumLocalParts();
1204 
1205  // weight of each vertex
1206  W_ = Teuchos::rcp( new Epetra_Vector(Matrix().RowMatrixRowMap()) );
1207  W_->PutScalar(0.0);
1208 
1209  for (int i = 0 ; i < NumLocalBlocks() ; ++i) {
1210 
1211  for (int j = 0 ; j < Partitioner_->NumRowsInPart(i) ; ++j) {
1212  int LID = (*Partitioner_)(i,j);
1213  (*W_)[LID]++;
1214  }
1215  }
1216  W_->Reciprocal(*W_);
1217 
1218  InitializeTime_ += Time_.ElapsedTime();
1219  IsInitialized_ = true;
1220  ++NumInitialize_;
1221 
1222  return(0);
1223 }
1224 
1225 //==============================================================================
1226 #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'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.
std::ostream & Print(std::ostream &os) const
Prints basic information on iostream. This function is used by operator<<.
int NumLocalBlocks() const
Returns the number local blocks.
virtual bool IsComputed() const
Returns true if the preconditioner has been successfully computed.
string Ifpack_toString(const int &x)
Convertes an integer to string.
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 Compute()
Computes the preconditioner.
Ifpack_METISPartitioner: A class to decompose Ifpack_Graph'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.
Ifpack_Graph_Epetra_RowMatrix: a class to define Ifpack_Graph as a light-weight conversion of Epetra_...
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 Condest(const Ifpack_CondestType CT=Ifpack_Cheap, const int MaxIters=1550, const double Tol=1e-9, Epetra_RowMatrix *Matrix_in=0)
Computes the condition number estimate, returns its value.
virtual double InitializeTime() const
Returns the time spent in Initialize().
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 NormInf() const
Returns the infinity norm of the global matrix (not implemented)
Ifpack_GreedyPartitioner: A class to decompose Ifpack_Graph's using a simple greedy algorithm...