1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * HYPRE_ParCSR_SuperLU interface
11 *
12 *****************************************************************************/
13
14 #include <stdlib.h>
15 #include <stdio.h>
16 #include <math.h>
17
18 #include "utilities/_hypre_utilities.h"
19 #include "HYPRE.h"
20 #include "parcsr_mv/_hypre_parcsr_mv.h"
21 #include "parcsr_ls/HYPRE_parcsr_ls.h"
22
23 #include "HYPRE_FEI.h"
24
25 /*---------------------------------------------------------------------------
26 * SUPERLU include files
27 *-------------------------------------------------------------------------*/
28
29 #ifdef HAVE_SUPERLU_20
30 #include "dsp_defs.h"
31 #include "superlu_util.h"
32
33 typedef struct HYPRE_SuperLU_Struct
34 {
35 int factorized_;
36 int *permR_;
37 int *permC_;
38 SuperMatrix SLU_Lmat;
39 SuperMatrix SLU_Umat;
40 int outputLevel_;
41 }
42 HYPRE_SuperLU;
43 #endif
44
45 #ifdef HAVE_SUPERLU
46 #include "slu_ddefs.h"
47 #include "slu_util.h"
48
49 typedef struct HYPRE_SuperLU_Struct
50 {
51 int factorized_;
52 int *permR_;
53 int *permC_;
54 SuperMatrix SLU_Lmat;
55 SuperMatrix SLU_Umat;
56 int outputLevel_;
57 }
58 HYPRE_SuperLU;
59 #endif
60
61 /***************************************************************************
62 * HYPRE_ParCSR_SuperLUCreate - Return a SuperLU object "solver".
63 *--------------------------------------------------------------------------*/
64
HYPRE_ParCSR_SuperLUCreate(MPI_Comm comm,HYPRE_Solver * solver)65 int HYPRE_ParCSR_SuperLUCreate( MPI_Comm comm, HYPRE_Solver *solver )
66 {
67 #ifdef HAVE_SUPERLU
68 int nprocs;
69 HYPRE_SuperLU *sluPtr;
70
71 MPI_Comm_size(comm, &nprocs);
72 if ( nprocs > 1 )
73 {
74 printf("HYPRE_ParCSR_SuperLUCreate ERROR - too many processors.\n");
75 return -1;
76 }
77 sluPtr = hypre_TAlloc(HYPRE_SuperLU, 1, HYPRE_MEMORY_HOST);
78 hypre_assert ( sluPtr != NULL );
79 sluPtr->factorized_ = 0;
80 sluPtr->permR_ = NULL;
81 sluPtr->permC_ = NULL;
82 sluPtr->outputLevel_ = 0;
83 *solver = (HYPRE_Solver) sluPtr;
84 return 0;
85 #else
86 printf("HYPRE_ParCSR_SuperLUCreate ERROR - SuperLU not enabled.\n");
87 *solver = (HYPRE_Solver) NULL;
88 return -1;
89 #endif
90 }
91
92 /***************************************************************************
93 * HYPRE_ParCSR_SuperLUDestroy - Destroy a SuperLU object.
94 *--------------------------------------------------------------------------*/
95
HYPRE_ParCSR_SuperLUDestroy(HYPRE_Solver solver)96 int HYPRE_ParCSR_SuperLUDestroy( HYPRE_Solver solver )
97 {
98 #ifdef HAVE_SUPERLU
99 HYPRE_SuperLU *sluPtr = (HYPRE_SuperLU *) solver;
100 hypre_assert ( sluPtr != NULL );
101 hypre_TFree(sluPtr->permR_, HYPRE_MEMORY_HOST);
102 hypre_TFree(sluPtr->permC_, HYPRE_MEMORY_HOST);
103 hypre_TFree(sluPtr, HYPRE_MEMORY_HOST);
104 return 0;
105 #else
106 printf("HYPRE_ParCSR_SuperLUDestroy ERROR - SuperLU not enabled.\n");
107 *solver = (HYPRE_Solver) NULL;
108 return -1;
109 #endif
110 }
111
112 /***************************************************************************
113 * HYPRE_ParCSR_SuperLUSetOutputLevel - Set debug level
114 *--------------------------------------------------------------------------*/
115
HYPRE_ParCSR_SuperLUSetOutputLevel(HYPRE_Solver solver,int level)116 int HYPRE_ParCSR_SuperLUSetOutputLevel(HYPRE_Solver solver, int level)
117 {
118 #ifdef HAVE_SUPERLU
119 HYPRE_SuperLU *sluPtr = (HYPRE_SuperLU *) solver;
120 hypre_assert ( sluPtr != NULL );
121 sluPtr->outputLevel_ = level;
122 return 0;
123 #else
124 printf("HYPRE_ParCSR_SuperLUSetOutputLevel ERROR - SuperLU not enabled.\n");
125 *solver = (HYPRE_Solver) NULL;
126 return -1;
127 #endif
128 }
129
130 /***************************************************************************
131 * HYPRE_ParCSR_SuperLUSetup - Set up function for SuperLU.
132 *--------------------------------------------------------------------------*/
133
HYPRE_ParCSR_SuperLUSetup(HYPRE_Solver solver,HYPRE_ParCSRMatrix A_csr,HYPRE_ParVector b,HYPRE_ParVector x)134 int HYPRE_ParCSR_SuperLUSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix A_csr,
135 HYPRE_ParVector b, HYPRE_ParVector x )
136 {
137 #ifdef HAVE_SUPERLU
138 int startRow, endRow, nrows, *partition, *AdiagI, *AdiagJ, nnz;
139 int irow, colNum, index, *cscI, *cscJ, jcol, *colLengs;
140 int *etree, permcSpec, lwork, panelSize, relax, info;
141 double *AdiagA, *cscA, diagPivotThresh, dropTol;
142 char refact[1];
143 hypre_CSRMatrix *Adiag;
144 HYPRE_SuperLU *sluPtr;
145 SuperMatrix sluAmat, auxAmat;
146 superlu_options_t slu_options;
147 SuperLUStat_t slu_stat;
148
149 /* ---------------------------------------------------------------- */
150 /* get matrix information */
151 /* ---------------------------------------------------------------- */
152
153 sluPtr = (HYPRE_SuperLU *) solver;
154 hypre_assert ( sluPtr != NULL );
155 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );
156 startRow = partition[0];
157 endRow = partition[1] - 1;
158 nrows = endRow - startRow + 1;
159 hypre_TFree(partition, HYPRE_MEMORY_HOST);
160 if ( startRow != 0 )
161 {
162 printf("HYPRE_ParCSR_SuperLUSetup ERROR - start row != 0.\n");
163 return -1;
164 }
165
166 /* ---------------------------------------------------------------- */
167 /* get hypre matrix */
168 /* ---------------------------------------------------------------- */
169
170 Adiag = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *) A_csr);
171 AdiagI = hypre_CSRMatrixI(Adiag);
172 AdiagJ = hypre_CSRMatrixJ(Adiag);
173 AdiagA = hypre_CSRMatrixData(Adiag);
174 nnz = AdiagI[nrows];
175
176 /* ---------------------------------------------------------------- */
177 /* convert the csr matrix into csc matrix */
178 /* ---------------------------------------------------------------- */
179
180 colLengs = hypre_TAlloc(int, nrows , HYPRE_MEMORY_HOST);
181 for ( irow = 0; irow < nrows; irow++ ) colLengs[irow] = 0;
182 for ( irow = 0; irow < nrows; irow++ )
183 for ( jcol = AdiagI[irow]; jcol < AdiagI[irow+1]; jcol++ )
184 colLengs[AdiagJ[jcol]]++;
185 cscJ = hypre_TAlloc(int, (nrows+1) , HYPRE_MEMORY_HOST);
186 cscI = hypre_TAlloc(int, nnz , HYPRE_MEMORY_HOST);
187 cscA = hypre_TAlloc(double, nnz , HYPRE_MEMORY_HOST);
188 cscJ[0] = 0;
189 nnz = 0;
190 for ( jcol = 1; jcol <= nrows; jcol++ )
191 {
192 nnz += colLengs[jcol-1];
193 cscJ[jcol] = nnz;
194 }
195 for ( irow = 0; irow < nrows; irow++ )
196 {
197 for ( jcol = AdiagI[irow]; jcol < AdiagI[irow+1]; jcol++ )
198 {
199 colNum = AdiagJ[jcol];
200 index = cscJ[colNum]++;
201 cscI[index] = irow;
202 cscA[index] = AdiagA[jcol];
203 }
204 }
205 cscJ[0] = 0;
206 nnz = 0;
207 for ( jcol = 1; jcol <= nrows; jcol++ )
208 {
209 nnz += colLengs[jcol-1];
210 cscJ[jcol] = nnz;
211 }
212 hypre_TFree(colLengs, HYPRE_MEMORY_HOST);
213
214 /* ---------------------------------------------------------------- */
215 /* create SuperMatrix */
216 /* ---------------------------------------------------------------- */
217
218 dCreate_CompCol_Matrix(&sluAmat,nrows,nrows,cscJ[nrows],cscA,cscI,
219 cscJ, SLU_NC, SLU_D, SLU_GE);
220 etree = hypre_TAlloc(int, nrows , HYPRE_MEMORY_HOST);
221 sluPtr->permC_ = hypre_TAlloc(int, nrows , HYPRE_MEMORY_HOST);
222 sluPtr->permR_ = hypre_TAlloc(int, nrows , HYPRE_MEMORY_HOST);
223 permcSpec = 0;
224 get_perm_c(permcSpec, &sluAmat, sluPtr->permC_);
225 slu_options.Fact = DOFACT;
226 slu_options.SymmetricMode = NO;
227 sp_preorder(&slu_options, &sluAmat, sluPtr->permC_, etree, &auxAmat);
228 diagPivotThresh = 1.0;
229 dropTol = 0.0;
230 panelSize = sp_ienv(1);
231 relax = sp_ienv(2);
232 StatInit(&slu_stat);
233 lwork = 0;
234 slu_options.ColPerm = MY_PERMC;
235 slu_options.DiagPivotThresh = diagPivotThresh;
236
237 dgstrf(&slu_options, &auxAmat, dropTol, relax, panelSize,
238 etree, NULL, lwork, sluPtr->permC_, sluPtr->permR_,
239 &(sluPtr->SLU_Lmat), &(sluPtr->SLU_Umat), &slu_stat, &info);
240 Destroy_CompCol_Permuted(&auxAmat);
241 Destroy_CompCol_Matrix(&sluAmat);
242 hypre_TFree(etree, HYPRE_MEMORY_HOST);
243 sluPtr->factorized_ = 1;
244 StatFree(&slu_stat);
245 return 0;
246 #else
247 printf("HYPRE_ParCSR_SuperLUSetup ERROR - SuperLU not enabled.\n");
248 *solver = (HYPRE_Solver) NULL;
249 return -1;
250 #endif
251 }
252
253 /***************************************************************************
254 * HYPRE_ParCSR_SuperLUSolve - Solve function for SuperLU.
255 *--------------------------------------------------------------------------*/
256
HYPRE_ParCSR_SuperLUSolve(HYPRE_Solver solver,HYPRE_ParCSRMatrix A,HYPRE_ParVector b,HYPRE_ParVector x)257 int HYPRE_ParCSR_SuperLUSolve(HYPRE_Solver solver, HYPRE_ParCSRMatrix A,
258 HYPRE_ParVector b, HYPRE_ParVector x )
259 {
260 #ifdef HAVE_SUPERLU
261 int nrows, i, info;
262 double *bData, *xData;
263 SuperMatrix B;
264 SuperLUStat_t slu_stat;
265 trans_t trans;
266 HYPRE_SuperLU *sluPtr = (HYPRE_SuperLU *) solver;
267
268 /* ---------------------------------------------------------------- */
269 /* make sure setup has been called */
270 /* ---------------------------------------------------------------- */
271
272 hypre_assert ( sluPtr != NULL );
273 if ( ! (sluPtr->factorized_) )
274 {
275 printf("HYPRE_ParCSR_SuperLUSolve ERROR - not factorized yet.\n");
276 return -1;
277 }
278
279 /* ---------------------------------------------------------------- */
280 /* fetch right hand side and solution vector */
281 /* ---------------------------------------------------------------- */
282
283 xData = hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *)x));
284 bData = hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *)b));
285 nrows = hypre_ParVectorGlobalSize((hypre_ParVector *)x);
286 for (i = 0; i < nrows; i++) xData[i] = bData[i];
287
288 /* ---------------------------------------------------------------- */
289 /* solve */
290 /* ---------------------------------------------------------------- */
291
292 dCreate_Dense_Matrix(&B, nrows, 1, bData, nrows, SLU_DN, SLU_D,SLU_GE);
293
294 /* -------------------------------------------------------------
295 * solve the problem
296 * -----------------------------------------------------------*/
297
298 trans = NOTRANS;
299 StatInit(&slu_stat);
300 dgstrs (trans, &(sluPtr->SLU_Lmat), &(sluPtr->SLU_Umat),
301 sluPtr->permC_, sluPtr->permR_, &B, &slu_stat, &info);
302 Destroy_SuperMatrix_Store(&B);
303 StatFree(&slu_stat);
304 return 0;
305 #else
306 printf("HYPRE_ParCSR_SuperLUSolve ERROR - SuperLU not enabled.\n");
307 *solver = (HYPRE_Solver) NULL;
308 return -1;
309 #endif
310 }
311
312