ShyLU  Version of the Day
shylu.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // ShyLU: Hybrid preconditioner package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
43 
56 #include <assert.h>
57 #include <iostream>
58 #include <sstream>
59 
60 /* This define will make S block diagonal, To make this happen even some
61  zero columns/rows of C and R will be stored.
62  */
63 #define BLOCK_DIAGONAL_Si
64 
65 #include "shylu_util.h"
66 
67 // Epetra includes
68 #ifdef HAVE_SHYLUCORE_MPI
69 # include "Epetra_MpiComm.h"
70 #endif // HAVE_SHYLUCORE_MPI
71 #include "Epetra_SerialComm.h"
72 #include "Epetra_Time.h"
73 #include "Epetra_CrsMatrix.h"
74 #include "Epetra_Map.h"
75 #include "Epetra_MultiVector.h"
76 #include "Epetra_LinearProblem.h"
77 #include "Epetra_Import.h"
78 #include "Epetra_Export.h"
79 
80 // Teuchos includes
81 #include "Teuchos_GlobalMPISession.hpp"
82 #include "Teuchos_XMLParameterListHelpers.hpp"
83 #include "Teuchos_LAPACK.hpp"
84 
85 // EpetraExt includes
86 #include "EpetraExt_RowMatrixOut.h"
87 #include "EpetraExt_MultiVectorOut.h"
88 #include "EpetraExt_CrsMatrixIn.h"
89 
90 // Amesos includes
91 #include "Amesos.h"
92 #include "Amesos_BaseSolver.h"
93 
94 
95 // Amesos2 includes
96 #ifdef HAVE_SHYLUCORE_AMESOS2
97 #include <Amesos2.hpp>
98 #endif
99 
100 // Tpetra includes
101 #ifdef HAVE_SHYLUCORE_TPETRA
102 #include <Tpetra_CrsMatrix_decl.hpp>
103 #include <Tpetra_CrsMatrix_def.hpp>
104 #endif
105 
106 #include "shylu.h"
107 
108 using namespace std;
109 
110 int main(int argc, char *argv[])
111 {
112 #ifdef HAVE_SHYLUCORE_MPI
113  Teuchos::GlobalMPISession mpiSession(&argc, &argv, 0);
114  Epetra_MpiComm Comm(MPI_COMM_WORLD);
115 #else
116  Epetra_SerialComm Comm;
117 #endif
118  int nProcs, myPID ;
119  Teuchos::ParameterList pLUList ; // ParaLU parameters
120  Teuchos::ParameterList isoList ; // Isorropia parameters
121  string ipFileName = "ShyLU.xml"; // TODO : Accept as i/p
122  int sym = 1; // TODO: Need rectangular factorization for unsymmetric case
123 
124  nProcs = mpiSession.getNProc();
125  myPID = Comm.MyPID();
126 
127  if (myPID == 0)
128  {
129  cout <<"Parallel execution: nProcs="<< nProcs << endl;
130  }
131 
132  // =================== Read input xml file =============================
133  Teuchos::updateParametersFromXmlFile(ipFileName, &pLUList);
134  isoList = pLUList.sublist("Isorropia Input");
135  // Get matrix market file name
136  string MMFileName = Teuchos::getParameter<string>(pLUList, "mm_file");
137 
138  if (myPID == 0)
139  {
140  cout << "Input :" << endl;
141  cout << "ParaLU params " << endl;
142  pLUList.print(std::cout, 2, true, true);
143  cout << "Matrix market file name: " << MMFileName << endl;
144  }
145 
146 
147 
148  // ==================== Read input Matrix ==============================
149  Epetra_CrsMatrix *tempA;
150  EpetraExt::MatrixMarketFileToCrsMatrix(MMFileName.c_str(), Comm, tempA);
151 
152  Epetra_CrsMatrix *A = balanceAndRedistribute(tempA, isoList);
153  delete tempA;
154 
155  Epetra_MultiVector *localS;
156  Epetra_MultiVector *CMV;
157  Epetra_LinearProblem *LP;
158  Amesos_BaseSolver *Solver;
159  int Dnr, Snr;
160  int *DRowElems, *SRowElems, *piv;
161 
162  shylu_factor(A, sym, localS, LP, Solver, Dnr, DRowElems, Snr, SRowElems,
163  piv, CMV);
164 
165  delete localS;
166  delete LP;
167  delete Solver;
168  delete[] DRowElems;
169  delete[] SRowElems;
170  delete[] piv;
171  delete CMV;
172 }
int shylu_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Main function call into ShylU.
Utilities for ShyLU.
Main header file of ShyLU (Include main user calls)