Thyra Package Browser (Single Doxygen Collection)  Version of the Day
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 
52 Teuchos::RCP<Epetra_Operator>
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.