Skip to content

Commit d965c7a

Browse files
jeremyltjedbrown
jeremylt
and
jedbrown
authored
CPU Impl of AssemblePointBlockDiagonal (CEED#503)
* Operator - add AssemblePointBlockDiagonal for CPU backends, with test * CUDA - add point block diagonal not supported message * make style * Operator - improve point block description * OCCA - explicitly remove OCCA fallback to CPU assembly functions, will update after new OCCA backend * Op - remove gap removal in point block diagonal * Op - update diagonal assembly documentation * Update backends/ref/ceed-ref-operator.c Co-authored-by: Jed Brown <[email protected]> * style - fix extra space in * with nopad Co-authored-by: Jed Brown <[email protected]>
1 parent 4eb74d6 commit d965c7a

17 files changed

+488
-21
lines changed

‎backends/cuda-gen/ceed-cuda-gen-operator.c

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -205,6 +205,17 @@ static int CeedOperatorAssembleLinearDiagonal_Cuda(CeedOperator op) {
205205
"Backend does not implement Operator diagonal assembly");
206206
}
207207

208+
//------------------------------------------------------------------------------
209+
// Assemble linear point block diagonal not supported
210+
//------------------------------------------------------------------------------
211+
static int CeedOperatorAssembleLinearPointBlockDiagonal_Cuda(CeedOperator op) {
212+
int ierr;
213+
Ceed ceed;
214+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
215+
return CeedError(ceed, 1,
216+
"Backend does not implement Operator point block diagonal assembly");
217+
}
218+
208219
//------------------------------------------------------------------------------
209220
// Create FDM element inverse not supported
210221
//------------------------------------------------------------------------------
@@ -233,6 +244,10 @@ int CeedOperatorCreate_Cuda_gen(CeedOperator op) {
233244
ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearDiagonal",
234245
CeedOperatorAssembleLinearDiagonal_Cuda);
235246
CeedChk(ierr);
247+
ierr = CeedSetBackendFunction(ceed, "Operator", op,
248+
"AssembleLinearPointBlockDiagonal",
249+
CeedOperatorAssembleLinearPointBlockDiagonal_Cuda);
250+
CeedChk(ierr);
236251
ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
237252
CeedOperatorCreateFDMElementInverse_Cuda);
238253
CeedChk(ierr);

‎backends/cuda/ceed-cuda-operator.c

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -444,6 +444,17 @@ static int CeedOperatorAssembleLinearDiagonal_Cuda(CeedOperator op) {
444444
"Backend does not implement Operator diagonal assembly");
445445
}
446446

447+
//------------------------------------------------------------------------------
448+
// Assemble linear point block diagonal not supported
449+
//------------------------------------------------------------------------------
450+
static int CeedOperatorAssembleLinearPointBlockDiagonal_Cuda(CeedOperator op) {
451+
int ierr;
452+
Ceed ceed;
453+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
454+
return CeedError(ceed, 1,
455+
"Backend does not implement Operator point block diagonal assembly");
456+
}
457+
447458
//------------------------------------------------------------------------------
448459
// Create FDM element inverse not supported
449460
//------------------------------------------------------------------------------
@@ -473,6 +484,10 @@ int CeedOperatorCreate_Cuda(CeedOperator op) {
473484
ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearDiagonal",
474485
CeedOperatorAssembleLinearDiagonal_Cuda);
475486
CeedChk(ierr);
487+
ierr = CeedSetBackendFunction(ceed, "Operator", op,
488+
"AssembleLinearPointBlockDiagonal",
489+
CeedOperatorAssembleLinearPointBlockDiagonal_Cuda);
490+
CeedChk(ierr);
476491
ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
477492
CeedOperatorCreateFDMElementInverse_Cuda);
478493
CeedChk(ierr);

‎backends/occa/ceed-occa-operator.c

Lines changed: 56 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -519,6 +519,48 @@ static int CeedOperatorApply_Occa(CeedOperator op,
519519
return 0;
520520
}
521521

522+
// *****************************************************************************
523+
// * Assemble linear QFunction not supported
524+
// *****************************************************************************
525+
static int CeedOperatorAssembleLinearQFunction_Occa(CeedOperator op) {
526+
int ierr;
527+
Ceed ceed;
528+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
529+
return CeedError(ceed, 1, "Backend does not implement QFunction assembly");
530+
}
531+
532+
// *****************************************************************************
533+
// * Assemble linear diagonal not supported
534+
// *****************************************************************************
535+
static int CeedOperatorAssembleLinearDiagonal_Occa(CeedOperator op) {
536+
int ierr;
537+
Ceed ceed;
538+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
539+
return CeedError(ceed, 1,
540+
"Backend does not implement Operator diagonal assembly");
541+
}
542+
543+
// *****************************************************************************
544+
// * Assemble linear point block diagonal not supported
545+
// *****************************************************************************
546+
static int CeedOperatorAssembleLinearPointBlockDiagonal_Occa(CeedOperator op) {
547+
int ierr;
548+
Ceed ceed;
549+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
550+
return CeedError(ceed, 1,
551+
"Backend does not implement Operator point block diagonal assembly");
552+
}
553+
554+
// *****************************************************************************
555+
// * Create FDM element inverse not supported
556+
// *****************************************************************************
557+
static int CeedOperatorCreateFDMElementInverse_Occa(CeedOperator op) {
558+
int ierr;
559+
Ceed ceed;
560+
ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
561+
return CeedError(ceed, 1, "Backend does not implement FDM inverse creation");
562+
}
563+
522564
// *****************************************************************************
523565
// * Create an operator
524566
// *****************************************************************************
@@ -532,6 +574,20 @@ int CeedOperatorCreate_Occa(CeedOperator op) {
532574
ierr = CeedCalloc(1, &impl); CeedChk(ierr);
533575
ierr = CeedOperatorSetData(op, (void *)&impl); CeedChk(ierr);
534576

577+
ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearQFunction",
578+
CeedOperatorAssembleLinearQFunction_Occa);
579+
CeedChk(ierr);
580+
ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearDiagonal",
581+
CeedOperatorAssembleLinearDiagonal_Occa);
582+
CeedChk(ierr);
583+
ierr = CeedSetBackendFunction(ceed, "Operator", op,
584+
"AssembleLinearPointBlockDiagonal",
585+
CeedOperatorAssembleLinearPointBlockDiagonal_Occa);
586+
CeedChk(ierr);
587+
ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
588+
CeedOperatorCreateFDMElementInverse_Occa);
589+
CeedChk(ierr);
590+
535591
ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd",
536592
CeedOperatorApply_Occa); CeedChk(ierr);
537593
ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy",

‎backends/ref/ceed-ref-operator.c

Lines changed: 98 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -566,7 +566,7 @@ static int CeedOperatorAssembleLinearQFunction_Ref(CeedOperator op,
566566
// LCOV_EXCL_STOP
567567

568568
// Create output restriction
569-
CeedInt strides[3] = {1, Q, numactivein *numactiveout*Q};
569+
CeedInt strides[3] = {1, Q, numactivein*numactiveout*Q}; /* *NOPAD* */
570570
ierr = CeedElemRestrictionCreateStrided(ceedparent, numelements, Q,
571571
numactivein*numactiveout,
572572
numactivein*numactiveout*numelements*Q,
@@ -664,11 +664,51 @@ static inline void CeedOperatorGetBasisPointer_Ref(const CeedScalar **basisptr,
664664
}
665665
}
666666

667+
//------------------------------------------------------------------------------
668+
// Create point block restriction
669+
//------------------------------------------------------------------------------
670+
static int CreatePBRestriction(CeedElemRestriction rstr,
671+
CeedElemRestriction *pbRstr) {
672+
int ierr;
673+
Ceed ceed;
674+
ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChk(ierr);
675+
const CeedInt *offsets;
676+
ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets);
677+
CeedChk(ierr);
678+
679+
// Expand offsets
680+
CeedInt nelem, ncomp, elemsize, compstride, max = 1, *pbOffsets;
681+
ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChk(ierr);
682+
ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChk(ierr);
683+
ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChk(ierr);
684+
ierr = CeedElemRestrictionGetCompStride(rstr, &compstride); CeedChk(ierr);
685+
CeedInt shift = ncomp;
686+
if (compstride != 1)
687+
shift *= ncomp;
688+
ierr = CeedCalloc(nelem*elemsize, &pbOffsets); CeedChk(ierr);
689+
for (CeedInt i = 0; i < nelem*elemsize; i++) {
690+
pbOffsets[i] = offsets[i]*shift;
691+
if (pbOffsets[i] > max)
692+
max = pbOffsets[i];
693+
}
694+
695+
// Create new restriction
696+
ierr = CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp*ncomp, 1,
697+
max + ncomp*ncomp, CEED_MEM_HOST,
698+
CEED_OWN_POINTER, pbOffsets, pbRstr);
699+
CeedChk(ierr);
700+
701+
// Cleanup
702+
ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChk(ierr);
703+
704+
return 0;
705+
}
706+
667707
//------------------------------------------------------------------------------
668708
// Assemble Linear Diagonal
669709
//------------------------------------------------------------------------------
670-
static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
671-
CeedVector *assembled, CeedRequest *request) {
710+
static int CeedOperatorAssembleDiagonalCore_Ref(CeedOperator op,
711+
CeedVector *assembled, CeedRequest *request, const bool pointBlock) {
672712
int ierr;
673713

674714
// Assemble QFunction
@@ -764,9 +804,15 @@ static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
764804
}
765805
}
766806

807+
// Assemble point-block diagonal restriction, if needed
808+
CeedElemRestriction diagrstr = rstrout;
809+
if (pointBlock) {
810+
ierr = CreatePBRestriction(rstrout, &diagrstr); CeedChk(ierr);
811+
}
812+
767813
// Create diagonal vector
768814
CeedVector elemdiag;
769-
ierr = CeedElemRestrictionCreateVector(rstrin, assembled, &elemdiag);
815+
ierr = CeedElemRestrictionCreateVector(diagrstr, assembled, &elemdiag);
770816
CeedChk(ierr);
771817

772818
// Assemble element operator diagonals
@@ -777,7 +823,7 @@ static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
777823
ierr = CeedVectorGetArray(assembledqf, CEED_MEM_HOST, &assembledqfarray);
778824
CeedChk(ierr);
779825
CeedInt nelem, nnodes, nqpts;
780-
ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); CeedChk(ierr);
826+
ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); CeedChk(ierr);
781827
ierr = CeedBasisGetNumNodes(basisin, &nnodes); CeedChk(ierr);
782828
ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChk(ierr);
783829
// Basis matrices
@@ -816,17 +862,30 @@ static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
816862
CeedOperatorGetBasisPointer_Ref(&b, emodein[ein], identity, interpin,
817863
&gradin[din*nqpts*nnodes]);
818864
// Each component
819-
for (CeedInt comp=0; comp<ncomp; comp++)
865+
for (CeedInt compOut=0; compOut<ncomp; compOut++)
820866
// Each qpoint/node pair
821-
for (CeedInt q=0; q<nqpts; q++) {
822-
const CeedScalar qfvalue =
823-
assembledqfarray[((((e*numemodein+ein)*ncomp+comp)*
824-
numemodeout+eout)*ncomp+comp)*nqpts+q];
825-
if (fabs(qfvalue) > maxnorm*1e-12)
826-
for (CeedInt n=0; n<nnodes; n++)
827-
elemdiagarray[(e*ncomp+comp)*nnodes+n] += bt[q*nnodes+n] *
828-
qfvalue * b[q*nnodes+n];
829-
}
867+
for (CeedInt q=0; q<nqpts; q++)
868+
if (pointBlock) {
869+
// Point Block Diagonal
870+
for (CeedInt compIn=0; compIn<ncomp; compIn++) {
871+
const CeedScalar qfvalue =
872+
assembledqfarray[((((e*numemodein+ein)*ncomp+compIn)*
873+
numemodeout+eout)*ncomp+compOut)*nqpts+q];
874+
if (fabs(qfvalue) > maxnorm*1e-12)
875+
for (CeedInt n=0; n<nnodes; n++)
876+
elemdiagarray[((e*ncomp+compOut)*ncomp+compIn)*nnodes+n] +=
877+
bt[q*nnodes+n] * qfvalue * b[q*nnodes+n];
878+
}
879+
} else {
880+
// Diagonal Only
881+
const CeedScalar qfvalue =
882+
assembledqfarray[((((e*numemodein+ein)*ncomp+compOut)*
883+
numemodeout+eout)*ncomp+compOut)*nqpts+q];
884+
if (fabs(qfvalue) > maxnorm*1e-12)
885+
for (CeedInt n=0; n<nnodes; n++)
886+
elemdiagarray[(e*ncomp+compOut)*nnodes+n] +=
887+
bt[q*nnodes+n] * qfvalue * b[q*nnodes+n];
888+
}
830889
}
831890
}
832891
}
@@ -835,10 +894,13 @@ static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
835894

836895
// Assemble local operator diagonal
837896
ierr = CeedVectorSetValue(*assembled, 0.0); CeedChk(ierr);
838-
ierr = CeedElemRestrictionApply(rstrout, CEED_TRANSPOSE, elemdiag,
897+
ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag,
839898
*assembled, request); CeedChk(ierr);
840899

841900
// Cleanup
901+
if (pointBlock) {
902+
ierr = CeedElemRestrictionDestroy(&diagrstr); CeedChk(ierr);
903+
}
842904
ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr);
843905
ierr = CeedVectorDestroy(&elemdiag); CeedChk(ierr);
844906
ierr = CeedFree(&emodein); CeedChk(ierr);
@@ -848,6 +910,22 @@ static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
848910
return 0;
849911
}
850912

913+
//------------------------------------------------------------------------------
914+
// Assemble Linear Diagonal
915+
//------------------------------------------------------------------------------
916+
static int CeedOperatorAssembleLinearDiagonal_Ref(CeedOperator op,
917+
CeedVector *assembled, CeedRequest *request) {
918+
return CeedOperatorAssembleDiagonalCore_Ref(op, assembled, request, false);
919+
}
920+
921+
//------------------------------------------------------------------------------
922+
// Assemble Linear Point Block Diagonal
923+
//------------------------------------------------------------------------------
924+
static int CeedOperatorAssembleLinearPointBlockDiagonal_Ref(CeedOperator op,
925+
CeedVector *assembled, CeedRequest *request) {
926+
return CeedOperatorAssembleDiagonalCore_Ref(op, assembled, request, true);
927+
}
928+
851929
//------------------------------------------------------------------------------
852930
// Create FDM Element Inverse
853931
//------------------------------------------------------------------------------
@@ -1101,6 +1179,10 @@ int CeedOperatorCreate_Ref(CeedOperator op) {
11011179
ierr = CeedSetBackendFunction(ceed, "Operator", op, "AssembleLinearDiagonal",
11021180
CeedOperatorAssembleLinearDiagonal_Ref);
11031181
CeedChk(ierr);
1182+
ierr = CeedSetBackendFunction(ceed, "Operator", op,
1183+
"AssembleLinearPointBlockDiagonal",
1184+
CeedOperatorAssembleLinearPointBlockDiagonal_Ref);
1185+
CeedChk(ierr);
11041186
ierr = CeedSetBackendFunction(ceed, "Operator", op, "CreateFDMElementInverse",
11051187
CeedOperatorCreateFDMElementInverse_Ref);
11061188
CeedChk(ierr);

‎include/ceed-backend.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ CEED_EXTERN int CeedElemRestrictionGetOffsets(CeedElemRestriction rstr,
9191
CeedMemType mtype, const CeedInt **offsets);
9292
CEED_EXTERN int CeedElemRestrictionRestoreOffsets(CeedElemRestriction rstr,
9393
const CeedInt **offsets);
94-
CEED_EXTERN int CeedElemRestrictionGetStridedStatus( CeedElemRestriction rstr,
94+
CEED_EXTERN int CeedElemRestrictionGetStridedStatus(CeedElemRestriction rstr,
9595
bool *status);
9696
CEED_EXTERN int CeedElemRestrictionGetBackendStridesStatus(
9797
CeedElemRestriction rstr, bool *status);

‎include/ceed-impl.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -271,6 +271,8 @@ struct CeedOperator_private {
271271
int (*AssembleLinearQFunction)(CeedOperator, CeedVector *,
272272
CeedElemRestriction *, CeedRequest *);
273273
int (*AssembleLinearDiagonal)(CeedOperator, CeedVector *, CeedRequest *);
274+
int (*AssembleLinearPointBlockDiagonal)(CeedOperator, CeedVector *,
275+
CeedRequest *);
274276
int (*CreateFDMElementInverse)(CeedOperator, CeedOperator *, CeedRequest *);
275277
int (*Apply)(CeedOperator, CeedVector, CeedVector, CeedRequest *);
276278
int (*ApplyComposite)(CeedOperator, CeedVector, CeedVector, CeedRequest *);

‎include/ceed.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -488,6 +488,8 @@ CEED_EXTERN int CeedOperatorAssembleLinearQFunction(CeedOperator op,
488488
CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request);
489489
CEED_EXTERN int CeedOperatorAssembleLinearDiagonal(CeedOperator op,
490490
CeedVector *assembled, CeedRequest *request);
491+
CEED_EXTERN int CeedOperatorAssembleLinearPointBlockDiagonal(CeedOperator op,
492+
CeedVector *assembled, CeedRequest *request);
491493
CEED_EXTERN int CeedOperatorCreateFDMElementInverse(CeedOperator op,
492494
CeedOperator *fdminv, CeedRequest *request);
493495
CEED_EXTERN int CeedOperatorView(CeedOperator op, FILE *stream);

0 commit comments

Comments
 (0)