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