source: CIVL/examples/mpi-omp/AMG2013/struct_mv/struct_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: 58.1 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 * Structured matrix-vector multiply routine
17 *
18 *****************************************************************************/
19
20#include "headers.h"
21
22/* this currently cannot be greater than 7 */
23#ifdef MAX_DEPTH
24#undef MAX_DEPTH
25#endif
26#define MAX_DEPTH 7
27
28/*--------------------------------------------------------------------------
29 * hypre_StructMatvecData data structure
30 *--------------------------------------------------------------------------*/
31
32typedef struct
33{
34 hypre_StructMatrix *A;
35 hypre_StructVector *x;
36 hypre_ComputePkg *compute_pkg;
37
38} hypre_StructMatvecData;
39
40/*--------------------------------------------------------------------------
41 * hypre_StructMatvecCreate
42 *--------------------------------------------------------------------------*/
43
44void *
45hypre_StructMatvecCreate( )
46{
47 hypre_StructMatvecData *matvec_data;
48
49 matvec_data = hypre_CTAlloc(hypre_StructMatvecData, 1);
50
51 return (void *) matvec_data;
52}
53
54/*--------------------------------------------------------------------------
55 * hypre_StructMatvecSetup
56 *--------------------------------------------------------------------------*/
57
58int
59hypre_StructMatvecSetup( void *matvec_vdata,
60 hypre_StructMatrix *A,
61 hypre_StructVector *x )
62{
63 int ierr = 0;
64
65 hypre_StructMatvecData *matvec_data = matvec_vdata;
66
67 hypre_StructGrid *grid;
68 hypre_StructStencil *stencil;
69 hypre_ComputeInfo *compute_info;
70 hypre_ComputePkg *compute_pkg;
71
72 /*----------------------------------------------------------
73 * Set up the compute package
74 *----------------------------------------------------------*/
75
76 grid = hypre_StructMatrixGrid(A);
77 stencil = hypre_StructMatrixStencil(A);
78
79 hypre_CreateComputeInfo(grid, stencil, &compute_info);
80 hypre_ComputePkgCreate(compute_info, hypre_StructVectorDataSpace(x), 1,
81 grid, &compute_pkg);
82
83 /*----------------------------------------------------------
84 * Set up the matvec data structure
85 *----------------------------------------------------------*/
86
87 (matvec_data -> A) = hypre_StructMatrixRef(A);
88 (matvec_data -> x) = hypre_StructVectorRef(x);
89 (matvec_data -> compute_pkg) = compute_pkg;
90
91 return ierr;
92}
93
94/*--------------------------------------------------------------------------
95 * hypre_StructMatvecCompute
96 *--------------------------------------------------------------------------*/
97
98int
99hypre_StructMatvecCompute( void *matvec_vdata,
100 double alpha,
101 hypre_StructMatrix *A,
102 hypre_StructVector *x,
103 double beta,
104 hypre_StructVector *y )
105{
106 int ierr = 0;
107
108 hypre_StructMatvecData *matvec_data = matvec_vdata;
109
110 hypre_ComputePkg *compute_pkg;
111
112 hypre_CommHandle *comm_handle;
113
114 hypre_BoxArrayArray *compute_box_aa;
115 hypre_Box *y_data_box;
116
117 int yi;
118
119 double *xp;
120 double *yp;
121
122 hypre_BoxArray *boxes;
123 hypre_Box *box;
124 hypre_Index loop_size;
125 hypre_IndexRef start;
126 hypre_IndexRef stride;
127
128 int constant_coefficient;
129
130 double temp;
131 int compute_i, i;
132 int loopi, loopj, loopk;
133
134 /*-----------------------------------------------------------------------
135 * Initialize some things
136 *-----------------------------------------------------------------------*/
137
138 constant_coefficient = hypre_StructMatrixConstantCoefficient(A);
139 if (constant_coefficient==1) hypre_StructVectorClearBoundGhostValues( x );
140
141 compute_pkg = (matvec_data -> compute_pkg);
142
143 stride = hypre_ComputePkgStride(compute_pkg);
144
145 /*-----------------------------------------------------------------------
146 * Do (alpha == 0.0) computation
147 *-----------------------------------------------------------------------*/
148
149 if (alpha == 0.0)
150 {
151 boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(A));
152 hypre_ForBoxI(i, boxes)
153 {
154 box = hypre_BoxArrayBox(boxes, i);
155 start = hypre_BoxIMin(box);
156
157 y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
158 yp = hypre_StructVectorBoxData(y, i);
159
160 hypre_BoxGetSize(box, loop_size);
161
162 hypre_BoxLoop1Begin(loop_size,
163 y_data_box, start, stride, yi);
164
165 hypre_BoxLoop1For(loopi, loopj, loopk, yi)
166 {
167 yp[yi] *= beta;
168 }
169 hypre_BoxLoop1End(yi);
170 }
171
172 return ierr;
173 }
174
175 /*-----------------------------------------------------------------------
176 * Do (alpha != 0.0) computation
177 *-----------------------------------------------------------------------*/
178
179 for (compute_i = 0; compute_i < 2; compute_i++)
180 {
181 switch(compute_i)
182 {
183 case 0:
184 {
185 xp = hypre_StructVectorData(x);
186 hypre_InitializeIndtComputations(compute_pkg, xp, &comm_handle);
187 compute_box_aa = hypre_ComputePkgIndtBoxes(compute_pkg);
188
189 /*--------------------------------------------------------------
190 * initialize y= (beta/alpha)*y normally (where everything
191 * is multiplied by alpha at the end),
192 * beta*y for constant coefficient (where only Ax gets multiplied by alpha)
193 *--------------------------------------------------------------*/
194
195 if ( constant_coefficient==1 )
196 {
197 temp = beta;
198 }
199 else
200 {
201 temp = beta / alpha;
202 }
203 if (temp != 1.0)
204 {
205 boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(A));
206 hypre_ForBoxI(i, boxes)
207 {
208 box = hypre_BoxArrayBox(boxes, i);
209 start = hypre_BoxIMin(box);
210
211 y_data_box =
212 hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
213 yp = hypre_StructVectorBoxData(y, i);
214
215 if (temp == 0.0)
216 {
217 hypre_BoxGetSize(box, loop_size);
218
219 hypre_BoxLoop1Begin(loop_size,
220 y_data_box, start, stride, yi);
221
222 hypre_BoxLoop1For(loopi, loopj, loopk, yi)
223 {
224 yp[yi] = 0.0;
225 }
226 hypre_BoxLoop1End(yi);
227 }
228 else
229 {
230 hypre_BoxGetSize(box, loop_size);
231
232 hypre_BoxLoop1Begin(loop_size,
233 y_data_box, start, stride, yi);
234
235 hypre_BoxLoop1For(loopi, loopj, loopk, yi)
236 {
237 yp[yi] *= temp;
238 }
239 hypre_BoxLoop1End(yi);
240 }
241 }
242 }
243 }
244 break;
245
246 case 1:
247 {
248 hypre_FinalizeIndtComputations(comm_handle);
249 compute_box_aa = hypre_ComputePkgDeptBoxes(compute_pkg);
250 }
251 break;
252 }
253
254 /*--------------------------------------------------------------------
255 * y += A*x
256 *--------------------------------------------------------------------*/
257
258 switch( constant_coefficient )
259 {
260 case 0:
261 {
262 ierr += hypre_StructMatvecCC0( alpha, A, x, y, compute_box_aa, stride );
263 break;
264 }
265 case 1:
266 {
267 ierr += hypre_StructMatvecCC1( alpha, A, x, y, compute_box_aa, stride );
268 break;
269 }
270 case 2:
271 {
272 ierr += hypre_StructMatvecCC2( alpha, A, x, y, compute_box_aa, stride );
273 break;
274 }
275 }
276
277 }
278
279 return ierr;
280}
281
282/*--------------------------------------------------------------------------
283 * hypre_StructMatvecCC0
284 * core of struct matvec computation, for the case constant_coefficient==0
285 * (all coefficients are variable)
286 *--------------------------------------------------------------------------*/
287
288int hypre_StructMatvecCC0( double alpha,
289 hypre_StructMatrix *A,
290 hypre_StructVector *x,
291 hypre_StructVector *y,
292 hypre_BoxArrayArray *compute_box_aa,
293 hypre_IndexRef stride
294 )
295{
296 int i, j, si;
297 int ierr = 0;
298 double *Ap0;
299 double *Ap1;
300 double *Ap2;
301 double *Ap3;
302 double *Ap4;
303 double *Ap5;
304 double *Ap6;
305 int xoff0;
306 int xoff1;
307 int xoff2;
308 int xoff3;
309 int xoff4;
310 int xoff5;
311 int xoff6;
312 int Ai;
313 int xi;
314 hypre_BoxArray *compute_box_a;
315 hypre_Box *compute_box;
316
317 hypre_Box *A_data_box;
318 hypre_Box *x_data_box;
319 hypre_StructStencil *stencil;
320 hypre_Index *stencil_shape;
321 int stencil_size;
322
323 hypre_Box *y_data_box;
324 double *xp;
325 double *yp;
326 int depth;
327 hypre_Index loop_size;
328 int loopi, loopj, loopk;
329 hypre_IndexRef start;
330 int yi;
331
332 stencil = hypre_StructMatrixStencil(A);
333 stencil_shape = hypre_StructStencilShape(stencil);
334 stencil_size = hypre_StructStencilSize(stencil);
335
336 hypre_ForBoxArrayI(i, compute_box_aa)
337 {
338 compute_box_a = hypre_BoxArrayArrayBoxArray(compute_box_aa, i);
339
340 A_data_box = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A), i);
341 x_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), i);
342 y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
343
344 xp = hypre_StructVectorBoxData(x, i);
345 yp = hypre_StructVectorBoxData(y, i);
346
347 hypre_ForBoxI(j, compute_box_a)
348 {
349 compute_box = hypre_BoxArrayBox(compute_box_a, j);
350
351 hypre_BoxGetSize(compute_box, loop_size);
352 start = hypre_BoxIMin(compute_box);
353
354 /* unroll up to depth MAX_DEPTH */
355 for (si = 0; si < stencil_size; si+= MAX_DEPTH)
356 {
357 depth = hypre_min(MAX_DEPTH, (stencil_size -si));
358 switch(depth)
359 {
360 case 7:
361 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
362 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
363 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
364 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
365 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
366 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
367 Ap6 = hypre_StructMatrixBoxData(A, i, si+6);
368
369 xoff0 = hypre_BoxOffsetDistance(x_data_box,
370 stencil_shape[si+0]);
371 xoff1 = hypre_BoxOffsetDistance(x_data_box,
372 stencil_shape[si+1]);
373 xoff2 = hypre_BoxOffsetDistance(x_data_box,
374 stencil_shape[si+2]);
375 xoff3 = hypre_BoxOffsetDistance(x_data_box,
376 stencil_shape[si+3]);
377 xoff4 = hypre_BoxOffsetDistance(x_data_box,
378 stencil_shape[si+4]);
379 xoff5 = hypre_BoxOffsetDistance(x_data_box,
380 stencil_shape[si+5]);
381 xoff6 = hypre_BoxOffsetDistance(x_data_box,
382 stencil_shape[si+6]);
383
384 hypre_BoxLoop3Begin(loop_size,
385 A_data_box, start, stride, Ai,
386 x_data_box, start, stride, xi,
387 y_data_box, start, stride, yi);
388
389 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
390 {
391 yp[yi] +=
392 Ap0[Ai] * xp[xi + xoff0] +
393 Ap1[Ai] * xp[xi + xoff1] +
394 Ap2[Ai] * xp[xi + xoff2] +
395 Ap3[Ai] * xp[xi + xoff3] +
396 Ap4[Ai] * xp[xi + xoff4] +
397 Ap5[Ai] * xp[xi + xoff5] +
398 Ap6[Ai] * xp[xi + xoff6];
399 }
400 hypre_BoxLoop3End(Ai, xi, yi);
401
402 break;
403
404 case 6:
405 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
406 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
407 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
408 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
409 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
410 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
411
412 xoff0 = hypre_BoxOffsetDistance(x_data_box,
413 stencil_shape[si+0]);
414 xoff1 = hypre_BoxOffsetDistance(x_data_box,
415 stencil_shape[si+1]);
416 xoff2 = hypre_BoxOffsetDistance(x_data_box,
417 stencil_shape[si+2]);
418 xoff3 = hypre_BoxOffsetDistance(x_data_box,
419 stencil_shape[si+3]);
420 xoff4 = hypre_BoxOffsetDistance(x_data_box,
421 stencil_shape[si+4]);
422 xoff5 = hypre_BoxOffsetDistance(x_data_box,
423 stencil_shape[si+5]);
424
425 hypre_BoxLoop3Begin(loop_size,
426 A_data_box, start, stride, Ai,
427 x_data_box, start, stride, xi,
428 y_data_box, start, stride, yi);
429
430 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
431 {
432 yp[yi] +=
433 Ap0[Ai] * xp[xi + xoff0] +
434 Ap1[Ai] * xp[xi + xoff1] +
435 Ap2[Ai] * xp[xi + xoff2] +
436 Ap3[Ai] * xp[xi + xoff3] +
437 Ap4[Ai] * xp[xi + xoff4] +
438 Ap5[Ai] * xp[xi + xoff5];
439 }
440 hypre_BoxLoop3End(Ai, xi, yi);
441
442 break;
443
444 case 5:
445 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
446 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
447 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
448 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
449 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
450
451 xoff0 = hypre_BoxOffsetDistance(x_data_box,
452 stencil_shape[si+0]);
453 xoff1 = hypre_BoxOffsetDistance(x_data_box,
454 stencil_shape[si+1]);
455 xoff2 = hypre_BoxOffsetDistance(x_data_box,
456 stencil_shape[si+2]);
457 xoff3 = hypre_BoxOffsetDistance(x_data_box,
458 stencil_shape[si+3]);
459 xoff4 = hypre_BoxOffsetDistance(x_data_box,
460 stencil_shape[si+4]);
461
462 hypre_BoxLoop3Begin(loop_size,
463 A_data_box, start, stride, Ai,
464 x_data_box, start, stride, xi,
465 y_data_box, start, stride, yi);
466
467 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
468 {
469 yp[yi] +=
470 Ap0[Ai] * xp[xi + xoff0] +
471 Ap1[Ai] * xp[xi + xoff1] +
472 Ap2[Ai] * xp[xi + xoff2] +
473 Ap3[Ai] * xp[xi + xoff3] +
474 Ap4[Ai] * xp[xi + xoff4];
475 }
476 hypre_BoxLoop3End(Ai, xi, yi);
477
478 break;
479
480 case 4:
481 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
482 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
483 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
484 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
485
486 xoff0 = hypre_BoxOffsetDistance(x_data_box,
487 stencil_shape[si+0]);
488 xoff1 = hypre_BoxOffsetDistance(x_data_box,
489 stencil_shape[si+1]);
490 xoff2 = hypre_BoxOffsetDistance(x_data_box,
491 stencil_shape[si+2]);
492 xoff3 = hypre_BoxOffsetDistance(x_data_box,
493 stencil_shape[si+3]);
494
495 hypre_BoxLoop3Begin(loop_size,
496 A_data_box, start, stride, Ai,
497 x_data_box, start, stride, xi,
498 y_data_box, start, stride, yi);
499
500 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
501 {
502 yp[yi] +=
503 Ap0[Ai] * xp[xi + xoff0] +
504 Ap1[Ai] * xp[xi + xoff1] +
505 Ap2[Ai] * xp[xi + xoff2] +
506 Ap3[Ai] * xp[xi + xoff3];
507 }
508 hypre_BoxLoop3End(Ai, xi, yi);
509
510 break;
511
512 case 3:
513 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
514 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
515 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
516
517 xoff0 = hypre_BoxOffsetDistance(x_data_box,
518 stencil_shape[si+0]);
519 xoff1 = hypre_BoxOffsetDistance(x_data_box,
520 stencil_shape[si+1]);
521 xoff2 = hypre_BoxOffsetDistance(x_data_box,
522 stencil_shape[si+2]);
523
524 hypre_BoxLoop3Begin(loop_size,
525 A_data_box, start, stride, Ai,
526 x_data_box, start, stride, xi,
527 y_data_box, start, stride, yi);
528
529 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
530 {
531 yp[yi] +=
532 Ap0[Ai] * xp[xi + xoff0] +
533 Ap1[Ai] * xp[xi + xoff1] +
534 Ap2[Ai] * xp[xi + xoff2];
535 }
536 hypre_BoxLoop3End(Ai, xi, yi);
537
538 break;
539
540 case 2:
541 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
542 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
543
544 xoff0 = hypre_BoxOffsetDistance(x_data_box,
545 stencil_shape[si+0]);
546 xoff1 = hypre_BoxOffsetDistance(x_data_box,
547 stencil_shape[si+1]);
548
549 hypre_BoxLoop3Begin(loop_size,
550 A_data_box, start, stride, Ai,
551 x_data_box, start, stride, xi,
552 y_data_box, start, stride, yi);
553
554 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
555 {
556 yp[yi] +=
557 Ap0[Ai] * xp[xi + xoff0] +
558 Ap1[Ai] * xp[xi + xoff1];
559 }
560 hypre_BoxLoop3End(Ai, xi, yi);
561
562 break;
563
564 case 1:
565 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
566
567 xoff0 = hypre_BoxOffsetDistance(x_data_box,
568 stencil_shape[si+0]);
569
570 hypre_BoxLoop3Begin(loop_size,
571 A_data_box, start, stride, Ai,
572 x_data_box, start, stride, xi,
573 y_data_box, start, stride, yi);
574
575 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
576 {
577 yp[yi] +=
578 Ap0[Ai] * xp[xi + xoff0];
579 }
580 hypre_BoxLoop3End(Ai, xi, yi);
581
582 break;
583 }
584 }
585
586 if (alpha != 1.0)
587 {
588 hypre_BoxLoop1Begin(loop_size,
589 y_data_box, start, stride, yi);
590
591 hypre_BoxLoop1For(loopi, loopj, loopk, yi)
592 {
593 yp[yi] *= alpha;
594 }
595 hypre_BoxLoop1End(yi);
596 }
597 }
598 }
599 return ierr;
600}
601
602
603/*--------------------------------------------------------------------------
604 * hypre_StructMatvecCC1
605 * core of struct matvec computation, for the case constant_coefficient==1
606 *--------------------------------------------------------------------------*/
607
608int hypre_StructMatvecCC1( double alpha,
609 hypre_StructMatrix *A,
610 hypre_StructVector *x,
611 hypre_StructVector *y,
612 hypre_BoxArrayArray *compute_box_aa,
613 hypre_IndexRef stride
614 )
615{
616 int i, j, si;
617 int ierr = 0;
618 double *Ap0;
619 double *Ap1;
620 double *Ap2;
621 double *Ap3;
622 double *Ap4;
623 double *Ap5;
624 double *Ap6;
625 double AAp0;
626 double AAp1;
627 double AAp2;
628 double AAp3;
629 double AAp4;
630 double AAp5;
631 double AAp6;
632 int xoff0;
633 int xoff1;
634 int xoff2;
635 int xoff3;
636 int xoff4;
637 int xoff5;
638 int xoff6;
639 int Ai;
640 int xi;
641 hypre_BoxArray *compute_box_a;
642 hypre_Box *compute_box;
643
644 hypre_Box *A_data_box;
645 hypre_Box *x_data_box;
646 hypre_StructStencil *stencil;
647 hypre_Index *stencil_shape;
648 int stencil_size;
649
650 hypre_Box *y_data_box;
651 double *xp;
652 double *yp;
653 int depth;
654 hypre_Index loop_size;
655 int loopi, loopj, loopk;
656 hypre_IndexRef start;
657 int yi;
658
659 stencil = hypre_StructMatrixStencil(A);
660 stencil_shape = hypre_StructStencilShape(stencil);
661 stencil_size = hypre_StructStencilSize(stencil);
662
663 hypre_ForBoxArrayI(i, compute_box_aa)
664 {
665 compute_box_a = hypre_BoxArrayArrayBoxArray(compute_box_aa, i);
666
667 A_data_box = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A), i);
668 x_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), i);
669 y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
670
671 xp = hypre_StructVectorBoxData(x, i);
672 yp = hypre_StructVectorBoxData(y, i);
673
674 hypre_ForBoxI(j, compute_box_a)
675 {
676 compute_box = hypre_BoxArrayBox(compute_box_a, j);
677
678 hypre_BoxGetSize(compute_box, loop_size);
679 start = hypre_BoxIMin(compute_box);
680
681 Ai = hypre_CCBoxIndexRank( A_data_box, start );
682
683 /* unroll up to depth MAX_DEPTH */
684 for (si = 0; si < stencil_size; si+= MAX_DEPTH)
685 {
686 depth = hypre_min(MAX_DEPTH, (stencil_size -si));
687 switch(depth)
688 {
689 case 7:
690 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
691 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
692 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
693 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
694 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
695 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
696 Ap6 = hypre_StructMatrixBoxData(A, i, si+6);
697 AAp0 = Ap0[Ai]*alpha;
698 AAp1 = Ap1[Ai]*alpha;
699 AAp2 = Ap2[Ai]*alpha;
700 AAp3 = Ap3[Ai]*alpha;
701 AAp4 = Ap4[Ai]*alpha;
702 AAp5 = Ap5[Ai]*alpha;
703 AAp6 = Ap6[Ai]*alpha;
704
705 xoff0 = hypre_BoxOffsetDistance(x_data_box,
706 stencil_shape[si+0]);
707 xoff1 = hypre_BoxOffsetDistance(x_data_box,
708 stencil_shape[si+1]);
709 xoff2 = hypre_BoxOffsetDistance(x_data_box,
710 stencil_shape[si+2]);
711 xoff3 = hypre_BoxOffsetDistance(x_data_box,
712 stencil_shape[si+3]);
713 xoff4 = hypre_BoxOffsetDistance(x_data_box,
714 stencil_shape[si+4]);
715 xoff5 = hypre_BoxOffsetDistance(x_data_box,
716 stencil_shape[si+5]);
717 xoff6 = hypre_BoxOffsetDistance(x_data_box,
718 stencil_shape[si+6]);
719
720 hypre_BoxLoop2Begin(loop_size,
721 x_data_box, start, stride, xi,
722 y_data_box, start, stride, yi);
723
724 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
725 {
726 yp[yi] +=
727 AAp0 * xp[xi + xoff0] +
728 AAp1 * xp[xi + xoff1] +
729 AAp2 * xp[xi + xoff2] +
730 AAp3 * xp[xi + xoff3] +
731 AAp4 * xp[xi + xoff4] +
732 AAp5 * xp[xi + xoff5] +
733 AAp6 * xp[xi + xoff6];
734 }
735 hypre_BoxLoop2End(xi, yi);
736 break;
737
738 case 6:
739 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
740 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
741 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
742 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
743 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
744 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
745 AAp0 = Ap0[Ai]*alpha;
746 AAp1 = Ap1[Ai]*alpha;
747 AAp2 = Ap2[Ai]*alpha;
748 AAp3 = Ap3[Ai]*alpha;
749 AAp4 = Ap4[Ai]*alpha;
750 AAp5 = Ap5[Ai]*alpha;
751
752 xoff0 = hypre_BoxOffsetDistance(x_data_box,
753 stencil_shape[si+0]);
754 xoff1 = hypre_BoxOffsetDistance(x_data_box,
755 stencil_shape[si+1]);
756 xoff2 = hypre_BoxOffsetDistance(x_data_box,
757 stencil_shape[si+2]);
758 xoff3 = hypre_BoxOffsetDistance(x_data_box,
759 stencil_shape[si+3]);
760 xoff4 = hypre_BoxOffsetDistance(x_data_box,
761 stencil_shape[si+4]);
762 xoff5 = hypre_BoxOffsetDistance(x_data_box,
763 stencil_shape[si+5]);
764
765 hypre_BoxLoop2Begin(loop_size,
766 x_data_box, start, stride, xi,
767 y_data_box, start, stride, yi);
768
769 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
770 {
771 yp[yi] +=
772 AAp0 * xp[xi + xoff0] +
773 AAp1 * xp[xi + xoff1] +
774 AAp2 * xp[xi + xoff2] +
775 AAp3 * xp[xi + xoff3] +
776 AAp4 * xp[xi + xoff4] +
777 AAp5 * xp[xi + xoff5];
778 }
779 hypre_BoxLoop2End(xi, yi);
780 break;
781
782 case 5:
783 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
784 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
785 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
786 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
787 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
788 AAp0 = Ap0[Ai]*alpha;
789 AAp1 = Ap1[Ai]*alpha;
790 AAp2 = Ap2[Ai]*alpha;
791 AAp3 = Ap3[Ai]*alpha;
792 AAp4 = Ap4[Ai]*alpha;
793
794 xoff0 = hypre_BoxOffsetDistance(x_data_box,
795 stencil_shape[si+0]);
796 xoff1 = hypre_BoxOffsetDistance(x_data_box,
797 stencil_shape[si+1]);
798 xoff2 = hypre_BoxOffsetDistance(x_data_box,
799 stencil_shape[si+2]);
800 xoff3 = hypre_BoxOffsetDistance(x_data_box,
801 stencil_shape[si+3]);
802 xoff4 = hypre_BoxOffsetDistance(x_data_box,
803 stencil_shape[si+4]);
804
805 hypre_BoxLoop2Begin(loop_size,
806 x_data_box, start, stride, xi,
807 y_data_box, start, stride, yi);
808
809 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
810 {
811 yp[yi] +=
812 AAp0 * xp[xi + xoff0] +
813 AAp1 * xp[xi + xoff1] +
814 AAp2 * xp[xi + xoff2] +
815 AAp3 * xp[xi + xoff3] +
816 AAp4 * xp[xi + xoff4];
817 }
818 hypre_BoxLoop2End(xi, yi);
819 break;
820
821 case 4:
822 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
823 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
824 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
825 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
826 AAp0 = Ap0[Ai]*alpha;
827 AAp1 = Ap1[Ai]*alpha;
828 AAp2 = Ap2[Ai]*alpha;
829 AAp3 = Ap3[Ai]*alpha;
830
831 xoff0 = hypre_BoxOffsetDistance(x_data_box,
832 stencil_shape[si+0]);
833 xoff1 = hypre_BoxOffsetDistance(x_data_box,
834 stencil_shape[si+1]);
835 xoff2 = hypre_BoxOffsetDistance(x_data_box,
836 stencil_shape[si+2]);
837 xoff3 = hypre_BoxOffsetDistance(x_data_box,
838 stencil_shape[si+3]);
839
840 hypre_BoxLoop2Begin(loop_size,
841 x_data_box, start, stride, xi,
842 y_data_box, start, stride, yi);
843
844 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
845 {
846 yp[yi] +=
847 AAp0 * xp[xi + xoff0] +
848 AAp1 * xp[xi + xoff1] +
849 AAp2 * xp[xi + xoff2] +
850 AAp3 * xp[xi + xoff3];
851 }
852 hypre_BoxLoop2End(xi, yi);
853 break;
854
855 case 3:
856 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
857 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
858 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
859 AAp0 = Ap0[Ai]*alpha;
860 AAp1 = Ap1[Ai]*alpha;
861 AAp2 = Ap2[Ai]*alpha;
862
863 xoff0 = hypre_BoxOffsetDistance(x_data_box,
864 stencil_shape[si+0]);
865 xoff1 = hypre_BoxOffsetDistance(x_data_box,
866 stencil_shape[si+1]);
867 xoff2 = hypre_BoxOffsetDistance(x_data_box,
868 stencil_shape[si+2]);
869
870 hypre_BoxLoop2Begin(loop_size,
871 x_data_box, start, stride, xi,
872 y_data_box, start, stride, yi);
873
874 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
875 {
876 yp[yi] +=
877 AAp0 * xp[xi + xoff0] +
878 AAp1 * xp[xi + xoff1] +
879 AAp2 * xp[xi + xoff2];
880 }
881 hypre_BoxLoop2End(xi, yi);
882 break;
883
884 case 2:
885 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
886 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
887 AAp0 = Ap0[Ai]*alpha;
888 AAp1 = Ap1[Ai]*alpha;
889
890 xoff0 = hypre_BoxOffsetDistance(x_data_box,
891 stencil_shape[si+0]);
892 xoff1 = hypre_BoxOffsetDistance(x_data_box,
893 stencil_shape[si+1]);
894
895 hypre_BoxLoop2Begin(loop_size,
896 x_data_box, start, stride, xi,
897 y_data_box, start, stride, yi);
898
899 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
900 {
901 yp[yi] +=
902 AAp0 * xp[xi + xoff0] +
903 AAp1 * xp[xi + xoff1];
904 }
905 hypre_BoxLoop2End(xi, yi);
906 break;
907
908 case 1:
909 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
910 AAp0 = Ap0[Ai]*alpha;
911
912 xoff0 = hypre_BoxOffsetDistance(x_data_box,
913 stencil_shape[si+0]);
914
915 hypre_BoxLoop2Begin(loop_size,
916 x_data_box, start, stride, xi,
917 y_data_box, start, stride, yi);
918
919 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
920 {
921 yp[yi] +=
922 AAp0 * xp[xi + xoff0];
923 }
924 hypre_BoxLoop2End(xi, yi);
925 }
926 }
927 }
928 }
929
930 return ierr;
931}
932
933
934/*--------------------------------------------------------------------------
935 * hypre_StructMatvecCC2
936 * core of struct matvec computation, for the case constant_coefficient==2
937 *--------------------------------------------------------------------------*/
938
939int hypre_StructMatvecCC2( double alpha,
940 hypre_StructMatrix *A,
941 hypre_StructVector *x,
942 hypre_StructVector *y,
943 hypre_BoxArrayArray *compute_box_aa,
944 hypre_IndexRef stride
945 )
946{
947 int i, j, si;
948 int ierr = 0;
949 double *Ap0;
950 double *Ap1;
951 double *Ap2;
952 double *Ap3;
953 double *Ap4;
954 double *Ap5;
955 double *Ap6;
956 double AAp0;
957 double AAp1;
958 double AAp2;
959 double AAp3;
960 double AAp4;
961 double AAp5;
962 double AAp6;
963 int xoff0;
964 int xoff1;
965 int xoff2;
966 int xoff3;
967 int xoff4;
968 int xoff5;
969 int xoff6;
970 int si_center, center_rank;
971 hypre_Index center_index;
972 int Ai, Ai_CC;
973 int xi;
974 hypre_BoxArray *compute_box_a;
975 hypre_Box *compute_box;
976
977 hypre_Box *A_data_box;
978 hypre_Box *x_data_box;
979 hypre_StructStencil *stencil;
980 hypre_Index *stencil_shape;
981 int stencil_size;
982
983 hypre_Box *y_data_box;
984 double *xp;
985 double *yp;
986 int depth;
987 hypre_Index loop_size;
988 int loopi, loopj, loopk;
989 hypre_IndexRef start;
990 int yi;
991
992 stencil = hypre_StructMatrixStencil(A);
993 stencil_shape = hypre_StructStencilShape(stencil);
994 stencil_size = hypre_StructStencilSize(stencil);
995
996
997 hypre_ForBoxArrayI(i, compute_box_aa)
998 {
999 compute_box_a = hypre_BoxArrayArrayBoxArray(compute_box_aa, i);
1000
1001 A_data_box = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(A), i);
1002 x_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(x), i);
1003 y_data_box = hypre_BoxArrayBox(hypre_StructVectorDataSpace(y), i);
1004
1005 xp = hypre_StructVectorBoxData(x, i);
1006 yp = hypre_StructVectorBoxData(y, i);
1007
1008 hypre_ForBoxI(j, compute_box_a)
1009 {
1010 compute_box = hypre_BoxArrayBox(compute_box_a, j);
1011
1012 hypre_BoxGetSize(compute_box, loop_size);
1013 start = hypre_BoxIMin(compute_box);
1014
1015 Ai_CC = hypre_CCBoxIndexRank( A_data_box, start );
1016
1017 /* Find the stencil index for the center of the stencil, which
1018 makes the matrix diagonal. This is the variable coefficient
1019 part of the matrix, so will get different treatment...*/
1020 hypre_SetIndex(center_index, 0, 0, 0);
1021 center_rank = hypre_StructStencilElementRank( stencil, center_index );
1022 si_center = center_rank;
1023
1024 /* unroll up to depth MAX_DEPTH
1025 Only the constant coefficient part of the matrix is referenced here,
1026 the center (variable) coefficient part is deferred. */
1027 for (si = 0; si < stencil_size; si+= MAX_DEPTH)
1028 {
1029 depth = hypre_min(MAX_DEPTH, (stencil_size -si));
1030 switch(depth)
1031 {
1032 case 7:
1033 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1034 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1035 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
1036 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
1037 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
1038 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
1039 Ap6 = hypre_StructMatrixBoxData(A, i, si+6);
1040 AAp0 = Ap0[Ai_CC];
1041 AAp1 = Ap1[Ai_CC];
1042 AAp2 = Ap2[Ai_CC];
1043 AAp3 = Ap3[Ai_CC];
1044 AAp4 = Ap4[Ai_CC];
1045 AAp5 = Ap5[Ai_CC];
1046 AAp6 = Ap6[Ai_CC];
1047 if ( (0 <= si_center-si) && (si_center-si < 7) )
1048 {
1049 switch ( si_center-si )
1050 {
1051 case 0: AAp0 = 0; break;
1052 case 1: AAp1 = 0; break;
1053 case 2: AAp2 = 0; break;
1054 case 3: AAp3 = 0; break;
1055 case 4: AAp4 = 0; break;
1056 case 5: AAp5 = 0; break;
1057 case 6: AAp6 = 0; break;
1058 }
1059 }
1060
1061 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1062 stencil_shape[si+0]);
1063 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1064 stencil_shape[si+0]);
1065 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1066 stencil_shape[si+1]);
1067 xoff2 = hypre_BoxOffsetDistance(x_data_box,
1068 stencil_shape[si+2]);
1069 xoff3 = hypre_BoxOffsetDistance(x_data_box,
1070 stencil_shape[si+3]);
1071 xoff4 = hypre_BoxOffsetDistance(x_data_box,
1072 stencil_shape[si+4]);
1073 xoff5 = hypre_BoxOffsetDistance(x_data_box,
1074 stencil_shape[si+5]);
1075 xoff6 = hypre_BoxOffsetDistance(x_data_box,
1076 stencil_shape[si+6]);
1077
1078 hypre_BoxLoop2Begin(loop_size,
1079 x_data_box, start, stride, xi,
1080 y_data_box, start, stride, yi);
1081
1082 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1083 {
1084 yp[yi] +=
1085 AAp0 * xp[xi + xoff0] +
1086 AAp1 * xp[xi + xoff1] +
1087 AAp2 * xp[xi + xoff2] +
1088 AAp3 * xp[xi + xoff3] +
1089 AAp4 * xp[xi + xoff4] +
1090 AAp5 * xp[xi + xoff5] +
1091 AAp6 * xp[xi + xoff6];
1092 }
1093 hypre_BoxLoop2End(xi, yi);
1094
1095 break;
1096
1097 case 6:
1098 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1099 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1100 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
1101 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
1102 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
1103 Ap5 = hypre_StructMatrixBoxData(A, i, si+5);
1104 AAp0 = Ap0[Ai_CC];
1105 AAp1 = Ap1[Ai_CC];
1106 AAp2 = Ap2[Ai_CC];
1107 AAp3 = Ap3[Ai_CC];
1108 AAp4 = Ap4[Ai_CC];
1109 AAp5 = Ap5[Ai_CC];
1110 if ( (0 <= si_center-si) && (si_center-si < 6) )
1111 {
1112 switch ( si_center-si )
1113 {
1114 case 0: AAp0 = 0; break;
1115 case 1: AAp1 = 0; break;
1116 case 2: AAp2 = 0; break;
1117 case 3: AAp3 = 0; break;
1118 case 4: AAp4 = 0; break;
1119 case 5: AAp5 = 0; break;
1120 }
1121 }
1122
1123 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1124 stencil_shape[si+0]);
1125 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1126 stencil_shape[si+1]);
1127 xoff2 = hypre_BoxOffsetDistance(x_data_box,
1128 stencil_shape[si+2]);
1129 xoff3 = hypre_BoxOffsetDistance(x_data_box,
1130 stencil_shape[si+3]);
1131 xoff4 = hypre_BoxOffsetDistance(x_data_box,
1132 stencil_shape[si+4]);
1133 xoff5 = hypre_BoxOffsetDistance(x_data_box,
1134 stencil_shape[si+5]);
1135
1136 hypre_BoxLoop2Begin(loop_size,
1137 x_data_box, start, stride, xi,
1138 y_data_box, start, stride, yi);
1139
1140 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1141 {
1142 yp[yi] +=
1143 AAp0 * xp[xi + xoff0] +
1144 AAp1 * xp[xi + xoff1] +
1145 AAp2 * xp[xi + xoff2] +
1146 AAp3 * xp[xi + xoff3] +
1147 AAp4 * xp[xi + xoff4] +
1148 AAp5 * xp[xi + xoff5];
1149 }
1150 hypre_BoxLoop2End(xi, yi);
1151 break;
1152
1153 case 5:
1154 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1155 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1156 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
1157 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
1158 Ap4 = hypre_StructMatrixBoxData(A, i, si+4);
1159 AAp0 = Ap0[Ai_CC];
1160 AAp1 = Ap1[Ai_CC];
1161 AAp2 = Ap2[Ai_CC];
1162 AAp3 = Ap3[Ai_CC];
1163 AAp4 = Ap4[Ai_CC];
1164 if ( (0 <= si_center-si) && (si_center-si < 5) )
1165 {
1166 switch ( si_center-si )
1167 {
1168 case 0: AAp0 = 0; break;
1169 case 1: AAp1 = 0; break;
1170 case 2: AAp2 = 0; break;
1171 case 3: AAp3 = 0; break;
1172 case 4: AAp4 = 0; break;
1173 }
1174 }
1175
1176 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1177 stencil_shape[si+0]);
1178 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1179 stencil_shape[si+1]);
1180 xoff2 = hypre_BoxOffsetDistance(x_data_box,
1181 stencil_shape[si+2]);
1182 xoff3 = hypre_BoxOffsetDistance(x_data_box,
1183 stencil_shape[si+3]);
1184 xoff4 = hypre_BoxOffsetDistance(x_data_box,
1185 stencil_shape[si+4]);
1186
1187 hypre_BoxLoop2Begin(loop_size,
1188 x_data_box, start, stride, xi,
1189 y_data_box, start, stride, yi);
1190
1191 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1192 {
1193 yp[yi] +=
1194 AAp0 * xp[xi + xoff0] +
1195 AAp1 * xp[xi + xoff1] +
1196 AAp2 * xp[xi + xoff2] +
1197 AAp3 * xp[xi + xoff3] +
1198 AAp4 * xp[xi + xoff4];
1199 }
1200 hypre_BoxLoop2End(xi, yi);
1201 break;
1202
1203 case 4:
1204 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1205 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1206 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
1207 Ap3 = hypre_StructMatrixBoxData(A, i, si+3);
1208 AAp0 = Ap0[Ai_CC];
1209 AAp1 = Ap1[Ai_CC];
1210 AAp2 = Ap2[Ai_CC];
1211 AAp3 = Ap3[Ai_CC];
1212 if ( (0 <= si_center-si) && (si_center-si < 4) )
1213 {
1214 switch ( si_center-si )
1215 {
1216 case 0: AAp0 = 0; break;
1217 case 1: AAp1 = 0; break;
1218 case 2: AAp2 = 0; break;
1219 case 3: AAp3 = 0; break;
1220 }
1221 }
1222
1223 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1224 stencil_shape[si+0]);
1225 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1226 stencil_shape[si+1]);
1227 xoff2 = hypre_BoxOffsetDistance(x_data_box,
1228 stencil_shape[si+2]);
1229 xoff3 = hypre_BoxOffsetDistance(x_data_box,
1230 stencil_shape[si+3]);
1231
1232 hypre_BoxLoop2Begin(loop_size,
1233 x_data_box, start, stride, xi,
1234 y_data_box, start, stride, yi);
1235
1236 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1237 {
1238 yp[yi] +=
1239 AAp0 * xp[xi + xoff0] +
1240 AAp1 * xp[xi + xoff1] +
1241 AAp2 * xp[xi + xoff2] +
1242 AAp3 * xp[xi + xoff3];
1243 }
1244 hypre_BoxLoop2End(xi, yi);
1245 break;
1246
1247 case 3:
1248 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1249 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1250 Ap2 = hypre_StructMatrixBoxData(A, i, si+2);
1251 AAp0 = Ap0[Ai_CC];
1252 AAp1 = Ap1[Ai_CC];
1253 AAp2 = Ap2[Ai_CC];
1254 if ( (0 <= si_center-si) && (si_center-si < 3) )
1255 {
1256 switch ( si_center-si )
1257 {
1258 case 0: AAp0 = 0; break;
1259 case 1: AAp1 = 0; break;
1260 case 2: AAp2 = 0; break;
1261 }
1262 }
1263
1264 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1265 stencil_shape[si+0]);
1266 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1267 stencil_shape[si+1]);
1268 xoff2 = hypre_BoxOffsetDistance(x_data_box,
1269 stencil_shape[si+2]);
1270
1271 hypre_BoxLoop2Begin(loop_size,
1272 x_data_box, start, stride, xi,
1273 y_data_box, start, stride, yi);
1274
1275 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1276 {
1277 yp[yi] +=
1278 AAp0 * xp[xi + xoff0] +
1279 AAp1 * xp[xi + xoff1] +
1280 AAp2 * xp[xi + xoff2];
1281 }
1282 hypre_BoxLoop2End(xi, yi);
1283 break;
1284
1285 case 2:
1286 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1287 Ap1 = hypre_StructMatrixBoxData(A, i, si+1);
1288 AAp0 = Ap0[Ai_CC];
1289 AAp1 = Ap1[Ai_CC];
1290 if ( (0 <= si_center-si) && (si_center-si < 2) )
1291 {
1292 switch ( si_center-si )
1293 {
1294 case 0: AAp0 = 0; break;
1295 case 1: AAp1 = 0; break;
1296 }
1297 }
1298
1299 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1300 stencil_shape[si+0]);
1301 xoff1 = hypre_BoxOffsetDistance(x_data_box,
1302 stencil_shape[si+1]);
1303
1304 hypre_BoxLoop2Begin(loop_size,
1305 x_data_box, start, stride, xi,
1306 y_data_box, start, stride, yi);
1307
1308 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1309 {
1310 yp[yi] +=
1311 AAp0 * xp[xi + xoff0] +
1312 AAp1 * xp[xi + xoff1];
1313 }
1314 hypre_BoxLoop2End(xi, yi);
1315 break;
1316
1317 case 1:
1318 Ap0 = hypre_StructMatrixBoxData(A, i, si+0);
1319 AAp0 = Ap0[Ai_CC];
1320 if ( si_center-si == 0 )
1321 {
1322 AAp0 = 0;
1323 }
1324
1325 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1326 stencil_shape[si+0]);
1327
1328 hypre_BoxLoop2Begin(loop_size,
1329 x_data_box, start, stride, xi,
1330 y_data_box, start, stride, yi);
1331
1332 hypre_BoxLoop2For(loopi, loopj, loopk, xi, yi)
1333 {
1334 yp[yi] +=
1335 AAp0 * xp[xi + xoff0];
1336 }
1337 hypre_BoxLoop2End(xi, yi);
1338
1339 break;
1340 }
1341 }
1342
1343 Ap0 = hypre_StructMatrixBoxData(A, i, si_center);
1344 xoff0 = hypre_BoxOffsetDistance(x_data_box,
1345 stencil_shape[si_center]);
1346 if (alpha!= 1.0 )
1347 {
1348 hypre_BoxLoop3Begin(loop_size,
1349 A_data_box, start, stride, Ai,
1350 x_data_box, start, stride, xi,
1351 y_data_box, start, stride, yi);
1352
1353 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
1354 {
1355 yp[yi] = alpha * ( yp[yi] +
1356 Ap0[Ai] * xp[xi + xoff0] );
1357 }
1358 hypre_BoxLoop3End(Ai, xi, yi);
1359 }
1360 else
1361 {
1362 hypre_BoxLoop3Begin(loop_size,
1363 A_data_box, start, stride, Ai,
1364 x_data_box, start, stride, xi,
1365 y_data_box, start, stride, yi);
1366
1367 hypre_BoxLoop3For(loopi, loopj, loopk, Ai, xi, yi)
1368 {
1369 yp[yi] +=
1370 Ap0[Ai] * xp[xi + xoff0];
1371 }
1372 hypre_BoxLoop3End(Ai, xi, yi);
1373
1374 }
1375
1376 }
1377 }
1378
1379 return ierr;
1380}
1381
1382
1383/*--------------------------------------------------------------------------
1384 * hypre_StructMatvecDestroy
1385 *--------------------------------------------------------------------------*/
1386
1387int
1388hypre_StructMatvecDestroy( void *matvec_vdata )
1389{
1390 int ierr = 0;
1391
1392 hypre_StructMatvecData *matvec_data = matvec_vdata;
1393
1394 if (matvec_data)
1395 {
1396 hypre_StructMatrixDestroy(matvec_data -> A);
1397 hypre_StructVectorDestroy(matvec_data -> x);
1398 hypre_ComputePkgDestroy(matvec_data -> compute_pkg );
1399 hypre_TFree(matvec_data);
1400 }
1401
1402 return ierr;
1403}
1404
1405/*--------------------------------------------------------------------------
1406 * hypre_StructMatvec
1407 *--------------------------------------------------------------------------*/
1408
1409int
1410hypre_StructMatvec( double alpha,
1411 hypre_StructMatrix *A,
1412 hypre_StructVector *x,
1413 double beta,
1414 hypre_StructVector *y )
1415{
1416 int ierr = 0;
1417
1418 void *matvec_data;
1419
1420 matvec_data = hypre_StructMatvecCreate();
1421 ierr = hypre_StructMatvecSetup(matvec_data, A, x);
1422 ierr = hypre_StructMatvecCompute(matvec_data, alpha, A, x, beta, y);
1423 ierr = hypre_StructMatvecDestroy(matvec_data);
1424
1425 return ierr;
1426}
Note: See TracBrowser for help on using the repository browser.