Amesos Package Browser (Single Doxygen Collection) Development
example_MC64.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                Amesos: Interface to Direct Solver Libraries
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 "Amesos_ConfigDefs.h"
00030 
00031 #ifdef HAVE_MPI
00032 #include "Epetra_MpiComm.h"
00033 #else
00034 #include "Epetra_SerialComm.h"
00035 #endif
00036 #include "Epetra_CrsMatrix.h"
00037 #include "Epetra_Map.h"
00038 #include "Amesos_MC64.h"
00039 
00040 class Ifpack_ReorderOperator 
00041 {
00042  public:
00043   Ifpack_ReorderOperator(const int size, int* perm, double* scale) :
00044     size_(size),
00045     perm_(perm),
00046     scale_(scale)
00047   { }
00048 
00049   int Apply(double* x, double* y)
00050   {
00051     for (int i = 0 ; i < size_ ; ++i)
00052       y[i] = scale_[i] * x[perm_[i]];
00053   }
00054  private:
00055   const int size_;
00056   int* perm_;
00057   double* scale_;
00058 };
00059 
00060 // Creates the following matrix:
00061 //
00062 // | 3.00  5.00           |
00063 // | 2.00  3.00      0.00 |
00064 // | 3.00       4.00 6.00 |
00065 // |            1.00 2.00 |
00066 //
00067 // as reported in the HSL MC64 version 1.3.0 user's guide. This simple driver
00068 // calls MC64 using JOB=5, then prints on screen the CPERM and DW vectors. The
00069 // reordered and scaled matrix looks like:
00070 //
00071 // | 1.00 1.00           |
00072 // | 0.90 1.00      0.00 |
00073 // |      1.00 1.00 1.00 |
00074 // |           0.75 1.00 |
00075 //
00076 // \warning The interface between Amesos and MC64 works for serial
00077 // computations only.
00078 //
00079 // \author Marzio Sala, ETHZ.
00080 //
00081 // \date Last updated on 05-Feb-06.
00082 
00083 int main(int argc, char *argv[])
00084 {
00085   // initialize MPI and Epetra communicator
00086 #ifdef HAVE_MPI
00087   MPI_Init(&argc,&argv);
00088   Epetra_MpiComm Comm(MPI_COMM_WORLD);
00089 #else
00090   Epetra_SerialComm Comm;
00091 #endif
00092   
00093   if (Comm.NumProc() != 1)
00094   {
00095     std::cerr << "This example can be run with one processor only!" << std::endl;
00096 #ifdef HAVE_MPI
00097     MPI_Finalize();
00098 #endif
00099     exit(EXIT_SUCCESS);
00100   }
00101 
00102   int NumMyRows = 4;
00103   Epetra_Map Map(NumMyRows, 0, Comm);
00104   Epetra_CrsMatrix A(Copy, Map, 0);
00105 
00106   int col;
00107   double val;
00108 
00109   col = 0; val = 3.0;
00110   A.InsertGlobalValues(0, 1, &val, &col);
00111 
00112   col = 1; val = 5.0;
00113   A.InsertGlobalValues(0, 1, &val, &col);
00114 
00115   col = 0; val = 2.0;
00116   A.InsertGlobalValues(1, 1, &val, &col);
00117 
00118   col = 1; val = 3.0;
00119   A.InsertGlobalValues(1, 1, &val, &col);
00120 
00121   col = 3; val = 0.0;
00122   A.InsertGlobalValues(1, 1, &val, &col);
00123 
00124   col = 0; val = 3.0;
00125   A.InsertGlobalValues(2, 1, &val, &col);
00126 
00127   col = 2; val = 4.0;
00128   A.InsertGlobalValues(2, 1, &val, &col);
00129 
00130   col = 3; val = 6.0;
00131   A.InsertGlobalValues(2, 1, &val, &col);
00132 
00133   col = 2; val = 1.0;
00134   A.InsertGlobalValues(3, 1, &val, &col);
00135 
00136   col = 3; val = 2.0;
00137   A.InsertGlobalValues(3, 1, &val, &col);
00138 
00139   A.FillComplete();
00140 
00141   // Creates an instance of the MC64 reordering and scaling interface
00142   // and computes the (column) reordering and scaling, using JOB=5
00143   Amesos_MC64 MC64(A, 5);
00144 
00145   // checks the return value
00146   std::cout << "INFO(1) = " << MC64.GetINFO(1) << std::endl;
00147 
00148   // Gets the pointer to reordering (CPERM) and scaling (DW). Both
00149   // vectors are allocated and free'd within Amesos_MC64.
00150   int*    CPERM = MC64.GetCPERM();
00151   double* DW    = MC64.GetDW();
00152 
00153   for (int i = 0 ; i < A.NumMyRows() ; ++i)
00154     std::cout << "CPERM[" << i << "] = " << CPERM[i] << std::endl;
00155 
00156   for (int i = 0 ; i < A.NumMyRows() * 2 ; ++i)
00157     std::cout << "DW[" << i << "] = " << DW[i] << std::endl;
00158 
00159 
00160   Ifpack_ReorderOperator RowPerm(4, MC64.GetRowPerm(), MC64.GetRowScaling());
00161   Ifpack_ReorderOperator ColPerm(4, MC64.GetColPerm(), MC64.GetColScaling());
00162 
00163 #ifdef HAVE_MPI
00164   MPI_Finalize() ; 
00165 #endif
00166 
00167   return(EXIT_SUCCESS);
00168 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines