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_, &parasails);
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