Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
test/Container/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#endif
48#include "Epetra_SerialComm.h"
49#include "Epetra_CrsMatrix.h"
50#include "Epetra_Vector.h"
51#include "Epetra_LinearProblem.h"
52#include "Galeri_Maps.h"
53#include "Galeri_CrsMatrices.h"
54#include "Galeri_Utils.h"
55#include "Teuchos_ParameterList.hpp"
56#include "Teuchos_RefCountPtr.hpp"
59#include "Ifpack_Amesos.h"
60#include "Ifpack_LocalFilter.h"
61
62static bool verbose = false;
63
64// ======================================================================
65bool TestContainer(std::string Type, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A)
66{
67 int NumVectors = 3;
68 int NumMyRows = A->NumMyRows();
69
70 Epetra_MultiVector LHS_exact(A->RowMatrixRowMap(), NumVectors);
71 Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
72 Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
73 LHS_exact.Random(); LHS.PutScalar(0.0);
74 A->Multiply(false, LHS_exact, RHS);
75
76 Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
77
78 if (verbose) {
79 cout << "Container type = " << Type << endl;
80 cout << "NumMyRows = " << NumMyRows << ", NumVectors = " << NumVectors << endl;
81 }
82 LHS.PutScalar(0.0);
83
84 Teuchos::RefCountPtr<Ifpack_Container> Container;
85
86 if (Type == "dense")
87 Container = Teuchos::rcp( new Ifpack_DenseContainer(A->NumMyRows(), NumVectors) );
88 else
89 Container = Teuchos::rcp( new Ifpack_SparseContainer<Ifpack_Amesos>(A->NumMyRows(), NumVectors) );
90
91 assert (Container != Teuchos::null);
92
93 IFPACK_CHK_ERR(Container->Initialize());
94 // set as ID all the local rows of A
95 for (int i = 0 ; i < A->NumMyRows() ; ++i)
96 Container->ID(i) = i;
97
98 // extract submatrix (in this case, the entire matrix)
99 // and complete setup
100 IFPACK_CHK_ERR(Container->Compute(*A));
101
102 // set the RHS and LHS
103 for (int i = 0 ; i < A->NumMyRows() ; ++i)
104 for (int j = 0 ; j < NumVectors ; ++j) {
105 Container->RHS(i,j) = RHS[j][i];
106 Container->LHS(i,j) = LHS[j][i];
107 }
108
109 // set parameters (empty for dense containers)
110 Teuchos::ParameterList List;
111 List.set("amesos: solver type", Type);
112 IFPACK_CHK_ERR(Container->SetParameters(List));
113
114 // solve the linear system
115 IFPACK_CHK_ERR(Container->ApplyInverse());
116
117 // get the computed solution, store it in LHS
118 for (int i = 0 ; i < A->NumMyRows() ; ++i)
119 for (int j = 0 ; j < NumVectors ; ++j) {
120 LHS[j][i] = Container->LHS(i,j);
121 }
122
123 double residual = Galeri::ComputeNorm(&LHS, &LHS_exact);
124
125 if (A->Comm().MyPID() == 0 && verbose) {
126 cout << "||x_exact - x||_2 = " << residual << endl;
127 cout << *Container;
128 }
129
130 bool passed = false;
131 if (residual < 1e-5)
132 passed = true;
133
134 return(passed);
135}
136
137// ======================================================================
138int main(int argc, char *argv[])
139{
140
141#ifdef HAVE_MPI
142 MPI_Init(&argc,&argv);
143 Epetra_MpiComm Comm(MPI_COMM_WORLD);
144#else
146#endif
147 Epetra_SerialComm SerialComm;
148
149 verbose = (Comm.MyPID() == 0);
150
151 const int nx = 30;
152 Teuchos::ParameterList GaleriList;
153 GaleriList.set("n", nx * nx);
154 GaleriList.set("nx", nx);
155 GaleriList.set("ny", nx);
156
157 Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Linear", Comm, GaleriList) );
158 Teuchos::RefCountPtr<Epetra_RowMatrix> Matrix = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
159
160 Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalMatrix = Teuchos::rcp( new Ifpack_LocalFilter(Matrix) );
161 int TestPassed = true;
162
163 if (!TestContainer("dense",LocalMatrix)) TestPassed = false;
164 if (!TestContainer("sparse",LocalMatrix)) TestPassed = false;
165
166 // add TriDi
167
168 if (TestPassed)
169 {
170 if (verbose)
171 cout << "Test `TestContainer.exe' passed!" << endl;
172 }
173 else {
174 cout << "Test `TestContainer.exe' FAILED!" << endl;
175 exit(EXIT_FAILURE);
176 }
177
178#ifdef HAVE_MPI
179 MPI_Finalize();
180#endif
181
182 return(EXIT_SUCCESS);
183}
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60
int MyPID() const
int PutScalar(double ScalarConstant)
Ifpack_DenseContainer: a class to define containers for dense matrices.
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
Ifpack_SparseContainer: a class for storing and solving linear systems using sparse matrices.
const int NumVectors
Definition: performance.cpp:71
bool TestContainer(std::string Type, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A)
int main(int argc, char *argv[])
static bool verbose