Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
test/OverlappingRowMatrix_LL/cxx_main.cpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack: Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2002) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #include "Ifpack_ConfigDefs.h"
44 
45 #ifdef HAVE_MPI
46 #include "Epetra_MpiComm.h"
47 #else
48 #include "Epetra_SerialComm.h"
49 #endif
50 #include "Epetra_CrsMatrix.h"
51 #include "Epetra_Vector.h"
52 #include "Epetra_LinearProblem.h"
53 #include "Epetra_Map.h"
54 #include "Epetra_Import.h"
55 #include "Epetra_Time.h"
56 #include "Galeri_Maps.h"
57 #include "Galeri_CrsMatrices.h"
58 #include "Teuchos_ParameterList.hpp"
59 #include "Teuchos_RefCountPtr.hpp"
61 #include "Ifpack_LocalFilter.h"
62 #include "Ifpack_Utils.h"
63 
64 int main(int argc, char *argv[])
65 {
66 #ifdef HAVE_MPI
67  MPI_Init(&argc,&argv);
68  Epetra_MpiComm Comm( MPI_COMM_WORLD );
69 #else
70  Epetra_SerialComm Comm;
71 #endif
72 
73  if (Comm.NumProc() == 1)
74  {
75 #ifdef HAVE_MPI
76  MPI_Finalize();
77 #endif
78  cout << "Test `TestOverlappingRowMatrix.exe' passed!" << endl;
79  exit(EXIT_SUCCESS);
80  }
81 
82  Teuchos::ParameterList GaleriList;
83  int nx = 100;
84  GaleriList.set("n", nx * nx);
85  GaleriList.set("nx", nx);
86  GaleriList.set("ny", nx);
87  Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap64("Linear", Comm, GaleriList) );
88  Teuchos::RefCountPtr<Epetra_CrsMatrix> A = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
89 
90  int OverlapLevel = 5;
91  Epetra_Time Time(Comm);
92 
93  // ======================================== //
94  // Build the overlapping matrix using class //
95  // Ifpack_OverlappingRowMatrix. //
96  // ======================================== //
97 
98  Time.ResetStartTime();
99  Ifpack_OverlappingRowMatrix B(A,OverlapLevel);
100  if (Comm.MyPID() == 0)
101  cout << "Time to create B = " << Time.ElapsedTime() << endl;
102 
103  long long NumGlobalRowsB = B.NumGlobalRows64();
104  long long NumGlobalNonzerosB = B.NumGlobalNonzeros64();
105 
106  Epetra_Vector X(A->RowMatrixRowMap());
107  Epetra_Vector Y(A->RowMatrixRowMap());
108  for (int i = 0 ; i < A->NumMyRows() ; ++i)
109  X[i] = 1.0* A->RowMatrixRowMap().GID64(i);
110  Y.PutScalar(0.0);
111 
112  Epetra_Vector ExtX_B(B.RowMatrixRowMap());
113  Epetra_Vector ExtY_B(B.RowMatrixRowMap());
114  ExtY_B.PutScalar(0.0);
115 
116  IFPACK_CHK_ERR(B.ImportMultiVector(X,ExtX_B));
117  IFPACK_CHK_ERR(B.Multiply(false,ExtX_B,ExtY_B));
118  IFPACK_CHK_ERR(B.ExportMultiVector(ExtY_B,Y,Add));
119 
120  double Norm_B;
121  Y.Norm2(&Norm_B);
122  if (Comm.MyPID() == 0)
123  cout << "Norm of Y using B = " << Norm_B << endl;
124 
125  // ================================================== //
126  //Build the overlapping matrix as an Epetra_CrsMatrix //
127  // ================================================== //
128 
129  Time.ResetStartTime();
131  *(Ifpack_CreateOverlappingCrsMatrix(&*A,OverlapLevel));
132  if (Comm.MyPID() == 0)
133  cout << "Time to create C = " << Time.ElapsedTime() << endl;
134 
135  // simple checks on global quantities
136  long long NumGlobalRowsC = C.NumGlobalRows64();
137  long long NumGlobalNonzerosC = C.NumGlobalNonzeros64();
138  if (NumGlobalRowsB != NumGlobalRowsC) {
139  std::ostringstream os;
140  os << "NumGlobalRowsB = " << NumGlobalRowsB
141  << " != NumGlobalRowsC = " << NumGlobalRowsC << ".";
142  throw std::logic_error (os.str ());
143  }
144  if (NumGlobalNonzerosB != NumGlobalNonzerosC) {
145  std::ostringstream os;
146  os << "NumGlobalNonzerosB = " << NumGlobalNonzerosB
147  << " != NumGlobalNonzerosC = " << NumGlobalNonzerosC << ".";
148  throw std::logic_error (os.str ());
149  }
150 
151  Epetra_Vector ExtX_C(C.RowMatrixRowMap());
152  Epetra_Vector ExtY_C(C.RowMatrixRowMap());
153  ExtY_C.PutScalar(0.0);
154  Y.PutScalar(0.0);
155 
156  IFPACK_CHK_ERR(C.Multiply(false,X,Y));
157 
158  double Norm_C;
159  Y.Norm2(&Norm_C);
160  if (Comm.MyPID() == 0)
161  cout << "Norm of Y using C = " << Norm_C << endl;
162 
163  if (IFPACK_ABS(Norm_B - Norm_C) > 1e-5)
164  IFPACK_CHK_ERR(-1);
165 
166  // ======================= //
167  // now localize the matrix //
168  // ======================= //
169 
170  Ifpack_LocalFilter D(Teuchos::rcp(&B, false));
171 
172 #ifdef HAVE_MPI
173  MPI_Finalize() ;
174 #endif
175 
176  if (Comm.MyPID() == 0)
177  cout << "Test `TestOverlappingRowMatrix.exe' passed!" << endl;
178 
179  return(EXIT_SUCCESS);
180 }
long long NumGlobalRows64() const
virtual long long NumGlobalNonzeros64() const
Returns the number of nonzero entries in the global matrix.
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
double ElapsedTime(void) const
const Epetra_Map & RowMatrixRowMap() const
int MyPID() const
Epetra_CrsMatrix * Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix *Matrix, const int OverlappingLevel)
Creates an overlapping Epetra_CrsMatrix. Returns 0 if OverlappingLevel is 0.
long long NumGlobalNonzeros64() const
virtual int Multiply(bool TransA, const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_RowMatrix multiplied by a Epetra_MultiVector X in Y.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
int main(int argc, char *argv[])
virtual long long NumGlobalRows64() const
Returns the number of global matrix rows.
int NumProc() const
virtual const Epetra_Map & RowMatrixRowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
#define IFPACK_ABS(x)
int ImportMultiVector(const Epetra_MultiVector &X, Epetra_MultiVector &OvX, Epetra_CombineMode CM=Insert)
int ExportMultiVector(const Epetra_MultiVector &OvX, Epetra_MultiVector &X, Epetra_CombineMode CM=Add)
#define IFPACK_CHK_ERR(ifpack_err)
void ResetStartTime(void)