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