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 // Date : Apr 26, 2002
10 //***************************************************************************
11 // system includes
12 //---------------------------------------------------------------------------
13
14 #include <stdlib.h>
15 #include <string.h>
16 #include <stdio.h>
17
18 #if 0 /* RDF: Not sure this is really needed */
19 #ifdef WIN32
20 #define strcmp _stricmp
21 #endif
22 #endif
23
24 //***************************************************************************
25 // local includes
26 //---------------------------------------------------------------------------
27
28 #include "HYPRE.h"
29 #include "HYPRE_LSI_UZAWA.h"
30
31 //---------------------------------------------------------------------------
32 // MLI include files
33 //---------------------------------------------------------------------------
34
35 #ifdef HAVE_MLI
36 #include "HYPRE_LSI_mli.h"
37 #endif
38
39 //***************************************************************************
40 // local defines and external functions
41 //---------------------------------------------------------------------------
42
43 extern "C"
44 {
45 int hypre_BoomerAMGBuildCoarseOperator(hypre_ParCSRMatrix*,
46 hypre_ParCSRMatrix*,
47 hypre_ParCSRMatrix*,
48 hypre_ParCSRMatrix**);
49 }
50
51 //***************************************************************************
52 //***************************************************************************
53 // C-Interface data structure
54 //---------------------------------------------------------------------------
55
56 typedef struct HYPRE_LSI_Uzawa_Struct
57 {
58 void *precon;
59 }
60 HYPRE_LSI_UzawaStruct;
61
62 //***************************************************************************
63 //***************************************************************************
64 // C-Interface functions to solver
65 //---------------------------------------------------------------------------
66
HYPRE_LSI_UzawaCreate(MPI_Comm mpi_comm,HYPRE_Solver * solver)67 extern "C" int HYPRE_LSI_UzawaCreate(MPI_Comm mpi_comm, HYPRE_Solver *solver)
68 {
69 (void) mpi_comm;
70 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *)
71 hypre_CTAlloc(HYPRE_LSI_UzawaStruct, 1, HYPRE_MEMORY_HOST);
72 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) new HYPRE_LSI_Uzawa(mpi_comm);
73 cprecon->precon = (void *) precon;
74 (*solver) = (HYPRE_Solver) cprecon;
75 return 0;
76 }
77
78 //***************************************************************************
79
HYPRE_LSI_UzawaDestroy(HYPRE_Solver solver)80 extern "C" int HYPRE_LSI_UzawaDestroy(HYPRE_Solver solver)
81 {
82 int err=0;
83
84 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
85 if ( cprecon == NULL ) err = 1;
86 else
87 {
88 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
89 if ( precon != NULL ) delete precon;
90 else err = 1;
91 free( cprecon );
92 }
93 return err;
94 }
95
96 //***************************************************************************
97
HYPRE_LSI_UzawaSetParams(HYPRE_Solver solver,char * params)98 extern "C" int HYPRE_LSI_UzawaSetParams(HYPRE_Solver solver, char *params)
99 {
100 int err=0;
101
102 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
103 if ( cprecon == NULL ) err = 1;
104 else
105 {
106 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
107 err = precon->setParams(params);
108 }
109 return err;
110 }
111
112 //***************************************************************************
113
HYPRE_LSI_UzawaSetMaxIterations(HYPRE_Solver solver,int iter)114 extern "C" int HYPRE_LSI_UzawaSetMaxIterations(HYPRE_Solver solver, int iter)
115 {
116 int err=0;
117
118 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
119 if ( cprecon == NULL ) err = 1;
120 else
121 {
122 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
123 err = precon->setMaxIterations(iter);
124 }
125 return err;
126 }
127
128 //***************************************************************************
129
HYPRE_LSI_UzawaSetTolerance(HYPRE_Solver solver,double tol)130 extern "C" int HYPRE_LSI_UzawaSetTolerance(HYPRE_Solver solver, double tol)
131 {
132 int err=0;
133
134 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
135 if ( cprecon == NULL ) err = 1;
136 else
137 {
138 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
139 err = precon->setTolerance(tol);
140 }
141 return err;
142 }
143
144 //***************************************************************************
145
HYPRE_LSI_UzawaGetNumIterations(HYPRE_Solver solver,int * iter)146 extern "C" int HYPRE_LSI_UzawaGetNumIterations(HYPRE_Solver solver, int *iter)
147 {
148 int err=0;
149
150 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
151 if ( cprecon == NULL ) err = 1;
152 else
153 {
154 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
155 err = precon->getNumIterations(*iter);
156 }
157 return err;
158 }
159
160 //***************************************************************************
161
162 extern "C"
HYPRE_LSI_UzawaSetup(HYPRE_Solver solver,HYPRE_ParCSRMatrix Amat,HYPRE_ParVector b,HYPRE_ParVector x)163 int HYPRE_LSI_UzawaSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix Amat,
164 HYPRE_ParVector b, HYPRE_ParVector x)
165 {
166 int err=0;
167
168 (void) b;
169 (void) x;
170 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
171 if ( cprecon == NULL ) err = 1;
172 else
173 {
174 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
175 err = precon->setup(Amat, x, b);
176 }
177 return err;
178 }
179
180 //***************************************************************************
181
182 extern "C"
HYPRE_LSI_UzawaSolve(HYPRE_Solver solver,HYPRE_ParCSRMatrix Amat,HYPRE_ParVector b,HYPRE_ParVector x)183 int HYPRE_LSI_UzawaSolve(HYPRE_Solver solver, HYPRE_ParCSRMatrix Amat,
184 HYPRE_ParVector b, HYPRE_ParVector x)
185 {
186 int err=0;
187
188 (void) Amat;
189 HYPRE_LSI_UzawaStruct *cprecon = (HYPRE_LSI_UzawaStruct *) solver;
190 if ( cprecon == NULL ) err = 1;
191 else
192 {
193 HYPRE_LSI_Uzawa *precon = (HYPRE_LSI_Uzawa *) cprecon->precon;
194 err = precon->solve(b, x);
195 }
196 return err;
197 }
198
199 //***************************************************************************
200 //***************************************************************************
201 // Constructor
202 //---------------------------------------------------------------------------
203
HYPRE_LSI_Uzawa(MPI_Comm comm)204 HYPRE_LSI_Uzawa::HYPRE_LSI_Uzawa(MPI_Comm comm)
205 {
206 Amat_ = NULL;
207 A11mat_ = NULL;
208 A12mat_ = NULL;
209 S22mat_ = NULL;
210 mpiComm_ = comm;
211 outputLevel_ = 2;
212 S22Scheme_ = 0;
213 maxIterations_ = 1;
214 tolerance_ = 1.0e-6;
215 numIterations_ = 0;
216 procA22Sizes_ = NULL;
217 A11Solver_ = NULL;
218 A11Precond_ = NULL;
219 S22Solver_ = NULL;
220 S22Precond_ = NULL;
221 modifiedScheme_ = 0;
222 S22SolverDampFactor_ = 1.0;
223 A11Params_.SolverID_ = 1; /* default : gmres */
224 S22Params_.SolverID_ = 1; /* default : cg */
225 A11Params_.PrecondID_ = 1; /* default : diagonal */
226 S22Params_.PrecondID_ = 1; /* default : diagonal */
227 A11Params_.Tol_ = 1.0e-3;
228 S22Params_.Tol_ = 1.0e-3;
229 A11Params_.MaxIter_ = 1000;
230 S22Params_.MaxIter_ = 1000;
231 A11Params_.PSNLevels_ = 1;
232 S22Params_.PSNLevels_ = 1;
233 A11Params_.PSThresh_ = 1.0e-1;
234 S22Params_.PSThresh_ = 1.0e-1;
235 A11Params_.PSFilter_ = 2.0e-1;
236 S22Params_.PSFilter_ = 2.0e-1;
237 A11Params_.AMGThresh_ = 7.5e-1;
238 S22Params_.AMGThresh_ = 7.5e-1;
239 A11Params_.AMGNSweeps_ = 2;
240 S22Params_.AMGNSweeps_ = 2;
241 A11Params_.AMGSystemSize_ = 1;
242 S22Params_.AMGSystemSize_ = 1;
243 A11Params_.PilutFillin_ = 100;
244 S22Params_.PilutFillin_ = 100;
245 A11Params_.PilutDropTol_ = 0.1;
246 S22Params_.PilutDropTol_ = 0.1;
247 A11Params_.EuclidNLevels_ = 1;
248 S22Params_.EuclidNLevels_ = 1;
249 A11Params_.EuclidThresh_ = 0.1;
250 S22Params_.EuclidThresh_ = 0.1;
251 A11Params_.MLIThresh_ = 0.08;
252 S22Params_.MLIThresh_ = 0.08;
253 A11Params_.MLINSweeps_ = 2;
254 S22Params_.MLINSweeps_ = 2;
255 A11Params_.MLIPweight_ = 0.0;
256 S22Params_.MLIPweight_ = 0.0;
257 A11Params_.MLINodeDOF_ = 3;
258 S22Params_.MLINodeDOF_ = 3;
259 A11Params_.MLINullDim_ = 3;
260 S22Params_.MLINullDim_ = 3;
261 }
262
263 //***************************************************************************
264 // destructor
265 //---------------------------------------------------------------------------
266
~HYPRE_LSI_Uzawa()267 HYPRE_LSI_Uzawa::~HYPRE_LSI_Uzawa()
268 {
269 Amat_ = NULL;
270 mpiComm_ = 0;
271 if ( procA22Sizes_ != NULL ) delete [] procA22Sizes_;
272 if ( A11mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(A11mat_);
273 if ( A12mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(A12mat_);
274 if ( S22mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(S22mat_);
275 }
276
277 //***************************************************************************
278 // set internal parameters
279 //---------------------------------------------------------------------------
280
setParams(char * params)281 int HYPRE_LSI_Uzawa::setParams(char *params)
282 {
283 char param1[256], param2[256], param3[256];
284
285 sscanf(params,"%s", param1);
286 if ( strcmp(param1, "Uzawa") )
287 {
288 printf("HYPRE_LSI_Uzawa::parameters not for me.\n");
289 return 1;
290 }
291 sscanf(params,"%s %s", param1, param2);
292 if ( !strcmp(param2, "help") )
293 {
294 printf("Available options for Uzawa are : \n");
295 printf(" outputLevel <d> \n");
296 printf(" A11Solver <cg,gmres> \n");
297 printf(" A11Tolerance <f> \n");
298 printf(" A11MaxIterations <d> \n");
299 printf(" A11Precon <pilut,boomeramg,euclid,parasails,ddilut,mli>\n");
300 printf(" A11PreconPSNlevels <d> \n");
301 printf(" A11PreconPSThresh <f> \n");
302 printf(" A11PreconPSFilter <f> \n");
303 printf(" A11PreconAMGThresh <f> \n");
304 printf(" A11PreconAMGNumSweeps <d> \n");
305 printf(" A11PreconAMGSystemSize <d> \n");
306 printf(" A11PreconEuclidNLevels <d> \n");
307 printf(" A11PreconEuclidThresh <f> \n");
308 printf(" A11PreconPilutFillin <d> \n");
309 printf(" A11PreconPilutDropTol <f> \n");
310 printf(" S22SolverDampingFactor <f> \n");
311 printf(" S22Solver <cg,gmres> \n");
312 printf(" S22Tolerance <f> \n");
313 printf(" S22MaxIterations <d> \n");
314 printf(" S22Precon <pilut,boomeramg,euclid,parasails,ddilut,mli>\n");
315 printf(" S22PreconPSNlevels <d> \n");
316 printf(" S22PreconPSThresh <f> \n");
317 printf(" S22PreconPSFilter <f> \n");
318 printf(" S22PreconAMGThresh <f> \n");
319 printf(" S22PreconAMGNumSweeps <d> \n");
320 printf(" S22PreconAMGSystemSize <d> \n");
321 printf(" S22PreconEuclidNLevels <d> \n");
322 printf(" S22PreconEuclidThresh <f> \n");
323 printf(" S22PreconPilutFillin <d> \n");
324 printf(" S22PreconPilutDropTol <f> \n");
325 }
326 else if ( !strcmp(param2, "outputLevel") )
327 {
328 sscanf(params,"%s %s %d", param1, param2, &outputLevel_);
329 if ( outputLevel_ > 0 )
330 printf("HYPRE_LSI_Uzawa::outputLevel = %d.\n", outputLevel_);
331 }
332 else if ( !strcmp(param2, "modified") )
333 {
334 modifiedScheme_ = 1;
335 if ( outputLevel_ > 0 ) printf("HYPRE_LSI_Uzawa::3 level scheme.\n");
336 }
337 else if ( !strcmp(param2, "A11Solver") )
338 {
339 sscanf(params,"%s %s %s", param1, param2, param3);
340 if ( !strcmp(param3, "none") )
341 {
342 A11Params_.SolverID_ = 0;
343 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11 solver = cg\n");
344 }
345 else if ( !strcmp(param3, "cg") )
346 {
347 A11Params_.SolverID_ = 1;
348 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11 solver = cg\n");
349 }
350 else if ( !strcmp(param3, "gmres") )
351 {
352 A11Params_.SolverID_ = 2;
353 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11 solver = gmres\n");
354 }
355 }
356 else if ( !strcmp(param2, "S22Solver") )
357 {
358 sscanf(params,"%s %s %s", param1, param2, param3);
359 if ( !strcmp(param3, "none") )
360 {
361 S22Params_.SolverID_ = 0;
362 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22 solver = cg\n");
363 }
364 else if ( !strcmp(param3, "cg") )
365 {
366 S22Params_.SolverID_ = 1;
367 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22 solver = cg\n");
368 }
369 else if ( !strcmp(param3, "gmres") )
370 {
371 S22Params_.SolverID_ = 2;
372 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22 solver = gmres\n");
373 }
374 }
375 else if ( !strcmp(param2, "S22SolverDampingFactor") )
376 {
377 sscanf(params,"%s %s %lg", param1, param2, &S22SolverDampFactor_);
378 if ( S22SolverDampFactor_ < 0.0 ) S22SolverDampFactor_ = 1.0;
379 }
380 else if ( !strcmp(param2, "A11Tolerance") )
381 {
382 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.Tol_));
383 if ( A11Params_.Tol_ >= 1.0 || A11Params_.Tol_ <= 0.0 )
384 A11Params_.Tol_ = 1.0e-12;
385 if (outputLevel_ > 0)
386 printf("HYPRE_LSI_Uzawa::A11 tol = %e\n", A11Params_.Tol_);
387 }
388 else if ( !strcmp(param2, "S22Tolerance") )
389 {
390 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.Tol_));
391 if ( S22Params_.Tol_ >= 1.0 || S22Params_.Tol_ <= 0.0 )
392 S22Params_.Tol_ = 1.0e-12;
393 if (outputLevel_ > 0)
394 printf("HYPRE_LSI_Uzawa::S22 tol = %e\n", S22Params_.Tol_);
395 }
396 else if ( !strcmp(param2, "A11MaxIterations") )
397 {
398 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.MaxIter_));
399 if ( A11Params_.MaxIter_ <= 0 ) A11Params_.MaxIter_ = 10;
400 if (outputLevel_ > 0)
401 printf("HYPRE_LSI_Uzawa::A11 maxiter = %d\n", A11Params_.MaxIter_);
402 }
403 else if ( !strcmp(param2, "S22MaxIterations") )
404 {
405 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.MaxIter_));
406 if ( S22Params_.MaxIter_ <= 0 ) S22Params_.MaxIter_ = 10;
407 if (outputLevel_ > 0)
408 printf("HYPRE_LSI_Uzawa::S22 maxiter = %d\n", S22Params_.MaxIter_);
409 }
410 else if ( !strcmp(param2, "A11Precon") )
411 {
412 sscanf(params,"%s %s %s", param1, param2, param3);
413 if ( !strcmp(param3, "diagonal") )
414 {
415 A11Params_.PrecondID_ = 1;
416 if (outputLevel_ > 0)
417 printf("HYPRE_LSI_Uzawa::A11 precon = diagonal\n");
418 }
419 else if ( !strcmp(param3, "parasails") )
420 {
421 A11Params_.PrecondID_ = 2;
422 if (outputLevel_ > 0)
423 printf("HYPRE_LSI_Uzawa::A11 precon = parasails\n");
424 }
425 else if ( !strcmp(param3, "boomeramg") )
426 {
427 A11Params_.PrecondID_ = 3;
428 if (outputLevel_ > 0)
429 printf("HYPRE_LSI_Uzawa::A11 precon = boomeramg\n");
430 }
431 else if ( !strcmp(param3, "pilut") )
432 {
433 A11Params_.PrecondID_ = 4;
434 if (outputLevel_ > 0)
435 printf("HYPRE_LSI_Uzawa::A11 precon = pilut\n");
436 }
437 else if ( !strcmp(param3, "euclid") )
438 {
439 A11Params_.PrecondID_ = 5;
440 if (outputLevel_ > 0)
441 printf("HYPRE_LSI_Uzawa::A11 precon = euclid\n");
442 }
443 else if ( !strcmp(param3, "mli") )
444 {
445 A11Params_.PrecondID_ = 6;
446 if (outputLevel_ > 0)
447 printf("HYPRE_LSI_Uzawa::A11 precon = MLISA\n");
448 }
449 }
450 else if ( !strcmp(param2, "S22Precon") )
451 {
452 sscanf(params,"%s %s %s", param1, param2, param3);
453 if ( !strcmp(param3, "diagonal") )
454 {
455 S22Params_.PrecondID_ = 1;
456 if (outputLevel_ > 0)
457 printf("HYPRE_LSI_Uzawa::S22 precon = diagonal\n");
458 }
459 else if ( !strcmp(param3, "parasails") )
460 {
461 S22Params_.PrecondID_ = 2;
462 if (outputLevel_ > 0)
463 printf("HYPRE_LSI_Uzawa::S22 precon = parasails\n");
464 }
465 else if ( !strcmp(param3, "boomeramg") )
466 {
467 S22Params_.PrecondID_ = 3;
468 if (outputLevel_ > 0)
469 printf("HYPRE_LSI_Uzawa::S22 precon = boomeramg\n");
470 }
471 else if ( !strcmp(param3, "pilut") )
472 {
473 S22Params_.PrecondID_ = 4;
474 if (outputLevel_ > 0)
475 printf("HYPRE_LSI_Uzawa::S22 precon = pilut\n");
476 }
477 else if ( !strcmp(param3, "euclid") )
478 {
479 S22Params_.PrecondID_ = 5;
480 if (outputLevel_ > 0)
481 printf("HYPRE_LSI_Uzawa::S22 precon = euclid\n");
482 }
483 else if ( !strcmp(param3, "mli") )
484 {
485 S22Params_.PrecondID_ = 6;
486 if (outputLevel_ > 0)
487 printf("HYPRE_LSI_Uzawa::S22 precon = MLISA\n");
488 }
489 }
490 else if ( !strcmp(param2, "A11PreconPSNlevels") )
491 {
492 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.PSNLevels_));
493 if ( A11Params_.PSNLevels_ < 0 ) A11Params_.PSNLevels_ = 0;
494 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconPSNLevels\n");
495 }
496 else if ( !strcmp(param2, "S22PreconPSNlevels") )
497 {
498 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.PSNLevels_));
499 if ( S22Params_.PSNLevels_ < 0 ) S22Params_.PSNLevels_ = 0;
500 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconPSNLevels\n");
501 }
502 else if ( !strcmp(param2, "A11PreconPSThresh") )
503 {
504 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.PSThresh_));
505 if ( A11Params_.PSThresh_ < 0 ) A11Params_.PSThresh_ = 0;
506 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconPSThresh\n");
507 }
508 else if ( !strcmp(param2, "S22PreconPSThresh") )
509 {
510 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.PSThresh_));
511 if ( S22Params_.PSThresh_ < 0 ) S22Params_.PSThresh_ = 0;
512 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconPSThresh\n");
513 }
514 else if ( !strcmp(param2, "A11PreconPSFilter") )
515 {
516 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.PSFilter_));
517 if ( A11Params_.PSFilter_ < 0 ) A11Params_.PSFilter_ = 0;
518 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconPSFilter\n");
519 }
520 else if ( !strcmp(param2, "S22PreconPSFilter") )
521 {
522 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.PSFilter_));
523 if ( S22Params_.PSFilter_ < 0 ) S22Params_.PSFilter_ = 0;
524 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconPSFilter\n");
525 }
526 else if ( !strcmp(param2, "A11PreconAMGThresh") )
527 {
528 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.AMGThresh_));
529 if ( A11Params_.AMGThresh_ < 0.0 ) A11Params_.AMGThresh_ = 0.0;
530 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconAMGThresh\n");
531 }
532 else if ( !strcmp(param2, "S22PreconAMGThresh") )
533 {
534 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.AMGThresh_));
535 if ( S22Params_.AMGThresh_ < 0.0 ) S22Params_.AMGThresh_ = 0.0;
536 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconAMGThresh\n");
537 }
538 else if ( !strcmp(param2, "A11PreconAMGNumSweeps") )
539 {
540 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.AMGNSweeps_));
541 if ( A11Params_.AMGNSweeps_ < 0 ) A11Params_.AMGNSweeps_ = 0;
542 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconAMGNSweeps\n");
543 }
544 else if ( !strcmp(param2, "S22PreconAMGNumSweeps") )
545 {
546 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.AMGNSweeps_));
547 if ( S22Params_.AMGNSweeps_ < 0 ) S22Params_.AMGNSweeps_ = 0;
548 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconAMGNSweeps\n");
549 }
550 else if ( !strcmp(param2, "A11PreconAMGSystemSize") )
551 {
552 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.AMGSystemSize_));
553 if ( A11Params_.AMGSystemSize_ < 1 ) A11Params_.AMGSystemSize_ = 1;
554 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconAMGSystemSize\n");
555 }
556 else if ( !strcmp(param2, "S22PreconAMGSystemSize") )
557 {
558 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.AMGSystemSize_));
559 if ( S22Params_.AMGSystemSize_ < 1 ) S22Params_.AMGSystemSize_ = 1;
560 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconAMGSystemSize\n");
561 }
562 else if ( !strcmp(param2, "A11PreconEuclidNLevels") )
563 {
564 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.EuclidNLevels_));
565 if ( A11Params_.EuclidNLevels_ < 0 ) A11Params_.EuclidNLevels_ = 0;
566 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconEuclidNLevels\n");
567 }
568 else if ( !strcmp(param2, "S22PreconEuclidNLevels") )
569 {
570 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.EuclidNLevels_));
571 if ( S22Params_.EuclidNLevels_ < 0 ) S22Params_.EuclidNLevels_ = 0;
572 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconEuclidNLevels\n");
573 }
574 else if ( !strcmp(param2, "A11PreconEuclidThresh") )
575 {
576 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.EuclidThresh_));
577 if ( A11Params_.EuclidThresh_ < 0 ) A11Params_.EuclidThresh_ = 0;
578 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconEuclidThresh\n");
579 }
580 else if ( !strcmp(param2, "S22PreconEuclidThresh") )
581 {
582 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.EuclidThresh_));
583 if ( S22Params_.EuclidThresh_ < 0 ) S22Params_.EuclidThresh_ = 0;
584 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconEuclidThresh\n");
585 }
586 else if ( !strcmp(param2, "A11PreconPilutFillin") )
587 {
588 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.PilutFillin_));
589 if ( A11Params_.PilutFillin_ < 0 ) A11Params_.PilutFillin_ = 0;
590 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconPilutFillin\n");
591 }
592 else if ( !strcmp(param2, "S22PreconPilutFillin") )
593 {
594 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.PilutFillin_));
595 if ( S22Params_.PilutFillin_ < 0 ) S22Params_.PilutFillin_ = 0;
596 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconPilutFillin\n");
597 }
598 else if ( !strcmp(param2, "A11PreconPilutDropTol") )
599 {
600 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.PilutDropTol_));
601 if ( A11Params_.PilutDropTol_ < 0 ) A11Params_.PilutDropTol_ = 0;
602 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconPilutDropTol\n");
603 }
604 else if ( !strcmp(param2, "S22PreconPilutDropTol") )
605 {
606 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.PilutDropTol_));
607 if ( S22Params_.PilutDropTol_ < 0 ) S22Params_.PilutDropTol_ = 0;
608 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconPilutDropTol\n");
609 }
610 else if ( !strcmp(param2, "A11PreconMLIThresh") )
611 {
612 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.MLIThresh_));
613 if ( A11Params_.MLIThresh_ < 0.0 ) A11Params_.MLIThresh_ = 0.0;
614 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconMLIThresh\n");
615 }
616 else if ( !strcmp(param2, "S22PreconMLIThresh") )
617 {
618 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.MLIThresh_));
619 if ( S22Params_.MLIThresh_ < 0.0 ) S22Params_.MLIThresh_ = 0.0;
620 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconMLIThresh\n");
621 }
622 else if ( !strcmp(param2, "A11PreconMLINumSweeps") )
623 {
624 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.MLINSweeps_));
625 if ( A11Params_.MLINSweeps_ < 0 ) A11Params_.MLINSweeps_ = 0;
626 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconMLINSweeps\n");
627 }
628 else if ( !strcmp(param2, "S22PreconMLINumSweeps") )
629 {
630 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.MLINSweeps_));
631 if ( S22Params_.MLINSweeps_ < 0 ) S22Params_.MLINSweeps_ = 0;
632 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconMLINSweeps\n");
633 }
634 else if ( !strcmp(param2, "A11PreconMLIPweight") )
635 {
636 sscanf(params,"%s %s %lg", param1, param2, &(A11Params_.MLIPweight_));
637 if ( A11Params_.MLIPweight_ < 0.0 ) A11Params_.MLIPweight_ = 0.0;
638 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconMLIPweight\n");
639 }
640 else if ( !strcmp(param2, "S22PreconMLIPweight") )
641 {
642 sscanf(params,"%s %s %lg", param1, param2, &(S22Params_.MLIPweight_));
643 if ( S22Params_.MLIPweight_ < 0.0 ) S22Params_.MLIPweight_ = 0.0;
644 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconMLIPweight\n");
645 }
646 else if ( !strcmp(param2, "A11PreconMLINodeDOF") )
647 {
648 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.MLINodeDOF_));
649 if ( A11Params_.MLINodeDOF_ < 1 ) A11Params_.MLINodeDOF_ = 1;
650 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconMLINodeDOF\n");
651 }
652 else if ( !strcmp(param2, "S22PreconMLINodeDOF") )
653 {
654 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.MLINodeDOF_));
655 if ( S22Params_.MLINodeDOF_ < 1 ) S22Params_.MLINodeDOF_ = 1;
656 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconMLINodeDOF\n");
657 }
658 else if ( !strcmp(param2, "A11PreconMLINullDim") )
659 {
660 sscanf(params,"%s %s %d", param1, param2, &(A11Params_.MLINullDim_));
661 if ( A11Params_.MLINullDim_ < 1 ) A11Params_.MLINullDim_ = 1;
662 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::A11PreconMLINullDim\n");
663 }
664 else if ( !strcmp(param2, "S22PreconMLINullDim") )
665 {
666 sscanf(params,"%s %s %d", param1, param2, &(S22Params_.MLINullDim_));
667 if ( S22Params_.MLINullDim_ < 1 ) S22Params_.MLINullDim_ = 1;
668 if (outputLevel_ > 0) printf("HYPRE_LSI_Uzawa::S22PreconMLINullDim\n");
669 }
670 else
671 {
672 printf("HYPRE_LSI_Uzawa:: string not recognized %s\n", params);
673 }
674 return 0;
675 }
676
677 //***************************************************************************
678 // set maximum number of iterations
679 //---------------------------------------------------------------------------
680
setMaxIterations(int niter)681 int HYPRE_LSI_Uzawa::setMaxIterations(int niter)
682 {
683 maxIterations_ = niter;
684 return 0;
685 }
686
687 //***************************************************************************
688 // set tolerance
689 //---------------------------------------------------------------------------
690
setTolerance(double tol)691 int HYPRE_LSI_Uzawa::setTolerance(double tol)
692 {
693 tolerance_ = tol;
694 return 0;
695 }
696
697 //***************************************************************************
698 // get number of iterations
699 //---------------------------------------------------------------------------
700
getNumIterations(int & iter)701 int HYPRE_LSI_Uzawa::getNumIterations(int &iter)
702 {
703 iter = numIterations_;
704 return 0;
705 }
706
707 //***************************************************************************
708 // Given the matrix (A) within the object, separate the blocks
709 //---------------------------------------------------------------------------
710
setup(HYPRE_ParCSRMatrix A,HYPRE_ParVector x,HYPRE_ParVector b)711 int HYPRE_LSI_Uzawa::setup(HYPRE_ParCSRMatrix A, HYPRE_ParVector x,
712 HYPRE_ParVector b)
713 {
714 int mypid;
715
716 //------------------------------------------------------------------
717 // initial set up
718 //------------------------------------------------------------------
719
720 MPI_Comm_rank( mpiComm_, &mypid );
721 if ( mypid == 0 && outputLevel_ >= 1 )
722 printf("%4d : HYPRE_LSI_Uzawa begins....\n", mypid);
723
724 Amat_ = A;
725
726 //------------------------------------------------------------------
727 // clean up first
728 //------------------------------------------------------------------
729
730 if ( procA22Sizes_ != NULL ) delete [] procA22Sizes_;
731 if ( A11mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(A11mat_);
732 if ( A12mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(A12mat_);
733 if ( S22mat_ != NULL ) HYPRE_ParCSRMatrixDestroy(S22mat_);
734 procA22Sizes_ = NULL;
735 A11mat_ = NULL;
736 A12mat_ = NULL;
737 S22mat_ = NULL;
738
739 //------------------------------------------------------------------
740 // find the size of A22 block in the local processor (procA22Sizes_)
741 //------------------------------------------------------------------
742
743 if ( findA22BlockSize() == 0 ) return 0;
744
745 //------------------------------------------------------------------
746 // build the reduced matrix
747 //------------------------------------------------------------------
748
749 buildBlockMatrices();
750
751 //------------------------------------------------------------------
752 // setup preconditioners
753 //------------------------------------------------------------------
754
755 setupPrecon(&A11Precond_, A11mat_, A11Params_);
756 setupPrecon(&S22Precond_, S22mat_, S22Params_);
757
758 //------------------------------------------------------------------
759 // return
760 //------------------------------------------------------------------
761
762 if ( mypid == 0 && outputLevel_ >= 1 )
763 printf("%4d : HYPRE_LSI_Uzawa ends.\n", mypid);
764 return 0;
765 }
766
767 //***************************************************************************
768 // Solve using the Uzawa algorithm
769 //---------------------------------------------------------------------------
770
solve(HYPRE_ParVector b,HYPRE_ParVector x)771 int HYPRE_LSI_Uzawa::solve(HYPRE_ParVector b, HYPRE_ParVector x)
772 {
773 int mypid, *procNRows, startRow, endRow, localNRows, ierr;
774 int irow, A22NRows;
775 double rnorm0, rnorm;
776 double *b_data, *x_data, *u1_data, *u2_data, *f1_data, *f2_data;
777 double *v1_data;
778 HYPRE_IJVector IJR, IJF1, IJF2, IJU1, IJU2, IJT1, IJT2, IJV1, IJV2;
779 HYPRE_ParVector r_csr, f1_csr, f2_csr, u1_csr, u2_csr, t1_csr, t2_csr;
780 HYPRE_ParVector v1_csr, v2_csr;
781 hypre_Vector *b_local, *f1_local, *f2_local;
782 hypre_Vector *x_local, *u1_local, *u2_local, *v1_local;
783
784 //------------------------------------------------------------------
785 // get machine information
786 //------------------------------------------------------------------
787
788 MPI_Comm_rank( mpiComm_, &mypid );
789
790 //------------------------------------------------------------------
791 // create temporary vectors
792 //------------------------------------------------------------------
793
794 HYPRE_ParCSRMatrixGetRowPartitioning( Amat_, &procNRows );
795 startRow = procNRows[mypid];
796 endRow = procNRows[mypid+1] - 1;
797 localNRows = endRow - startRow + 1;
798 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJR);
799 ierr = HYPRE_IJVectorSetObjectType(IJR, HYPRE_PARCSR);
800 ierr = HYPRE_IJVectorInitialize(IJR);
801 ierr = HYPRE_IJVectorAssemble(IJR);
802 hypre_assert(!ierr);
803 HYPRE_IJVectorGetObject(IJR, (void **) &r_csr);
804
805 startRow = procNRows[mypid] - procA22Sizes_[mypid];
806 endRow = procNRows[mypid+1] - procA22Sizes_[mypid+1] - 1;
807 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJF1);
808 ierr = HYPRE_IJVectorSetObjectType(IJF1, HYPRE_PARCSR);
809 ierr = HYPRE_IJVectorInitialize(IJF1);
810 ierr = HYPRE_IJVectorAssemble(IJF1);
811 hypre_assert(!ierr);
812 HYPRE_IJVectorGetObject(IJF1, (void **) &f1_csr);
813 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJU1);
814 ierr = HYPRE_IJVectorSetObjectType(IJU1, HYPRE_PARCSR);
815 ierr = HYPRE_IJVectorInitialize(IJU1);
816 ierr = HYPRE_IJVectorAssemble(IJU1);
817 hypre_assert(!ierr);
818 HYPRE_IJVectorGetObject(IJU1, (void **) &u1_csr);
819 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJT1);
820 ierr = HYPRE_IJVectorSetObjectType(IJT1, HYPRE_PARCSR);
821 ierr = HYPRE_IJVectorInitialize(IJT1);
822 ierr = HYPRE_IJVectorAssemble(IJT1);
823 hypre_assert(!ierr);
824 HYPRE_IJVectorGetObject(IJT1, (void **) &t1_csr);
825 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJV1);
826 ierr = HYPRE_IJVectorSetObjectType(IJV1, HYPRE_PARCSR);
827 ierr = HYPRE_IJVectorInitialize(IJV1);
828 ierr = HYPRE_IJVectorAssemble(IJV1);
829 hypre_assert(!ierr);
830 HYPRE_IJVectorGetObject(IJV1, (void **) &v1_csr);
831
832 startRow = procA22Sizes_[mypid];
833 endRow = procA22Sizes_[mypid+1] - 1;
834 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJF2);
835 ierr = HYPRE_IJVectorSetObjectType(IJF2, HYPRE_PARCSR);
836 ierr = HYPRE_IJVectorInitialize(IJF2);
837 ierr = HYPRE_IJVectorAssemble(IJF2);
838 hypre_assert(!ierr);
839 HYPRE_IJVectorGetObject(IJF2, (void **) &f2_csr);
840 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJU2);
841 ierr = HYPRE_IJVectorSetObjectType(IJU2, HYPRE_PARCSR);
842 ierr = HYPRE_IJVectorInitialize(IJU2);
843 ierr = HYPRE_IJVectorAssemble(IJU2);
844 hypre_assert(!ierr);
845 HYPRE_IJVectorGetObject(IJU2, (void **) &u2_csr);
846 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJT2);
847 ierr = HYPRE_IJVectorSetObjectType(IJT2, HYPRE_PARCSR);
848 ierr = HYPRE_IJVectorInitialize(IJT2);
849 ierr = HYPRE_IJVectorAssemble(IJT2);
850 hypre_assert(!ierr);
851 HYPRE_IJVectorGetObject(IJT2, (void **) &t2_csr);
852 ierr = HYPRE_IJVectorCreate(mpiComm_, startRow, endRow, &IJV2);
853 ierr = HYPRE_IJVectorSetObjectType(IJV2, HYPRE_PARCSR);
854 ierr = HYPRE_IJVectorInitialize(IJV2);
855 ierr = HYPRE_IJVectorAssemble(IJV2);
856 hypre_assert(!ierr);
857 HYPRE_IJVectorGetObject(IJV2, (void **) &v2_csr);
858 free( procNRows );
859
860 //------------------------------------------------------------------
861 // compute initial residual
862 //------------------------------------------------------------------
863
864 if ( maxIterations_ > 1 )
865 {
866 HYPRE_ParVectorCopy( b, r_csr );
867 HYPRE_ParCSRMatrixMatvec( -1.0, Amat_, x, 1.0, r_csr );
868 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
869 rnorm = sqrt( rnorm );
870 rnorm0 = rnorm * tolerance_;
871 if ( rnorm < rnorm0 ) return 0;
872 if ( mypid == 0 ) printf("Uzawa : initial rnorm = %e\n",rnorm);
873 }
874 else rnorm = rnorm0 = 1.0;
875
876 //------------------------------------------------------------------
877 // set up solvers
878 //------------------------------------------------------------------
879
880 if ( A11Solver_ == NULL )
881 setupSolver(&A11Solver_,A11mat_,f1_csr,u1_csr,A11Precond_,A11Params_);
882 if ( S22Params_.SolverID_ != 0 && S22Solver_ == NULL )
883 setupSolver(&S22Solver_,S22mat_,f2_csr,u2_csr,S22Precond_,S22Params_);
884
885 //------------------------------------------------------------------
886 // distribute the vectors
887 //------------------------------------------------------------------
888
889 A22NRows = procA22Sizes_[mypid+1] - procA22Sizes_[mypid];
890 b_local = hypre_ParVectorLocalVector((hypre_ParVector *) b);
891 b_data = (double *) hypre_VectorData(b_local);
892 f1_local = hypre_ParVectorLocalVector((hypre_ParVector *) f1_csr);
893 f1_data = (double *) hypre_VectorData(f1_local);
894 f2_local = hypre_ParVectorLocalVector((hypre_ParVector *) f2_csr);
895 f2_data = (double *) hypre_VectorData(f2_local);
896 x_local = hypre_ParVectorLocalVector((hypre_ParVector *) x);
897 x_data = (double *) hypre_VectorData(x_local);
898 u1_local = hypre_ParVectorLocalVector((hypre_ParVector *) u1_csr);
899 u1_data = (double *) hypre_VectorData(u1_local);
900 u2_local = hypre_ParVectorLocalVector((hypre_ParVector *) u2_csr);
901 u2_data = (double *) hypre_VectorData(u2_local);
902 v1_local = hypre_ParVectorLocalVector((hypre_ParVector *) v1_csr);
903 v1_data = (double *) hypre_VectorData(v1_local);
904
905 for ( irow = 0; irow < localNRows-A22NRows; irow++ )
906 f1_data[irow] = b_data[irow];
907 for ( irow = localNRows-A22NRows; irow < localNRows; irow++ )
908 f2_data[irow-localNRows+A22NRows] = b_data[irow];
909
910 //------------------------------------------------------------------
911 // iterate (using the Cheng-Zou method)
912 //------------------------------------------------------------------
913
914 numIterations_ = 0;
915
916 while ( numIterations_ < maxIterations_ && rnorm >= rnorm0 )
917 {
918 numIterations_++;
919
920 //----------------------------------------------------------------
921 // copy x to split vectors
922 //---------------------------------------------------------------
923
924 for ( irow = 0; irow < localNRows-A22NRows; irow++ )
925 v1_data[irow] = u1_data[irow] = x_data[irow];
926 for ( irow = localNRows-A22NRows; irow < localNRows; irow++ )
927 u2_data[irow-localNRows+A22NRows] = x_data[irow];
928
929 //----------------------------------------------------------------
930 // x_{i+1/2} = x_i + Q_A^{-1} (f1 - A x_i - A12 y_i)
931 //---------------------------------------------------------------
932
933 HYPRE_ParVectorCopy( f1_csr, t1_csr );
934 HYPRE_ParCSRMatrixMatvec( -1.0, A11mat_, u1_csr, 1.0, t1_csr );
935 HYPRE_ParCSRMatrixMatvec( -1.0, A12mat_, u2_csr, 1.0, t1_csr );
936
937 if ( A11Params_.SolverID_ == 1 )
938 HYPRE_ParCSRPCGSolve(A11Solver_, A11mat_, t1_csr, u1_csr);
939 else if ( A11Params_.SolverID_ == 2 )
940 HYPRE_ParCSRGMRESSolve(A11Solver_, A11mat_, t1_csr, u1_csr);
941
942 hypre_ParVectorAxpy( 1.0, (hypre_ParVector*)v1_csr ,
943 (hypre_ParVector*)u1_csr );
944
945 //----------------------------------------------------------------
946 // y_{i+1/2} = y_i + Q_B^{-1} (A21 x_{i+1/2} - f2)
947 //---------------------------------------------------------------
948
949 if ( modifiedScheme_ >= 1 )
950 {
951 HYPRE_ParVectorCopy( f2_csr, t2_csr );
952 HYPRE_ParCSRMatrixMatvecT( 1.0, A12mat_, u1_csr, -1.0, t2_csr );
953
954 if ( S22Params_.SolverID_ == 1 )
955 HYPRE_ParCSRPCGSolve(S22Solver_, S22mat_, t2_csr, v2_csr);
956 else if ( S22Params_.SolverID_ == 2 )
957 HYPRE_ParCSRGMRESSolve(S22Solver_, S22mat_, t2_csr, v2_csr);
958 else
959 {
960 HYPRE_ParVectorCopy( t2_csr, v2_csr );
961 HYPRE_ParVectorScale( S22SolverDampFactor_, v2_csr );
962 }
963
964 hypre_ParVectorAxpy( 1.0, (hypre_ParVector*)v2_csr ,
965 (hypre_ParVector*)u2_csr );
966
967 //----------------------------------------------------------------
968 // x_{i+1} = x_i + Q_A^{-1} (f - A x_i - A12 y_{i+1/2})
969 //---------------------------------------------------------------
970
971 HYPRE_ParVectorCopy( f1_csr, t1_csr );
972 HYPRE_ParCSRMatrixMatvec( -1.0, A11mat_, v1_csr, 1.0, t1_csr );
973 HYPRE_ParCSRMatrixMatvec( -1.0, A12mat_, u2_csr, 1.0, t1_csr );
974
975 if ( A11Params_.SolverID_ == 1 )
976 HYPRE_ParCSRPCGSolve(A11Solver_, A11mat_, t1_csr, u1_csr);
977 else if ( A11Params_.SolverID_ == 2 )
978 HYPRE_ParCSRGMRESSolve(A11Solver_, A11mat_, t1_csr, u1_csr);
979
980 hypre_ParVectorAxpy( 1.0, (hypre_ParVector*)v1_csr ,
981 (hypre_ParVector*)u1_csr );
982
983 //----------------------------------------------------------------
984 // y_{i+1} = y_{i+1/2} + Q_B^{-1} (A21 x_{i+1} - f2)
985 //---------------------------------------------------------------
986
987 HYPRE_ParVectorCopy( f2_csr, t2_csr );
988 HYPRE_ParCSRMatrixMatvecT( 1.0, A12mat_, u1_csr, -1.0, t2_csr );
989
990 if ( S22Params_.SolverID_ == 1 )
991 HYPRE_ParCSRPCGSolve(S22Solver_, S22mat_, t2_csr, v2_csr);
992 else if ( S22Params_.SolverID_ == 2 )
993 HYPRE_ParCSRGMRESSolve(S22Solver_, S22mat_, t2_csr, v2_csr);
994 else
995 {
996 HYPRE_ParVectorCopy( t2_csr, v2_csr );
997 HYPRE_ParVectorScale( S22SolverDampFactor_, v2_csr );
998 }
999
1000 hypre_ParVectorAxpy( 1.0, (hypre_ParVector*)v2_csr ,
1001 (hypre_ParVector*)u2_csr );
1002 }
1003
1004 //----------------------------------------------------------------
1005 // merge solution vector and compute residual norm
1006 //---------------------------------------------------------------
1007
1008 for ( irow = 0; irow < localNRows-A22NRows; irow++ )
1009 x_data[irow] = u1_data[irow];
1010 for ( irow = localNRows-A22NRows; irow < localNRows; irow++ )
1011 x_data[irow] = u2_data[irow-localNRows+A22NRows];
1012
1013 if ( maxIterations_ > 1 )
1014 {
1015 HYPRE_ParVectorCopy( b, r_csr );
1016 HYPRE_ParCSRMatrixMatvec( -1.0, Amat_, x, 1.0, r_csr );
1017 HYPRE_ParVectorInnerProd( r_csr, r_csr, &rnorm);
1018 rnorm = sqrt( rnorm );
1019 if ( mypid == 0 )
1020 printf("Uzawa : iteration = %5d, rnorm = %e\n",numIterations_,
1021 rnorm);
1022 }
1023 }
1024 HYPRE_IJVectorDestroy(IJR);
1025 HYPRE_IJVectorDestroy(IJF1);
1026 HYPRE_IJVectorDestroy(IJF2);
1027 HYPRE_IJVectorDestroy(IJU1);
1028 HYPRE_IJVectorDestroy(IJU2);
1029 HYPRE_IJVectorDestroy(IJT1);
1030 HYPRE_IJVectorDestroy(IJT2);
1031 HYPRE_IJVectorDestroy(IJV2);
1032 return 0;
1033 }
1034
1035 //***************************************************************************
1036 // search for the separator for
1037 // A = |A_11 A_12|
1038 // |A_12^T A_22|
1039 //---------------------------------------------------------------------------
1040
findA22BlockSize()1041 int HYPRE_LSI_Uzawa::findA22BlockSize()
1042 {
1043 int mypid, nprocs, *procNRows, startRow, endRow;
1044 int A22LocalSize, irow, zeroDiag, jcol, rowSize, *colInd;
1045 int *iTempList, ip, ncnt, A22GlobalSize;
1046 double *colVal;
1047
1048 //------------------------------------------------------------------
1049 // get matrix information
1050 //------------------------------------------------------------------
1051
1052 MPI_Comm_rank(mpiComm_, &mypid);
1053 MPI_Comm_size(mpiComm_, &nprocs);
1054 HYPRE_ParCSRMatrixGetRowPartitioning( Amat_, &procNRows );
1055 startRow = procNRows[mypid];
1056 endRow = procNRows[mypid+1] - 1;
1057 free( procNRows );
1058
1059 //------------------------------------------------------------------
1060 // search for dimension of A_22
1061 //------------------------------------------------------------------
1062
1063 A22LocalSize = 0;
1064 for ( irow = endRow; irow >= startRow; irow-- )
1065 {
1066 HYPRE_ParCSRMatrixGetRow(Amat_,irow,&rowSize,&colInd,&colVal);
1067 zeroDiag = 1;
1068 for ( jcol = 0; jcol < rowSize; jcol++ )
1069 {
1070 if ( colInd[jcol] == irow && colVal[jcol] != 0.0 )
1071 {
1072 zeroDiag = 0;
1073 break;
1074 }
1075 }
1076 HYPRE_ParCSRMatrixRestoreRow(Amat_,irow,&rowSize,&colInd,&colVal);
1077 if ( zeroDiag ) A22LocalSize++;
1078 else break;
1079 }
1080 if ( outputLevel_ >= 1 )
1081 printf("%4d : findA22BlockSize - local nrows = %d\n",mypid,A22LocalSize);
1082
1083 //------------------------------------------------------------------
1084 // gather the block size information on all processors
1085 //------------------------------------------------------------------
1086
1087 iTempList = new int[nprocs];
1088 if ( procA22Sizes_ != NULL ) delete [] procA22Sizes_;
1089 procA22Sizes_ = new int[nprocs+1];
1090 for ( ip = 0; ip < nprocs; ip++ ) iTempList[ip] = 0;
1091 iTempList[mypid] = A22LocalSize;
1092 MPI_Allreduce(iTempList,procA22Sizes_,nprocs,MPI_INT,MPI_SUM,mpiComm_);
1093 delete [] iTempList;
1094 A22GlobalSize = 0;
1095 ncnt = 0;
1096 for ( ip = 0; ip < nprocs; ip++ )
1097 {
1098 ncnt = procA22Sizes_[ip];
1099 procA22Sizes_[ip] = A22GlobalSize;
1100 A22GlobalSize += ncnt;
1101 }
1102 procA22Sizes_[nprocs] = A22GlobalSize;
1103 return A22GlobalSize;
1104 }
1105
1106 //****************************************************************************
1107 // build the block matrices A11, A21, and S22
1108 //----------------------------------------------------------------------------
1109
buildBlockMatrices()1110 int HYPRE_LSI_Uzawa::buildBlockMatrices()
1111 {
1112 int ierr=0;
1113
1114 ierr += buildA11A12Mat();
1115 ierr += buildS22Mat();
1116 return ierr;
1117 }
1118
1119 //****************************************************************************
1120 // build A11 and A12 matrix
1121 //----------------------------------------------------------------------------
1122
buildA11A12Mat()1123 int HYPRE_LSI_Uzawa::buildA11A12Mat()
1124 {
1125 int mypid, nprocs, *procNRows, startRow, endRow, newEndRow, ierr;
1126 int A12NCols, A11NRows, A11StartRow, A12StartCol, *A11MatSize, ip;
1127 int *A12MatSize, irow, jcol, colIndex, uBound, A11RowSize, A12RowSize;
1128 int *A11ColInd, *A12ColInd, rowIndex, rowSize, *colInd, ncnt;
1129 int localNRows, maxA11RowSize, maxA12RowSize;
1130 double *colVal, *A11ColVal, *A12ColVal;
1131 HYPRE_IJMatrix IJA11, IJA12;
1132
1133 //------------------------------------------------------------------
1134 // get matrix information
1135 //------------------------------------------------------------------
1136
1137 MPI_Comm_rank(mpiComm_, &mypid);
1138 MPI_Comm_size(mpiComm_, &nprocs);
1139 HYPRE_ParCSRMatrixGetRowPartitioning( Amat_, &procNRows );
1140 startRow = procNRows[mypid];
1141 endRow = procNRows[mypid+1] - 1;
1142 localNRows = endRow - startRow + 1;
1143 newEndRow = endRow - (procA22Sizes_[mypid+1] - procA22Sizes_[mypid]);
1144
1145 //------------------------------------------------------------------
1146 // calculate the dimension of A11 and A12
1147 //------------------------------------------------------------------
1148
1149 A12NCols = procA22Sizes_[mypid+1] - procA22Sizes_[mypid];
1150 A11NRows = localNRows - A12NCols;
1151 A11StartRow = procNRows[mypid] - procA22Sizes_[mypid];
1152 A12StartCol = procA22Sizes_[mypid];
1153
1154 if ( outputLevel_ >= 1 )
1155 {
1156 printf("%4d : buildA11A12Mat - A11StartRow = %d\n", mypid, A11StartRow);
1157 printf("%4d : buildA11A12Mat - A11LocalDim = %d\n", mypid, A11NRows);
1158 printf("%4d : buildA11A12Mat - A12StartRow = %d\n", mypid, A12StartCol);
1159 printf("%4d : buildA11A12Mat - A12LocalCol = %d\n", mypid, A12NCols);
1160 }
1161
1162 //------------------------------------------------------------------
1163 // create a matrix context for A11 and A12
1164 //------------------------------------------------------------------
1165
1166 ierr = HYPRE_IJMatrixCreate(mpiComm_,A11StartRow,A11StartRow+A11NRows-1,
1167 A11StartRow,A11StartRow+A11NRows-1,&IJA11);
1168 ierr += HYPRE_IJMatrixSetObjectType(IJA11, HYPRE_PARCSR);
1169 hypre_assert(!ierr);
1170 ierr = HYPRE_IJMatrixCreate(mpiComm_,A11StartRow,A11StartRow+A11NRows-1,
1171 A12StartCol,A12StartCol+A12NCols-1,&IJA12);
1172 ierr += HYPRE_IJMatrixSetObjectType(IJA12, HYPRE_PARCSR);
1173 hypre_assert(!ierr);
1174
1175 //------------------------------------------------------------------
1176 // compute the number of nonzeros in each matrix
1177 //------------------------------------------------------------------
1178
1179 A11MatSize = new int[A11NRows];
1180 A12MatSize = new int[A11NRows];
1181 maxA11RowSize = maxA12RowSize = 0;
1182
1183 for ( irow = startRow; irow <= newEndRow ; irow++ )
1184 {
1185 A11RowSize = A12RowSize = 0;
1186 HYPRE_ParCSRMatrixGetRow(Amat_,irow,&rowSize,&colInd,NULL);
1187 for ( jcol = 0; jcol < rowSize; jcol++ )
1188 {
1189 colIndex = colInd[jcol];
1190 for ( ip = 1; ip <= nprocs; ip++ )
1191 if ( procNRows[ip] > colIndex ) break;
1192 uBound = procNRows[ip] - (procA22Sizes_[ip] - procA22Sizes_[ip-1]);
1193 if ( colIndex < uBound ) A11RowSize++;
1194 else A12RowSize++;
1195 }
1196 A11MatSize[irow-startRow] = A11RowSize;
1197 A12MatSize[irow-startRow] = A12RowSize;
1198 maxA11RowSize = (A11RowSize > maxA11RowSize) ? A11RowSize : maxA11RowSize;
1199 maxA12RowSize = (A12RowSize > maxA12RowSize) ? A12RowSize : maxA12RowSize;
1200 HYPRE_ParCSRMatrixRestoreRow(Amat_,irow,&rowSize,&colInd,NULL);
1201 }
1202
1203 //------------------------------------------------------------------
1204 // after fetching the row sizes, initialize the matrices
1205 //------------------------------------------------------------------
1206
1207 ierr = HYPRE_IJMatrixSetRowSizes(IJA11, A11MatSize);
1208 ierr += HYPRE_IJMatrixInitialize(IJA11);
1209 hypre_assert(!ierr);
1210 ierr = HYPRE_IJMatrixSetRowSizes(IJA12, A12MatSize);
1211 ierr += HYPRE_IJMatrixInitialize(IJA12);
1212 hypre_assert(!ierr);
1213
1214 //------------------------------------------------------------------
1215 // next load the matrices
1216 //------------------------------------------------------------------
1217
1218 A11ColInd = new int[maxA11RowSize+1];
1219 A11ColVal = new double[maxA11RowSize+1];
1220 A12ColInd = new int[maxA12RowSize+1];
1221 A12ColVal = new double[maxA12RowSize+1];
1222
1223 for ( irow = startRow; irow <= newEndRow ; irow++ )
1224 {
1225 A11RowSize = A12RowSize = 0;
1226 HYPRE_ParCSRMatrixGetRow(Amat_,irow,&rowSize,&colInd,&colVal);
1227 for ( jcol = 0; jcol < rowSize; jcol++ )
1228 {
1229 colIndex = colInd[jcol];
1230 for ( ip = 1; ip <= nprocs; ip++ )
1231 if ( procNRows[ip] > colIndex ) break;
1232 uBound = procNRows[ip] - (procA22Sizes_[ip] - procA22Sizes_[ip-1]);
1233 if ( colIndex < uBound )
1234 {
1235 A11ColInd[A11RowSize] = colIndex - procA22Sizes_[ip-1];
1236 A11ColVal[A11RowSize++] = colVal[jcol];
1237 }
1238 else
1239 {
1240 A12ColInd[A12RowSize] = colIndex - uBound + procA22Sizes_[ip-1];
1241 A12ColVal[A12RowSize++] = colVal[jcol];
1242 }
1243 }
1244 HYPRE_ParCSRMatrixRestoreRow(Amat_,irow,&rowSize,&colInd,&colVal);
1245 rowIndex = irow - procA22Sizes_[mypid];
1246 ierr = HYPRE_IJMatrixSetValues(IJA11, 1, &A11RowSize,
1247 (const int *) &rowIndex, (const int *) A11ColInd,
1248 (const double *) A11ColVal);
1249 hypre_assert( !ierr );
1250 ierr = HYPRE_IJMatrixSetValues(IJA12, 1, &A12RowSize,
1251 (const int *) &rowIndex, (const int *) A12ColInd,
1252 (const double *) A12ColVal);
1253 hypre_assert( !ierr );
1254 }
1255
1256 //------------------------------------------------------------------
1257 // finally assemble the matrix and sanitize
1258 //------------------------------------------------------------------
1259
1260 HYPRE_IJMatrixAssemble(IJA11);
1261 HYPRE_IJMatrixGetObject(IJA11, (void **) &A11mat_);
1262 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A11mat_);
1263 HYPRE_IJMatrixAssemble(IJA12);
1264 HYPRE_IJMatrixGetObject(IJA12, (void **) &A12mat_);
1265 hypre_MatvecCommPkgCreate((hypre_ParCSRMatrix *) A11mat_);
1266
1267 HYPRE_IJMatrixSetObjectType(IJA11, -1);
1268 HYPRE_IJMatrixDestroy(IJA11);
1269 HYPRE_IJMatrixSetObjectType(IJA12, -1);
1270 HYPRE_IJMatrixDestroy(IJA12);
1271 delete [] A11MatSize;
1272 delete [] A12MatSize;
1273 delete [] A11ColInd;
1274 delete [] A11ColVal;
1275 delete [] A12ColInd;
1276 delete [] A12ColVal;
1277 free( procNRows );
1278
1279 //------------------------------------------------------------------
1280 // diagnostics
1281 //------------------------------------------------------------------
1282
1283 if ( outputLevel_ >= 3 )
1284 {
1285 ncnt = 0;
1286 MPI_Barrier(mpiComm_);
1287 while ( ncnt < nprocs )
1288 {
1289 if ( mypid == ncnt )
1290 {
1291 printf("====================================================\n");
1292 printf("%4d : Printing A11 matrix... \n", mypid);
1293 fflush(stdout);
1294 for (irow = A11StartRow;irow < A11StartRow+A11NRows;irow++)
1295 {
1296 HYPRE_ParCSRMatrixGetRow(A11mat_,irow,&rowSize,&colInd,&colVal);
1297 for ( jcol = 0; jcol < rowSize; jcol++ )
1298 if ( colVal[jcol] != 0.0 )
1299 printf("%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
1300 colVal[jcol]);
1301 HYPRE_ParCSRMatrixRestoreRow(A11mat_, irow, &rowSize,
1302 &colInd, &colVal);
1303 }
1304 printf("====================================================\n");
1305 }
1306 ncnt++;
1307 MPI_Barrier(mpiComm_);
1308 }
1309 }
1310 if ( outputLevel_ >= 3 )
1311 {
1312 ncnt = 0;
1313 MPI_Barrier(mpiComm_);
1314 while ( ncnt < nprocs )
1315 {
1316 if ( mypid == ncnt )
1317 {
1318 printf("====================================================\n");
1319 printf("%4d : Printing A12 matrix... \n", mypid);
1320 fflush(stdout);
1321 for (irow = A11StartRow;irow < A11StartRow+A11NRows;irow++)
1322 {
1323 HYPRE_ParCSRMatrixGetRow(A12mat_,irow,&rowSize,&colInd,&colVal);
1324 for ( jcol = 0; jcol < rowSize; jcol++ )
1325 if ( colVal[jcol] != 0.0 )
1326 printf("%6d %6d %25.16e \n",irow+1,colInd[jcol]+1,
1327 colVal[jcol]);
1328 HYPRE_ParCSRMatrixRestoreRow(A12mat_, irow, &rowSize,
1329 &colInd, &colVal);
1330 }
1331 printf("====================================================\n");
1332 }
1333 ncnt++;
1334 MPI_Barrier(mpiComm_);
1335 }
1336 }
1337 return 0;
1338 }
1339
1340 //****************************************************************************
1341 // build the block matrix S22 (B diag(A)^{-1} B^T)
1342 //----------------------------------------------------------------------------
1343
buildS22Mat()1344 int HYPRE_LSI_Uzawa::buildS22Mat()
1345 {
1346 int mypid, nprocs, *procNRows, one=1, A11StartRow, A11NRows, irow, jcol;
1347 int rowSize, *colInd, *A11MatSize, ierr;
1348 double *colVal, ddata;
1349 HYPRE_ParCSRMatrix ainvA11_csr;
1350 HYPRE_IJMatrix ainvA11;
1351 HYPRE_Solver parasails;
1352
1353 //------------------------------------------------------------------
1354 // get machine information
1355 //------------------------------------------------------------------
1356
1357 MPI_Comm_rank(mpiComm_, &mypid);
1358 MPI_Comm_size(mpiComm_, &nprocs);
1359
1360 //------------------------------------------------------------------
1361 // choose between diagonal or approximate inverse
1362 //------------------------------------------------------------------
1363
1364 if ( S22Scheme_ == 1 )
1365 {
1366 //---------------------------------------------------------------
1367 // build approximate inverse of A11
1368 //---------------------------------------------------------------
1369
1370 HYPRE_ParaSailsCreate(mpiComm_, ¶sails);
1371 HYPRE_ParaSailsSetParams(parasails, 0.1, 1);
1372 HYPRE_ParaSailsSetFilter(parasails, 0.1);
1373 HYPRE_ParaSailsSetLogging(parasails, 1);
1374 HYPRE_ParaSailsSetup(parasails, A11mat_, NULL, NULL);
1375 HYPRE_ParaSailsBuildIJMatrix(parasails, &ainvA11);
1376 }
1377 else
1378 {
1379 //---------------------------------------------------------------
1380 // build inverse of diagonal of A11
1381 //---------------------------------------------------------------
1382
1383 HYPRE_ParCSRMatrixGetRowPartitioning( A11mat_, &procNRows );
1384 A11StartRow = procNRows[mypid];
1385 A11NRows = procNRows[mypid+1] - A11StartRow;
1386
1387 ierr = HYPRE_IJMatrixCreate(mpiComm_,A11StartRow,
1388 A11StartRow+A11NRows-1, A11StartRow,
1389 A11StartRow+A11NRows-1,&ainvA11);
1390 ierr += HYPRE_IJMatrixSetObjectType(ainvA11, HYPRE_PARCSR);
1391 hypre_assert(!ierr);
1392
1393 A11MatSize = new int[A11NRows];
1394 for ( irow = 0; irow < A11NRows; irow++ ) A11MatSize[irow] = 1;
1395 ierr = HYPRE_IJMatrixSetRowSizes(ainvA11, A11MatSize);
1396 ierr += HYPRE_IJMatrixInitialize(ainvA11);
1397 hypre_assert(!ierr);
1398
1399 for ( irow = A11StartRow; irow < A11StartRow+A11NRows; irow++ )
1400 {
1401 HYPRE_ParCSRMatrixGetRow(A11mat_,irow,&rowSize,&colInd,&colVal);
1402 ddata = 0.0;
1403 for ( jcol = 0; jcol < rowSize; jcol++ )
1404 {
1405 if ( colInd[jcol] == irow )
1406 {
1407 ddata = 1.0 / colVal[jcol];
1408 break;
1409 }
1410 }
1411 HYPRE_ParCSRMatrixRestoreRow(A11mat_,irow,&rowSize,&colInd,&colVal);
1412 ierr = HYPRE_IJMatrixSetValues(ainvA11, 1, &one, (const int *) &irow,
1413 (const int *) &irow, (const double *) &ddata);
1414 hypre_assert( !ierr );
1415 }
1416 HYPRE_IJMatrixAssemble(ainvA11);
1417 free( procNRows );
1418 delete [] A11MatSize;
1419 }
1420
1421 //------------------------------------------------------------------
1422 // perform the triple matrix product A12' * diagA11 * A12
1423 //------------------------------------------------------------------
1424
1425 HYPRE_IJMatrixGetObject(ainvA11, (void **) &ainvA11_csr);
1426 hypre_BoomerAMGBuildCoarseOperator((hypre_ParCSRMatrix *) A12mat_,
1427 (hypre_ParCSRMatrix *) ainvA11_csr,
1428 (hypre_ParCSRMatrix *) A12mat_,
1429 (hypre_ParCSRMatrix **) &S22mat_);
1430
1431 //------------------------------------------------------------------
1432 // clean up and return
1433 //------------------------------------------------------------------
1434
1435 HYPRE_IJMatrixDestroy(ainvA11);
1436 return 0;
1437 }
1438
1439 //***************************************************************************
1440 // setup preconditioner
1441 //---------------------------------------------------------------------------
1442
setupPrecon(HYPRE_Solver * precon,HYPRE_ParCSRMatrix Amat,HYPRE_Uzawa_PARAMS paramPtr)1443 int HYPRE_LSI_Uzawa::setupPrecon(HYPRE_Solver *precon,HYPRE_ParCSRMatrix Amat,
1444 HYPRE_Uzawa_PARAMS paramPtr)
1445 {
1446 int i, *nsweeps, *relaxType;
1447 char **targv;
1448 #ifdef HAVE_MLI
1449 char paramString[100];
1450 #endif
1451
1452 (void) Amat;
1453 if ( paramPtr.SolverID_ == 0 ) return 0;
1454
1455 //------------------------------------------------------------------
1456 // set up the solvers and preconditioners
1457 //------------------------------------------------------------------
1458
1459 switch( paramPtr.PrecondID_ )
1460 {
1461 case 2 :
1462 HYPRE_ParCSRParaSailsCreate( mpiComm_, precon );
1463 if (paramPtr.SolverID_ == 0) HYPRE_ParCSRParaSailsSetSym(*precon,1);
1464 else HYPRE_ParCSRParaSailsSetSym(*precon,0);
1465 HYPRE_ParCSRParaSailsSetParams(*precon, paramPtr.PSThresh_,
1466 paramPtr.PSNLevels_);
1467 HYPRE_ParCSRParaSailsSetFilter(*precon, paramPtr.PSFilter_);
1468 break;
1469 case 3 :
1470 HYPRE_BoomerAMGCreate(precon);
1471 HYPRE_BoomerAMGSetMaxIter(*precon, 1);
1472 HYPRE_BoomerAMGSetCycleType(*precon, 1);
1473 HYPRE_BoomerAMGSetPrintLevel(*precon, outputLevel_);
1474 HYPRE_BoomerAMGSetMaxLevels(*precon, 25);
1475 HYPRE_BoomerAMGSetMeasureType(*precon, 0);
1476 HYPRE_BoomerAMGSetCoarsenType(*precon, 0);
1477 HYPRE_BoomerAMGSetStrongThreshold(*precon,paramPtr.AMGThresh_);
1478 if ( paramPtr.AMGSystemSize_ > 1 )
1479 HYPRE_BoomerAMGSetNumFunctions(*precon,paramPtr.AMGSystemSize_);
1480 nsweeps = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
1481 for ( i = 0; i < 4; i++ ) nsweeps[i] = paramPtr.AMGNSweeps_;
1482 HYPRE_BoomerAMGSetNumGridSweeps(*precon, nsweeps);
1483 relaxType = hypre_CTAlloc(int,4,HYPRE_MEMORY_HOST);
1484 for ( i = 0; i < 4; i++ ) relaxType[i] = 6;
1485 HYPRE_BoomerAMGSetGridRelaxType(*precon, relaxType);
1486 break;
1487 case 4 :
1488 HYPRE_ParCSRPilutCreate( mpiComm_, precon );
1489 HYPRE_ParCSRPilutSetMaxIter( *precon, 1 );
1490 HYPRE_ParCSRPilutSetFactorRowSize(*precon,paramPtr.PilutFillin_);
1491 HYPRE_ParCSRPilutSetDropTolerance(*precon,paramPtr.PilutDropTol_);
1492 break;
1493 case 5 :
1494 HYPRE_EuclidCreate( mpiComm_, precon );
1495 targv = hypre_TAlloc(char*, 4 , HYPRE_MEMORY_HOST);
1496 for ( i = 0; i < 4; i++ ) targv[i] = hypre_TAlloc(char, 50, HYPRE_MEMORY_HOST);
1497 strcpy(targv[0], "-level");
1498 sprintf(targv[1], "%1d", paramPtr.EuclidNLevels_);
1499 strcpy(targv[2], "-sparseA");
1500 sprintf(targv[3], "%f", paramPtr.EuclidThresh_);
1501 HYPRE_EuclidSetParams(*precon, 4, targv);
1502 for ( i = 0; i < 4; i++ ) free(targv[i]);
1503 free(targv);
1504 break;
1505 case 6 :
1506 #ifdef HAVE_MLI
1507 HYPRE_LSI_MLICreate(mpiComm_, precon);
1508 sprintf(paramString, "MLI outputLevel %d", outputLevel_);
1509 HYPRE_LSI_MLISetParams(*precon, paramString);
1510 sprintf(paramString, "MLI strengthThreshold %e",paramPtr.MLIThresh_);
1511 HYPRE_LSI_MLISetParams(*precon, paramString);
1512 sprintf(paramString, "MLI method AMGSA");
1513 HYPRE_LSI_MLISetParams(*precon, paramString);
1514 sprintf(paramString, "MLI smoother SGS");
1515 HYPRE_LSI_MLISetParams(*precon, paramString);
1516 sprintf(paramString, "MLI numSweeps %d",paramPtr.MLINSweeps_);
1517 HYPRE_LSI_MLISetParams(*precon, paramString);
1518 sprintf(paramString, "MLI Pweight %e",paramPtr.MLIPweight_);
1519 HYPRE_LSI_MLISetParams(*precon, paramString);
1520 sprintf(paramString, "MLI nodeDOF %d",paramPtr.MLINodeDOF_);
1521 HYPRE_LSI_MLISetParams(*precon, paramString);
1522 sprintf(paramString, "MLI nullSpaceDim %d",paramPtr.MLINullDim_);
1523 HYPRE_LSI_MLISetParams(*precon, paramString);
1524 #else
1525 printf("Uzawa setupPrecon ERROR : mli not available.\n");
1526 exit(1);
1527 #endif
1528 break;
1529 }
1530 return 0;
1531 }
1532
1533 //***************************************************************************
1534 // setup solver
1535 //---------------------------------------------------------------------------
1536
setupSolver(HYPRE_Solver * solver,HYPRE_ParCSRMatrix Amat,HYPRE_ParVector fvec,HYPRE_ParVector xvec,HYPRE_Solver precon,HYPRE_Uzawa_PARAMS paramPtr)1537 int HYPRE_LSI_Uzawa::setupSolver(HYPRE_Solver *solver,HYPRE_ParCSRMatrix Amat,
1538 HYPRE_ParVector fvec, HYPRE_ParVector xvec,
1539 HYPRE_Solver precon, HYPRE_Uzawa_PARAMS paramPtr)
1540 {
1541 //------------------------------------------------------------------
1542 // create solver context
1543 //------------------------------------------------------------------
1544
1545 switch ( paramPtr.SolverID_ )
1546 {
1547 case 1 :
1548 HYPRE_ParCSRPCGCreate(mpiComm_, solver);
1549 HYPRE_ParCSRPCGSetMaxIter(*solver, paramPtr.MaxIter_ );
1550 HYPRE_ParCSRPCGSetTol(*solver, paramPtr.Tol_);
1551 HYPRE_ParCSRPCGSetLogging(*solver, outputLevel_);
1552 HYPRE_ParCSRPCGSetRelChange(*solver, 0);
1553 HYPRE_ParCSRPCGSetTwoNorm(*solver, 1);
1554 switch ( paramPtr.PrecondID_ )
1555 {
1556 case 1 :
1557 HYPRE_ParCSRPCGSetPrecond(*solver, HYPRE_ParCSRDiagScale,
1558 HYPRE_ParCSRDiagScaleSetup,precon);
1559 break;
1560 case 2 :
1561 HYPRE_ParCSRPCGSetPrecond(*solver,HYPRE_ParCSRParaSailsSolve,
1562 HYPRE_ParCSRParaSailsSetup,precon);
1563 break;
1564 case 3 :
1565 HYPRE_ParCSRPCGSetPrecond(*solver, HYPRE_BoomerAMGSolve,
1566 HYPRE_BoomerAMGSetup, precon);
1567 break;
1568 case 4 :
1569 HYPRE_ParCSRPCGSetPrecond(*solver, HYPRE_ParCSRPilutSolve,
1570 HYPRE_ParCSRPilutSetup, precon);
1571 break;
1572 case 5 :
1573 HYPRE_ParCSRPCGSetPrecond(*solver, HYPRE_EuclidSolve,
1574 HYPRE_EuclidSetup, precon);
1575 break;
1576 case 6 :
1577 #ifdef HAVE_MLI
1578 HYPRE_ParCSRPCGSetPrecond(*solver,HYPRE_LSI_MLISolve,
1579 HYPRE_LSI_MLISetup, precon);
1580 #else
1581 printf("Uzawa setupSolver ERROR : mli not available.\n");
1582 exit(1);
1583 #endif
1584 break;
1585 }
1586 HYPRE_ParCSRPCGSetup(*solver, Amat, fvec, xvec);
1587 break;
1588
1589 case 2 :
1590 HYPRE_ParCSRGMRESCreate(mpiComm_, solver);
1591 HYPRE_ParCSRGMRESSetMaxIter(*solver, paramPtr.MaxIter_ );
1592 HYPRE_ParCSRGMRESSetTol(*solver, paramPtr.Tol_);
1593 HYPRE_ParCSRGMRESSetLogging(*solver, outputLevel_);
1594 HYPRE_ParCSRGMRESSetKDim(*solver, 50);
1595 switch ( paramPtr.PrecondID_ )
1596 {
1597 case 1 :
1598 HYPRE_ParCSRGMRESSetPrecond(*solver, HYPRE_ParCSRDiagScale,
1599 HYPRE_ParCSRDiagScaleSetup,precon);
1600 break;
1601 case 2 :
1602 HYPRE_ParCSRGMRESSetPrecond(*solver,HYPRE_ParCSRParaSailsSolve,
1603 HYPRE_ParCSRParaSailsSetup,precon);
1604 break;
1605 case 3 :
1606 HYPRE_ParCSRGMRESSetPrecond(*solver, HYPRE_BoomerAMGSolve,
1607 HYPRE_BoomerAMGSetup, precon);
1608 break;
1609 case 4 :
1610 HYPRE_ParCSRGMRESSetPrecond(*solver, HYPRE_ParCSRPilutSolve,
1611 HYPRE_ParCSRPilutSetup, precon);
1612 break;
1613 case 5 :
1614 HYPRE_ParCSRGMRESSetPrecond(*solver, HYPRE_EuclidSolve,
1615 HYPRE_EuclidSetup, precon);
1616 break;
1617 case 6 :
1618 #ifdef HAVEL_MLI
1619 HYPRE_ParCSRGMRESSetPrecond(*solver,HYPRE_LSI_MLISolve,
1620 HYPRE_LSI_MLISetup, precon);
1621 #else
1622 printf("Uzawa setupSolver ERROR : mli not available.\n");
1623 exit(1);
1624 #endif
1625 break;
1626 }
1627 HYPRE_ParCSRGMRESSetup(*solver, Amat, fvec, xvec);
1628 break;
1629 }
1630 return 0;
1631 }
1632
1633