Thyra Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)