source: CIVL/examples/mpi-omp/AMG2013/sstruct_mv/sstruct_matvec.c@ 7d77e64

main test-branch
Last change on this file since 7d77e64 was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 12.5 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15/******************************************************************************
16 *
17 * SStruct matrix-vector multiply routine
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23/*==========================================================================
24 * PMatvec routines
25 *==========================================================================*/
26
27/*--------------------------------------------------------------------------
28 * hypre_SStructPMatvecData data structure
29 *--------------------------------------------------------------------------*/
30
31typedef struct
32{
33 int nvars;
34 void ***smatvec_data;
35
36} hypre_SStructPMatvecData;
37
38/*--------------------------------------------------------------------------
39 * hypre_SStructPMatvecCreate
40 *--------------------------------------------------------------------------*/
41
42int
43hypre_SStructPMatvecCreate( void **pmatvec_vdata_ptr )
44{
45 int ierr = 0;
46 hypre_SStructPMatvecData *pmatvec_data;
47
48 pmatvec_data = hypre_CTAlloc(hypre_SStructPMatvecData, 1);
49 *pmatvec_vdata_ptr = (void *) pmatvec_data;
50
51 return ierr;
52}
53
54/*--------------------------------------------------------------------------
55 * hypre_SStructPMatvecSetup
56 *--------------------------------------------------------------------------*/
57
58int
59hypre_SStructPMatvecSetup( void *pmatvec_vdata,
60 hypre_SStructPMatrix *pA,
61 hypre_SStructPVector *px )
62{
63 int ierr = 0;
64 hypre_SStructPMatvecData *pmatvec_data = pmatvec_vdata;
65 int nvars;
66 void ***smatvec_data;
67 hypre_StructMatrix *sA;
68 hypre_StructVector *sx;
69 int vi, vj;
70
71 nvars = hypre_SStructPMatrixNVars(pA);
72 smatvec_data = hypre_TAlloc(void **, nvars);
73 for (vi = 0; vi < nvars; vi++)
74 {
75 smatvec_data[vi] = hypre_TAlloc(void *, nvars);
76 for (vj = 0; vj < nvars; vj++)
77 {
78 sA = hypre_SStructPMatrixSMatrix(pA, vi, vj);
79 sx = hypre_SStructPVectorSVector(px, vj);
80 smatvec_data[vi][vj] = NULL;
81 if (sA != NULL)
82 {
83 smatvec_data[vi][vj] = hypre_StructMatvecCreate();
84 hypre_StructMatvecSetup(smatvec_data[vi][vj], sA, sx);
85 }
86 }
87 }
88 (pmatvec_data -> nvars) = nvars;
89 (pmatvec_data -> smatvec_data) = smatvec_data;
90
91 return ierr;
92}
93
94/*--------------------------------------------------------------------------
95 * hypre_SStructPMatvecCompute
96 *--------------------------------------------------------------------------*/
97
98int
99hypre_SStructPMatvecCompute( void *pmatvec_vdata,
100 double alpha,
101 hypre_SStructPMatrix *pA,
102 hypre_SStructPVector *px,
103 double beta,
104 hypre_SStructPVector *py )
105{
106 int ierr = 0;
107
108 hypre_SStructPMatvecData *pmatvec_data = pmatvec_vdata;
109 int nvars = (pmatvec_data -> nvars);
110 void ***smatvec_data = (pmatvec_data -> smatvec_data);
111
112 void *sdata;
113 hypre_StructMatrix *sA;
114 hypre_StructVector *sx;
115 hypre_StructVector *sy;
116
117 int vi, vj;
118
119 for (vi = 0; vi < nvars; vi++)
120 {
121 sy = hypre_SStructPVectorSVector(py, vi);
122
123 /* diagonal block computation */
124 if (smatvec_data[vi][vi] != NULL)
125 {
126 sdata = smatvec_data[vi][vi];
127 sA = hypre_SStructPMatrixSMatrix(pA, vi, vi);
128 sx = hypre_SStructPVectorSVector(px, vi);
129 hypre_StructMatvecCompute(sdata, alpha, sA, sx, beta, sy);
130 }
131 else
132 {
133 hypre_StructScale(beta, sy);
134 }
135
136 /* off-diagonal block computation */
137 for (vj = 0; vj < nvars; vj++)
138 {
139 if ((smatvec_data[vi][vj] != NULL) && (vj != vi))
140 {
141 sdata = smatvec_data[vi][vj];
142 sA = hypre_SStructPMatrixSMatrix(pA, vi, vj);
143 sx = hypre_SStructPVectorSVector(px, vj);
144 hypre_StructMatvecCompute(sdata, alpha, sA, sx, 1.0, sy);
145 }
146 }
147 }
148
149 return ierr;
150}
151
152/*--------------------------------------------------------------------------
153 * hypre_SStructPMatvecDestroy
154 *--------------------------------------------------------------------------*/
155
156int
157hypre_SStructPMatvecDestroy( void *pmatvec_vdata )
158{
159 int ierr = 0;
160 hypre_SStructPMatvecData *pmatvec_data = pmatvec_vdata;
161 int nvars;
162 void ***smatvec_data;
163 int vi, vj;
164
165 if (pmatvec_data)
166 {
167 nvars = (pmatvec_data -> nvars);
168 smatvec_data = (pmatvec_data -> smatvec_data);
169 for (vi = 0; vi < nvars; vi++)
170 {
171 for (vj = 0; vj < nvars; vj++)
172 {
173 if (smatvec_data[vi][vj] != NULL)
174 {
175 hypre_StructMatvecDestroy(smatvec_data[vi][vj]);
176 }
177 }
178 hypre_TFree(smatvec_data[vi]);
179 }
180 hypre_TFree(smatvec_data);
181 hypre_TFree(pmatvec_data);
182 }
183
184 return ierr;
185}
186
187/*--------------------------------------------------------------------------
188 * hypre_SStructPMatvec
189 *--------------------------------------------------------------------------*/
190
191int
192hypre_SStructPMatvec( double alpha,
193 hypre_SStructPMatrix *pA,
194 hypre_SStructPVector *px,
195 double beta,
196 hypre_SStructPVector *py )
197{
198 int ierr = 0;
199
200 void *pmatvec_data;
201
202 hypre_SStructPMatvecCreate(&pmatvec_data);
203 ierr = hypre_SStructPMatvecSetup(pmatvec_data, pA, px);
204 ierr = hypre_SStructPMatvecCompute(pmatvec_data, alpha, pA, px, beta, py);
205 ierr = hypre_SStructPMatvecDestroy(pmatvec_data);
206
207 return ierr;
208}
209
210/*==========================================================================
211 * Matvec routines
212 *==========================================================================*/
213
214/*--------------------------------------------------------------------------
215 * hypre_SStructMatvecData data structure
216 *--------------------------------------------------------------------------*/
217
218typedef struct
219{
220 int nparts;
221 void **pmatvec_data;
222
223} hypre_SStructMatvecData;
224
225/*--------------------------------------------------------------------------
226 * hypre_SStructMatvecCreate
227 *--------------------------------------------------------------------------*/
228
229int
230hypre_SStructMatvecCreate( void **matvec_vdata_ptr )
231{
232 int ierr = 0;
233 hypre_SStructMatvecData *matvec_data;
234
235 matvec_data = hypre_CTAlloc(hypre_SStructMatvecData, 1);
236 *matvec_vdata_ptr = (void *) matvec_data;
237
238 return ierr;
239}
240
241/*--------------------------------------------------------------------------
242 * hypre_SStructMatvecSetup
243 *--------------------------------------------------------------------------*/
244
245int
246hypre_SStructMatvecSetup( void *matvec_vdata,
247 hypre_SStructMatrix *A,
248 hypre_SStructVector *x )
249{
250 int ierr = 0;
251 hypre_SStructMatvecData *matvec_data = matvec_vdata;
252 int nparts;
253 void **pmatvec_data;
254 hypre_SStructPMatrix *pA;
255 hypre_SStructPVector *px;
256 int part;
257
258 nparts = hypre_SStructMatrixNParts(A);
259 pmatvec_data = hypre_TAlloc(void *, nparts);
260 for (part = 0; part < nparts; part++)
261 {
262 hypre_SStructPMatvecCreate(&pmatvec_data[part]);
263 pA = hypre_SStructMatrixPMatrix(A, part);
264 px = hypre_SStructVectorPVector(x, part);
265 hypre_SStructPMatvecSetup(pmatvec_data[part], pA, px);
266 }
267 (matvec_data -> nparts) = nparts;
268 (matvec_data -> pmatvec_data) = pmatvec_data;
269
270 return ierr;
271}
272
273/*--------------------------------------------------------------------------
274 * hypre_SStructMatvecCompute
275 *--------------------------------------------------------------------------*/
276
277int
278hypre_SStructMatvecCompute( void *matvec_vdata,
279 double alpha,
280 hypre_SStructMatrix *A,
281 hypre_SStructVector *x,
282 double beta,
283 hypre_SStructVector *y )
284{
285 int ierr = 0;
286
287 hypre_SStructMatvecData *matvec_data = matvec_vdata;
288 int nparts = (matvec_data -> nparts);
289 void **pmatvec_data = (matvec_data -> pmatvec_data);
290
291 void *pdata;
292 hypre_SStructPMatrix *pA;
293 hypre_SStructPVector *px;
294 hypre_SStructPVector *py;
295
296 hypre_ParCSRMatrix *parcsrA = hypre_SStructMatrixParCSRMatrix(A);
297 hypre_ParVector *parx;
298 hypre_ParVector *pary;
299
300 int part;
301 int x_object_type= hypre_SStructVectorObjectType(x);
302 int A_object_type= hypre_SStructMatrixObjectType(A);
303
304 if (x_object_type != A_object_type)
305 {
306 printf("possible error: A and x are different object types\n");
307 }
308
309 if (x_object_type == HYPRE_SSTRUCT)
310 {
311 /* do S-matrix computations */
312 for (part = 0; part < nparts; part++)
313 {
314 pdata = pmatvec_data[part];
315 pA = hypre_SStructMatrixPMatrix(A, part);
316 px = hypre_SStructVectorPVector(x, part);
317 py = hypre_SStructVectorPVector(y, part);
318 hypre_SStructPMatvecCompute(pdata, alpha, pA, px, beta, py);
319 }
320
321 /* do U-matrix computations */
322
323 /* GEC1002 the data chunk pointed by the local-parvectors
324 * inside the semistruct vectors x and y is now identical to the
325 * data chunk of the structure vectors x and y. The role of the function
326 * convert is to pass the addresses of the data chunk
327 * to the parx and pary. */
328
329 hypre_SStructVectorConvert(x, &parx);
330 hypre_SStructVectorConvert(y, &pary);
331
332 hypre_ParCSRMatrixMatvec(alpha, parcsrA, parx, 1.0, pary);
333
334 /* dummy functions since there is nothing to restore */
335
336 hypre_SStructVectorRestore(x, NULL);
337 hypre_SStructVectorRestore(y, pary);
338
339 parx = NULL;
340
341 }
342
343 else if (x_object_type == HYPRE_PARCSR)
344 {
345 hypre_SStructVectorConvert(x, &parx);
346 hypre_SStructVectorConvert(y, &pary);
347
348 hypre_ParCSRMatrixMatvec(alpha, parcsrA, parx, beta, pary);
349
350 hypre_SStructVectorRestore(x, NULL);
351 hypre_SStructVectorRestore(y, pary);
352
353 parx = NULL;
354 }
355
356 return ierr;
357}
358
359/*--------------------------------------------------------------------------
360 * hypre_SStructMatvecDestroy
361 *--------------------------------------------------------------------------*/
362
363int
364hypre_SStructMatvecDestroy( void *matvec_vdata )
365{
366 int ierr = 0;
367 hypre_SStructMatvecData *matvec_data = matvec_vdata;
368 int nparts;
369 void **pmatvec_data;
370 int part;
371
372 if (matvec_data)
373 {
374 nparts = (matvec_data -> nparts);
375 pmatvec_data = (matvec_data -> pmatvec_data);
376 for (part = 0; part < nparts; part++)
377 {
378 hypre_SStructPMatvecDestroy(pmatvec_data[part]);
379 }
380 hypre_TFree(pmatvec_data);
381 hypre_TFree(matvec_data);
382 }
383
384 return ierr;
385}
386
387/*--------------------------------------------------------------------------
388 * hypre_SStructMatvec
389 *--------------------------------------------------------------------------*/
390
391int
392hypre_SStructMatvec( double alpha,
393 hypre_SStructMatrix *A,
394 hypre_SStructVector *x,
395 double beta,
396 hypre_SStructVector *y )
397{
398 int ierr = 0;
399
400 void *matvec_data;
401
402 hypre_SStructMatvecCreate(&matvec_data);
403 ierr = hypre_SStructMatvecSetup(matvec_data, A, x);
404 ierr = hypre_SStructMatvecCompute(matvec_data, alpha, A, x, beta, y);
405 ierr = hypre_SStructMatvecDestroy(matvec_data);
406
407 return ierr;
408}
Note: See TracBrowser for help on using the repository browser.