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 #include <math.h>
9 #include <string.h>
10 #include "_hypre_parcsr_mv.h"
11 #include "_hypre_parcsr_ls.h"
12 #include "par_amg.h"
13 #include "mli_solver_hsgs.h"
14
15 /******************************************************************************
16 * symmetric Gauss-Seidel relaxation scheme in BoomerAMG
17 *****************************************************************************/
18
19 /******************************************************************************
20 * constructor
21 *---------------------------------------------------------------------------*/
22
MLI_Solver_HSGS(char * name)23 MLI_Solver_HSGS::MLI_Solver_HSGS(char *name) : MLI_Solver(name)
24 {
25 Amat_ = NULL;
26 nSweeps_ = 1;
27 relaxWeights_ = 1.0;
28 relaxOmega_ = 1.0;
29 mliVec_ = NULL;
30 calcOmega_ = 1;
31 }
32
33 /******************************************************************************
34 * destructor
35 *---------------------------------------------------------------------------*/
36
~MLI_Solver_HSGS()37 MLI_Solver_HSGS::~MLI_Solver_HSGS()
38 {
39 if (mliVec_ != NULL) delete mliVec_;
40 mliVec_ = NULL;
41 }
42
43 /******************************************************************************
44 * set up the smoother
45 *---------------------------------------------------------------------------*/
46
setup(MLI_Matrix * mat)47 int MLI_Solver_HSGS::setup(MLI_Matrix *mat)
48 {
49 Amat_ = mat;
50 if (mliVec_ != NULL) delete mliVec_;
51 mliVec_ = Amat_->createVector();
52 if (calcOmega_ == 1) calcOmega();
53 return 0;
54 }
55
56 /******************************************************************************
57 * apply function
58 *---------------------------------------------------------------------------*/
59
solve(MLI_Vector * fIn,MLI_Vector * uIn)60 int MLI_Solver_HSGS::solve(MLI_Vector *fIn, MLI_Vector *uIn)
61 {
62 int relaxType=6, relaxPts=0, iS;
63 hypre_ParCSRMatrix *A;
64 hypre_ParVector *f, *u, *vTemp;
65 hypre_ParVector *zTemp = NULL;
66
67 //int mypid;
68 //double rnorm;
69 //MPI_Comm comm;
70
71 /*-----------------------------------------------------------------
72 * fetch machine and smoother parameters
73 *-----------------------------------------------------------------*/
74
75 A = (hypre_ParCSRMatrix *) Amat_->getMatrix();
76 u = (hypre_ParVector *) uIn->getVector();
77 f = (hypre_ParVector *) fIn->getVector();
78 vTemp = (hypre_ParVector *) mliVec_->getVector();
79
80 /* AB: need an extra vector for threading */
81 if (hypre_NumThreads() > 1)
82 {
83 zTemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A),
84 hypre_ParCSRMatrixGlobalNumRows(A),
85 hypre_ParCSRMatrixRowStarts(A));
86 hypre_ParVectorInitialize(zTemp);
87 }
88
89 //comm = hypre_ParCSRMatrixComm(A);
90 //MPI_Comm_rank(comm, &mypid);
91 for (iS = 0; iS < nSweeps_; iS++)
92 {
93 hypre_BoomerAMGRelax(A,f,NULL,relaxType,relaxPts,relaxWeights_,
94 relaxOmega_,NULL,u,vTemp, zTemp);
95 //hypre_ParVectorCopy( f, vTemp );
96 //hypre_ParCSRMatrixMatvec( -1.0, A, u, 1.0, vTemp );
97 //rnorm = sqrt(hypre_ParVectorInnerProd( vTemp, vTemp ));
98 //if ( mypid == 0 )
99 // printf("\tMLI_Solver_HSGS iter = %4d, rnorm = %e (omega=%e)\n",
100 // iS, rnorm, relaxWeights_);
101 }
102
103 if (zTemp)
104 hypre_ParVectorDestroy(zTemp);
105
106 return 0;
107 }
108
109 /******************************************************************************
110 * set SGS parameters
111 *---------------------------------------------------------------------------*/
112
setParams(char * paramString,int argc,char ** argv)113 int MLI_Solver_HSGS::setParams(char *paramString, int argc, char **argv)
114 {
115 double *weights=NULL;
116 char param1[100];
117
118 sscanf(paramString, "%s", param1);
119 if (!strcmp(param1, "numSweeps"))
120 {
121 if ( argc != 1 )
122 {
123 printf("MLI_Solver_HSGS::setParams ERROR : needs 1 arg.\n");
124 return 1;
125 }
126 nSweeps_ = *(int*) argv[0];
127 if ( nSweeps_ < 1 ) nSweeps_ = 1;
128 return 0;
129 }
130 else if ( !strcmp(param1, "relaxWeight") )
131 {
132 if ( argc != 2 && argc != 1 )
133 {
134 printf("MLI_Solver_HSGS::setParams ERROR : needs 1 or 2 args.\n");
135 return 1;
136 }
137 if ( argc >= 1 ) nSweeps_ = *(int*) argv[0];
138 if ( argc == 2 ) weights = (double*) argv[1];
139 if ( nSweeps_ < 1 ) nSweeps_ = 1;
140 if ( weights != NULL ) relaxWeights_ = weights[0];
141 }
142 else if ( !strcmp(param1, "calcOmega") )
143 {
144 calcOmega_ = 1;
145 }
146 else return 1;
147 return 0;
148 }
149
150 /******************************************************************************
151 * calculate relax weight
152 *---------------------------------------------------------------------------*/
153
calcOmega()154 int MLI_Solver_HSGS::calcOmega()
155 {
156 int relaxType=6, relaxTypes[2], level=0, numCGSweeps=10;
157 hypre_ParCSRMatrix *A;
158 hypre_ParVector *vTemp;
159 hypre_ParAMGData *amgData;
160
161 A = (hypre_ParCSRMatrix *) Amat_->getMatrix();
162 amgData = (hypre_ParAMGData *) hypre_BoomerAMGCreate();
163 amgData->CF_marker_array = new int*[1];
164 amgData->CF_marker_array[0] = NULL;
165 amgData->A_array = new hypre_ParCSRMatrix*[1];
166 amgData->A_array[0] = A;
167 vTemp = (hypre_ParVector *) mliVec_->getVector();
168 amgData->Vtemp = vTemp;
169 relaxTypes[0] = 0;
170 relaxTypes[1] = relaxType;
171 amgData->grid_relax_type = relaxTypes;
172 amgData->smooth_num_levels = 0;
173 amgData->smooth_type = 0;
174 hypre_BoomerAMGCGRelaxWt((void *) amgData,level,numCGSweeps,&relaxOmega_);
175 //printf("HYPRE/FEI/MLI HSGS : relaxOmega = %e\n", relaxOmega_);
176 delete [] amgData->A_array;
177 delete [] amgData->CF_marker_array;
178 hypre_TFree(amgData, HYPRE_MEMORY_HOST);
179 return 0;
180 }
181