Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/OverlappingRowMatrix/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"
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_RCP.hpp"
61#include "Ifpack_LocalFilter.h"
62#include "Ifpack_Utils.h"
63#include <sstream>
64#include <stdexcept>
65
66int main(int argc, char *argv[])
67{
68#ifdef HAVE_MPI
69 MPI_Init(&argc,&argv);
70 Epetra_MpiComm Comm( MPI_COMM_WORLD );
71#else
73#endif
74
75 if (Comm.NumProc() == 1)
76 {
77#ifdef HAVE_MPI
78 MPI_Finalize();
79#endif
80 cout << "Test `TestOverlappingRowMatrix.exe' passed!" << endl;
81 exit(EXIT_SUCCESS);
82 }
83
84 Teuchos::ParameterList GaleriList;
85 int nx = 100;
86 GaleriList.set("n", nx * nx);
87 GaleriList.set("nx", nx);
88 GaleriList.set("ny", nx);
89 Teuchos::RCP<Epetra_Map> Map
90 (Galeri::CreateMap ("Linear", Comm, GaleriList));
91 Teuchos::RCP<Epetra_CrsMatrix> A
92 (Galeri::CreateCrsMatrix ("Laplace2D", Map.getRawPtr (), GaleriList));
93
94 int OverlapLevel = 5;
95 Epetra_Time Time(Comm);
96
97 // ======================================== //
98 // Build the overlapping matrix using class //
99 // Ifpack_OverlappingRowMatrix. //
100 // ======================================== //
101
102 Time.ResetStartTime();
103 Ifpack_OverlappingRowMatrix B(A,OverlapLevel);
104 if (Comm.MyPID() == 0)
105 cout << "Time to create B = " << Time.ElapsedTime() << endl;
106
107 int NumGlobalRowsB = B.NumGlobalRows();
108 int NumGlobalNonzerosB = B.NumGlobalNonzeros();
109
110 Epetra_Vector X(A->RowMatrixRowMap());
111 Epetra_Vector Y(A->RowMatrixRowMap());
112 for (int i = 0 ; i < A->NumMyRows() ; ++i)
113 X[i] = 1.0* A->RowMatrixRowMap().GID(i);
114 Y.PutScalar(0.0);
115
116 Epetra_Vector ExtX_B(B.RowMatrixRowMap());
117 Epetra_Vector ExtY_B(B.RowMatrixRowMap());
118 ExtY_B.PutScalar(0.0);
119
121 IFPACK_CHK_ERR(B.Multiply(false,ExtX_B,ExtY_B));
123
124 double Norm_B;
125 Y.Norm2(&Norm_B);
126 if (Comm.MyPID() == 0)
127 cout << "Norm of Y using B = " << Norm_B << endl;
128
129 // ================================================== //
130 //Build the overlapping matrix as an Epetra_CrsMatrix //
131 // ================================================== //
132
133 Time.ResetStartTime();
134 Epetra_CrsMatrix& C =
135 *(Ifpack_CreateOverlappingCrsMatrix(&*A,OverlapLevel));
136 if (Comm.MyPID() == 0)
137 cout << "Time to create C = " << Time.ElapsedTime() << endl;
138
139 // simple checks on global quantities
140 int NumGlobalRowsC = C.NumGlobalRows();
141 int NumGlobalNonzerosC = C.NumGlobalNonzeros();
142
143 if (NumGlobalRowsB != NumGlobalRowsC) {
144 std::ostringstream os;
145 os << "NumGlobalRowsB = " << NumGlobalRowsB
146 << " != NumGlobalRowsC = " << NumGlobalRowsC
147 << "." << std::endl;
148 throw std::logic_error (os.str ());
149 }
150 if (NumGlobalNonzerosB != NumGlobalNonzerosC) {
151 std::ostringstream os;
152 os << "NumGlobalNonzerosB = " << NumGlobalNonzerosB
153 << " != NumGlobalNonzerosC = " << NumGlobalNonzerosC
154 << "." << std::endl;
155 throw std::logic_error (os.str ());
156 }
157
158 Epetra_Vector ExtX_C(C.RowMatrixRowMap());
159 Epetra_Vector ExtY_C(C.RowMatrixRowMap());
160 ExtY_C.PutScalar(0.0);
161 Y.PutScalar(0.0);
162
163 IFPACK_CHK_ERR(C.Multiply(false,X,Y));
164
165 double Norm_C;
166 Y.Norm2(&Norm_C);
167 if (Comm.MyPID() == 0)
168 cout << "Norm of Y using C = " << Norm_C << endl;
169
170 if (IFPACK_ABS(Norm_B - Norm_C) > 1e-5)
171 IFPACK_CHK_ERR(-1);
172
173 // ======================= //
174 // now localize the matrix //
175 // ======================= //
176
177 Ifpack_LocalFilter D(Teuchos::rcp(&B, false));
178
179#ifdef HAVE_MPI
180 MPI_Finalize() ;
181#endif
182
183 if (Comm.MyPID() == 0)
184 cout << "Test `TestOverlappingRowMatrix.exe' passed!" << endl;
185
186 return(EXIT_SUCCESS);
187}
Add
#define IFPACK_CHK_ERR(ifpack_err)
#define IFPACK_ABS(x)
Epetra_CrsMatrix * Ifpack_CreateOverlappingCrsMatrix(const Epetra_RowMatrix *Matrix, const int OverlappingLevel)
Creates an overlapping Epetra_CrsMatrix. Returns 0 if OverlappingLevel is 0.
int NumGlobalNonzeros() const
int NumGlobalRows() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_Map & RowMatrixRowMap() const
int NumProc() const
int MyPID() const
int PutScalar(double ScalarConstant)
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
Ifpack_OverlappingRowMatrix: matrix with ghost rows, based on Epetra_RowMatrix.
int ExportMultiVector(const Epetra_MultiVector &OvX, Epetra_MultiVector &X, Epetra_CombineMode CM=Add)
int ImportMultiVector(const Epetra_MultiVector &X, Epetra_MultiVector &OvX, Epetra_CombineMode CM=Insert)
virtual int NumGlobalNonzeros() const
Returns the number of nonzero entries in the global matrix.
virtual int NumGlobalRows() const
Returns the number of global matrix rows.
virtual const Epetra_Map & RowMatrixRowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
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.
int main(int argc, char *argv[])