1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed/ceed.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-tools.h> 11 #include <hip/hip_runtime.h> 12 #include <assert.h> 13 #include <stdbool.h> 14 #include <string.h> 15 #include "ceed-hip-ref.h" 16 #include "../hip/ceed-hip-compile.h" 17 18 //------------------------------------------------------------------------------ 19 // Destroy operator 20 //------------------------------------------------------------------------------ 21 static int CeedOperatorDestroy_Hip(CeedOperator op) { 22 int ierr; 23 CeedOperator_Hip *impl; 24 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 25 26 // Apply data 27 for (CeedInt i = 0; i < impl->numein + impl->numeout; i++) { 28 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChkBackend(ierr); 29 } 30 ierr = CeedFree(&impl->evecs); CeedChkBackend(ierr); 31 32 for (CeedInt i = 0; i < impl->numein; i++) { 33 ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChkBackend(ierr); 34 } 35 ierr = CeedFree(&impl->qvecsin); CeedChkBackend(ierr); 36 37 for (CeedInt i = 0; i < impl->numeout; i++) { 38 ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChkBackend(ierr); 39 } 40 ierr = CeedFree(&impl->qvecsout); CeedChkBackend(ierr); 41 42 // QFunction diagonal assembly data 43 for (CeedInt i=0; i<impl->qfnumactivein; i++) { 44 ierr = CeedVectorDestroy(&impl->qfactivein[i]); CeedChkBackend(ierr); 45 } 46 ierr = CeedFree(&impl->qfactivein); CeedChkBackend(ierr); 47 48 // Diag data 49 if (impl->diag) { 50 Ceed ceed; 51 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 52 CeedChk_Hip(ceed, hipModuleUnload(impl->diag->module)); 53 ierr = CeedFree(&impl->diag->h_emodein); CeedChkBackend(ierr); 54 ierr = CeedFree(&impl->diag->h_emodeout); CeedChkBackend(ierr); 55 ierr = hipFree(impl->diag->d_emodein); CeedChk_Hip(ceed, ierr); 56 ierr = hipFree(impl->diag->d_emodeout); CeedChk_Hip(ceed, ierr); 57 ierr = hipFree(impl->diag->d_identity); CeedChk_Hip(ceed, ierr); 58 ierr = hipFree(impl->diag->d_interpin); CeedChk_Hip(ceed, ierr); 59 ierr = hipFree(impl->diag->d_interpout); CeedChk_Hip(ceed, ierr); 60 ierr = hipFree(impl->diag->d_gradin); CeedChk_Hip(ceed, ierr); 61 ierr = hipFree(impl->diag->d_gradout); CeedChk_Hip(ceed, ierr); 62 ierr = CeedElemRestrictionDestroy(&impl->diag->pbdiagrstr); 63 CeedChkBackend(ierr); 64 ierr = CeedVectorDestroy(&impl->diag->elemdiag); CeedChkBackend(ierr); 65 ierr = CeedVectorDestroy(&impl->diag->pbelemdiag); CeedChkBackend(ierr); 66 } 67 ierr = CeedFree(&impl->diag); CeedChkBackend(ierr); 68 69 if (impl->asmb) { 70 Ceed ceed; 71 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 72 CeedChk_Hip(ceed, hipModuleUnload(impl->asmb->module)); 73 ierr = hipFree(impl->asmb->d_B_in); CeedChk_Hip(ceed, ierr); 74 ierr = hipFree(impl->asmb->d_B_out); CeedChk_Hip(ceed, ierr); 75 } 76 ierr = CeedFree(&impl->asmb); CeedChkBackend(ierr); 77 78 ierr = CeedFree(&impl); CeedChkBackend(ierr); 79 return CEED_ERROR_SUCCESS; 80 } 81 82 //------------------------------------------------------------------------------ 83 // Setup infields or outfields 84 //------------------------------------------------------------------------------ 85 static int CeedOperatorSetupFields_Hip(CeedQFunction qf, CeedOperator op, 86 bool isinput, CeedVector *evecs, 87 CeedVector *qvecs, CeedInt starte, 88 CeedInt numfields, CeedInt Q, 89 CeedInt numelements) { 90 CeedInt dim, ierr, size; 91 CeedSize q_size; 92 Ceed ceed; 93 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 94 CeedBasis basis; 95 CeedElemRestriction Erestrict; 96 CeedOperatorField *opfields; 97 CeedQFunctionField *qffields; 98 CeedVector fieldvec; 99 bool strided; 100 bool skiprestrict; 101 102 if (isinput) { 103 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 104 CeedChkBackend(ierr); 105 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 106 CeedChkBackend(ierr); 107 } else { 108 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 109 CeedChkBackend(ierr); 110 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 111 CeedChkBackend(ierr); 112 } 113 114 // Loop over fields 115 for (CeedInt i = 0; i < numfields; i++) { 116 CeedEvalMode emode; 117 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 118 119 strided = false; 120 skiprestrict = false; 121 if (emode != CEED_EVAL_WEIGHT) { 122 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &Erestrict); 123 CeedChkBackend(ierr); 124 125 // Check whether this field can skip the element restriction: 126 // must be passive input, with emode NONE, and have a strided restriction with 127 // CEED_STRIDES_BACKEND. 128 129 // First, check whether the field is input or output: 130 if (isinput) { 131 // Check for passive input: 132 ierr = CeedOperatorFieldGetVector(opfields[i], &fieldvec); CeedChkBackend(ierr); 133 if (fieldvec != CEED_VECTOR_ACTIVE) { 134 // Check emode 135 if (emode == CEED_EVAL_NONE) { 136 // Check for strided restriction 137 ierr = CeedElemRestrictionIsStrided(Erestrict, &strided); 138 CeedChkBackend(ierr); 139 if (strided) { 140 // Check if vector is already in preferred backend ordering 141 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, 142 &skiprestrict); CeedChkBackend(ierr); 143 } 144 } 145 } 146 } 147 if (skiprestrict) { 148 // We do not need an E-Vector, but will use the input field vector's data 149 // directly in the operator application. 150 evecs[i + starte] = NULL; 151 } else { 152 ierr = CeedElemRestrictionCreateVector(Erestrict, NULL, 153 &evecs[i + starte]); 154 CeedChkBackend(ierr); 155 } 156 } 157 158 switch (emode) { 159 case CEED_EVAL_NONE: 160 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 161 q_size = (CeedSize)numelements * Q * size; 162 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 163 break; 164 case CEED_EVAL_INTERP: 165 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 166 q_size = (CeedSize)numelements * Q * size; 167 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 168 break; 169 case CEED_EVAL_GRAD: 170 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 171 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 172 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 173 q_size = (CeedSize)numelements * Q * size; 174 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 175 break; 176 case CEED_EVAL_WEIGHT: // Only on input fields 177 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 178 q_size = (CeedSize)numelements * Q; 179 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 180 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 181 CEED_EVAL_WEIGHT, NULL, qvecs[i]); CeedChkBackend(ierr); 182 break; 183 case CEED_EVAL_DIV: 184 break; // TODO: Not implemented 185 case CEED_EVAL_CURL: 186 break; // TODO: Not implemented 187 } 188 } 189 return CEED_ERROR_SUCCESS; 190 } 191 192 //------------------------------------------------------------------------------ 193 // CeedOperator needs to connect all the named fields (be they active or passive) 194 // to the named inputs and outputs of its CeedQFunction. 195 //------------------------------------------------------------------------------ 196 static int CeedOperatorSetup_Hip(CeedOperator op) { 197 int ierr; 198 bool setupdone; 199 ierr = CeedOperatorIsSetupDone(op, &setupdone); CeedChkBackend(ierr); 200 if (setupdone) 201 return CEED_ERROR_SUCCESS; 202 Ceed ceed; 203 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 204 CeedOperator_Hip *impl; 205 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 206 CeedQFunction qf; 207 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 208 CeedInt Q, numelements, numinputfields, numoutputfields; 209 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 210 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 211 CeedOperatorField *opinputfields, *opoutputfields; 212 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 213 &numoutputfields, &opoutputfields); 214 CeedChkBackend(ierr); 215 CeedQFunctionField *qfinputfields, *qfoutputfields; 216 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 217 CeedChkBackend(ierr); 218 219 // Allocate 220 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs); 221 CeedChkBackend(ierr); 222 223 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsin); CeedChkBackend(ierr); 224 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsout); CeedChkBackend(ierr); 225 226 impl->numein = numinputfields; impl->numeout = numoutputfields; 227 228 // Set up infield and outfield evecs and qvecs 229 // Infields 230 ierr = CeedOperatorSetupFields_Hip(qf, op, true, 231 impl->evecs, impl->qvecsin, 0, 232 numinputfields, Q, numelements); 233 CeedChkBackend(ierr); 234 235 // Outfields 236 ierr = CeedOperatorSetupFields_Hip(qf, op, false, 237 impl->evecs, impl->qvecsout, 238 numinputfields, numoutputfields, Q, 239 numelements); CeedChkBackend(ierr); 240 241 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 242 return CEED_ERROR_SUCCESS; 243 } 244 245 //------------------------------------------------------------------------------ 246 // Setup Operator Inputs 247 //------------------------------------------------------------------------------ 248 static inline int CeedOperatorSetupInputs_Hip(CeedInt numinputfields, 249 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 250 CeedVector invec, const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 251 CeedOperator_Hip *impl, CeedRequest *request) { 252 CeedInt ierr; 253 CeedEvalMode emode; 254 CeedVector vec; 255 CeedElemRestriction Erestrict; 256 257 for (CeedInt i = 0; i < numinputfields; i++) { 258 // Get input vector 259 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 260 if (vec == CEED_VECTOR_ACTIVE) { 261 if (skipactive) 262 continue; 263 else 264 vec = invec; 265 } 266 267 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 268 CeedChkBackend(ierr); 269 if (emode == CEED_EVAL_WEIGHT) { // Skip 270 } else { 271 // Get input vector 272 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 273 // Get input element restriction 274 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 275 CeedChkBackend(ierr); 276 if (vec == CEED_VECTOR_ACTIVE) 277 vec = invec; 278 // Restrict, if necessary 279 if (!impl->evecs[i]) { 280 // No restriction for this field; read data directly from vec. 281 ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, 282 (const CeedScalar **) &edata[i]); 283 CeedChkBackend(ierr); 284 } else { 285 ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE, vec, 286 impl->evecs[i], request); CeedChkBackend(ierr); 287 // Get evec 288 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_DEVICE, 289 (const CeedScalar **) &edata[i]); 290 CeedChkBackend(ierr); 291 } 292 } 293 } 294 return CEED_ERROR_SUCCESS; 295 } 296 297 //------------------------------------------------------------------------------ 298 // Input Basis Action 299 //------------------------------------------------------------------------------ 300 static inline int CeedOperatorInputBasis_Hip(CeedInt numelements, 301 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 302 CeedInt numinputfields, const bool skipactive, 303 CeedScalar *edata[2*CEED_FIELD_MAX], CeedOperator_Hip *impl) { 304 CeedInt ierr; 305 CeedInt elemsize, size; 306 CeedElemRestriction Erestrict; 307 CeedEvalMode emode; 308 CeedBasis basis; 309 310 for (CeedInt i=0; i<numinputfields; i++) { 311 // Skip active input 312 if (skipactive) { 313 CeedVector vec; 314 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 315 if (vec == CEED_VECTOR_ACTIVE) 316 continue; 317 } 318 // Get elemsize, emode, size 319 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 320 CeedChkBackend(ierr); 321 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 322 CeedChkBackend(ierr); 323 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 324 CeedChkBackend(ierr); 325 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 326 // Basis action 327 switch (emode) { 328 case CEED_EVAL_NONE: 329 ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_DEVICE, 330 CEED_USE_POINTER, edata[i]); CeedChkBackend(ierr); 331 break; 332 case CEED_EVAL_INTERP: 333 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 334 CeedChkBackend(ierr); 335 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 336 CEED_EVAL_INTERP, impl->evecs[i], 337 impl->qvecsin[i]); CeedChkBackend(ierr); 338 break; 339 case CEED_EVAL_GRAD: 340 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 341 CeedChkBackend(ierr); 342 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 343 CEED_EVAL_GRAD, impl->evecs[i], 344 impl->qvecsin[i]); CeedChkBackend(ierr); 345 break; 346 case CEED_EVAL_WEIGHT: 347 break; // No action 348 case CEED_EVAL_DIV: 349 break; // TODO: Not implemented 350 case CEED_EVAL_CURL: 351 break; // TODO: Not implemented 352 } 353 } 354 return CEED_ERROR_SUCCESS; 355 } 356 357 //------------------------------------------------------------------------------ 358 // Restore Input Vectors 359 //------------------------------------------------------------------------------ 360 static inline int CeedOperatorRestoreInputs_Hip(CeedInt numinputfields, 361 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 362 const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 363 CeedOperator_Hip *impl) { 364 CeedInt ierr; 365 CeedEvalMode emode; 366 CeedVector vec; 367 368 for (CeedInt i = 0; i < numinputfields; i++) { 369 // Skip active input 370 if (skipactive) { 371 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 372 if (vec == CEED_VECTOR_ACTIVE) 373 continue; 374 } 375 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 376 CeedChkBackend(ierr); 377 if (emode == CEED_EVAL_WEIGHT) { // Skip 378 } else { 379 if (!impl->evecs[i]) { // This was a skiprestrict case 380 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 381 ierr = CeedVectorRestoreArrayRead(vec, 382 (const CeedScalar **)&edata[i]); 383 CeedChkBackend(ierr); 384 } else { 385 ierr = CeedVectorRestoreArrayRead(impl->evecs[i], 386 (const CeedScalar **) &edata[i]); 387 CeedChkBackend(ierr); 388 } 389 } 390 } 391 return CEED_ERROR_SUCCESS; 392 } 393 394 //------------------------------------------------------------------------------ 395 // Apply and add to output 396 //------------------------------------------------------------------------------ 397 static int CeedOperatorApplyAdd_Hip(CeedOperator op, CeedVector invec, 398 CeedVector outvec, CeedRequest *request) { 399 int ierr; 400 CeedOperator_Hip *impl; 401 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 402 CeedQFunction qf; 403 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 404 CeedInt Q, numelements, elemsize, numinputfields, numoutputfields, size; 405 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 406 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 407 CeedOperatorField *opinputfields, *opoutputfields; 408 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 409 &numoutputfields, &opoutputfields); 410 CeedChkBackend(ierr); 411 CeedQFunctionField *qfinputfields, *qfoutputfields; 412 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 413 CeedChkBackend(ierr); 414 CeedEvalMode emode; 415 CeedVector vec; 416 CeedBasis basis; 417 CeedElemRestriction Erestrict; 418 CeedScalar *edata[2*CEED_FIELD_MAX]; 419 420 // Setup 421 ierr = CeedOperatorSetup_Hip(op); CeedChkBackend(ierr); 422 423 // Input Evecs and Restriction 424 ierr = CeedOperatorSetupInputs_Hip(numinputfields, qfinputfields, 425 opinputfields, invec, false, edata, 426 impl, request); CeedChkBackend(ierr); 427 428 // Input basis apply if needed 429 ierr = CeedOperatorInputBasis_Hip(numelements, qfinputfields, opinputfields, 430 numinputfields, false, edata, impl); 431 CeedChkBackend(ierr); 432 433 // Output pointers, as necessary 434 for (CeedInt i = 0; i < numoutputfields; i++) { 435 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 436 CeedChkBackend(ierr); 437 if (emode == CEED_EVAL_NONE) { 438 // Set the output Q-Vector to use the E-Vector data directly. 439 ierr = CeedVectorGetArrayWrite(impl->evecs[i + impl->numein], CEED_MEM_DEVICE, 440 &edata[i + numinputfields]); CeedChkBackend(ierr); 441 ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_DEVICE, 442 CEED_USE_POINTER, edata[i + numinputfields]); 443 CeedChkBackend(ierr); 444 } 445 } 446 447 // Q function 448 ierr = CeedQFunctionApply(qf, numelements * Q, impl->qvecsin, impl->qvecsout); 449 CeedChkBackend(ierr); 450 451 // Output basis apply if needed 452 for (CeedInt i = 0; i < numoutputfields; i++) { 453 // Get elemsize, emode, size 454 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 455 CeedChkBackend(ierr); 456 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 457 CeedChkBackend(ierr); 458 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 459 CeedChkBackend(ierr); 460 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 461 CeedChkBackend(ierr); 462 // Basis action 463 switch (emode) { 464 case CEED_EVAL_NONE: 465 break; 466 case CEED_EVAL_INTERP: 467 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 468 CeedChkBackend(ierr); 469 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 470 CEED_EVAL_INTERP, impl->qvecsout[i], 471 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 472 break; 473 case CEED_EVAL_GRAD: 474 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 475 CeedChkBackend(ierr); 476 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 477 CEED_EVAL_GRAD, impl->qvecsout[i], 478 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 479 break; 480 // LCOV_EXCL_START 481 case CEED_EVAL_WEIGHT: { 482 Ceed ceed; 483 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 484 return CeedError(ceed, CEED_ERROR_BACKEND, 485 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 486 break; // Should not occur 487 } 488 case CEED_EVAL_DIV: 489 break; // TODO: Not implemented 490 case CEED_EVAL_CURL: 491 break; // TODO: Not implemented 492 // LCOV_EXCL_STOP 493 } 494 } 495 496 // Output restriction 497 for (CeedInt i = 0; i < numoutputfields; i++) { 498 // Restore evec 499 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 500 CeedChkBackend(ierr); 501 if (emode == CEED_EVAL_NONE) { 502 ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein], 503 &edata[i + numinputfields]); 504 CeedChkBackend(ierr); 505 } 506 // Get output vector 507 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 508 CeedChkBackend(ierr); 509 // Restrict 510 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 511 CeedChkBackend(ierr); 512 // Active 513 if (vec == CEED_VECTOR_ACTIVE) 514 vec = outvec; 515 516 ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE, 517 impl->evecs[i + impl->numein], vec, 518 request); CeedChkBackend(ierr); 519 } 520 521 // Restore input arrays 522 ierr = CeedOperatorRestoreInputs_Hip(numinputfields, qfinputfields, 523 opinputfields, false, edata, impl); 524 CeedChkBackend(ierr); 525 return CEED_ERROR_SUCCESS; 526 } 527 528 //------------------------------------------------------------------------------ 529 // Core code for assembling linear QFunction 530 //------------------------------------------------------------------------------ 531 static inline int CeedOperatorLinearAssembleQFunctionCore_Hip(CeedOperator op, 532 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 533 CeedRequest *request) { 534 int ierr; 535 CeedOperator_Hip *impl; 536 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 537 CeedQFunction qf; 538 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 539 CeedInt Q, numelements, numinputfields, numoutputfields, size; 540 CeedSize q_size; 541 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 542 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 543 CeedOperatorField *opinputfields, *opoutputfields; 544 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 545 &numoutputfields, &opoutputfields); 546 CeedChkBackend(ierr); 547 CeedQFunctionField *qfinputfields, *qfoutputfields; 548 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 549 CeedChkBackend(ierr); 550 CeedVector vec; 551 CeedInt numactivein = impl->qfnumactivein, numactiveout = impl->qfnumactiveout; 552 CeedVector *activein = impl->qfactivein; 553 CeedScalar *a, *tmp; 554 Ceed ceed, ceedparent; 555 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 556 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceedparent); 557 CeedChkBackend(ierr); 558 ceedparent = ceedparent ? ceedparent : ceed; 559 CeedScalar *edata[2*CEED_FIELD_MAX]; 560 561 // Setup 562 ierr = CeedOperatorSetup_Hip(op); CeedChkBackend(ierr); 563 564 // Check for identity 565 bool identityqf; 566 ierr = CeedQFunctionIsIdentity(qf, &identityqf); CeedChkBackend(ierr); 567 if (identityqf) 568 // LCOV_EXCL_START 569 return CeedError(ceed, CEED_ERROR_BACKEND, 570 "Assembling identity QFunctions not supported"); 571 // LCOV_EXCL_STOP 572 573 // Input Evecs and Restriction 574 ierr = CeedOperatorSetupInputs_Hip(numinputfields, qfinputfields, 575 opinputfields, NULL, true, edata, 576 impl, request); CeedChkBackend(ierr); 577 578 // Count number of active input fields 579 if (!numactivein) { 580 for (CeedInt i=0; i<numinputfields; i++) { 581 // Get input vector 582 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 583 // Check if active input 584 if (vec == CEED_VECTOR_ACTIVE) { 585 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 586 ierr = CeedVectorSetValue(impl->qvecsin[i], 0.0); CeedChkBackend(ierr); 587 ierr = CeedVectorGetArray(impl->qvecsin[i], CEED_MEM_DEVICE, &tmp); 588 CeedChkBackend(ierr); 589 ierr = CeedRealloc(numactivein + size, &activein); CeedChkBackend(ierr); 590 for (CeedInt field = 0; field < size; field++) { 591 q_size = (CeedSize)Q*numelements; 592 ierr = CeedVectorCreate(ceed, q_size, &activein[numactivein+field]); 593 CeedChkBackend(ierr); 594 ierr = CeedVectorSetArray(activein[numactivein+field], CEED_MEM_DEVICE, 595 CEED_USE_POINTER, &tmp[field*Q*numelements]); 596 CeedChkBackend(ierr); 597 } 598 numactivein += size; 599 ierr = CeedVectorRestoreArray(impl->qvecsin[i], &tmp); CeedChkBackend(ierr); 600 } 601 } 602 impl->qfnumactivein = numactivein; 603 impl->qfactivein = activein; 604 } 605 606 // Count number of active output fields 607 if (!numactiveout) { 608 for (CeedInt i=0; i<numoutputfields; i++) { 609 // Get output vector 610 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 611 CeedChkBackend(ierr); 612 // Check if active output 613 if (vec == CEED_VECTOR_ACTIVE) { 614 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 615 CeedChkBackend(ierr); 616 numactiveout += size; 617 } 618 } 619 impl->qfnumactiveout = numactiveout; 620 } 621 622 // Check sizes 623 if (!numactivein || !numactiveout) 624 // LCOV_EXCL_START 625 return CeedError(ceed, CEED_ERROR_BACKEND, 626 "Cannot assemble QFunction without active inputs " 627 "and outputs"); 628 // LCOV_EXCL_STOP 629 630 // Build objects if needed 631 if (build_objects) { 632 // Create output restriction 633 CeedInt strides[3] = {1, numelements*Q, Q}; /* *NOPAD* */ 634 ierr = CeedElemRestrictionCreateStrided(ceedparent, numelements, Q, 635 numactivein*numactiveout, 636 numactivein*numactiveout*numelements*Q, 637 strides, rstr); CeedChkBackend(ierr); 638 // Create assembled vector 639 CeedSize l_size = (CeedSize)numelements*Q*numactivein*numactiveout; 640 ierr = CeedVectorCreate(ceedparent, l_size, assembled); CeedChkBackend(ierr); 641 } 642 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 643 ierr = CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &a); 644 CeedChkBackend(ierr); 645 646 // Input basis apply 647 ierr = CeedOperatorInputBasis_Hip(numelements, qfinputfields, opinputfields, 648 numinputfields, true, edata, impl); 649 CeedChkBackend(ierr); 650 651 // Assemble QFunction 652 for (CeedInt in=0; in<numactivein; in++) { 653 // Set Inputs 654 ierr = CeedVectorSetValue(activein[in], 1.0); CeedChkBackend(ierr); 655 if (numactivein > 1) { 656 ierr = CeedVectorSetValue(activein[(in+numactivein-1)%numactivein], 657 0.0); CeedChkBackend(ierr); 658 } 659 // Set Outputs 660 for (CeedInt out=0; out<numoutputfields; out++) { 661 // Get output vector 662 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 663 CeedChkBackend(ierr); 664 // Check if active output 665 if (vec == CEED_VECTOR_ACTIVE) { 666 CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_DEVICE, 667 CEED_USE_POINTER, a); CeedChkBackend(ierr); 668 ierr = CeedQFunctionFieldGetSize(qfoutputfields[out], &size); 669 CeedChkBackend(ierr); 670 a += size*Q*numelements; // Advance the pointer by the size of the output 671 } 672 } 673 // Apply QFunction 674 ierr = CeedQFunctionApply(qf, Q*numelements, impl->qvecsin, impl->qvecsout); 675 CeedChkBackend(ierr); 676 } 677 678 // Un-set output Qvecs to prevent accidental overwrite of Assembled 679 for (CeedInt out=0; out<numoutputfields; out++) { 680 // Get output vector 681 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 682 CeedChkBackend(ierr); 683 // Check if active output 684 if (vec == CEED_VECTOR_ACTIVE) { 685 ierr = CeedVectorTakeArray(impl->qvecsout[out], CEED_MEM_DEVICE, NULL); 686 CeedChkBackend(ierr); 687 } 688 } 689 690 // Restore input arrays 691 ierr = CeedOperatorRestoreInputs_Hip(numinputfields, qfinputfields, 692 opinputfields, true, edata, impl); 693 CeedChkBackend(ierr); 694 695 // Restore output 696 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 697 698 return CEED_ERROR_SUCCESS; 699 } 700 701 //------------------------------------------------------------------------------ 702 // Assemble Linear QFunction 703 //------------------------------------------------------------------------------ 704 static int CeedOperatorLinearAssembleQFunction_Hip(CeedOperator op, 705 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 706 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, true, assembled, rstr, 707 request); 708 } 709 710 //------------------------------------------------------------------------------ 711 // Assemble Linear QFunction 712 //------------------------------------------------------------------------------ 713 static int CeedOperatorLinearAssembleQFunctionUpdate_Hip(CeedOperator op, 714 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 715 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, false, &assembled, &rstr, 716 request); 717 } 718 719 //------------------------------------------------------------------------------ 720 // Create point block restriction 721 //------------------------------------------------------------------------------ 722 static int CreatePBRestriction(CeedElemRestriction rstr, 723 CeedElemRestriction *pbRstr) { 724 int ierr; 725 Ceed ceed; 726 ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr); 727 const CeedInt *offsets; 728 ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets); 729 CeedChkBackend(ierr); 730 731 // Expand offsets 732 CeedInt nelem, ncomp, elemsize, compstride, max = 1, *pbOffsets; 733 ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChkBackend(ierr); 734 ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChkBackend(ierr); 735 ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChkBackend(ierr); 736 ierr = CeedElemRestrictionGetCompStride(rstr, &compstride); 737 CeedChkBackend(ierr); 738 CeedInt shift = ncomp; 739 if (compstride != 1) 740 shift *= ncomp; 741 ierr = CeedCalloc(nelem*elemsize, &pbOffsets); CeedChkBackend(ierr); 742 for (CeedInt i = 0; i < nelem*elemsize; i++) { 743 pbOffsets[i] = offsets[i]*shift; 744 if (pbOffsets[i] > max) 745 max = pbOffsets[i]; 746 } 747 748 // Create new restriction 749 ierr = CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp*ncomp, 1, 750 max + ncomp*ncomp, CEED_MEM_HOST, 751 CEED_OWN_POINTER, pbOffsets, pbRstr); 752 CeedChkBackend(ierr); 753 754 // Cleanup 755 ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChkBackend(ierr); 756 757 return CEED_ERROR_SUCCESS; 758 } 759 760 //------------------------------------------------------------------------------ 761 // Assemble diagonal setup 762 //------------------------------------------------------------------------------ 763 static inline int CeedOperatorAssembleDiagonalSetup_Hip(CeedOperator op, 764 const bool pointBlock) { 765 int ierr; 766 Ceed ceed; 767 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 768 CeedQFunction qf; 769 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 770 CeedInt numinputfields, numoutputfields; 771 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 772 CeedChkBackend(ierr); 773 774 // Determine active input basis 775 CeedOperatorField *opfields; 776 CeedQFunctionField *qffields; 777 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 778 CeedChkBackend(ierr); 779 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 780 CeedChkBackend(ierr); 781 CeedInt numemodein = 0, ncomp = 0, dim = 1; 782 CeedEvalMode *emodein = NULL; 783 CeedBasis basisin = NULL; 784 CeedElemRestriction rstrin = NULL; 785 for (CeedInt i = 0; i < numinputfields; i++) { 786 CeedVector vec; 787 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 788 if (vec == CEED_VECTOR_ACTIVE) { 789 CeedElemRestriction rstr; 790 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisin); CeedChkBackend(ierr); 791 ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChkBackend(ierr); 792 ierr = CeedBasisGetDimension(basisin, &dim); CeedChkBackend(ierr); 793 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 794 CeedChkBackend(ierr); 795 if (rstrin && rstrin != rstr) 796 // LCOV_EXCL_START 797 return CeedError(ceed, CEED_ERROR_BACKEND, 798 "Multi-field non-composite operator diagonal assembly not supported"); 799 // LCOV_EXCL_STOP 800 rstrin = rstr; 801 CeedEvalMode emode; 802 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); 803 CeedChkBackend(ierr); 804 switch (emode) { 805 case CEED_EVAL_NONE: 806 case CEED_EVAL_INTERP: 807 ierr = CeedRealloc(numemodein + 1, &emodein); CeedChkBackend(ierr); 808 emodein[numemodein] = emode; 809 numemodein += 1; 810 break; 811 case CEED_EVAL_GRAD: 812 ierr = CeedRealloc(numemodein + dim, &emodein); CeedChkBackend(ierr); 813 for (CeedInt d = 0; d < dim; d++) 814 emodein[numemodein+d] = emode; 815 numemodein += dim; 816 break; 817 case CEED_EVAL_WEIGHT: 818 case CEED_EVAL_DIV: 819 case CEED_EVAL_CURL: 820 break; // Caught by QF Assembly 821 } 822 } 823 } 824 825 // Determine active output basis 826 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 827 CeedChkBackend(ierr); 828 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 829 CeedChkBackend(ierr); 830 CeedInt numemodeout = 0; 831 CeedEvalMode *emodeout = NULL; 832 CeedBasis basisout = NULL; 833 CeedElemRestriction rstrout = NULL; 834 for (CeedInt i = 0; i < numoutputfields; i++) { 835 CeedVector vec; 836 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 837 if (vec == CEED_VECTOR_ACTIVE) { 838 CeedElemRestriction rstr; 839 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisout); CeedChkBackend(ierr); 840 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 841 CeedChkBackend(ierr); 842 if (rstrout && rstrout != rstr) 843 // LCOV_EXCL_START 844 return CeedError(ceed, CEED_ERROR_BACKEND, 845 "Multi-field non-composite operator diagonal assembly not supported"); 846 // LCOV_EXCL_STOP 847 rstrout = rstr; 848 CeedEvalMode emode; 849 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 850 switch (emode) { 851 case CEED_EVAL_NONE: 852 case CEED_EVAL_INTERP: 853 ierr = CeedRealloc(numemodeout + 1, &emodeout); CeedChkBackend(ierr); 854 emodeout[numemodeout] = emode; 855 numemodeout += 1; 856 break; 857 case CEED_EVAL_GRAD: 858 ierr = CeedRealloc(numemodeout + dim, &emodeout); CeedChkBackend(ierr); 859 for (CeedInt d = 0; d < dim; d++) 860 emodeout[numemodeout+d] = emode; 861 numemodeout += dim; 862 break; 863 case CEED_EVAL_WEIGHT: 864 case CEED_EVAL_DIV: 865 case CEED_EVAL_CURL: 866 break; // Caught by QF Assembly 867 } 868 } 869 } 870 871 // Operator data struct 872 CeedOperator_Hip *impl; 873 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 874 ierr = CeedCalloc(1, &impl->diag); CeedChkBackend(ierr); 875 CeedOperatorDiag_Hip *diag = impl->diag; 876 diag->basisin = basisin; 877 diag->basisout = basisout; 878 diag->h_emodein = emodein; 879 diag->h_emodeout = emodeout; 880 diag->numemodein = numemodein; 881 diag->numemodeout = numemodeout; 882 883 // Assemble kernel 884 885 char *diagonal_kernel_path, *diagonal_kernel_source; 886 ierr = CeedGetJitAbsolutePath(ceed, 887 "ceed/jit-source/hip/hip-ref-operator-assemble-diagonal.h", 888 &diagonal_kernel_path); CeedChkBackend(ierr); 889 CeedDebug256(ceed, 2, "----- Loading Diagonal Assembly Kernel Source -----\n"); 890 ierr = CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, 891 &diagonal_kernel_source); 892 CeedChkBackend(ierr); 893 CeedDebug256(ceed, 2, 894 "----- Loading Diagonal Assembly Source Complete! -----\n"); 895 CeedInt nnodes, nqpts; 896 ierr = CeedBasisGetNumNodes(basisin, &nnodes); CeedChkBackend(ierr); 897 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChkBackend(ierr); 898 diag->nnodes = nnodes; 899 ierr = CeedCompileHip(ceed, diagonal_kernel_source, &diag->module, 5, 900 "NUMEMODEIN", numemodein, 901 "NUMEMODEOUT", numemodeout, 902 "NNODES", nnodes, 903 "NQPTS", nqpts, 904 "NCOMP", ncomp 905 ); CeedChk_Hip(ceed, ierr); 906 ierr = CeedGetKernelHip(ceed, diag->module, "linearDiagonal", 907 &diag->linearDiagonal); CeedChk_Hip(ceed, ierr); 908 ierr = CeedGetKernelHip(ceed, diag->module, "linearPointBlockDiagonal", 909 &diag->linearPointBlock); 910 CeedChk_Hip(ceed, ierr); 911 ierr = CeedFree(&diagonal_kernel_path); CeedChkBackend(ierr); 912 ierr = CeedFree(&diagonal_kernel_source); CeedChkBackend(ierr); 913 914 // Basis matrices 915 const CeedInt qBytes = nqpts * sizeof(CeedScalar); 916 const CeedInt iBytes = qBytes * nnodes; 917 const CeedInt gBytes = qBytes * nnodes * dim; 918 const CeedInt eBytes = sizeof(CeedEvalMode); 919 const CeedScalar *interpin, *interpout, *gradin, *gradout; 920 921 // CEED_EVAL_NONE 922 CeedScalar *identity = NULL; 923 bool evalNone = false; 924 for (CeedInt i=0; i<numemodein; i++) 925 evalNone = evalNone || (emodein[i] == CEED_EVAL_NONE); 926 for (CeedInt i=0; i<numemodeout; i++) 927 evalNone = evalNone || (emodeout[i] == CEED_EVAL_NONE); 928 if (evalNone) { 929 ierr = CeedCalloc(nqpts*nnodes, &identity); CeedChkBackend(ierr); 930 for (CeedInt i=0; i<(nnodes<nqpts?nnodes:nqpts); i++) 931 identity[i*nnodes+i] = 1.0; 932 ierr = hipMalloc((void **)&diag->d_identity, iBytes); CeedChk_Hip(ceed, ierr); 933 ierr = hipMemcpy(diag->d_identity, identity, iBytes, 934 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 935 } 936 937 // CEED_EVAL_INTERP 938 ierr = CeedBasisGetInterp(basisin, &interpin); CeedChkBackend(ierr); 939 ierr = hipMalloc((void **)&diag->d_interpin, iBytes); CeedChk_Hip(ceed, ierr); 940 ierr = hipMemcpy(diag->d_interpin, interpin, iBytes, 941 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 942 ierr = CeedBasisGetInterp(basisout, &interpout); CeedChkBackend(ierr); 943 ierr = hipMalloc((void **)&diag->d_interpout, iBytes); CeedChk_Hip(ceed, ierr); 944 ierr = hipMemcpy(diag->d_interpout, interpout, iBytes, 945 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 946 947 // CEED_EVAL_GRAD 948 ierr = CeedBasisGetGrad(basisin, &gradin); CeedChkBackend(ierr); 949 ierr = hipMalloc((void **)&diag->d_gradin, gBytes); CeedChk_Hip(ceed, ierr); 950 ierr = hipMemcpy(diag->d_gradin, gradin, gBytes, 951 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 952 ierr = CeedBasisGetGrad(basisout, &gradout); CeedChkBackend(ierr); 953 ierr = hipMalloc((void **)&diag->d_gradout, gBytes); CeedChk_Hip(ceed, ierr); 954 ierr = hipMemcpy(diag->d_gradout, gradout, gBytes, 955 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 956 957 // Arrays of emodes 958 ierr = hipMalloc((void **)&diag->d_emodein, numemodein * eBytes); 959 CeedChk_Hip(ceed, ierr); 960 ierr = hipMemcpy(diag->d_emodein, emodein, numemodein * eBytes, 961 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 962 ierr = hipMalloc((void **)&diag->d_emodeout, numemodeout * eBytes); 963 CeedChk_Hip(ceed, ierr); 964 ierr = hipMemcpy(diag->d_emodeout, emodeout, numemodeout * eBytes, 965 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 966 967 // Restriction 968 diag->diagrstr = rstrout; 969 970 return CEED_ERROR_SUCCESS; 971 } 972 973 //------------------------------------------------------------------------------ 974 // Assemble diagonal common code 975 //------------------------------------------------------------------------------ 976 static inline int CeedOperatorAssembleDiagonalCore_Hip(CeedOperator op, 977 CeedVector assembled, CeedRequest *request, const bool pointBlock) { 978 int ierr; 979 Ceed ceed; 980 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 981 CeedOperator_Hip *impl; 982 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 983 984 // Assemble QFunction 985 CeedVector assembledqf; 986 CeedElemRestriction rstr; 987 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembledqf, 988 &rstr, request); CeedChkBackend(ierr); 989 ierr = CeedElemRestrictionDestroy(&rstr); CeedChkBackend(ierr); 990 CeedScalar maxnorm = 0; 991 ierr = CeedVectorNorm(assembledqf, CEED_NORM_MAX, &maxnorm); 992 CeedChkBackend(ierr); 993 994 // Setup 995 if (!impl->diag) { 996 ierr = CeedOperatorAssembleDiagonalSetup_Hip(op, pointBlock); 997 CeedChkBackend(ierr); 998 } 999 CeedOperatorDiag_Hip *diag = impl->diag; 1000 assert(diag != NULL); 1001 1002 // Restriction 1003 if (pointBlock && !diag->pbdiagrstr) { 1004 CeedElemRestriction pbdiagrstr; 1005 ierr = CreatePBRestriction(diag->diagrstr, &pbdiagrstr); CeedChkBackend(ierr); 1006 diag->pbdiagrstr = pbdiagrstr; 1007 } 1008 CeedElemRestriction diagrstr = pointBlock ? diag->pbdiagrstr : diag->diagrstr; 1009 1010 // Create diagonal vector 1011 CeedVector elemdiag = pointBlock ? diag->pbelemdiag : diag->elemdiag; 1012 if (!elemdiag) { 1013 // Element diagonal vector 1014 ierr = CeedElemRestrictionCreateVector(diagrstr, NULL, &elemdiag); 1015 CeedChkBackend(ierr); 1016 if (pointBlock) 1017 diag->pbelemdiag = elemdiag; 1018 else 1019 diag->elemdiag = elemdiag; 1020 } 1021 ierr = CeedVectorSetValue(elemdiag, 0.0); CeedChkBackend(ierr); 1022 1023 // Assemble element operator diagonals 1024 CeedScalar *elemdiagarray; 1025 const CeedScalar *assembledqfarray; 1026 ierr = CeedVectorGetArray(elemdiag, CEED_MEM_DEVICE, &elemdiagarray); 1027 CeedChkBackend(ierr); 1028 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_DEVICE, &assembledqfarray); 1029 CeedChkBackend(ierr); 1030 CeedInt nelem; 1031 ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); 1032 CeedChkBackend(ierr); 1033 1034 // Compute the diagonal of B^T D B 1035 int elemsPerBlock = 1; 1036 int grid = nelem/elemsPerBlock+((nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1037 void *args[] = {(void *) &nelem, (void *) &maxnorm, &diag->d_identity, 1038 &diag->d_interpin, &diag->d_gradin, &diag->d_interpout, 1039 &diag->d_gradout, &diag->d_emodein, &diag->d_emodeout, 1040 &assembledqfarray, &elemdiagarray 1041 }; 1042 if (pointBlock) { 1043 ierr = CeedRunKernelDimHip(ceed, diag->linearPointBlock, grid, 1044 diag->nnodes, 1, elemsPerBlock, args); 1045 CeedChkBackend(ierr); 1046 } else { 1047 ierr = CeedRunKernelDimHip(ceed, diag->linearDiagonal, grid, 1048 diag->nnodes, 1, elemsPerBlock, args); 1049 CeedChkBackend(ierr); 1050 } 1051 1052 // Restore arrays 1053 ierr = CeedVectorRestoreArray(elemdiag, &elemdiagarray); CeedChkBackend(ierr); 1054 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1055 CeedChkBackend(ierr); 1056 1057 // Assemble local operator diagonal 1058 ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag, 1059 assembled, request); CeedChkBackend(ierr); 1060 1061 // Cleanup 1062 ierr = CeedVectorDestroy(&assembledqf); CeedChkBackend(ierr); 1063 1064 return CEED_ERROR_SUCCESS; 1065 } 1066 1067 //------------------------------------------------------------------------------ 1068 // Assemble composite diagonal common code 1069 //------------------------------------------------------------------------------ 1070 static inline int CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip( 1071 CeedOperator op, CeedVector assembled, CeedRequest *request, 1072 const bool pointBlock) { 1073 int ierr; 1074 CeedInt numSub; 1075 CeedOperator *subOperators; 1076 ierr = CeedOperatorGetNumSub(op, &numSub); CeedChkBackend(ierr); 1077 ierr = CeedOperatorGetSubList(op, &subOperators); CeedChkBackend(ierr); 1078 for (CeedInt i = 0; i < numSub; i++) { 1079 ierr = CeedOperatorAssembleDiagonalCore_Hip(subOperators[i], assembled, 1080 request, pointBlock); CeedChkBackend(ierr); 1081 } 1082 return CEED_ERROR_SUCCESS; 1083 } 1084 1085 //------------------------------------------------------------------------------ 1086 // Assemble Linear Diagonal 1087 //------------------------------------------------------------------------------ 1088 static int CeedOperatorLinearAssembleAddDiagonal_Hip(CeedOperator op, 1089 CeedVector assembled, CeedRequest *request) { 1090 int ierr; 1091 bool isComposite; 1092 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1093 if (isComposite) { 1094 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip(op, assembled, 1095 request, false); 1096 } else { 1097 return CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, false); 1098 } 1099 } 1100 1101 //------------------------------------------------------------------------------ 1102 // Assemble Linear Point Block Diagonal 1103 //------------------------------------------------------------------------------ 1104 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip(CeedOperator op, 1105 CeedVector assembled, CeedRequest *request) { 1106 int ierr; 1107 bool isComposite; 1108 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1109 if (isComposite) { 1110 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip(op, assembled, 1111 request, true); 1112 } else { 1113 return CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, true); 1114 } 1115 } 1116 1117 //------------------------------------------------------------------------------ 1118 // Single operator assembly setup 1119 //------------------------------------------------------------------------------ 1120 static int CeedSingleOperatorAssembleSetup_Hip(CeedOperator op) { 1121 int ierr; 1122 Ceed ceed; 1123 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1124 CeedOperator_Hip *impl; 1125 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1126 1127 // Get intput and output fields 1128 CeedInt num_input_fields, num_output_fields; 1129 CeedOperatorField *input_fields; 1130 CeedOperatorField *output_fields; 1131 ierr = CeedOperatorGetFields(op, &num_input_fields, &input_fields, 1132 &num_output_fields, &output_fields); CeedChkBackend(ierr); 1133 1134 // Determine active input basis eval mode 1135 CeedQFunction qf; 1136 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 1137 CeedQFunctionField *qf_fields; 1138 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 1139 CeedChkBackend(ierr); 1140 // Note that the kernel will treat each dimension of a gradient action separately; 1141 // i.e., when an active input has a CEED_EVAL_GRAD mode, num_emode_in will increment 1142 // by dim. However, for the purposes of loading the B matrices, it will be treated 1143 // as one mode, and we will load/copy the entire gradient matrix at once, so 1144 // num_B_in_mats_to_load will be incremented by 1. 1145 CeedInt num_emode_in = 0, dim = 1, num_B_in_mats_to_load = 0, size_B_in = 0; 1146 CeedEvalMode *eval_mode_in = NULL; //will be of size num_B_in_mats_load 1147 CeedBasis basis_in = NULL; 1148 CeedInt nqpts = 0, esize = 0; 1149 CeedElemRestriction rstr_in = NULL; 1150 for (CeedInt i=0; i<num_input_fields; i++) { 1151 CeedVector vec; 1152 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChkBackend(ierr); 1153 if (vec == CEED_VECTOR_ACTIVE) { 1154 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1155 CeedChkBackend(ierr); 1156 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr); 1157 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &nqpts); CeedChkBackend(ierr); 1158 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1159 CeedChkBackend(ierr); 1160 ierr = CeedElemRestrictionGetElementSize(rstr_in, &esize); CeedChkBackend(ierr); 1161 CeedEvalMode eval_mode; 1162 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1163 CeedChkBackend(ierr); 1164 if (eval_mode != CEED_EVAL_NONE) { 1165 ierr = CeedRealloc(num_B_in_mats_to_load + 1, &eval_mode_in); 1166 CeedChkBackend(ierr); 1167 eval_mode_in[num_B_in_mats_to_load] = eval_mode; 1168 num_B_in_mats_to_load += 1; 1169 if (eval_mode == CEED_EVAL_GRAD) { 1170 num_emode_in += dim; 1171 size_B_in += dim * esize * nqpts; 1172 } else { 1173 num_emode_in +=1; 1174 size_B_in += esize * nqpts; 1175 } 1176 } 1177 } 1178 } 1179 1180 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1181 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 1182 CeedChkBackend(ierr); 1183 CeedInt num_emode_out = 0, num_B_out_mats_to_load = 0, size_B_out = 0; 1184 CeedEvalMode *eval_mode_out = NULL; 1185 CeedBasis basis_out = NULL; 1186 CeedElemRestriction rstr_out = NULL; 1187 for (CeedInt i=0; i<num_output_fields; i++) { 1188 CeedVector vec; 1189 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChkBackend(ierr); 1190 if (vec == CEED_VECTOR_ACTIVE) { 1191 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1192 CeedChkBackend(ierr); 1193 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1194 CeedChkBackend(ierr); 1195 if (rstr_out && rstr_out != rstr_in) 1196 // LCOV_EXCL_START 1197 return CeedError(ceed, CEED_ERROR_BACKEND, 1198 "Multi-field non-composite operator assembly not supported"); 1199 // LCOV_EXCL_STOP 1200 CeedEvalMode eval_mode; 1201 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1202 CeedChkBackend(ierr); 1203 if (eval_mode != CEED_EVAL_NONE) { 1204 ierr = CeedRealloc(num_B_out_mats_to_load + 1, &eval_mode_out); 1205 CeedChkBackend(ierr); 1206 eval_mode_out[num_B_out_mats_to_load] = eval_mode; 1207 num_B_out_mats_to_load += 1; 1208 if (eval_mode == CEED_EVAL_GRAD) { 1209 num_emode_out += dim; 1210 size_B_out += dim * esize * nqpts; 1211 } else { 1212 num_emode_out +=1; 1213 size_B_out += esize * nqpts; 1214 } 1215 } 1216 } 1217 } 1218 1219 if (num_emode_in == 0 || num_emode_out == 0) 1220 // LCOV_EXCL_START 1221 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1222 "Cannot assemble operator without inputs/outputs"); 1223 // LCOV_EXCL_STOP 1224 1225 CeedInt nelem, ncomp; 1226 ierr = CeedElemRestrictionGetNumElements(rstr_in, &nelem); CeedChkBackend(ierr); 1227 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &ncomp); 1228 CeedChkBackend(ierr); 1229 1230 ierr = CeedCalloc(1, &impl->asmb); CeedChkBackend(ierr); 1231 CeedOperatorAssemble_Hip *asmb = impl->asmb; 1232 asmb->nelem = nelem; 1233 1234 // Compile kernels 1235 int elemsPerBlock = 1; 1236 asmb->elemsPerBlock = elemsPerBlock; 1237 CeedInt block_size = esize * esize * elemsPerBlock; 1238 char *assembly_kernel_path, *assembly_kernel_source; 1239 ierr = CeedGetJitAbsolutePath(ceed, 1240 "ceed/jit-source/hip/hip-ref-operator-assemble.h", 1241 &assembly_kernel_path); CeedChkBackend(ierr); 1242 CeedDebug256(ceed, 2, "----- Loading Assembly Kernel Source -----\n"); 1243 ierr = CeedLoadSourceToBuffer(ceed, assembly_kernel_path, 1244 &assembly_kernel_source); 1245 CeedChkBackend(ierr); 1246 CeedDebug256(ceed, 2, "----- Loading Assembly Source Complete! -----\n"); 1247 bool fallback = block_size > 1024; 1248 if (fallback) { // Use fallback kernel with 1D threadblock 1249 block_size = esize * elemsPerBlock; 1250 asmb->block_size_x = esize; 1251 asmb->block_size_y = 1; 1252 } else { // Use kernel with 2D threadblock 1253 asmb->block_size_x = esize; 1254 asmb->block_size_y = esize; 1255 } 1256 ierr = CeedCompileHip(ceed, assembly_kernel_source, &asmb->module, 7, 1257 "NELEM", nelem, 1258 "NUMEMODEIN", num_emode_in, 1259 "NUMEMODEOUT", num_emode_out, 1260 "NQPTS", nqpts, 1261 "NNODES", esize, 1262 "BLOCK_SIZE", block_size, 1263 "NCOMP", ncomp 1264 ); CeedChk_Hip(ceed, ierr); 1265 ierr = CeedGetKernelHip(ceed, asmb->module, 1266 fallback ? "linearAssembleFallback" : "linearAssemble", 1267 &asmb->linearAssemble); CeedChk_Hip(ceed, ierr); 1268 ierr = CeedFree(&assembly_kernel_path); CeedChkBackend(ierr); 1269 ierr = CeedFree(&assembly_kernel_source); CeedChkBackend(ierr); 1270 1271 // Build 'full' B matrices (not 1D arrays used for tensor-product matrices) 1272 const CeedScalar *interp_in, *grad_in; 1273 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr); 1274 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr); 1275 1276 // Load into B_in, in order that they will be used in eval_mode 1277 const CeedInt inBytes = size_B_in * sizeof(CeedScalar); 1278 CeedInt mat_start = 0; 1279 ierr = hipMalloc((void **) &asmb->d_B_in, inBytes); CeedChk_Hip(ceed, ierr); 1280 for (int i = 0; i < num_B_in_mats_to_load; i++) { 1281 CeedEvalMode eval_mode = eval_mode_in[i]; 1282 if (eval_mode == CEED_EVAL_INTERP) { 1283 ierr = hipMemcpy(&asmb->d_B_in[mat_start], interp_in, 1284 esize * nqpts * sizeof(CeedScalar), 1285 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1286 mat_start += esize * nqpts; 1287 } else if (eval_mode == CEED_EVAL_GRAD) { 1288 ierr = hipMemcpy(asmb->d_B_in, grad_in, 1289 dim * esize * nqpts * sizeof(CeedScalar), 1290 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1291 mat_start += dim * esize * nqpts; 1292 } 1293 } 1294 1295 const CeedScalar *interp_out, *grad_out; 1296 // Note that this function currently assumes 1 basis, so this should always be true 1297 // for now 1298 if (basis_out == basis_in) { 1299 interp_out = interp_in; 1300 grad_out = grad_in; 1301 } else { 1302 ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr); 1303 ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr); 1304 } 1305 1306 // Load into B_out, in order that they will be used in eval_mode 1307 const CeedInt outBytes = size_B_out * sizeof(CeedScalar); 1308 mat_start = 0; 1309 ierr = hipMalloc((void **) &asmb->d_B_out, outBytes); CeedChk_Hip(ceed, ierr); 1310 for (int i = 0; i < num_B_out_mats_to_load; i++) { 1311 CeedEvalMode eval_mode = eval_mode_out[i]; 1312 if (eval_mode == CEED_EVAL_INTERP) { 1313 ierr = hipMemcpy(&asmb->d_B_out[mat_start], interp_out, 1314 esize * nqpts * sizeof(CeedScalar), 1315 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1316 mat_start += esize * nqpts; 1317 } else if (eval_mode == CEED_EVAL_GRAD) { 1318 ierr = hipMemcpy(&asmb->d_B_out[mat_start], grad_out, 1319 dim * esize * nqpts * sizeof(CeedScalar), 1320 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1321 mat_start += dim * esize * nqpts; 1322 } 1323 } 1324 return CEED_ERROR_SUCCESS; 1325 } 1326 1327 //------------------------------------------------------------------------------ 1328 // Single operator assembly 1329 //------------------------------------------------------------------------------ 1330 static int CeedSingleOperatorAssemble_Hip(CeedOperator op, CeedInt offset, 1331 CeedVector values) { 1332 1333 int ierr; 1334 Ceed ceed; 1335 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1336 CeedOperator_Hip *impl; 1337 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1338 1339 // Setup 1340 if (!impl->asmb) { 1341 ierr = CeedSingleOperatorAssembleSetup_Hip(op); 1342 CeedChkBackend(ierr); 1343 assert(impl->asmb != NULL); 1344 } 1345 1346 // Assemble QFunction 1347 CeedVector assembled_qf; 1348 CeedElemRestriction rstr_q; 1349 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate( 1350 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChkBackend(ierr); 1351 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChkBackend(ierr); 1352 CeedScalar *values_array; 1353 ierr = CeedVectorGetArrayWrite(values, CEED_MEM_DEVICE, &values_array); 1354 CeedChkBackend(ierr); 1355 values_array += offset; 1356 const CeedScalar *qf_array; 1357 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &qf_array); 1358 CeedChkBackend(ierr); 1359 1360 // Compute B^T D B 1361 const CeedInt nelem = impl->asmb->nelem; // to satisfy clang-tidy 1362 const CeedInt elemsPerBlock = impl->asmb->elemsPerBlock; 1363 const CeedInt grid = nelem/elemsPerBlock+(( 1364 nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1365 void *args[] = {&impl->asmb->d_B_in, &impl->asmb->d_B_out, 1366 &qf_array, &values_array 1367 }; 1368 ierr = CeedRunKernelDimHip(ceed, impl->asmb->linearAssemble, grid, 1369 impl->asmb->block_size_x, impl->asmb->block_size_y, 1370 elemsPerBlock, args); 1371 CeedChkBackend(ierr); 1372 1373 1374 // Restore arrays 1375 ierr = CeedVectorRestoreArray(values, &values_array); CeedChkBackend(ierr); 1376 ierr = CeedVectorRestoreArrayRead(assembled_qf, &qf_array); 1377 CeedChkBackend(ierr); 1378 1379 // Cleanup 1380 ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr); 1381 1382 return CEED_ERROR_SUCCESS; 1383 } 1384 1385 //------------------------------------------------------------------------------ 1386 // Assemble matrix data for COO matrix of assembled operator. 1387 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1388 // 1389 // Note that this (and other assembly routines) currently assume only one 1390 // active input restriction/basis per operator (could have multiple basis eval 1391 // modes). 1392 // TODO: allow multiple active input restrictions/basis objects 1393 //------------------------------------------------------------------------------ 1394 int CeedOperatorLinearAssemble_Hip(CeedOperator op, CeedVector values) { 1395 1396 // As done in the default implementation, loop through suboperators 1397 // for composite operators, or call single operator assembly otherwise 1398 bool is_composite; 1399 CeedInt ierr; 1400 ierr = CeedOperatorIsComposite(op, &is_composite); CeedChkBackend(ierr); 1401 1402 CeedElemRestriction rstr; 1403 CeedInt num_elem, elem_size, num_comp; 1404 1405 CeedInt offset = 0; 1406 if (is_composite) { 1407 CeedInt num_suboperators; 1408 ierr = CeedOperatorGetNumSub(op, &num_suboperators); CeedChkBackend(ierr); 1409 CeedOperator *sub_operators; 1410 ierr = CeedOperatorGetSubList(op, &sub_operators); CeedChkBackend(ierr); 1411 for (int k = 0; k < num_suboperators; ++k) { 1412 ierr = CeedSingleOperatorAssemble_Hip(sub_operators[k], offset, values); 1413 CeedChkBackend(ierr); 1414 ierr = CeedOperatorGetActiveElemRestriction(sub_operators[k], &rstr); 1415 CeedChkBackend(ierr); 1416 ierr = CeedElemRestrictionGetNumElements(rstr, &num_elem); CeedChkBackend(ierr); 1417 ierr = CeedElemRestrictionGetElementSize(rstr, &elem_size); 1418 CeedChkBackend(ierr); 1419 ierr = CeedElemRestrictionGetNumComponents(rstr, &num_comp); 1420 CeedChkBackend(ierr); 1421 offset += elem_size*num_comp * elem_size*num_comp * num_elem; 1422 } 1423 } else { 1424 ierr = CeedSingleOperatorAssemble_Hip(op, offset, values); CeedChkBackend(ierr); 1425 } 1426 1427 return CEED_ERROR_SUCCESS; 1428 } 1429 1430 //------------------------------------------------------------------------------ 1431 // Create operator 1432 //------------------------------------------------------------------------------ 1433 int CeedOperatorCreate_Hip(CeedOperator op) { 1434 int ierr; 1435 Ceed ceed; 1436 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1437 CeedOperator_Hip *impl; 1438 1439 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1440 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1441 1442 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1443 CeedOperatorLinearAssembleQFunction_Hip); 1444 CeedChkBackend(ierr); 1445 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1446 "LinearAssembleQFunctionUpdate", 1447 CeedOperatorLinearAssembleQFunctionUpdate_Hip); 1448 CeedChkBackend(ierr); 1449 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1450 CeedOperatorLinearAssembleAddDiagonal_Hip); 1451 CeedChkBackend(ierr); 1452 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1453 "LinearAssembleAddPointBlockDiagonal", 1454 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1455 CeedChkBackend(ierr); 1456 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1457 "LinearAssemble", CeedOperatorLinearAssemble_Hip); 1458 CeedChkBackend(ierr); 1459 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1460 CeedOperatorApplyAdd_Hip); CeedChkBackend(ierr); 1461 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1462 CeedOperatorDestroy_Hip); CeedChkBackend(ierr); 1463 return CEED_ERROR_SUCCESS; 1464 } 1465 1466 //------------------------------------------------------------------------------ 1467 // Composite Operator Create 1468 //------------------------------------------------------------------------------ 1469 int CeedCompositeOperatorCreate_Hip(CeedOperator op) { 1470 int ierr; 1471 Ceed ceed; 1472 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1473 1474 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1475 CeedOperatorLinearAssembleAddDiagonal_Hip); 1476 CeedChkBackend(ierr); 1477 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1478 "LinearAssembleAddPointBlockDiagonal", 1479 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1480 CeedChkBackend(ierr); 1481 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1482 "LinearAssemble", CeedOperatorLinearAssemble_Hip); 1483 CeedChkBackend(ierr); 1484 return CEED_ERROR_SUCCESS; 1485 } 1486 //------------------------------------------------------------------------------ 1487