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 // This file holds the other functions for HYPRE_LinSysCore
10 //***************************************************************************
11
12 //***************************************************************************
13 // include files
14 //---------------------------------------------------------------------------
15
16 //---------------------------------------------------------------------------
17 // system include files
18 //---------------------------------------------------------------------------
19
20 #include <stdlib.h>
21 #include <string.h>
22 #include <stdio.h>
23 #include <math.h>
24
25 #if 0 /* RDF: Not sure this is really needed */
26 #ifdef WIN32
27 #define strcmp _stricmp
28 #endif
29 #endif
30
31 //#define HAVE_SYSPDE
32
33 //---------------------------------------------------------------------------
34 // HYPRE include files
35 //---------------------------------------------------------------------------
36
37 #include "HYPRE.h"
38 #include "utilities/_hypre_utilities.h"
39 #include "IJ_mv/HYPRE_IJ_mv.h"
40 #include "parcsr_mv/HYPRE_parcsr_mv.h"
41 #include "parcsr_mv/_hypre_parcsr_mv.h"
42 #include "parcsr_ls/HYPRE_parcsr_ls.h"
43 #include "HYPRE_parcsr_bicgstabl.h"
44 #include "HYPRE_parcsr_lsicg.h"
45 #include "HYPRE_parcsr_TFQmr.h"
46 #include "HYPRE_parcsr_bicgs.h"
47 #include "HYPRE_parcsr_symqmr.h"
48 #include "HYPRE_parcsr_fgmres.h"
49 #include "HYPRE_LinSysCore.h"
50 #include "parcsr_mv/_hypre_parcsr_mv.h"
51 #include "HYPRE_LSI_schwarz.h"
52 #include "HYPRE_LSI_ddilut.h"
53 #include "HYPRE_LSI_ddict.h"
54 #include "HYPRE_LSI_poly.h"
55 #include "HYPRE_LSI_block.h"
56 #include "HYPRE_LSI_Uzawa_c.h"
57 #include "HYPRE_LSI_Dsuperlu.h"
58 #include "HYPRE_MLMaxwell.h"
59 #include "HYPRE_FEI.h"
60 //---------------------------------------------------------------------------
61 // SUPERLU include files
62 //---------------------------------------------------------------------------
63
64 #ifdef HAVE_SUPERLU
65 #include "slu_ddefs.h"
66 #include "slu_util.h"
67 #endif
68
69 //---------------------------------------------------------------------------
70 // MLI include files
71 //---------------------------------------------------------------------------
72
73 #ifdef HAVE_MLI
74 #include "HYPRE_LSI_mli.h"
75 #endif
76
77 //***************************************************************************
78 // external functions needed here
79 //---------------------------------------------------------------------------
80
81 extern "C" {
82
83 /*-------------------------------------------------------------------------*
84 * ML functions
85 *-------------------------------------------------------------------------*/
86
87 #ifdef HAVE_ML
88 int HYPRE_LSI_MLCreate( MPI_Comm, HYPRE_Solver *);
89 int HYPRE_LSI_MLDestroy( HYPRE_Solver );
90 int HYPRE_LSI_MLSetup( HYPRE_Solver, HYPRE_ParCSRMatrix,
91 HYPRE_ParVector, HYPRE_ParVector );
92 int HYPRE_LSI_MLSolve( HYPRE_Solver, HYPRE_ParCSRMatrix,
93 HYPRE_ParVector, HYPRE_ParVector );
94 int HYPRE_LSI_MLSetStrongThreshold( HYPRE_Solver, double );
95 int HYPRE_LSI_MLSetNumPreSmoothings( HYPRE_Solver, int );
96 int HYPRE_LSI_MLSetNumPostSmoothings( HYPRE_Solver, int );
97 int HYPRE_LSI_MLSetPreSmoother( HYPRE_Solver, int );
98 int HYPRE_LSI_MLSetPostSmoother( HYPRE_Solver, int );
99 int HYPRE_LSI_MLSetDampingFactor( HYPRE_Solver, double );
100 int HYPRE_LSI_MLSetMethod( HYPRE_Solver, int );
101 int HYPRE_LSI_MLSetCoarsenScheme( HYPRE_Solver , int );
102 int HYPRE_LSI_MLSetCoarseSolver( HYPRE_Solver, int );
103 int HYPRE_LSI_MLSetNumPDEs( HYPRE_Solver, int );
104 #endif
105
106 /*-------------------------------------------------------------------------*
107 * MLMaxwell functions
108 *-------------------------------------------------------------------------*/
109
110 #ifdef HAVE_MLMAXWELL
111 int HYPRE_LSI_MLMaxwellCreate(MPI_Comm, HYPRE_Solver *);
112 int HYPRE_LSI_MLMaxwellDestroy(HYPRE_Solver);
113 int HYPRE_LSI_MLMaxwellSetup(HYPRE_Solver, HYPRE_ParCSRMatrix,
114 HYPRE_ParVector, HYPRE_ParVector);
115 int HYPRE_LSI_MLMaxwellSolve(HYPRE_Solver, HYPRE_ParCSRMatrix,
116 HYPRE_ParVector, HYPRE_ParVector);
117 int HYPRE_LSI_MLMaxwellSetGMatrix(HYPRE_Solver, HYPRE_ParCSRMatrix);
118 int HYPRE_LSI_MLMaxwellSetANNMatrix(HYPRE_Solver, HYPRE_ParCSRMatrix);
119 int HYPRE_LSI_MLMaxwellSetStrongThreshold(HYPRE_Solver, double);
120 #endif
121
122 /*-------------------------------------------------------------------------*
123 * other functions
124 *-------------------------------------------------------------------------*/
125
126 void hypre_qsort1(int *, double *, int, int);
HYPRE_DummyFunction(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector)127 int HYPRE_DummyFunction(HYPRE_Solver, HYPRE_ParCSRMatrix,
128 HYPRE_ParVector, HYPRE_ParVector) {return 0;}
129 int HYPRE_LSI_Search(int*, int, int);
130 int HYPRE_LSI_SolveIdentity(HYPRE_Solver, HYPRE_ParCSRMatrix,
131 HYPRE_ParVector , HYPRE_ParVector);
132 int HYPRE_LSI_GetParCSRMatrix(HYPRE_IJMatrix,int nrows,int nnz,int*,
133 int*,double*);
134
135 /*-------------------------------------------------------------------------*
136 * Y12 functions (obsolete)
137 *-------------------------------------------------------------------------*/
138
139 #ifdef Y12M
140 void y12maf_(int*,int*,double*,int*,int*,int*,int*,double*,
141 int*,int*, double*,int*,double*,int*);
142 #endif
143
144 }
145
146 //***************************************************************************
147 // this function takes parameters for setting internal things like solver
148 // and preconditioner choice, etc.
149 //---------------------------------------------------------------------------
150
parameters(int numParams,char ** params)151 int HYPRE_LinSysCore::parameters(int numParams, char **params)
152 {
153 int i, k, nsweeps, rtype, olevel, reuse=0, precon_override=0;
154 int solver_override=0, solver_index=-1, precon_index=-1;
155 double weight, dtemp;
156 char param[256], param1[256], param2[80], param3[80];
157 int recognized;
158
159 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
160 {
161 printf("%4d : HYPRE_LSC::entering parameters function.\n",mypid_);
162 if ( mypid_ == 0 )
163 {
164 printf("%4d : HYPRE_LSC::parameters - numParams = %d\n", mypid_,
165 numParams);
166 for ( i = 0; i < numParams; i++ )
167 printf(" param %d = %s \n", i, params[i]);
168 }
169 }
170 if ( numParams <= 0 ) return (0);
171
172 //-------------------------------------------------------------------
173 // process the solver and preconditioner selection first
174 //-------------------------------------------------------------------
175
176 for ( i = 0; i < numParams; i++ )
177 {
178 sscanf(params[i],"%s", param1);
179 strcpy(param3, "invalid");
180 if ( !strcmp(param1, "solver") && (!solver_override) )
181 {
182 sscanf(params[i],"%s %s %s", param, param2, param3);
183 solver_index = i;
184 if (!strcmp(param3, "override")) solver_override = 1;
185 }
186 if ( !strcmp(param1, "preconditioner") && (!precon_override) )
187 {
188 sscanf(params[i],"%s %s %s", param, param2, param3);
189 if ( strcmp(param2, "reuse") )
190 {
191 precon_index = i;
192 if (!strcmp(param3, "override")) precon_override = 1;
193 }
194 }
195 }
196
197 //-------------------------------------------------------------------
198 // select solver : cg, gmres, superlu, superlux, y12m
199 //-------------------------------------------------------------------
200
201 if ( solver_index >= 0 )
202 {
203 sscanf(params[solver_index],"%s %s", param, HYSolverName_);
204 selectSolver(HYSolverName_);
205 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 && mypid_ == 0 )
206 printf(" HYPRE_LSC::parameters solver = %s\n",HYSolverName_);
207 }
208
209 //-------------------------------------------------------------------
210 // select preconditioner : diagonal, pilut, boomeramg, parasails
211 //-------------------------------------------------------------------
212
213 if ( precon_index >= 0 )
214 {
215 sscanf(params[precon_index],"%s %s", param, param1);
216 selectPreconditioner(param1);
217 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 && mypid_ == 0 )
218 printf(" HYPRE_LSC::parameters preconditioner = %s\n",
219 param1);
220 }
221
222 //-------------------------------------------------------------------
223 // parse all parameters
224 //-------------------------------------------------------------------
225
226 for ( i = 0; i < numParams; i++ )
227 {
228 sscanf(params[i],"%s", param1);
229
230 //----------------------------------------------------------------
231 // help menu
232 //----------------------------------------------------------------
233
234 recognized = 1;
235 if ( !strcmp(param1, "help") )
236 {
237 printf("%4d : HYPRE_LinSysCore::parameters - available ones : \n",
238 mypid_);
239 printf(" - outputLevel <d> \n");
240 printf(" - optimizeMemory \n");
241 printf(" - setDebug <slideReduction1,amgDebug,printFEInfo>\n");
242 printf(" - haveFEData <0,1>\n");
243 printf(" - schurReduction\n");
244 printf(" - slideReduction, slideReduction2, slideReduction3\n");
245 printf(" - slideReductionMinNorm <f>\n");
246 printf(" - matrixPartition\n");
247 printf(" - AConjugateProjection <dsize>\n");
248 printf(" - minResProjection <dsize>\n");
249 printf(" - solver <cg,gmres,bicgstab,boomeramg,superlux,..>\n");
250 printf(" - maxIterations <d>\n");
251 printf(" - tolerance <f>\n");
252 printf(" - gmresDim <d>\n");
253 printf(" - stopCrit <absolute,relative>\n");
254 printf(" - pcgRecomputeResiudal\n");
255 printf(" - preconditioner <identity,diagonal,pilut,parasails,\n");
256 printf(" - boomeramg,ddilut,schwarz,ddict,poly,euclid,...\n");
257 printf(" - blockP,ml,mli,reuse,parasails_reuse> <override>\n");
258 printf(" - pilutFillin or pilutRowSize <d>\n");
259 printf(" - pilutDropTol <f>\n");
260 printf(" - ddilutFillin <f>\n");
261 printf(" - ddilutDropTol <f> (f*sparsity(A))\n");
262 printf(" - ddilutOverlap\n");
263 printf(" - ddilutReorder\n");
264 printf(" - ddictFillin <f>\n");
265 printf(" - ddictDropTol <f> (f*sparsity(A))\n");
266 printf(" - schwarzNBlocks <d>\n");
267 printf(" - schwarzBlockSize <d>\n");
268 printf(" - polyorder <d>\n");
269 printf(" - superluOrdering <natural,mmd>\n");
270 printf(" - superluScale <y,n>\n");
271 printf(" - amgMaxLevels <d>\n");
272 printf(" - amgCoarsenType <cljp,falgout,ruge,ruge3c,pmis,hmis>\n");
273 printf(" - amgMeasureType <global,local>\n");
274 printf(" - amgRelaxType <jacobi,gsFast,hybrid,hybridsym,l1gs>\n");
275 printf(" - amgNumSweeps <d>\n");
276 printf(" - amgRelaxWeight <f>\n");
277 printf(" - amgRelaxOmega <f>\n");
278 printf(" - amgStrongThreshold <f>\n");
279 printf(" - amgSystemSize <d>\n");
280 printf(" - amgMaxIterations <d>\n");
281 printf(" - amgSmoothType <d>\n");
282 printf(" - amgSmoothNumLevels <d>\n");
283 printf(" - amgSmoothNumSweeps <d>\n");
284 printf(" - amgSchwarzRelaxWt <d>\n");
285 printf(" - amgSchwarzVariant <d>\n");
286 printf(" - amgSchwarzOverlap <d>\n");
287 printf(" - amgSchwarzDomainType <d>\n");
288 printf(" - amgUseGSMG\n");
289 printf(" - amgGSMGNumSamples\n");
290 printf(" - amgAggLevels <d>\n");
291 printf(" - amgInterpType <d>\n");
292 printf(" - amgPmax <d>\n");
293 printf(" - parasailsThreshold <f>\n");
294 printf(" - parasailsNlevels <d>\n");
295 printf(" - parasailsFilter <f>\n");
296 printf(" - parasailsLoadbal <f>\n");
297 printf(" - parasailsSymmetric\n");
298 printf(" - parasailsUnSymmetric\n");
299 printf(" - parasailsReuse <0,1>\n");
300 printf(" - euclidNlevels <d>\n");
301 printf(" - euclidThreshold <f>\n");
302 printf(" - blockP help (to get blockP options) \n");
303 printf(" - amsNumPDEs <d>\n");
304 printf(" - MLI help (to get MLI options) \n");
305 printf(" - syspdeNVars <d> \n");
306 #ifdef HAVE_ML
307 printf(" - mlNumSweeps <d>\n");
308 printf(" - mlRelaxType <jacobi,sgs,vbsgs>\n");
309 printf(" - mlRelaxWeight <f>\n");
310 printf(" - mlStrongThreshold <f>\n");
311 printf(" - mlMethod <amg>\n");
312 printf(" - mlCoarseSolver <superlu,aggregation,GS>\n");
313 printf(" - mlCoarsenScheme <uncoupled,coupled,mis>\n");
314 printf(" - mlNumPDEs <d>\n");
315 #endif
316 }
317
318 //----------------------------------------------------------------
319 // output level
320 //----------------------------------------------------------------
321
322 else if ( !strcmp(param1, "outputLevel") )
323 {
324 sscanf(params[i],"%s %d", param, &olevel);
325 if ( olevel < 0 ) olevel = 0;
326 if ( olevel > 7 ) olevel = 7;
327 HYOutputLevel_ = ( HYOutputLevel_ & HYFEI_HIGHMASK ) + olevel;
328 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
329 printf(" HYPRE_LSC::parameters outputLevel = %d\n",
330 HYOutputLevel_);
331 }
332
333 //----------------------------------------------------------------
334 // turn on memory optimizer
335 //----------------------------------------------------------------
336
337 else if ( !strcmp(param1, "optimizeMemory") )
338 {
339 memOptimizerFlag_ = 1;
340 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
341 printf(" HYPRE_LSC::parameters optimizeMemory on\n");
342 }
343
344 //----------------------------------------------------------------
345 // put no boundary condition on the matrix (for diagnostics only)
346 //----------------------------------------------------------------
347
348 else if ( !strcmp(param1, "imposeNoBC") )
349 {
350 HYOutputLevel_ |= HYFEI_IMPOSENOBC;
351 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
352 printf(" HYPRE_LSC::parameters imposeNoBC on\n");
353 }
354
355 //----------------------------------------------------------------
356 // turn on multiple right hand side
357 //----------------------------------------------------------------
358
359 else if ( !strcmp(param1, "mRHS") )
360 {
361 mRHSFlag_ = 1;
362 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
363 printf(" HYPRE_LSC::parameters multiple rhs on\n");
364 }
365
366 //----------------------------------------------------------------
367 // set matrix trunction threshold
368 //----------------------------------------------------------------
369
370 else if ( !strcmp(param1, "setTruncationThreshold") )
371 {
372 sscanf(params[i],"%s %lg", param, &truncThresh_);
373 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
374 printf(" HYPRE_LSC::parameters truncThresh = %e\n",
375 truncThresh_);
376 }
377
378 //----------------------------------------------------------------
379 // turn on fetching diagonal
380 //----------------------------------------------------------------
381
382 else if ( !strcmp(param1, "set_mixed_diag") )
383 {
384 FEI_mixedDiagFlag_ = 1;
385 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
386 printf(" HYPRE_LSC::parameters set mixed diagonal\n");
387 }
388
389 //----------------------------------------------------------------
390 // scale the matrix
391 //----------------------------------------------------------------
392
393 else if ( !strcmp(param1, "slideReductionScaleMatrix") )
394 {
395 slideReductionScaleMatrix_ = 1;
396 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
397 printf(" HYPRE_LSC::parameters slideReduction scaleMat\n");
398 }
399
400 //----------------------------------------------------------------
401 // special output level
402 //----------------------------------------------------------------
403
404 else if ( !strcmp(param1, "setDebug") )
405 {
406 sscanf(params[i],"%s %s", param, param2);
407 if (!strcmp(param2, "slideReduction1"))
408 HYOutputLevel_ |= HYFEI_SLIDEREDUCE1;
409 else if (!strcmp(param2, "slideReduction2"))
410 HYOutputLevel_ |= HYFEI_SLIDEREDUCE2;
411 else if (!strcmp(param2, "slideReduction3"))
412 HYOutputLevel_ |= HYFEI_SLIDEREDUCE3;
413 else if (!strcmp(param2, "schurReduction1"))
414 HYOutputLevel_ |= HYFEI_SCHURREDUCE1;
415 else if (!strcmp(param2, "schurReduction2"))
416 HYOutputLevel_ |= HYFEI_SCHURREDUCE2;
417 else if (!strcmp(param2, "schurReduction3"))
418 HYOutputLevel_ |= HYFEI_SCHURREDUCE3;
419 else if (!strcmp(param2, "amgDebug"))
420 HYOutputLevel_ |= HYFEI_AMGDEBUG;
421 else if (!strcmp(param2, "printMat"))
422 HYOutputLevel_ |= HYFEI_PRINTMAT;
423 else if (!strcmp(param2, "printSol"))
424 HYOutputLevel_ |= HYFEI_PRINTSOL;
425 else if (!strcmp(param2, "printReducedMat"))
426 HYOutputLevel_ |= HYFEI_PRINTREDMAT;
427 else if (!strcmp(param2, "printParCSRMat"))
428 HYOutputLevel_ |= HYFEI_PRINTPARCSRMAT;
429 else if (!strcmp(param2, "printFEInfo"))
430 HYOutputLevel_ |= HYFEI_PRINTFEINFO;
431 else if (!strcmp(param2, "ddilut"))
432 HYOutputLevel_ |= HYFEI_DDILUT;
433 else if (!strcmp(param2, "stopAfterPrint"))
434 HYOutputLevel_ |= HYFEI_STOPAFTERPRINT;
435 else if (!strcmp(param2, "off"))
436 HYOutputLevel_ = 0;
437 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
438 printf(" HYPRE_LSC::parameters setDebug %s.\n", param2);
439 }
440
441 //----------------------------------------------------------------
442 // turn on MLI's FEData module
443 //----------------------------------------------------------------
444
445 else if ( !strcmp(param1, "haveFEData") )
446 {
447 sscanf(params[i],"%s %d", param, &haveFEData_);
448 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
449 printf(" HYPRE_LSC::parameters haveFEData = %d\n",
450 haveFEData_);
451 }
452 else if ( !strcmp(param1, "haveSFEI") )
453 {
454 haveFEData_ = 2;
455 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
456 printf(" HYPRE_LSC::parameters haveSFEI\n");
457 }
458
459 //----------------------------------------------------------------
460 // turn on normal equation option
461 //----------------------------------------------------------------
462
463 else if ( !strcmp(param1, "normalEquation") )
464 {
465 normalEqnFlag_ = 1;
466 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
467 printf(" HYPRE_LSC::parameters - normal equation on.\n");
468 }
469
470 //----------------------------------------------------------------
471 // perform Schur complement reduction
472 //----------------------------------------------------------------
473
474 else if ( !strcmp(param1, "schurReduction") )
475 {
476 schurReduction_ = 1;
477 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
478 printf(" HYPRE_LSC::parameters - schur reduction.\n");
479 }
480
481 //----------------------------------------------------------------
482 // perform slide reduction
483 //----------------------------------------------------------------
484
485 else if ( !strcmp(param1, "slideReduction") )
486 {
487 slideReduction_ = 1;
488 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
489 printf(" HYPRE_LSC::parameters - slideReduction.\n");
490 }
491 else if ( !strcmp(param1, "slideReduction2") )
492 {
493 slideReduction_ = 2;
494 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
495 printf(" HYPRE_LSC::parameters - slideReduction2.\n");
496 }
497 else if ( !strcmp(param1, "slideReduction3") )
498 {
499 slideReduction_ = 3;
500 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
501 printf(" HYPRE_LSC::parameters - slideReduction3.\n");
502 }
503 else if ( !strcmp(param1, "slideReduction4") )
504 {
505 slideReduction_ = 4;
506 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
507 printf(" HYPRE_LSC::parameters - slideReduction4.\n");
508 }
509 else if ( !strcmp(param1, "slideReductionMinNorm") )
510 {
511 sscanf(params[i],"%s %lg", param, &slideReductionMinNorm_);
512 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
513 printf(" HYPRE_LSC::parameters - slideReductionMinNorm.\n");
514 }
515 else if ( !strcmp(param1, "matrixPartition") )
516 {
517 matrixPartition_ = 1;
518 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
519 printf(" HYPRE_LSC::parameters - matrixPartition.\n");
520 }
521
522 //----------------------------------------------------------------
523 // perform A-conjugate projection
524 //----------------------------------------------------------------
525
526 else if ( !strcmp(param1, "AConjugateProjection") )
527 {
528 if ( HYpbs_ != NULL )
529 {
530 for ( k = 0; k <= projectSize_; k++ )
531 if ( HYpbs_[k] != NULL ) HYPRE_IJVectorDestroy(HYpbs_[k]);
532 delete [] HYpbs_;
533 HYpbs_ = NULL;
534 }
535 if ( HYpxs_ != NULL )
536 {
537 for ( k = 0; k <= projectSize_; k++ )
538 if ( HYpxs_[k] != NULL ) HYPRE_IJVectorDestroy(HYpxs_[k]);
539 delete [] HYpxs_;
540 HYpxs_ = NULL;
541 }
542 sscanf(params[i],"%s %d", param, &k);
543 if ( k > 0 && k < 100 ) projectSize_ = k; else projectSize_ = 10;
544 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
545 printf(" HYPRE_LSC::parameters AConjugateProjection = %d\n",
546 projectSize_);
547 projectionScheme_ = 1;
548 }
549
550 //----------------------------------------------------------------
551 // perform minimal residual projection
552 //----------------------------------------------------------------
553
554 else if ( !strcmp(param1, "minResProjection") )
555 {
556 if ( HYpbs_ != NULL )
557 {
558 for ( k = 0; k <= projectSize_; k++ )
559 if ( HYpbs_[k] != NULL ) HYPRE_IJVectorDestroy(HYpbs_[k]);
560 delete [] HYpbs_;
561 HYpbs_ = NULL;
562 }
563 if ( HYpxs_ != NULL )
564 {
565 for ( k = 0; k <= projectSize_; k++ )
566 if ( HYpxs_[k] != NULL ) HYPRE_IJVectorDestroy(HYpxs_[k]);
567 delete [] HYpxs_;
568 HYpxs_ = NULL;
569 }
570 sscanf(params[i],"%s %d", param, &k);
571 if ( k > 0 && k < 100 ) projectSize_ = k; else projectSize_ = 10;
572 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
573 printf(" HYPRE_LSC::parameters minResProjection = %d\n",
574 projectSize_);
575 projectionScheme_ = 2;
576 }
577
578 //----------------------------------------------------------------
579 // which solver to pick : cg, gmres, superlu, superlux, y12m
580 //----------------------------------------------------------------
581
582 else if ( !strcmp(param1, "solver") )
583 {
584 solver_override = 0;
585 }
586
587 //----------------------------------------------------------------
588 // for GMRES, the restart size
589 //----------------------------------------------------------------
590
591 else if ( !strcmp(param1, "gmresDim") )
592 {
593 sscanf(params[i],"%s %d", param, &gmresDim_);
594 if ( gmresDim_ < 1 ) gmresDim_ = 100;
595 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
596 printf(" HYPRE_LSC::parameters gmresDim = %d\n",
597 gmresDim_);
598 }
599
600 //----------------------------------------------------------------
601 // for FGMRES, tell it to update the convergence criterion of its
602 // preconditioner, if it is Block preconditioning
603 //----------------------------------------------------------------
604
605 else if ( !strcmp(param1, "fgmresUpdateTol") )
606 {
607 fgmresUpdateTol_ = 1;
608 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
609 printf(" HYPRE_LSC::parameters fgmresUpdateTol on\n");
610 }
611
612 //----------------------------------------------------------------
613 // for GMRES, the convergence criterion
614 //----------------------------------------------------------------
615
616 else if ( !strcmp(param1, "gmresStopCrit") )
617 {
618 sscanf(params[i],"%s %s", param, param2);
619 if ( !strcmp(param2, "absolute" ) ) normAbsRel_ = 1;
620 else if ( !strcmp(param2, "relative" ) ) normAbsRel_ = 0;
621 else normAbsRel_ = 0;
622 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
623 printf(" HYPRE_LSC::parameters gmresStopCrit = %s\n",
624 param2);
625 }
626
627 else if ( !strcmp(param1, "stopCrit") )
628 {
629 sscanf(params[i],"%s %s", param, param2);
630 if ( !strcmp(param2, "absolute") ) normAbsRel_ = 1;
631 else if ( !strcmp(param2, "relative") ) normAbsRel_ = 0;
632 else normAbsRel_ = 0;
633
634 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
635 printf(" HYPRE_LSC::parameters stopCrit = %s\n",
636 param2);
637 }
638
639 //----------------------------------------------------------------
640 // for PCG only
641 //----------------------------------------------------------------
642
643 else if ( !strcmp(param1, "pcgRecomputeResidual") )
644 {
645 pcgRecomputeRes_ = 1;
646 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
647 printf(" HYPRE_LSC::parameters pcgRecomputeResidual\n");
648 }
649
650 //----------------------------------------------------------------
651 // preconditioner reuse
652 //----------------------------------------------------------------
653
654 else if ( !strcmp(param1, "precond_reuse") )
655 {
656 sscanf(params[i],"%s %s", param, param2);
657 if ( !strcmp(param2, "on") ) HYPreconReuse_ = reuse = 1;
658 else HYPreconReuse_ = 0;
659 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
660 printf(" HYPRE_LSC::parameters precond_reuse = %s\n",
661 param2);
662 }
663
664 //----------------------------------------------------------------
665 // which preconditioner : diagonal, pilut, boomeramg, parasails
666 //----------------------------------------------------------------
667
668 else if ( !strcmp(param1, "preconditioner") )
669 {
670 sscanf(params[i],"%s %s", param, param2);
671 if ( !strcmp(param2, "reuse") ) HYPreconReuse_ = reuse = 1;
672 else if ( !strcmp(param2, "parasails_reuse") ) parasailsReuse_ = 1;
673 }
674
675 //----------------------------------------------------------------
676 // maximum number of iterations for pcg or gmres
677 //----------------------------------------------------------------
678
679 else if ( !strcmp(param1, "maxIterations") )
680 {
681 sscanf(params[i],"%s %d", param, &maxIterations_);
682 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
683 printf(" HYPRE_LSC::parameters maxIterations = %d\n",
684 maxIterations_);
685 }
686
687 //----------------------------------------------------------------
688 // tolerance as termination criterion
689 //----------------------------------------------------------------
690
691 else if ( !strcmp(param1, "tolerance") )
692 {
693 sscanf(params[i],"%s %lg", param, &tolerance_);
694 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
695 printf(" HYPRE_LSC::parameters tolerance = %e\n",
696 tolerance_);
697 }
698
699 //----------------------------------------------------------------
700 // pilut preconditioner : max number of nonzeros to keep per row
701 //----------------------------------------------------------------
702
703 else if ( !strcmp(param1, "pilutFillin") )
704 {
705 sscanf(params[i],"%s %d", param, &pilutFillin_);
706 if ( pilutFillin_ < 1 ) pilutFillin_ = 50;
707 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
708 printf(" HYPRE_LSC::parameters pilutFillin_ = %d\n",
709 pilutFillin_);
710 }
711 else if ( !strcmp(param1, "pilutRowSize") )
712 {
713 sscanf(params[i],"%s %d", param, &pilutFillin_);
714 if ( pilutFillin_ < 1 ) pilutFillin_ = 50;
715 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
716 printf(" HYPRE_LSC::parameters pilutFillin_ = %d\n",
717 pilutFillin_);
718 }
719
720 //----------------------------------------------------------------
721 // pilut preconditioner : threshold to drop small nonzeros
722 //----------------------------------------------------------------
723
724 else if ( !strcmp(param1, "pilutDropTol") )
725 {
726 sscanf(params[i],"%s %lg", param, &pilutDropTol_);
727 if (pilutDropTol_<0.0 || pilutDropTol_ >=1.0) pilutDropTol_ = 0.0;
728 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
729 printf(" HYPRE_LSC::parameters pilutDropTol = %e\n",
730 pilutDropTol_);
731 }
732
733 //----------------------------------------------------------------
734 // DDILUT preconditioner : amount of fillin (0 == same as A)
735 //----------------------------------------------------------------
736
737 else if ( !strcmp(param1, "ddilutFillin") )
738 {
739 sscanf(params[i],"%s %lg", param, &ddilutFillin_);
740 if ( ddilutFillin_ < 0.0 ) ddilutFillin_ = 0.0;
741 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
742 printf(" HYPRE_LSC::parameters ddilutFillin = %e\n",
743 ddilutFillin_);
744 }
745
746 //----------------------------------------------------------------
747 // DDILUT preconditioner : threshold to drop small nonzeros
748 //----------------------------------------------------------------
749
750 else if ( !strcmp(param1, "ddilutDropTol") )
751 {
752 sscanf(params[i],"%s %lg", param, &ddilutDropTol_);
753 if (ddilutDropTol_<0.0 || ddilutDropTol_ >=1.0) ddilutDropTol_ = 0.0;
754 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
755 printf(" HYPRE_LSC::parameters ddilutDropTol = %e\n",
756 ddilutDropTol_);
757 }
758
759 //----------------------------------------------------------------
760 // DDILUT preconditioner : turn on processor overlap
761 //----------------------------------------------------------------
762
763 else if ( !strcmp(param1, "ddilutOverlap") )
764 {
765 ddilutOverlap_ = 1;
766 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
767 printf(" HYPRE_LSC::parameters ddilutOverlap = on\n");
768 }
769
770 //----------------------------------------------------------------
771 // DDILUT preconditioner : reorder based on Cuthill McKee
772 //----------------------------------------------------------------
773
774 else if ( !strcmp(param1, "ddilutReorder") )
775 {
776 ddilutReorder_ = 1;
777 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
778 printf(" HYPRE_LSC::parameters ddilutReorder = on\n");
779 }
780
781 //----------------------------------------------------------------
782 // DDICT preconditioner : amount of fillin (0 == same as A)
783 //----------------------------------------------------------------
784
785 else if ( !strcmp(param1, "ddictFillin") )
786 {
787 sscanf(params[i],"%s %lg", param, &ddictFillin_);
788 if ( ddictFillin_ < 0.0 ) ddictFillin_ = 0.0;
789 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
790 printf(" HYPRE_LSC::parameters ddictFillin = %e\n",
791 ddictFillin_);
792 }
793
794 //----------------------------------------------------------------
795 // DDICT preconditioner : threshold to drop small nonzeros
796 //----------------------------------------------------------------
797
798 else if ( !strcmp(param1, "ddictDropTol") )
799 {
800 sscanf(params[i],"%s %lg", param, &ddictDropTol_);
801 if (ddictDropTol_<0.0 || ddictDropTol_ >=1.0) ddictDropTol_ = 0.0;
802 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
803 printf(" HYPRE_LSC::parameters ddictDropTol = %e\n",
804 ddictDropTol_);
805 }
806
807 //----------------------------------------------------------------
808 // Schwarz preconditioner : Fillin
809 //----------------------------------------------------------------
810
811 else if ( !strcmp(param1, "schwarzFillin") )
812 {
813 sscanf(params[i],"%s %lg", param, &schwarzFillin_);
814 if ( schwarzFillin_ < 0.0 ) schwarzFillin_ = 0.0;
815 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
816 printf(" HYPRE_LSC::parameters schwarzFillin = %e\n",
817 schwarzFillin_);
818 }
819
820 //----------------------------------------------------------------
821 // Schwarz preconditioner : block size
822 //----------------------------------------------------------------
823
824 else if ( !strcmp(param1, "schwarzNBlocks") )
825 {
826 sscanf(params[i],"%s %d", param, &schwarzNblocks_);
827 if ( schwarzNblocks_ <= 0 ) schwarzNblocks_ = 1;
828 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
829 printf(" HYPRE_LSC::parameters schwarzNblocks = %d\n",
830 schwarzNblocks_);
831 }
832
833 //----------------------------------------------------------------
834 // Schwarz preconditioner : block size
835 //----------------------------------------------------------------
836
837 else if ( !strcmp(param1, "schwarzBlockSize") )
838 {
839 sscanf(params[i],"%s %d", param, &schwarzBlksize_);
840 if ( schwarzBlksize_ <= 0 ) schwarzBlksize_ = 1000;
841 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
842 printf(" HYPRE_LSC::parameters schwarzBlockSize = %d\n",
843 schwarzBlksize_);
844 }
845
846 //----------------------------------------------------------------
847 // Polynomial preconditioner : order
848 //----------------------------------------------------------------
849
850 else if ( !strcmp(param1, "polyOrder") )
851 {
852 sscanf(params[i],"%s %d", param, &polyOrder_);
853 if ( polyOrder_ < 0 ) polyOrder_ = 0;
854 if ( polyOrder_ > 8 ) polyOrder_ = 8;
855 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
856 printf(" HYPRE_LSC::parameters polyOrder = %d\n",
857 polyOrder_);
858 }
859
860 //----------------------------------------------------------------
861 // superlu : ordering to use (natural, mmd)
862 //----------------------------------------------------------------
863
864 else if ( !strcmp(param1, "superluOrdering") )
865 {
866 sscanf(params[i],"%s %s", param, param2);
867 if ( !strcmp(param2, "natural" ) ) superluOrdering_ = 0;
868 else if ( !strcmp(param2, "mmd") ) superluOrdering_ = 2;
869 else superluOrdering_ = 0;
870 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
871 printf(" HYPRE_LSC::parameters superluOrdering = %s\n",
872 param2);
873 }
874
875 //----------------------------------------------------------------
876 // superlu : scaling none ('N') or both col/row ('B')
877 //----------------------------------------------------------------
878
879 else if ( !strcmp(param1, "superluScale") )
880 {
881 sscanf(params[i],"%s %s", param, param2);
882 if ( !strcmp(param2, "y" ) ) superluScale_[0] = 'B';
883 else superluScale_[0] = 'N';
884 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
885 printf(" HYPRE_LSC::parameters superluScale = %s\n",
886 param2);
887 }
888
889 else recognized = 0;
890
891 //----------------------------------------------------------------
892 // amg preconditoner : coarsening type
893 //----------------------------------------------------------------
894
895 if (!recognized)
896 {
897 recognized = 1;
898 if ( !strcmp(param1, "amgMaxLevels") )
899 {
900 sscanf(params[i],"%s %d", param, &amgMaxLevels_);
901 if ( amgMaxLevels_ <= 0 ) amgMaxLevels_ = 30;
902 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
903 printf(" HYPRE_LSC::parameters amgMaxLevels = %d\n",
904 amgMaxLevels_);
905 }
906
907 //----------------------------------------------------------------
908 // amg preconditoner : coarsening type
909 //----------------------------------------------------------------
910
911 else if ( !strcmp(param1, "amgCoarsenType") )
912 {
913 sscanf(params[i],"%s %s", param, param2);
914 if ( !strcmp(param2, "cljp" ) ) amgCoarsenType_ = 0;
915 else if ( !strcmp(param2, "ruge" ) ) amgCoarsenType_ = 1;
916 else if ( !strcmp(param2, "ruge3c" ) ) amgCoarsenType_ = 4;
917 else if ( !strcmp(param2, "falgout" ) ) amgCoarsenType_ = 6;
918 else if ( !strcmp(param2, "pmis" ) ) amgCoarsenType_ = 8;
919 else if ( !strcmp(param2, "hmis" ) ) amgCoarsenType_ = 10;
920 else amgCoarsenType_ = 0;
921 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
922 printf(" HYPRE_LSC::parameters amgCoarsenType = %s\n",
923 param2);
924 }
925
926 //----------------------------------------------------------------
927 // amg preconditoner : measure
928 //----------------------------------------------------------------
929
930 else if ( !strcmp(param1, "amgMeasureType") )
931 {
932 sscanf(params[i],"%s %s", param, param2);
933 if ( !strcmp(param2, "local" ) ) amgMeasureType_ = 0;
934 else if ( !strcmp(param2, "global" ) ) amgMeasureType_ = 1;
935 else amgMeasureType_ = 0;
936 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
937 printf(" HYPRE_LSC::parameters amgCoarsenType = %s\n",
938 param2);
939 }
940
941 //----------------------------------------------------------------
942 // amg preconditoner : no of relaxation sweeps per level
943 //----------------------------------------------------------------
944
945 else if ( !strcmp(param1, "amgNumSweeps") )
946 {
947 sscanf(params[i],"%s %d", param, &nsweeps);
948 if ( nsweeps < 1 ) nsweeps = 1;
949 for ( k = 0; k < 4; k++ ) amgNumSweeps_[k] = nsweeps;
950 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
951 printf(" HYPRE_LSC::parameters amgNumSweeps = %d\n",
952 nsweeps);
953 }
954
955 //---------------------------------------------------------------
956 // amg preconditoner : which smoother to use
957 //----------------------------------------------------------------
958
959 else if ( !strcmp(param1, "amgRelaxType") )
960 {
961 sscanf(params[i],"%s %s", param, param2);
962 if ( !strcmp(param2, "jacobi" ) ) rtype = 0;
963 else if ( !strcmp(param2, "CFjacobi" ) )
964 {rtype = 0; amgGridRlxType_ = 1;}
965 else if ( !strcmp(param2, "gsSlow") ) rtype = 1;
966 else if ( !strcmp(param2, "gsFast") ) rtype = 4;
967 else if ( !strcmp(param2, "hybrid" ) ) rtype = 3;
968 else if ( !strcmp(param2, "CFhybrid" ) )
969 {rtype = 3; amgGridRlxType_ = 1;}
970 else if ( !strcmp(param2, "hybridsym" ) ) rtype = 6;
971 else if ( !strcmp(param2, "l1gs" ) ) rtype = 8;
972 else if ( !strcmp(param2, "CFhybridsym" ) )
973 {rtype = 6; amgGridRlxType_ = 1;}
974 else rtype = 4;
975 for ( k = 0; k < 3; k++ ) amgRelaxType_[k] = rtype;
976 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
977 printf(" HYPRE_LSC::parameters amgRelaxType = %s\n",
978 param2);
979 }
980
981 //---------------------------------------------------------------
982 // amg preconditoner : damping factor for Jacobi smoother
983 //---------------------------------------------------------------
984
985 else if ( !strcmp(param1, "amgRelaxWeight") )
986 {
987 sscanf(params[i],"%s %lg", param, &weight);
988 for ( k = 0; k < 25; k++ ) amgRelaxWeight_[k] = weight;
989 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
990 printf(" HYPRE_LSC::parameters amgRelaxWeight = %e\n",
991 weight);
992 }
993
994 //---------------------------------------------------------------
995 // amg preconditoner : relaxation factor for hybrid smoother
996 //---------------------------------------------------------------
997
998 else if ( !strcmp(param1, "amgRelaxOmega") )
999 {
1000 sscanf(params[i],"%s %lg", param, &weight);
1001 for ( k = 0; k < 25; k++ ) amgRelaxOmega_[k] = weight;
1002 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1003 printf(" HYPRE_LSC::parameters amgRelaxOmega = %e\n",
1004 weight);
1005 }
1006
1007 //---------------------------------------------------------------
1008 // amg preconditoner : threshold to determine strong coupling
1009 //---------------------------------------------------------------
1010
1011 else if ( !strcmp(param1, "amgStrongThreshold") )
1012 {
1013 sscanf(params[i],"%s %lg", param, &amgStrongThreshold_);
1014 if ( amgStrongThreshold_ < 0.0 || amgStrongThreshold_ > 1.0 )
1015 amgStrongThreshold_ = 0.25;
1016 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1017 printf(" HYPRE_LSC::parameters amgStrongThreshold = %e\n",
1018 amgStrongThreshold_);
1019 }
1020 else if ( !strcmp(param1, "amgStrengthThreshold") )
1021 {
1022 sscanf(params[i],"%s %lg", param, &amgStrongThreshold_);
1023 if ( amgStrongThreshold_ < 0.0 || amgStrongThreshold_ > 1.0 )
1024 amgStrongThreshold_ = 0.25;
1025 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1026 printf(" HYPRE_LSC::parameters amgStrengthThreshold = %e\n",
1027 amgStrongThreshold_);
1028 }
1029
1030 //---------------------------------------------------------------
1031 // amg preconditoner : choose system size
1032 //---------------------------------------------------------------
1033
1034 else if ( !strcmp(param1, "amgSystemSize") )
1035 {
1036 sscanf(params[i],"%s %d", param, &amgSystemSize_);
1037 if ( amgSystemSize_ <= 0 ) amgSystemSize_ = 1;
1038 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1039 printf(" HYPRE_LSC::parameters amgSystemSize = %d\n",
1040 amgSystemSize_);
1041 }
1042
1043 //---------------------------------------------------------------
1044 // amg preconditoner : choose max iterations
1045 //---------------------------------------------------------------
1046
1047 else if ( !strcmp(param1, "amgMaxIterations") )
1048 {
1049 sscanf(params[i],"%s %d", param, &amgMaxIter_);
1050 if ( amgMaxIter_ <= 0 ) amgMaxIter_ = 1;
1051 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1052 printf(" HYPRE_LSC::parameters amgMaxIter = %d\n",
1053 amgMaxIter_);
1054 }
1055
1056 //---------------------------------------------------------------
1057 // amg preconditoner : choose more complex smoother
1058 //---------------------------------------------------------------
1059
1060 else if ( !strcmp(param1, "amgSmoothType") )
1061 {
1062 sscanf(params[i],"%s %d", param, &amgSmoothType_);
1063 if ( amgSmoothType_ < 0 ) amgSmoothType_ = 0;
1064 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1065 printf(" HYPRE_LSC::parameters amgSmoothType = %d\n",
1066 amgSmoothType_);
1067 }
1068
1069 //---------------------------------------------------------------
1070 // amg preconditoner : choose no. of levels for complex smoother
1071 //---------------------------------------------------------------
1072
1073 else if ( !strcmp(param1, "amgSmoothNumLevels") )
1074 {
1075 sscanf(params[i],"%s %d", param, &amgSmoothNumLevels_);
1076 if ( amgSmoothNumLevels_ < 0 ) amgSmoothNumLevels_ = 0;
1077 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1078 printf(" HYPRE_LSC::parameters amgSmoothNumLevels = %d\n",
1079 amgSmoothNumLevels_);
1080 }
1081
1082 //---------------------------------------------------------------
1083 // amg preconditoner : choose no. of sweeps for complex smoother
1084 //---------------------------------------------------------------
1085
1086 else if ( !strcmp(param1, "amgSmoothNumSweeps") )
1087 {
1088 sscanf(params[i],"%s %d", param, &amgSmoothNumSweeps_);
1089 if ( amgSmoothNumSweeps_ < 0 ) amgSmoothNumSweeps_ = 1;
1090 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1091 printf(" HYPRE_LSC::parameters amgSmoothNumSweeps = %d\n",
1092 amgSmoothNumSweeps_);
1093 }
1094
1095 //---------------------------------------------------------------
1096 // amg preconditoner : choose relaxation weight for Schwarz smoother
1097 //---------------------------------------------------------------
1098
1099 else if ( !strcmp(param1, "amgSchwarzRelaxWt") )
1100 {
1101 sscanf(params[i],"%s %lg", param, &amgSchwarzRelaxWt_);
1102 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1103 printf(" HYPRE_LSC::parameters amgSchwarzRelaxWt = %e\n",
1104 amgSchwarzRelaxWt_);
1105 }
1106
1107 //---------------------------------------------------------------
1108 // amg preconditoner : choose Schwarz smoother variant
1109 //---------------------------------------------------------------
1110
1111 else if ( !strcmp(param1, "amgSchwarzVariant") )
1112 {
1113 sscanf(params[i],"%s %d", param, &amgSchwarzVariant_);
1114 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1115 printf(" HYPRE_LSC::parameters amgSchwarzVariant = %d\n",
1116 amgSchwarzVariant_);
1117 }
1118
1119 //---------------------------------------------------------------
1120 // amg preconditoner : choose Schwarz smoother overlap
1121 //---------------------------------------------------------------
1122
1123 else if ( !strcmp(param1, "amgSchwarzOverlap") )
1124 {
1125 sscanf(params[i],"%s %d", param, &amgSchwarzOverlap_);
1126 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1127 printf(" HYPRE_LSC::parameters amgSchwarzOverlap = %d\n",
1128 amgSchwarzOverlap_);
1129 }
1130
1131 //----------------------------------------------------------------
1132 //---------------------------------------------------------------
1133 // amg preconditoner : choose Schwarz smoother domain type
1134 //---------------------------------------------------------------
1135
1136 else if ( !strcmp(param1, "amgSchwarzDomainType") )
1137 {
1138 sscanf(params[i],"%s %d", param, &amgSchwarzDomainType_);
1139 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1140 printf(" HYPRE_LSC::parameters amgSchwarzDomainType = %d\n",
1141 amgSchwarzDomainType_);
1142 }
1143
1144 //----------------------------------------------------------------
1145 //----------------------------------------------------------------
1146 // amg preconditoner : use gsmg
1147 //----------------------------------------------------------------
1148
1149 else if ( !strcmp(param1, "amgUseGSMG") )
1150 {
1151 amgUseGSMG_ = 1;
1152 if ( amgGSMGNSamples_ == 0 ) amgGSMGNSamples_ = 5;
1153 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1154 printf(" HYPRE_LSC::parameters amgUseGSMG.\n");
1155 }
1156
1157 //----------------------------------------------------------------
1158 // amg preconditoner : levels of aggresive coarsening
1159 //----------------------------------------------------------------
1160
1161 else if ( !strcmp(param1, "amgAggLevels") )
1162 {
1163 sscanf(params[i],"%s %d", param, &amgAggLevels_);
1164 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1165 printf(" HYPRE_LSC::parameters %s = %d\n",
1166 param1, amgAggLevels_);
1167 }
1168
1169 //----------------------------------------------------------------
1170 // amg preconditoner : interpolation type
1171 //----------------------------------------------------------------
1172
1173 else if ( !strcmp(param1, "amgInterpType") )
1174 {
1175 sscanf(params[i],"%s %d", param, &amgInterpType_);
1176 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1177 printf(" HYPRE_LSC::parameters %s = %d\n",
1178 param1, amgInterpType_);
1179 }
1180
1181 //----------------------------------------------------------------
1182 // amg preconditoner : interpolation truncation
1183 //----------------------------------------------------------------
1184
1185 else if ( !strcmp(param1, "amgPmax") )
1186 {
1187 sscanf(params[i],"%s %d", param, &amgPmax_);
1188 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1189 printf(" HYPRE_LSC::parameters %s = %d\n",
1190 param1, amgPmax_);
1191 }
1192
1193 //---------------------------------------------------------------
1194 // parasails preconditoner : gsmg number of samples
1195 //---------------------------------------------------------------
1196
1197 else if ( !strcmp(param1, "amgGSMGNumSamples") )
1198 {
1199 sscanf(params[i],"%s %d", param, &amgGSMGNSamples_);
1200 if ( amgGSMGNSamples_ < 0 ) amgGSMGNSamples_ = 5;
1201 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1202 printf(" HYPRE_LSC::parameters amgGSMGNumSamples = %d\n",
1203 amgGSMGNSamples_);
1204 }
1205
1206 else recognized = 0;
1207 }
1208
1209 //---------------------------------------------------------------
1210 // parasails preconditoner : threshold ( >= 0.0 )
1211 //---------------------------------------------------------------
1212
1213 if (!recognized)
1214 {
1215 recognized = 1;
1216 if ( !strcmp(param1, "parasailsThreshold") )
1217 {
1218 sscanf(params[i],"%s %lg", param, ¶sailsThreshold_);
1219 if ( parasailsThreshold_ < 0.0 ) parasailsThreshold_ = 0.1;
1220 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1221 printf(" HYPRE_LSC::parameters parasailsThreshold = %e\n",
1222 parasailsThreshold_);
1223 }
1224
1225 //---------------------------------------------------------------
1226 // parasails preconditoner : nlevels ( >= 0)
1227 //---------------------------------------------------------------
1228
1229 else if ( !strcmp(param1, "parasailsNlevels") )
1230 {
1231 sscanf(params[i],"%s %d", param, ¶sailsNlevels_);
1232 if ( parasailsNlevels_ < 0 ) parasailsNlevels_ = 1;
1233 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1234 printf(" HYPRE_LSC::parameters parasailsNlevels = %d\n",
1235 parasailsNlevels_);
1236 }
1237
1238 //---------------------------------------------------------------
1239 // parasails preconditoner : filter
1240 //---------------------------------------------------------------
1241
1242 else if ( !strcmp(param1, "parasailsFilter") )
1243 {
1244 sscanf(params[i],"%s %lg", param, ¶sailsFilter_);
1245 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1246 printf(" HYPRE_LSC::parameters parasailsFilter = %e\n",
1247 parasailsFilter_);
1248 }
1249
1250 //---------------------------------------------------------------
1251 // parasails preconditoner : loadbal
1252 //---------------------------------------------------------------
1253
1254 else if ( !strcmp(param1, "parasailsLoadbal") )
1255 {
1256 sscanf(params[i],"%s %lg", param, ¶sailsLoadbal_);
1257 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1258 printf(" HYPRE_LSC::parameters parasailsLoadbal = %e\n",
1259 parasailsLoadbal_);
1260 }
1261
1262 //---------------------------------------------------------------
1263 // parasails preconditoner : symmetry flag (1 - symm, 0 - nonsym)
1264 //---------------------------------------------------------------
1265
1266 else if ( !strcmp(param1, "parasailsSymmetric") )
1267 {
1268 parasailsSym_ = 1;
1269 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1270 printf(" HYPRE_LSC::parameters parasailsSym = sym\n");
1271 }
1272 else if ( !strcmp(param1, "parasailsUnSymmetric") )
1273 {
1274 parasailsSym_ = 0;
1275 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1276 printf(" HYPRE_LSC::parameters parasailsSym = nonsym\n");
1277 }
1278
1279 //---------------------------------------------------------------
1280 // parasails preconditoner : reuse flag
1281 //---------------------------------------------------------------
1282
1283 else if ( !strcmp(param1, "parasailsReuse") )
1284 {
1285 sscanf(params[i],"%s %d", param, ¶sailsReuse_);
1286 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1287 printf(" HYPRE_LSC::parameters parasailsReuse = %d\n",
1288 parasailsReuse_);
1289 }
1290
1291 else recognized = 0;
1292 }
1293
1294 //---------------------------------------------------------------
1295 // Euclid preconditoner : fill-in
1296 //---------------------------------------------------------------
1297
1298 if (!recognized)
1299 {
1300 recognized = 1;
1301 if ( !strcmp(param1, "euclidNlevels") )
1302 {
1303 sscanf(params[i],"%s %d", param, &olevel);
1304 if ( olevel < 0 ) olevel = 0;
1305 sprintf( euclidargv_[1], "%d", olevel);
1306 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1307 printf(" HYPRE_LSC::parameters euclidNlevels = %d\n",
1308 olevel);
1309 }
1310
1311 //---------------------------------------------------------------
1312 // Euclid preconditoner : threshold
1313 //---------------------------------------------------------------
1314
1315 else if ( !strcmp(param1, "euclidThreshold") )
1316 {
1317 sscanf(params[i],"%s %lg", param, &dtemp);
1318 if ( dtemp < 0.0 ) dtemp = 0.0;
1319 sprintf( euclidargv_[3], "%e", dtemp);
1320 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1321 printf(" HYPRE_LSC::parameters euclidThreshold = %e\n",
1322 dtemp);
1323 }
1324
1325 //---------------------------------------------------------------
1326 // block preconditoner (hold this until this end)
1327 //---------------------------------------------------------------
1328
1329 else if ( !strcmp(param1, "blockP") )
1330 {
1331 if ( HYPreconID_ == HYBLOCK )
1332 HYPRE_LSI_BlockPrecondSetParams(HYPrecon_, params[i]);
1333 }
1334
1335 //---------------------------------------------------------------
1336 // MLI preconditoners (hold this until the end)
1337 //---------------------------------------------------------------
1338
1339 else if ( !strcmp(param1, "MLI_Hybrid") )
1340 {
1341 MLI_Hybrid_GSA_ = 1;
1342 }
1343 else if ( !strcmp(param1, "MLI_Hybrid_NSIncr") )
1344 {
1345 sscanf(params[i],"%s %d", param, &MLI_Hybrid_NSIncr_);
1346 if ( MLI_Hybrid_NSIncr_ <= 0 ) MLI_Hybrid_NSIncr_ = 1;
1347 if ( MLI_Hybrid_NSIncr_ > 10 ) MLI_Hybrid_NSIncr_ = 10;
1348 }
1349 else if ( !strcmp(param1, "MLI_Hybrid_MaxIter") )
1350 {
1351 sscanf(params[i],"%s %d", param, &MLI_Hybrid_MaxIter_);
1352 if ( MLI_Hybrid_MaxIter_ <= 0 ) MLI_Hybrid_MaxIter_ = 10;
1353 }
1354 else if ( !strcmp(param1, "MLI_Hybrid_NTrials") )
1355 {
1356 sscanf(params[i],"%s %d", param, &MLI_Hybrid_NTrials_);
1357 if ( MLI_Hybrid_NTrials_ <= 0 ) MLI_Hybrid_NTrials_ = 1;
1358 }
1359 else if ( !strcmp(param1, "MLI_Hybrid_ConvRate") )
1360 {
1361 sscanf(params[i],"%s %lg", param, &MLI_Hybrid_ConvRate_);
1362 if ( MLI_Hybrid_ConvRate_ >= 1.0 ) MLI_Hybrid_ConvRate_ = 1.0;
1363 if ( MLI_Hybrid_ConvRate_ <= 0.0 ) MLI_Hybrid_ConvRate_ = 0.0;
1364 }
1365 else if ( !strcmp(param1, "MLI") )
1366 {
1367 #ifdef HAVE_MLI
1368 if ( HYPreconID_ == HYMLI )
1369 HYPRE_LSI_MLISetParams(HYPrecon_, params[i]);
1370 #else
1371 // if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 && mypid_ == 0 )
1372 // printf(" HYPRE_LSC::MLI SetParams - MLI unavailable.\n");
1373 #endif
1374 }
1375
1376 //----------------------------------------------------------------
1377 // for Uzawa, the various parameters
1378 //----------------------------------------------------------------
1379
1380 else if ( !strcmp(param1, "Uzawa") )
1381 {
1382 if ( HYPreconID_ == HYUZAWA )
1383 HYPRE_LSI_UzawaSetParams(HYPrecon_, params[i]);
1384 }
1385
1386 else recognized = 0;
1387 }
1388
1389 //---------------------------------------------------------------
1390 // mlpack preconditoner : no of relaxation sweeps per level
1391 //---------------------------------------------------------------
1392
1393 if (!recognized)
1394 {
1395 recognized = 1;
1396 if ( !strcmp(param1, "mlNumPresweeps") )
1397 {
1398 sscanf(params[i],"%s %d", param, &nsweeps);
1399 if ( nsweeps < 1 ) nsweeps = 1;
1400 mlNumPreSweeps_ = nsweeps;
1401 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1402 printf(" HYPRE_LSC::parameters mlNumPresweeps = %d\n",
1403 nsweeps);
1404 }
1405 else if ( !strcmp(param1, "mlNumPostsweeps") )
1406 {
1407 sscanf(params[i],"%s %d", param, &nsweeps);
1408 if ( nsweeps < 1 ) nsweeps = 1;
1409 mlNumPostSweeps_ = nsweeps;
1410 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1411 printf(" HYPRE_LSC::parameters mlNumPostsweeps = %d\n",
1412 nsweeps);
1413 }
1414 else if ( !strcmp(param1, "mlNumSweeps") )
1415 {
1416 sscanf(params[i],"%s %d", param, &nsweeps);
1417 if ( nsweeps < 1 ) nsweeps = 1;
1418 mlNumPreSweeps_ = nsweeps;
1419 mlNumPostSweeps_ = nsweeps;
1420 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1421 printf(" HYPRE_LSC::parameters mlNumSweeps = %d\n",
1422 nsweeps);
1423 }
1424
1425 //---------------------------------------------------------------
1426 // mlpack preconditoner : which smoother to use
1427 //---------------------------------------------------------------
1428
1429 else if ( !strcmp(param1, "mlPresmootherType") )
1430 {
1431 sscanf(params[i],"%s %s", param, param2);
1432 rtype = 1;
1433 if ( !strcmp(param2, "jacobi" ) ) rtype = 0;
1434 else if ( !strcmp(param2, "sgs") ) rtype = 1;
1435 else if ( !strcmp(param2, "sgsseq") ) rtype = 2;
1436 else if ( !strcmp(param2, "vbjacobi")) rtype = 3;
1437 else if ( !strcmp(param2, "vbsgs") ) rtype = 4;
1438 else if ( !strcmp(param2, "vbsgsseq")) rtype = 5;
1439 else if ( !strcmp(param2, "ilut") ) rtype = 6;
1440 else if ( !strcmp(param2, "aSchwarz") ) rtype = 7;
1441 else if ( !strcmp(param2, "mSchwarz") ) rtype = 8;
1442 mlPresmootherType_ = rtype;
1443 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1444 printf(" HYPRE_LSC::parameters mlPresmootherType = %s\n",
1445 param2);
1446 }
1447 else if ( !strcmp(param1, "mlPostsmootherType") )
1448 {
1449 sscanf(params[i],"%s %s", param, param2);
1450 rtype = 1;
1451 if ( !strcmp(param2, "jacobi" ) ) rtype = 0;
1452 else if ( !strcmp(param2, "sgs") ) rtype = 1;
1453 else if ( !strcmp(param2, "sgsseq") ) rtype = 2;
1454 else if ( !strcmp(param2, "vbjacobi")) rtype = 3;
1455 else if ( !strcmp(param2, "vbsgs") ) rtype = 4;
1456 else if ( !strcmp(param2, "vbsgsseq")) rtype = 5;
1457 mlPostsmootherType_ = rtype;
1458 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1459 printf(" HYPRE_LSC::parameters mlPostsmootherType = %s\n",
1460 param2);
1461 }
1462 else if ( !strcmp(param1, "mlRelaxType") )
1463 {
1464 sscanf(params[i],"%s %s", param, param2);
1465 rtype = 1;
1466 if ( !strcmp(param2, "jacobi" ) ) rtype = 0;
1467 else if ( !strcmp(param2, "sgs") ) rtype = 1;
1468 else if ( !strcmp(param2, "sgsseq") ) rtype = 2;
1469 else if ( !strcmp(param2, "vbjacobi")) rtype = 3;
1470 else if ( !strcmp(param2, "vbsgs") ) rtype = 4;
1471 else if ( !strcmp(param2, "vbsgsseq")) rtype = 5;
1472 mlPresmootherType_ = rtype;
1473 mlPostsmootherType_ = rtype;
1474 if ( rtype == 6 ) mlPostsmootherType_ = 1;
1475 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1476 printf(" HYPRE_LSC::parameters mlRelaxType = %s\n",
1477 param2);
1478 }
1479
1480 //---------------------------------------------------------------
1481 // mlpack preconditoner : damping factor for Jacobi smoother
1482 //---------------------------------------------------------------
1483
1484 else if ( !strcmp(param1, "mlRelaxWeight") )
1485 {
1486 sscanf(params[i],"%s %lg", param, &weight);
1487 if ( weight < 0.0 || weight > 1.0 ) weight = 0.5;
1488 mlRelaxWeight_ = weight;
1489 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1490 printf(" HYPRE_LSC::parameters mlRelaxWeight = %e\n",
1491 weight);
1492 }
1493
1494 //---------------------------------------------------------------
1495 // mlpack preconditoner : threshold to determine strong coupling
1496 //---------------------------------------------------------------
1497
1498 else if ( !strcmp(param1, "mlStrongThreshold") )
1499 {
1500 sscanf(params[i],"%s %lg", param, &mlStrongThreshold_);
1501 if ( mlStrongThreshold_ < 0.0 || mlStrongThreshold_ > 1.0 )
1502 mlStrongThreshold_ = 0.08;
1503 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1504 printf(" HYPRE_LSC::parameters mlStrongThreshold = %e\n",
1505 mlStrongThreshold_);
1506 }
1507
1508 //---------------------------------------------------------------
1509 // mlpack preconditoner : method to use
1510 //---------------------------------------------------------------
1511
1512 else if ( !strcmp(param1, "mlMethod") )
1513 {
1514 sscanf(params[i],"%s %s", param, param2);
1515 if ( !strcmp(param2, "amg" ) ) mlMethod_ = 0;
1516 else mlMethod_ = 1;
1517 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1518 printf(" HYPRE_LSC::parameters mlMethod = %d\n",mlMethod_);
1519 }
1520
1521 //---------------------------------------------------------------
1522 // mlpack preconditoner : coarse solver to use
1523 //---------------------------------------------------------------
1524
1525 else if ( !strcmp(param1, "mlCoarseSolver") )
1526 {
1527 sscanf(params[i],"%s %s", param, param2);
1528 if ( !strcmp(param2, "superlu" ) ) mlCoarseSolver_ = 0;
1529 else if ( !strcmp(param2, "aggregation" ) ) mlCoarseSolver_ = 1;
1530 else if ( !strcmp(param2, "GS" ) ) mlCoarseSolver_ = 2;
1531 else mlCoarseSolver_ = 1;
1532 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1533 printf(" HYPRE_LSC::parameters mlCoarseSolver = %d\n",
1534 mlCoarseSolver_);
1535 }
1536
1537 //---------------------------------------------------------------
1538 // mlpack preconditoner : coarsening scheme to use
1539 //---------------------------------------------------------------
1540
1541 else if ( !strcmp(param1, "mlCoarsenScheme") )
1542 {
1543 sscanf(params[i],"%s %s", param, param2);
1544 if ( !strcmp(param2, "uncoupled" ) ) mlCoarsenScheme_ = 1;
1545 else if ( !strcmp(param2, "coupled" ) ) mlCoarsenScheme_ = 2;
1546 else if ( !strcmp(param2, "mis" ) ) mlCoarsenScheme_ = 3;
1547 else if ( !strcmp(param2, "hybridum" ) ) mlCoarsenScheme_ = 5;
1548 else if ( !strcmp(param2, "hybriduc" ) ) mlCoarsenScheme_ = 6;
1549 else mlCoarsenScheme_ = 1;
1550 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1551 printf(" HYPRE_LSC::parameters mlCoarsenScheme = %d\n",
1552 mlCoarsenScheme_);
1553 }
1554
1555 //---------------------------------------------------------------
1556 // mlpack preconditoner : no of PDEs (block size)
1557 //---------------------------------------------------------------
1558
1559 else if ( !strcmp(param1, "mlNumPDEs") )
1560 {
1561 sscanf(params[i],"%s %d", param, &mlNumPDEs_);
1562 if ( mlNumPDEs_ < 1 ) mlNumPDEs_ = 1;
1563 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1564 printf(" HYPRE_LSC::parameters mlNumPDEs = %d\n",
1565 mlNumPDEs_);
1566 }
1567
1568 else recognized = 0;
1569 }
1570
1571 //---------------------------------------------------------------
1572 // ams preconditoner : no of PDEs (block size)
1573 //---------------------------------------------------------------
1574
1575 if (!recognized)
1576 {
1577 recognized = 1;
1578 if ( !strcmp(param1, "amsNumPDEs") )
1579 {
1580 sscanf(params[i],"%s %d", param, &amsNumPDEs_);
1581 if ( amsNumPDEs_ < 1 ) amsNumPDEs_ = 1;
1582 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1583 printf(" HYPRE_LSC::parameters amsNumPDEs = %d\n",
1584 amsNumPDEs_);
1585 }
1586
1587 //---------------------------------------------------------------
1588 // ams preconditoner : which smoother to use
1589 //----------------------------------------------------------------
1590
1591 else if ( !strcmp(param1, "amsRelaxType") )
1592 {
1593 sscanf(params[i],"%s %s", param, param2);
1594 if (!strcmp(param2, "jacobi")) amsRelaxType_ = 0;
1595 else if (!strcmp(param2, "scjacobi")) amsRelaxType_ = 1;
1596 else if (!strcmp(param2, "scgs")) amsRelaxType_ = 2;
1597 else if (!strcmp(param2, "kaczmarz")) amsRelaxType_ = 3;
1598 else if (!strcmp(param2, "l1gs*")) amsRelaxType_ = 4;
1599 else if (!strcmp(param2, "hybridsym")) amsRelaxType_ = 6;
1600 else if (!strcmp(param2, "cheby")) amsRelaxType_ = 16;
1601 else amsRelaxType_ = 2;
1602 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1603 printf(" HYPRE_LSC::parameters amsRelaxType = %s\n",
1604 param2);
1605 }
1606
1607 //----------------------------------------------------------------
1608 // ams preconditoner : no of relaxation sweeps per level
1609 //----------------------------------------------------------------
1610
1611 else if (!strcmp(param1, "amsRelaxTimes"))
1612 {
1613 sscanf(params[i],"%s %d", param, &amsRelaxTimes_);
1614 if (amsRelaxTimes_ < 1) amsRelaxTimes_ = 1;
1615 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1616 printf(" HYPRE_LSC::parameters amgRelaxTimes = %d\n",
1617 amsRelaxTimes_);
1618 }
1619
1620 //---------------------------------------------------------------
1621 // ams preconditoner : damping factor for Jacobi smoother
1622 //---------------------------------------------------------------
1623
1624 else if (!strcmp(param1, "amsRelaxWeight"))
1625 {
1626 sscanf(params[i],"%s %lg", param, &amsRelaxWt_);
1627 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1628 printf(" HYPRE_LSC::parameters amgRelaxWeight = %e\n",
1629 amsRelaxWt_);
1630 }
1631
1632 //---------------------------------------------------------------
1633 // ams preconditoner : omega
1634 //---------------------------------------------------------------
1635
1636 else if (!strcmp(param1, "amsRelaxOmega"))
1637 {
1638 sscanf(params[i],"%s %lg", param, &amsRelaxOmega_);
1639 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1640 printf(" HYPRE_LSC::parameters amgRelaxWeight = %e\n",
1641 amsRelaxOmega_);
1642 }
1643
1644 //---------------------------------------------------------------
1645 // ams preconditoner : cycle type
1646 //---------------------------------------------------------------
1647
1648 else if (!strcmp(param1, "amsCycleType"))
1649 {
1650 sscanf(params[i],"%s %d", param, &amsCycleType_);
1651 if (amsCycleType_ < 1) amsCycleType_ = 1;
1652 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1653 printf(" HYPRE_LSC::parameters amgCycleType = %s\n",
1654 param2);
1655 }
1656
1657 //---------------------------------------------------------------
1658 // ams preconditoner : max iterations
1659 //---------------------------------------------------------------
1660
1661 else if (!strcmp(param1, "amsMaxIterations"))
1662 {
1663 sscanf(params[i],"%s %d", param, &amsMaxIter_);
1664 if (amsMaxIter_ < 1) amsMaxIter_ = 1;
1665 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1666 printf(" HYPRE_LSC::parameters amsMaxIterations = %d\n",
1667 amsMaxIter_);
1668 }
1669
1670 //---------------------------------------------------------------
1671 // ams preconditoner : tolerance
1672 //---------------------------------------------------------------
1673
1674 else if (!strcmp(param1, "amsTolerance"))
1675 {
1676 sscanf(params[i],"%s %lg", param, &amsTol_);
1677 if (amsTol_ < 0.0) amsTol_ = 0.0;
1678 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1679 printf(" HYPRE_LSC::parameters amsTolerance = %e\n",
1680 amsTol_);
1681 }
1682
1683 //---------------------------------------------------------------
1684 // ams preconditoner : print level
1685 //---------------------------------------------------------------
1686
1687 else if (!strcmp(param1, "amsPrintLevel"))
1688 {
1689 sscanf(params[i],"%s %d", param, &amsPrintLevel_);
1690 if (amsPrintLevel_ < 0) amsPrintLevel_ = 0;
1691 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1692 printf(" HYPRE_LSC::parameters amsPrintLevel = %d\n",
1693 amsPrintLevel_);
1694 }
1695
1696 //----------------------------------------------------------------
1697 // amg preconditoner : alpha coarsening type
1698 //----------------------------------------------------------------
1699
1700 else if ( !strcmp(param1, "amsAlphaCoarsenType") )
1701 {
1702 sscanf(params[i],"%s %s", param, param2);
1703 if (!strcmp(param2, "cljp")) amsAlphaCoarsenType_ = 0;
1704 else if (!strcmp(param2, "ruge")) amsAlphaCoarsenType_ = 1;
1705 else if (!strcmp(param2, "ruge3c")) amsAlphaCoarsenType_ = 4;
1706 else if (!strcmp(param2, "falgout")) amsAlphaCoarsenType_ = 6;
1707 else if (!strcmp(param2, "pmis")) amsAlphaCoarsenType_ = 8;
1708 else if (!strcmp(param2, "hmis")) amsAlphaCoarsenType_ = 10;
1709 else amsAlphaCoarsenType_ = 0;
1710 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1711 printf(" HYPRE_LSC::parameters amsAlphaCoarsenType = %s\n",
1712 param2);
1713 }
1714
1715 //----------------------------------------------------------------
1716 // amg preconditoner : coarsening type
1717 //----------------------------------------------------------------
1718
1719 else if ( !strcmp(param1, "amsBetaCoarsenType") )
1720 {
1721 sscanf(params[i],"%s %s", param, param2);
1722 if (!strcmp(param2, "cljp")) amsBetaCoarsenType_ = 0;
1723 else if (!strcmp(param2, "ruge")) amsBetaCoarsenType_ = 1;
1724 else if (!strcmp(param2, "ruge3c")) amsBetaCoarsenType_ = 4;
1725 else if (!strcmp(param2, "falgout")) amsBetaCoarsenType_ = 6;
1726 else if (!strcmp(param2, "pmis")) amsBetaCoarsenType_ = 8;
1727 else if (!strcmp(param2, "hmis")) amsBetaCoarsenType_ = 10;
1728 else amsBetaCoarsenType_ = 0;
1729 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1730 printf(" HYPRE_LSC::parameters amsBetaCoarsenType = %s\n",
1731 param2);
1732 }
1733
1734 //---------------------------------------------------------------
1735 // ams preconditoner : levels of aggresive coarseinig
1736 //---------------------------------------------------------------
1737
1738 else if (!strcmp(param1, "amsAlphaAggLevels"))
1739 {
1740 sscanf(params[i],"%s %d", param, &amsAlphaAggLevels_);
1741 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1742 printf(" HYPRE_LSC::parameters amsAlphaAggLevels = %d\n",
1743 amsAlphaAggLevels_);
1744 }
1745
1746 //----------------------------------------------------------------
1747 // ams preconditoner : interpolation type
1748 //----------------------------------------------------------------
1749
1750 else if ( !strcmp(param1, "amsAlphaInterpType") )
1751 {
1752 sscanf(params[i],"%s %d", param, &amsAlphaInterpType_);
1753 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1754 printf(" HYPRE_LSC::parameters amsAlphaInterpType = %d\n",
1755 amsAlphaInterpType_);
1756 }
1757
1758 //----------------------------------------------------------------
1759 // ams preconditoner : interpolation truncation
1760 //----------------------------------------------------------------
1761
1762 else if ( !strcmp(param1, "amsAlphaPmax") )
1763 {
1764 sscanf(params[i],"%s %d", param, &amsAlphaPmax_);
1765 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1766 printf(" HYPRE_LSC::parameters amsAlphaPmax = %d\n",
1767 amsAlphaPmax_);
1768 }
1769
1770 //---------------------------------------------------------------
1771 // ams preconditoner : levels of aggresive coarseinig
1772 //---------------------------------------------------------------
1773
1774 else if (!strcmp(param1, "amsBetaAggLevels"))
1775 {
1776 sscanf(params[i],"%s %d", param, &amsBetaAggLevels_);
1777 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1778 printf(" HYPRE_LSC::parameters amsBetaAggLevels = %d\n",
1779 amsAlphaAggLevels_);
1780 }
1781
1782 //----------------------------------------------------------------
1783 // ams preconditoner : interpolation type
1784 //----------------------------------------------------------------
1785
1786 else if ( !strcmp(param1, "amsBetaInterpType") )
1787 {
1788 sscanf(params[i],"%s %d", param, &amsBetaInterpType_);
1789 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1790 printf(" HYPRE_LSC::parameters amsBetaInterpType = %d\n",
1791 amsBetaInterpType_);
1792 }
1793
1794 //----------------------------------------------------------------
1795 // ams preconditoner : interpolation truncation
1796 //----------------------------------------------------------------
1797
1798 else if ( !strcmp(param1, "amsBetaPmax") )
1799 {
1800 sscanf(params[i],"%s %d", param, &amsBetaPmax_);
1801 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0)
1802 printf(" HYPRE_LSC::parameters amsBetaPmax = %d\n",
1803 amsBetaPmax_);
1804 }
1805
1806 //---------------------------------------------------------------
1807 // ams preconditoner : which smoother to use
1808 //----------------------------------------------------------------
1809
1810 else if ( !strcmp(param1, "amsAlphaRelaxType") )
1811 {
1812 sscanf(params[i],"%s %s", param, param2);
1813 if (!strcmp(param2, "jacobi")) amsAlphaRelaxType_ = 0;
1814 else if (!strcmp(param2, "gsSlow")) amsAlphaRelaxType_ = 1;
1815 else if (!strcmp(param2, "gsFast")) amsAlphaRelaxType_ = 4;
1816 else if (!strcmp(param2, "hybrid")) amsAlphaRelaxType_ = 3;
1817 else if (!strcmp(param2, "hybridsym")) amsAlphaRelaxType_ = 6;
1818 else if (!strcmp(param2, "l1gs" ) ) amsAlphaRelaxType_ = 8;
1819 else amsAlphaRelaxType_ = 4;
1820 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1821 printf(" HYPRE_LSC::parameters amsAlphaRelaxType = %s\n",
1822 param2);
1823 }
1824
1825 //---------------------------------------------------------------
1826 // ams preconditoner : which smoother to use
1827 //----------------------------------------------------------------
1828
1829 else if ( !strcmp(param1, "amsBetaRelaxType") )
1830 {
1831 sscanf(params[i],"%s %s", param, param2);
1832 if (!strcmp(param2, "jacobi")) amsBetaRelaxType_ = 0;
1833 else if (!strcmp(param2, "gsSlow")) amsBetaRelaxType_ = 1;
1834 else if (!strcmp(param2, "gsFast")) amsBetaRelaxType_ = 4;
1835 else if (!strcmp(param2, "hybrid")) amsBetaRelaxType_ = 3;
1836 else if (!strcmp(param2, "hybridsym")) amsBetaRelaxType_ = 6;
1837 else if (!strcmp(param2, "l1gs" ) ) amsBetaRelaxType_ = 8;
1838 else amsBetaRelaxType_ = 4;
1839 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1840 printf(" HYPRE_LSC::parameters amsBetaRelaxType = %s\n",
1841 param2);
1842 }
1843
1844 //---------------------------------------------------------------
1845 // amg preconditoner : threshold to determine strong coupling
1846 //---------------------------------------------------------------
1847
1848 else if ( !strcmp(param1, "amsAlphaStrengthThreshold") )
1849 {
1850 sscanf(params[i],"%s %lg", param, &amsAlphaStrengthThresh_);
1851 if (amsAlphaStrengthThresh_<0.0 || amsAlphaStrengthThresh_>1.0)
1852 amsAlphaStrengthThresh_ = 0.25;
1853 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1854 printf(" HYPRE_LSC::parameters amsAlphaStrengthThresh = %e\n",
1855 amsAlphaStrengthThresh_);
1856 }
1857
1858 //---------------------------------------------------------------
1859 // amg preconditoner : threshold to determine strong coupling
1860 //---------------------------------------------------------------
1861
1862 else if ( !strcmp(param1, "amsBetaStrengthThreshold") )
1863 {
1864 sscanf(params[i],"%s %lg", param, &amsBetaStrengthThresh_);
1865 if (amsBetaStrengthThresh_<0.0 || amsBetaStrengthThresh_>1.0)
1866 amsBetaStrengthThresh_ = 0.25;
1867 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1868 printf(" HYPRE_LSC::parameters amsBetaStrengthThresh = %e\n",
1869 amsBetaStrengthThresh_);
1870 }
1871
1872 //---------------------------------------------------------------
1873 // syspde preconditoner : nvars
1874 //---------------------------------------------------------------
1875
1876 else if ( !strcmp(param1, "syspdeNVars") )
1877 {
1878 sscanf(params[i],"%s %d", param, &sysPDENVars_);
1879 if (sysPDENVars_<0 || sysPDENVars_>10)
1880 sysPDENVars_ = 1;
1881 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 && mypid_ == 0 )
1882 printf(" HYPRE_LSC::parameters syspdeNVars = %d\n",
1883 sysPDENVars_);
1884 }
1885
1886 else recognized = 0;
1887 }
1888
1889 //---------------------------------------------------------------
1890 // error
1891 //---------------------------------------------------------------
1892
1893 if (!recognized)
1894 {
1895 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 && mypid_ == 0 )
1896 printf("HYPRE_LSC::parameters WARNING : %s not recognized\n",
1897 params[i]);
1898 }
1899 }
1900
1901 //-------------------------------------------------------------------
1902 // if reuse is requested, set preconditioner reuse flag
1903 //-------------------------------------------------------------------
1904
1905 if ( reuse == 1 ) HYPreconReuse_ = 1;
1906 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
1907 printf("%4d : HYPRE_LSC::leaving parameters function.\n",mypid_);
1908 return(0);
1909 }
1910
1911 //***************************************************************************
1912 // set up preconditioners for PCG
1913 //---------------------------------------------------------------------------
1914
setupPCGPrecon()1915 void HYPRE_LinSysCore::setupPCGPrecon()
1916 {
1917 //-------------------------------------------------------------------
1918 // if matrix has been reloaded, reset preconditioner
1919 //-------------------------------------------------------------------
1920
1921 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
1922 selectPreconditioner( HYPreconName_ );
1923
1924 //-------------------------------------------------------------------
1925 // set up preconditioners
1926 //-------------------------------------------------------------------
1927
1928 switch ( HYPreconID_ )
1929 {
1930 case HYIDENTITY :
1931 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
1932 printf("No preconditioning \n");
1933 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
1934 HYPRE_DummyFunction, HYPrecon_);
1935 break;
1936
1937 case HYDIAGONAL :
1938 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
1939 printf("Diagonal preconditioning \n");
1940 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
1941 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
1942 HYPRE_DummyFunction, HYPrecon_);
1943 else
1944 {
1945 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
1946 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
1947 HYPreconSetup_ = 1;
1948 }
1949 break;
1950
1951 case HYPILUT :
1952 if ( mypid_ == 0 )
1953 printf("HYPRE_LSI : PCG does not work with pilut.\n");
1954 exit(1);
1955 break;
1956
1957 case HYDDILUT :
1958 if ( mypid_ == 0 )
1959 printf("HYPRE_LSI : PCG does not work with ddilut.\n");
1960 exit(1);
1961 break;
1962
1963 case HYDDICT :
1964 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
1965 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
1966 HYPRE_DummyFunction, HYPrecon_);
1967 else
1968 {
1969 setupPreconDDICT();
1970 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
1971 HYPRE_LSI_DDICTSetup, HYPrecon_);
1972 HYPreconSetup_ = 1;
1973 }
1974 break;
1975
1976 case HYSCHWARZ :
1977 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
1978 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
1979 HYPRE_DummyFunction, HYPrecon_);
1980 else
1981 {
1982 setupPreconSchwarz();
1983 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
1984 HYPRE_LSI_SchwarzSetup, HYPrecon_);
1985 HYPreconSetup_ = 1;
1986 }
1987 break;
1988
1989 case HYPOLY :
1990 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
1991 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
1992 HYPRE_DummyFunction, HYPrecon_);
1993 else
1994 {
1995 setupPreconPoly();
1996 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
1997 HYPRE_LSI_PolySetup, HYPrecon_);
1998 HYPreconSetup_ = 1;
1999 }
2000 break;
2001
2002 case HYPARASAILS :
2003 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2004 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRParaSailsSolve,
2005 HYPRE_DummyFunction, HYPrecon_);
2006 else
2007 {
2008 setupPreconParaSails();
2009 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRParaSailsSolve,
2010 HYPRE_ParCSRParaSailsSetup, HYPrecon_);
2011 HYPreconSetup_ = 1;
2012 }
2013 break;
2014
2015 case HYBOOMERAMG :
2016 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2017 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2018 HYPRE_DummyFunction, HYPrecon_);
2019 else
2020 {
2021 setupPreconBoomerAMG();
2022 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2023 HYPRE_BoomerAMGSetup, HYPrecon_);
2024 HYPreconSetup_ = 1;
2025 }
2026 break;
2027
2028 case HYEUCLID :
2029 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2030 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2031 HYPRE_DummyFunction, HYPrecon_);
2032 else
2033 {
2034 setupPreconEuclid();
2035 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2036 HYPRE_EuclidSetup, HYPrecon_);
2037 HYPreconSetup_ = 1;
2038 }
2039 break;
2040
2041 case HYBLOCK :
2042 printf("PCG : block preconditioning not available.\n");
2043 exit(1);
2044 break;
2045
2046 case HYML :
2047 #ifdef HAVE_ML
2048 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2049 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
2050 HYPRE_DummyFunction, HYPrecon_);
2051 else
2052 {
2053 setupPreconML();
2054 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
2055 HYPRE_LSI_MLSetup, HYPrecon_);
2056 HYPreconSetup_ = 1;
2057 }
2058 #else
2059 printf("PCG : ML preconditioning not available.\n");
2060 #endif
2061 break;
2062
2063 case HYMLMAXWELL :
2064 #ifdef HAVE_MLMAXWELL
2065 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2066 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_MLMaxwellSolve,
2067 HYPRE_DummyFunction, HYPrecon_);
2068 else
2069 {
2070 setupPreconMLMaxwell();
2071 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_MLMaxwellSolve,
2072 HYPRE_LSI_MLMaxwellSetup, HYPrecon_);
2073 HYPreconSetup_ = 1;
2074 }
2075 #else
2076 printf("PCG : ML preconditioning not available.\n");
2077 #endif
2078 break;
2079
2080 case HYMLI :
2081 #ifdef HAVE_MLI
2082 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2083 printf("MLI preconditioning\n");
2084 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2085 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
2086 HYPRE_DummyFunction, HYPrecon_);
2087 else
2088 {
2089 HYPRE_ParCSRPCGSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
2090 HYPRE_LSI_MLISetup, HYPrecon_);
2091 HYPreconSetup_ = 1;
2092 }
2093 #else
2094 printf("PCG : MLI preconditioning not available.\n");
2095 #endif
2096 break;
2097
2098 case HYAMS :
2099 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2100 printf("AMS preconditioning\n");
2101 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2102 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_AMSSolve,
2103 HYPRE_DummyFunction, HYPrecon_);
2104 else
2105 {
2106 setupPreconAMS();
2107 HYPRE_ParCSRPCGSetPrecond(HYSolver_,HYPRE_AMSSolve,
2108 HYPRE_AMSSetup, HYPrecon_);
2109 HYPreconSetup_ = 1;
2110 }
2111 break;
2112
2113 case HYUZAWA :
2114 printf("PCG : Uzawa preconditioning not available.\n");
2115 exit(1);
2116 break;
2117
2118 case HYSYSPDE :
2119 #ifdef HAVE_SYSPDE
2120 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2121 printf("SysPDe preconditioning\n");
2122 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2123 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2124 HYPRE_DummyFunction, HYPrecon_);
2125 else
2126 {
2127 setupPreconSysPDE();
2128 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2129 HYPRE_ParCSRSysPDESetup, HYPrecon_);
2130 HYPreconSetup_ = 1;
2131 }
2132 #else
2133 printf("PCG : SysPDE preconditioning not available.\n");
2134 #endif
2135 break;
2136
2137 case HYDSLU :
2138 #ifdef HYPRE_USING_DSUPERLU
2139 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2140 printf("DSuperLU preconditioning\n");
2141 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2142 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2143 HYPRE_DummyFunction, HYPrecon_);
2144 else
2145 {
2146 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
2147 HYPRE_ParCSRPCGSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2148 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
2149 HYPreconSetup_ = 1;
2150 }
2151 #else
2152 printf("PCG : DSuperLU preconditioning not available.\n");
2153 #endif
2154 }
2155 return;
2156 }
2157
2158 //***************************************************************************
2159 // set up preconditioners for LSICG
2160 //---------------------------------------------------------------------------
2161
setupLSICGPrecon()2162 void HYPRE_LinSysCore::setupLSICGPrecon()
2163 {
2164 //-------------------------------------------------------------------
2165 // if matrix has been reloaded, reset preconditioner
2166 //-------------------------------------------------------------------
2167
2168 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
2169 selectPreconditioner( HYPreconName_ );
2170
2171 //-------------------------------------------------------------------
2172 // set up preconditioners
2173 //-------------------------------------------------------------------
2174
2175 switch ( HYPreconID_ )
2176 {
2177 case HYIDENTITY :
2178 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2179 printf("No preconditioning \n");
2180 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
2181 HYPRE_DummyFunction, HYPrecon_);
2182 break;
2183
2184 case HYDIAGONAL :
2185 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2186 printf("Diagonal preconditioning \n");
2187 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2188 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2189 HYPRE_DummyFunction, HYPrecon_);
2190 else
2191 {
2192 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2193 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
2194 HYPreconSetup_ = 1;
2195 }
2196 break;
2197
2198 case HYPILUT :
2199 if ( mypid_ == 0 )
2200 printf("HYPRE_LSI : LSICG does not work with pilut.\n");
2201 exit(1);
2202 break;
2203
2204 case HYDDILUT :
2205 if ( mypid_ == 0 )
2206 printf("HYPRE_LSI : LSICG does not work with ddilut.\n");
2207 exit(1);
2208 break;
2209
2210 case HYDDICT :
2211 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2212 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2213 HYPRE_DummyFunction, HYPrecon_);
2214 else
2215 {
2216 setupPreconDDICT();
2217 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2218 HYPRE_LSI_DDICTSetup, HYPrecon_);
2219 HYPreconSetup_ = 1;
2220 }
2221 break;
2222
2223 case HYSCHWARZ :
2224 if ( mypid_ == 0 )
2225 printf("HYPRE_LSI : LSICG does not work with Schwarz.\n");
2226 exit(1);
2227 break;
2228
2229 case HYPOLY :
2230 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2231 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2232 HYPRE_DummyFunction, HYPrecon_);
2233 else
2234 {
2235 setupPreconPoly();
2236 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2237 HYPRE_LSI_PolySetup, HYPrecon_);
2238 HYPreconSetup_ = 1;
2239 }
2240 break;
2241
2242 case HYPARASAILS :
2243 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2244 HYPRE_ParCSRLSICGSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
2245 HYPRE_DummyFunction, HYPrecon_);
2246 else
2247 {
2248 setupPreconParaSails();
2249 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_ParCSRParaSailsSolve,
2250 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
2251 HYPreconSetup_ = 1;
2252 }
2253 break;
2254
2255 case HYBOOMERAMG :
2256 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2257 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2258 HYPRE_DummyFunction, HYPrecon_);
2259 else
2260 {
2261 setupPreconBoomerAMG();
2262 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2263 HYPRE_BoomerAMGSetup, HYPrecon_);
2264 HYPreconSetup_ = 1;
2265 }
2266 break;
2267
2268 case HYEUCLID :
2269 if ( mypid_ == 0 )
2270 printf("HYPRE_LSI : LSICG does not work with Euclid.\n");
2271 exit(1);
2272 break;
2273
2274 case HYBLOCK :
2275 if ( mypid_ == 0 )
2276 printf("HYPRE_LSI : LSICG does not work with blkprec.\n");
2277 exit(1);
2278 break;
2279
2280 case HYML :
2281 printf("HYPRE_LSI : LSICG - MLI preconditioning not available.\n");
2282 break;
2283
2284 case HYMLMAXWELL :
2285 printf("HYPRE_LSI : LSICG - MLMAXWELL not available.\n");
2286 break;
2287
2288 case HYMLI :
2289 #ifdef HAVE_MLI
2290 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2291 printf("MLI preconditioning\n");
2292 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2293 HYPRE_ParCSRLSICGSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
2294 HYPRE_DummyFunction, HYPrecon_);
2295 else
2296 {
2297 HYPRE_ParCSRLSICGSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
2298 HYPRE_LSI_MLISetup, HYPrecon_);
2299 HYPreconSetup_ = 1;
2300 }
2301 #else
2302 printf("HYPRE_LSI LSICG : MLI preconditioning not available.\n");
2303 #endif
2304 break;
2305
2306 case HYUZAWA :
2307 if ( mypid_ == 0 )
2308 printf("HYPRE_LSI : LSICG does not work with Uzawa.\n");
2309 exit(1);
2310 break;
2311
2312 default :
2313 printf("CG : preconditioner unknown.\n");
2314 exit(1);
2315 break;
2316 }
2317 return;
2318 }
2319
2320 //***************************************************************************
2321 // set up preconditioners for GMRES
2322 //---------------------------------------------------------------------------
2323
setupGMRESPrecon()2324 void HYPRE_LinSysCore::setupGMRESPrecon()
2325 {
2326 //-------------------------------------------------------------------
2327 // if matrix has been reloaded, reset preconditioner
2328 //-------------------------------------------------------------------
2329
2330 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
2331 selectPreconditioner( HYPreconName_ );
2332
2333 //-------------------------------------------------------------------
2334 // set up preconditioners
2335 //-------------------------------------------------------------------
2336
2337 switch ( HYPreconID_ )
2338 {
2339 case HYIDENTITY :
2340 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2341 printf("No preconditioning \n");
2342 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
2343 HYPRE_DummyFunction, HYPrecon_);
2344 break;
2345
2346 case HYDIAGONAL :
2347 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2348 printf("Diagonal preconditioning \n");
2349 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2350 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2351 HYPRE_DummyFunction, HYPrecon_);
2352 else
2353 {
2354 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2355 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
2356 HYPreconSetup_ = 1;
2357 }
2358 break;
2359
2360 case HYPILUT :
2361 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2362 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2363 HYPRE_DummyFunction, HYPrecon_);
2364 else
2365 {
2366 setupPreconPILUT();
2367 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2368 HYPRE_ParCSRPilutSetup, HYPrecon_);
2369 HYPreconSetup_ = 1;
2370 }
2371 break;
2372
2373 case HYDDILUT :
2374 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2375 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2376 HYPRE_DummyFunction, HYPrecon_);
2377 else
2378 {
2379 setupPreconDDILUT();
2380 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2381 HYPRE_LSI_DDIlutSetup, HYPrecon_);
2382 HYPreconSetup_ = 1;
2383 }
2384 break;
2385
2386 case HYDDICT :
2387 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2388 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2389 HYPRE_DummyFunction, HYPrecon_);
2390 else
2391 {
2392 setupPreconDDICT();
2393 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2394 HYPRE_LSI_DDICTSetup, HYPrecon_);
2395 HYPreconSetup_ = 1;
2396 }
2397 break;
2398
2399 case HYSCHWARZ :
2400 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2401 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2402 HYPRE_DummyFunction, HYPrecon_);
2403 else
2404 {
2405 setupPreconSchwarz();
2406 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2407 HYPRE_LSI_SchwarzSetup, HYPrecon_);
2408 HYPreconSetup_ = 1;
2409 }
2410 break;
2411
2412 case HYPOLY :
2413 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2414 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2415 HYPRE_DummyFunction, HYPrecon_);
2416 else
2417 {
2418 setupPreconPoly();
2419 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2420 HYPRE_LSI_PolySetup, HYPrecon_);
2421 HYPreconSetup_ = 1;
2422 }
2423 break;
2424
2425 case HYPARASAILS :
2426 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2427 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
2428 HYPRE_DummyFunction, HYPrecon_);
2429 else
2430 {
2431 setupPreconParaSails();
2432 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
2433 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
2434 HYPreconSetup_ = 1;
2435 }
2436 break;
2437
2438 case HYBOOMERAMG :
2439 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2440 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
2441 HYPRE_DummyFunction, HYPrecon_);
2442 else
2443 {
2444 setupPreconBoomerAMG();
2445 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2446 HYPRE_BoomerAMGSetup, HYPrecon_);
2447 HYPreconSetup_ = 1;
2448 }
2449 break;
2450
2451 case HYEUCLID :
2452 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2453 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2454 HYPRE_DummyFunction, HYPrecon_);
2455 else
2456 {
2457 setupPreconEuclid();
2458 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2459 HYPRE_EuclidSetup, HYPrecon_);
2460 HYPreconSetup_ = 1;
2461 }
2462 break;
2463
2464 case HYBLOCK :
2465 printf("GMRES : block preconditioning not available.\n");
2466 exit(1);
2467 break;
2468
2469 case HYML :
2470 #ifdef HAVE_ML
2471 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2472 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
2473 HYPRE_DummyFunction, HYPrecon_);
2474 else
2475 {
2476 setupPreconML();
2477 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
2478 HYPRE_LSI_MLSetup, HYPrecon_);
2479 HYPreconSetup_ = 1;
2480 }
2481 #else
2482 printf("GMRES : ML preconditioning not available.\n");
2483 #endif
2484 break;
2485
2486 case HYMLMAXWELL :
2487 #ifdef HAVE_MLMAXWELL
2488 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2489 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_MLMaxwellSolve,
2490 HYPRE_DummyFunction, HYPrecon_);
2491 else
2492 {
2493 setupPreconMLMaxwell();
2494 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_LSI_MLMaxwellSolve,
2495 HYPRE_LSI_MLMaxwellSetup, HYPrecon_);
2496 HYPreconSetup_ = 1;
2497 }
2498 #else
2499 printf("GMRES : ML preconditioning not available.\n");
2500 #endif
2501 break;
2502
2503 case HYMLI :
2504 #ifdef HAVE_MLI
2505 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2506 printf("MLI preconditioning \n");
2507 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2508 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
2509 HYPRE_DummyFunction, HYPrecon_);
2510 else
2511 {
2512 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
2513 HYPRE_LSI_MLISetup, HYPrecon_);
2514 HYPreconSetup_ = 1;
2515 }
2516 #else
2517 printf("GMRES : MLI preconditioning not available.\n");
2518 #endif
2519 break;
2520
2521 case HYUZAWA :
2522 printf("GMRES : Uzawa preconditioning not available.\n");
2523 exit(1);
2524 break;
2525
2526 case HYAMS :
2527 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2528 printf("AMS preconditioning\n");
2529 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2530 HYPRE_ParCSRGMRESSetPrecond(HYSolver_,HYPRE_AMSSolve,
2531 HYPRE_DummyFunction, HYPrecon_);
2532 else
2533 {
2534 setupPreconAMS();
2535 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_AMSSolve,
2536 HYPRE_AMSSetup, HYPrecon_);
2537 HYPreconSetup_ = 1;
2538 }
2539 break;
2540
2541 case HYSYSPDE :
2542 #ifdef HY_SYSPDE
2543 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2544 printf("SysPDe preconditioning\n");
2545 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2546 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2547 HYPRE_DummyFunction, HYPrecon_);
2548 else
2549 {
2550 setupPreconSysPDE();
2551 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2552 HYPRE_ParCSRSysPDESetup, HYPrecon_);
2553 HYPreconSetup_ = 1;
2554 }
2555 #else
2556 printf("GMRES : SysPDe preconditioning not available.\n");
2557 #endif
2558 break;
2559
2560 case HYDSLU :
2561 #ifdef HYPRE_USING_DSUPERLU
2562 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2563 printf("DSuperLU preconditioning\n");
2564 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2565 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2566 HYPRE_DummyFunction, HYPrecon_);
2567 else
2568 {
2569 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
2570 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2571 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
2572 HYPreconSetup_ = 1;
2573 }
2574 #else
2575 printf("GMRES : DSuperLU preconditioning not available.\n");
2576 #endif
2577 break;
2578 }
2579 return;
2580 }
2581
2582 //***************************************************************************
2583 // set up preconditioners for FGMRES
2584 //---------------------------------------------------------------------------
2585
setupFGMRESPrecon()2586 void HYPRE_LinSysCore::setupFGMRESPrecon()
2587 {
2588 //-------------------------------------------------------------------
2589 // if matrix has been reloaded, reset preconditioner
2590 //-------------------------------------------------------------------
2591
2592 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
2593 selectPreconditioner( HYPreconName_ );
2594
2595 //-------------------------------------------------------------------
2596 // set up preconditioners
2597 //-------------------------------------------------------------------
2598
2599 switch ( HYPreconID_ )
2600 {
2601 case HYIDENTITY :
2602 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2603 printf("No preconditioning \n");
2604 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
2605 HYPRE_DummyFunction, HYPrecon_);
2606 break;
2607
2608 case HYDIAGONAL :
2609 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2610 printf("Diagonal preconditioning \n");
2611 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2612 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2613 HYPRE_DummyFunction, HYPrecon_);
2614 else
2615 {
2616 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2617 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
2618 HYPreconSetup_ = 1;
2619 }
2620 break;
2621
2622 case HYPILUT :
2623 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2624 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2625 HYPRE_DummyFunction, HYPrecon_);
2626 else
2627 {
2628 setupPreconPILUT();
2629 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2630 HYPRE_ParCSRPilutSetup, HYPrecon_);
2631 HYPreconSetup_ = 1;
2632 }
2633 break;
2634
2635 case HYDDILUT :
2636 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2637 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2638 HYPRE_DummyFunction, HYPrecon_);
2639 else
2640 {
2641 setupPreconDDILUT();
2642 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2643 HYPRE_LSI_DDIlutSetup, HYPrecon_);
2644 HYPreconSetup_ = 1;
2645 }
2646 break;
2647
2648 case HYDDICT :
2649 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2650 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2651 HYPRE_DummyFunction, HYPrecon_);
2652 else
2653 {
2654 setupPreconDDICT();
2655 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2656 HYPRE_LSI_DDICTSetup, HYPrecon_);
2657 HYPreconSetup_ = 1;
2658 }
2659 break;
2660
2661 case HYSCHWARZ :
2662 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2663 HYPRE_ParCSRGMRESSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2664 HYPRE_DummyFunction, HYPrecon_);
2665 else
2666 {
2667 setupPreconSchwarz();
2668 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2669 HYPRE_LSI_SchwarzSetup, HYPrecon_);
2670 HYPreconSetup_ = 1;
2671 }
2672 break;
2673
2674 case HYPOLY :
2675 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2676 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2677 HYPRE_DummyFunction, HYPrecon_);
2678 else
2679 {
2680 setupPreconPoly();
2681 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2682 HYPRE_LSI_PolySetup, HYPrecon_);
2683 HYPreconSetup_ = 1;
2684 }
2685 break;
2686
2687 case HYPARASAILS :
2688 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2689 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
2690 HYPRE_DummyFunction, HYPrecon_);
2691 else
2692 {
2693 setupPreconParaSails();
2694 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
2695 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
2696 HYPreconSetup_ = 1;
2697 }
2698 break;
2699
2700 case HYBOOMERAMG :
2701 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2702 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
2703 HYPRE_DummyFunction, HYPrecon_);
2704 else
2705 {
2706 setupPreconBoomerAMG();
2707 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2708 HYPRE_BoomerAMGSetup, HYPrecon_);
2709 HYPreconSetup_ = 1;
2710 }
2711 break;
2712
2713 case HYEUCLID :
2714 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2715 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2716 HYPRE_DummyFunction, HYPrecon_);
2717 else
2718 {
2719 setupPreconEuclid();
2720 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2721 HYPRE_EuclidSetup, HYPrecon_);
2722 HYPreconSetup_ = 1;
2723 }
2724 break;
2725
2726 case HYBLOCK :
2727 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2728 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,
2729 HYPRE_LSI_BlockPrecondSolve, HYPRE_DummyFunction,
2730 HYPrecon_);
2731 else
2732 {
2733 setupPreconBlock();
2734 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,
2735 HYPRE_LSI_BlockPrecondSolve,
2736 HYPRE_LSI_BlockPrecondSetup, HYPrecon_);
2737 HYPreconSetup_ = 1;
2738 }
2739 break;
2740
2741 case HYML :
2742 #ifdef HAVE_ML
2743 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2744 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
2745 HYPRE_DummyFunction, HYPrecon_);
2746 else
2747 {
2748 setupPreconML();
2749 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
2750 HYPRE_LSI_MLSetup, HYPrecon_);
2751 HYPreconSetup_ = 1;
2752 }
2753 #else
2754 printf("FGMRES : ML preconditioning not available.\n");
2755 #endif
2756 break;
2757
2758 case HYMLMAXWELL :
2759 printf("FGMRES : MLMaxwell preconditioning not available.\n");
2760 break;
2761
2762 case HYMLI :
2763 #ifdef HAVE_MLI
2764 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2765 printf("MLI preconditioning \n");
2766 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2767 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
2768 HYPRE_DummyFunction, HYPrecon_);
2769 else
2770 {
2771 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
2772 HYPRE_LSI_MLISetup, HYPrecon_);
2773 HYPreconSetup_ = 1;
2774 }
2775 #else
2776 printf("FGMRES : ML preconditioning not available.\n");
2777 #endif
2778 break;
2779
2780 case HYUZAWA :
2781 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2782 printf("Uzawa preconditioning \n");
2783 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2784 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_UzawaSolve,
2785 HYPRE_DummyFunction, HYPrecon_);
2786 else
2787 {
2788 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_LSI_UzawaSolve,
2789 HYPRE_LSI_UzawaSetup, HYPrecon_);
2790 HYPreconSetup_ = 1;
2791 }
2792 break;
2793
2794 case HYAMS :
2795 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2796 printf("AMS preconditioning\n");
2797 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2798 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_,HYPRE_AMSSolve,
2799 HYPRE_DummyFunction, HYPrecon_);
2800 else
2801 {
2802 setupPreconAMS();
2803 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_AMSSolve,
2804 HYPRE_AMSSetup, HYPrecon_);
2805 HYPreconSetup_ = 1;
2806 }
2807 break;
2808
2809 case HYSYSPDE :
2810 #ifdef HY_SYSPDE
2811 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2812 printf("SysPDe preconditioning\n");
2813 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2814 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2815 HYPRE_DummyFunction, HYPrecon_);
2816 else
2817 {
2818 setupPreconSysPDE();
2819 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
2820 HYPRE_ParCSRSysPDESetup, HYPrecon_);
2821 HYPreconSetup_ = 1;
2822 }
2823 #else
2824 printf("FGMRES : SysPDe preconditioning not available.\n");
2825 #endif
2826 break;
2827
2828 case HYDSLU :
2829 #ifdef HYPRE_USING_DSUPERLU
2830 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
2831 printf("DSuperLU preconditioning\n");
2832 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2833 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2834 HYPRE_DummyFunction, HYPrecon_);
2835 else
2836 {
2837 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
2838 HYPRE_ParCSRFGMRESSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
2839 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
2840 HYPreconSetup_ = 1;
2841 }
2842 #else
2843 printf("FGMRES : DSuperLU preconditioning not available.\n");
2844 #endif
2845 break;
2846 }
2847 return;
2848 }
2849
2850 //***************************************************************************
2851 // set up preconditioners for BiCGSTAB
2852 //---------------------------------------------------------------------------
2853
setupBiCGSTABPrecon()2854 void HYPRE_LinSysCore::setupBiCGSTABPrecon()
2855 {
2856 //-------------------------------------------------------------------
2857 // if matrix has been reloaded, reset preconditioner
2858 //-------------------------------------------------------------------
2859
2860 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
2861 selectPreconditioner( HYPreconName_ );
2862
2863 //-------------------------------------------------------------------
2864 // set up preconditioners
2865 //-------------------------------------------------------------------
2866
2867 switch ( HYPreconID_ )
2868 {
2869 case HYIDENTITY :
2870 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2871 printf("No preconditioning \n");
2872 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
2873 HYPRE_DummyFunction, HYPrecon_);
2874 break;
2875
2876 case HYDIAGONAL :
2877 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
2878 printf("Diagonal preconditioning \n");
2879 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2880 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2881 HYPRE_DummyFunction, HYPrecon_);
2882 else
2883 {
2884 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
2885 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
2886 HYPreconSetup_ = 1;
2887 }
2888 break;
2889
2890 case HYPILUT :
2891 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2892 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2893 HYPRE_DummyFunction, HYPrecon_);
2894 else
2895 {
2896 setupPreconPILUT();
2897 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
2898 HYPRE_ParCSRPilutSetup,HYPrecon_);
2899 HYPreconSetup_ = 1;
2900 }
2901 break;
2902
2903 case HYDDILUT :
2904 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2905 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2906 HYPRE_DummyFunction, HYPrecon_);
2907 else
2908 {
2909 setupPreconDDILUT();
2910 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
2911 HYPRE_LSI_DDIlutSetup, HYPrecon_);
2912 HYPreconSetup_ = 1;
2913 }
2914 break;
2915
2916 case HYDDICT :
2917 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2918 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2919 HYPRE_DummyFunction, HYPrecon_);
2920 else
2921 {
2922 setupPreconDDICT();
2923 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
2924 HYPRE_LSI_DDICTSetup, HYPrecon_);
2925 HYPreconSetup_ = 1;
2926 }
2927 break;
2928
2929 case HYSCHWARZ :
2930 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2931 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2932 HYPRE_DummyFunction, HYPrecon_);
2933 else
2934 {
2935 setupPreconSchwarz();
2936 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
2937 HYPRE_LSI_SchwarzSetup,HYPrecon_);
2938 HYPreconSetup_ = 1;
2939 }
2940 break;
2941
2942 case HYPOLY :
2943 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2944 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2945 HYPRE_DummyFunction, HYPrecon_);
2946 else
2947 {
2948 setupPreconPoly();
2949 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
2950 HYPRE_LSI_PolySetup, HYPrecon_);
2951 HYPreconSetup_ = 1;
2952 }
2953 break;
2954
2955 case HYPARASAILS :
2956 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2957 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,
2958 HYPRE_ParCSRParaSailsSolve,
2959 HYPRE_DummyFunction, HYPrecon_);
2960 else
2961 {
2962 setupPreconParaSails();
2963 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,
2964 HYPRE_ParCSRParaSailsSolve,
2965 HYPRE_ParCSRParaSailsSetup,
2966 HYPrecon_);
2967 HYPreconSetup_ = 1;
2968 }
2969 break;
2970
2971 case HYBOOMERAMG :
2972 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2973 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
2974 HYPRE_DummyFunction, HYPrecon_);
2975 else
2976 {
2977 setupPreconBoomerAMG();
2978 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
2979 HYPRE_BoomerAMGSetup, HYPrecon_);
2980 HYPreconSetup_ = 1;
2981 }
2982 break;
2983
2984 case HYEUCLID :
2985 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
2986 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2987 HYPRE_DummyFunction, HYPrecon_);
2988 else
2989 {
2990 setupPreconEuclid();
2991 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_EuclidSolve,
2992 HYPRE_EuclidSetup,HYPrecon_);
2993 HYPreconSetup_ = 1;
2994 }
2995 break;
2996
2997 case HYBLOCK :
2998 printf("BiCGSTAB : block preconditioning not available.\n");
2999 exit(1);
3000 break;
3001
3002 case HYML :
3003 #ifdef HAVE_ML
3004 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3005 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
3006 HYPRE_DummyFunction, HYPrecon_);
3007 else
3008 {
3009 setupPreconML();
3010 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
3011 HYPRE_LSI_MLSetup, HYPrecon_);
3012 HYPreconSetup_ = 1;
3013 }
3014 #else
3015 printf("BiCGSTAB : ML preconditioning not available.\n");
3016 #endif
3017 break;
3018
3019 case HYMLMAXWELL :
3020 printf("BiCGSTAB : MLMaxwell preconditioning not available.\n");
3021 break;
3022
3023 case HYMLI :
3024 #ifdef HAVE_MLI
3025 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3026 printf("MLI preconditioning\n");
3027 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3028 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
3029 HYPRE_DummyFunction, HYPrecon_);
3030 else
3031 {
3032 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
3033 HYPRE_LSI_MLISetup, HYPrecon_);
3034 HYPreconSetup_ = 1;
3035 }
3036 #else
3037 printf("BiCGSTAB : MLI preconditioning not available.\n");
3038 #endif
3039 break;
3040
3041 case HYUZAWA :
3042 printf("BiCGSTAB : Uzawa preconditioning not available.\n");
3043 exit(1);
3044 break;
3045
3046 case HYAMS :
3047 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3048 printf("AMS preconditioning\n");
3049 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3050 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_,HYPRE_AMSSolve,
3051 HYPRE_DummyFunction, HYPrecon_);
3052 else
3053 {
3054 setupPreconAMS();
3055 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_AMSSolve,
3056 HYPRE_AMSSetup, HYPrecon_);
3057 HYPreconSetup_ = 1;
3058 }
3059 break;
3060
3061 case HYSYSPDE :
3062 #ifdef HY_SYSPDE
3063 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3064 printf("SysPDe preconditioning\n");
3065 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3066 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3067 HYPRE_DummyFunction, HYPrecon_);
3068 else
3069 {
3070 setupPreconSysPDE();
3071 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3072 HYPRE_ParCSRSysPDESetup, HYPrecon_);
3073 HYPreconSetup_ = 1;
3074 }
3075 #else
3076 printf("BiCGSTAB : SysPDe preconditioning not available.\n");
3077 #endif
3078 break;
3079
3080 case HYDSLU :
3081 #ifdef HYPRE_USING_DSUPERLU
3082 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3083 printf("DSuperLU preconditioning\n");
3084 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3085 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3086 HYPRE_DummyFunction, HYPrecon_);
3087 else
3088 {
3089 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
3090 HYPRE_ParCSRBiCGSTABSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3091 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
3092 HYPreconSetup_ = 1;
3093 }
3094 #else
3095 printf("BiCGSTAB : DSuperLU preconditioning not available.\n");
3096 #endif
3097 break;
3098 }
3099 return;
3100 }
3101
3102 //***************************************************************************
3103 // set up preconditioners for BiCGSTABL
3104 //---------------------------------------------------------------------------
3105
setupBiCGSTABLPrecon()3106 void HYPRE_LinSysCore::setupBiCGSTABLPrecon()
3107 {
3108 //-------------------------------------------------------------------
3109 // if matrix has been reloaded, reset preconditioner
3110 //-------------------------------------------------------------------
3111
3112 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
3113 selectPreconditioner( HYPreconName_ );
3114
3115 //-------------------------------------------------------------------
3116 // set up preconditioners
3117 //-------------------------------------------------------------------
3118
3119 switch ( HYPreconID_ )
3120 {
3121 case HYIDENTITY :
3122 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3123 printf("No preconditioning \n");
3124 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
3125 HYPRE_DummyFunction, HYPrecon_);
3126 break;
3127
3128 case HYDIAGONAL :
3129 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3130 printf("Diagonal preconditioning \n");
3131 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3132 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3133 HYPRE_DummyFunction, HYPrecon_);
3134 else
3135 {
3136 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3137 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
3138 HYPreconSetup_ = 1;
3139 }
3140 break;
3141
3142 case HYPILUT :
3143 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3144 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_ParCSRPilutSolve,
3145 HYPRE_DummyFunction, HYPrecon_);
3146 else
3147 {
3148 setupPreconPILUT();
3149 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_ParCSRPilutSolve,
3150 HYPRE_ParCSRPilutSetup, HYPrecon_);
3151 HYPreconSetup_ = 1;
3152 }
3153 break;
3154
3155 case HYDDILUT :
3156 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3157 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3158 HYPRE_DummyFunction, HYPrecon_);
3159 else
3160 {
3161 setupPreconDDILUT();
3162 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3163 HYPRE_LSI_DDIlutSetup,HYPrecon_);
3164 HYPreconSetup_ = 1;
3165 }
3166 break;
3167
3168 case HYDDICT :
3169 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3170 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3171 HYPRE_DummyFunction, HYPrecon_);
3172 else
3173 {
3174 setupPreconDDICT();
3175 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3176 HYPRE_LSI_DDICTSetup, HYPrecon_);
3177 HYPreconSetup_ = 1;
3178 }
3179 break;
3180
3181 case HYSCHWARZ :
3182 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3183 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_LSI_SchwarzSolve,
3184 HYPRE_DummyFunction, HYPrecon_);
3185 else
3186 {
3187 setupPreconSchwarz();
3188 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_LSI_SchwarzSolve,
3189 HYPRE_LSI_SchwarzSetup,HYPrecon_);
3190 HYPreconSetup_ = 1;
3191 }
3192 break;
3193
3194 case HYPOLY :
3195 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3196 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3197 HYPRE_DummyFunction, HYPrecon_);
3198 else
3199 {
3200 setupPreconPoly();
3201 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3202 HYPRE_LSI_PolySetup, HYPrecon_);
3203 HYPreconSetup_ = 1;
3204 }
3205 break;
3206
3207 case HYPARASAILS :
3208 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3209 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,
3210 HYPRE_ParCSRParaSailsSolve,
3211 HYPRE_DummyFunction, HYPrecon_);
3212 else
3213 {
3214 setupPreconParaSails();
3215 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,
3216 HYPRE_ParCSRParaSailsSolve,
3217 HYPRE_ParCSRParaSailsSetup,
3218 HYPrecon_);
3219 HYPreconSetup_ = 1;
3220 }
3221 break;
3222
3223 case HYBOOMERAMG :
3224 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3225 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
3226 HYPRE_DummyFunction, HYPrecon_);
3227 else
3228 {
3229 setupPreconBoomerAMG();
3230 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
3231 HYPRE_BoomerAMGSetup, HYPrecon_);
3232 HYPreconSetup_ = 1;
3233 }
3234 break;
3235
3236 case HYEUCLID :
3237 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3238 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_EuclidSolve,
3239 HYPRE_DummyFunction, HYPrecon_);
3240 else
3241 {
3242 setupPreconEuclid();
3243 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_EuclidSolve,
3244 HYPRE_EuclidSetup,
3245 HYPrecon_);
3246 HYPreconSetup_ = 1;
3247 }
3248 break;
3249
3250 case HYBLOCK :
3251 printf("BiCGSTABL : block preconditioning not available.\n");
3252 exit(1);
3253 break;
3254
3255 case HYML :
3256 #ifdef HAVE_ML
3257 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3258 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
3259 HYPRE_DummyFunction, HYPrecon_);
3260 else
3261 {
3262 setupPreconML();
3263 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
3264 HYPRE_LSI_MLSetup, HYPrecon_);
3265 HYPreconSetup_ = 1;
3266 }
3267 #else
3268 printf("BiCGSTABL : ML preconditioning not available.\n");
3269 #endif
3270 break;
3271
3272 case HYMLMAXWELL :
3273 printf("BiCGSTABL : MLMaxwell preconditioning not available.\n");
3274 break;
3275
3276 case HYMLI :
3277 #ifdef HAVE_MLI
3278 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3279 printf("MLI preconditioning \n");
3280 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3281 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
3282 HYPRE_DummyFunction, HYPrecon_);
3283 else
3284 {
3285 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
3286 HYPRE_LSI_MLISetup, HYPrecon_);
3287 HYPreconSetup_ = 1;
3288 }
3289 #else
3290 printf("BiCGSTABL : ML preconditioning not available.\n");
3291 #endif
3292 break;
3293
3294 case HYUZAWA :
3295 printf("BiCGSTABL : Uzawa preconditioning not available.\n");
3296 exit(1);
3297 break;
3298
3299 case HYAMS :
3300 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3301 printf("AMS preconditioning\n");
3302 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3303 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_,HYPRE_AMSSolve,
3304 HYPRE_DummyFunction, HYPrecon_);
3305 else
3306 {
3307 setupPreconAMS();
3308 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_AMSSolve,
3309 HYPRE_AMSSetup, HYPrecon_);
3310 HYPreconSetup_ = 1;
3311 }
3312 break;
3313
3314 case HYSYSPDE :
3315 #ifdef HY_SYSPDE
3316 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3317 printf("SysPDe preconditioning\n");
3318 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3319 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3320 HYPRE_DummyFunction, HYPrecon_);
3321 else
3322 {
3323 setupPreconSysPDE();
3324 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3325 HYPRE_ParCSRSysPDESetup, HYPrecon_);
3326 HYPreconSetup_ = 1;
3327 }
3328 #else
3329 printf("BiCGSTABL : SysPDe preconditioning not available.\n");
3330 #endif
3331 break;
3332
3333 case HYDSLU :
3334 #ifdef HYPRE_USING_DSUPERLU
3335 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3336 printf("DSuperLU preconditioning\n");
3337 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3338 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3339 HYPRE_DummyFunction, HYPrecon_);
3340 else
3341 {
3342 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
3343 HYPRE_ParCSRBiCGSTABLSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3344 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
3345 HYPreconSetup_ = 1;
3346 }
3347 #else
3348 printf("BiCGSTABL : DSuperLU preconditioning not available.\n");
3349 #endif
3350 break;
3351 }
3352 return;
3353 }
3354
3355 //***************************************************************************
3356 // set up preconditioners for TFQMR
3357 //---------------------------------------------------------------------------
3358
setupTFQmrPrecon()3359 void HYPRE_LinSysCore::setupTFQmrPrecon()
3360 {
3361 //-------------------------------------------------------------------
3362 // if matrix has been reloaded, reset preconditioner
3363 //-------------------------------------------------------------------
3364
3365 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
3366 selectPreconditioner( HYPreconName_ );
3367
3368 //-------------------------------------------------------------------
3369 // set up preconditioners
3370 //-------------------------------------------------------------------
3371
3372 switch ( HYPreconID_ )
3373 {
3374 case HYIDENTITY :
3375 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3376 printf("No preconditioning \n");
3377 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
3378 HYPRE_DummyFunction, HYPrecon_);
3379 break;
3380
3381 case HYDIAGONAL :
3382 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3383 printf("Diagonal preconditioning \n");
3384 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3385 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3386 HYPRE_DummyFunction, HYPrecon_);
3387 else
3388 {
3389 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3390 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
3391 HYPreconSetup_ = 1;
3392 }
3393 break;
3394
3395 case HYPILUT :
3396 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3397 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
3398 HYPRE_DummyFunction, HYPrecon_);
3399 else
3400 {
3401 setupPreconPILUT();
3402 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
3403 HYPRE_ParCSRPilutSetup, HYPrecon_);
3404 HYPreconSetup_ = 1;
3405 }
3406 break;
3407
3408 case HYDDILUT :
3409 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3410 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3411 HYPRE_DummyFunction, HYPrecon_);
3412 else
3413 {
3414 setupPreconDDILUT();
3415 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3416 HYPRE_LSI_DDIlutSetup, HYPrecon_);
3417 HYPreconSetup_ = 1;
3418 }
3419 break;
3420
3421 case HYDDICT :
3422 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3423 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3424 HYPRE_DummyFunction, HYPrecon_);
3425 else
3426 {
3427 setupPreconDDICT();
3428 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3429 HYPRE_LSI_DDICTSetup, HYPrecon_);
3430 HYPreconSetup_ = 1;
3431 }
3432 break;
3433
3434 case HYSCHWARZ :
3435 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3436 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
3437 HYPRE_DummyFunction, HYPrecon_);
3438 else
3439 {
3440 setupPreconSchwarz();
3441 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
3442 HYPRE_LSI_SchwarzSetup, HYPrecon_);
3443 HYPreconSetup_ = 1;
3444 }
3445 break;
3446
3447 case HYPOLY :
3448 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3449 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3450 HYPRE_DummyFunction, HYPrecon_);
3451 else
3452 {
3453 setupPreconPoly();
3454 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3455 HYPRE_LSI_PolySetup, HYPrecon_);
3456 HYPreconSetup_ = 1;
3457 }
3458 break;
3459
3460 case HYPARASAILS :
3461 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3462 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3463 HYPRE_DummyFunction, HYPrecon_);
3464 else
3465 {
3466 setupPreconParaSails();
3467 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3468 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
3469 HYPreconSetup_ = 1;
3470 }
3471 break;
3472
3473 case HYBOOMERAMG :
3474 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3475 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
3476 HYPRE_DummyFunction, HYPrecon_);
3477 else
3478 {
3479 setupPreconBoomerAMG();
3480 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
3481 HYPRE_BoomerAMGSetup, HYPrecon_);
3482 HYPreconSetup_ = 1;
3483 }
3484 break;
3485
3486 case HYEUCLID :
3487 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3488 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_EuclidSolve,
3489 HYPRE_DummyFunction, HYPrecon_);
3490 else
3491 {
3492 setupPreconEuclid();
3493 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_EuclidSolve,
3494 HYPRE_EuclidSetup, HYPrecon_);
3495 HYPreconSetup_ = 1;
3496 }
3497 break;
3498
3499 case HYBLOCK :
3500 printf("TFQMR : block preconditioning not available.\n");
3501 exit(1);
3502 break;
3503
3504 case HYML :
3505 #ifdef HAVE_ML
3506 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3507 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
3508 HYPRE_DummyFunction, HYPrecon_);
3509 else
3510 {
3511 setupPreconML();
3512 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
3513 HYPRE_LSI_MLSetup, HYPrecon_);
3514 HYPreconSetup_ = 1;
3515 }
3516 #else
3517 printf("TFQMR : ML preconditioning not available.\n");
3518 #endif
3519 break;
3520
3521 case HYMLMAXWELL :
3522 printf("TFQMR : MLMaxwell preconditioning not available.\n");
3523 break;
3524
3525 case HYMLI :
3526 #ifdef HAVE_MLI
3527 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3528 printf("MLI preconditioning \n");
3529 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3530 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
3531 HYPRE_DummyFunction, HYPrecon_);
3532 else
3533 {
3534 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
3535 HYPRE_LSI_MLISetup, HYPrecon_);
3536 HYPreconSetup_ = 1;
3537 }
3538 #else
3539 printf("TFQMR : MLI preconditioning not available.\n");
3540 #endif
3541 break;
3542
3543 case HYUZAWA :
3544 printf("TFQMR : Uzawa preconditioning not available.\n");
3545 exit(1);
3546 break;
3547
3548 case HYAMS :
3549 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3550 printf("AMS preconditioning\n");
3551 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3552 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_,HYPRE_AMSSolve,
3553 HYPRE_DummyFunction, HYPrecon_);
3554 else
3555 {
3556 setupPreconAMS();
3557 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_AMSSolve,
3558 HYPRE_AMSSetup, HYPrecon_);
3559 HYPreconSetup_ = 1;
3560 }
3561 break;
3562
3563 case HYSYSPDE :
3564 #ifdef HY_SYSPDE
3565 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3566 printf("SysPDe preconditioning\n");
3567 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3568 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3569 HYPRE_DummyFunction, HYPrecon_);
3570 else
3571 {
3572 setupPreconSysPDE();
3573 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3574 HYPRE_ParCSRSysPDESetup, HYPrecon_);
3575 HYPreconSetup_ = 1;
3576 }
3577 #else
3578 printf("TFQMR : SysPDe preconditioning not available.\n");
3579 #endif
3580 break;
3581
3582 case HYDSLU :
3583 #ifdef HYPRE_USING_DSUPERLU
3584 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3585 printf("DSuperLU preconditioning\n");
3586 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3587 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3588 HYPRE_DummyFunction, HYPrecon_);
3589 else
3590 {
3591 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
3592 HYPRE_ParCSRTFQmrSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3593 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
3594 HYPreconSetup_ = 1;
3595 }
3596 #else
3597 printf("TFQMR : DSuperLU preconditioning not available.\n");
3598 #endif
3599 break;
3600 }
3601 return;
3602 }
3603
3604 //***************************************************************************
3605 // set up preconditioners for BiCGS
3606 //---------------------------------------------------------------------------
3607
setupBiCGSPrecon()3608 void HYPRE_LinSysCore::setupBiCGSPrecon()
3609 {
3610 //-------------------------------------------------------------------
3611 // if matrix has been reloaded, reset preconditioner
3612 //-------------------------------------------------------------------
3613
3614 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
3615 selectPreconditioner( HYPreconName_ );
3616
3617 //-------------------------------------------------------------------
3618 // set up preconditioners
3619 //-------------------------------------------------------------------
3620
3621 switch ( HYPreconID_ )
3622 {
3623 case HYIDENTITY :
3624 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3625 printf("No preconditioning \n");
3626 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
3627 HYPRE_DummyFunction, HYPrecon_);
3628 break;
3629
3630 case HYDIAGONAL :
3631 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3632 printf("Diagonal preconditioning \n");
3633 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3634 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3635 HYPRE_DummyFunction, HYPrecon_);
3636 else
3637 {
3638 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3639 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
3640 HYPreconSetup_ = 1;
3641 }
3642 break;
3643
3644 case HYPILUT :
3645 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3646 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
3647 HYPRE_DummyFunction, HYPrecon_);
3648 else
3649 {
3650 setupPreconPILUT();
3651 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRPilutSolve,
3652 HYPRE_ParCSRPilutSetup, HYPrecon_);
3653 HYPreconSetup_ = 1;
3654 }
3655 break;
3656
3657 case HYDDILUT :
3658 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3659 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3660 HYPRE_DummyFunction, HYPrecon_);
3661 else
3662 {
3663 setupPreconDDILUT();
3664 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DDIlutSolve,
3665 HYPRE_LSI_DDIlutSetup, HYPrecon_);
3666 HYPreconSetup_ = 1;
3667 }
3668 break;
3669
3670 case HYDDICT :
3671 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3672 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3673 HYPRE_DummyFunction, HYPrecon_);
3674 else
3675 {
3676 setupPreconDDICT();
3677 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3678 HYPRE_LSI_DDICTSetup, HYPrecon_);
3679 HYPreconSetup_ = 1;
3680 }
3681 break;
3682
3683 case HYSCHWARZ :
3684 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3685 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
3686 HYPRE_DummyFunction, HYPrecon_);
3687 else
3688 {
3689 setupPreconSchwarz();
3690 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_SchwarzSolve,
3691 HYPRE_LSI_SchwarzSetup, HYPrecon_);
3692 HYPreconSetup_ = 1;
3693 }
3694 break;
3695
3696 case HYPOLY :
3697 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3698 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3699 HYPRE_DummyFunction, HYPrecon_);
3700 else
3701 {
3702 setupPreconPoly();
3703 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3704 HYPRE_LSI_PolySetup, HYPrecon_);
3705 HYPreconSetup_ = 1;
3706 }
3707 break;
3708
3709 case HYPARASAILS :
3710 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3711 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3712 HYPRE_DummyFunction, HYPrecon_);
3713 else
3714 {
3715 setupPreconParaSails();
3716 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3717 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
3718 HYPreconSetup_ = 1;
3719 }
3720 break;
3721
3722 case HYBOOMERAMG :
3723 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3724 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
3725 HYPRE_DummyFunction, HYPrecon_);
3726 else
3727 {
3728 setupPreconBoomerAMG();
3729 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
3730 HYPRE_BoomerAMGSetup, HYPrecon_);
3731 HYPreconSetup_ = 1;
3732 }
3733 break;
3734
3735 case HYEUCLID :
3736 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3737 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_EuclidSolve,
3738 HYPRE_DummyFunction, HYPrecon_);
3739 else
3740 {
3741 setupPreconEuclid();
3742 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_EuclidSolve,
3743 HYPRE_EuclidSetup, HYPrecon_);
3744 HYPreconSetup_ = 1;
3745 }
3746 break;
3747
3748 case HYBLOCK :
3749 printf("BiCGS : block preconditioning not available.\n");
3750 exit(1);
3751 break;
3752
3753 case HYML :
3754 #ifdef HAVE_ML
3755 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3756 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
3757 HYPRE_DummyFunction, HYPrecon_);
3758 else
3759 {
3760 setupPreconML();
3761 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
3762 HYPRE_LSI_MLSetup, HYPrecon_);
3763 HYPreconSetup_ = 1;
3764 }
3765 #else
3766 printf("BiCGS : ML preconditioning not available.\n");
3767 #endif
3768 break;
3769
3770 case HYMLMAXWELL :
3771 printf("BiCGS : MLMaxwell preconditioning not available.\n");
3772 break;
3773
3774 case HYMLI :
3775 #ifdef HAVE_MLI
3776 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3777 printf("MLI preconditioning \n");
3778 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3779 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
3780 HYPRE_DummyFunction, HYPrecon_);
3781 else
3782 {
3783 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
3784 HYPRE_LSI_MLISetup, HYPrecon_);
3785 HYPreconSetup_ = 1;
3786 }
3787 #else
3788 printf("BiCGS : MLI preconditioning not available.\n");
3789 #endif
3790 break;
3791
3792 case HYUZAWA :
3793 printf("BiCGS : Uzawa preconditioning not available.\n");
3794 exit(1);
3795 break;
3796
3797 case HYAMS :
3798 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3799 printf("AMS preconditioning\n");
3800 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3801 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_,HYPRE_AMSSolve,
3802 HYPRE_DummyFunction, HYPrecon_);
3803 else
3804 {
3805 setupPreconAMS();
3806 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_AMSSolve,
3807 HYPRE_AMSSetup, HYPrecon_);
3808 HYPreconSetup_ = 1;
3809 }
3810 break;
3811
3812 case HYSYSPDE :
3813 #ifdef HY_SYSPDE
3814 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3815 printf("SysPDe preconditioning\n");
3816 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3817 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3818 HYPRE_DummyFunction, HYPrecon_);
3819 else
3820 {
3821 setupPreconSysPDE();
3822 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
3823 HYPRE_ParCSRSysPDESetup, HYPrecon_);
3824 HYPreconSetup_ = 1;
3825 }
3826 #else
3827 printf("BiCGS : SysPDe preconditioning not available.\n");
3828 #endif
3829 break;
3830
3831 case HYDSLU :
3832 #ifdef HYPRE_USING_DSUPERLU
3833 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
3834 printf("DSuperLU preconditioning\n");
3835 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3836 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3837 HYPRE_DummyFunction, HYPrecon_);
3838 else
3839 {
3840 HYPRE_LSI_DSuperLUSetOutputLevel(HYPrecon_, HYOutputLevel_);
3841 HYPRE_ParCSRBiCGSSetPrecond(HYSolver_, HYPRE_LSI_DSuperLUSolve,
3842 HYPRE_LSI_DSuperLUSetup, HYPrecon_);
3843 HYPreconSetup_ = 1;
3844 }
3845 #else
3846 printf("BiCGS : DSuperLU preconditioning not available.\n");
3847 #endif
3848 break;
3849 }
3850 return;
3851 }
3852
3853 //***************************************************************************
3854 // set up preconditioners for Symmetric QMR
3855 //---------------------------------------------------------------------------
3856
setupSymQMRPrecon()3857 void HYPRE_LinSysCore::setupSymQMRPrecon()
3858 {
3859 //-------------------------------------------------------------------
3860 // if matrix has been reloaded, reset preconditioner
3861 //-------------------------------------------------------------------
3862
3863 if ( HYPreconReuse_ == 0 && HYPreconSetup_ == 1 )
3864 selectPreconditioner( HYPreconName_ );
3865
3866 //-------------------------------------------------------------------
3867 // set up preconditioners
3868 //-------------------------------------------------------------------
3869
3870 switch ( HYPreconID_ )
3871 {
3872 case HYIDENTITY :
3873 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3874 printf("No preconditioning \n");
3875 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_SolveIdentity,
3876 HYPRE_DummyFunction, HYPrecon_);
3877 break;
3878
3879 case HYDIAGONAL :
3880 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
3881 printf("Diagonal preconditioning \n");
3882 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3883 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3884 HYPRE_DummyFunction, HYPrecon_);
3885 else
3886 {
3887 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_ParCSRDiagScale,
3888 HYPRE_ParCSRDiagScaleSetup,HYPrecon_);
3889 HYPreconSetup_ = 1;
3890 }
3891 break;
3892
3893 case HYPILUT :
3894 printf("ERROR : PILUT does not match SymQMR in general.\n");
3895 exit(1);
3896 break;
3897
3898 case HYDDILUT :
3899 printf("ERROR : DDILUT does not match SymQMR in general.\n");
3900 exit(1);
3901 break;
3902
3903 case HYDDICT :
3904 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3905 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3906 HYPRE_DummyFunction, HYPrecon_);
3907 else
3908 {
3909 setupPreconDDICT();
3910 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_DDICTSolve,
3911 HYPRE_LSI_DDICTSetup, HYPrecon_);
3912 HYPreconSetup_ = 1;
3913 }
3914 break;
3915
3916 case HYSCHWARZ :
3917 printf("ERROR : Schwarz does not match SymQMR in general.\n");
3918 exit(1);
3919 break;
3920
3921 case HYPOLY :
3922 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3923 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3924 HYPRE_DummyFunction, HYPrecon_);
3925 else
3926 {
3927 setupPreconPoly();
3928 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_PolySolve,
3929 HYPRE_LSI_PolySetup, HYPrecon_);
3930 HYPreconSetup_ = 1;
3931 }
3932 break;
3933
3934 case HYPARASAILS :
3935 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
3936 HYPRE_ParCSRParaSailsSetLogging(HYPrecon_, 1);
3937 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3938 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3939 HYPRE_DummyFunction, HYPrecon_);
3940 else
3941 {
3942 setupPreconParaSails();
3943 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_ParCSRParaSailsSolve,
3944 HYPRE_ParCSRParaSailsSetup,HYPrecon_);
3945 HYPreconSetup_ = 1;
3946 }
3947 break;
3948
3949 case HYBOOMERAMG :
3950 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3951 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_BoomerAMGSolve,
3952 HYPRE_DummyFunction, HYPrecon_);
3953 else
3954 {
3955 setupPreconBoomerAMG();
3956 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_BoomerAMGSolve,
3957 HYPRE_BoomerAMGSetup, HYPrecon_);
3958 HYPreconSetup_ = 1;
3959 }
3960 break;
3961
3962 case HYEUCLID :
3963 printf("ERROR : Euclid does not match SymQMR in general.\n");
3964 exit(1);
3965 break;
3966
3967 case HYBLOCK :
3968 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3969 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,
3970 HYPRE_LSI_BlockPrecondSolve, HYPRE_DummyFunction,
3971 HYPrecon_);
3972 else
3973 {
3974 setupPreconBlock();
3975 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,
3976 HYPRE_LSI_BlockPrecondSolve,
3977 HYPRE_LSI_BlockPrecondSetup, HYPrecon_);
3978 HYPreconSetup_ = 1;
3979 }
3980 break;
3981
3982 case HYML :
3983 #ifdef HAVE_ML
3984 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
3985 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_MLSolve,
3986 HYPRE_DummyFunction, HYPrecon_);
3987 else
3988 {
3989 setupPreconML();
3990 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_LSI_MLSolve,
3991 HYPRE_LSI_MLSetup, HYPrecon_);
3992 HYPreconSetup_ = 1;
3993 }
3994 #else
3995 printf("SymQMR : ML preconditioning not available.\n");
3996 #endif
3997 break;
3998
3999 case HYMLMAXWELL :
4000 printf("SymQMR : MLMaxwell preconditioning not available.\n");
4001 break;
4002
4003 case HYMLI :
4004 #ifdef HAVE_MLI
4005 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4006 printf("MLI preconditioning \n");
4007 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
4008 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_LSI_MLISolve,
4009 HYPRE_DummyFunction, HYPrecon_);
4010 else
4011 {
4012 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_LSI_MLISolve,
4013 HYPRE_LSI_MLISetup, HYPrecon_);
4014 HYPreconSetup_ = 1;
4015 }
4016 #else
4017 printf("SymQMR : MLI preconditioning not available.\n");
4018 #endif
4019 break;
4020
4021 case HYUZAWA :
4022 printf("SymQMR : Uzawa preconditioning not available.\n");
4023 exit(1);
4024 break;
4025
4026 case HYAMS :
4027 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4028 printf("AMS preconditioning\n");
4029 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
4030 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_,HYPRE_AMSSolve,
4031 HYPRE_DummyFunction, HYPrecon_);
4032 else
4033 {
4034 setupPreconAMS();
4035 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_AMSSolve,
4036 HYPRE_AMSSetup, HYPrecon_);
4037 HYPreconSetup_ = 1;
4038 }
4039 break;
4040
4041 case HYSYSPDE :
4042 #ifdef HY_SYSPDE
4043 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4044 printf("SysPDe preconditioning\n");
4045 if ( HYPreconReuse_ == 1 && HYPreconSetup_ == 1 )
4046 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
4047 HYPRE_DummyFunction, HYPrecon_);
4048 else
4049 {
4050 setupPreconSysPDE();
4051 HYPRE_ParCSRSymQMRSetPrecond(HYSolver_, HYPRE_ParCSRSysPDESolve,
4052 HYPRE_ParCSRSysPDESetup, HYPrecon_);
4053 HYPreconSetup_ = 1;
4054 }
4055 #else
4056 printf("SymQMR : SysPDe preconditioning not available.\n");
4057 #endif
4058 break;
4059
4060 case HYDSLU :
4061 printf("BiCGS : DSuperLU preconditioning not an option.\n");
4062 break;
4063 }
4064 return;
4065 }
4066
4067 //***************************************************************************
4068 // this function sets up BOOMERAMG preconditioner
4069 //---------------------------------------------------------------------------
4070
setupPreconBoomerAMG()4071 void HYPRE_LinSysCore::setupPreconBoomerAMG()
4072 {
4073 int i, j, *num_sweeps, *relax_type, **relax_points;
4074 double *relax_wt, *relax_omega;
4075
4076 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4077 {
4078 printf("AMG max levels = %d\n", amgMaxLevels_);
4079 printf("AMG coarsen type = %d\n", amgCoarsenType_);
4080 printf("AMG measure type = %d\n", amgMeasureType_);
4081 printf("AMG threshold = %e\n", amgStrongThreshold_);
4082 printf("AMG numsweeps = %d\n", amgNumSweeps_[0]);
4083 printf("AMG relax type = %d\n", amgRelaxType_[0]);
4084 if (amgGridRlxType_) printf("AMG CF smoothing \n");
4085 printf("AMG relax weight = %e\n", amgRelaxWeight_[0]);
4086 printf("AMG relax omega = %e\n", amgRelaxOmega_[0]);
4087 printf("AMG system size = %d\n", amgSystemSize_);
4088 printf("AMG smooth type = %d\n", amgSmoothType_);
4089 printf("AMG smooth numlevels = %d\n", amgSmoothNumLevels_);
4090 printf("AMG smooth numsweeps = %d\n", amgSmoothNumSweeps_);
4091 printf("AMG Schwarz variant = %d\n", amgSchwarzVariant_);
4092 printf("AMG Schwarz overlap = %d\n", amgSchwarzOverlap_);
4093 printf("AMG Schwarz domain type = %d\n", amgSchwarzDomainType_);
4094 printf("AMG Schwarz relax weight = %e\n", amgSchwarzRelaxWt_);
4095 }
4096 if ( HYOutputLevel_ & HYFEI_AMGDEBUG )
4097 {
4098 HYPRE_BoomerAMGSetDebugFlag(HYPrecon_, 0);
4099 HYPRE_BoomerAMGSetPrintLevel(HYPrecon_, 1);
4100 }
4101 if ( amgSystemSize_ > 1 )
4102 HYPRE_BoomerAMGSetNumFunctions(HYPrecon_, amgSystemSize_);
4103 HYPRE_BoomerAMGSetMaxLevels(HYPrecon_, amgMaxLevels_);
4104 HYPRE_BoomerAMGSetCoarsenType(HYPrecon_, amgCoarsenType_);
4105 HYPRE_BoomerAMGSetMeasureType(HYPrecon_, amgMeasureType_);
4106 HYPRE_BoomerAMGSetStrongThreshold(HYPrecon_,amgStrongThreshold_);
4107 HYPRE_BoomerAMGSetTol(HYPrecon_, 0.0e0);
4108 HYPRE_BoomerAMGSetMaxIter(HYPrecon_, 1);
4109 num_sweeps = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
4110 for ( i = 0; i < 4; i++ ) num_sweeps[i] = amgNumSweeps_[i];
4111
4112 HYPRE_BoomerAMGSetNumGridSweeps(HYPrecon_, num_sweeps);
4113 relax_type = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
4114 for ( i = 0; i < 4; i++ ) relax_type[i] = amgRelaxType_[i];
4115
4116 HYPRE_BoomerAMGSetGridRelaxType(HYPrecon_, relax_type);
4117 relax_wt = hypre_CTAlloc(double,amgMaxLevels_,HYPRE_MEMORY_HOST);
4118 for ( i = 0; i < amgMaxLevels_; i++ ) relax_wt[i] = amgRelaxWeight_[i];
4119 HYPRE_BoomerAMGSetRelaxWeight(HYPrecon_, relax_wt);
4120
4121 relax_omega = hypre_CTAlloc(double,amgMaxLevels_,HYPRE_MEMORY_HOST);
4122 for ( i = 0; i < amgMaxLevels_; i++ ) relax_omega[i] = amgRelaxOmega_[i];
4123 HYPRE_BoomerAMGSetOmega(HYPrecon_, relax_omega);
4124
4125 if (amgGridRlxType_)
4126 {
4127 relax_points = hypre_CTAlloc(int*,4,HYPRE_MEMORY_HOST);
4128 relax_points[0] = hypre_CTAlloc(int,num_sweeps[0],HYPRE_MEMORY_HOST);
4129 for ( j = 0; j < num_sweeps[0]; j++ ) relax_points[0][j] = 0;
4130 relax_points[1] = hypre_CTAlloc(int,2*num_sweeps[1],HYPRE_MEMORY_HOST);
4131 for ( j = 0; j < num_sweeps[1]; j+=2 )
4132 {relax_points[1][j] = -1;relax_points[1][j+1] = 1;}
4133 relax_points[2] = hypre_CTAlloc(int,2*num_sweeps[2],HYPRE_MEMORY_HOST);
4134 for ( j = 0; j < num_sweeps[2]; j+=2 )
4135 {relax_points[2][j] = -1;relax_points[2][j+1] = 1;}
4136 relax_points[3] = hypre_CTAlloc(int,num_sweeps[3],HYPRE_MEMORY_HOST);
4137 for ( j = 0; j < num_sweeps[3]; j++ ) relax_points[3][j] = 0;
4138 }
4139 else
4140 {
4141 relax_points = hypre_CTAlloc(int*,4,HYPRE_MEMORY_HOST);
4142 for ( i = 0; i < 4; i++ )
4143 {
4144 relax_points[i] = hypre_CTAlloc(int,num_sweeps[i],HYPRE_MEMORY_HOST);
4145 for ( j = 0; j < num_sweeps[i]; j++ ) relax_points[i][j] = 0;
4146 }
4147 }
4148 HYPRE_BoomerAMGSetGridRelaxPoints(HYPrecon_, relax_points);
4149
4150 if ( amgSmoothNumLevels_ > 0 )
4151 {
4152 HYPRE_BoomerAMGSetSmoothType(HYPrecon_, amgSmoothType_);
4153 HYPRE_BoomerAMGSetSmoothNumLevels(HYPrecon_, amgSmoothNumLevels_);
4154 HYPRE_BoomerAMGSetSmoothNumSweeps(HYPrecon_, amgSmoothNumSweeps_);
4155 HYPRE_BoomerAMGSetSchwarzRlxWeight(HYPrecon_, amgSchwarzRelaxWt_);
4156 HYPRE_BoomerAMGSetVariant(HYPrecon_, amgSchwarzVariant_);
4157 HYPRE_BoomerAMGSetOverlap(HYPrecon_, amgSchwarzOverlap_);
4158 HYPRE_BoomerAMGSetDomainType(HYPrecon_, amgSchwarzDomainType_);
4159 }
4160
4161 if ( amgUseGSMG_ == 1 )
4162 {
4163 HYPRE_BoomerAMGSetGSMG(HYPrecon_, 4);
4164 HYPRE_BoomerAMGSetNumSamples(HYPrecon_,amgGSMGNSamples_);
4165 }
4166
4167 HYPRE_BoomerAMGSetAggNumLevels(HYPrecon_, amgAggLevels_);
4168 HYPRE_BoomerAMGSetInterpType(HYPrecon_, amgInterpType_);
4169 HYPRE_BoomerAMGSetPMaxElmts(HYPrecon_, amgPmax_);
4170 }
4171
4172 //***************************************************************************
4173 // this function sets up ML preconditioner
4174 //---------------------------------------------------------------------------
4175
setupPreconML()4176 void HYPRE_LinSysCore::setupPreconML()
4177 {
4178 #ifdef HAVE_ML
4179 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4180 {
4181 printf("ML strong threshold = %e\n", mlStrongThreshold_);
4182 printf("ML numsweeps(pre) = %d\n", mlNumPreSweeps_);
4183 printf("ML numsweeps(post) = %d\n", mlNumPostSweeps_);
4184 printf("ML smoother (pre) = %d\n", mlPresmootherType_);
4185 printf("ML smoother (post) = %d\n", mlPostsmootherType_);
4186 printf("ML relax weight = %e\n", mlRelaxWeight_);
4187 }
4188 HYPRE_LSI_MLSetMethod(HYPrecon_,mlMethod_);
4189 HYPRE_LSI_MLSetCoarseSolver(HYPrecon_,mlCoarseSolver_);
4190 HYPRE_LSI_MLSetCoarsenScheme(HYPrecon_,mlCoarsenScheme_);
4191 HYPRE_LSI_MLSetStrongThreshold(HYPrecon_,mlStrongThreshold_);
4192 HYPRE_LSI_MLSetNumPreSmoothings(HYPrecon_,mlNumPreSweeps_);
4193 HYPRE_LSI_MLSetNumPostSmoothings(HYPrecon_,mlNumPostSweeps_);
4194 HYPRE_LSI_MLSetPreSmoother(HYPrecon_,mlPresmootherType_);
4195 HYPRE_LSI_MLSetPostSmoother(HYPrecon_,mlPostsmootherType_);
4196 HYPRE_LSI_MLSetDampingFactor(HYPrecon_,mlRelaxWeight_);
4197 HYPRE_LSI_MLSetNumPDEs(HYPrecon_,mlNumPDEs_);
4198 #else
4199 return;
4200 #endif
4201 }
4202
4203 //***************************************************************************
4204 // this function sets up MLMaxwell preconditioner
4205 //---------------------------------------------------------------------------
4206
setupPreconMLMaxwell()4207 void HYPRE_LinSysCore::setupPreconMLMaxwell()
4208 {
4209 #ifdef HAVE_MLMAXWELL
4210 HYPRE_ParCSRMatrix A_csr;
4211
4212 if (maxwellGEN_ != NULL)
4213 HYPRE_LSI_MLMaxwellSetGMatrix(HYPrecon_,maxwellGEN_);
4214 else
4215 {
4216 printf("HYPRE_LSC::setupPreconMLMaxwell ERROR - no G matrix.\n");
4217 exit(1);
4218 }
4219 if (maxwellANN_ == NULL)
4220 {
4221 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4222 hypre_BoomerAMGBuildCoarseOperator((hypre_ParCSRMatrix *) maxwellGEN_,
4223 (hypre_ParCSRMatrix *) A_csr,
4224 (hypre_ParCSRMatrix *) maxwellGEN_,
4225 (hypre_ParCSRMatrix **) &maxwellANN_);
4226 }
4227 HYPRE_LSI_MLMaxwellSetANNMatrix(HYPrecon_,maxwellANN_);
4228 #else
4229 return;
4230 #endif
4231 }
4232
4233 //***************************************************************************
4234 // this function sets up AMS preconditioner
4235 //---------------------------------------------------------------------------
4236
setupPreconAMS()4237 void HYPRE_LinSysCore::setupPreconAMS()
4238 {
4239 #if 0
4240 HYPRE_ParVector parVecX, parVecY, parVecZ;
4241 #endif
4242
4243 /* Set AMS parameters */
4244 HYPRE_AMSSetDimension(HYPrecon_, amsNumPDEs_);
4245 HYPRE_AMSSetMaxIter(HYPrecon_, amsMaxIter_);
4246 HYPRE_AMSSetTol(HYPrecon_, amsTol_);
4247 HYPRE_AMSSetCycleType(HYPrecon_, amsCycleType_);
4248 HYPRE_AMSSetPrintLevel(HYPrecon_, amsPrintLevel_);
4249 HYPRE_AMSSetSmoothingOptions(HYPrecon_, amsRelaxType_, amsRelaxTimes_,
4250 amsRelaxWt_, amsRelaxOmega_);
4251
4252 if (amsBetaPoisson_ != NULL)
4253 HYPRE_AMSSetBetaPoissonMatrix(HYPrecon_, amsBetaPoisson_);
4254
4255 HYPRE_AMSSetAlphaAMGOptions(HYPrecon_, amsAlphaCoarsenType_,
4256 amsAlphaAggLevels_, amsAlphaRelaxType_, amsAlphaStrengthThresh_,
4257 amsAlphaInterpType_, amsAlphaPmax_);
4258 HYPRE_AMSSetBetaAMGOptions(HYPrecon_, amsBetaCoarsenType_,
4259 amsBetaAggLevels_, amsBetaRelaxType_, amsBetaStrengthThresh_,
4260 amsBetaInterpType_, amsBetaPmax_);
4261
4262 #if 0
4263 if (maxwellGEN_ != NULL)
4264 HYPRE_AMSSetDiscreteGradient(HYPrecon_, maxwellGEN_);
4265 else
4266 {
4267 printf("HYPRE_LSC::setupPreconAMS ERROR - no G matrix.\n");
4268 exit(1);
4269 }
4270 if (amsX_ == NULL && amsY_ != NULL)
4271 {
4272 HYPRE_IJVectorGetObject(amsX_, (void **) &parVecX);
4273 HYPRE_IJVectorGetObject(amsY_, (void **) &parVecY);
4274 HYPRE_IJVectorGetObject(amsZ_, (void **) &parVecZ);
4275 HYPRE_AMSSetCoordinateVectors(HYPrecon_,parVecX,parVecY,parVecZ);
4276 }
4277 #endif
4278
4279 // Call AMS to construct the discrete gradient matrix G
4280 // and the nodal coordinate vectors
4281 {
4282 HYPRE_ParCSRMatrix A_csr;
4283 HYPRE_ParVector b_csr;
4284 HYPRE_ParVector x_csr;
4285
4286 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4287 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
4288 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
4289
4290 if( amsG_ == NULL ) {
4291
4292 //Old way of doing things
4293 //only works for 1 domain per processor (in ALE3D)
4294 //not compatible with contact
4295 HYPRE_AMSFEISetup(HYPrecon_,
4296 A_csr,
4297 b_csr,
4298 x_csr,
4299 AMSData_.EdgeNodeList_,
4300 AMSData_.NodeNumbers_,
4301 AMSData_.numEdges_,
4302 AMSData_.numLocalNodes_,
4303 AMSData_.numNodes_,
4304 AMSData_.NodalCoord_);
4305 } else {
4306 //New Code//
4307 HYPRE_ParCSRMatrix G_csr;
4308 HYPRE_ParVector X_csr;
4309 HYPRE_ParVector Y_csr;
4310 HYPRE_ParVector Z_csr;
4311 HYPRE_IJMatrixGetObject(amsG_, (void **) &G_csr);
4312 HYPRE_IJVectorGetObject(amsX_, (void **) &X_csr);
4313 HYPRE_IJVectorGetObject(amsY_, (void **) &Y_csr);
4314 HYPRE_IJVectorGetObject(amsZ_, (void **) &Z_csr);
4315 HYPRE_AMSSetCoordinateVectors(HYPrecon_,X_csr,Y_csr,Z_csr);
4316 bool debugprint = false;
4317 if( debugprint ) {
4318 HYPRE_ParCSRMatrixPrint( G_csr, "G.parcsr" );
4319 HYPRE_ParCSRMatrixPrint( A_csr, "A.parcsr" );
4320 HYPRE_ParVectorPrint( b_csr, "B.parvector" );
4321 HYPRE_ParVectorPrint( X_csr, "X.parvector" );
4322 HYPRE_ParVectorPrint( Y_csr, "Y.parvector" );
4323 HYPRE_ParVectorPrint( Z_csr, "Z.parvector" );
4324 }
4325 HYPRE_AMSSetDiscreteGradient(HYPrecon_,G_csr);
4326 }
4327 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4328 printf("AMSprecon: finished building auxiliary info, calling AMSSetup\n");
4329 //int ierr = HYPRE_AMSSetup(HYPrecon_,A_csr,b_csr,x_csr);
4330 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4331 printf("AMSprecon: finished with AMSSetup\n");
4332 }
4333
4334 }
4335
4336 //***************************************************************************
4337 // this function sets up DDICT preconditioner
4338 //---------------------------------------------------------------------------
4339
setupPreconDDICT()4340 void HYPRE_LinSysCore::setupPreconDDICT()
4341 {
4342 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4343 {
4344 printf("DDICT - fillin = %e\n", ddictFillin_);
4345 printf("DDICT - drop tol = %e\n", ddictDropTol_);
4346 }
4347 if ( HYOutputLevel_ & HYFEI_DDILUT )
4348 HYPRE_LSI_DDICTSetOutputLevel(HYPrecon_,2);
4349 HYPRE_LSI_DDICTSetFillin(HYPrecon_,ddictFillin_);
4350 HYPRE_LSI_DDICTSetDropTolerance(HYPrecon_,ddictDropTol_);
4351 }
4352
4353 //***************************************************************************
4354 // this function sets up DDILUT preconditioner
4355 //---------------------------------------------------------------------------
4356
setupPreconDDILUT()4357 void HYPRE_LinSysCore::setupPreconDDILUT()
4358 {
4359 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4360 {
4361 printf("DDILUT - fillin = %e\n", ddilutFillin_);
4362 printf("DDILUT - drop tol = %e\n", ddilutDropTol_);
4363 }
4364 if ( HYOutputLevel_ & HYFEI_DDILUT )
4365 HYPRE_LSI_DDIlutSetOutputLevel(HYPrecon_,2);
4366 if ( ddilutReorder_ ) HYPRE_LSI_DDIlutSetReorder(HYPrecon_);
4367 HYPRE_LSI_DDIlutSetFillin(HYPrecon_,ddilutFillin_);
4368 HYPRE_LSI_DDIlutSetDropTolerance(HYPrecon_,ddilutDropTol_);
4369 if ( ddilutOverlap_ == 1 ) HYPRE_LSI_DDIlutSetOverlap(HYPrecon_);
4370 if ( ddilutReorder_ == 1 ) HYPRE_LSI_DDIlutSetReorder(HYPrecon_);
4371 }
4372
4373 //***************************************************************************
4374 // this function sets up Schwarz preconditioner
4375 //---------------------------------------------------------------------------
4376
setupPreconSchwarz()4377 void HYPRE_LinSysCore::setupPreconSchwarz()
4378 {
4379 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4380 {
4381 printf("Schwarz - ILU fillin = %e\n", schwarzFillin_);
4382 printf("Schwarz - nBlocks = %d\n", schwarzNblocks_);
4383 printf("Schwarz - blockSize = %d\n", schwarzBlksize_);
4384 }
4385 if ( HYOutputLevel_ & HYFEI_DDILUT )
4386 HYPRE_LSI_SchwarzSetOutputLevel(HYPrecon_,2);
4387 HYPRE_LSI_SchwarzSetILUTFillin(HYPrecon_,schwarzFillin_);
4388 HYPRE_LSI_SchwarzSetNBlocks(HYPrecon_, schwarzNblocks_);
4389 HYPRE_LSI_SchwarzSetBlockSize(HYPrecon_, schwarzBlksize_);
4390 }
4391
4392 //***************************************************************************
4393 // this function sets up Polynomial preconditioner
4394 //---------------------------------------------------------------------------
4395
setupPreconPoly()4396 void HYPRE_LinSysCore::setupPreconPoly()
4397 {
4398 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4399 printf("Polynomial preconditioning - order = %d\n",polyOrder_ );
4400 HYPRE_LSI_PolySetOrder(HYPrecon_, polyOrder_);
4401 }
4402
4403 //***************************************************************************
4404 // this function sets up ParaSails preconditioner
4405 //---------------------------------------------------------------------------
4406
setupPreconParaSails()4407 void HYPRE_LinSysCore::setupPreconParaSails()
4408 {
4409 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4410 {
4411 printf("ParaSails - nlevels = %d\n",parasailsNlevels_);
4412 printf("ParaSails - threshold = %e\n",parasailsThreshold_);
4413 printf("ParaSails - filter = %e\n",parasailsFilter_);
4414 printf("ParaSails - sym = %d\n",parasailsSym_);
4415 printf("ParaSails - loadbal = %e\n",parasailsLoadbal_);
4416 }
4417 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
4418 HYPRE_ParCSRParaSailsSetLogging(HYPrecon_, 1);
4419 HYPRE_ParCSRParaSailsSetSym(HYPrecon_,parasailsSym_);
4420 HYPRE_ParCSRParaSailsSetParams(HYPrecon_, parasailsThreshold_,
4421 parasailsNlevels_);
4422 HYPRE_ParCSRParaSailsSetFilter(HYPrecon_, parasailsFilter_);
4423 HYPRE_ParCSRParaSailsSetLoadbal(HYPrecon_, parasailsLoadbal_);
4424 HYPRE_ParCSRParaSailsSetReuse(HYPrecon_, parasailsReuse_);
4425 }
4426
4427 //***************************************************************************
4428 // this function sets up Euclid preconditioner
4429 //---------------------------------------------------------------------------
4430
setupPreconEuclid()4431 void HYPRE_LinSysCore::setupPreconEuclid()
4432 {
4433 if ((HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0)
4434 {
4435 for ( int i = 0; i < euclidargc_; i++ )
4436 printf("Euclid parameter : %s %s\n", euclidargv_[2*i],
4437 euclidargv_[2*i+1]);
4438 }
4439 HYPRE_EuclidSetParams(HYPrecon_,euclidargc_*2,euclidargv_);
4440 }
4441
4442 //***************************************************************************
4443 // this function sets up Pilut preconditioner
4444 //---------------------------------------------------------------------------
4445
setupPreconPILUT()4446 void HYPRE_LinSysCore::setupPreconPILUT()
4447 {
4448 if (pilutFillin_ == 0) pilutFillin_ = pilutMaxNnzPerRow_;
4449 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4450 {
4451 printf("PILUT - row size = %d\n", pilutFillin_);
4452 printf("PILUT - drop tol = %e\n", pilutDropTol_);
4453 }
4454 HYPRE_ParCSRPilutSetFactorRowSize(HYPrecon_,pilutFillin_);
4455 HYPRE_ParCSRPilutSetDropTolerance(HYPrecon_,pilutDropTol_);
4456 }
4457
4458 //***************************************************************************
4459 // this function sets up block preconditioner
4460 //---------------------------------------------------------------------------
4461
setupPreconBlock()4462 void HYPRE_LinSysCore::setupPreconBlock()
4463 {
4464 HYPRE_Lookup *newLookup;
4465 newLookup = hypre_TAlloc(HYPRE_Lookup, 1, HYPRE_MEMORY_HOST);
4466 newLookup->object = (void *) lookup_;
4467 HYPRE_LSI_BlockPrecondSetLookup( HYPrecon_, newLookup );
4468 free( newLookup );
4469 }
4470
4471 //***************************************************************************
4472 // this function sets up system pde preconditioner
4473 //---------------------------------------------------------------------------
4474
setupPreconSysPDE()4475 void HYPRE_LinSysCore::setupPreconSysPDE()
4476 {
4477 #ifdef HAVE_SYSPDE
4478 if (sysPDEMethod_ >= 0)
4479 HYPRE_ParCSRSysPDESetMethod(HYPrecon_, sysPDEMethod_);
4480 if (sysPDEFormat_ >= 0)
4481 HYPRE_ParCSRFormat(HYPrecon_, sysPDEFormat_);
4482 if (sysPDETol_ > 0.0)
4483 HYPRE_ParCSRSysPDESetTol(HYPrecon_, sysPDETol_);
4484 if (sysPDEMaxIter_ > 0)
4485 HYPRE_ParCSRSysPDESetMaxIter(HYPrecon_, sysPDEMaxIter_);
4486 if (sysPDENumPre_ > 0)
4487 HYPRE_ParCSRSysPDESetNPre(HYPrecon_, sysPDENumPre_);
4488 if (sysPDENumPost_ > 0)
4489 HYPRE_ParCSRSysPDESetNPost(HYPrecon_, sysPDENumPost_);
4490 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
4491 HYPRE_ParCSRSysPDESetLogging(HYPrecon_, 1);
4492 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
4493 HYPRE_ParCSRSysPDESetPrintLevel(HYPrecon_, 1);
4494 #endif
4495 }
4496
4497 //***************************************************************************
4498 // this function solve the incoming linear system using Boomeramg
4499 //---------------------------------------------------------------------------
4500
solveUsingBoomeramg(int & status)4501 void HYPRE_LinSysCore::solveUsingBoomeramg(int& status)
4502 {
4503 int i, j, *relax_type, *num_sweeps, **relax_points;
4504 double *relax_wt, *relax_omega;
4505 HYPRE_ParCSRMatrix A_csr;
4506 HYPRE_ParVector b_csr;
4507 HYPRE_ParVector x_csr;
4508
4509 //-------------------------------------------------------------------
4510 // get matrix and vectors in ParCSR format
4511 //-------------------------------------------------------------------
4512
4513 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4514 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
4515 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
4516
4517 //-------------------------------------------------------------------
4518 // set BoomerAMG parameters
4519 //-------------------------------------------------------------------
4520
4521 HYPRE_BoomerAMGSetCoarsenType(HYSolver_, amgCoarsenType_);
4522 HYPRE_BoomerAMGSetMeasureType(HYSolver_, amgMeasureType_);
4523 HYPRE_BoomerAMGSetStrongThreshold(HYSolver_, amgStrongThreshold_);
4524
4525 num_sweeps = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
4526 for ( i = 0; i < 4; i++ ) num_sweeps[i] = amgNumSweeps_[i];
4527 HYPRE_BoomerAMGSetNumGridSweeps(HYSolver_, num_sweeps);
4528
4529 relax_type = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
4530 for ( i = 0; i < 4; i++ ) relax_type[i] = amgRelaxType_[i];
4531 HYPRE_BoomerAMGSetGridRelaxType(HYSolver_, relax_type);
4532
4533 HYPRE_BoomerAMGSetMaxLevels(HYPrecon_, amgMaxLevels_);
4534 relax_wt = hypre_CTAlloc(double, amgMaxLevels_,HYPRE_MEMORY_HOST);
4535 for ( i = 0; i < amgMaxLevels_; i++ ) relax_wt[i] = amgRelaxWeight_[i];
4536 HYPRE_BoomerAMGSetRelaxWeight(HYSolver_, relax_wt);
4537
4538 relax_omega = hypre_CTAlloc(double, amgMaxLevels_,HYPRE_MEMORY_HOST);
4539 for ( i = 0; i < amgMaxLevels_; i++ ) relax_omega[i] = amgRelaxOmega_[i];
4540 HYPRE_BoomerAMGSetOmega(HYPrecon_, relax_omega);
4541
4542 relax_points = hypre_CTAlloc(int*,4,HYPRE_MEMORY_HOST);
4543 for ( i = 0; i < 4; i++ )
4544 {
4545 relax_points[i] = hypre_CTAlloc(int,num_sweeps[i],HYPRE_MEMORY_HOST);
4546 for ( j = 0; j < num_sweeps[i]; j++ ) relax_points[i][j] = 0;
4547 }
4548 HYPRE_BoomerAMGSetGridRelaxPoints(HYPrecon_, relax_points);
4549 if ( amgSmoothNumLevels_ > 0 )
4550 {
4551 HYPRE_BoomerAMGSetSmoothType(HYPrecon_, amgSmoothType_);
4552 HYPRE_BoomerAMGSetSmoothNumLevels(HYPrecon_, amgSmoothNumLevels_);
4553 HYPRE_BoomerAMGSetSmoothNumSweeps(HYPrecon_, amgSmoothNumSweeps_);
4554 HYPRE_BoomerAMGSetSchwarzRlxWeight(HYPrecon_, amgSchwarzRelaxWt_);
4555 HYPRE_BoomerAMGSetVariant(HYPrecon_, amgSchwarzVariant_);
4556 HYPRE_BoomerAMGSetOverlap(HYPrecon_, amgSchwarzOverlap_);
4557 HYPRE_BoomerAMGSetDomainType(HYPrecon_, amgSchwarzDomainType_);
4558 }
4559
4560 if ( amgUseGSMG_ == 1 )
4561 {
4562 HYPRE_BoomerAMGSetGSMG(HYPrecon_, 4);
4563 HYPRE_BoomerAMGSetNumSamples(HYPrecon_,amgGSMGNSamples_);
4564 }
4565 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 && mypid_ == 0 )
4566 {
4567 printf("***************************************************\n");
4568 printf("* Boomeramg (AMG) solver \n");
4569 printf("* coarsen type = %d\n", amgCoarsenType_);
4570 printf("* measure type = %d\n", amgMeasureType_);
4571 printf("* threshold = %e\n", amgStrongThreshold_);
4572 printf("* numsweeps = %d\n", amgNumSweeps_[0]);
4573 printf("* relax type = %d\n", amgRelaxType_[0]);
4574 printf("* relax weight = %e\n", amgRelaxWeight_[0]);
4575 printf("* maximum iterations = %d\n", maxIterations_);
4576 printf("* smooth type = %d\n", amgSmoothType_);
4577 printf("* smooth numlevels = %d\n", amgSmoothNumLevels_);
4578 printf("* smooth numsweeps = %d\n", amgSmoothNumSweeps_);
4579 printf("* Schwarz variant = %d\n", amgSchwarzVariant_);
4580 printf("* Schwarz overlap = %d\n", amgSchwarzOverlap_);
4581 printf("* Schwarz domain type = %d\n", amgSchwarzDomainType_);
4582 printf("* Schwarz relax weight = %e\n", amgSchwarzRelaxWt_);
4583 printf("* convergence tolerance = %e\n", tolerance_);
4584 printf("*--------------------------------------------------\n");
4585 }
4586 if ( HYOutputLevel_ & HYFEI_AMGDEBUG )
4587 {
4588 HYPRE_BoomerAMGSetDebugFlag(HYSolver_, 0);
4589 HYPRE_BoomerAMGSetPrintLevel(HYSolver_, 1);
4590 }
4591 HYPRE_BoomerAMGSetMaxIter(HYSolver_, maxIterations_);
4592 HYPRE_BoomerAMGSetMeasureType(HYSolver_, 0);
4593 HYPRE_BoomerAMGSetup( HYSolver_, A_csr, b_csr, x_csr );
4594
4595 //-------------------------------------------------------------------
4596 // BoomerAMG solve
4597 //-------------------------------------------------------------------
4598
4599 HYPRE_BoomerAMGSolve( HYSolver_, A_csr, b_csr, x_csr );
4600
4601 status = 0;
4602 }
4603
4604 //***************************************************************************
4605 // this function solve the incoming linear system using SuperLU
4606 //---------------------------------------------------------------------------
4607
solveUsingSuperLU(int & status)4608 double HYPRE_LinSysCore::solveUsingSuperLU(int& status)
4609 {
4610 double rnorm=-1.0;
4611 #ifdef HAVE_SUPERLU
4612 int i, nnz, nrows, ierr;
4613 int rowSize, *colInd, *new_ia, *new_ja, *ind_array;
4614 int nz_ptr, *partition, start_row, end_row;
4615 double *colVal, *new_a;
4616 HYPRE_ParCSRMatrix A_csr;
4617 HYPRE_ParVector r_csr;
4618 HYPRE_ParVector b_csr;
4619 HYPRE_ParVector x_csr;
4620
4621 int info=0, permc_spec;
4622 int *perm_r, *perm_c;
4623 double *rhs, *soln;
4624 superlu_options_t slu_options;
4625 SuperLUStat_t slu_stat;
4626 SuperMatrix A2, B, L, U;
4627 NRformat *Ustore;
4628 SCformat *Lstore;
4629
4630 //-------------------------------------------------------------------
4631 // available for sequential processing only for now
4632 //-------------------------------------------------------------------
4633
4634 if ( numProcs_ > 1 )
4635 {
4636 printf("solveUsingSuperLU ERROR - too many processors.\n");
4637 status = -1;
4638 return rnorm;
4639 }
4640
4641 //-------------------------------------------------------------------
4642 // need to construct a CSR matrix, and the column indices should
4643 // have been stored in colIndices and rowLengths
4644 //-------------------------------------------------------------------
4645
4646 if ( localStartRow_ != 1 )
4647 {
4648 printf("solveUsingSuperLU ERROR - row does not start at 1\n");
4649 status = -1;
4650 return rnorm;
4651 }
4652
4653 //-------------------------------------------------------------------
4654 // get information about the current matrix
4655 //-------------------------------------------------------------------
4656
4657 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
4658 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );
4659 start_row = partition[0];
4660 end_row = partition[1] - 1;
4661 nrows = end_row - start_row + 1;
4662 free( partition );
4663
4664 nnz = 0;
4665 for ( i = start_row; i <= end_row; i++ )
4666 {
4667 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
4668 nnz += rowSize;
4669 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
4670 }
4671
4672 new_ia = new int[nrows+1];
4673 new_ja = new int[nnz];
4674 new_a = new double[nnz];
4675 nz_ptr = HYPRE_LSI_GetParCSRMatrix(currA_,nrows,nnz,new_ia,new_ja,new_a);
4676 nnz = nz_ptr;
4677
4678 //-------------------------------------------------------------------
4679 // set up SuperLU CSR matrix and the corresponding rhs
4680 //-------------------------------------------------------------------
4681
4682 dCreate_CompRow_Matrix(&A2,nrows,nrows,nnz,new_a,new_ja,new_ia,
4683 SLU_NR,SLU_D,SLU_GE);
4684 ind_array = new int[nrows];
4685 for ( i = 0; i < nrows; i++ ) ind_array[i] = i;
4686 rhs = new double[nrows];
4687
4688 ierr = HYPRE_IJVectorGetValues(currB_, nrows, ind_array, rhs);
4689
4690 hypre_assert(!ierr);
4691 dCreate_Dense_Matrix(&B, nrows, 1, rhs, nrows, SLU_DN, SLU_D, SLU_GE);
4692
4693 //-------------------------------------------------------------------
4694 // set up the rest and solve (permc_spec=0 : natural ordering)
4695 //-------------------------------------------------------------------
4696
4697 perm_r = new int[nrows];
4698 perm_c = new int[nrows];
4699 permc_spec = superluOrdering_;
4700 get_perm_c(permc_spec, &A2, perm_c);
4701 for ( i = 0; i < nrows; i++ ) perm_r[i] = 0;
4702
4703 set_default_options(&slu_options);
4704 slu_options.ColPerm = MY_PERMC;
4705 slu_options.Fact = DOFACT;
4706 StatInit(&slu_stat);
4707 dgssv(&slu_options, &A2, perm_c, perm_r, &L, &U, &B, &slu_stat, &info);
4708
4709 //-------------------------------------------------------------------
4710 // postprocessing of the return status information
4711 //-------------------------------------------------------------------
4712
4713 if ( info == 0 )
4714 {
4715 status = 1;
4716 Lstore = (SCformat *) L.Store;
4717 Ustore = (NRformat *) U.Store;
4718 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4719 {
4720 printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
4721 printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
4722 printf("SuperLU : NNZ in L+U = %d\n",Lstore->nnz+Ustore->nnz-nrows);
4723 }
4724 }
4725 else
4726 {
4727 status = 0;
4728 printf("HYPRE_LinSysCore::solveUsingSuperLU - dgssv error = %d\n",info);
4729 }
4730
4731 //-------------------------------------------------------------------
4732 // fetch the solution and find residual norm
4733 //-------------------------------------------------------------------
4734
4735 if ( info == 0 )
4736 {
4737 soln = (double *) ((DNformat *) B.Store)->nzval;
4738 ierr = HYPRE_IJVectorSetValues(currX_, nrows, (const int *) ind_array,
4739 (const double *) soln);
4740 hypre_assert(!ierr);
4741
4742 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
4743 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
4744 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
4745
4746 ierr = HYPRE_ParVectorCopy( b_csr, r_csr );
4747 hypre_assert(!ierr);
4748 HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
4749 ierr = HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
4750 hypre_assert(!ierr);
4751 rnorm = sqrt( rnorm );
4752 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
4753 printf("HYPRE_LSC::solveUsingSuperLU - FINAL NORM = %e.\n",rnorm);
4754 }
4755
4756 //-------------------------------------------------------------------
4757 // clean up
4758 //-------------------------------------------------------------------
4759
4760 delete [] ind_array;
4761 delete [] rhs;
4762 delete [] perm_c;
4763 delete [] perm_r;
4764 delete [] new_ia;
4765 delete [] new_ja;
4766 delete [] new_a;
4767 Destroy_SuperMatrix_Store(&B);
4768 Destroy_SuperNode_Matrix(&L);
4769 SUPERLU_FREE( A2.Store );
4770 SUPERLU_FREE( ((NRformat *) U.Store)->colind);
4771 SUPERLU_FREE( ((NRformat *) U.Store)->rowptr);
4772 SUPERLU_FREE( ((NRformat *) U.Store)->nzval);
4773 SUPERLU_FREE( U.Store );
4774 StatFree(&slu_stat);
4775 #else
4776 status = -1;
4777 printf("HYPRE_LSC::solveUsingSuperLU : not available.\n");
4778 #endif
4779 return rnorm;
4780 }
4781
4782 //***************************************************************************
4783 // this function solve the incoming linear system using SuperLU
4784 // using expert mode
4785 //---------------------------------------------------------------------------
4786
solveUsingSuperLUX(int & status)4787 double HYPRE_LinSysCore::solveUsingSuperLUX(int& status)
4788 {
4789 double rnorm=-1.0;
4790 #ifdef HAVE_SUPERLU
4791 int i, nnz, nrows, ierr;
4792 int rowSize, *colInd, *new_ia, *new_ja, *ind_array;
4793 int nz_ptr;
4794 int *partition, start_row, end_row;
4795 double *colVal, *new_a;
4796 HYPRE_ParCSRMatrix A_csr;
4797 HYPRE_ParVector r_csr;
4798 HYPRE_ParVector b_csr;
4799 HYPRE_ParVector x_csr;
4800
4801 int info, permc_spec;
4802 int *perm_r, *perm_c, *etree, lwork;
4803 double *rhs, *soln, *sol2;
4804 double *R, *C;
4805 double *ferr, *berr;
4806 double rpg, rcond;
4807 void *work=NULL;
4808 char equed[1];
4809 GlobalLU_t Glu;
4810 mem_usage_t mem_usage;
4811 superlu_options_t slu_options;
4812 SuperLUStat_t slu_stat;
4813 SuperMatrix A2, B, X, L, U;
4814 NRformat *Ustore;
4815 SCformat *Lstore;
4816
4817 //-------------------------------------------------------------------
4818 // available for sequential processing only for now
4819 //-------------------------------------------------------------------
4820
4821 if ( numProcs_ > 1 )
4822 {
4823 printf("solveUsingSuperLUX ERROR - too many processors.\n");
4824 status = -1;
4825 return rnorm;
4826 }
4827
4828 //-------------------------------------------------------------------
4829 // need to construct a CSR matrix, and the column indices should
4830 // have been stored in colIndices and rowLengths
4831 //-------------------------------------------------------------------
4832
4833 if ( localStartRow_ != 1 )
4834 {
4835 printf("solveUsingSuperLUX ERROR - row not start at 1\n");
4836 status = -1;
4837 return rnorm;
4838 }
4839
4840 //-------------------------------------------------------------------
4841 // get information about the current matrix
4842 //-------------------------------------------------------------------
4843
4844 HYPRE_IJMatrixGetObject(currA_, (void**) &A_csr);
4845 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );
4846 start_row = partition[0];
4847 end_row = partition[1] - 1;
4848 nrows = end_row - start_row + 1;
4849 free( partition );
4850
4851 nnz = 0;
4852 for ( i = 0; i < nrows; i++ )
4853 {
4854 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
4855 nnz += rowSize;
4856 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
4857 }
4858
4859 new_ia = new int[nrows+1];
4860 new_ja = new int[nnz];
4861 new_a = new double[nnz];
4862 nz_ptr = HYPRE_LSI_GetParCSRMatrix(currA_,nrows,nnz,new_ia,new_ja,new_a);
4863 nnz = nz_ptr;
4864
4865 //-------------------------------------------------------------------
4866 // set up SuperLU CSR matrix and the corresponding rhs
4867 //-------------------------------------------------------------------
4868
4869 dCreate_CompRow_Matrix(&A2,nrows,nrows,nnz,new_a,new_ja,new_ia,SLU_NR,
4870 SLU_D,SLU_GE);
4871 ind_array = new int[nrows];
4872 for ( i = 0; i < nrows; i++ ) ind_array[i] = i;
4873
4874 rhs = new double[nrows];
4875 ierr = HYPRE_IJVectorGetValues(currB_, nrows, ind_array, rhs);
4876 hypre_assert(!ierr);
4877 dCreate_Dense_Matrix(&B, nrows, 1, rhs, nrows, SLU_DN, SLU_D, SLU_GE);
4878
4879 soln = new double[nrows];
4880 for ( i = 0; i < nrows; i++ ) soln[i] = 0.0;
4881 dCreate_Dense_Matrix(&X, nrows, 1, soln, nrows, SLU_DN, SLU_D, SLU_GE);
4882
4883 //-------------------------------------------------------------------
4884 // set up the other parameters (permc_spec=0 : natural ordering)
4885 //-------------------------------------------------------------------
4886
4887 perm_r = new int[nrows];
4888 for ( i = 0; i < nrows; i++ ) perm_r[i] = 0;
4889 perm_c = new int[nrows];
4890 etree = new int[nrows];
4891 permc_spec = superluOrdering_;
4892 get_perm_c(permc_spec, &A2, perm_c);
4893 lwork = 0;
4894 set_default_options(&slu_options);
4895 slu_options.ColPerm = MY_PERMC;
4896 slu_options.Equil = YES;
4897 slu_options.Trans = NOTRANS;
4898 slu_options.Fact = DOFACT;
4899 slu_options.IterRefine = SLU_DOUBLE;
4900 slu_options.DiagPivotThresh = 1.0;
4901 slu_options.PivotGrowth = YES;
4902 slu_options.ConditionNumber = YES;
4903
4904 StatInit(&slu_stat);
4905 *equed = 'N';
4906 R = (double *) SUPERLU_MALLOC(A2.nrow * sizeof(double));
4907 C = (double *) SUPERLU_MALLOC(A2.ncol * sizeof(double));
4908 ferr = (double *) SUPERLU_MALLOC(sizeof(double));
4909 berr = (double *) SUPERLU_MALLOC(sizeof(double));
4910
4911 //-------------------------------------------------------------------
4912 // solve
4913 //-------------------------------------------------------------------
4914
4915 // dgssvx(&slu_options, &A2, perm_c, perm_r, etree,
4916 // equed, R, C, &L, &U, work, lwork, &B, &X,
4917 // &rpg, &rcond, ferr, berr, &mem_usage, &slu_stat, &info);
4918 dgssvx(&slu_options, &A2, perm_c, perm_r, etree,
4919 equed, R, C, &L, &U, work, lwork, &B, &X,
4920 &rpg, &rcond, ferr, berr, &Glu, &mem_usage, &slu_stat, &info);
4921
4922 //-------------------------------------------------------------------
4923 // print SuperLU internal information at the first step
4924 //-------------------------------------------------------------------
4925
4926 if ( info == 0 || info == nrows+1 )
4927 {
4928 status = 1;
4929 Lstore = (SCformat *) L.Store;
4930 Ustore = (NRformat *) U.Store;
4931 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
4932 {
4933 printf("Recip. pivot growth = %e\n", rpg);
4934 printf("%8s%16s%16s\n", "rhs", "FERR", "BERR");
4935 printf("%8d%16e%16e\n", 1, ferr[0], berr[0]);
4936 if ( rcond != 0.0 )
4937 printf(" SuperLU : condition number = %e\n", 1.0/rcond);
4938 else
4939 printf(" SuperLU : Recip. condition number = %e\n", rcond);
4940 printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
4941 printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
4942 printf("SuperLUX : NNZ in L+U = %d\n", Lstore->nnz+Ustore->nnz-nrows);
4943 }
4944 }
4945 else
4946 {
4947 printf("solveUsingSuperLUX - dgssvx error code = %d\n",info);
4948 status = 0;
4949 }
4950
4951 //-------------------------------------------------------------------
4952 // fetch the solution and find residual norm
4953 //-------------------------------------------------------------------
4954
4955 if ( status == 1 )
4956 {
4957 sol2 = (double *) ((DNformat *) X.Store)->nzval;
4958
4959 ierr = HYPRE_IJVectorSetValues(currX_, nrows, (const int *) ind_array,
4960 (const double *) sol2);
4961 hypre_assert(!ierr);
4962
4963 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
4964 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
4965 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
4966 ierr = HYPRE_ParVectorCopy( b_csr, r_csr );
4967 hypre_assert(!ierr);
4968 ierr = HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
4969 hypre_assert(!ierr);
4970 ierr = HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
4971 hypre_assert(!ierr);
4972 rnorm = sqrt( rnorm );
4973 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 2 )
4974 printf("HYPRE_LSC::solveUsingSuperLUX - FINAL NORM = %e.\n",rnorm);
4975 }
4976
4977 //-------------------------------------------------------------------
4978 // clean up
4979 //-------------------------------------------------------------------
4980
4981 delete [] ind_array;
4982 delete [] perm_c;
4983 delete [] perm_r;
4984 delete [] etree;
4985 delete [] rhs;
4986 delete [] soln;
4987 delete [] new_ia;
4988 delete [] new_ja;
4989 delete [] new_a;
4990 Destroy_SuperMatrix_Store(&B);
4991 Destroy_SuperMatrix_Store(&X);
4992 Destroy_SuperNode_Matrix(&L);
4993 SUPERLU_FREE( A2.Store );
4994 SUPERLU_FREE( ((NRformat *) U.Store)->colind);
4995 SUPERLU_FREE( ((NRformat *) U.Store)->rowptr);
4996 SUPERLU_FREE( ((NRformat *) U.Store)->nzval);
4997 SUPERLU_FREE( U.Store );
4998 SUPERLU_FREE (R);
4999 SUPERLU_FREE (C);
5000 SUPERLU_FREE (ferr);
5001 SUPERLU_FREE (berr);
5002 StatFree(&slu_stat);
5003 #else
5004 status = -1;
5005 printf("HYPRE_LSC::solveUsingSuperLUX : not available.\n");
5006 #endif
5007 return rnorm;
5008 }
5009
5010 //***************************************************************************
5011 // this function solve the incoming linear system using distributed SuperLU
5012 //---------------------------------------------------------------------------
5013
solveUsingDSuperLU(int & status)5014 double HYPRE_LinSysCore::solveUsingDSuperLU(int& status)
5015 {
5016 double rnorm=1.0;
5017 #ifdef HYPRE_USING_DSUPERLU
5018 int ierr;
5019 HYPRE_ParCSRMatrix A_csr;
5020 HYPRE_ParVector x_csr, b_csr, r_csr;
5021
5022 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5023 HYPRE_IJVectorGetObject(currX_, (void **) &x_csr);
5024 HYPRE_IJVectorGetObject(currB_, (void **) &b_csr);
5025 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5026
5027 HYPRE_LSI_DSuperLUCreate(comm_, &HYSolver_);
5028 HYPRE_LSI_DSuperLUSetOutputLevel(HYSolver_, HYOutputLevel_);
5029 HYPRE_LSI_DSuperLUSetup(HYSolver_, A_csr, b_csr, x_csr);
5030 HYPRE_LSI_DSuperLUSolve(HYSolver_, A_csr, b_csr, x_csr);
5031 HYPRE_LSI_DSuperLUDestroy(HYSolver_);
5032 ierr = HYPRE_ParVectorCopy( b_csr, r_csr );
5033 hypre_assert(!ierr);
5034 ierr = HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5035 hypre_assert(!ierr);
5036 ierr = HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5037 hypre_assert(!ierr);
5038 rnorm = sqrt( rnorm );
5039 //if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5040 // printf("HYPRE_LSC::solveUsingDSuperLU - FINAL NORM = %e.\n",rnorm);
5041 #endif
5042 return rnorm;
5043 }
5044
5045 //***************************************************************************
5046 // this function solve the incoming linear system using Y12M
5047 //---------------------------------------------------------------------------
5048
solveUsingY12M(int & status)5049 void HYPRE_LinSysCore::solveUsingY12M(int& status)
5050 {
5051 #ifdef Y12M
5052 int i, k, nnz, nrows, ierr;
5053 int rowSize, *colInd, *ind_array;
5054 int j, nz_ptr, *colLengths, count, maxRowSize;
5055 double *colVal, rnorm;
5056 double upperSum, lowerSum, *accuSoln, *origRhs;
5057 HYPRE_ParCSRMatrix A_csr;
5058 HYPRE_ParVector r_csr;
5059 HYPRE_ParVector b_csr;
5060 HYPRE_ParVector x_csr;
5061
5062 int n, nn, nn1, *rnr, *snr, *ha, iha, iflag[10], ifail;
5063 double *pivot, *val, *rhs, aflag[8];
5064
5065 //-------------------------------------------------------------------
5066 // available for sequential processing only for now
5067 //-------------------------------------------------------------------
5068
5069 if ( numProcs_ > 1 )
5070 {
5071 printf("solveUsingY12M ERROR - too many processors.\n");
5072 status = 0;
5073 return;
5074 }
5075
5076 //-------------------------------------------------------------------
5077 // need to construct a CSR matrix, and the column indices should
5078 // have been stored in colIndices and rowLengths
5079 //-------------------------------------------------------------------
5080
5081 if ( localStartRow_ != 1 )
5082 {
5083 printf("solveUsingY12M ERROR - row does not start at 1.\n");
5084 status = -1;
5085 return;
5086 }
5087 if (slideReduction_ == 1)
5088 nrows = localEndRow_ - 2 * nConstraints_;
5089 else if (slideReduction_ == 2 || slideReduction_ == 3)
5090 nrows = localEndRow_ - nConstraints_;
5091 else if (schurReduction_ == 1)
5092 nrows = localEndRow_ - localStartRow_ + 1 - A21NRows_;
5093 else nrows = localEndRow_;
5094
5095 colLengths = new int[nrows];
5096 for ( i = 0; i < nrows; i++ ) colLengths[i] = 0;
5097
5098 maxRowSize = 0;
5099 HYPRE_IJMatrixGetObject(currA_, (void**) &A_csr);
5100
5101 for ( i = 0; i < nrows; i++ )
5102 {
5103 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
5104 maxRowSize = ( rowSize > maxRowSize ) ? rowSize : maxRowSize;
5105 for ( j = 0; j < rowSize; j++ )
5106 if ( colVal[j] != 0.0 ) colLengths[colInd[j]]++;
5107 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
5108 }
5109 nnz = 0;
5110 for ( i = 0; i < nrows; i++ ) nnz += colLengths[i];
5111
5112 nn = 2 * nnz;
5113 nn1 = 2 * nnz;
5114 snr = new int[nn];
5115 rnr = new int[nn1];
5116 val = new double[nn];
5117 pivot = new double[nrows];
5118 iha = nrows;
5119 ha = new int[iha*11];
5120
5121 nz_ptr = 0;
5122 for ( i = 0; i < nrows; i++ )
5123 {
5124 HYPRE_ParCSRMatrixGetRow(A_csr,i,&rowSize,&colInd,&colVal);
5125 for ( j = 0; j < rowSize; j++ )
5126 {
5127 if ( colVal[j] != 0.0 )
5128 {
5129 rnr[nz_ptr] = i + 1;
5130 snr[nz_ptr] = colInd[j] + 1;
5131 val[nz_ptr] = colVal[j];
5132 nz_ptr++;
5133 }
5134 }
5135 HYPRE_ParCSRMatrixRestoreRow(A_csr,i,&rowSize,&colInd,&colVal);
5136 }
5137
5138 nnz = nz_ptr;
5139
5140 //-------------------------------------------------------------------
5141 // set up other parameters and the right hand side
5142 //-------------------------------------------------------------------
5143
5144 aflag[0] = 16.0;
5145 aflag[1] = 0.0;
5146 aflag[2] = 1.0e8;
5147 aflag[3] = 1.0e-12;
5148 iflag[0] = 1;
5149 iflag[1] = 3;
5150 iflag[2] = 1;
5151 iflag[3] = 0;
5152 iflag[4] = 2;
5153 ind_array = new int[nrows];
5154 for ( i = 0; i < nrows; i++ ) ind_array[i] = i;
5155 rhs = new double[nrows];
5156
5157 ierr = HYPRE_IJVectorGetValues(currB_, nrows, ind_array, rhs);
5158 hypre_assert(!ierr);
5159
5160 //-------------------------------------------------------------------
5161 // call Y12M to solve the linear system
5162 //-------------------------------------------------------------------
5163
5164 y12maf_(&nrows,&nnz,val,snr,&nn,rnr,&nn1,pivot,ha,&iha,aflag,iflag,
5165 rhs,&ifail);
5166 if ( ifail != 0 && (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5167 {
5168 printf("solveUsingY12M WARNING - ifail = %d\n", ifail);
5169 }
5170
5171 //-------------------------------------------------------------------
5172 // postprocessing
5173 //-------------------------------------------------------------------
5174
5175 if ( ifail == 0 )
5176 {
5177 ierr = HYPRE_IJVectorSetValues(currX_, nrows, (const int *) &ind_array,
5178 (const double *) rhs);
5179 hypre_assert(!ierr);
5180
5181 HYPRE_IJVectorGetObject(currX_, (void**) &x_csr);
5182 HYPRE_IJVectorGetObject(currR_, (void**) &r_csr);
5183 HYPRE_IJVectorGetObject(currB_, (void**) &b_csr);
5184 ierr = HYPRE_ParVectorCopy( b_csr, r_csr );
5185 hypre_assert(!ierr);
5186 ierr = HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5187 hypre_assert(!ierr);
5188 ierr = HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5189 hypre_assert(!ierr);
5190 rnorm = sqrt( rnorm );
5191 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5192 printf("HYPRE_LSC::solveUsingY12M - final norm = %e.\n", rnorm);
5193 }
5194
5195 //-------------------------------------------------------------------
5196 // clean up
5197 //-------------------------------------------------------------------
5198
5199 delete [] ind_array;
5200 delete [] rhs;
5201 delete [] val;
5202 delete [] snr;
5203 delete [] rnr;
5204 delete [] ha;
5205 delete [] pivot;
5206 #else
5207 status = -1;
5208 printf("HYPRE_LSC::solveUsingY12M - not available.\n");
5209 #endif
5210 }
5211
5212 //***************************************************************************
5213 // this function solve the incoming linear system using Y12M
5214 //---------------------------------------------------------------------------
5215
solveUsingAMGe(int & iterations)5216 void HYPRE_LinSysCore::solveUsingAMGe(int &iterations)
5217 {
5218 #ifdef HAVE_AMGE
5219 int i, nrows, ierr, *ind_array, status;
5220 double rnorm, *rhs, *sol;
5221 HYPRE_ParCSRMatrix A_csr;
5222 HYPRE_ParVector r_csr;
5223 HYPRE_ParVector b_csr;
5224 HYPRE_ParVector x_csr;
5225
5226 //-------------------------------------------------------------------
5227 // available for sequential processing only for now
5228 //-------------------------------------------------------------------
5229
5230 if ( numProcs_ > 1 )
5231 {
5232 printf("solveUsingAMGE ERROR - too many processors.\n");
5233 iterations = 0;
5234 return;
5235 }
5236
5237 //-------------------------------------------------------------------
5238 // need to construct a CSR matrix, and the column indices should
5239 // have been stored in colIndices and rowLengths
5240 //-------------------------------------------------------------------
5241
5242 if ( localStartRow_ != 1 )
5243 {
5244 printf("solveUsingAMGe ERROR - row does not start at 1.\n");
5245 status = -1;
5246 return;
5247 }
5248 if (slideReduction_ == 1)
5249 nrows = localEndRow_ - 2 * nConstraints_;
5250 else if (slideReduction_ == 2 || slideReduction_ == 3)
5251 nrows = localEndRow_ - nConstraints_;
5252 else if (schurReduction_ == 1)
5253 nrows = localEndRow_ - localStartRow_ + 1 - A21NRows_;
5254 else nrows = localEndRow_;
5255
5256 //-------------------------------------------------------------------
5257 // set up the right hand side
5258 //-------------------------------------------------------------------
5259
5260 ind_array = new int[nrows];
5261 for ( i = 0; i < nrows; i++ ) ind_array[i] = i;
5262 rhs = new double[nrows];
5263
5264 ierr = HYPRE_IJVectorGetValues(currB_, nrows, ind_array, rhs);
5265 hypre_assert(!ierr);
5266
5267 //-------------------------------------------------------------------
5268 // call Y12M to solve the linear system
5269 //-------------------------------------------------------------------
5270
5271 sol = new double[nrows];
5272 status = HYPRE_LSI_AMGeSolve( rhs, sol );
5273
5274 //-------------------------------------------------------------------
5275 // postprocessing
5276 //-------------------------------------------------------------------
5277
5278 ierr = HYPRE_IJVectorSetValues(currX_, nrows, (const int *) &ind_array,
5279 (const double *) sol);
5280 hypre_assert(!ierr);
5281
5282 HYPRE_IJVectorGetObject(currX_, (void**) &x_csr);
5283 HYPRE_IJVectorGetObject(currR_, (void**) &r_csr);
5284 HYPRE_IJVectorGetObject(currB_, (void**) &b_csr);
5285
5286 ierr = HYPRE_ParVectorCopy( b_csr, r_csr );
5287 hypre_assert(!ierr);
5288 ierr = HYPRE_ParCSRMatrixMatvec( -1.0, A_csr, x_csr, 1.0, r_csr );
5289 hypre_assert(!ierr);
5290 ierr = HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
5291 hypre_assert(!ierr);
5292 rnorm = sqrt( rnorm );
5293 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 1 )
5294 printf("HYPRE_LSC::solveUsingAMGe - final norm = %e.\n", rnorm);
5295
5296 //-------------------------------------------------------------------
5297 // clean up
5298 //-------------------------------------------------------------------
5299
5300 delete [] ind_array;
5301 delete [] rhs;
5302 delete [] sol;
5303 #else
5304 iterations = 0;
5305 printf("HYPRE_LSC::solveUsingAMGe - not available.\n");
5306 #endif
5307 }
5308
5309 //***************************************************************************
5310 // this function loads in the constraint numbers for reduction
5311 // (to activate automatic slave search, constrList should be NULL)
5312 //---------------------------------------------------------------------------
5313
loadConstraintNumbers(int nConstr,int * constrList)5314 void HYPRE_LinSysCore::loadConstraintNumbers(int nConstr, int *constrList)
5315 {
5316 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
5317 printf("%4d : HYPRE_LSC::loadConstraintNumbers - size = %d\n",
5318 mypid_, nConstr);
5319 nConstraints_ = nConstr;
5320 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
5321 printf("%4d : HYPRE_LSC::leaving loadConstraintNumbers\n", mypid_);
5322 }
5323
5324 //***************************************************************************
5325 // this function extracts the the version number from HYPRE
5326 //---------------------------------------------------------------------------
5327
getVersion()5328 char *HYPRE_LinSysCore::getVersion()
5329 {
5330 static char extVersion[100];
5331 char hypre[200], hypreVersion[50], ctmp[50];
5332 sprintf(hypre, "%s", HYPRE_VERSION);
5333 sscanf(hypre, "%s %s", ctmp, hypreVersion);
5334 sprintf(extVersion, "%s-%s", HYPRE_FEI_Version(), hypreVersion);
5335 return extVersion;
5336 }
5337
5338 //***************************************************************************
5339 // create a node to equation map from the solution vector
5340 //---------------------------------------------------------------------------
5341
beginCreateMapFromSoln()5342 void HYPRE_LinSysCore::beginCreateMapFromSoln()
5343 {
5344 mapFromSolnFlag_ = 1;
5345 mapFromSolnLengMax_ = 10;
5346 mapFromSolnLeng_ = 0;
5347 mapFromSolnList_ = new int[mapFromSolnLengMax_];
5348 mapFromSolnList2_ = new int[mapFromSolnLengMax_];
5349 return;
5350 }
5351
5352 //***************************************************************************
5353 // create a node to equation map from the solution vector
5354 //---------------------------------------------------------------------------
5355
endCreateMapFromSoln()5356 void HYPRE_LinSysCore::endCreateMapFromSoln()
5357 {
5358 int i, *iarray;
5359 double *darray;
5360
5361 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
5362 printf("%4d : HYPRE_LSC::entering endCreateMapFromSoln.\n",mypid_);
5363
5364 mapFromSolnFlag_ = 0;
5365 if ( mapFromSolnLeng_ > 0 )
5366 darray = new double[mapFromSolnLeng_];
5367 for ( i = 0; i < mapFromSolnLeng_; i++ )
5368 darray[i] = (double) mapFromSolnList_[i];
5369
5370 hypre_qsort1(mapFromSolnList2_, darray, 0, mapFromSolnLeng_-1);
5371 iarray = mapFromSolnList2_;
5372 mapFromSolnList2_ = mapFromSolnList_;
5373 mapFromSolnList_ = iarray;
5374 for ( i = 0; i < mapFromSolnLeng_; i++ )
5375 mapFromSolnList2_[i] = (int) darray[i];
5376 delete [] darray;
5377
5378 for ( i = 0; i < mapFromSolnLeng_; i++ )
5379 printf("HYPRE_LSC::mapFromSoln %d = %d\n",mapFromSolnList_[i],
5380 mapFromSolnList2_[i]);
5381
5382 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 3 )
5383 printf("%4d : HYPRE_LSC::leaving endCreateMapFromSoln.\n",mypid_);
5384 }
5385
5386 //***************************************************************************
5387 // add extra nonzero entries into the matrix data structure
5388 //---------------------------------------------------------------------------
5389
putIntoMappedMatrix(int row,int numValues,const double * values,const int * scatterIndices)5390 void HYPRE_LinSysCore::putIntoMappedMatrix(int row, int numValues,
5391 const double* values, const int* scatterIndices)
5392 {
5393 int i, index, colIndex, localRow, mappedRow, mappedCol, newLeng;
5394 int *tempInd, ind2;
5395 double *tempVal;
5396
5397 //-------------------------------------------------------------------
5398 // error checking
5399 //-------------------------------------------------------------------
5400
5401 if ( systemAssembled_ == 1 )
5402 {
5403 printf("putIntoMappedMatrix ERROR : matrix already assembled\n");
5404 exit(1);
5405 }
5406 if ( (row+1) < localStartRow_ || (row+1) > localEndRow_ )
5407 {
5408 printf("putIntoMappedMatrix ERROR : invalid row number %d.\n",row);
5409 exit(1);
5410 }
5411 index = HYPRE_LSI_Search(mapFromSolnList_, row, mapFromSolnLeng_);
5412
5413 if ( index >= 0 ) mappedRow = mapFromSolnList2_[index];
5414 else mappedRow = row;
5415 localRow = mappedRow - localStartRow_ + 1;
5416
5417 //-------------------------------------------------------------------
5418 // load the local matrix
5419 //-------------------------------------------------------------------
5420
5421 newLeng = rowLengths_[localRow] + numValues;
5422 tempInd = new int[newLeng];
5423 tempVal = new double[newLeng];
5424 for ( i = 0; i < rowLengths_[localRow]; i++ )
5425 {
5426 tempVal[i] = colValues_[localRow][i];
5427 tempInd[i] = colIndices_[localRow][i];
5428 }
5429 delete [] colValues_[localRow];
5430 delete [] colIndices_[localRow];
5431 colValues_[localRow] = tempVal;
5432 colIndices_[localRow] = tempInd;
5433
5434 index = rowLengths_[localRow];
5435
5436 for ( i = 0; i < numValues; i++ )
5437 {
5438 colIndex = scatterIndices[i];
5439
5440 ind2 = HYPRE_LSI_Search(mapFromSolnList_,colIndex,mapFromSolnLeng_);
5441 if ( mapFromSolnList_ != NULL ) mappedCol = mapFromSolnList2_[ind2];
5442 else mappedCol = colIndex;
5443
5444 ind2 = HYPRE_LSI_Search(colIndices_[localRow],mappedCol+1,index);
5445 if ( ind2 >= 0 )
5446 {
5447 newLeng--;
5448 colValues_[localRow][ind2] = values[i];
5449 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5450 printf("%4d : putIntoMappedMatrix (add) : row, col = %8d %8d %e \n",
5451 mypid_, localRow, colIndices_[localRow][ind2]-1,
5452 colValues_[localRow][ind2]);
5453 }
5454 else
5455 {
5456 ind2 = index;
5457 colIndices_[localRow][index] = mappedCol + 1;
5458 colValues_[localRow][index++] = values[i];
5459 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5460 printf("%4d : putIntoMappedMatrix : row, col = %8d %8d %e \n",
5461 mypid_, localRow, colIndices_[localRow][ind2]-1,
5462 colValues_[localRow][ind2]);
5463 hypre_qsort1(colIndices_[localRow],colValues_[localRow],0,index-1);
5464 }
5465 }
5466 rowLengths_[localRow] = newLeng;
5467 }
5468
5469 //***************************************************************************
5470 // project the initial guess into the previous solutions (x + X inv(R) Q^T b)
5471 // Given r and B (a collection of right hand vectors such that A X = B)
5472 //
5473 // min || r - B v ||
5474 // v
5475 //
5476 // = min (trans(r) r - trans(r) B v-trans(v) trans(B) r+trans(v) trans(B) B v)
5477 //
5478 // ==> trans(B) r = trans(B) B v ==> v = inv(trans(B) B) trans(B) r
5479 //
5480 // Once v is computed, x = x + X v
5481 //---------------------------------------------------------------------------
5482
computeMinResProjection(HYPRE_ParCSRMatrix A_csr,HYPRE_ParVector x_csr,HYPRE_ParVector b_csr)5483 void HYPRE_LinSysCore::computeMinResProjection(HYPRE_ParCSRMatrix A_csr,
5484 HYPRE_ParVector x_csr, HYPRE_ParVector b_csr)
5485 {
5486 int i;
5487 double alpha;
5488 HYPRE_ParVector r_csr, v_csr, w_csr;
5489
5490 //-----------------------------------------------------------------------
5491 // diagnostic message
5492 //-----------------------------------------------------------------------
5493
5494 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5495 printf("%4d : HYPRE_LSC::entering computeMinResProjection %d\n",mypid_,
5496 projectCurrSize_);
5497 if ( projectCurrSize_ == 0 && HYpxs_ == NULL ) return;
5498
5499 //-----------------------------------------------------------------------
5500 // compute r = b - A x, save Ax (to w)
5501 //-----------------------------------------------------------------------
5502
5503 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5504 HYPRE_IJVectorGetObject(HYpbs_[projectSize_], (void **) &w_csr);
5505 HYPRE_ParCSRMatrixMatvec( 1.0, A_csr, x_csr, 0.0, w_csr );
5506 HYPRE_ParVectorCopy( b_csr, r_csr );
5507 alpha = -1.0;
5508 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)w_csr,(hypre_ParVector*)r_csr);
5509
5510 //-----------------------------------------------------------------------
5511 // compute x + X v, accumulate offset to b (in w)
5512 //-----------------------------------------------------------------------
5513
5514 for ( i = 0; i < projectCurrSize_; i++ )
5515 {
5516 HYPRE_IJVectorGetObject(HYpbs_[i], (void **) &v_csr);
5517 HYPRE_ParVectorInnerProd( r_csr, v_csr, &alpha);
5518 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5519 (hypre_ParVector*)w_csr);
5520 HYPRE_IJVectorGetObject(HYpxs_[i], (void **) &v_csr);
5521 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5522 (hypre_ParVector*)x_csr);
5523 }
5524
5525 //-----------------------------------------------------------------------
5526 // save x and b away (and adjust b)
5527 //-----------------------------------------------------------------------
5528
5529 alpha = - 1.0;
5530 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)w_csr,(hypre_ParVector*)b_csr);
5531
5532 HYPRE_IJVectorGetObject(HYpxs_[projectSize_], (void **) &v_csr);
5533 HYPRE_ParVectorCopy( x_csr, v_csr );
5534 hypre_ParVectorScale(0.0,(hypre_ParVector*)x_csr);
5535
5536 //-----------------------------------------------------------------------
5537 // diagnostic message
5538 //-----------------------------------------------------------------------
5539
5540 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5541 printf("%4d : HYPRE_LSC:: leaving computeMinResProjection n", mypid_);
5542 return;
5543 }
5544
5545 //***************************************************************************
5546 // add a new pair of (x,b) vectors to the projection space
5547 //---------------------------------------------------------------------------
5548
addToMinResProjectionSpace(HYPRE_IJVector xvec,HYPRE_IJVector bvec)5549 void HYPRE_LinSysCore::addToMinResProjectionSpace(HYPRE_IJVector xvec,
5550 HYPRE_IJVector bvec)
5551 {
5552 int i, ierr, *partition, start_row, end_row;
5553 double alpha;
5554 HYPRE_ParVector v_csr, x_csr, xn_csr, b_csr, r_csr, bn_csr;
5555 HYPRE_ParCSRMatrix A_csr;
5556
5557 //-----------------------------------------------------------------------
5558 // diagnostic message
5559 //-----------------------------------------------------------------------
5560
5561 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5562 printf("%4d : HYPRE_LSC::addToProjectionSpace %d\n",mypid_,
5563 projectCurrSize_);
5564
5565 //-----------------------------------------------------------------------
5566 // fetch the matrix and vectors
5567 //-----------------------------------------------------------------------
5568
5569 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5570 HYPRE_IJVectorGetObject(xvec, (void **) &x_csr);
5571 HYPRE_IJVectorGetObject(bvec, (void **) &b_csr);
5572 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5573
5574 //-----------------------------------------------------------------------
5575 // initially, allocate space for B's and X's and R
5576 //-----------------------------------------------------------------------
5577
5578 if ( projectCurrSize_ == 0 && HYpbs_ == NULL )
5579 {
5580 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );
5581 start_row = partition[mypid_];
5582 end_row = partition[mypid_+1] - 1;
5583 free( partition );
5584 HYpxs_ = new HYPRE_IJVector[projectSize_+1];
5585 HYpbs_ = new HYPRE_IJVector[projectSize_+1];
5586
5587 for ( i = 0; i <= projectSize_; i++ )
5588 {
5589 ierr = HYPRE_IJVectorCreate(comm_, start_row, end_row, &(HYpbs_[i]));
5590 ierr = HYPRE_IJVectorSetObjectType(HYpbs_[i], HYPRE_PARCSR);
5591 ierr = HYPRE_IJVectorInitialize(HYpbs_[i]);
5592 ierr = HYPRE_IJVectorAssemble(HYpbs_[i]);
5593 hypre_assert( !ierr );
5594 }
5595 for ( i = 0; i <= projectSize_; i++ )
5596 {
5597 ierr = HYPRE_IJVectorCreate(comm_, start_row, end_row, &(HYpxs_[i]));
5598 ierr = HYPRE_IJVectorSetObjectType(HYpxs_[i], HYPRE_PARCSR);
5599 ierr = HYPRE_IJVectorInitialize(HYpxs_[i]);
5600 ierr = HYPRE_IJVectorAssemble(HYpxs_[i]);
5601 hypre_assert(!ierr);
5602 }
5603 }
5604
5605 //-----------------------------------------------------------------------
5606 // if buffer has been filled, move things up (but for now, restart)
5607 //-----------------------------------------------------------------------
5608
5609 if ( projectCurrSize_ >= projectSize_ )
5610 {
5611 //projectCurrSize_--;
5612 //tmpxvec = HYpxs_[0];
5613 //tmpbvec = HYpbs_[0];
5614 //for ( i = 0; i < projectCurrSize_; i++ )
5615 //{
5616 // HYpbs_[i] = HYpbs_[i+1];
5617 // HYpxs_[i] = HYpxs_[i+1];
5618 //}
5619 //HYpxs_[projectCurrSize_] = tmpxvec;
5620 //HYpbs_[projectCurrSize_] = tmpbvec;
5621 projectCurrSize_ = 0;
5622 }
5623
5624 //-----------------------------------------------------------------------
5625 // fetch projection vectors
5626 //-----------------------------------------------------------------------
5627
5628 HYPRE_IJVectorGetObject(HYpxs_[projectCurrSize_], (void **) &xn_csr);
5629 HYPRE_IJVectorGetObject(HYpbs_[projectCurrSize_], (void **) &bn_csr);
5630
5631 //-----------------------------------------------------------------------
5632 // copy incoming initial guess to buffer
5633 //-----------------------------------------------------------------------
5634
5635 HYPRE_ParVectorCopy( x_csr, xn_csr );
5636
5637 //-----------------------------------------------------------------------
5638 // compute bn = A * x
5639 //-----------------------------------------------------------------------
5640
5641 HYPRE_ParCSRMatrixMatvec( 1.0, A_csr, x_csr, 0.0, bn_csr );
5642 HYPRE_ParVectorCopy( bn_csr, r_csr );
5643
5644 //-----------------------------------------------------------------------
5645 // compute new vectors
5646 //-----------------------------------------------------------------------
5647
5648 for ( i = 0; i < projectCurrSize_; i++ )
5649 {
5650 HYPRE_IJVectorGetObject(HYpbs_[i], (void **) &v_csr);
5651 HYPRE_ParVectorInnerProd(r_csr, v_csr, &alpha);
5652 alpha = - alpha;
5653 if ( alpha != 0.0 )
5654 {
5655 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5656 (hypre_ParVector*)bn_csr);
5657 HYPRE_IJVectorGetObject(HYpxs_[i], (void **) &v_csr);
5658 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5659 (hypre_ParVector*)xn_csr);
5660 }
5661 }
5662 HYPRE_ParVectorInnerProd( bn_csr, bn_csr, &alpha);
5663 alpha = sqrt( alpha );
5664 if ( alpha != 0.0 )
5665 {
5666 alpha = 1.0 / alpha;
5667 hypre_ParVectorScale(alpha,(hypre_ParVector*)bn_csr);
5668 hypre_ParVectorScale(alpha,(hypre_ParVector*)xn_csr);
5669 projectCurrSize_++;
5670 }
5671
5672 //-----------------------------------------------------------------------
5673 // update final solution
5674 //-----------------------------------------------------------------------
5675
5676 if ( alpha != 0.0 )
5677 {
5678 HYPRE_IJVectorGetObject(HYpxs_[projectSize_], (void **) &v_csr);
5679 hypre_ParVectorAxpy(1.0,(hypre_ParVector*)v_csr,(hypre_ParVector*)x_csr);
5680
5681 HYPRE_IJVectorGetObject(HYpbs_[projectSize_], (void **) &v_csr);
5682 hypre_ParVectorAxpy(1.0,(hypre_ParVector*)v_csr,(hypre_ParVector*)b_csr);
5683 }
5684
5685 //-----------------------------------------------------------------------
5686 // diagnostic message
5687 //-----------------------------------------------------------------------
5688
5689 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5690 printf("%4d : HYPRE_LSC::leaving addToProjectionSpace %d\n",mypid_,
5691 projectCurrSize_);
5692 }
5693
5694 //***************************************************************************
5695 // project the initial guess into the previous solution space
5696 //
5697 // min || trans(x - xbar) A (x - xbar) ||
5698 //
5699 // solutions (phi_i) at previous steps
5700 //
5701 // (1) compute r = b - A * x_0
5702 // (2) compute alpha_i = (r, phi_i) for all previous stored vectors
5703 // (3) x_stored = x_0 + sum (alpha_i * phi_i)
5704 // (4) b_stored = A * x_0 + sum (alpha_i * psi_i)
5705 // (5) b = b - b_stored, x = 0
5706 //
5707 //---------------------------------------------------------------------------
5708
computeAConjProjection(HYPRE_ParCSRMatrix A_csr,HYPRE_ParVector x_csr,HYPRE_ParVector b_csr)5709 void HYPRE_LinSysCore::computeAConjProjection(HYPRE_ParCSRMatrix A_csr,
5710 HYPRE_ParVector x_csr, HYPRE_ParVector b_csr)
5711 {
5712 int i;
5713 double alpha;
5714 HYPRE_ParVector r_csr, v_csr, w_csr;
5715
5716 //-----------------------------------------------------------------------
5717 // diagnostic message
5718 //-----------------------------------------------------------------------
5719
5720 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5721 printf("%4d : HYPRE_LSC::entering computeAConjProjection %d\n",mypid_,
5722 projectCurrSize_);
5723 if ( projectCurrSize_ == 0 && HYpxs_ == NULL ) return;
5724
5725 //-----------------------------------------------------------------------
5726 // fetch vectors
5727 //-----------------------------------------------------------------------
5728
5729 HYPRE_IJVectorGetObject(currR_, (void **) &r_csr);
5730 HYPRE_IJVectorGetObject(HYpbs_[projectSize_], (void **) &w_csr);
5731
5732 //-----------------------------------------------------------------------
5733 // compute r = b - A x_0, save A * x_0 (to w)
5734 //-----------------------------------------------------------------------
5735
5736 HYPRE_ParCSRMatrixMatvec( 1.0, A_csr, x_csr, 0.0, w_csr );
5737 HYPRE_ParVectorCopy( b_csr, r_csr );
5738 alpha = -1.0;
5739 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)w_csr,(hypre_ParVector*)r_csr);
5740
5741 //-----------------------------------------------------------------------
5742 // compute alpha_i = (phi_i, r)
5743 // then x = x + alpha_i * phi_i for all i
5744 // then w = w + alpha_i * psi_i for all i
5745 //-----------------------------------------------------------------------
5746
5747 for ( i = 0; i < projectCurrSize_; i++ )
5748 {
5749 HYPRE_IJVectorGetObject(HYpxs_[i], (void **) &v_csr);
5750 HYPRE_ParVectorInnerProd(r_csr, v_csr, &alpha);
5751 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5752 (hypre_ParVector*)x_csr);
5753
5754 HYPRE_IJVectorGetObject(HYpbs_[i], (void **) &v_csr);
5755 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5756 (hypre_ParVector*)w_csr);
5757 }
5758
5759 //-----------------------------------------------------------------------
5760 // store x away
5761 //-----------------------------------------------------------------------
5762
5763 HYPRE_IJVectorGetObject(HYpxs_[projectSize_], (void **) &v_csr);
5764 HYPRE_ParVectorCopy( x_csr, v_csr );
5765 hypre_ParVectorScale(0.0,(hypre_ParVector*)x_csr);
5766
5767 //-----------------------------------------------------------------------
5768 // compute new residual b = b - w
5769 //-----------------------------------------------------------------------
5770
5771 alpha = -1.0;
5772 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)w_csr,(hypre_ParVector*)b_csr);
5773
5774 //-----------------------------------------------------------------------
5775 // diagnostic message
5776 //-----------------------------------------------------------------------
5777
5778 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5779 printf("%4d : HYPRE_LSC:: leaving computeAConjProjection n", mypid_);
5780 return;
5781 }
5782
5783 //***************************************************************************
5784 // add x to the projection space
5785 //
5786 // (1) compute alpha_i = (x, psi_i) for all previous stored vectors
5787 // (2) phi_n = x - sum(alpha_i * phi_i)
5788 // (3) phi_n = phi_n / norm(phi_n)_A
5789 // (4) psi_n = A * phi_n
5790 //---------------------------------------------------------------------------
5791
addToAConjProjectionSpace(HYPRE_IJVector xvec,HYPRE_IJVector bvec)5792 void HYPRE_LinSysCore::addToAConjProjectionSpace(HYPRE_IJVector xvec,
5793 HYPRE_IJVector bvec)
5794 {
5795 int i, ierr, *partition, start_row, end_row;
5796 double alpha;
5797 HYPRE_ParVector v_csr, x_csr, b_csr, bn_csr, xn_csr;
5798 //HYPRE_IJVector tmpxvec;
5799 HYPRE_ParCSRMatrix A_csr;
5800
5801 //-----------------------------------------------------------------------
5802 // diagnostic message
5803 //-----------------------------------------------------------------------
5804
5805 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5806 printf("%4d : HYPRE_LSC::addToAConjProjectionSpace %d\n",mypid_,
5807 projectCurrSize_);
5808
5809 //-----------------------------------------------------------------------
5810 // fetch the matrix and vectors
5811 //-----------------------------------------------------------------------
5812
5813 HYPRE_IJMatrixGetObject(currA_, (void **) &A_csr);
5814 HYPRE_IJVectorGetObject(xvec, (void **) &x_csr);
5815 HYPRE_IJVectorGetObject(bvec, (void **) &b_csr);
5816
5817 //-----------------------------------------------------------------------
5818 // initially, allocate space for the phi's and psi's
5819 //-----------------------------------------------------------------------
5820
5821 if ( projectCurrSize_ == 0 && HYpxs_ == NULL )
5822 {
5823 HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );
5824 start_row = partition[mypid_];
5825 end_row = partition[mypid_+1] - 1;
5826 free( partition );
5827 HYpxs_ = new HYPRE_IJVector[projectSize_+1];
5828 HYpbs_ = new HYPRE_IJVector[projectSize_+1];
5829
5830 for ( i = 0; i <= projectSize_; i++ )
5831 {
5832 ierr = HYPRE_IJVectorCreate(comm_, start_row, end_row, &(HYpbs_[i]));
5833 ierr = HYPRE_IJVectorSetObjectType(HYpbs_[i], HYPRE_PARCSR);
5834 ierr = HYPRE_IJVectorInitialize(HYpbs_[i]);
5835 ierr = HYPRE_IJVectorAssemble(HYpbs_[i]);
5836 hypre_assert( !ierr );
5837 }
5838 for ( i = 0; i <= projectSize_; i++ )
5839 {
5840 ierr = HYPRE_IJVectorCreate(comm_, start_row, end_row, &(HYpxs_[i]));
5841 ierr = HYPRE_IJVectorSetObjectType(HYpxs_[i], HYPRE_PARCSR);
5842 ierr = HYPRE_IJVectorInitialize(HYpxs_[i]);
5843 ierr = HYPRE_IJVectorAssemble(HYpxs_[i]);
5844 hypre_assert(!ierr);
5845 }
5846 }
5847
5848 //-----------------------------------------------------------------------
5849 // if buffer has been filled, move things up (but for now, restart)
5850 //-----------------------------------------------------------------------
5851
5852 if ( projectCurrSize_ >= projectSize_ )
5853 {
5854 //projectCurrSize_--;
5855 //tmpxvec = HYpxs_[0];
5856 //for ( i = 0; i < projectCurrSize_; i++ ) HYpxs_[i] = HYpxs_[i+1];
5857 //HYpxs_[projectCurrSize_] = tmpxvec;
5858 projectCurrSize_ = 0;
5859 }
5860
5861 //-----------------------------------------------------------------------
5862 // fetch the projection vectors
5863 //-----------------------------------------------------------------------
5864
5865 HYPRE_IJVectorGetObject(HYpxs_[projectCurrSize_], (void **) &xn_csr);
5866 HYPRE_IJVectorGetObject(HYpbs_[projectCurrSize_], (void **) &bn_csr);
5867
5868 //-----------------------------------------------------------------------
5869 // compute the new A-conjugate vector and its A-norm
5870 //-----------------------------------------------------------------------
5871
5872 HYPRE_ParVectorCopy( x_csr, xn_csr );
5873 for ( i = 0; i < projectCurrSize_; i++ )
5874 {
5875 HYPRE_IJVectorGetObject(HYpbs_[i], (void **) &v_csr);
5876 HYPRE_ParVectorInnerProd( x_csr, v_csr, &alpha);
5877 if ( alpha != 0.0 )
5878 {
5879 alpha = - alpha;
5880 HYPRE_IJVectorGetObject(HYpxs_[i], (void **) &v_csr);
5881 hypre_ParVectorAxpy(alpha,(hypre_ParVector*)v_csr,
5882 (hypre_ParVector*)xn_csr);
5883 }
5884 }
5885 HYPRE_ParCSRMatrixMatvec( 1.0, A_csr, xn_csr, 0.0, bn_csr );
5886 HYPRE_ParVectorInnerProd( xn_csr, bn_csr, &alpha);
5887 if ( alpha != 0.0 )
5888 {
5889 alpha = 1.0 / sqrt( alpha );
5890 hypre_ParVectorScale(alpha,(hypre_ParVector*)xn_csr);
5891 hypre_ParVectorScale(alpha,(hypre_ParVector*)bn_csr);
5892 projectCurrSize_++;
5893 }
5894
5895 //-----------------------------------------------------------------------
5896 // update final solution
5897 //-----------------------------------------------------------------------
5898
5899 if ( alpha != 0.0 )
5900 {
5901 HYPRE_IJVectorGetObject(HYpxs_[projectSize_], (void **) &v_csr);
5902 hypre_ParVectorAxpy(1.0,(hypre_ParVector*)v_csr,(hypre_ParVector*)x_csr);
5903
5904 HYPRE_IJVectorGetObject(HYpbs_[projectSize_], (void **) &v_csr);
5905 hypre_ParVectorAxpy(1.0,(hypre_ParVector*)v_csr,(hypre_ParVector*)b_csr);
5906 }
5907
5908 //-----------------------------------------------------------------------
5909 // diagnostic message
5910 //-----------------------------------------------------------------------
5911
5912 if ( (HYOutputLevel_ & HYFEI_SPECIALMASK) >= 4 )
5913 printf("%4d : HYPRE_LSC::leaving addToAConjProjectionSpace %d\n",mypid_,
5914 projectCurrSize_);
5915 }
5916
5917 //***************************************************************************
5918 //***************************************************************************
5919 // MLI specific functions
5920 //***************************************************************************
5921 //***************************************************************************
5922
5923 //***************************************************************************
5924 // initialize field information
5925 //---------------------------------------------------------------------------
5926
FE_initFields(int nFields,int * fieldSizes,int * fieldIDs)5927 void HYPRE_LinSysCore::FE_initFields(int nFields, int *fieldSizes,
5928 int *fieldIDs)
5929 {
5930 #ifdef HAVE_MLI
5931 if ( haveFEData_ == 1 && feData_ != NULL )
5932 HYPRE_LSI_MLIFEDataInitFields(feData_,nFields,fieldSizes,fieldIDs);
5933 #else
5934 (void) nFields;
5935 (void) fieldSizes;
5936 (void) fieldIDs;
5937 #endif
5938 return;
5939 }
5940
5941 //***************************************************************************
5942 // initialize element block
5943 //---------------------------------------------------------------------------
5944
FE_initElemBlock(int nElems,int nNodesPerElem,int numNodeFields,int * nodeFieldIDs)5945 void HYPRE_LinSysCore::FE_initElemBlock(int nElems, int nNodesPerElem,
5946 int numNodeFields, int *nodeFieldIDs)
5947 {
5948 #ifdef HAVE_MLI
5949 int status;
5950 if ( haveFEData_ == 1 && feData_ != NULL )
5951 {
5952 status = HYPRE_LSI_MLIFEDataInitElemBlock(feData_, nElems,
5953 nNodesPerElem, numNodeFields, nodeFieldIDs);
5954 if ( status )
5955 {
5956 if (haveFEData_ == 1) HYPRE_LSI_MLIFEDataDestroy(feData_);
5957 else if (haveFEData_ == 2) HYPRE_LSI_MLISFEIDestroy(feData_);
5958 feData_ = NULL;
5959 haveFEData_ = 0;
5960 }
5961 }
5962 #else
5963 (void) nElems;
5964 (void) nNodesPerElem;
5965 (void) numNodeFields;
5966 (void) nodeFieldIDs;
5967 #endif
5968 return;
5969 }
5970
5971 //***************************************************************************
5972 // initialize element node list
5973 //---------------------------------------------------------------------------
5974
FE_initElemNodeList(int elemID,int nNodesPerElem,int * nodeIDs)5975 void HYPRE_LinSysCore::FE_initElemNodeList(int elemID, int nNodesPerElem,
5976 int *nodeIDs)
5977 {
5978 #ifdef HAVE_MLI
5979 if ( haveFEData_ == 1 && feData_ != NULL )
5980 HYPRE_LSI_MLIFEDataInitElemNodeList(feData_, elemID, nNodesPerElem,
5981 nodeIDs);
5982 #else
5983 (void) elemID;
5984 (void) nNodesPerElem;
5985 (void) nodeIDs;
5986 #endif
5987 return;
5988 }
5989
5990 //***************************************************************************
5991 // initialize shared nodes
5992 //---------------------------------------------------------------------------
5993
FE_initSharedNodes(int nShared,int * sharedIDs,int * sharedPLengs,int ** sharedProcs)5994 void HYPRE_LinSysCore::FE_initSharedNodes(int nShared, int *sharedIDs,
5995 int *sharedPLengs, int **sharedProcs)
5996 {
5997 #ifdef HAVE_MLI
5998 if ( haveFEData_ == 1 && feData_ != NULL )
5999 HYPRE_LSI_MLIFEDataInitSharedNodes(feData_, nShared, sharedIDs,
6000 sharedPLengs, sharedProcs);
6001 #else
6002 (void) nShared;
6003 (void) sharedIDs;
6004 (void) sharedPLengs;
6005 (void) sharedProcs;
6006 #endif
6007 return;
6008 }
6009
6010 //***************************************************************************
6011 // initialize complete
6012 //---------------------------------------------------------------------------
6013
FE_initComplete()6014 void HYPRE_LinSysCore::FE_initComplete()
6015 {
6016 #ifdef HAVE_MLI
6017 if ( haveFEData_ == 1 && feData_ != NULL )
6018 HYPRE_LSI_MLIFEDataInitComplete(feData_);
6019 #endif
6020 return;
6021 }
6022
6023 //***************************************************************************
6024 // load element matrix
6025 //---------------------------------------------------------------------------
6026
FE_loadElemMatrix(int elemID,int nNodes,int * elemNodeList,int matDim,double ** elemMat)6027 void HYPRE_LinSysCore::FE_loadElemMatrix(int elemID, int nNodes,
6028 int *elemNodeList, int matDim, double **elemMat)
6029 {
6030 #ifdef HAVE_MLI
6031 if ( haveFEData_ == 1 && feData_ != NULL )
6032 HYPRE_LSI_MLIFEDataLoadElemMatrix(feData_, elemID, nNodes, elemNodeList,
6033 matDim, elemMat);
6034 #else
6035 (void) elemID;
6036 (void) nNodes;
6037 (void) elemNodeList;
6038 (void) matDim;
6039 (void) elemMat;
6040 #endif
6041 return;
6042 }
6043
6044 //***************************************************************************
6045 // build nodal coordinates
6046 // (to be used by AMS preconditioner, but it has been replaced with better
6047 // scheme. So when time is ripe, this function should be deleted.)
6048 //---------------------------------------------------------------------------
6049
HYPRE_LSI_BuildNodalCoordinates()6050 void HYPRE_LinSysCore::HYPRE_LSI_BuildNodalCoordinates()
6051 {
6052 int localNRows, *procNRows, *iTempArray, iP, iN, iS, *nodeProcMap;
6053 int *procList, nSends, *sendProcs, *sendLengs, **iSendBufs;
6054 int nRecvs, *recvLengs, *recvProcs, **iRecvBufs, iR, numNodes, eqnInd;
6055 int *flags, arrayLeng, coordLength, ierr, procIndex, iD;
6056 double **dSendBufs, **dRecvBufs, *nCoords, *vecData;
6057 MPI_Request *mpiRequests;
6058 MPI_Status mpiStatus;
6059 HYPRE_ParVector parVec;
6060
6061 /* -------------------------------------------------------- */
6062 /* construct procNRows array */
6063 /* -------------------------------------------------------- */
6064
6065 localNRows = localEndRow_ - localStartRow_ + 1;
6066 procNRows = new int[numProcs_+1];
6067 iTempArray = new int[numProcs_];
6068 for (iP = 0; iP <= numProcs_; iP++) procNRows[iP] = 0;
6069 procNRows[mypid_] = localNRows;
6070 MPI_Allreduce(procNRows,iTempArray,numProcs_,MPI_INT,MPI_SUM,comm_);
6071 procNRows[0] = 0;
6072 for (iP = 1; iP <= numProcs_; iP++)
6073 procNRows[iP] = procNRows[iP-1] + iTempArray[iP-1];
6074
6075 /* -------------------------------------------------------- */
6076 /* construct node to processor map */
6077 /* -------------------------------------------------------- */
6078
6079 nodeProcMap = new int[MLI_NumNodes_];
6080 for (iN = 0; iN < MLI_NumNodes_; iN++)
6081 {
6082 nodeProcMap[iN] = -1;
6083 if (MLI_EqnNumbers_[iN] < procNRows[mypid_] ||
6084 MLI_EqnNumbers_[iN] >= procNRows[mypid_+1])
6085 {
6086 for (iP = 0; iP < numProcs_; iP++)
6087 if (MLI_EqnNumbers_[iN] < procNRows[iP]) break;
6088 nodeProcMap[iN] = iP - 1;
6089 }
6090 }
6091
6092 /* -------------------------------------------------------- */
6093 /* construct send information */
6094 /* -------------------------------------------------------- */
6095
6096 procList = new int[numProcs_];
6097 for (iP = 0; iP < numProcs_; iP++) procList[iP] = 0;
6098 for (iN = 0; iN < numProcs_; iN++)
6099 if (nodeProcMap[iN] >= 0) procList[nodeProcMap[iN]]++;
6100 nSends = 0;
6101 for (iP = 0; iP < numProcs_; iP++) if (procList[iP] > 0) nSends++;
6102 if (nSends > 0)
6103 {
6104 sendProcs = new int[nSends];
6105 sendLengs = new int[nSends];
6106 iSendBufs = new int*[nSends];
6107 dSendBufs = new double*[nSends];
6108 }
6109 nSends = 0;
6110 for (iP = 0; iP < numProcs_; iP++)
6111 {
6112 if (procList[iP] > 0)
6113 {
6114 sendLengs[nSends] = procList[iP];
6115 sendProcs[nSends++] = iP;
6116 }
6117 }
6118
6119 /* -------------------------------------------------------- */
6120 /* construct recv information */
6121 /* -------------------------------------------------------- */
6122
6123 for (iP = 0; iP < numProcs_; iP++) procList[iP] = 0;
6124 for (iP = 0; iP < nSends; iP++) procList[sendProcs[iP]]++;
6125 MPI_Allreduce(procList,iTempArray,numProcs_,MPI_INT,MPI_SUM,comm_);
6126 nRecvs = iTempArray[mypid_];
6127 if (nRecvs > 0)
6128 {
6129 recvLengs = new int[nRecvs];
6130 recvProcs = new int[nRecvs];
6131 iRecvBufs = new int*[nRecvs];
6132 dRecvBufs = new double*[nRecvs];
6133 mpiRequests = new MPI_Request[nRecvs];
6134 }
6135 for (iP = 0; iP < nRecvs; iP++)
6136 MPI_Irecv(&(recvLengs[iP]), 1, MPI_INT, MPI_ANY_SOURCE, 29421,
6137 comm_, &(mpiRequests[iP]));
6138 for (iP = 0; iP < nSends; iP++)
6139 MPI_Send(&(sendLengs[iP]), 1, MPI_INT, sendProcs[iP], 29421, comm_);
6140 for (iP = 0; iP < nRecvs; iP++)
6141 {
6142 MPI_Wait(&(mpiRequests[iP]), &mpiStatus);
6143 recvProcs[iP] = mpiStatus.MPI_SOURCE;
6144 }
6145
6146 /* -------------------------------------------------------- */
6147 /* communicate equation numbers information */
6148 /* -------------------------------------------------------- */
6149
6150 for (iP = 0; iP < nRecvs; iP++)
6151 {
6152 iRecvBufs[iP] = new int[recvLengs[iP]];
6153 MPI_Irecv(iRecvBufs[iP], recvLengs[iP], MPI_INT, recvProcs[iP],
6154 29422, comm_, &(mpiRequests[iP]));
6155 }
6156 for (iP = 0; iP < nSends; iP++)
6157 {
6158 iSendBufs[iP] = new int[sendLengs[iP]];
6159 sendLengs[iP] = 0;
6160 }
6161 for (iN = 0; iN < MLI_NumNodes_; iN++)
6162 {
6163 if (nodeProcMap[iN] >= 0)
6164 {
6165 procIndex = nodeProcMap[iN];
6166 for (iP = 0; iP < nSends; iP++)
6167 if (procIndex == sendProcs[iP]) break;
6168 iSendBufs[iP][sendLengs[iP]++] = MLI_EqnNumbers_[iN];
6169 }
6170 }
6171 for (iP = 0; iP < nSends; iP++)
6172 {
6173 MPI_Send(iSendBufs[iP], sendLengs[iP], MPI_INT, sendProcs[iP],
6174 29422, comm_);
6175 }
6176 for (iP = 0; iP < nRecvs; iP++) MPI_Wait(&(mpiRequests[iP]),&mpiStatus);
6177
6178 /* -------------------------------------------------------- */
6179 /* communicate coordinate information */
6180 /* -------------------------------------------------------- */
6181
6182 for (iP = 0; iP < nRecvs; iP++)
6183 {
6184 dRecvBufs[iP] = new double[recvLengs[iP]*MLI_FieldSize_];
6185 MPI_Irecv(dRecvBufs[iP], recvLengs[iP]*MLI_FieldSize_, MPI_DOUBLE,
6186 recvProcs[iP], 29425, comm_, &(mpiRequests[iP]));
6187 }
6188 for (iP = 0; iP < nSends; iP++)
6189 {
6190 dSendBufs[iP] = new double[sendLengs[iP]*MLI_FieldSize_];
6191 sendLengs[iP] = 0;
6192 }
6193 for (iN = 0; iN < MLI_NumNodes_; iN++)
6194 {
6195 if (nodeProcMap[iN] >= 0)
6196 {
6197 procIndex = nodeProcMap[iN];
6198 for (iP = 0; iP < nSends; iP++)
6199 if (procIndex == sendProcs[iP]) break;
6200 for (iD = 0; iD < MLI_FieldSize_; iD++)
6201 dSendBufs[iP][sendLengs[iP]++] =
6202 MLI_NodalCoord_[iN*MLI_FieldSize_+iD];
6203 }
6204 }
6205 for (iP = 0; iP < nSends; iP++)
6206 {
6207 sendLengs[iP] /= MLI_FieldSize_;
6208 MPI_Send(dSendBufs[iP], sendLengs[iP]*MLI_FieldSize_, MPI_DOUBLE,
6209 sendProcs[iP], 29425, comm_);
6210 }
6211 for (iP = 0; iP < nRecvs; iP++) MPI_Wait(&(mpiRequests[iP]),&mpiStatus);
6212
6213 /* -------------------------------------------------------- */
6214 /* check any duplicate coordinate information */
6215 /* -------------------------------------------------------- */
6216
6217 arrayLeng = MLI_NumNodes_;
6218 for (iP = 0; iP < nRecvs; iP++) arrayLeng += recvLengs[iP];
6219 flags = new int[arrayLeng];
6220 for (iN = 0; iN < arrayLeng; iN++) flags[iN] = 0;
6221 for (iN = 0; iN < MLI_NumNodes_; iN++)
6222 {
6223 if (nodeProcMap[iN] < 0)
6224 {
6225 eqnInd = (MLI_EqnNumbers_[iN] - procNRows[mypid_]) / MLI_FieldSize_;
6226 if (eqnInd >= arrayLeng)
6227 {
6228 printf("%d : LoadNodalCoordinates - ERROR(1).\n", mypid_);
6229 exit(1);
6230 }
6231 flags[eqnInd] = 1;
6232 }
6233 }
6234 for (iP = 0; iP < nRecvs; iP++)
6235 {
6236 for (iR = 0; iR < recvLengs[iP]; iR++)
6237 {
6238 eqnInd = (iRecvBufs[iP][iR] - procNRows[mypid_]) / MLI_FieldSize_;
6239 if (eqnInd >= arrayLeng)
6240 {
6241 printf("%d : LoadNodalCoordinates - ERROR(2).\n", mypid_);
6242 exit(1);
6243 }
6244 flags[eqnInd] = 1;
6245 }
6246 }
6247 numNodes = 0;
6248 for (iN = 0; iN < arrayLeng; iN++)
6249 {
6250 if ( flags[iN] == 0 ) break;
6251 else numNodes++;
6252 }
6253 delete [] flags;
6254
6255 /* -------------------------------------------------------- */
6256 /* set up nodal coordinate information in correct order */
6257 /* -------------------------------------------------------- */
6258
6259 coordLength = MLI_NumNodes_ * MLI_FieldSize_;
6260 nCoords = new double[coordLength];
6261
6262 arrayLeng = MLI_NumNodes_ * MLI_FieldSize_;
6263 for (iN = 0; iN < MLI_NumNodes_; iN++)
6264 {
6265 if (nodeProcMap[iN] < 0)
6266 {
6267 eqnInd = (MLI_EqnNumbers_[iN] - procNRows[mypid_]) / MLI_FieldSize_;
6268 if (eqnInd >= 0 && eqnInd < arrayLeng)
6269 for (iD = 0; iD < MLI_FieldSize_; iD++)
6270 nCoords[eqnInd*MLI_FieldSize_+iD] =
6271 MLI_NodalCoord_[iN*MLI_FieldSize_+iD];
6272 }
6273 }
6274 for (iP = 0; iP < nRecvs; iP++)
6275 {
6276 for (iR = 0; iR < recvLengs[iP]; iR++)
6277 {
6278 eqnInd = (iRecvBufs[iP][iR] - procNRows[mypid_]) / MLI_FieldSize_;
6279 if (eqnInd >= 0 && eqnInd < arrayLeng)
6280 for (iD = 0; iD < MLI_FieldSize_; iD++)
6281 nCoords[eqnInd*MLI_FieldSize_+iD] =
6282 dRecvBufs[iP][iR*MLI_FieldSize_+iD];
6283 }
6284 }
6285
6286 /* -------------------------------------------------------- */
6287 /* create AMS vectors */
6288 /* -------------------------------------------------------- */
6289
6290 localNRows = localEndRow_ - localStartRow_ + 1;
6291 ierr = HYPRE_IJVectorCreate(comm_,(localStartRow_-1)/MLI_FieldSize_,
6292 localEndRow_/MLI_FieldSize_-1, &amsX_);
6293 ierr += HYPRE_IJVectorSetObjectType(amsX_, HYPRE_PARCSR);
6294 ierr += HYPRE_IJVectorInitialize(amsX_);
6295 ierr += HYPRE_IJVectorAssemble(amsX_);
6296 hypre_assert(!ierr);
6297 HYPRE_IJVectorGetObject(amsX_, (void **) &parVec);
6298 vecData = (double *) hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *) parVec));
6299 for (iN = 0; iN < localNRows/MLI_FieldSize_; iN++)
6300 vecData[iN] = nCoords[iN*MLI_FieldSize_];
6301 ierr = HYPRE_IJVectorCreate(comm_,(localStartRow_-1)/MLI_FieldSize_,
6302 localEndRow_/MLI_FieldSize_-1, &amsY_);
6303 ierr += HYPRE_IJVectorSetObjectType(amsY_, HYPRE_PARCSR);
6304 ierr += HYPRE_IJVectorInitialize(amsY_);
6305 ierr += HYPRE_IJVectorAssemble(amsY_);
6306 hypre_assert(!ierr);
6307 HYPRE_IJVectorGetObject(amsY_, (void **) &parVec);
6308 vecData = (double *) hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *) parVec));
6309 for (iN = 0; iN < localNRows/MLI_FieldSize_; iN++)
6310 vecData[iN] = nCoords[iN*MLI_FieldSize_+1];
6311 ierr = HYPRE_IJVectorCreate(comm_,(localStartRow_-1)/MLI_FieldSize_,
6312 localEndRow_/MLI_FieldSize_-1, &amsZ_);
6313 ierr += HYPRE_IJVectorSetObjectType(amsZ_, HYPRE_PARCSR);
6314 ierr += HYPRE_IJVectorInitialize(amsZ_);
6315 ierr += HYPRE_IJVectorAssemble(amsZ_);
6316 hypre_assert(!ierr);
6317 HYPRE_IJVectorGetObject(amsZ_, (void **) &parVec);
6318 vecData = (double *) hypre_VectorData(hypre_ParVectorLocalVector((hypre_ParVector *) parVec));
6319 for (iN = 0; iN < localNRows/MLI_FieldSize_; iN++)
6320 vecData[iN] = nCoords[iN*MLI_FieldSize_+2];
6321
6322 /* -------------------------------------------------------- */
6323 /* clean up */
6324 /* -------------------------------------------------------- */
6325
6326 delete [] procList;
6327 delete [] iTempArray;
6328 delete [] nodeProcMap;
6329 delete [] procNRows;
6330 delete [] nCoords;
6331 if (nSends > 0)
6332 {
6333 delete [] sendProcs;
6334 delete [] sendLengs;
6335 for (iS = 0; iS < nSends; iS++) delete [] iSendBufs[iS];
6336 for (iS = 0; iS < nSends; iS++) delete [] dSendBufs[iS];
6337 delete [] dSendBufs;
6338 delete [] iSendBufs;
6339 }
6340 if (nRecvs > 0)
6341 {
6342 delete [] recvProcs;
6343 delete [] recvLengs;
6344 for (iR = 0; iR < nRecvs; iR++) delete [] iRecvBufs[iR];
6345 for (iR = 0; iR < nRecvs; iR++) delete [] dRecvBufs[iR];
6346 delete [] iRecvBufs;
6347 delete [] dRecvBufs;
6348 delete [] mpiRequests;
6349 }
6350 return;
6351 }
6352
6353