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 /* HYPRE_LSI_MLI interface */
10 /*--------------------------------------------------------------------------*/
11 /* local functions :
12 *
13 * HYPRE_LSI_MLICreate
14 * HYPRE_LSI_MLIDestroy
15 * HYPRE_LSI_MLISetup
16 * HYPRE_LSI_MLISolve
17 * HYPRE_LSI_MLISetParams
18 *--------------------------------------------------------------------------
19 * HYPRE_LSI_MLICreateNodeEqnMap
20 * HYPRE_LSI_MLIAdjustNodeEqnMap
21 * HYPRE_LSI_MLIAdjustNullSpace
22 * HYPRE_LSI_MLISetFEData
23 * HYPRE_LSI_MLILoadNodalCoordinates
24 * HYPRE_LSI_MLILoadMatrixScalings
25 * HYPRE_LSI_MLILoadMaterialLabels
26 *--------------------------------------------------------------------------
27 * HYPRE_LSI_MLIFEDataCreate
28 * HYPRE_LSI_MLIFEDataDestroy
29 * HYPRE_LSI_MLIFEDataInitFields
30 * HYPRE_LSI_MLIFEDataInitElemBlock
31 * HYPRE_LSI_MLIFEDataInitElemNodeList
32 * HYPRE_LSI_MLIFEDataInitSharedNodes
33 * HYPRE_LSI_MLIFEDataInitComplete
34 * HYPRE_LSI_MLIFEDataLoadElemMatrix
35 * HYPRE_LSI_MLIFEDataWriteToFile
36 *--------------------------------------------------------------------------
37 * HYPRE_LSI_MLISFEICreate
38 * HYPRE_LSI_MLISFEIDestroy
39 * HYPRE_LSI_MLISFEILoadElemMatrices
40 * HYPRE_LSI_MLISFEIAddNumElems
41 ****************************************************************************/
42
43 /****************************************************************************/
44 /* system include files */
45 /*--------------------------------------------------------------------------*/
46
47 #include "HYPRE_utilities.h"
48 #include <string.h>
49 #include <stdlib.h>
50 #include <stdio.h>
51 #include <math.h>
52
53 #if 0 /* RDF: Not sure this is really needed */
54 #ifdef WIN32
55 #define strcmp _stricmp
56 #endif
57 #endif
58
59 /****************************************************************************/
60 /* MLI include files */
61 /*--------------------------------------------------------------------------*/
62
63 #ifdef HAVE_MLI
64 #include "mli.h"
65 #include "mli_utils.h"
66 #include "mli_method.h"
67 #endif
68 #include "HYPRE_LSI_mli.h"
69
70 /****************************************************************************/
71 /* HYPRE_LSI_MLI data structure */
72 /*--------------------------------------------------------------------------*/
73
74 typedef struct HYPRE_LSI_MLI_Struct
75 {
76 #ifdef HAVE_MLI
77 MLI *mli_;
78 MLI_FEData *feData_; /* holds FE information */
79 MLI_SFEI *sfei_; /* holds FE information */
80 MLI_Mapper *mapper_; /* holds mapping information */
81 #endif
82 MPI_Comm mpiComm_;
83 int outputLevel_; /* for diagnostics */
84 int nLevels_; /* max number of levels */
85 int cycleType_; /* 1 for V and 2 for W */
86 int maxIterations_; /* default - 1 iteration */
87 char method_[20]; /* default - smoothed aggregation */
88 char coarsenScheme_[20]; /* default - default in MLI */
89 char preSmoother_[20]; /* default - symmetric Gauss Seidel */
90 char postSmoother_ [20]; /* default - symmetric Gauss Seidel */
91 int preNSweeps_; /* default - 2 smoothing steps */
92 int postNSweeps_; /* default - 2 smoothing steps */
93 double *preSmootherWts_; /* relaxation weights */
94 double *postSmootherWts_; /* relaxation weights */
95 int smootherPrintRNorm_; /* for smoother diagnostics */
96 int smootherFindOmega_; /* for SGS smoother */
97 double strengthThreshold_; /* strength threshold */
98 char coarseSolver_[20]; /* default = SuperLU */
99 int coarseSolverNSweeps_; /* number of sweeps (if iterative used) */
100 double *coarseSolverWts_; /* relaxation weight (if Jacobi used) */
101 int minCoarseSize_; /* minimum coarse grid size */
102 int scalar_; /* aggregate in a scalar manner */
103 int nodeDOF_; /* node degree of freedom */
104 int spaceDim_; /* 2D or 3D */
105 int nSpaceDim_; /* number of null vectors */
106 int localNEqns_; /* number of equations locally */
107 int nCoordAccept_; /* flag to accept nodal coordinate or not */
108 double *nCoordinates_; /* for storing nodal coordinates */
109 double *nullScales_; /* scaling vector for null space */
110 int calibrationSize_; /* for calibration smoothed aggregation */
111 double Pweight_;
112 int SPLevel_;
113 char paramFile_[50];
114 int adjustNullSpace_;
115 int numResetNull_;
116 int *resetNullIndices_;
117 int numMatLabels_; /* for controlling aggregation */
118 int *matLabels_;
119 int printNullSpace_;
120 int symmetric_;
121 int injectionForR_;
122 HYPRE_ParCSRMatrix correctionMatrix_; /* for nullspace correction */
123 int numSmoothVecs_;
124 int smoothVecSteps_;
125 double arpackTol_;
126 }
127 HYPRE_LSI_MLI;
128
129 /****************************************************************************/
130 /* HYPRE_MLI_FEData data structure */
131 /*--------------------------------------------------------------------------*/
132
133 typedef struct HYPRE_MLI_FEData_Struct
134 {
135 #ifdef HAVE_MLI
136 MPI_Comm comm_; /* MPI communicator */
137 MLI_FEData *fedata_; /* holds FE information */
138 int fedataOwn_; /* flag to indicate ownership */
139 int computeNull_; /* flag - compute null space or not */
140 int nullDim_; /* number of null space vectors */
141 #endif
142 }
143 HYPRE_MLI_FEData;
144
145 /****************************************************************************/
146 /* HYPRE_MLI_SFEI data structure */
147 /*--------------------------------------------------------------------------*/
148
149 typedef struct HYPRE_MLI_SFEI_Struct
150 {
151 #ifdef HAVE_MLI
152 MPI_Comm comm_; /* MPI communicator */
153 MLI_SFEI *sfei_; /* holds FE information */
154 int sfeiOwn_; /* flag to indicate ownership */
155 #endif
156 }
157 HYPRE_MLI_SFEI;
158
159 /****************************************************************************/
160 /* HYPRE_LSI_MLICreate */
161 /*--------------------------------------------------------------------------*/
162
163 extern "C"
HYPRE_LSI_MLICreate(MPI_Comm comm,HYPRE_Solver * solver)164 int HYPRE_LSI_MLICreate( MPI_Comm comm, HYPRE_Solver *solver )
165 {
166 HYPRE_LSI_MLI *mli_object = hypre_TAlloc(HYPRE_LSI_MLI, 1, HYPRE_MEMORY_HOST);
167 *solver = (HYPRE_Solver) mli_object;
168 mli_object->mpiComm_ = comm;
169 mli_object->outputLevel_ = 0;
170 mli_object->nLevels_ = 0;
171 mli_object->maxIterations_ = 1;
172 mli_object->cycleType_ = 1;
173 strcpy(mli_object->method_ , "AMGSA");
174 strcpy(mli_object->coarsenScheme_ , "default"); /* default in MLI */
175 strcpy(mli_object->preSmoother_, "default");
176 strcpy(mli_object->postSmoother_, "default");
177 mli_object->preNSweeps_ = 1;
178 mli_object->postNSweeps_ = 1;
179 mli_object->preSmootherWts_ = NULL;
180 mli_object->postSmootherWts_ = NULL;
181 mli_object->smootherPrintRNorm_ = 0;
182 mli_object->smootherFindOmega_ = 0;
183 mli_object->strengthThreshold_ = 0.0;
184 strcpy(mli_object->coarseSolver_, "default"); /* default in MLI */
185 mli_object->coarseSolverNSweeps_ = 0;
186 mli_object->coarseSolverWts_ = NULL;
187 mli_object->minCoarseSize_ = 0;
188 mli_object->scalar_ = 0;
189 mli_object->nodeDOF_ = 1;
190 mli_object->spaceDim_ = 1;
191 mli_object->nSpaceDim_ = 1;
192 mli_object->localNEqns_ = 0;
193 mli_object->nCoordinates_ = NULL;
194 mli_object->nCoordAccept_ = 0;
195 mli_object->nullScales_ = NULL;
196 mli_object->calibrationSize_ = 0;
197 mli_object->Pweight_ = -1.0; /* default in MLI */
198 mli_object->SPLevel_ = 0; /* default in MLI */
199 mli_object->adjustNullSpace_ = 0;
200 mli_object->numResetNull_ = 0;
201 mli_object->resetNullIndices_ = NULL;
202 mli_object->correctionMatrix_ = NULL;
203 strcpy(mli_object->paramFile_, "empty");
204 mli_object->numMatLabels_ = 0;
205 mli_object->matLabels_ = NULL;
206 mli_object->printNullSpace_ = 0;
207 mli_object->symmetric_ = 1;
208 mli_object->injectionForR_ = 0;
209 mli_object->numSmoothVecs_ = 0;
210 mli_object->smoothVecSteps_ = 0;
211 mli_object->arpackTol_ = 0.0;
212 #ifdef HAVE_MLI
213 mli_object->mli_ = NULL;
214 mli_object->sfei_ = NULL;
215 mli_object->feData_ = NULL;
216 mli_object->mapper_ = NULL;
217 return 0;
218 #else
219 printf("MLI not available.\n");
220 exit(1);
221 return -1;
222 #endif
223 }
224
225 /****************************************************************************/
226 /* HYPRE_LSI_MLIDestroy */
227 /*--------------------------------------------------------------------------*/
228
229 extern "C"
HYPRE_LSI_MLIDestroy(HYPRE_Solver solver)230 int HYPRE_LSI_MLIDestroy( HYPRE_Solver solver )
231 {
232 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
233
234 if ( mli_object->preSmootherWts_ != NULL )
235 delete [] mli_object->preSmootherWts_;
236 if ( mli_object->postSmootherWts_ != NULL )
237 delete [] mli_object->postSmootherWts_;
238 if ( mli_object->coarseSolverWts_ != NULL )
239 delete [] mli_object->coarseSolverWts_;
240 if ( mli_object->nCoordinates_ != NULL )
241 delete [] mli_object->nCoordinates_;
242 if ( mli_object->nullScales_ != NULL )
243 delete [] mli_object->nullScales_;
244 if ( mli_object->resetNullIndices_ != NULL )
245 delete [] mli_object->resetNullIndices_;
246 if ( mli_object->correctionMatrix_ != NULL )
247 HYPRE_ParCSRMatrixDestroy(mli_object->correctionMatrix_);
248 if ( mli_object->matLabels_ != NULL ) delete [] mli_object->matLabels_;
249 #ifdef HAVE_MLI
250 if ( mli_object->feData_ != NULL ) delete mli_object->feData_;
251 if ( mli_object->mli_ != NULL ) delete mli_object->mli_;
252 free( mli_object );
253 return 0;
254 #else
255 free( mli_object );
256 printf("MLI not available.\n");
257 return -1;
258 #endif
259 }
260
261 /****************************************************************************/
262 /* HYPRE_LSI_MLISetup */
263 /*--------------------------------------------------------------------------*/
264
265 extern "C"
HYPRE_LSI_MLISetup(HYPRE_Solver solver,HYPRE_ParCSRMatrix A,HYPRE_ParVector b,HYPRE_ParVector x)266 int HYPRE_LSI_MLISetup( HYPRE_Solver solver, HYPRE_ParCSRMatrix A,
267 HYPRE_ParVector b, HYPRE_ParVector x )
268 {
269 #ifdef HAVE_MLI
270 int targc, nNodes, iZero=0;
271 double tol=1.0e-8;
272 char *targv[6], paramString[100];;
273 HYPRE_LSI_MLI *mli_object;
274 MLI_Matrix *mli_mat;
275 MLI_Method *method;
276 MPI_Comm mpiComm;
277 MLI *mli;
278
279 /* -------------------------------------------------------- */
280 /* create object */
281 /* -------------------------------------------------------- */
282
283 mli_object = (HYPRE_LSI_MLI *) solver;
284 mpiComm = mli_object->mpiComm_;
285 mli = new MLI( mpiComm );
286 if ( mli_object->mli_ != NULL ) delete mli_object->mli_;
287 mli_object->mli_ = mli;
288
289 /* -------------------------------------------------------- */
290 /* set general parameters */
291 /* -------------------------------------------------------- */
292
293 if (!strcmp(mli_object->method_,"AMGSADD") ||
294 !strcmp(mli_object->method_,"AMGSADDe")) mli_object->nLevels_ = 2;
295
296 mli->setNumLevels( mli_object->nLevels_ );
297 mli->setTolerance( tol );
298
299 /* -------------------------------------------------------- */
300 /* set method specific parameters */
301 /* -------------------------------------------------------- */
302
303 method = MLI_Method_CreateFromName(mli_object->method_,mpiComm);
304 if ( mli_object->outputLevel_ > 0 )
305 {
306 sprintf(paramString, "setOutputLevel %d", mli_object->outputLevel_);
307 method->setParams( paramString, 0, NULL );
308 }
309 if ( mli_object->nLevels_ > 0 )
310 {
311 sprintf(paramString, "setNumLevels %d", mli_object->nLevels_);
312 method->setParams( paramString, 0, NULL );
313 }
314 if ( mli_object->strengthThreshold_ > 0.0 )
315 {
316 sprintf(paramString, "setStrengthThreshold %f",
317 mli_object->strengthThreshold_);
318 method->setParams( paramString, 0, NULL );
319 }
320 if ( mli_object->scalar_ == 1 )
321 {
322 strcpy(paramString, "scalar");
323 method->setParams( paramString, 0, NULL );
324 }
325 if ( mli_object->symmetric_ == 0 )
326 {
327 strcpy(paramString, "nonsymmetric");
328 method->setParams( paramString, 0, NULL );
329 }
330 if ( mli_object->injectionForR_ == 1 )
331 {
332 strcpy(paramString, "useInjectionForR");
333 method->setParams( paramString, 0, NULL );
334 }
335 if ( mli_object->smootherPrintRNorm_ == 1 )
336 {
337 strcpy(paramString, "setSmootherPrintRNorm");
338 method->setParams( paramString, 0, NULL );
339 }
340 if ( mli_object->smootherFindOmega_ == 1 )
341 {
342 strcpy(paramString, "setSmootherFindOmega");
343 method->setParams( paramString, 0, NULL );
344 }
345 if ( mli_object->numSmoothVecs_ > 0 )
346 {
347 sprintf(paramString, "setSmoothVec %d", mli_object->numSmoothVecs_);
348 method->setParams( paramString, 0, NULL );
349 if ( mli_object->smoothVecSteps_ > 0 )
350 sprintf(paramString, "setSmoothVecSteps %d",
351 mli_object->smoothVecSteps_);
352 else
353 sprintf(paramString, "setSmoothVecSteps 5");
354 method->setParams( paramString, 0, NULL );
355 }
356 if ( mli_object->arpackTol_ > 0.0 )
357 {
358 sprintf(paramString, "arpackTol %e", mli_object->arpackTol_);
359 method->setParams( paramString, 0, NULL );
360 }
361
362 /* -------------------------------------------------------- */
363 /* set up presmoother */
364 /* -------------------------------------------------------- */
365
366 if ( strcmp(mli_object->preSmoother_, "default") )
367 {
368 targc = 2;
369 targv[0] = (char *) &mli_object->preNSweeps_;
370 targv[1] = (char *) mli_object->preSmootherWts_;
371 sprintf(paramString, "setPreSmoother %s", mli_object->preSmoother_);
372 method->setParams( paramString, targc, targv );
373 }
374
375 /* -------------------------------------------------------- */
376 /* set up postsmoother */
377 /* -------------------------------------------------------- */
378
379 if ( strcmp(mli_object->preSmoother_, "default") )
380 {
381 targc = 2;
382 targv[0] = (char *) &mli_object->postNSweeps_;
383 targv[1] = (char *) mli_object->postSmootherWts_;
384 sprintf(paramString, "setPostSmoother %s", mli_object->postSmoother_);
385 method->setParams( paramString, targc, targv );
386 }
387
388 /* -------------------------------------------------------- */
389 /* set up coarse solver */
390 /* -------------------------------------------------------- */
391
392 if ( strcmp(mli_object->coarseSolver_, "default") )
393 {
394 targc = 2;
395 targv[0] = (char *) &(mli_object->coarseSolverNSweeps_);
396 targv[1] = (char *) mli_object->coarseSolverWts_;
397 sprintf(paramString, "setCoarseSolver %s", mli_object->coarseSolver_);
398 method->setParams( paramString, targc, targv );
399 }
400
401 /* -------------------------------------------------------- */
402 /* load minimum coarse grid size and others */
403 /* -------------------------------------------------------- */
404
405 if ( mli_object->minCoarseSize_ != 0 )
406 {
407 sprintf(paramString, "setMinCoarseSize %d",mli_object->minCoarseSize_);
408 method->setParams( paramString, 0, NULL );
409 }
410 if ( mli_object->Pweight_ >= 0.0 )
411 {
412 sprintf( paramString, "setPweight %e", mli_object->Pweight_ );
413 method->setParams( paramString, 0, NULL );
414 if ( mli_object->SPLevel_ > 0 )
415 {
416 sprintf( paramString, "setSPLevel %d", mli_object->SPLevel_ );
417 method->setParams( paramString, 0, NULL );
418 }
419 }
420 if ( strcmp(mli_object->coarsenScheme_, "default") )
421 {
422 sprintf(paramString, "setCoarsenScheme %s",mli_object->coarsenScheme_);
423 method->setParams( paramString, 0, NULL );
424 }
425
426 /* -------------------------------------------------------- */
427 /* load calibration size */
428 /* -------------------------------------------------------- */
429
430 if ( mli_object->calibrationSize_ > 0 )
431 {
432 sprintf(paramString,"setCalibrationSize %d",mli_object->calibrationSize_);
433 method->setParams( paramString, 0, NULL );
434 }
435
436 /* -------------------------------------------------------- */
437 /* load FEData, if there is any */
438 /* -------------------------------------------------------- */
439
440 if ( mli_object->feData_ != NULL )
441 mli->setFEData( 0, mli_object->feData_, mli_object->mapper_ );
442 if ( mli_object->sfei_ != NULL ) mli->setSFEI(0, mli_object->sfei_);
443 //mli_object->mapper_ = NULL;
444 //mli_object->feData_ = NULL;
445
446 /* -------------------------------------------------------- */
447 /* load null space, if there is any */
448 /* -------------------------------------------------------- */
449
450 if ((mli_object->printNullSpace_ & 1) == 1 )
451 {
452 strcpy( paramString, "printNullSpace" );
453 method->setParams( paramString, 0, NULL );
454 }
455 if ((mli_object->printNullSpace_ & 2) == 2 )
456 {
457 strcpy( paramString, "printElemNodeList" );
458 method->setParams( paramString, 0, NULL );
459 }
460 if ((mli_object->printNullSpace_ & 4) == 4 )
461 {
462 strcpy( paramString, "printNodalCoord" );
463 method->setParams( paramString, 0, NULL );
464 }
465 if ( mli_object->nCoordinates_ != NULL )
466 {
467 nNodes = mli_object->localNEqns_ / mli_object->nodeDOF_;
468 targc = 6;
469 targv[0] = (char *) &nNodes;
470 targv[1] = (char *) &(mli_object->nodeDOF_);
471 targv[2] = (char *) &(mli_object->spaceDim_);
472 targv[3] = (char *) mli_object->nCoordinates_;
473 targv[4] = (char *) &(mli_object->nSpaceDim_);
474 targv[5] = (char *) mli_object->nullScales_;
475 strcpy( paramString, "setNodalCoord" );
476 method->setParams( paramString, targc, targv );
477 #if 0
478 if ( mli_object->adjustNullSpace_ == 1 &&
479 mli_object->correctionMatrix_ != NULL )
480 {
481 int iV, irow, numNS, length, *partition, mypid;
482 double *NSpace, *vecInData, *vecOutData, *nullCorrection;
483 HYPRE_IJVector IJVecIn, IJVecOut;
484 HYPRE_ParCSRMatrix hypreA;
485 hypre_ParVector *hypreVecIn, *hypreVecOut;
486 HYPRE_Solver solver;
487
488 strcpy( paramString, "getNullSpace" );
489 method->getParams( paramString, &targc, targv );
490 numNS = *(int *) targv[1];
491 NSpace = (double *) targv[2];
492 length = *(int *) targv[3];
493 hypreA = mli_object->correctionMatrix_;
494 HYPRE_ParCSRMatrixGetRowPartitioning( hypreA, &partition );
495 MPI_Comm_rank( mpiComm , &mypid );
496 HYPRE_IJVectorCreate(mpiComm, partition[mypid],
497 partition[mypid+1]-1, &IJVecIn);
498 HYPRE_IJVectorSetObjectType(IJVecIn, HYPRE_PARCSR);
499 HYPRE_IJVectorInitialize(IJVecIn);
500 HYPRE_IJVectorAssemble(IJVecIn);
501 HYPRE_IJVectorCreate(mpiComm, partition[mypid],
502 partition[mypid+1]-1, &IJVecOut);
503 HYPRE_IJVectorSetObjectType(IJVecOut, HYPRE_PARCSR);
504 HYPRE_IJVectorInitialize(IJVecOut);
505 HYPRE_IJVectorAssemble(IJVecOut);
506 HYPRE_IJVectorGetObject(IJVecIn, (void **) &hypreVecIn);
507 HYPRE_IJVectorGetObject(IJVecOut, (void **) &hypreVecOut);
508 vecInData = hypre_VectorData(hypre_ParVectorLocalVector(hypreVecIn));
509 vecOutData = hypre_VectorData(hypre_ParVectorLocalVector(hypreVecOut));
510 nullCorrection = new double[numNS*length];
511 HYPRE_ParCSRGMRESCreate(mpiComm, &solver);
512 HYPRE_ParCSRGMRESSetKDim(solver, 100);
513 HYPRE_ParCSRGMRESSetMaxIter(solver, 100);
514 HYPRE_ParCSRGMRESSetTol(solver, 1.0e-8);
515 HYPRE_ParCSRGMRESSetPrecond(solver, HYPRE_ParCSRDiagScale,
516 HYPRE_ParCSRDiagScaleSetup,NULL);
517 HYPRE_ParCSRGMRESSetLogging(solver, mli_object->outputLevel_);
518 HYPRE_ParCSRGMRESSetup(solver, hypreA, (HYPRE_ParVector) hypreVecIn,
519 (HYPRE_ParVector) hypreVecOut);
520
521 for ( iV = 0; iV < numNS; iV++ )
522 {
523 for ( irow = 0; irow < length; irow++ )
524 vecOutData[irow] = NSpace[length*iV+irow];
525 HYPRE_ParCSRMatrixMatvec(1.0,hypreA,(HYPRE_ParVector) hypreVecOut,
526 0.0, (HYPRE_ParVector) hypreVecIn);
527 for ( irow = 0; irow < length; irow++ ) vecOutData[irow] = 0.0;
528 HYPRE_ParCSRGMRESSolve(solver,A,(HYPRE_ParVector) hypreVecIn,
529 (HYPRE_ParVector) hypreVecOut);
530 for ( irow = 0; irow < length; irow++ )
531 nullCorrection[length*iV+irow] = - vecOutData[irow];
532 }
533 strcpy( paramString, "adjustNullSpace" );
534 targc = 1;
535 targv[0] = (char *) nullCorrection;
536 method->setParams( paramString, targc, targv );
537 HYPRE_ParCSRGMRESDestroy(solver);
538 HYPRE_IJVectorDestroy( IJVecIn );
539 HYPRE_IJVectorDestroy( IJVecOut );
540 delete [] nullCorrection;
541 free( partition );
542 }
543 if ( mli_object->numResetNull_ != 0 &&
544 mli_object->resetNullIndices_ != NULL )
545 {
546 int *rowPartition, my_id;
547 HYPRE_ParCSRMatrixGetRowPartitioning( A, &rowPartition );
548 MPI_Comm_rank( mpiComm , &my_id );
549 targc = 3;
550 targv[0] = (char *) &(mli_object->numResetNull_);
551 targv[1] = (char *) &(rowPartition[my_id]);
552 targv[2] = (char *) (mli_object->resetNullIndices_);
553 strcpy( paramString, "resetNullSpaceComponents" );
554 method->setParams( paramString, targc, targv );
555 free( rowPartition );
556 }
557 #endif
558 }
559 else
560 {
561 targc = 4;
562 targv[0] = (char *) &(mli_object->nodeDOF_);
563 //if ( mli_object->nSpaceDim_ > mli_object->nodeDOF_ )
564 // mli_object->nSpaceDim_ = mli_object->nodeDOF_;
565 targv[1] = (char *) &(mli_object->nSpaceDim_);
566 targv[2] = (char *) NULL;
567 targv[3] = (char *) &iZero;
568 strcpy( paramString, "setNullSpace" );
569 method->setParams( paramString, targc, targv );
570 }
571 if ( mli_object->correctionMatrix_ != NULL )
572 {
573 HYPRE_ParCSRMatrixDestroy( mli_object->correctionMatrix_ );
574 mli_object->correctionMatrix_ = NULL;
575 }
576 if (!strcmp(mli_object->method_,"AMGRS") )
577 {
578 sprintf( paramString, "setNodeDOF %d", mli_object->nodeDOF_ );
579 method->setParams( paramString, 0, NULL );
580 }
581
582 /* -------------------------------------------------------- */
583 /* load material labels, if there is any */
584 /* -------------------------------------------------------- */
585
586 if ( mli_object->matLabels_ != NULL )
587 {
588 strcpy( paramString, "setLabels" );
589 targc = 3;
590 targv[0] = (char *) &(mli_object->numMatLabels_);
591 targv[1] = (char *) &iZero;
592 targv[2] = (char *) mli_object->matLabels_;
593 method->setParams( paramString, targc, targv );
594 }
595
596 /* -------------------------------------------------------- */
597 /* set parameter file */
598 /* -------------------------------------------------------- */
599
600 if ( strcmp(mli_object->paramFile_, "empty") )
601 {
602 targc = 1;
603 targv[0] = (char *) mli_object->paramFile_;
604 strcpy( paramString, "setParamFile" );
605 method->setParams( paramString, targc, targv );
606 }
607 if ( mli_object->outputLevel_ >= 1 )
608 {
609 strcpy( paramString, "print" );
610 method->setParams( paramString, 0, NULL );
611 }
612
613 /* -------------------------------------------------------- */
614 /* finally, set up */
615 /* -------------------------------------------------------- */
616
617 strcpy( paramString, "HYPRE_ParCSR" );
618 mli_mat = new MLI_Matrix((void*) A, paramString, NULL);
619 mli->setMethod( method );
620 mli->setSystemMatrix( 0, mli_mat );
621 mli->setOutputLevel( mli_object->outputLevel_ );
622 mli->setup();
623 mli->setMaxIterations( mli_object->maxIterations_ );
624 mli->setCyclesAtLevel( -1, mli_object->cycleType_ );
625 return 0;
626 #else
627 printf("MLI not available.\n");
628 return -1;
629 #endif
630 }
631
632 /****************************************************************************/
633 /* HYPRE_LSI_MLISolve */
634 /*--------------------------------------------------------------------------*/
635
636 extern "C"
HYPRE_LSI_MLISolve(HYPRE_Solver solver,HYPRE_ParCSRMatrix A,HYPRE_ParVector b,HYPRE_ParVector x)637 int HYPRE_LSI_MLISolve( HYPRE_Solver solver, HYPRE_ParCSRMatrix A,
638 HYPRE_ParVector b, HYPRE_ParVector x )
639 {
640 #ifdef HAVE_MLI
641 HYPRE_LSI_MLI *mli_object;
642 MLI_Vector *sol, *rhs;
643 char paramString[100];
644
645 strcpy(paramString, "HYPRE_ParVector");
646 sol = new MLI_Vector( (void *) x, paramString, NULL);
647 rhs = new MLI_Vector( (void *) b, paramString, NULL);
648
649 mli_object = (HYPRE_LSI_MLI *) solver;
650 if ( mli_object->mli_ == NULL )
651 {
652 printf("HYPRE_LSI_MLISolve ERROR : mli not instantiated.\n");
653 exit(1);
654 }
655 mli_object->mli_->solve( sol, rhs);
656
657 return 0;
658 #else
659 printf("MLI not available.\n");
660 return -1;
661 #endif
662 }
663
664 /****************************************************************************/
665 /* HYPRE_LSI_MLISetParams */
666 /*--------------------------------------------------------------------------*/
667
668 extern "C"
HYPRE_LSI_MLISetParams(HYPRE_Solver solver,char * paramString)669 int HYPRE_LSI_MLISetParams( HYPRE_Solver solver, char *paramString )
670 {
671 int i, mypid;
672 double weight;
673 HYPRE_LSI_MLI *mli_object;
674 char param1[256], param2[256], param3[256];
675
676 mli_object = (HYPRE_LSI_MLI *) solver;
677 sscanf(paramString,"%s", param1);
678 if ( strcmp(param1, "MLI") )
679 {
680 printf("HYPRE_LSI_MLI::parameters not for me.\n");
681 return 1;
682 }
683 MPI_Comm_rank( mli_object->mpiComm_, &mypid );
684 sscanf(paramString,"%s %s", param1, param2);
685 if ( !strcmp(param2, "help") )
686 {
687 if ( mypid == 0 )
688 {
689 printf("%4d : Available options for MLI are : \n", mypid);
690 printf("\t outputLevel <d> \n");
691 printf("\t numLevels <d> \n");
692 printf("\t maxIterations <d> \n");
693 printf("\t cycleType <'V','W'> \n");
694 printf("\t strengthThreshold <f> \n");
695 printf("\t method <AMGSA, AMGSAe> \n");
696 printf("\t coarsenScheme <local, hybrid, cljp, falgout> \n");
697 printf("\t smoother <Jacobi,GS,...> \n");
698 printf("\t coarseSolver <Jacobi,GS,...> \n");
699 printf("\t numSweeps <d> \n");
700 printf("\t smootherWeight <f> \n");
701 printf("\t smootherPrintRNorm \n");
702 printf("\t smootherFindOmega \n");
703 printf("\t minCoarseSize <d> \n");
704 printf("\t Pweight <f> \n");
705 printf("\t SPLevel <d> \n");
706 printf("\t scalar\n");
707 printf("\t nodeDOF <d> \n");
708 printf("\t nullSpaceDim <d> \n");
709 printf("\t useNodalCoord <on,off> \n");
710 printf("\t saAMGCalibrationSize <d> \n");
711 printf("\t rsAMGSymmetric <d> \n");
712 printf("\t rsAMGInjectionForR\n");
713 printf("\t printNullSpace\n");
714 printf("\t printElemNodeList \n");
715 printf("\t printNodalCoord \n");
716 printf("\t paramFile <s> \n");
717 printf("\t numSmoothVecs <d> \n");
718 printf("\t smoothVecSteps <d> \n");
719 printf("\t arpackTol <f> \n");
720 }
721 }
722 else if ( !strcmp(param2, "outputLevel") )
723 {
724 sscanf(paramString,"%s %s %d",param1,param2,&(mli_object->outputLevel_));
725 }
726 else if ( !strcmp(param2, "numLevels") )
727 {
728 sscanf(paramString,"%s %s %d", param1, param2, &(mli_object->nLevels_));
729 if ( mli_object->nLevels_ <= 0 ) mli_object->nLevels_ = 1;
730 }
731 else if ( !strcmp(param2, "maxIterations") )
732 {
733 sscanf(paramString,"%s %s %d", param1, param2,
734 &(mli_object->maxIterations_));
735 if ( mli_object->maxIterations_ <= 0 ) mli_object->maxIterations_ = 1;
736 }
737 else if ( !strcmp(param2, "cycleType") )
738 {
739 sscanf(paramString,"%s %s %s", param1, param2, param3);
740 if ( ! strcmp( param3, "V" ) ) mli_object->cycleType_ = 1;
741 else if ( ! strcmp( param3, "W" ) ) mli_object->cycleType_ = 2;
742 }
743 else if ( !strcmp(param2, "strengthThreshold") )
744 {
745 sscanf(paramString,"%s %s %lg", param1, param2,
746 &(mli_object->strengthThreshold_));
747 if ( mli_object->strengthThreshold_ < 0.0 )
748 mli_object->strengthThreshold_ = 0.0;
749 }
750 else if ( !strcmp(param2, "method") )
751 {
752 sscanf(paramString,"%s %s %s", param1, param2, param3);
753 strcpy( mli_object->method_, param3 );
754 }
755 else if ( !strcmp(param2, "coarsenScheme") )
756 {
757 sscanf(paramString,"%s %s %s", param1, param2, param3);
758 strcpy( mli_object->coarsenScheme_, param3 );
759 }
760 else if ( !strcmp(param2, "smoother") )
761 {
762 sscanf(paramString,"%s %s %s", param1, param2, param3);
763 strcpy( mli_object->preSmoother_, param3 );
764 strcpy( mli_object->postSmoother_, param3 );
765 }
766 else if ( !strcmp(param2, "coarseSolver") )
767 {
768 sscanf(paramString,"%s %s %s", param1, param2, param3);
769 strcpy( mli_object->coarseSolver_, param3 );
770 }
771 else if ( !strcmp(param2, "coarseSolverNumSweeps") )
772 {
773 sscanf(paramString,"%s %s %d", param1, param2,
774 &(mli_object->coarseSolverNSweeps_));
775 if ( mli_object->coarseSolverNSweeps_ < 1 )
776 mli_object->coarseSolverNSweeps_ = 1;
777 }
778 else if ( !strcmp(param2, "numSweeps") )
779 {
780 sscanf(paramString,"%s %s %d",param1,param2,&(mli_object->preNSweeps_));
781 if ( mli_object->preNSweeps_ <= 0 ) mli_object->preNSweeps_ = 1;
782 mli_object->postNSweeps_ = mli_object->preNSweeps_;
783 if ( mli_object->preSmootherWts_ != NULL )
784 {
785 weight = mli_object->preSmootherWts_[0];
786 delete [] mli_object->preSmootherWts_;
787 mli_object->preSmootherWts_ = new double[mli_object->preNSweeps_];
788 for ( i = 0; i < mli_object->preNSweeps_; i++ )
789 mli_object->preSmootherWts_[i] = weight;
790 }
791 if ( mli_object->postSmootherWts_ != NULL )
792 {
793 weight = mli_object->postSmootherWts_[0];
794 delete [] mli_object->postSmootherWts_;
795 mli_object->postSmootherWts_ = new double[mli_object->postNSweeps_];
796 for ( i = 0; i < mli_object->postNSweeps_; i++ )
797 mli_object->postSmootherWts_[i] = weight;
798 }
799 }
800 else if ( !strcmp(param2, "smootherWeight") )
801 {
802 sscanf(paramString,"%s %s %lg",param1,param2,&weight);
803 if ( weight < 0.0 || weight > 2.0 ) weight = 1.0;
804 if ( mli_object->preNSweeps_ > 0 )
805 {
806 if ( mli_object->preSmootherWts_ != NULL )
807 delete [] mli_object->preSmootherWts_;
808 mli_object->preSmootherWts_ = new double[mli_object->preNSweeps_];
809 for ( i = 0; i < mli_object->preNSweeps_; i++ )
810 mli_object->preSmootherWts_[i] = weight;
811 mli_object->postNSweeps_ = mli_object->preNSweeps_;
812 if ( mli_object->postSmootherWts_ != NULL )
813 delete [] mli_object->postSmootherWts_;
814 mli_object->postSmootherWts_ = new double[mli_object->preNSweeps_];
815 for ( i = 0; i < mli_object->preNSweeps_; i++ )
816 mli_object->postSmootherWts_[i] = weight;
817 }
818 }
819 else if ( !strcmp(param2, "smootherPrintRNorm") )
820 {
821 mli_object->smootherPrintRNorm_ = 1;
822 }
823 else if ( !strcmp(param2, "smootherFindOmega") )
824 {
825 mli_object->smootherFindOmega_ = 1;
826 }
827 else if ( !strcmp(param2, "minCoarseSize") )
828 {
829 sscanf(paramString,"%s %s %d",param1,param2,
830 &(mli_object->minCoarseSize_));
831 if ( mli_object->minCoarseSize_ <= 0 ) mli_object->minCoarseSize_ = 20;
832 }
833 else if ( !strcmp(param2, "Pweight") )
834 {
835 sscanf(paramString,"%s %s %lg",param1,param2, &(mli_object->Pweight_));
836 if ( mli_object->Pweight_ < 0. ) mli_object->Pweight_ = 1.333;
837 }
838 else if ( !strcmp(param2, "SPLevel") )
839 {
840 sscanf(paramString,"%s %s %d",param1,param2, &(mli_object->SPLevel_));
841 if ( mli_object->SPLevel_ < 0 ) mli_object->SPLevel_ = 0;
842 }
843 else if ( !strcmp(param2, "scalar") )
844 {
845 mli_object->scalar_ = 1;
846 }
847 else if ( !strcmp(param2, "nodeDOF") )
848 {
849 sscanf(paramString,"%s %s %d",param1,param2, &(mli_object->nodeDOF_));
850 if ( mli_object->nodeDOF_ <= 0 ) mli_object->nodeDOF_ = 1;
851 }
852 else if ( !strcmp(param2, "nullSpaceDim") )
853 {
854 sscanf(paramString,"%s %s %d",param1,param2, &(mli_object->nSpaceDim_));
855 if ( mli_object->nSpaceDim_ <= 0 ) mli_object->nSpaceDim_ = 1;
856 }
857 else if ( !strcmp(param2, "useNodalCoord") )
858 {
859 sscanf(paramString,"%s %s %s",param1,param2,param3);
860 if ( !strcmp(param3, "on") ) mli_object->nCoordAccept_ = 1;
861 else mli_object->nCoordAccept_ = 0;
862 }
863 else if ( !strcmp(param2, "saAMGCalibrationSize") )
864 {
865 sscanf(paramString,"%s %s %d",param1,param2,
866 &(mli_object->calibrationSize_));
867 if (mli_object->calibrationSize_ < 0) mli_object->calibrationSize_ = 0;
868 }
869 else if ( !strcmp(param2, "rsAMGSymmetric") )
870 {
871 sscanf(paramString,"%s %s %d",param1,param2, &(mli_object->symmetric_));
872 if ( mli_object->symmetric_ < 0 ) mli_object->symmetric_ = 0;
873 if ( mli_object->symmetric_ > 1 ) mli_object->symmetric_ = 1;
874 }
875 else if ( !strcmp(param2, "rsAMGInjectionForR") )
876 {
877 mli_object->injectionForR_ = 1;
878 }
879 else if ( !strcmp(param2, "printNullSpace") )
880 {
881 mli_object->printNullSpace_ |= 1;
882 }
883 else if ( !strcmp(param2, "printElemNodeList") )
884 {
885 mli_object->printNullSpace_ |= 2;
886 }
887 else if ( !strcmp(param2, "printNodalCoord") )
888 {
889 mli_object->printNullSpace_ |= 4;
890 }
891 else if ( !strcmp(param2, "paramFile") )
892 {
893 sscanf(paramString,"%s %s %s",param1,param2,mli_object->paramFile_);
894 }
895 else if ( !strcmp(param2, "numSmoothVecs") )
896 {
897 sscanf(paramString,"%s %s %d",param1,param2,
898 &(mli_object->numSmoothVecs_));
899 if ( mli_object->numSmoothVecs_ < 0 ) mli_object->numSmoothVecs_ = 0;
900 }
901 else if ( !strcmp(param2, "smoothVecSteps") )
902 {
903 sscanf(paramString,"%s %s %d",param1,param2,
904 &(mli_object->smoothVecSteps_));
905 if ( mli_object->smoothVecSteps_ < 0 ) mli_object->smoothVecSteps_ = 0;
906 }
907 else if ( !strcmp(param2, "arpackTol") )
908 {
909 sscanf(paramString,"%s %s %lg",param1,param2,
910 &(mli_object->arpackTol_));
911 if ( mli_object->arpackTol_ <= 0.0 ) mli_object->arpackTol_ = 0.0;
912 }
913 else if ( !strcmp(param2, "incrNullSpaceDim") )
914 {
915 sscanf(paramString,"%s %s %d",param1,param2, &i);
916 mli_object->nSpaceDim_ += i;
917 }
918 else
919 {
920 if ( mypid == 0 )
921 {
922 printf("%4d : HYPRE_LSI_MLISetParams ERROR : unrecognized request.\n",
923 mypid);
924 printf("\t offending request = %s.\n", paramString);
925 printf("\tAvailable options for MLI are : \n");
926 printf("\t outputLevel <d> \n");
927 printf("\t numLevels <d> \n");
928 printf("\t maxIterations <d> \n");
929 printf("\t cycleType <'V','W'> \n");
930 printf("\t strengthThreshold <f> \n");
931 printf("\t method <AMGSA, AMGSAe> \n");
932 printf("\t smoother <Jacobi,GS,...> \n");
933 printf("\t coarseSolver <Jacobi,GS,...> \n");
934 printf("\t numSweeps <d> \n");
935 printf("\t smootherWeight <f> \n");
936 printf("\t smootherPrintRNorm\n");
937 printf("\t smootherFindOmega\n");
938 printf("\t minCoarseSize <d> \n");
939 printf("\t Pweight <f> \n");
940 printf("\t SPLevel <d> \n");
941 printf("\t nodeDOF <d> \n");
942 printf("\t nullSpaceDim <d> \n");
943 printf("\t useNodalCoord <on,off> \n");
944 printf("\t saAMGCalibrationSize <d> \n");
945 printf("\t rsAMGSymmetric <d> \n");
946 printf("\t rsAMGInjectionForR\n");
947 printf("\t printNullSpace\n");
948 printf("\t printElemNodeList\n");
949 printf("\t printNodalCoord\n");
950 printf("\t paramFile <s> \n");
951 printf("\t numSmoothVecs <d> \n");
952 printf("\t smoothVecSteps <d> \n");
953 printf("\t arpackTol <f> \n");
954 exit(1);
955 }
956 }
957 return 0;
958 }
959
960 /****************************************************************************/
961 /* HYPRE_LSI_MLICreateNodeEqnMap */
962 /*--------------------------------------------------------------------------*/
963
964 extern "C"
HYPRE_LSI_MLICreateNodeEqnMap(HYPRE_Solver solver,int nNodes,int * nodeNumbers,int * eqnNumbers,int * procNRows)965 int HYPRE_LSI_MLICreateNodeEqnMap(HYPRE_Solver solver, int nNodes,
966 int *nodeNumbers, int *eqnNumbers,
967 int *procNRows)
968 {
969 #ifdef HAVE_MLI
970 int iN, iP, mypid, nprocs, *procMapArray, *iTempArray;
971 int iS, nSends, *sendLengs, *sendProcs, **iSendBufs;
972 int iR, nRecvs, *recvLengs, *recvProcs, **iRecvBufs, *procList;
973 int newNumNodes, *newNodeNumbers, *newEqnNumbers, procIndex;
974 MPI_Comm mpiComm;
975 MPI_Request *mpiRequests;
976 MPI_Status mpiStatus;
977 MLI_Mapper *mapper;
978 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
979
980 /* -------------------------------------------------------- */
981 /* fetch processor information */
982 /* -------------------------------------------------------- */
983
984 if ( mli_object == NULL ) return 1;
985 if ( mli_object->mapper_ != NULL ) delete mli_object->mapper_;
986 mpiComm = mli_object->mpiComm_;
987 MPI_Comm_rank( mpiComm, &mypid );
988 MPI_Comm_size( mpiComm, &nprocs );
989
990 /* -------------------------------------------------------- */
991 /* construct node to processor map array */
992 /* -------------------------------------------------------- */
993
994 procMapArray = new int[nNodes];
995 for ( iN = 0; iN < nNodes; iN++ )
996 {
997 procMapArray[iN] = -1;
998 if ( eqnNumbers[iN] < procNRows[mypid] ||
999 eqnNumbers[iN] >= procNRows[mypid+1] )
1000 {
1001 for ( iP = 0; iP < nprocs; iP++ )
1002 if ( eqnNumbers[iN] < procNRows[iP] ) break;
1003 procMapArray[iN] = iP - 1;
1004 }
1005 }
1006
1007 /* -------------------------------------------------------- */
1008 /* construct send information */
1009 /* -------------------------------------------------------- */
1010
1011 procList = new int[nprocs];
1012 for ( iP = 0; iP < nprocs; iP++ ) procList[iP] = 0;
1013 for ( iN = 0; iN < nNodes; iN++ )
1014 if ( procMapArray[iN] >= 0 ) procList[procMapArray[iN]]++;
1015 nSends = 0;
1016 for ( iP = 0; iP < nprocs; iP++ ) if ( procList[iP] > 0 ) nSends++;
1017 if ( nSends > 0 )
1018 {
1019 sendProcs = new int[nSends];
1020 sendLengs = new int[nSends];
1021 iSendBufs = new int*[nSends];
1022 }
1023 nSends = 0;
1024 for ( iP = 0; iP < nprocs; iP++ )
1025 {
1026 if ( procList[iP] > 0 )
1027 {
1028 sendLengs[nSends] = procList[iP];
1029 sendProcs[nSends++] = iP;
1030 }
1031 }
1032
1033 /* -------------------------------------------------------- */
1034 /* construct recv information */
1035 /* -------------------------------------------------------- */
1036
1037 for ( iP = 0; iP < nprocs; iP++ ) procList[iP] = 0;
1038 for ( iP = 0; iP < nSends; iP++ ) procList[sendProcs[iP]]++;
1039 iTempArray = new int[nprocs];
1040 MPI_Allreduce(procList,iTempArray,nprocs,MPI_INT,MPI_SUM,mpiComm);
1041 nRecvs = iTempArray[mypid];
1042 delete [] procList;
1043 delete [] iTempArray;
1044 if ( nRecvs > 0 )
1045 {
1046 recvLengs = new int[nRecvs];
1047 recvProcs = new int[nRecvs];
1048 iRecvBufs = new int*[nRecvs];
1049 mpiRequests = new MPI_Request[nRecvs];
1050 }
1051 for ( iP = 0; iP < nRecvs; iP++ )
1052 MPI_Irecv(&(recvLengs[iP]), 1, MPI_INT, MPI_ANY_SOURCE, 29421,
1053 mpiComm, &(mpiRequests[iP]));
1054 for ( iP = 0; iP < nSends; iP++ )
1055 MPI_Send(&(sendLengs[iP]), 1, MPI_INT, sendProcs[iP], 29421, mpiComm);
1056 for ( iP = 0; iP < nRecvs; iP++ )
1057 {
1058 MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
1059 recvProcs[iP] = mpiStatus.MPI_SOURCE;
1060 }
1061
1062 /* -------------------------------------------------------- */
1063 /* communicate node and equation information */
1064 /* -------------------------------------------------------- */
1065
1066 for ( iP = 0; iP < nRecvs; iP++ )
1067 {
1068 iRecvBufs[iP] = new int[recvLengs[iP]*2];
1069 MPI_Irecv(iRecvBufs[iP], recvLengs[iP]*2, MPI_INT, recvProcs[iP],
1070 29422, mpiComm, &(mpiRequests[iP]));
1071 }
1072 for ( iP = 0; iP < nSends; iP++ )
1073 {
1074 iSendBufs[iP] = new int[sendLengs[iP]*2];
1075 sendLengs[iP] = 0;
1076 }
1077 for ( iN = 0; iN < nNodes; iN++ )
1078 {
1079 if ( procMapArray[iN] >= 0 )
1080 {
1081 procIndex = procMapArray[iN];
1082 for ( iP = 0; iP < nSends; iP++ )
1083 if ( sendProcs[iP] == procIndex ) break;
1084 iSendBufs[iP][sendLengs[iP]++] = nodeNumbers[iN];
1085 iSendBufs[iP][sendLengs[iP]++] = eqnNumbers[iN];
1086 }
1087 }
1088 for ( iP = 0; iP < nSends; iP++ )
1089 {
1090 sendLengs[iP] /= 2;
1091 MPI_Send(iSendBufs[iP], sendLengs[iP]*2, MPI_INT, sendProcs[iP],
1092 29422, mpiComm);
1093 }
1094 for ( iP = 0; iP < nRecvs; iP++ ) MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
1095
1096 /* -------------------------------------------------------- */
1097 /* form node and equation map */
1098 /* -------------------------------------------------------- */
1099
1100 newNumNodes = nNodes;
1101 for ( iP = 0; iP < nRecvs; iP++ ) newNumNodes += recvLengs[iP];
1102 newNodeNumbers = new int[newNumNodes];
1103 newEqnNumbers = new int[newNumNodes];
1104 newNumNodes = 0;
1105 for (iN = 0; iN < nNodes; iN++)
1106 {
1107 newNodeNumbers[newNumNodes] = nodeNumbers[iN];
1108 newEqnNumbers[newNumNodes++] = eqnNumbers[iN];
1109 }
1110 for ( iP = 0; iP < nRecvs; iP++ )
1111 {
1112 for ( iR = 0; iR < recvLengs[iP]; iR++ )
1113 {
1114 newNodeNumbers[newNumNodes] = iRecvBufs[iP][iR*2];
1115 newEqnNumbers[newNumNodes++] = iRecvBufs[iP][iR*2+1];
1116 }
1117 }
1118 mapper = new MLI_Mapper();
1119 mapper->setMap( newNumNodes, newNodeNumbers, newEqnNumbers );
1120 mli_object->mapper_ = mapper;
1121
1122 /* -------------------------------------------------------- */
1123 /* clean up and return */
1124 /* -------------------------------------------------------- */
1125
1126 delete [] procMapArray;
1127 if ( nSends > 0 )
1128 {
1129 delete [] sendProcs;
1130 delete [] sendLengs;
1131 for ( iS = 0; iS < nSends; iS++ ) delete [] iSendBufs[iS];
1132 delete [] iSendBufs;
1133 }
1134 if ( nRecvs > 0 )
1135 {
1136 delete [] recvProcs;
1137 delete [] recvLengs;
1138 for ( iR = 0; iR < nRecvs; iR++ ) delete [] iRecvBufs[iR];
1139 delete [] iRecvBufs;
1140 delete [] mpiRequests;
1141 }
1142 delete [] newNodeNumbers;
1143 delete [] newEqnNumbers;
1144 return 0;
1145 #else
1146 (void) solver;
1147 (void) nNodes;
1148 (void) nodeNumbers;
1149 (void) eqnNumbers;
1150 (void) procNRows;
1151 return 1;
1152 #endif
1153 }
1154
1155 /****************************************************************************/
1156 /* HYPRE_LSI_MLIAdjustNodeEqnMap */
1157 /*--------------------------------------------------------------------------*/
1158
1159 extern "C"
HYPRE_LSI_MLIAdjustNodeEqnMap(HYPRE_Solver solver,int * procNRows,int * procOffsets)1160 int HYPRE_LSI_MLIAdjustNodeEqnMap(HYPRE_Solver solver, int *procNRows,
1161 int *procOffsets)
1162 {
1163 #ifdef HAVE_MLI
1164 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1165 if ( mli_object == NULL ) return 1;
1166 if ( mli_object->mapper_ == NULL ) return 1;
1167 mli_object->mapper_->adjustMapOffset( mli_object->mpiComm_, procNRows,
1168 procOffsets );
1169 return 0;
1170 #else
1171 (void) solver;
1172 (void) procNRows;
1173 (void) procOffsets;
1174 return 1;
1175 #endif
1176 }
1177
1178 /****************************************************************************/
1179 /* HYPRE_LSI_MLIAdjustNullSpace */
1180 /*--------------------------------------------------------------------------*/
1181
1182 extern "C"
HYPRE_LSI_MLIAdjustNullSpace(HYPRE_Solver solver,int nConstraints,int * slaveIndices,HYPRE_ParCSRMatrix hypreA)1183 int HYPRE_LSI_MLIAdjustNullSpace(HYPRE_Solver solver, int nConstraints,
1184 int *slaveIndices,
1185 HYPRE_ParCSRMatrix hypreA)
1186 {
1187 #ifdef HAVE_MLI
1188 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1189 if ( mli_object == NULL ) return 1;
1190 mli_object->adjustNullSpace_ = 1;
1191 mli_object->numResetNull_ = nConstraints;
1192 if ( nConstraints > 0 )
1193 mli_object->resetNullIndices_ = new int[nConstraints];
1194 for ( int i = 0; i < nConstraints; i++ )
1195 mli_object->resetNullIndices_[i] = slaveIndices[i];
1196 mli_object->correctionMatrix_ = hypreA;
1197 return 0;
1198 #else
1199 (void) solver;
1200 (void) nConstraints;
1201 (void) slaveIndices;
1202 (void) hypreA;
1203 return 1;
1204 #endif
1205 }
1206
1207 /****************************************************************************/
1208 /* HYPRE_LSI_MLISetFEData */
1209 /*--------------------------------------------------------------------------*/
1210
1211 extern "C"
HYPRE_LSI_MLISetFEData(HYPRE_Solver solver,void * object)1212 int HYPRE_LSI_MLISetFEData(HYPRE_Solver solver, void *object)
1213 {
1214 #ifdef HAVE_MLI
1215 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1216 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1217 mli_object->feData_ = hypre_fedata->fedata_;
1218 hypre_fedata->fedata_ = NULL;
1219 hypre_fedata->fedataOwn_ = 0;
1220 return 0;
1221 #else
1222 (void) solver;
1223 (void) object;
1224 return 1;
1225 #endif
1226 }
1227
1228 /****************************************************************************/
1229 /* HYPRE_LSI_MLISetSFEI */
1230 /*--------------------------------------------------------------------------*/
1231
1232 extern "C"
HYPRE_LSI_MLISetSFEI(HYPRE_Solver solver,void * object)1233 int HYPRE_LSI_MLISetSFEI(HYPRE_Solver solver, void *object)
1234 {
1235 #ifdef HAVE_MLI
1236 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1237 HYPRE_MLI_SFEI *hypre_sfei = (HYPRE_MLI_SFEI *) object;
1238 mli_object->sfei_ = hypre_sfei->sfei_;
1239 hypre_sfei->sfei_ = NULL;
1240 hypre_sfei->sfeiOwn_ = 0;
1241 return 0;
1242 #else
1243 (void) solver;
1244 (void) object;
1245 return 1;
1246 #endif
1247 }
1248
1249 /****************************************************************************/
1250 /* HYPRE_LSI_MLILoadNodalCoordinates */
1251 /* (The nodal coordinates loaded in here conforms to the nodal labeling in */
1252 /* FEI, so the lookup object can be used to find the equation number. */
1253 /* In addition, node numbers and coordinates need to be shuffled between */
1254 /* processors) */
1255 /*--------------------------------------------------------------------------*/
1256
1257 /*
1258 #define FEI_250
1259 */
1260
1261 extern "C"
HYPRE_LSI_MLILoadNodalCoordinates(HYPRE_Solver solver,int nNodes,int nodeDOF,int * eqnNumbers,int nDim,double * coords,int localNRows)1262 int HYPRE_LSI_MLILoadNodalCoordinates(HYPRE_Solver solver, int nNodes,
1263 int nodeDOF, int *eqnNumbers, int nDim, double *coords,
1264 int localNRows)
1265 {
1266 int iN, iD, eqnInd, mypid, *flags, arrayLeng;
1267 double *nCoords;
1268 MPI_Comm mpiComm;
1269 #ifdef FEI_250
1270 int iMax, iMin, offFlag;
1271 #else
1272 int iP, nprocs, *nodeProcMap, *iTempArray;
1273 int iS, nSends, *sendLengs, *sendProcs, **iSendBufs, procIndex;
1274 int iR, nRecvs, *recvLengs, *recvProcs, **iRecvBufs, *procList;
1275 int *procNRows, numNodes, coordLength;
1276 double **dSendBufs, **dRecvBufs;
1277 MPI_Request *mpiRequests;
1278 MPI_Status mpiStatus;
1279 #endif
1280 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1281
1282 /* -------------------------------------------------------- */
1283 /* if nodal coordinates flag off, do not take coordinates */
1284 /* -------------------------------------------------------- */
1285
1286 if ( ! mli_object->nCoordAccept_ ) return 1;
1287 mpiComm = mli_object->mpiComm_;
1288 MPI_Comm_rank( mpiComm, &mypid );
1289
1290 /* -------------------------------------------------------- */
1291 /* clean up previously allocated space */
1292 /* -------------------------------------------------------- */
1293
1294 if ( mli_object->nCoordinates_ != NULL )
1295 delete [] mli_object->nCoordinates_;
1296 if ( mli_object->nullScales_ != NULL )
1297 delete [] mli_object->nullScales_;
1298 mli_object->nCoordinates_ = NULL;
1299 mli_object->nullScales_ = NULL;
1300
1301 /* -------------------------------------------------------- */
1302 /* This code is used in place of the 'else' block in view */
1303 /* of the changes made to FEI 2.5.0 */
1304 /* -------------------------------------------------------- */
1305
1306 #ifdef FEI_250
1307 mli_object->spaceDim_ = nDim;
1308 mli_object->nodeDOF_ = nodeDOF;
1309 iMin = 1000000000; iMax = 0;
1310 for ( iN = 0; iN < nNodes; iN++ )
1311 {
1312 iMin = ( eqnNumbers[iN] < iMin ) ? eqnNumbers[iN] : iMin;
1313 iMax = ( eqnNumbers[iN] > iMax ) ? eqnNumbers[iN] : iMax;
1314 }
1315 mli_object->localNEqns_ = ( iMax/nDim - iMin/nDim + 1 ) * nDim;
1316 offFlag = 0;
1317 if ( mli_object->localNEqns_ != nNodes*nDim ) offFlag = 1;
1318 iN = offFlag;
1319 MPI_Allreduce(&iN,&offFlag,1,MPI_INT,MPI_SUM,mpiComm);
1320 if ( offFlag != 0 )
1321 {
1322 if ( mypid == 0 )
1323 printf("HYPRE_LSI_MLILoadNodalCoordinates - turned off.\n");
1324 mli_object->nCoordAccept_ = 0;
1325 mli_object->localNEqns_ = 0;
1326 return 1;
1327 }
1328 mli_object->nCoordinates_ = new double[nNodes*nDim];
1329 nCoords = mli_object->nCoordinates_;
1330 for ( iN = 0; iN < nNodes; iN++ )
1331 {
1332 eqnInd = (eqnNumbers[iN] - iMin) / nodeDOF;
1333 for ( iD = 0; iD < nDim; iD++ )
1334 nCoords[eqnInd*nDim+iD] = coords[iN*nDim+iD];
1335 }
1336
1337 #else
1338 /* -------------------------------------------------------- */
1339 /* fetch machine information */
1340 /* -------------------------------------------------------- */
1341
1342 MPI_Comm_size( mpiComm, &nprocs );
1343
1344 /* -------------------------------------------------------- */
1345 /* construct procNRows array */
1346 /* -------------------------------------------------------- */
1347
1348 procNRows = new int[nprocs+1];
1349 iTempArray = new int[nprocs];
1350 for ( iP = 0; iP <= nprocs; iP++ ) procNRows[iP] = 0;
1351 procNRows[mypid] = localNRows;
1352 MPI_Allreduce(procNRows,iTempArray,nprocs,MPI_INT,MPI_SUM,mpiComm);
1353 procNRows[0] = 0;
1354 for ( iP = 1; iP <= nprocs; iP++ )
1355 procNRows[iP] = procNRows[iP-1] + iTempArray[iP-1];
1356
1357 /* -------------------------------------------------------- */
1358 /* construct node to processor map */
1359 /* -------------------------------------------------------- */
1360
1361 nodeProcMap = new int[nNodes];
1362 for ( iN = 0; iN < nNodes; iN++ )
1363 {
1364 nodeProcMap[iN] = -1;
1365 if ( eqnNumbers[iN] < procNRows[mypid] ||
1366 eqnNumbers[iN] >= procNRows[mypid+1] )
1367 {
1368 for ( iP = 0; iP < nprocs; iP++ )
1369 if ( eqnNumbers[iN] < procNRows[iP] ) break;
1370 nodeProcMap[iN] = iP - 1;
1371 }
1372 }
1373
1374 /* -------------------------------------------------------- */
1375 /* construct send information */
1376 /* -------------------------------------------------------- */
1377
1378 procList = new int[nprocs];
1379 for ( iP = 0; iP < nprocs; iP++ ) procList[iP] = 0;
1380 for ( iN = 0; iN < nNodes; iN++ )
1381 if ( nodeProcMap[iN] >= 0 ) procList[nodeProcMap[iN]]++;
1382 nSends = 0;
1383 for ( iP = 0; iP < nprocs; iP++ ) if ( procList[iP] > 0 ) nSends++;
1384 if ( nSends > 0 )
1385 {
1386 sendProcs = new int[nSends];
1387 sendLengs = new int[nSends];
1388 iSendBufs = new int*[nSends];
1389 dSendBufs = new double*[nSends];
1390 }
1391 nSends = 0;
1392 for ( iP = 0; iP < nprocs; iP++ )
1393 {
1394 if ( procList[iP] > 0 )
1395 {
1396 sendLengs[nSends] = procList[iP];
1397 sendProcs[nSends++] = iP;
1398 }
1399 }
1400
1401 /* -------------------------------------------------------- */
1402 /* construct recv information */
1403 /* -------------------------------------------------------- */
1404
1405 for ( iP = 0; iP < nprocs; iP++ ) procList[iP] = 0;
1406 for ( iP = 0; iP < nSends; iP++ ) procList[sendProcs[iP]]++;
1407 MPI_Allreduce(procList,iTempArray,nprocs,MPI_INT,MPI_SUM,mpiComm);
1408 nRecvs = iTempArray[mypid];
1409 if ( nRecvs > 0 )
1410 {
1411 recvLengs = new int[nRecvs];
1412 recvProcs = new int[nRecvs];
1413 iRecvBufs = new int*[nRecvs];
1414 dRecvBufs = new double*[nRecvs];
1415 mpiRequests = new MPI_Request[nRecvs];
1416 }
1417 for ( iP = 0; iP < nRecvs; iP++ )
1418 MPI_Irecv(&(recvLengs[iP]), 1, MPI_INT, MPI_ANY_SOURCE, 29421,
1419 mpiComm, &(mpiRequests[iP]));
1420 for ( iP = 0; iP < nSends; iP++ )
1421 MPI_Send(&(sendLengs[iP]), 1, MPI_INT, sendProcs[iP], 29421, mpiComm);
1422 for ( iP = 0; iP < nRecvs; iP++ )
1423 {
1424 MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
1425 recvProcs[iP] = mpiStatus.MPI_SOURCE;
1426 }
1427
1428 /* -------------------------------------------------------- */
1429 /* communicate equation numbers information */
1430 /* -------------------------------------------------------- */
1431
1432 for ( iP = 0; iP < nRecvs; iP++ )
1433 {
1434 iRecvBufs[iP] = new int[recvLengs[iP]];
1435 MPI_Irecv(iRecvBufs[iP], recvLengs[iP], MPI_INT, recvProcs[iP],
1436 29422, mpiComm, &(mpiRequests[iP]));
1437 }
1438 for ( iP = 0; iP < nSends; iP++ )
1439 {
1440 iSendBufs[iP] = new int[sendLengs[iP]];
1441 sendLengs[iP] = 0;
1442 }
1443 for ( iN = 0; iN < nNodes; iN++ )
1444 {
1445 if ( nodeProcMap[iN] >= 0 )
1446 {
1447 procIndex = nodeProcMap[iN];
1448 for ( iP = 0; iP < nSends; iP++ )
1449 if ( procIndex == sendProcs[iP] ) break;
1450 iSendBufs[iP][sendLengs[iP]++] = eqnNumbers[iN];
1451 }
1452 }
1453 for ( iP = 0; iP < nSends; iP++ )
1454 {
1455 MPI_Send(iSendBufs[iP], sendLengs[iP], MPI_INT, sendProcs[iP],
1456 29422, mpiComm);
1457 }
1458 for ( iP = 0; iP < nRecvs; iP++ ) MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
1459
1460 /* -------------------------------------------------------- */
1461 /* communicate coordinate information */
1462 /* -------------------------------------------------------- */
1463
1464 for ( iP = 0; iP < nRecvs; iP++ )
1465 {
1466 dRecvBufs[iP] = new double[recvLengs[iP]*nDim];
1467 MPI_Irecv(dRecvBufs[iP], recvLengs[iP]*nDim, MPI_DOUBLE,
1468 recvProcs[iP], 29425, mpiComm, &(mpiRequests[iP]));
1469 }
1470 for ( iP = 0; iP < nSends; iP++ )
1471 {
1472 dSendBufs[iP] = new double[sendLengs[iP]*nDim];
1473 sendLengs[iP] = 0;
1474 }
1475 for ( iN = 0; iN < nNodes; iN++ )
1476 {
1477 if ( nodeProcMap[iN] >= 0 )
1478 {
1479 procIndex = nodeProcMap[iN];
1480 for ( iP = 0; iP < nSends; iP++ )
1481 if ( procIndex == sendProcs[iP] ) break;
1482 for ( iD = 0; iD < nDim; iD++ )
1483 dSendBufs[iP][sendLengs[iP]++]=coords[iN*nDim+iD];
1484 }
1485 }
1486 for ( iP = 0; iP < nSends; iP++ )
1487 {
1488 sendLengs[iP] /= nDim;
1489 MPI_Send(dSendBufs[iP], sendLengs[iP]*nDim, MPI_DOUBLE,
1490 sendProcs[iP], 29425, mpiComm);
1491 }
1492 for ( iP = 0; iP < nRecvs; iP++ ) MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
1493
1494 /* -------------------------------------------------------- */
1495 /* check any duplicate coordinate information */
1496 /* -------------------------------------------------------- */
1497
1498 arrayLeng = nNodes;
1499 for ( iP = 0; iP < nRecvs; iP++ ) arrayLeng += recvLengs[iP];
1500 flags = new int[arrayLeng];
1501 for ( iN = 0; iN < arrayLeng; iN++ ) flags[iN] = 0;
1502 for ( iN = 0; iN < nNodes; iN++ )
1503 {
1504 if ( nodeProcMap[iN] < 0 )
1505 {
1506 eqnInd = (eqnNumbers[iN] - procNRows[mypid]) / nodeDOF;
1507 if ( eqnInd >= arrayLeng )
1508 {
1509 printf("%d : HYPRE_LSI_MLILoadNodalCoordinates - ERROR(1).\n",
1510 mypid);
1511 exit(1);
1512 }
1513 flags[eqnInd] = 1;
1514 }
1515 }
1516 for ( iP = 0; iP < nRecvs; iP++ )
1517 {
1518 for ( iR = 0; iR < recvLengs[iP]; iR++ )
1519 {
1520 eqnInd = (iRecvBufs[iP][iR] - procNRows[mypid]) / nodeDOF;
1521 if ( eqnInd >= arrayLeng )
1522 {
1523 printf("%d : HYPRE_LSI_MLILoadNodalCoordinates - ERROR(2).\n",
1524 mypid);
1525 exit(1);
1526 }
1527 flags[eqnInd] = 1;
1528 }
1529 }
1530 numNodes = 0;
1531 for ( iN = 0; iN < arrayLeng; iN++ )
1532 if ( flags[iN] == 0 ) break;
1533 else numNodes++;
1534 delete [] flags;
1535
1536 /* -------------------------------------------------------- */
1537 /* set up nodal coordinate information in correct order */
1538 /* -------------------------------------------------------- */
1539
1540 mli_object->spaceDim_ = nDim;
1541 mli_object->nodeDOF_ = nodeDOF;
1542 mli_object->localNEqns_ = numNodes * nodeDOF;
1543 coordLength = numNodes * nodeDOF * nDim;
1544 mli_object->nCoordinates_ = new double[coordLength];
1545 nCoords = mli_object->nCoordinates_;
1546
1547 arrayLeng = numNodes * nodeDOF;
1548 for ( iN = 0; iN < nNodes; iN++ )
1549 {
1550 if ( nodeProcMap[iN] < 0 )
1551 {
1552 eqnInd = (eqnNumbers[iN] - procNRows[mypid]) / nodeDOF;
1553 if ( eqnInd >= 0 && eqnInd < arrayLeng )
1554 for ( iD = 0; iD < nDim; iD++ )
1555 nCoords[eqnInd*nDim+iD] = coords[iN*nDim+iD];
1556 }
1557 }
1558 for ( iP = 0; iP < nRecvs; iP++ )
1559 {
1560 for ( iR = 0; iR < recvLengs[iP]; iR++ )
1561 {
1562 eqnInd = (iRecvBufs[iP][iR] - procNRows[mypid]) / nodeDOF;
1563 if ( eqnInd >= 0 && eqnInd < arrayLeng )
1564 for ( iD = 0; iD < nDim; iD++ )
1565 nCoords[eqnInd*nDim+iD] = dRecvBufs[iP][iR*nDim+iD];
1566 }
1567 }
1568 for ( iN = 0; iN < numNodes*nodeDOF; iN++ )
1569 if (nCoords[iN] == -99999.0) printf("%d : LSI_mli error %d\n",mypid,iN);
1570
1571 /* -------------------------------------------------------- */
1572 /* clean up */
1573 /* -------------------------------------------------------- */
1574
1575 delete [] procList;
1576 delete [] iTempArray;
1577 delete [] nodeProcMap;
1578 delete [] procNRows;
1579 if ( nSends > 0 )
1580 {
1581 delete [] sendProcs;
1582 delete [] sendLengs;
1583 for ( iS = 0; iS < nSends; iS++ ) delete [] iSendBufs[iS];
1584 for ( iS = 0; iS < nSends; iS++ ) delete [] dSendBufs[iS];
1585 delete [] dSendBufs;
1586 delete [] iSendBufs;
1587 }
1588 if ( nRecvs > 0 )
1589 {
1590 delete [] recvProcs;
1591 delete [] recvLengs;
1592 for ( iR = 0; iR < nRecvs; iR++ ) delete [] iRecvBufs[iR];
1593 for ( iR = 0; iR < nRecvs; iR++ ) delete [] dRecvBufs[iR];
1594 delete [] iRecvBufs;
1595 delete [] dRecvBufs;
1596 delete [] mpiRequests;
1597 }
1598 #endif
1599 return 0;
1600 }
1601
1602 /****************************************************************************/
1603 /* HYPRE_LSI_MLILoadMatrixScalings */
1604 /*--------------------------------------------------------------------------*/
1605
1606 extern "C"
HYPRE_LSI_MLILoadMatrixScalings(HYPRE_Solver solver,int nEqns,double * scalings)1607 int HYPRE_LSI_MLILoadMatrixScalings(HYPRE_Solver solver, int nEqns,
1608 double *scalings)
1609 {
1610 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1611 if ( scalings != NULL )
1612 {
1613 mli_object->nullScales_ = new double[nEqns];
1614 for ( int i = 0; i < nEqns; i++ )
1615 mli_object->nullScales_[i] = scalings[i];
1616 }
1617 return 0;
1618 }
1619
1620 /****************************************************************************/
1621 /* HYPRE_LSI_MLILoadMaterialLabels */
1622 /*--------------------------------------------------------------------------*/
1623
1624 extern "C"
HYPRE_LSI_MLILoadMaterialLabels(HYPRE_Solver solver,int nLabels,int * labels)1625 int HYPRE_LSI_MLILoadMaterialLabels(HYPRE_Solver solver, int nLabels,
1626 int *labels)
1627 {
1628 HYPRE_LSI_MLI *mli_object = (HYPRE_LSI_MLI *) solver;
1629 if ( labels != NULL )
1630 {
1631 mli_object->matLabels_ = new int[nLabels];
1632 for ( int i = 0; i < nLabels; i++ )
1633 mli_object->matLabels_[i] = labels[i];
1634 mli_object->numMatLabels_ = nLabels;
1635 }
1636 return 0;
1637 }
1638
1639 /****************************************************************************/
1640 /* HYPRE_LSI_MLIFEDataCreate */
1641 /*--------------------------------------------------------------------------*/
1642
1643 extern "C"
HYPRE_LSI_MLIFEDataCreate(MPI_Comm mpi_comm)1644 void *HYPRE_LSI_MLIFEDataCreate( MPI_Comm mpi_comm )
1645 {
1646 #ifdef HAVE_MLI
1647 HYPRE_MLI_FEData *hypre_fedata;
1648 hypre_fedata = hypre_TAlloc(HYPRE_MLI_FEData, 1, HYPRE_MEMORY_HOST);
1649 hypre_fedata->comm_ = mpi_comm;
1650 hypre_fedata->fedata_ = NULL;
1651 hypre_fedata->fedataOwn_ = 0;
1652 hypre_fedata->computeNull_ = 0;
1653 hypre_fedata->nullDim_ = 1;
1654 return ((void *) hypre_fedata);
1655 #else
1656 return NULL;
1657 #endif
1658 }
1659
1660 /****************************************************************************/
1661 /* HYPRE_LSI_MLIFEDataDestroy */
1662 /*--------------------------------------------------------------------------*/
1663
1664 extern "C"
HYPRE_LSI_MLIFEDataDestroy(void * object)1665 int HYPRE_LSI_MLIFEDataDestroy( void *object )
1666 {
1667 #ifdef HAVE_MLI
1668 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1669 if ( hypre_fedata == NULL ) return 1;
1670 if ( hypre_fedata->fedataOwn_ && hypre_fedata->fedata_ != NULL )
1671 delete hypre_fedata->fedata_;
1672 hypre_fedata->fedata_ = NULL;
1673 free( hypre_fedata );
1674 return 0;
1675 #else
1676 return 1;
1677 #endif
1678 }
1679
1680 /****************************************************************************/
1681 /* HYPRE_LSI_MLIFEDataInitFields */
1682 /*--------------------------------------------------------------------------*/
1683
1684 extern "C"
HYPRE_LSI_MLIFEDataInitFields(void * object,int nFields,int * fieldSizes,int * fieldIDs)1685 int HYPRE_LSI_MLIFEDataInitFields( void *object, int nFields, int *fieldSizes,
1686 int *fieldIDs )
1687 {
1688 #ifdef HAVE_MLI
1689 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1690 MLI_FEData *fedata;
1691
1692 if ( hypre_fedata == NULL ) return 1;
1693 if ( hypre_fedata->fedata_ != NULL ) delete hypre_fedata->fedata_;
1694 hypre_fedata->fedata_ = new MLI_FEData(hypre_fedata->comm_);
1695 hypre_fedata->fedataOwn_ = 1;
1696 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1697 fedata->initFields( nFields, fieldSizes, fieldIDs );
1698 return 0;
1699 #else
1700 (void) object;
1701 (void) nFields;
1702 (void) fieldSizes;
1703 (void) fieldIDs;
1704 return 1;
1705 #endif
1706 }
1707
1708 /****************************************************************************/
1709 /* HYPRE_LSI_MLIFEDataInitElemBlock */
1710 /*--------------------------------------------------------------------------*/
1711
1712 extern "C"
HYPRE_LSI_MLIFEDataInitElemBlock(void * object,int nElems,int nNodesPerElem,int numNodeFields,int * nodeFieldIDs)1713 int HYPRE_LSI_MLIFEDataInitElemBlock(void *object, int nElems,
1714 int nNodesPerElem, int numNodeFields,
1715 int *nodeFieldIDs)
1716 {
1717 #ifdef HAVE_MLI
1718 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1719 MLI_FEData *fedata;
1720
1721 if ( hypre_fedata == NULL ) return 1;
1722 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1723 if ( fedata == NULL ) return 1;
1724 if ( numNodeFields != 1 ) return 1;
1725 hypre_fedata->fedata_->initElemBlock(nElems,nNodesPerElem,numNodeFields,
1726 nodeFieldIDs,0,NULL);
1727 return 0;
1728 #else
1729 (void) object;
1730 (void) nElems;
1731 (void) nNodesPerElem;
1732 (void) numNodeFields;
1733 (void) nodeFieldIDs;
1734 return 1;
1735 #endif
1736 }
1737
1738 /****************************************************************************/
1739 /* HYPRE_LSI_MLIFEDataInitElemNodeList */
1740 /*--------------------------------------------------------------------------*/
1741
1742 extern "C"
HYPRE_LSI_MLIFEDataInitElemNodeList(void * object,int elemID,int nNodesPerElem,int * elemNodeList)1743 int HYPRE_LSI_MLIFEDataInitElemNodeList( void *object, int elemID,
1744 int nNodesPerElem, int *elemNodeList)
1745 {
1746 #ifdef HAVE_MLI
1747 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1748 MLI_FEData *fedata;
1749
1750 if ( hypre_fedata == NULL ) return 1;
1751 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1752 if ( fedata == NULL ) return 1;
1753 fedata->initElemNodeList(elemID,nNodesPerElem,elemNodeList,3,NULL);
1754 return 0;
1755 #else
1756 return 1;
1757 #endif
1758 }
1759
1760 /****************************************************************************/
1761 /* HYPRE_LSI_MLIFEDataInitSharedNodes */
1762 /*--------------------------------------------------------------------------*/
1763
1764 extern "C"
HYPRE_LSI_MLIFEDataInitSharedNodes(void * object,int nSharedNodes,int * sharedNodeIDs,int * sharedProcLengs,int ** sharedProcIDs)1765 int HYPRE_LSI_MLIFEDataInitSharedNodes(void *object, int nSharedNodes,
1766 int *sharedNodeIDs, int *sharedProcLengs,
1767 int **sharedProcIDs)
1768 {
1769 #ifdef HAVE_MLI
1770 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1771 MLI_FEData *fedata;
1772
1773 if ( hypre_fedata == NULL ) return 1;
1774 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1775 if ( fedata == NULL ) return 1;
1776 if ( nSharedNodes > 0 )
1777 fedata->initSharedNodes(nSharedNodes, sharedNodeIDs, sharedProcLengs,
1778 sharedProcIDs);
1779 return 0;
1780 #else
1781 (void) object;
1782 (void) nSharedNodes;
1783 (void) sharedNodeIDs;
1784 (void) sharedProcLengs;
1785 (void) sharedProcIDs;
1786 return 1;
1787 #endif
1788 }
1789
1790 /****************************************************************************/
1791 /* HYPRE_LSI_MLIFEDataInitComplete */
1792 /*--------------------------------------------------------------------------*/
1793
1794 extern "C"
HYPRE_LSI_MLIFEDataInitComplete(void * object)1795 int HYPRE_LSI_MLIFEDataInitComplete( void *object )
1796 {
1797 #ifdef HAVE_MLI
1798 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1799 MLI_FEData *fedata;
1800
1801 if ( hypre_fedata == NULL ) return 1;
1802 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1803 if ( fedata == NULL ) return 1;
1804 fedata->initComplete();
1805 return 0;
1806 #else
1807 (void) object;
1808 return 1;
1809 #endif
1810 }
1811
1812 /****************************************************************************/
1813 /* HYPRE_LSI_MLIFEDataLoadElemMatrix */
1814 /*--------------------------------------------------------------------------*/
1815
1816 extern "C"
HYPRE_LSI_MLIFEDataLoadElemMatrix(void * object,int elemID,int nNodes,int * nodeList,int matDim,double ** inMat)1817 int HYPRE_LSI_MLIFEDataLoadElemMatrix(void *object, int elemID, int nNodes,
1818 int *nodeList, int matDim, double **inMat)
1819 {
1820 (void) nNodes;
1821 (void) nodeList;
1822 #ifdef HAVE_MLI
1823 int i, j;
1824 double *elemMat;
1825 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1826 MLI_FEData *fedata;
1827
1828 /* -------------------------------------------------------- */
1829 /* error checking */
1830 /* -------------------------------------------------------- */
1831
1832 if ( hypre_fedata == NULL ) return 1;
1833 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1834 if ( fedata == NULL ) return 1;
1835
1836 /* -------------------------------------------------------- */
1837 /* load the element matrix */
1838 /* -------------------------------------------------------- */
1839
1840 elemMat = new double[matDim*matDim];
1841 for ( i = 0; i < matDim; i++ )
1842 for ( j = 0; j < matDim; j++ )
1843 elemMat[i+j*matDim] = inMat[i][j];
1844 fedata->loadElemMatrix(elemID, matDim, elemMat);
1845 delete [] elemMat;
1846 return 0;
1847
1848 #else
1849 (void) object;
1850 (void) elemID;
1851 (void) matDim;
1852 (void) inMat;
1853 return 1;
1854 #endif
1855 }
1856
1857 /****************************************************************************/
1858 /* HYPRE_LSI_MLIFEDataWriteToFile */
1859 /*--------------------------------------------------------------------------*/
1860
1861 extern "C"
HYPRE_LSI_MLIFEDataWriteToFile(void * object,char * filename)1862 int HYPRE_LSI_MLIFEDataWriteToFile( void *object, char *filename )
1863 {
1864 #ifdef HAVE_MLI
1865 MLI_FEData *fedata;
1866 HYPRE_MLI_FEData *hypre_fedata = (HYPRE_MLI_FEData *) object;
1867
1868 /* -------------------------------------------------------- */
1869 /* error checking */
1870 /* -------------------------------------------------------- */
1871
1872 if ( hypre_fedata == NULL ) return 1;
1873 fedata = (MLI_FEData *) hypre_fedata->fedata_;
1874 if ( fedata == NULL ) return 1;
1875
1876 /* -------------------------------------------------------- */
1877 /* write to file */
1878 /* -------------------------------------------------------- */
1879
1880 fedata->writeToFile( filename );
1881 return 0;
1882 #else
1883 (void) object;
1884 (void) filename;
1885 return 1;
1886 #endif
1887 }
1888
1889 /****************************************************************************/
1890 /* HYPRE_LSI_MLISFEICreate */
1891 /*--------------------------------------------------------------------------*/
1892
1893 extern "C"
HYPRE_LSI_MLISFEICreate(MPI_Comm mpiComm)1894 void *HYPRE_LSI_MLISFEICreate( MPI_Comm mpiComm )
1895 {
1896 #ifdef HAVE_MLI
1897 HYPRE_MLI_SFEI *hypre_sfei;
1898 hypre_sfei = hypre_TAlloc(HYPRE_MLI_SFEI, 1, HYPRE_MEMORY_HOST);
1899 hypre_sfei->comm_ = mpiComm;
1900 hypre_sfei->sfei_ = new MLI_SFEI(mpiComm);;
1901 hypre_sfei->sfeiOwn_ = 1;
1902 return ((void *) hypre_sfei);
1903 #else
1904 return NULL;
1905 #endif
1906 }
1907
1908 /****************************************************************************/
1909 /* HYPRE_LSI_MLISFEIDestroy */
1910 /*--------------------------------------------------------------------------*/
1911
1912 extern "C"
HYPRE_LSI_MLISFEIDestroy(void * object)1913 int HYPRE_LSI_MLISFEIDestroy( void *object )
1914 {
1915 #ifdef HAVE_MLI
1916 HYPRE_MLI_SFEI *hypre_sfei = (HYPRE_MLI_SFEI *) object;
1917 if ( hypre_sfei == NULL ) return 1;
1918 if ( hypre_sfei->sfeiOwn_ && hypre_sfei->sfei_ != NULL )
1919 delete hypre_sfei->sfei_;
1920 hypre_sfei->sfei_ = NULL;
1921 free( hypre_sfei );
1922 return 0;
1923 #else
1924 return 1;
1925 #endif
1926 }
1927
1928 /****************************************************************************/
1929 /* HYPRE_LSI_MLISFEILoadElemMatrices */
1930 /*--------------------------------------------------------------------------*/
1931
1932 extern "C"
HYPRE_LSI_MLISFEILoadElemMatrices(void * object,int elemBlk,int nElems,int * elemIDs,double *** inMat,int elemNNodes,int ** nodeLists)1933 int HYPRE_LSI_MLISFEILoadElemMatrices(void *object, int elemBlk, int nElems,
1934 int *elemIDs, double ***inMat, int elemNNodes, int **nodeLists)
1935 {
1936 #ifdef HAVE_MLI
1937 HYPRE_MLI_SFEI *hypre_sfei = (HYPRE_MLI_SFEI *) object;
1938 MLI_SFEI *sfei;
1939
1940 /* -------------------------------------------------------- */
1941 /* error checking */
1942 /* -------------------------------------------------------- */
1943
1944 if ( hypre_sfei == NULL ) return 1;
1945 sfei = (MLI_SFEI *) hypre_sfei->sfei_;
1946 if ( sfei == NULL ) return 1;
1947
1948 /* -------------------------------------------------------- */
1949 /* load the element matrix */
1950 /* -------------------------------------------------------- */
1951
1952 sfei->loadElemBlock(elemBlk,nElems,elemIDs,inMat,elemNNodes,nodeLists);
1953 return 0;
1954 #else
1955 (void) object;
1956 (void) elemBlk;
1957 (void) nElems;
1958 (void) elemIDs;
1959 (void) inMat;
1960 (void) elemNNodes;
1961 (void) nodeLists;
1962 return 1;
1963 #endif
1964 }
1965
1966 /****************************************************************************/
1967 /* HYPRE_LSI_MLISFEIAddNumElems */
1968 /*--------------------------------------------------------------------------*/
1969
1970 extern "C"
HYPRE_LSI_MLISFEIAddNumElems(void * object,int elemBlk,int nElems,int elemNNodes)1971 int HYPRE_LSI_MLISFEIAddNumElems(void *object, int elemBlk, int nElems,
1972 int elemNNodes)
1973 {
1974 #ifdef HAVE_MLI
1975 HYPRE_MLI_SFEI *hypre_sfei = (HYPRE_MLI_SFEI *) object;
1976 MLI_SFEI *sfei;
1977
1978 /* -------------------------------------------------------- */
1979 /* error checking */
1980 /* -------------------------------------------------------- */
1981
1982 if ( hypre_sfei == NULL ) return 1;
1983 sfei = (MLI_SFEI *) hypre_sfei->sfei_;
1984 if ( sfei == NULL ) return 1;
1985
1986 /* -------------------------------------------------------- */
1987 /* send information to sfei object */
1988 /* -------------------------------------------------------- */
1989
1990 sfei->addNumElems(elemBlk,nElems,elemNNodes);
1991 return 0;
1992 #else
1993 (void) object;
1994 (void) elemBlk;
1995 (void) nElems;
1996 (void) elemNNodes;
1997 return 1;
1998 #endif
1999 }
2000
2001