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

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