Thyra Package Browser (Single Doxygen Collection) Version of the Day
createTridiagEpetraLinearOp.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include "createTridiagEpetraLinearOp.hpp"
00043 #include "Thyra_EpetraLinearOp.hpp"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_CrsMatrix.h"
00046 #ifdef HAVE_MPI
00047 # include "Epetra_MpiComm.h"
00048 #else
00049 # include "Epetra_SerialComm.h"
00050 #endif
00051 
00052 Teuchos::RCP<Epetra_Operator>
00053 createTridiagEpetraLinearOp(
00054   const int      globalDim
00055 #ifdef HAVE_MPI
00056   ,MPI_Comm      mpiComm
00057 #endif
00058   ,const double  diagScale
00059   ,const bool    verbose
00060   ,std::ostream  &out
00061   )
00062 {
00063   using Teuchos::RCP;
00064   using Teuchos::rcp;
00065 
00066   //
00067   // (A) Create Epetra_Map
00068   //
00069 
00070 #ifdef HAVE_MPI
00071   if(verbose) out << "\nCreating Epetra_MpiComm ...\n";
00072   Epetra_MpiComm epetra_comm(mpiComm); // Note, Epetra_MpiComm is just a handle class to Epetra_MpiCommData object!
00073 #else
00074   if(verbose) out << "\nCreating Epetra_SerialComm ...\n";
00075   Epetra_SerialComm epetra_comm; // Note, Epetra_SerialComm is just a handle class to Epetra_SerialCommData object!
00076 #endif
00077   // Create the Epetra_Map object giving it the handle to the Epetra_Comm object
00078   const Epetra_Map epetra_map(globalDim,0,epetra_comm); // Note, Epetra_Map is just a handle class to an Epetra_BlockMapData object!
00079   // Note that the above Epetra_Map object "copies" the Epetra_Comm object in some
00080   // since so that memory mangaement is performed safely.
00081 
00082   //
00083   // (B) Create the tridiagonal Epetra object
00084   //
00085   //       [  2  -1             ]
00086   //       [ -1   2  -1         ]
00087   //  A =  [      .   .   .     ]
00088   //       [          -1  2  -1 ]
00089   //       [             -1   2 ]
00090   //
00091   
00092   // (B.1) Allocate the Epetra_CrsMatrix object.
00093   RCP<Epetra_CrsMatrix> A_epetra = rcp(new Epetra_CrsMatrix(::Copy,epetra_map,3));
00094   // Note that Epetra_CrsMatrix is *not* a handle object so have to use RCP above.
00095   // Also note that the Epetra_Map object is copied in some sence by the Epetra_CrsMatrix
00096   // object so the memory managment is handled safely.
00097 
00098   // (B.2) Get the indexes of the rows on this processor
00099   const int numMyElements = epetra_map.NumMyElements();
00100   std::vector<int> myGlobalElements(numMyElements);
00101   epetra_map.MyGlobalElements(&myGlobalElements[0]);
00102 
00103   // (B.3) Fill the local matrix entries one row at a time.
00104   // Note, we set up Epetra_Map above to use zero-based indexing and that is what we must use below:
00105   const double offDiag = -1.0, diag = 2.0*diagScale;
00106   int numEntries; double values[3]; int indexes[3];
00107   for( int k = 0; k < numMyElements; ++k ) {
00108     const int rowIndex = myGlobalElements[k];
00109     if( rowIndex == 0 ) {                     // First row
00110       numEntries = 2;
00111       values[0]  = diag;             values[1]  = offDiag;
00112       indexes[0] = 0;                indexes[1] = 1; 
00113     }
00114     else if( rowIndex == globalDim - 1 ) {    // Last row
00115       numEntries = 2;
00116       values[0]  = offDiag;         values[1]  = diag;
00117       indexes[0] = globalDim-2;     indexes[1] = globalDim-1; 
00118     }
00119     else {                                    // Middle rows
00120       numEntries = 3;
00121       values[0]  = offDiag;         values[1]  = diag;          values[2]  = offDiag;
00122       indexes[0] = rowIndex-1;      indexes[1] = rowIndex;      indexes[2] = rowIndex+1; 
00123     }
00124     TEUCHOS_TEST_FOR_EXCEPT( 0!=A_epetra->InsertGlobalValues(rowIndex,numEntries,values,indexes) );
00125   }
00126 
00127   // (B.4) Finish the construction of the Epetra_CrsMatrix
00128   TEUCHOS_TEST_FOR_EXCEPT( 0!=A_epetra->FillComplete() );
00129 
00130   // Return the Epetra_Operator object
00131   return A_epetra;
00132 
00133   // Note that when this function returns the returned RCP-wrapped
00134   // Epetra_Operator object will own all of the Epetra objects that went into
00135   // its construction and these objects will stay around until all of the
00136   // RCP objects to the allocated Epetra_Operator object are removed
00137   // and destruction occurs!
00138 
00139 } // end createTridiagLinearOp()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines