Skip to content

Commit c2849f3

Browse files
jedbrownjedbrown
committed
Merge branch 'jeremy/vec-ops' [PR #744]
* jeremy/vec-ops: vec - add tests for ceed compatibility for AXPY and PointwiseMult vec - add PointwiseMult and AXPY convenience functions
2 parents 49dd03b + 2d04630 commit c2849f3

6 files changed

Lines changed: 242 additions & 1 deletion

File tree

‎include/ceed-impl.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -138,6 +138,8 @@ struct CeedVector_private {
138138
int (*RestoreArray)(CeedVector);
139139
int (*RestoreArrayRead)(CeedVector);
140140
int (*Norm)(CeedVector, CeedNormType, CeedScalar *);
141+
int (*AXPY)(CeedVector, CeedScalar, CeedVector);
142+
int (*PointwiseMult)(CeedVector, CeedVector, CeedVector);
141143
int (*Reciprocal)(CeedVector);
142144
int (*Destroy)(CeedVector);
143145
int ref_count;

‎include/ceed/ceed.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -236,7 +236,7 @@ CEED_EXTERN int CeedGetVersion(int *major, int *minor, int *patch,
236236
/// data is likely modified or corrupted.
237237
/// @ingroup Ceed
238238
typedef enum {
239-
/// Sucess error code
239+
/// Success error code
240240
CEED_ERROR_SUCCESS = 0,
241241
/// Minor error, generic
242242
CEED_ERROR_MINOR = 1,
@@ -321,6 +321,8 @@ CEED_EXTERN int CeedVectorRestoreArrayRead(CeedVector vec,
321321
const CeedScalar **array);
322322
CEED_EXTERN int CeedVectorNorm(CeedVector vec, CeedNormType type,
323323
CeedScalar *norm);
324+
CEED_EXTERN int CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x);
325+
CEED_EXTERN int CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y);
324326
CEED_EXTERN int CeedVectorReciprocal(CeedVector vec);
325327
CEED_EXTERN int CeedVectorView(CeedVector vec, const char *fp_fmt, FILE *stream);
326328
CEED_EXTERN int CeedVectorGetLength(CeedVector vec, CeedInt *length);

‎interface/ceed-vector.c‎

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -531,6 +531,132 @@ int CeedVectorNorm(CeedVector vec, CeedNormType norm_type, CeedScalar *norm) {
531531
return CEED_ERROR_SUCCESS;
532532
}
533533

534+
/**
535+
@brief Compute y = alpha x + y
536+
537+
@param y[in,out] target vector for sum
538+
@param alpha[in] scaling factor
539+
@param x[in] second vector, must be different than y
540+
541+
@return An error code: 0 - success, otherwise - failure
542+
543+
@ref User
544+
**/
545+
int CeedVectorAXPY(CeedVector y, CeedScalar alpha, CeedVector x) {
546+
int ierr;
547+
CeedScalar *y_array;
548+
CeedScalar const *x_array;
549+
CeedInt n_x, n_y;
550+
551+
ierr = CeedVectorGetLength(y, &n_y); CeedChk(ierr);
552+
ierr = CeedVectorGetLength(x, &n_x); CeedChk(ierr);
553+
if (n_x != n_y)
554+
// LCOV_EXCL_START
555+
return CeedError(y->ceed, CEED_ERROR_UNSUPPORTED,
556+
"Cannot add vector of different lengths");
557+
// LCOV_EXCL_STOP
558+
if (x == y)
559+
// LCOV_EXCL_START
560+
return CeedError(y->ceed, CEED_ERROR_UNSUPPORTED,
561+
"Cannot use same vector for x and y in CeedVectorAXPY");
562+
// LCOV_EXCL_STOP
563+
564+
Ceed ceed_parent_x, ceed_parent_y;
565+
ierr = CeedGetParent(x->ceed, &ceed_parent_x); CeedChk(ierr);
566+
ierr = CeedGetParent(y->ceed, &ceed_parent_y); CeedChk(ierr);
567+
if (ceed_parent_x != ceed_parent_y)
568+
// LCOV_EXCL_START
569+
return CeedError(y->ceed, CEED_ERROR_INCOMPATIBLE,
570+
"Vectors x and y must be created by the same Ceed context");
571+
// LCOV_EXCL_STOP
572+
573+
// Backend implementation
574+
if (y->AXPY)
575+
return y->AXPY(y, alpha, x);
576+
577+
// Default implementation
578+
ierr = CeedVectorGetArray(y, CEED_MEM_HOST, &y_array); CeedChk(ierr);
579+
ierr = CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); CeedChk(ierr);
580+
581+
for (CeedInt i=0; i<n_y; i++)
582+
y_array[i] += alpha * x_array[i];
583+
584+
ierr = CeedVectorRestoreArray(y, &y_array); CeedChk(ierr);
585+
ierr = CeedVectorRestoreArrayRead(x, &x_array); CeedChk(ierr);
586+
587+
return CEED_ERROR_SUCCESS;
588+
}
589+
590+
/**
591+
@brief Compute the pointwise multiplication w = x * y. Any
592+
subset of x, y, and w may be the same vector.
593+
594+
@param w[out] target vector for the product
595+
@param x[in] first vector for product
596+
@param y[in] second vector for the product
597+
598+
@return An error code: 0 - success, otherwise - failure
599+
600+
@ ref User
601+
**/
602+
int CeedVectorPointwiseMult(CeedVector w, CeedVector x, CeedVector y) {
603+
int ierr;
604+
CeedScalar *w_array;
605+
CeedScalar const *x_array, *y_array;
606+
CeedInt n_w, n_x, n_y;
607+
608+
ierr = CeedVectorGetLength(w, &n_w); CeedChk(ierr);
609+
ierr = CeedVectorGetLength(x, &n_x); CeedChk(ierr);
610+
ierr = CeedVectorGetLength(y, &n_y); CeedChk(ierr);
611+
if (n_w != n_x || n_w != n_y)
612+
// LCOV_EXCL_START
613+
return CeedError(w->ceed, CEED_ERROR_UNSUPPORTED,
614+
"Cannot multiply vectors of different lengths");
615+
// LCOV_EXCL_STOP
616+
617+
Ceed ceed_parent_w, ceed_parent_x, ceed_parent_y;
618+
ierr = CeedGetParent(w->ceed, &ceed_parent_w); CeedChk(ierr);
619+
ierr = CeedGetParent(x->ceed, &ceed_parent_x); CeedChk(ierr);
620+
ierr = CeedGetParent(y->ceed, &ceed_parent_y); CeedChk(ierr);
621+
if ((ceed_parent_w != ceed_parent_y) ||
622+
(ceed_parent_w != ceed_parent_y))
623+
// LCOV_EXCL_START
624+
return CeedError(w->ceed, CEED_ERROR_INCOMPATIBLE,
625+
"Vectors w, x, and y must be created by the same Ceed context");
626+
// LCOV_EXCL_STOP
627+
628+
// Backend implementation
629+
if (w->PointwiseMult)
630+
return w->PointwiseMult(w, x, y);
631+
632+
// Default implementation
633+
ierr = CeedVectorGetArray(w, CEED_MEM_HOST, &w_array); CeedChk(ierr);
634+
if (x != w) {
635+
ierr = CeedVectorGetArrayRead(x, CEED_MEM_HOST, &x_array); CeedChk(ierr);
636+
} else {
637+
x_array = w_array;
638+
}
639+
if (y != w && y != x) {
640+
ierr = CeedVectorGetArrayRead(y, CEED_MEM_HOST, &y_array); CeedChk(ierr);
641+
} else if (y != x) {
642+
y_array = w_array;
643+
} else {
644+
y_array = x_array;
645+
}
646+
647+
for (CeedInt i=0; i<n_w; i++)
648+
w_array[i] = x_array[i] * y_array[i];
649+
650+
if (y != w && y != x) {
651+
ierr = CeedVectorRestoreArrayRead(y, &y_array); CeedChk(ierr);
652+
}
653+
if (x != w) {
654+
ierr = CeedVectorRestoreArrayRead(x, &x_array); CeedChk(ierr);
655+
}
656+
ierr = CeedVectorRestoreArray(w, &w_array); CeedChk(ierr);
657+
return CEED_ERROR_SUCCESS;
658+
}
659+
534660
/**
535661
@brief Take the reciprocal of a CeedVector.
536662

‎interface/ceed.c‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -770,6 +770,8 @@ int CeedInit(const char *resource, Ceed *ceed) {
770770
CEED_FTABLE_ENTRY(CeedVector, RestoreArray),
771771
CEED_FTABLE_ENTRY(CeedVector, RestoreArrayRead),
772772
CEED_FTABLE_ENTRY(CeedVector, Norm),
773+
CEED_FTABLE_ENTRY(CeedVector, AXPY),
774+
CEED_FTABLE_ENTRY(CeedVector, PointwiseMult),
773775
CEED_FTABLE_ENTRY(CeedVector, Reciprocal),
774776
CEED_FTABLE_ENTRY(CeedVector, Destroy),
775777
CEED_FTABLE_ENTRY(CeedElemRestriction, Apply),

‎tests/t121-vector.c‎

Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
/// @file
2+
/// Test summing of a pair of vectors
3+
/// \test Test summing of a pair of vectors
4+
#include <ceed.h>
5+
#include <math.h>
6+
7+
int main(int argc, char **argv) {
8+
Ceed ceed;
9+
CeedVector x, y;
10+
CeedInt n;
11+
CeedScalar a[10];
12+
const CeedScalar *b;
13+
14+
CeedInit(argv[1], &ceed);
15+
16+
n = 10;
17+
CeedVectorCreate(ceed, n, &x);
18+
CeedVectorCreate(ceed, n, &y);
19+
for (CeedInt i=0; i<n; i++)
20+
a[i] = 10 + i;
21+
CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, a);
22+
CeedVectorSetArray(y, CEED_MEM_HOST, CEED_COPY_VALUES, a);
23+
24+
CeedVectorAXPY(y, -0.5, x);
25+
26+
CeedVectorGetArrayRead(y, CEED_MEM_HOST, &b);
27+
for (CeedInt i=0; i<n; i++)
28+
if (fabs(b[i] - (10.0 + i)/2 ) > 1e-14)
29+
// LCOV_EXCL_START
30+
printf("Error in alpha x + y, computed: %f actual: %f\n", b[i],
31+
(10.0 + i)/2);
32+
// LCOV_EXCL_STOP
33+
CeedVectorRestoreArrayRead(y, &b);
34+
35+
CeedVectorDestroy(&x);
36+
CeedVectorDestroy(&y);
37+
CeedDestroy(&ceed);
38+
return 0;
39+
}

‎tests/t122-vector.c‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
/// @file
2+
/// Test poinwise muliplication of a pair of vectors
3+
/// \test Test poinwise muliplication of a pair of vectors
4+
#include <ceed.h>
5+
#include <math.h>
6+
7+
int main(int argc, char **argv) {
8+
Ceed ceed;
9+
CeedVector x, y, w;
10+
CeedInt n;
11+
CeedScalar a[10];
12+
const CeedScalar *b;
13+
14+
CeedInit(argv[1], &ceed);
15+
16+
n = 10;
17+
CeedVectorCreate(ceed, n, &x);
18+
CeedVectorCreate(ceed, n, &y);
19+
CeedVectorCreate(ceed, n, &w);
20+
for (CeedInt i=0; i<n; i++)
21+
a[i] = i;
22+
CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, a);
23+
CeedVectorSetArray(y, CEED_MEM_HOST, CEED_COPY_VALUES, a);
24+
25+
// Test multiplying two vectors into third
26+
CeedVectorPointwiseMult(w, x, y);
27+
CeedVectorGetArrayRead(w, CEED_MEM_HOST, &b);
28+
for (CeedInt i=0; i<n; i++)
29+
if (fabs(b[i] - i*i ) > 1e-14)
30+
// LCOV_EXCL_START
31+
printf("Error in alpha x + y, computed: %f actual: %f\n", b[i], 1.0*i*i);
32+
// LCOV_EXCL_STOP
33+
CeedVectorRestoreArrayRead(w, &b);
34+
35+
// Test multiplying two vectors into one of the two
36+
CeedVectorPointwiseMult(w, w, y);
37+
CeedVectorGetArrayRead(w, CEED_MEM_HOST, &b);
38+
for (CeedInt i=0; i<n; i++)
39+
if (fabs(b[i] - i*i*i ) > 1e-14)
40+
// LCOV_EXCL_START
41+
printf("Error in alpha x + y, computed: %f actual: %f\n", b[i], 1.0*i*i*i);
42+
// LCOV_EXCL_STOP
43+
CeedVectorRestoreArrayRead(w, &b);
44+
45+
// Test multiplying two vectors into one of the two
46+
CeedVectorPointwiseMult(w, x, w);
47+
CeedVectorGetArrayRead(w, CEED_MEM_HOST, &b);
48+
for (CeedInt i=0; i<n; i++)
49+
if (fabs(b[i] - i*i*i*i ) > 1e-14)
50+
// LCOV_EXCL_START
51+
printf("Error in alpha x + y, computed: %f actual: %f\n", b[i], 1.0*i*i*i*i);
52+
// LCOV_EXCL_STOP
53+
CeedVectorRestoreArrayRead(w, &b);
54+
55+
// Test multiplying vector by itself and putting product into self
56+
CeedVectorPointwiseMult(y, y, y);
57+
CeedVectorGetArrayRead(y, CEED_MEM_HOST, &b);
58+
for (CeedInt i=0; i<n; i++)
59+
if (fabs(b[i] - i*i ) > 1e-14)
60+
// LCOV_EXCL_START
61+
printf("Error in alpha x + y, computed: %f actual: %f\n", b[i], 1.0*i*i);
62+
// LCOV_EXCL_STOP
63+
CeedVectorRestoreArrayRead(y, &b);
64+
65+
CeedVectorDestroy(&x);
66+
CeedVectorDestroy(&y);
67+
CeedVectorDestroy(&w);
68+
CeedDestroy(&ceed);
69+
return 0;
70+
}

0 commit comments

Comments
 (0)