Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
57 #include "Ifpack_DenseContainer.h"
58 #include "Ifpack_SparseContainer.h"
59 #include "Ifpack_Amesos.h"
60 #include "Ifpack_LocalFilter.h"
61 
62 static bool verbose = false;
63 
64 // ======================================================================
65 bool 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)
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 // ======================================================================
141 int 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
150  Epetra_SerialComm Comm;
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 }
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
const int NumVectors
Definition: performance.cpp:71
int MyPID() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Ifpack_SparseContainer: a class for storing and solving linear systems using sparse matrices...
Ifpack_LocalFilter a class for light-weight extraction of the submatrix corresponding to local rows a...
int main(int argc, char *argv[])
Ifpack_DenseContainer: a class to define containers for dense matrices.
bool TestContainer(std::string Type, const Teuchos::RefCountPtr< Epetra_RowMatrix > &A)
#define IFPACK_CHK_ERR(ifpack_err)
#define RHS(a)
Definition: MatGenFD.c:60