source: CIVL/examples/mpi-omp/AMG2013/parcsr_mv/par_csr_matvec.c@ e898536

1.23 2.0 acw/focus-triggers main test-branch
Last change on this file since e898536 was 2aa6644, checked in by Michael Rogers <mrogers08@…>, 11 years ago

Added AMG2013 (hybrid program of mpi/omp) example to mpi-omp directory

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@2362 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 17.8 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 * Matvec functions for hypre_CSRMatrix class.
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22#include <assert.h>
23
24/*--------------------------------------------------------------------------
25 * hypre_ParCSRMatrixMatvec
26 *--------------------------------------------------------------------------*/
27
28int
29hypre_ParCSRMatrixMatvec( double alpha,
30 hypre_ParCSRMatrix *A,
31 hypre_ParVector *x,
32 double beta,
33 hypre_ParVector *y )
34{
35 hypre_ParCSRCommHandle **comm_handle;
36 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
37 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
38 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
39 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
40 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
41 HYPRE_BigInt num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
42 HYPRE_BigInt num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
43
44 hypre_Vector *x_tmp;
45 HYPRE_BigInt x_size = hypre_ParVectorGlobalSize(x);
46 HYPRE_BigInt y_size = hypre_ParVectorGlobalSize(y);
47 int num_vectors = 1;
48 int num_cols_offd = hypre_CSRMatrixNumCols(offd);
49 int ierr = 0;
50 int num_sends, i, j, jv, index, start;
51
52 /*int vecstride = hypre_VectorVectorStride( x_local );
53 int idxstride = hypre_VectorIndexStride( x_local );*/
54
55 double *x_tmp_data, **x_buf_data;
56 double *x_local_data = hypre_VectorData(x_local);
57 /*---------------------------------------------------------------------
58 * Check for size compatibility. ParMatvec returns ierr = 11 if
59 * length of X doesn't equal the number of columns of A,
60 * ierr = 12 if the length of Y doesn't equal the number of rows
61 * of A, and ierr = 13 if both are true.
62 *
63 * Because temporary vectors are often used in ParMatvec, none of
64 * these conditions terminates processing, and the ierr flag
65 * is informational only.
66 *--------------------------------------------------------------------*/
67
68 /*hypre_assert( idxstride>0 );*/
69
70 if (num_cols != x_size)
71 ierr = 11;
72
73 if (num_rows != y_size)
74 ierr = 12;
75
76 if (num_cols != x_size && num_rows != y_size)
77 ierr = 13;
78
79 x_tmp = hypre_SeqVectorCreate( num_cols_offd );
80
81 hypre_SeqVectorInitialize(x_tmp);
82 x_tmp_data = hypre_VectorData(x_tmp);
83
84 comm_handle = hypre_CTAlloc(hypre_ParCSRCommHandle*,num_vectors);
85
86 /*---------------------------------------------------------------------
87 * If there exists no CommPkg for A, a CommPkg is generated using
88 * equally load balanced partitionings
89 *--------------------------------------------------------------------*/
90 if (!comm_pkg)
91 {
92#ifdef HYPRE_NO_GLOBAL_PARTITION
93 hypre_NewCommPkgCreate(A);
94#else
95 hypre_MatvecCommPkgCreate(A);
96#endif
97 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
98 }
99
100 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
101 x_buf_data = hypre_CTAlloc( double*, num_vectors );
102 for ( jv=0; jv<num_vectors; ++jv )
103 x_buf_data[jv] = hypre_CTAlloc(double, hypre_ParCSRCommPkgSendMapStart
104 (comm_pkg, num_sends));
105
106 /*if ( num_vectors==1 )*/
107 {
108 index = 0;
109 for (i = 0; i < num_sends; i++)
110 {
111 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
112 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
113 x_buf_data[0][index++]
114 = x_local_data[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
115 }
116 }
117 /*else
118 for ( jv=0; jv<num_vectors; ++jv )
119 {
120 index = 0;
121 for (i = 0; i < num_sends; i++)
122 {
123 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
124 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
125 x_buf_data[jv][index++]
126 = x_local_data[
127 jv*vecstride +
128 idxstride*hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j) ];
129 }
130 }
131
132 hypre_assert( idxstride==1 );*/
133 /* >>> ... The assert is because the following loop only works for 'column' storage of a multivector <<<
134 >>> This needs to be fixed to work more generally, at least for 'row' storage. <<<
135 >>> This in turn, means either change CommPkg so num_sends is no.zones*no.vectors (not no.zones)
136 >>> or, less dangerously, put a stride in the logic of CommHandleCreate (stride either from a
137 >>> new arg or a new variable inside CommPkg). Or put the num_vector iteration inside
138 >>> CommHandleCreate (perhaps a new multivector variant of it).
139 */
140 for ( jv=0; jv<num_vectors; ++jv )
141 {
142 comm_handle[jv] = hypre_ParCSRCommHandleCreate
143 ( 1, comm_pkg, x_buf_data[jv], &(x_tmp_data[jv*num_cols_offd]) );
144 }
145
146 hypre_CSRMatrixMatvec( alpha, diag, x_local, beta, y_local);
147
148 for ( jv=0; jv<num_vectors; ++jv )
149 {
150 hypre_ParCSRCommHandleDestroy(comm_handle[jv]);
151 comm_handle[jv] = NULL;
152 }
153 hypre_TFree(comm_handle);
154
155 if (num_cols_offd) hypre_CSRMatrixMatvec( alpha, offd, x_tmp, 1.0, y_local);
156
157 hypre_SeqVectorDestroy(x_tmp);
158 x_tmp = NULL;
159 for ( jv=0; jv<num_vectors; ++jv ) hypre_TFree(x_buf_data[jv]);
160 hypre_TFree(x_buf_data);
161
162 return ierr;
163}
164
165/*--------------------------------------------------------------------------
166 * hypre_ParCSRMatrixMatvecT
167 *
168 * Performs y <- alpha * A^T * x + beta * y
169 *
170 *--------------------------------------------------------------------------*/
171
172int
173hypre_ParCSRMatrixMatvecT( double alpha,
174 hypre_ParCSRMatrix *A,
175 hypre_ParVector *x,
176 double beta,
177 hypre_ParVector *y )
178{
179 hypre_ParCSRCommHandle **comm_handle;
180 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
181 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
182 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
183 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
184 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
185 hypre_Vector *y_tmp;
186 int vecstride = hypre_VectorVectorStride( y_local );
187 int idxstride = hypre_VectorIndexStride( y_local );
188 double *y_tmp_data, **y_buf_data;
189 double *y_local_data = hypre_VectorData(y_local);
190
191 HYPRE_BigInt num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
192 HYPRE_BigInt num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
193 int num_cols_offd = hypre_CSRMatrixNumCols(offd);
194 HYPRE_BigInt x_size = hypre_ParVectorGlobalSize(x);
195 HYPRE_BigInt y_size = hypre_ParVectorGlobalSize(y);
196 int num_vectors = hypre_VectorNumVectors(y_local);
197
198 int i, j, jv, index, start, num_sends;
199
200 int ierr = 0;
201
202 /*---------------------------------------------------------------------
203 * Check for size compatibility. MatvecT returns ierr = 1 if
204 * length of X doesn't equal the number of rows of A,
205 * ierr = 2 if the length of Y doesn't equal the number of
206 * columns of A, and ierr = 3 if both are true.
207 *
208 * Because temporary vectors are often used in MatvecT, none of
209 * these conditions terminates processing, and the ierr flag
210 * is informational only.
211 *--------------------------------------------------------------------*/
212
213 if (num_rows != x_size)
214 ierr = 1;
215
216 if (num_cols != y_size)
217 ierr = 2;
218
219 if (num_rows != x_size && num_cols != y_size)
220 ierr = 3;
221 /*-----------------------------------------------------------------------
222 *-----------------------------------------------------------------------*/
223
224 comm_handle = hypre_CTAlloc(hypre_ParCSRCommHandle*,num_vectors);
225
226 /*if ( num_vectors==1 )
227 {*/
228 y_tmp = hypre_SeqVectorCreate(num_cols_offd);
229 /*}
230 else
231 {
232 y_tmp = hypre_SeqMultiVectorCreate(num_cols_offd,num_vectors);
233 }*/
234 hypre_SeqVectorInitialize(y_tmp);
235
236 /*---------------------------------------------------------------------
237 * If there exists no CommPkg for A, a CommPkg is generated using
238 * equally load balanced partitionings
239 *--------------------------------------------------------------------*/
240 if (!comm_pkg)
241 {
242#ifdef HYPRE_NO_GLOBAL_PARTITION
243 hypre_NewCommPkgCreate(A);
244#else
245 hypre_MatvecCommPkgCreate(A);
246#endif
247 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
248 }
249
250 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
251 y_buf_data = hypre_CTAlloc( double*, num_vectors );
252 for ( jv=0; jv<num_vectors; ++jv )
253 y_buf_data[jv] = hypre_CTAlloc(double, hypre_ParCSRCommPkgSendMapStart
254 (comm_pkg, num_sends));
255 y_tmp_data = hypre_VectorData(y_tmp);
256 y_local_data = hypre_VectorData(y_local);
257
258 hypre_assert( idxstride==1 ); /* >>> only 'column' storage of multivectors implemented so far */
259
260 if (num_cols_offd) hypre_CSRMatrixMatvecT(alpha, offd, x_local, 0.0, y_tmp);
261
262 for ( jv=0; jv<num_vectors; ++jv )
263 {
264 /* >>> this is where we assume multivectors are 'column' storage */
265 comm_handle[jv] = hypre_ParCSRCommHandleCreate
266 ( 2, comm_pkg, &(y_tmp_data[jv*num_cols_offd]), y_buf_data[jv] );
267 }
268
269 hypre_CSRMatrixMatvecT(alpha, diag, x_local, beta, y_local);
270
271 for ( jv=0; jv<num_vectors; ++jv )
272 {
273 hypre_ParCSRCommHandleDestroy(comm_handle[jv]);
274 comm_handle[jv] = NULL;
275 }
276 hypre_TFree(comm_handle);
277
278 if ( num_vectors==1 )
279 {
280 index = 0;
281 for (i = 0; i < num_sends; i++)
282 {
283 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
284 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
285 y_local_data[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)]
286 += y_buf_data[0][index++];
287 }
288 }
289 else
290 for ( jv=0; jv<num_vectors; ++jv )
291 {
292 index = 0;
293 for (i = 0; i < num_sends; i++)
294 {
295 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
296 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
297 y_local_data[ jv*vecstride +
298 idxstride*hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j) ]
299 += y_buf_data[jv][index++];
300 }
301 }
302
303 hypre_SeqVectorDestroy(y_tmp);
304 y_tmp = NULL;
305 for ( jv=0; jv<num_vectors; ++jv ) hypre_TFree(y_buf_data[jv]);
306 hypre_TFree(y_buf_data);
307
308 return ierr;
309}
310/*--------------------------------------------------------------------------
311 * hypre_ParCSRMatrixMatvec_FF
312 *--------------------------------------------------------------------------*/
313
314int
315hypre_ParCSRMatrixMatvec_FF( double alpha,
316 hypre_ParCSRMatrix *A,
317 hypre_ParVector *x,
318 double beta,
319 hypre_ParVector *y,
320 int *CF_marker,
321 int fpt )
322{
323 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
324 hypre_ParCSRCommHandle *comm_handle;
325 hypre_ParCSRCommPkg *comm_pkg = hypre_ParCSRMatrixCommPkg(A);
326 hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
327 hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
328 hypre_Vector *x_local = hypre_ParVectorLocalVector(x);
329 hypre_Vector *y_local = hypre_ParVectorLocalVector(y);
330 HYPRE_BigInt num_rows = hypre_ParCSRMatrixGlobalNumRows(A);
331 HYPRE_BigInt num_cols = hypre_ParCSRMatrixGlobalNumCols(A);
332
333 hypre_Vector *x_tmp;
334 HYPRE_BigInt x_size = hypre_ParVectorGlobalSize(x);
335 HYPRE_BigInt y_size = hypre_ParVectorGlobalSize(y);
336 HYPRE_BigInt num_cols_offd = hypre_CSRMatrixNumCols(offd);
337 int ierr = 0;
338 int num_sends, i, j, index, start, num_procs;
339 int *int_buf_data = NULL;
340 int *CF_marker_offd = NULL;
341
342
343 double *x_tmp_data = NULL;
344 double *x_buf_data = NULL;
345 double *x_local_data = hypre_VectorData(x_local);
346 /*---------------------------------------------------------------------
347 * Check for size compatibility. ParMatvec returns ierr = 11 if
348 * length of X doesn't equal the number of columns of A,
349 * ierr = 12 if the length of Y doesn't equal the number of rows
350 * of A, and ierr = 13 if both are true.
351 *
352 * Because temporary vectors are often used in ParMatvec, none of
353 * these conditions terminates processing, and the ierr flag
354 * is informational only.
355 *--------------------------------------------------------------------*/
356
357 MPI_Comm_size(comm,&num_procs);
358
359 if (num_cols != x_size)
360 ierr = 11;
361
362 if (num_rows != y_size)
363 ierr = 12;
364
365 if (num_cols != x_size && num_rows != y_size)
366 ierr = 13;
367
368 if (num_procs > 1)
369 {
370 if (num_cols_offd)
371 {
372 x_tmp = hypre_SeqVectorCreate( num_cols_offd );
373 hypre_SeqVectorInitialize(x_tmp);
374 x_tmp_data = hypre_VectorData(x_tmp);
375 }
376
377 /*---------------------------------------------------------------------
378 * If there exists no CommPkg for A, a CommPkg is generated using
379 * equally load balanced partitionings
380 *--------------------------------------------------------------------*/
381 if (!comm_pkg)
382 {
383#ifdef HYPRE_NO_GLOBAL_PARTITION
384 hypre_NewCommPkgCreate(A);
385#else
386 hypre_MatvecCommPkgCreate(A);
387#endif
388 comm_pkg = hypre_ParCSRMatrixCommPkg(A);
389 }
390
391 num_sends = hypre_ParCSRCommPkgNumSends(comm_pkg);
392 if (num_sends)
393 x_buf_data = hypre_CTAlloc(double, hypre_ParCSRCommPkgSendMapStart
394 (comm_pkg, num_sends));
395
396 index = 0;
397 for (i = 0; i < num_sends; i++)
398 {
399 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
400 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
401 x_buf_data[index++]
402 = x_local_data[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
403 }
404 comm_handle = hypre_ParCSRCommHandleCreate ( 1, comm_pkg, x_buf_data, x_tmp_data );
405 }
406 hypre_CSRMatrixMatvec_FF( alpha, diag, x_local, beta, y_local, CF_marker, CF_marker, fpt);
407
408 if (num_procs > 1)
409 {
410 hypre_ParCSRCommHandleDestroy(comm_handle);
411 comm_handle = NULL;
412
413 if (num_sends)
414 int_buf_data = hypre_CTAlloc(int, hypre_ParCSRCommPkgSendMapStart
415 (comm_pkg, num_sends));
416 if (num_cols_offd) CF_marker_offd = hypre_CTAlloc(int, num_cols_offd);
417 index = 0;
418 for (i = 0; i < num_sends; i++)
419 {
420 start = hypre_ParCSRCommPkgSendMapStart(comm_pkg, i);
421 for (j = start; j < hypre_ParCSRCommPkgSendMapStart(comm_pkg, i+1); j++)
422 int_buf_data[index++]
423 = CF_marker[hypre_ParCSRCommPkgSendMapElmt(comm_pkg,j)];
424 }
425 comm_handle = hypre_ParCSRCommHandleCreate(11,comm_pkg,int_buf_data,CF_marker_offd );
426
427 hypre_ParCSRCommHandleDestroy(comm_handle);
428 comm_handle = NULL;
429
430 if (num_cols_offd) hypre_CSRMatrixMatvec_FF( alpha, offd, x_tmp, 1.0, y_local,
431 CF_marker, CF_marker_offd, fpt);
432
433 hypre_SeqVectorDestroy(x_tmp);
434 x_tmp = NULL;
435 hypre_TFree(x_buf_data);
436 hypre_TFree(int_buf_data);
437 hypre_TFree(CF_marker_offd);
438 }
439
440 return ierr;
441}
Note: See TracBrowser for help on using the repository browser.