Thyra Package Browser (Single Doxygen Collection) Version of the Day
createTridiagEpetraLinearOp.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //               Epetra: Linear Algebra Services Package 
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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "createTridiagEpetraLinearOp.hpp"
00030 #include "Thyra_EpetraLinearOp.hpp"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_CrsMatrix.h"
00033 #ifdef HAVE_MPI
00034 # include "Epetra_MpiComm.h"
00035 #else
00036 # include "Epetra_SerialComm.h"
00037 #endif
00038 
00039 Teuchos::RCP<Epetra_Operator>
00040 createTridiagEpetraLinearOp(
00041   const int      globalDim
00042 #ifdef HAVE_MPI
00043   ,MPI_Comm      mpiComm
00044 #endif
00045   ,const double  diagScale
00046   ,const bool    verbose
00047   ,std::ostream  &out
00048   )
00049 {
00050   using Teuchos::RCP;
00051   using Teuchos::rcp;
00052 
00053   //
00054   // (A) Create Epetra_Map
00055   //
00056 
00057 #ifdef HAVE_MPI
00058   if(verbose) out << "\nCreating Epetra_MpiComm ...\n";
00059   Epetra_MpiComm epetra_comm(mpiComm); // Note, Epetra_MpiComm is just a handle class to Epetra_MpiCommData object!
00060 #else
00061   if(verbose) out << "\nCreating Epetra_SerialComm ...\n";
00062   Epetra_SerialComm epetra_comm; // Note, Epetra_SerialComm is just a handle class to Epetra_SerialCommData object!
00063 #endif
00064   // Create the Epetra_Map object giving it the handle to the Epetra_Comm object
00065   const Epetra_Map epetra_map(globalDim,0,epetra_comm); // Note, Epetra_Map is just a handle class to an Epetra_BlockMapData object!
00066   // Note that the above Epetra_Map object "copies" the Epetra_Comm object in some
00067   // since so that memory mangaement is performed safely.
00068 
00069   //
00070   // (B) Create the tridiagonal Epetra object
00071   //
00072   //       [  2  -1             ]
00073   //       [ -1   2  -1         ]
00074   //  A =  [      .   .   .     ]
00075   //       [          -1  2  -1 ]
00076   //       [             -1   2 ]
00077   //
00078   
00079   // (B.1) Allocate the Epetra_CrsMatrix object.
00080   RCP<Epetra_CrsMatrix> A_epetra = rcp(new Epetra_CrsMatrix(::Copy,epetra_map,3));
00081   // Note that Epetra_CrsMatrix is *not* a handle object so have to use RCP above.
00082   // Also note that the Epetra_Map object is copied in some sence by the Epetra_CrsMatrix
00083   // object so the memory managment is handled safely.
00084 
00085   // (B.2) Get the indexes of the rows on this processor
00086   const int numMyElements = epetra_map.NumMyElements();
00087   std::vector<int> myGlobalElements(numMyElements);
00088   epetra_map.MyGlobalElements(&myGlobalElements[0]);
00089 
00090   // (B.3) Fill the local matrix entries one row at a time.
00091   // Note, we set up Epetra_Map above to use zero-based indexing and that is what we must use below:
00092   const double offDiag = -1.0, diag = 2.0*diagScale;
00093   int numEntries; double values[3]; int indexes[3];
00094   for( int k = 0; k < numMyElements; ++k ) {
00095     const int rowIndex = myGlobalElements[k];
00096     if( rowIndex == 0 ) {                     // First row
00097       numEntries = 2;
00098       values[0]  = diag;             values[1]  = offDiag;
00099       indexes[0] = 0;                indexes[1] = 1; 
00100     }
00101     else if( rowIndex == globalDim - 1 ) {    // Last row
00102       numEntries = 2;
00103       values[0]  = offDiag;         values[1]  = diag;
00104       indexes[0] = globalDim-2;     indexes[1] = globalDim-1; 
00105     }
00106     else {                                    // Middle rows
00107       numEntries = 3;
00108       values[0]  = offDiag;         values[1]  = diag;          values[2]  = offDiag;
00109       indexes[0] = rowIndex-1;      indexes[1] = rowIndex;      indexes[2] = rowIndex+1; 
00110     }
00111     TEST_FOR_EXCEPT( 0!=A_epetra->InsertGlobalValues(rowIndex,numEntries,values,indexes) );
00112   }
00113 
00114   // (B.4) Finish the construction of the Epetra_CrsMatrix
00115   TEST_FOR_EXCEPT( 0!=A_epetra->FillComplete() );
00116 
00117   // Return the Epetra_Operator object
00118   return A_epetra;
00119 
00120   // Note that when this function returns the returned RCP-wrapped
00121   // Epetra_Operator object will own all of the Epetra objects that went into
00122   // its construction and these objects will stay around until all of the
00123   // RCP objects to the allocated Epetra_Operator object are removed
00124   // and destruction occurs!
00125 
00126 } // end createTridiagLinearOp()
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines