Thyra Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
createTridiagEpetraLinearOp.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) 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 (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
43 #include "Thyra_EpetraLinearOp.hpp"
44 #include "Epetra_Map.h"
45 #include "Epetra_CrsMatrix.h"
46 #ifdef HAVE_MPI
47 # include "Epetra_MpiComm.h"
48 #else
49 # include "Epetra_SerialComm.h"
50 #endif
51 
54  const int globalDim
55 #ifdef HAVE_MPI
56  ,MPI_Comm mpiComm
57 #endif
58  ,const double diagScale
59  ,const bool verbose
60  ,std::ostream &out
61  )
62 {
63  using Teuchos::RCP;
64  using Teuchos::rcp;
65 
66  //
67  // (A) Create Epetra_Map
68  //
69 
70 #ifdef HAVE_MPI
71  if(verbose) out << "\nCreating Epetra_MpiComm ...\n";
72  Epetra_MpiComm epetra_comm(mpiComm); // Note, Epetra_MpiComm is just a handle class to Epetra_MpiCommData object!
73 #else
74  if(verbose) out << "\nCreating Epetra_SerialComm ...\n";
75  Epetra_SerialComm epetra_comm; // Note, Epetra_SerialComm is just a handle class to Epetra_SerialCommData object!
76 #endif
77  // Create the Epetra_Map object giving it the handle to the Epetra_Comm object
78  const Epetra_Map epetra_map(globalDim,0,epetra_comm); // Note, Epetra_Map is just a handle class to an Epetra_BlockMapData object!
79  // Note that the above Epetra_Map object "copies" the Epetra_Comm object in some
80  // since so that memory mangaement is performed safely.
81 
82  //
83  // (B) Create the tridiagonal Epetra object
84  //
85  // [ 2 -1 ]
86  // [ -1 2 -1 ]
87  // A = [ . . . ]
88  // [ -1 2 -1 ]
89  // [ -1 2 ]
90  //
91 
92  // (B.1) Allocate the Epetra_CrsMatrix object.
93  RCP<Epetra_CrsMatrix> A_epetra = rcp(new Epetra_CrsMatrix(::Copy,epetra_map,3));
94  // Note that Epetra_CrsMatrix is *not* a handle object so have to use RCP above.
95  // Also note that the Epetra_Map object is copied in some sence by the Epetra_CrsMatrix
96  // object so the memory managment is handled safely.
97 
98  // (B.2) Get the indexes of the rows on this processor
99  const int numMyElements = epetra_map.NumMyElements();
100  std::vector<int> myGlobalElements(numMyElements);
101  epetra_map.MyGlobalElements(&myGlobalElements[0]);
102 
103  // (B.3) Fill the local matrix entries one row at a time.
104  // Note, we set up Epetra_Map above to use zero-based indexing and that is what we must use below:
105  const double offDiag = -1.0, diag = 2.0*diagScale;
106  int numEntries; double values[3]; int indexes[3];
107  for( int k = 0; k < numMyElements; ++k ) {
108  const int rowIndex = myGlobalElements[k];
109  if( rowIndex == 0 ) { // First row
110  numEntries = 2;
111  values[0] = diag; values[1] = offDiag;
112  indexes[0] = 0; indexes[1] = 1;
113  }
114  else if( rowIndex == globalDim - 1 ) { // Last row
115  numEntries = 2;
116  values[0] = offDiag; values[1] = diag;
117  indexes[0] = globalDim-2; indexes[1] = globalDim-1;
118  }
119  else { // Middle rows
120  numEntries = 3;
121  values[0] = offDiag; values[1] = diag; values[2] = offDiag;
122  indexes[0] = rowIndex-1; indexes[1] = rowIndex; indexes[2] = rowIndex+1;
123  }
124  TEUCHOS_TEST_FOR_EXCEPT( 0!=A_epetra->InsertGlobalValues(rowIndex,numEntries,values,indexes) );
125  }
126 
127  // (B.4) Finish the construction of the Epetra_CrsMatrix
128  TEUCHOS_TEST_FOR_EXCEPT( 0!=A_epetra->FillComplete() );
129 
130  // Return the Epetra_Operator object
131  return A_epetra;
132 
133  // Note that when this function returns the returned RCP-wrapped
134  // Epetra_Operator object will own all of the Epetra objects that went into
135  // its construction and these objects will stay around until all of the
136  // RCP objects to the allocated Epetra_Operator object are removed
137  // and destruction occurs!
138 
139 } // end createTridiagLinearOp()
Teuchos::RCP< Epetra_Operator > createTridiagEpetraLinearOp(const int globalDim, const double diagScale, const bool verbose, std::ostream &out)
This function generates a tridiagonal linear operator using Epetra.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)