Ifpack Package Browser (Single Doxygen Collection) Development
test/Container/cxx_main.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                IFPACK
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 "Ifpack_ConfigDefs.h"
00030 
00031 #ifdef HAVE_MPI
00032 #include "Epetra_MpiComm.h"
00033 #endif
00034 #include "Epetra_SerialComm.h"
00035 #include "Epetra_CrsMatrix.h"
00036 #include "Epetra_Vector.h"
00037 #include "Epetra_LinearProblem.h"
00038 #include "Galeri_Maps.h"
00039 #include "Galeri_CrsMatrices.h"
00040 #include "Galeri_Utils.h"
00041 #include "Teuchos_ParameterList.hpp"
00042 #include "Teuchos_RefCountPtr.hpp"
00043 #include "Ifpack_DenseContainer.h"
00044 #include "Ifpack_SparseContainer.h"
00045 #include "Ifpack_Amesos.h"
00046 #include "Ifpack_LocalFilter.h"
00047 
00048 static bool verbose = false;
00049 
00050 // ======================================================================
00051 bool TestContainer(string Type, const Teuchos::RefCountPtr<Epetra_RowMatrix>& A)
00052 {
00053   int NumVectors = 3;
00054   int NumMyRows = A->NumMyRows();
00055 
00056   Epetra_MultiVector LHS_exact(A->RowMatrixRowMap(), NumVectors);
00057   Epetra_MultiVector LHS(A->RowMatrixRowMap(), NumVectors);
00058   Epetra_MultiVector RHS(A->RowMatrixRowMap(), NumVectors);
00059   LHS_exact.Random(); LHS.PutScalar(0.0); 
00060   A->Multiply(false, LHS_exact, RHS);
00061 
00062   Epetra_LinearProblem Problem(&*A, &LHS, &RHS);
00063 
00064   if (verbose) {
00065     cout << "Container type = " << Type << endl;
00066     cout << "NumMyRows = " << NumMyRows << ", NumVectors = " << NumVectors << endl;
00067   }
00068   LHS.PutScalar(0.0);
00069   
00070   Teuchos::RefCountPtr<Ifpack_Container> Container;
00071 
00072   if (Type == "dense")
00073     Container = Teuchos::rcp( new Ifpack_DenseContainer(A->NumMyRows(), NumVectors) );
00074   else
00075     Container = Teuchos::rcp( new Ifpack_SparseContainer<Ifpack_Amesos>(A->NumMyRows(), NumVectors) );
00076 
00077   assert (Container != Teuchos::null);
00078 
00079   IFPACK_CHK_ERR(Container->Initialize());
00080   // set as ID all the local rows of A
00081   for (int i = 0 ; i < A->NumMyRows() ; ++i)
00082     Container->ID(i) = i;
00083 
00084   // extract submatrix (in this case, the entire matrix)
00085   // and complete setup
00086   IFPACK_CHK_ERR(Container->Compute(*A));
00087 
00088   // set the RHS and LHS
00089   for (int i = 0 ; i < A->NumMyRows() ; ++i)
00090     for (int j = 0 ; j < NumVectors ; ++j) {
00091       Container->RHS(i,j) = RHS[j][i];
00092       Container->LHS(i,j) = LHS[j][i];
00093     }
00094   
00095   // set parameters (empty for dense containers)
00096   Teuchos::ParameterList List;
00097   List.set("amesos: solver type", Type);
00098   IFPACK_CHK_ERR(Container->SetParameters(List));
00099 
00100   // solve the linear system
00101   IFPACK_CHK_ERR(Container->ApplyInverse());
00102 
00103   // get the computed solution, store it in LHS
00104   for (int i = 0 ; i < A->NumMyRows() ; ++i)
00105     for (int j = 0 ; j < NumVectors ; ++j) {
00106        LHS[j][i] = Container->LHS(i,j);
00107     }
00108 
00109   double residual = Galeri::ComputeNorm(&LHS, &LHS_exact);
00110 
00111   if (A->Comm().MyPID() == 0 && verbose) {
00112     cout << "||x_exact - x||_2 = " << residual << endl;
00113     cout << *Container;
00114   }
00115 
00116   bool passed = false;
00117   if (residual < 1e-5)
00118     passed = true;
00119 
00120   return(passed);
00121 }
00122 
00123 // ======================================================================
00124 int main(int argc, char *argv[])
00125 {
00126 
00127 #ifdef HAVE_MPI
00128   MPI_Init(&argc,&argv);
00129   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00130 #else
00131   Epetra_SerialComm Comm;
00132 #endif
00133   Epetra_SerialComm SerialComm;
00134 
00135   verbose = (Comm.MyPID() == 0);
00136 
00137   const int nx = 30;
00138   Teuchos::ParameterList GaleriList;
00139   GaleriList.set("n", nx * nx);
00140   GaleriList.set("nx", nx);
00141   GaleriList.set("ny", nx);
00142 
00143   Teuchos::RefCountPtr<Epetra_Map> Map = Teuchos::rcp( Galeri::CreateMap("Linear", Comm, GaleriList) );
00144   Teuchos::RefCountPtr<Epetra_RowMatrix> Matrix = Teuchos::rcp( Galeri::CreateCrsMatrix("Laplace2D", &*Map, GaleriList) );
00145   
00146   Teuchos::RefCountPtr<Ifpack_LocalFilter> LocalMatrix = Teuchos::rcp( new Ifpack_LocalFilter(Matrix) );
00147   int TestPassed = true;
00148 
00149   if (!TestContainer("dense",LocalMatrix))  TestPassed = false;
00150   if (!TestContainer("sparse",LocalMatrix)) TestPassed = false;
00151 
00152   if (TestPassed)
00153   {
00154     if (verbose)
00155       cout << "Test `TestContainer.exe' passed!" << endl;
00156   }
00157   else {
00158     cout << "Test `TestContainer.exe' FAILED!" << endl;
00159     exit(EXIT_FAILURE);
00160   }
00161   
00162 #ifdef HAVE_MPI
00163   MPI_Finalize(); 
00164 #endif
00165 
00166   return(EXIT_SUCCESS);
00167 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines