Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cxx_main.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stratimikos: Thyra-based strategies for linear solvers
5// Copyright (2006) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41//
42// This test uses the MVOPTester.hpp functions to test the Belos adapters
43// to Epetra and Thyra.
44//
45
46#include "Epetra_Map.h"
47#include "Epetra_CrsMatrix.h"
48#ifdef HAVE_MPI
49#include "mpi.h"
50#include "Epetra_MpiComm.h"
51#endif
52#ifndef __cplusplus
53#define __cplusplus
54#endif
55#include "Epetra_Comm.h"
56#include "Epetra_SerialComm.h"
57
58#include "BelosConfigDefs.hpp"
59#include "BelosMVOPTester.hpp"
60#include "BelosEpetraAdapter.hpp"
61
62#ifdef HAVE_EPETRA_THYRA
63#include "BelosThyraAdapter.hpp"
64#include "Thyra_EpetraThyraWrappers.hpp"
65#include "Thyra_EpetraLinearOp.hpp"
66#endif
67
68
69
70int main(int argc, char *argv[])
71{
72
73 using Teuchos::rcp_implicit_cast;
74
75 int i, ierr, gerr;
76 gerr = 0;
77#ifdef NDEBUG
78 (void)ierr; // Eliminate unused warning
79#endif
80
81#ifdef HAVE_MPI
82 // Initialize MPI and setup an Epetra communicator
83 MPI_Init(&argc,&argv);
84 Teuchos::RCP<Epetra_MpiComm> Comm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
85#else
86 // If we aren't using MPI, then setup a serial communicator.
87 Teuchos::RCP<Epetra_SerialComm> Comm = Teuchos::rcp( new Epetra_SerialComm() );
88#endif
89
90
91 // number of global elements
92 int dim = 100;
93 int blockSize = 3;
94
95 // PID info
96 int MyPID = Comm->MyPID();
97 bool verbose = 0;
98
99 if (argc>1) {
100 if (argv[1][0]=='-' && argv[1][1]=='v') {
101 verbose = true;
102 }
103 }
104
105 // Construct a Map that puts approximately the same number of
106 // equations on each processor.
107 Teuchos::RCP<Epetra_Map> Map = Teuchos::rcp( new Epetra_Map(dim, 0, *Comm) );
108
109 // Get update list and number of local equations from newly created Map.
110 int NumMyElements = Map->NumMyElements();
111 std::vector<int> MyGlobalElements(NumMyElements);
112 Map->MyGlobalElements(&MyGlobalElements[0]);
113
114 // Create an integer std::vector NumNz that is used to build the Petra Matrix.
115 // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation
116 // on this processor
117 std::vector<int> NumNz(NumMyElements);
118
119 // We are building a tridiagonal matrix where each row has (-1 2 -1)
120 // So we need 2 off-diagonal terms (except for the first and last equation)
121 for (i=0; i<NumMyElements; i++) {
122 if (MyGlobalElements[i]==0 || MyGlobalElements[i] == dim-1) {
123 NumNz[i] = 2;
124 }
125 else {
126 NumNz[i] = 3;
127 }
128 }
129
130 // Create an Epetra_Matrix
131 Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Map, &NumNz[0]) );
132
133 // Add rows one-at-a-time
134 // Need some vectors to help
135 // Off diagonal Values will always be -1
136 std::vector<double> Values(2);
137 Values[0] = -1.0; Values[1] = -1.0;
138 std::vector<int> Indices(2);
139 double two = 2.0;
140 int NumEntries;
141 for (i=0; i<NumMyElements; i++) {
142 if (MyGlobalElements[i]==0) {
143 Indices[0] = 1;
144 NumEntries = 1;
145 }
146 else if (MyGlobalElements[i] == dim-1) {
147 Indices[0] = dim-2;
148 NumEntries = 1;
149 }
150 else {
151 Indices[0] = MyGlobalElements[i]-1;
152 Indices[1] = MyGlobalElements[i]+1;
153 NumEntries = 2;
154 }
155 ierr = A->InsertGlobalValues(MyGlobalElements[i],NumEntries,&Values[0],&Indices[0]);
156 assert(ierr==0);
157 // Put in the diagonal entry
158 ierr = A->InsertGlobalValues(MyGlobalElements[i],1,&two,&MyGlobalElements[i]);
159 assert(ierr==0);
160 }
161
162 // Finish building the epetra matrix A
163 ierr = A->FillComplete();
164 assert(ierr==0);
165
166 // Create an Belos::EpetraOp from this Epetra_CrsMatrix
167 Teuchos::RCP<Belos::EpetraOp> op = Teuchos::rcp(new Belos::EpetraOp(A));
168
169 // Issue several useful typedefs;
170 // typedef Belos::MultiVec<double> EMV; // unused
171 // typedef Belos::Operator<double> EOP; // unused
172
173 // Create an Epetra_MultiVector for an initial std::vector to start the solver.
174 // Note that this needs to have the same number of columns as the blocksize.
175 Teuchos::RCP<Belos::EpetraMultiVec> ivec = Teuchos::rcp( new Belos::EpetraMultiVec(*Map, blockSize) );
176 ivec->Random();
177
178 // Create an output manager to handle the I/O from the solver
179 Teuchos::RCP<Belos::OutputManager<double> > MyOM = Teuchos::rcp( new Belos::OutputManager<double>( MyPID ) );
180 if (verbose) {
181 MyOM->setVerbosity( Belos::Errors + Belos::Warnings );
182 }
183
184#ifdef HAVE_EPETRA_THYRA
185 typedef Thyra::MultiVectorBase<double> TMVB;
186 typedef Thyra::LinearOpBase<double> TLOB;
187 // create thyra objects from the epetra objects
188
189 // first, a Thyra::VectorSpaceBase
190 Teuchos::RCP<const Thyra::VectorSpaceBase<double> > epetra_vs =
191 Thyra::create_VectorSpace(Map);
192
193 // then, a MultiVectorBase (from the Epetra_MultiVector)
194 Teuchos::RCP<Thyra::MultiVectorBase<double> > thyra_ivec =
195 Thyra::create_MultiVector(rcp_implicit_cast<Epetra_MultiVector>(ivec),epetra_vs);
196
197 // then, a LinearOpBase (from the Epetra_CrsMatrix)
198 Teuchos::RCP<Thyra::LinearOpBase<double> > thyra_op =
199 Teuchos::rcp( new Thyra::EpetraLinearOp(A) );
200
201
202 // test the Thyra adapter multivector
203 ierr = Belos::TestMultiVecTraits<double,TMVB>(MyOM,thyra_ivec);
204 gerr |= ierr;
205 switch (ierr) {
206 case Belos::Ok:
207 if ( verbose && MyPID==0 ) {
208 std::cout << "*** ThyraAdapter PASSED TestMultiVecTraits()" << std::endl;
209 }
210 break;
211 case Belos::Error:
212 if ( verbose && MyPID==0 ) {
213 std::cout << "*** ThyraAdapter FAILED TestMultiVecTraits() ***"
214 << std::endl << std::endl;
215 }
216 break;
217 }
218
219 // test the Thyra adapter operator
220 ierr = Belos::TestOperatorTraits<double,TMVB,TLOB>(MyOM,thyra_ivec,thyra_op);
221 gerr |= ierr;
222 switch (ierr) {
223 case Belos::Ok:
224 if ( verbose && MyPID==0 ) {
225 std::cout << "*** ThyraAdapter PASSED TestOperatorTraits()" << std::endl;
226 }
227 break;
228 case Belos::Error:
229 if ( verbose && MyPID==0 ) {
230 std::cout << "*** ThyraAdapter FAILED TestOperatorTraits() ***"
231 << std::endl << std::endl;
232 }
233 break;
234 }
235#endif
236
237#ifdef HAVE_MPI
238 MPI_Finalize();
239#endif
240
241 if (gerr) {
242 if (verbose && MyPID==0)
243 std::cout << "End Result: TEST FAILED" << std::endl;
244 return -1;
245 }
246 //
247 // Default return value
248 //
249 if (verbose && MyPID==0)
250 std::cout << "End Result: TEST PASSED" << std::endl;
251 return 0;
252
253}
Thyra specializations of MultiVecTraits and OperatorTraits.
int main(int argc, char *argv[])
Definition cxx_main.cpp:70