1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <hip/hip_runtime.h> 20 #include <assert.h> 21 #include <stdbool.h> 22 #include <string.h> 23 #include "ceed-hip-ref.h" 24 #include "../hip/ceed-hip-compile.h" 25 26 //------------------------------------------------------------------------------ 27 // Destroy operator 28 //------------------------------------------------------------------------------ 29 static int CeedOperatorDestroy_Hip(CeedOperator op) { 30 int ierr; 31 CeedOperator_Hip *impl; 32 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 33 34 // Apply data 35 for (CeedInt i = 0; i < impl->numein + impl->numeout; i++) { 36 ierr = CeedVectorDestroy(&impl->evecs[i]); CeedChkBackend(ierr); 37 } 38 ierr = CeedFree(&impl->evecs); CeedChkBackend(ierr); 39 40 for (CeedInt i = 0; i < impl->numein; i++) { 41 ierr = CeedVectorDestroy(&impl->qvecsin[i]); CeedChkBackend(ierr); 42 } 43 ierr = CeedFree(&impl->qvecsin); CeedChkBackend(ierr); 44 45 for (CeedInt i = 0; i < impl->numeout; i++) { 46 ierr = CeedVectorDestroy(&impl->qvecsout[i]); CeedChkBackend(ierr); 47 } 48 ierr = CeedFree(&impl->qvecsout); CeedChkBackend(ierr); 49 50 // QFunction diagonal assembly data 51 for (CeedInt i=0; i<impl->qfnumactivein; i++) { 52 ierr = CeedVectorDestroy(&impl->qfactivein[i]); CeedChkBackend(ierr); 53 } 54 ierr = CeedFree(&impl->qfactivein); CeedChkBackend(ierr); 55 56 // Diag data 57 if (impl->diag) { 58 Ceed ceed; 59 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 60 CeedChk_Hip(ceed, hipModuleUnload(impl->diag->module)); 61 ierr = CeedFree(&impl->diag->h_emodein); CeedChkBackend(ierr); 62 ierr = CeedFree(&impl->diag->h_emodeout); CeedChkBackend(ierr); 63 ierr = hipFree(impl->diag->d_emodein); CeedChk_Hip(ceed, ierr); 64 ierr = hipFree(impl->diag->d_emodeout); CeedChk_Hip(ceed, ierr); 65 ierr = hipFree(impl->diag->d_identity); CeedChk_Hip(ceed, ierr); 66 ierr = hipFree(impl->diag->d_interpin); CeedChk_Hip(ceed, ierr); 67 ierr = hipFree(impl->diag->d_interpout); CeedChk_Hip(ceed, ierr); 68 ierr = hipFree(impl->diag->d_gradin); CeedChk_Hip(ceed, ierr); 69 ierr = hipFree(impl->diag->d_gradout); CeedChk_Hip(ceed, ierr); 70 ierr = CeedElemRestrictionDestroy(&impl->diag->pbdiagrstr); 71 CeedChkBackend(ierr); 72 ierr = CeedVectorDestroy(&impl->diag->elemdiag); CeedChkBackend(ierr); 73 ierr = CeedVectorDestroy(&impl->diag->pbelemdiag); CeedChkBackend(ierr); 74 } 75 ierr = CeedFree(&impl->diag); CeedChkBackend(ierr); 76 77 ierr = CeedFree(&impl); CeedChkBackend(ierr); 78 return CEED_ERROR_SUCCESS; 79 } 80 81 //------------------------------------------------------------------------------ 82 // Setup infields or outfields 83 //------------------------------------------------------------------------------ 84 static int CeedOperatorSetupFields_Hip(CeedQFunction qf, CeedOperator op, 85 bool isinput, CeedVector *evecs, 86 CeedVector *qvecs, CeedInt starte, 87 CeedInt numfields, CeedInt Q, 88 CeedInt numelements) { 89 CeedInt dim, ierr, size; 90 CeedSize q_size; 91 Ceed ceed; 92 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 93 CeedBasis basis; 94 CeedElemRestriction Erestrict; 95 CeedOperatorField *opfields; 96 CeedQFunctionField *qffields; 97 CeedVector fieldvec; 98 bool strided; 99 bool skiprestrict; 100 101 if (isinput) { 102 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 103 CeedChkBackend(ierr); 104 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 105 CeedChkBackend(ierr); 106 } else { 107 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 108 CeedChkBackend(ierr); 109 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 110 CeedChkBackend(ierr); 111 } 112 113 // Loop over fields 114 for (CeedInt i = 0; i < numfields; i++) { 115 CeedEvalMode emode; 116 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 117 118 strided = false; 119 skiprestrict = false; 120 if (emode != CEED_EVAL_WEIGHT) { 121 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &Erestrict); 122 CeedChkBackend(ierr); 123 124 // Check whether this field can skip the element restriction: 125 // must be passive input, with emode NONE, and have a strided restriction with 126 // CEED_STRIDES_BACKEND. 127 128 // First, check whether the field is input or output: 129 if (isinput) { 130 // Check for passive input: 131 ierr = CeedOperatorFieldGetVector(opfields[i], &fieldvec); CeedChkBackend(ierr); 132 if (fieldvec != CEED_VECTOR_ACTIVE) { 133 // Check emode 134 if (emode == CEED_EVAL_NONE) { 135 // Check for strided restriction 136 ierr = CeedElemRestrictionIsStrided(Erestrict, &strided); 137 CeedChkBackend(ierr); 138 if (strided) { 139 // Check if vector is already in preferred backend ordering 140 ierr = CeedElemRestrictionHasBackendStrides(Erestrict, 141 &skiprestrict); CeedChkBackend(ierr); 142 } 143 } 144 } 145 } 146 if (skiprestrict) { 147 // We do not need an E-Vector, but will use the input field vector's data 148 // directly in the operator application. 149 evecs[i + starte] = NULL; 150 } else { 151 ierr = CeedElemRestrictionCreateVector(Erestrict, NULL, 152 &evecs[i + starte]); 153 CeedChkBackend(ierr); 154 } 155 } 156 157 switch (emode) { 158 case CEED_EVAL_NONE: 159 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 160 q_size = (CeedSize)numelements * Q * size; 161 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 162 break; 163 case CEED_EVAL_INTERP: 164 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 165 q_size = (CeedSize)numelements * Q * size; 166 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 167 break; 168 case CEED_EVAL_GRAD: 169 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 170 ierr = CeedQFunctionFieldGetSize(qffields[i], &size); CeedChkBackend(ierr); 171 ierr = CeedBasisGetDimension(basis, &dim); CeedChkBackend(ierr); 172 q_size = (CeedSize)numelements * Q * size; 173 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 174 break; 175 case CEED_EVAL_WEIGHT: // Only on input fields 176 ierr = CeedOperatorFieldGetBasis(opfields[i], &basis); CeedChkBackend(ierr); 177 q_size = (CeedSize)numelements * Q; 178 ierr = CeedVectorCreate(ceed, q_size, &qvecs[i]); CeedChkBackend(ierr); 179 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 180 CEED_EVAL_WEIGHT, NULL, qvecs[i]); CeedChkBackend(ierr); 181 break; 182 case CEED_EVAL_DIV: 183 break; // TODO: Not implemented 184 case CEED_EVAL_CURL: 185 break; // TODO: Not implemented 186 } 187 } 188 return CEED_ERROR_SUCCESS; 189 } 190 191 //------------------------------------------------------------------------------ 192 // CeedOperator needs to connect all the named fields (be they active or passive) 193 // to the named inputs and outputs of its CeedQFunction. 194 //------------------------------------------------------------------------------ 195 static int CeedOperatorSetup_Hip(CeedOperator op) { 196 int ierr; 197 bool setupdone; 198 ierr = CeedOperatorIsSetupDone(op, &setupdone); CeedChkBackend(ierr); 199 if (setupdone) 200 return CEED_ERROR_SUCCESS; 201 Ceed ceed; 202 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 203 CeedOperator_Hip *impl; 204 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 205 CeedQFunction qf; 206 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 207 CeedInt Q, numelements, numinputfields, numoutputfields; 208 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 209 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 210 CeedOperatorField *opinputfields, *opoutputfields; 211 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 212 &numoutputfields, &opoutputfields); 213 CeedChkBackend(ierr); 214 CeedQFunctionField *qfinputfields, *qfoutputfields; 215 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 216 CeedChkBackend(ierr); 217 218 // Allocate 219 ierr = CeedCalloc(numinputfields + numoutputfields, &impl->evecs); 220 CeedChkBackend(ierr); 221 222 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsin); CeedChkBackend(ierr); 223 ierr = CeedCalloc(CEED_FIELD_MAX, &impl->qvecsout); CeedChkBackend(ierr); 224 225 impl->numein = numinputfields; impl->numeout = numoutputfields; 226 227 // Set up infield and outfield evecs and qvecs 228 // Infields 229 ierr = CeedOperatorSetupFields_Hip(qf, op, true, 230 impl->evecs, impl->qvecsin, 0, 231 numinputfields, Q, numelements); 232 CeedChkBackend(ierr); 233 234 // Outfields 235 ierr = CeedOperatorSetupFields_Hip(qf, op, false, 236 impl->evecs, impl->qvecsout, 237 numinputfields, numoutputfields, Q, 238 numelements); CeedChkBackend(ierr); 239 240 ierr = CeedOperatorSetSetupDone(op); CeedChkBackend(ierr); 241 return CEED_ERROR_SUCCESS; 242 } 243 244 //------------------------------------------------------------------------------ 245 // Setup Operator Inputs 246 //------------------------------------------------------------------------------ 247 static inline int CeedOperatorSetupInputs_Hip(CeedInt numinputfields, 248 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 249 CeedVector invec, const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 250 CeedOperator_Hip *impl, CeedRequest *request) { 251 CeedInt ierr; 252 CeedEvalMode emode; 253 CeedVector vec; 254 CeedElemRestriction Erestrict; 255 256 for (CeedInt i = 0; i < numinputfields; i++) { 257 // Get input vector 258 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 259 if (vec == CEED_VECTOR_ACTIVE) { 260 if (skipactive) 261 continue; 262 else 263 vec = invec; 264 } 265 266 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 267 CeedChkBackend(ierr); 268 if (emode == CEED_EVAL_WEIGHT) { // Skip 269 } else { 270 // Get input vector 271 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 272 // Get input element restriction 273 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 274 CeedChkBackend(ierr); 275 if (vec == CEED_VECTOR_ACTIVE) 276 vec = invec; 277 // Restrict, if necessary 278 if (!impl->evecs[i]) { 279 // No restriction for this field; read data directly from vec. 280 ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, 281 (const CeedScalar **) &edata[i]); 282 CeedChkBackend(ierr); 283 } else { 284 ierr = CeedElemRestrictionApply(Erestrict, CEED_NOTRANSPOSE, vec, 285 impl->evecs[i], request); CeedChkBackend(ierr); 286 // Get evec 287 ierr = CeedVectorGetArrayRead(impl->evecs[i], CEED_MEM_DEVICE, 288 (const CeedScalar **) &edata[i]); 289 CeedChkBackend(ierr); 290 } 291 } 292 } 293 return CEED_ERROR_SUCCESS; 294 } 295 296 //------------------------------------------------------------------------------ 297 // Input Basis Action 298 //------------------------------------------------------------------------------ 299 static inline int CeedOperatorInputBasis_Hip(CeedInt numelements, 300 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 301 CeedInt numinputfields, const bool skipactive, 302 CeedScalar *edata[2*CEED_FIELD_MAX], CeedOperator_Hip *impl) { 303 CeedInt ierr; 304 CeedInt elemsize, size; 305 CeedElemRestriction Erestrict; 306 CeedEvalMode emode; 307 CeedBasis basis; 308 309 for (CeedInt i=0; i<numinputfields; i++) { 310 // Skip active input 311 if (skipactive) { 312 CeedVector vec; 313 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 314 if (vec == CEED_VECTOR_ACTIVE) 315 continue; 316 } 317 // Get elemsize, emode, size 318 ierr = CeedOperatorFieldGetElemRestriction(opinputfields[i], &Erestrict); 319 CeedChkBackend(ierr); 320 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 321 CeedChkBackend(ierr); 322 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 323 CeedChkBackend(ierr); 324 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 325 // Basis action 326 switch (emode) { 327 case CEED_EVAL_NONE: 328 ierr = CeedVectorSetArray(impl->qvecsin[i], CEED_MEM_DEVICE, 329 CEED_USE_POINTER, edata[i]); CeedChkBackend(ierr); 330 break; 331 case CEED_EVAL_INTERP: 332 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 333 CeedChkBackend(ierr); 334 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 335 CEED_EVAL_INTERP, impl->evecs[i], 336 impl->qvecsin[i]); CeedChkBackend(ierr); 337 break; 338 case CEED_EVAL_GRAD: 339 ierr = CeedOperatorFieldGetBasis(opinputfields[i], &basis); 340 CeedChkBackend(ierr); 341 ierr = CeedBasisApply(basis, numelements, CEED_NOTRANSPOSE, 342 CEED_EVAL_GRAD, impl->evecs[i], 343 impl->qvecsin[i]); CeedChkBackend(ierr); 344 break; 345 case CEED_EVAL_WEIGHT: 346 break; // No action 347 case CEED_EVAL_DIV: 348 break; // TODO: Not implemented 349 case CEED_EVAL_CURL: 350 break; // TODO: Not implemented 351 } 352 } 353 return CEED_ERROR_SUCCESS; 354 } 355 356 //------------------------------------------------------------------------------ 357 // Restore Input Vectors 358 //------------------------------------------------------------------------------ 359 static inline int CeedOperatorRestoreInputs_Hip(CeedInt numinputfields, 360 CeedQFunctionField *qfinputfields, CeedOperatorField *opinputfields, 361 const bool skipactive, CeedScalar *edata[2*CEED_FIELD_MAX], 362 CeedOperator_Hip *impl) { 363 CeedInt ierr; 364 CeedEvalMode emode; 365 CeedVector vec; 366 367 for (CeedInt i = 0; i < numinputfields; i++) { 368 // Skip active input 369 if (skipactive) { 370 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 371 if (vec == CEED_VECTOR_ACTIVE) 372 continue; 373 } 374 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 375 CeedChkBackend(ierr); 376 if (emode == CEED_EVAL_WEIGHT) { // Skip 377 } else { 378 if (!impl->evecs[i]) { // This was a skiprestrict case 379 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 380 ierr = CeedVectorRestoreArrayRead(vec, 381 (const CeedScalar **)&edata[i]); 382 CeedChkBackend(ierr); 383 } else { 384 ierr = CeedVectorRestoreArrayRead(impl->evecs[i], 385 (const CeedScalar **) &edata[i]); 386 CeedChkBackend(ierr); 387 } 388 } 389 } 390 return CEED_ERROR_SUCCESS; 391 } 392 393 //------------------------------------------------------------------------------ 394 // Apply and add to output 395 //------------------------------------------------------------------------------ 396 static int CeedOperatorApplyAdd_Hip(CeedOperator op, CeedVector invec, 397 CeedVector outvec, CeedRequest *request) { 398 int ierr; 399 CeedOperator_Hip *impl; 400 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 401 CeedQFunction qf; 402 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 403 CeedInt Q, numelements, elemsize, numinputfields, numoutputfields, size; 404 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 405 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 406 CeedOperatorField *opinputfields, *opoutputfields; 407 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 408 &numoutputfields, &opoutputfields); 409 CeedChkBackend(ierr); 410 CeedQFunctionField *qfinputfields, *qfoutputfields; 411 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 412 CeedChkBackend(ierr); 413 CeedEvalMode emode; 414 CeedVector vec; 415 CeedBasis basis; 416 CeedElemRestriction Erestrict; 417 CeedScalar *edata[2*CEED_FIELD_MAX]; 418 419 // Setup 420 ierr = CeedOperatorSetup_Hip(op); CeedChkBackend(ierr); 421 422 // Input Evecs and Restriction 423 ierr = CeedOperatorSetupInputs_Hip(numinputfields, qfinputfields, 424 opinputfields, invec, false, edata, 425 impl, request); CeedChkBackend(ierr); 426 427 // Input basis apply if needed 428 ierr = CeedOperatorInputBasis_Hip(numelements, qfinputfields, opinputfields, 429 numinputfields, false, edata, impl); 430 CeedChkBackend(ierr); 431 432 // Output pointers, as necessary 433 for (CeedInt i = 0; i < numoutputfields; i++) { 434 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 435 CeedChkBackend(ierr); 436 if (emode == CEED_EVAL_NONE) { 437 // Set the output Q-Vector to use the E-Vector data directly. 438 ierr = CeedVectorGetArrayWrite(impl->evecs[i + impl->numein], CEED_MEM_DEVICE, 439 &edata[i + numinputfields]); CeedChkBackend(ierr); 440 ierr = CeedVectorSetArray(impl->qvecsout[i], CEED_MEM_DEVICE, 441 CEED_USE_POINTER, edata[i + numinputfields]); 442 CeedChkBackend(ierr); 443 } 444 } 445 446 // Q function 447 ierr = CeedQFunctionApply(qf, numelements * Q, impl->qvecsin, impl->qvecsout); 448 CeedChkBackend(ierr); 449 450 // Output basis apply if needed 451 for (CeedInt i = 0; i < numoutputfields; i++) { 452 // Get elemsize, emode, size 453 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 454 CeedChkBackend(ierr); 455 ierr = CeedElemRestrictionGetElementSize(Erestrict, &elemsize); 456 CeedChkBackend(ierr); 457 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 458 CeedChkBackend(ierr); 459 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 460 CeedChkBackend(ierr); 461 // Basis action 462 switch (emode) { 463 case CEED_EVAL_NONE: 464 break; 465 case CEED_EVAL_INTERP: 466 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 467 CeedChkBackend(ierr); 468 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 469 CEED_EVAL_INTERP, impl->qvecsout[i], 470 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 471 break; 472 case CEED_EVAL_GRAD: 473 ierr = CeedOperatorFieldGetBasis(opoutputfields[i], &basis); 474 CeedChkBackend(ierr); 475 ierr = CeedBasisApply(basis, numelements, CEED_TRANSPOSE, 476 CEED_EVAL_GRAD, impl->qvecsout[i], 477 impl->evecs[i + impl->numein]); CeedChkBackend(ierr); 478 break; 479 // LCOV_EXCL_START 480 case CEED_EVAL_WEIGHT: { 481 Ceed ceed; 482 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 483 return CeedError(ceed, CEED_ERROR_BACKEND, 484 "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 485 break; // Should not occur 486 } 487 case CEED_EVAL_DIV: 488 break; // TODO: Not implemented 489 case CEED_EVAL_CURL: 490 break; // TODO: Not implemented 491 // LCOV_EXCL_STOP 492 } 493 } 494 495 // Output restriction 496 for (CeedInt i = 0; i < numoutputfields; i++) { 497 // Restore evec 498 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 499 CeedChkBackend(ierr); 500 if (emode == CEED_EVAL_NONE) { 501 ierr = CeedVectorRestoreArray(impl->evecs[i+impl->numein], 502 &edata[i + numinputfields]); 503 CeedChkBackend(ierr); 504 } 505 // Get output vector 506 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 507 CeedChkBackend(ierr); 508 // Restrict 509 ierr = CeedOperatorFieldGetElemRestriction(opoutputfields[i], &Erestrict); 510 CeedChkBackend(ierr); 511 // Active 512 if (vec == CEED_VECTOR_ACTIVE) 513 vec = outvec; 514 515 ierr = CeedElemRestrictionApply(Erestrict, CEED_TRANSPOSE, 516 impl->evecs[i + impl->numein], vec, 517 request); CeedChkBackend(ierr); 518 } 519 520 // Restore input arrays 521 ierr = CeedOperatorRestoreInputs_Hip(numinputfields, qfinputfields, 522 opinputfields, false, edata, impl); 523 CeedChkBackend(ierr); 524 return CEED_ERROR_SUCCESS; 525 } 526 527 //------------------------------------------------------------------------------ 528 // Core code for assembling linear QFunction 529 //------------------------------------------------------------------------------ 530 static inline int CeedOperatorLinearAssembleQFunctionCore_Hip(CeedOperator op, 531 bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 532 CeedRequest *request) { 533 int ierr; 534 CeedOperator_Hip *impl; 535 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 536 CeedQFunction qf; 537 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 538 CeedInt Q, numelements, numinputfields, numoutputfields, size; 539 CeedSize q_size; 540 ierr = CeedOperatorGetNumQuadraturePoints(op, &Q); CeedChkBackend(ierr); 541 ierr = CeedOperatorGetNumElements(op, &numelements); CeedChkBackend(ierr); 542 CeedOperatorField *opinputfields, *opoutputfields; 543 ierr = CeedOperatorGetFields(op, &numinputfields, &opinputfields, 544 &numoutputfields, &opoutputfields); 545 CeedChkBackend(ierr); 546 CeedQFunctionField *qfinputfields, *qfoutputfields; 547 ierr = CeedQFunctionGetFields(qf, NULL, &qfinputfields, NULL, &qfoutputfields); 548 CeedChkBackend(ierr); 549 CeedVector vec; 550 CeedInt numactivein = impl->qfnumactivein, numactiveout = impl->qfnumactiveout; 551 CeedVector *activein = impl->qfactivein; 552 CeedScalar *a, *tmp; 553 Ceed ceed, ceedparent; 554 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 555 ierr = CeedGetOperatorFallbackParentCeed(ceed, &ceedparent); 556 CeedChkBackend(ierr); 557 ceedparent = ceedparent ? ceedparent : ceed; 558 CeedScalar *edata[2*CEED_FIELD_MAX]; 559 560 // Setup 561 ierr = CeedOperatorSetup_Hip(op); CeedChkBackend(ierr); 562 563 // Check for identity 564 bool identityqf; 565 ierr = CeedQFunctionIsIdentity(qf, &identityqf); CeedChkBackend(ierr); 566 if (identityqf) 567 // LCOV_EXCL_START 568 return CeedError(ceed, CEED_ERROR_BACKEND, 569 "Assembling identity QFunctions not supported"); 570 // LCOV_EXCL_STOP 571 572 // Input Evecs and Restriction 573 ierr = CeedOperatorSetupInputs_Hip(numinputfields, qfinputfields, 574 opinputfields, NULL, true, edata, 575 impl, request); CeedChkBackend(ierr); 576 577 // Count number of active input fields 578 if (!numactivein) { 579 for (CeedInt i=0; i<numinputfields; i++) { 580 // Get input vector 581 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 582 // Check if active input 583 if (vec == CEED_VECTOR_ACTIVE) { 584 ierr = CeedQFunctionFieldGetSize(qfinputfields[i], &size); CeedChkBackend(ierr); 585 ierr = CeedVectorSetValue(impl->qvecsin[i], 0.0); CeedChkBackend(ierr); 586 ierr = CeedVectorGetArray(impl->qvecsin[i], CEED_MEM_DEVICE, &tmp); 587 CeedChkBackend(ierr); 588 ierr = CeedRealloc(numactivein + size, &activein); CeedChkBackend(ierr); 589 for (CeedInt field = 0; field < size; field++) { 590 q_size = (CeedSize)Q*numelements; 591 ierr = CeedVectorCreate(ceed, q_size, &activein[numactivein+field]); 592 CeedChkBackend(ierr); 593 ierr = CeedVectorSetArray(activein[numactivein+field], CEED_MEM_DEVICE, 594 CEED_USE_POINTER, &tmp[field*Q*numelements]); 595 CeedChkBackend(ierr); 596 } 597 numactivein += size; 598 ierr = CeedVectorRestoreArray(impl->qvecsin[i], &tmp); CeedChkBackend(ierr); 599 } 600 } 601 impl->qfnumactivein = numactivein; 602 impl->qfactivein = activein; 603 } 604 605 // Count number of active output fields 606 if (!numactiveout) { 607 for (CeedInt i=0; i<numoutputfields; i++) { 608 // Get output vector 609 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 610 CeedChkBackend(ierr); 611 // Check if active output 612 if (vec == CEED_VECTOR_ACTIVE) { 613 ierr = CeedQFunctionFieldGetSize(qfoutputfields[i], &size); 614 CeedChkBackend(ierr); 615 numactiveout += size; 616 } 617 } 618 impl->qfnumactiveout = numactiveout; 619 } 620 621 // Check sizes 622 if (!numactivein || !numactiveout) 623 // LCOV_EXCL_START 624 return CeedError(ceed, CEED_ERROR_BACKEND, 625 "Cannot assemble QFunction without active inputs " 626 "and outputs"); 627 // LCOV_EXCL_STOP 628 629 // Build objects if needed 630 if (build_objects) { 631 // Create output restriction 632 CeedInt strides[3] = {1, numelements*Q, Q}; /* *NOPAD* */ 633 ierr = CeedElemRestrictionCreateStrided(ceedparent, numelements, Q, 634 numactivein*numactiveout, 635 numactivein*numactiveout*numelements*Q, 636 strides, rstr); CeedChkBackend(ierr); 637 // Create assembled vector 638 CeedSize l_size = (CeedSize)numelements*Q*numactivein*numactiveout; 639 ierr = CeedVectorCreate(ceedparent, l_size, assembled); CeedChkBackend(ierr); 640 } 641 ierr = CeedVectorSetValue(*assembled, 0.0); CeedChkBackend(ierr); 642 ierr = CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &a); 643 CeedChkBackend(ierr); 644 645 // Input basis apply 646 ierr = CeedOperatorInputBasis_Hip(numelements, qfinputfields, opinputfields, 647 numinputfields, true, edata, impl); 648 CeedChkBackend(ierr); 649 650 // Assemble QFunction 651 for (CeedInt in=0; in<numactivein; in++) { 652 // Set Inputs 653 ierr = CeedVectorSetValue(activein[in], 1.0); CeedChkBackend(ierr); 654 if (numactivein > 1) { 655 ierr = CeedVectorSetValue(activein[(in+numactivein-1)%numactivein], 656 0.0); CeedChkBackend(ierr); 657 } 658 // Set Outputs 659 for (CeedInt out=0; out<numoutputfields; out++) { 660 // Get output vector 661 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 662 CeedChkBackend(ierr); 663 // Check if active output 664 if (vec == CEED_VECTOR_ACTIVE) { 665 CeedVectorSetArray(impl->qvecsout[out], CEED_MEM_DEVICE, 666 CEED_USE_POINTER, a); CeedChkBackend(ierr); 667 ierr = CeedQFunctionFieldGetSize(qfoutputfields[out], &size); 668 CeedChkBackend(ierr); 669 a += size*Q*numelements; // Advance the pointer by the size of the output 670 } 671 } 672 // Apply QFunction 673 ierr = CeedQFunctionApply(qf, Q*numelements, impl->qvecsin, impl->qvecsout); 674 CeedChkBackend(ierr); 675 } 676 677 // Un-set output Qvecs to prevent accidental overwrite of Assembled 678 for (CeedInt out=0; out<numoutputfields; out++) { 679 // Get output vector 680 ierr = CeedOperatorFieldGetVector(opoutputfields[out], &vec); 681 CeedChkBackend(ierr); 682 // Check if active output 683 if (vec == CEED_VECTOR_ACTIVE) { 684 ierr = CeedVectorTakeArray(impl->qvecsout[out], CEED_MEM_DEVICE, NULL); 685 CeedChkBackend(ierr); 686 } 687 } 688 689 // Restore input arrays 690 ierr = CeedOperatorRestoreInputs_Hip(numinputfields, qfinputfields, 691 opinputfields, true, edata, impl); 692 CeedChkBackend(ierr); 693 694 // Restore output 695 ierr = CeedVectorRestoreArray(*assembled, &a); CeedChkBackend(ierr); 696 697 return CEED_ERROR_SUCCESS; 698 } 699 700 //------------------------------------------------------------------------------ 701 // Assemble Linear QFunction 702 //------------------------------------------------------------------------------ 703 static int CeedOperatorLinearAssembleQFunction_Hip(CeedOperator op, 704 CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 705 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, true, assembled, rstr, 706 request); 707 } 708 709 //------------------------------------------------------------------------------ 710 // Assemble Linear QFunction 711 //------------------------------------------------------------------------------ 712 static int CeedOperatorLinearAssembleQFunctionUpdate_Hip(CeedOperator op, 713 CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 714 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, false, &assembled, &rstr, 715 request); 716 } 717 718 //------------------------------------------------------------------------------ 719 // Diagonal assembly kernels 720 //------------------------------------------------------------------------------ 721 // *INDENT-OFF* 722 static const char *diagonalkernels = QUOTE( 723 724 typedef enum { 725 /// Perform no evaluation (either because there is no data or it is already at 726 /// quadrature points) 727 CEED_EVAL_NONE = 0, 728 /// Interpolate from nodes to quadrature points 729 CEED_EVAL_INTERP = 1, 730 /// Evaluate gradients at quadrature points from input in a nodal basis 731 CEED_EVAL_GRAD = 2, 732 /// Evaluate divergence at quadrature points from input in a nodal basis 733 CEED_EVAL_DIV = 4, 734 /// Evaluate curl at quadrature points from input in a nodal basis 735 CEED_EVAL_CURL = 8, 736 /// Using no input, evaluate quadrature weights on the reference element 737 CEED_EVAL_WEIGHT = 16, 738 } CeedEvalMode; 739 740 //------------------------------------------------------------------------------ 741 // Get Basis Emode Pointer 742 //------------------------------------------------------------------------------ 743 extern "C" __device__ void CeedOperatorGetBasisPointer_Hip(const CeedScalar **basisptr, 744 CeedEvalMode emode, const CeedScalar *identity, const CeedScalar *interp, 745 const CeedScalar *grad) { 746 switch (emode) { 747 case CEED_EVAL_NONE: 748 *basisptr = identity; 749 break; 750 case CEED_EVAL_INTERP: 751 *basisptr = interp; 752 break; 753 case CEED_EVAL_GRAD: 754 *basisptr = grad; 755 break; 756 case CEED_EVAL_WEIGHT: 757 case CEED_EVAL_DIV: 758 case CEED_EVAL_CURL: 759 break; // Caught by QF Assembly 760 } 761 } 762 763 //------------------------------------------------------------------------------ 764 // Core code for diagonal assembly 765 //------------------------------------------------------------------------------ 766 __device__ void diagonalCore(const CeedInt nelem, 767 const CeedScalar maxnorm, const bool pointBlock, 768 const CeedScalar *identity, 769 const CeedScalar *interpin, const CeedScalar *gradin, 770 const CeedScalar *interpout, const CeedScalar *gradout, 771 const CeedEvalMode *emodein, const CeedEvalMode *emodeout, 772 const CeedScalar *__restrict__ assembledqfarray, 773 CeedScalar *__restrict__ elemdiagarray) { 774 const int tid = threadIdx.x; // running with P threads, tid is evec node 775 const CeedScalar qfvaluebound = maxnorm*1e-12; 776 777 // Compute the diagonal of B^T D B 778 // Each element 779 for (CeedInt e = blockIdx.x*blockDim.z + threadIdx.z; e < nelem; 780 e += gridDim.x*blockDim.z) { 781 CeedInt dout = -1; 782 // Each basis eval mode pair 783 for (CeedInt eout = 0; eout < NUMEMODEOUT; eout++) { 784 const CeedScalar *bt = NULL; 785 if (emodeout[eout] == CEED_EVAL_GRAD) 786 dout += 1; 787 CeedOperatorGetBasisPointer_Hip(&bt, emodeout[eout], identity, interpout, 788 &gradout[dout*NQPTS*NNODES]); 789 CeedInt din = -1; 790 for (CeedInt ein = 0; ein < NUMEMODEIN; ein++) { 791 const CeedScalar *b = NULL; 792 if (emodein[ein] == CEED_EVAL_GRAD) 793 din += 1; 794 CeedOperatorGetBasisPointer_Hip(&b, emodein[ein], identity, interpin, 795 &gradin[din*NQPTS*NNODES]); 796 // Each component 797 for (CeedInt compOut = 0; compOut < NCOMP; compOut++) { 798 // Each qpoint/node pair 799 if (pointBlock) { 800 // Point Block Diagonal 801 for (CeedInt compIn = 0; compIn < NCOMP; compIn++) { 802 CeedScalar evalue = 0.; 803 for (CeedInt q = 0; q < NQPTS; q++) { 804 const CeedScalar qfvalue = 805 assembledqfarray[((((ein*NCOMP+compIn)*NUMEMODEOUT+eout)* 806 NCOMP+compOut)*nelem+e)*NQPTS+q]; 807 if (abs(qfvalue) > qfvaluebound) 808 evalue += bt[q*NNODES+tid] * qfvalue * b[q*NNODES+tid]; 809 } 810 elemdiagarray[((compOut*NCOMP+compIn)*nelem+e)*NNODES+tid] += evalue; 811 } 812 } else { 813 // Diagonal Only 814 CeedScalar evalue = 0.; 815 for (CeedInt q = 0; q < NQPTS; q++) { 816 const CeedScalar qfvalue = 817 assembledqfarray[((((ein*NCOMP+compOut)*NUMEMODEOUT+eout)* 818 NCOMP+compOut)*nelem+e)*NQPTS+q]; 819 if (abs(qfvalue) > qfvaluebound) 820 evalue += bt[q*NNODES+tid] * qfvalue * b[q*NNODES+tid]; 821 } 822 elemdiagarray[(compOut*nelem+e)*NNODES+tid] += evalue; 823 } 824 } 825 } 826 } 827 } 828 } 829 830 //------------------------------------------------------------------------------ 831 // Linear diagonal 832 //------------------------------------------------------------------------------ 833 extern "C" __global__ void linearDiagonal(const CeedInt nelem, 834 const CeedScalar maxnorm, const CeedScalar *identity, 835 const CeedScalar *interpin, const CeedScalar *gradin, 836 const CeedScalar *interpout, const CeedScalar *gradout, 837 const CeedEvalMode *emodein, const CeedEvalMode *emodeout, 838 const CeedScalar *__restrict__ assembledqfarray, 839 CeedScalar *__restrict__ elemdiagarray) { 840 diagonalCore(nelem, maxnorm, false, identity, interpin, gradin, interpout, 841 gradout, emodein, emodeout, assembledqfarray, elemdiagarray); 842 } 843 844 //------------------------------------------------------------------------------ 845 // Linear point block diagonal 846 //------------------------------------------------------------------------------ 847 extern "C" __global__ void linearPointBlockDiagonal(const CeedInt nelem, 848 const CeedScalar maxnorm, const CeedScalar *identity, 849 const CeedScalar *interpin, const CeedScalar *gradin, 850 const CeedScalar *interpout, const CeedScalar *gradout, 851 const CeedEvalMode *emodein, const CeedEvalMode *emodeout, 852 const CeedScalar *__restrict__ assembledqfarray, 853 CeedScalar *__restrict__ elemdiagarray) { 854 diagonalCore(nelem, maxnorm, true, identity, interpin, gradin, interpout, 855 gradout, emodein, emodeout, assembledqfarray, elemdiagarray); 856 } 857 858 ); 859 // *INDENT-ON* 860 861 //------------------------------------------------------------------------------ 862 // Create point block restriction 863 //------------------------------------------------------------------------------ 864 static int CreatePBRestriction(CeedElemRestriction rstr, 865 CeedElemRestriction *pbRstr) { 866 int ierr; 867 Ceed ceed; 868 ierr = CeedElemRestrictionGetCeed(rstr, &ceed); CeedChkBackend(ierr); 869 const CeedInt *offsets; 870 ierr = CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets); 871 CeedChkBackend(ierr); 872 873 // Expand offsets 874 CeedInt nelem, ncomp, elemsize, compstride, max = 1, *pbOffsets; 875 ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChkBackend(ierr); 876 ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChkBackend(ierr); 877 ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChkBackend(ierr); 878 ierr = CeedElemRestrictionGetCompStride(rstr, &compstride); 879 CeedChkBackend(ierr); 880 CeedInt shift = ncomp; 881 if (compstride != 1) 882 shift *= ncomp; 883 ierr = CeedCalloc(nelem*elemsize, &pbOffsets); CeedChkBackend(ierr); 884 for (CeedInt i = 0; i < nelem*elemsize; i++) { 885 pbOffsets[i] = offsets[i]*shift; 886 if (pbOffsets[i] > max) 887 max = pbOffsets[i]; 888 } 889 890 // Create new restriction 891 ierr = CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp*ncomp, 1, 892 max + ncomp*ncomp, CEED_MEM_HOST, 893 CEED_OWN_POINTER, pbOffsets, pbRstr); 894 CeedChkBackend(ierr); 895 896 // Cleanup 897 ierr = CeedElemRestrictionRestoreOffsets(rstr, &offsets); CeedChkBackend(ierr); 898 899 return CEED_ERROR_SUCCESS; 900 } 901 902 //------------------------------------------------------------------------------ 903 // Assemble diagonal setup 904 //------------------------------------------------------------------------------ 905 static inline int CeedOperatorAssembleDiagonalSetup_Hip(CeedOperator op, 906 const bool pointBlock) { 907 int ierr; 908 Ceed ceed; 909 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 910 CeedQFunction qf; 911 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 912 CeedInt numinputfields, numoutputfields; 913 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 914 CeedChkBackend(ierr); 915 916 // Determine active input basis 917 CeedOperatorField *opfields; 918 CeedQFunctionField *qffields; 919 ierr = CeedOperatorGetFields(op, NULL, &opfields, NULL, NULL); 920 CeedChkBackend(ierr); 921 ierr = CeedQFunctionGetFields(qf, NULL, &qffields, NULL, NULL); 922 CeedChkBackend(ierr); 923 CeedInt numemodein = 0, ncomp = 0, dim = 1; 924 CeedEvalMode *emodein = NULL; 925 CeedBasis basisin = NULL; 926 CeedElemRestriction rstrin = NULL; 927 for (CeedInt i = 0; i < numinputfields; i++) { 928 CeedVector vec; 929 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 930 if (vec == CEED_VECTOR_ACTIVE) { 931 CeedElemRestriction rstr; 932 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisin); CeedChkBackend(ierr); 933 ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChkBackend(ierr); 934 ierr = CeedBasisGetDimension(basisin, &dim); CeedChkBackend(ierr); 935 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 936 CeedChkBackend(ierr); 937 if (rstrin && rstrin != rstr) 938 // LCOV_EXCL_START 939 return CeedError(ceed, CEED_ERROR_BACKEND, 940 "Multi-field non-composite operator diagonal assembly not supported"); 941 // LCOV_EXCL_STOP 942 rstrin = rstr; 943 CeedEvalMode emode; 944 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); 945 CeedChkBackend(ierr); 946 switch (emode) { 947 case CEED_EVAL_NONE: 948 case CEED_EVAL_INTERP: 949 ierr = CeedRealloc(numemodein + 1, &emodein); CeedChkBackend(ierr); 950 emodein[numemodein] = emode; 951 numemodein += 1; 952 break; 953 case CEED_EVAL_GRAD: 954 ierr = CeedRealloc(numemodein + dim, &emodein); CeedChkBackend(ierr); 955 for (CeedInt d = 0; d < dim; d++) 956 emodein[numemodein+d] = emode; 957 numemodein += dim; 958 break; 959 case CEED_EVAL_WEIGHT: 960 case CEED_EVAL_DIV: 961 case CEED_EVAL_CURL: 962 break; // Caught by QF Assembly 963 } 964 } 965 } 966 967 // Determine active output basis 968 ierr = CeedOperatorGetFields(op, NULL, NULL, NULL, &opfields); 969 CeedChkBackend(ierr); 970 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qffields); 971 CeedChkBackend(ierr); 972 CeedInt numemodeout = 0; 973 CeedEvalMode *emodeout = NULL; 974 CeedBasis basisout = NULL; 975 CeedElemRestriction rstrout = NULL; 976 for (CeedInt i = 0; i < numoutputfields; i++) { 977 CeedVector vec; 978 ierr = CeedOperatorFieldGetVector(opfields[i], &vec); CeedChkBackend(ierr); 979 if (vec == CEED_VECTOR_ACTIVE) { 980 CeedElemRestriction rstr; 981 ierr = CeedOperatorFieldGetBasis(opfields[i], &basisout); CeedChkBackend(ierr); 982 ierr = CeedOperatorFieldGetElemRestriction(opfields[i], &rstr); 983 CeedChkBackend(ierr); 984 if (rstrout && rstrout != rstr) 985 // LCOV_EXCL_START 986 return CeedError(ceed, CEED_ERROR_BACKEND, 987 "Multi-field non-composite operator diagonal assembly not supported"); 988 // LCOV_EXCL_STOP 989 rstrout = rstr; 990 CeedEvalMode emode; 991 ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode); CeedChkBackend(ierr); 992 switch (emode) { 993 case CEED_EVAL_NONE: 994 case CEED_EVAL_INTERP: 995 ierr = CeedRealloc(numemodeout + 1, &emodeout); CeedChkBackend(ierr); 996 emodeout[numemodeout] = emode; 997 numemodeout += 1; 998 break; 999 case CEED_EVAL_GRAD: 1000 ierr = CeedRealloc(numemodeout + dim, &emodeout); CeedChkBackend(ierr); 1001 for (CeedInt d = 0; d < dim; d++) 1002 emodeout[numemodeout+d] = emode; 1003 numemodeout += dim; 1004 break; 1005 case CEED_EVAL_WEIGHT: 1006 case CEED_EVAL_DIV: 1007 case CEED_EVAL_CURL: 1008 break; // Caught by QF Assembly 1009 } 1010 } 1011 } 1012 1013 // Operator data struct 1014 CeedOperator_Hip *impl; 1015 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1016 ierr = CeedCalloc(1, &impl->diag); CeedChkBackend(ierr); 1017 CeedOperatorDiag_Hip *diag = impl->diag; 1018 diag->basisin = basisin; 1019 diag->basisout = basisout; 1020 diag->h_emodein = emodein; 1021 diag->h_emodeout = emodeout; 1022 diag->numemodein = numemodein; 1023 diag->numemodeout = numemodeout; 1024 1025 // Assemble kernel 1026 CeedInt nnodes, nqpts; 1027 ierr = CeedBasisGetNumNodes(basisin, &nnodes); CeedChkBackend(ierr); 1028 ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChkBackend(ierr); 1029 diag->nnodes = nnodes; 1030 ierr = CeedCompileHip(ceed, diagonalkernels, &diag->module, 5, 1031 "NUMEMODEIN", numemodein, 1032 "NUMEMODEOUT", numemodeout, 1033 "NNODES", nnodes, 1034 "NQPTS", nqpts, 1035 "NCOMP", ncomp 1036 ); CeedChk_Hip(ceed, ierr); 1037 ierr = CeedGetKernelHip(ceed, diag->module, "linearDiagonal", 1038 &diag->linearDiagonal); CeedChk_Hip(ceed, ierr); 1039 ierr = CeedGetKernelHip(ceed, diag->module, "linearPointBlockDiagonal", 1040 &diag->linearPointBlock); 1041 CeedChk_Hip(ceed, ierr); 1042 1043 // Basis matrices 1044 const CeedInt qBytes = nqpts * sizeof(CeedScalar); 1045 const CeedInt iBytes = qBytes * nnodes; 1046 const CeedInt gBytes = qBytes * nnodes * dim; 1047 const CeedInt eBytes = sizeof(CeedEvalMode); 1048 const CeedScalar *interpin, *interpout, *gradin, *gradout; 1049 1050 // CEED_EVAL_NONE 1051 CeedScalar *identity = NULL; 1052 bool evalNone = false; 1053 for (CeedInt i=0; i<numemodein; i++) 1054 evalNone = evalNone || (emodein[i] == CEED_EVAL_NONE); 1055 for (CeedInt i=0; i<numemodeout; i++) 1056 evalNone = evalNone || (emodeout[i] == CEED_EVAL_NONE); 1057 if (evalNone) { 1058 ierr = CeedCalloc(nqpts*nnodes, &identity); CeedChkBackend(ierr); 1059 for (CeedInt i=0; i<(nnodes<nqpts?nnodes:nqpts); i++) 1060 identity[i*nnodes+i] = 1.0; 1061 ierr = hipMalloc((void **)&diag->d_identity, iBytes); CeedChk_Hip(ceed, ierr); 1062 ierr = hipMemcpy(diag->d_identity, identity, iBytes, 1063 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1064 } 1065 1066 // CEED_EVAL_INTERP 1067 ierr = CeedBasisGetInterp(basisin, &interpin); CeedChkBackend(ierr); 1068 ierr = hipMalloc((void **)&diag->d_interpin, iBytes); CeedChk_Hip(ceed, ierr); 1069 ierr = hipMemcpy(diag->d_interpin, interpin, iBytes, 1070 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1071 ierr = CeedBasisGetInterp(basisout, &interpout); CeedChkBackend(ierr); 1072 ierr = hipMalloc((void **)&diag->d_interpout, iBytes); CeedChk_Hip(ceed, ierr); 1073 ierr = hipMemcpy(diag->d_interpout, interpout, iBytes, 1074 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1075 1076 // CEED_EVAL_GRAD 1077 ierr = CeedBasisGetGrad(basisin, &gradin); CeedChkBackend(ierr); 1078 ierr = hipMalloc((void **)&diag->d_gradin, gBytes); CeedChk_Hip(ceed, ierr); 1079 ierr = hipMemcpy(diag->d_gradin, gradin, gBytes, 1080 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1081 ierr = CeedBasisGetGrad(basisout, &gradout); CeedChkBackend(ierr); 1082 ierr = hipMalloc((void **)&diag->d_gradout, gBytes); CeedChk_Hip(ceed, ierr); 1083 ierr = hipMemcpy(diag->d_gradout, gradout, gBytes, 1084 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1085 1086 // Arrays of emodes 1087 ierr = hipMalloc((void **)&diag->d_emodein, numemodein * eBytes); 1088 CeedChk_Hip(ceed, ierr); 1089 ierr = hipMemcpy(diag->d_emodein, emodein, numemodein * eBytes, 1090 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1091 ierr = hipMalloc((void **)&diag->d_emodeout, numemodeout * eBytes); 1092 CeedChk_Hip(ceed, ierr); 1093 ierr = hipMemcpy(diag->d_emodeout, emodeout, numemodeout * eBytes, 1094 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1095 1096 // Restriction 1097 diag->diagrstr = rstrout; 1098 1099 return CEED_ERROR_SUCCESS; 1100 } 1101 1102 //------------------------------------------------------------------------------ 1103 // Assemble diagonal common code 1104 //------------------------------------------------------------------------------ 1105 static inline int CeedOperatorAssembleDiagonalCore_Hip(CeedOperator op, 1106 CeedVector assembled, CeedRequest *request, const bool pointBlock) { 1107 int ierr; 1108 Ceed ceed; 1109 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1110 CeedOperator_Hip *impl; 1111 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1112 1113 // Assemble QFunction 1114 CeedVector assembledqf; 1115 CeedElemRestriction rstr; 1116 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembledqf, 1117 &rstr, request); CeedChkBackend(ierr); 1118 ierr = CeedElemRestrictionDestroy(&rstr); CeedChkBackend(ierr); 1119 CeedScalar maxnorm = 0; 1120 ierr = CeedVectorNorm(assembledqf, CEED_NORM_MAX, &maxnorm); 1121 CeedChkBackend(ierr); 1122 1123 // Setup 1124 if (!impl->diag) { 1125 ierr = CeedOperatorAssembleDiagonalSetup_Hip(op, pointBlock); 1126 CeedChkBackend(ierr); 1127 } 1128 CeedOperatorDiag_Hip *diag = impl->diag; 1129 assert(diag != NULL); 1130 1131 // Restriction 1132 if (pointBlock && !diag->pbdiagrstr) { 1133 CeedElemRestriction pbdiagrstr; 1134 ierr = CreatePBRestriction(diag->diagrstr, &pbdiagrstr); CeedChkBackend(ierr); 1135 diag->pbdiagrstr = pbdiagrstr; 1136 } 1137 CeedElemRestriction diagrstr = pointBlock ? diag->pbdiagrstr : diag->diagrstr; 1138 1139 // Create diagonal vector 1140 CeedVector elemdiag = pointBlock ? diag->pbelemdiag : diag->elemdiag; 1141 if (!elemdiag) { 1142 // Element diagonal vector 1143 ierr = CeedElemRestrictionCreateVector(diagrstr, NULL, &elemdiag); 1144 CeedChkBackend(ierr); 1145 if (pointBlock) 1146 diag->pbelemdiag = elemdiag; 1147 else 1148 diag->elemdiag = elemdiag; 1149 } 1150 ierr = CeedVectorSetValue(elemdiag, 0.0); CeedChkBackend(ierr); 1151 1152 // Assemble element operator diagonals 1153 CeedScalar *elemdiagarray; 1154 const CeedScalar *assembledqfarray; 1155 ierr = CeedVectorGetArray(elemdiag, CEED_MEM_DEVICE, &elemdiagarray); 1156 CeedChkBackend(ierr); 1157 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_DEVICE, &assembledqfarray); 1158 CeedChkBackend(ierr); 1159 CeedInt nelem; 1160 ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); 1161 CeedChkBackend(ierr); 1162 1163 // Compute the diagonal of B^T D B 1164 int elemsPerBlock = 1; 1165 int grid = nelem/elemsPerBlock+((nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1166 void *args[] = {(void *) &nelem, (void *) &maxnorm, &diag->d_identity, 1167 &diag->d_interpin, &diag->d_gradin, &diag->d_interpout, 1168 &diag->d_gradout, &diag->d_emodein, &diag->d_emodeout, 1169 &assembledqfarray, &elemdiagarray 1170 }; 1171 if (pointBlock) { 1172 ierr = CeedRunKernelDimHip(ceed, diag->linearPointBlock, grid, 1173 diag->nnodes, 1, elemsPerBlock, args); 1174 CeedChkBackend(ierr); 1175 } else { 1176 ierr = CeedRunKernelDimHip(ceed, diag->linearDiagonal, grid, 1177 diag->nnodes, 1, elemsPerBlock, args); 1178 CeedChkBackend(ierr); 1179 } 1180 1181 // Restore arrays 1182 ierr = CeedVectorRestoreArray(elemdiag, &elemdiagarray); CeedChkBackend(ierr); 1183 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1184 CeedChkBackend(ierr); 1185 1186 // Assemble local operator diagonal 1187 ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag, 1188 assembled, request); CeedChkBackend(ierr); 1189 1190 // Cleanup 1191 ierr = CeedVectorDestroy(&assembledqf); CeedChkBackend(ierr); 1192 1193 return CEED_ERROR_SUCCESS; 1194 } 1195 1196 //------------------------------------------------------------------------------ 1197 // Assemble composite diagonal common code 1198 //------------------------------------------------------------------------------ 1199 static inline int CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip( 1200 CeedOperator op, CeedVector assembled, CeedRequest *request, 1201 const bool pointBlock) { 1202 int ierr; 1203 CeedInt numSub; 1204 CeedOperator *subOperators; 1205 ierr = CeedOperatorGetNumSub(op, &numSub); CeedChkBackend(ierr); 1206 ierr = CeedOperatorGetSubList(op, &subOperators); CeedChkBackend(ierr); 1207 for (CeedInt i = 0; i < numSub; i++) { 1208 ierr = CeedOperatorAssembleDiagonalCore_Hip(subOperators[i], assembled, 1209 request, pointBlock); CeedChkBackend(ierr); 1210 } 1211 return CEED_ERROR_SUCCESS; 1212 } 1213 1214 //------------------------------------------------------------------------------ 1215 // Assemble Linear Diagonal 1216 //------------------------------------------------------------------------------ 1217 static int CeedOperatorLinearAssembleAddDiagonal_Hip(CeedOperator op, 1218 CeedVector assembled, CeedRequest *request) { 1219 int ierr; 1220 bool isComposite; 1221 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1222 if (isComposite) { 1223 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip(op, assembled, 1224 request, false); 1225 } else { 1226 return CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, false); 1227 } 1228 } 1229 1230 //------------------------------------------------------------------------------ 1231 // Assemble Linear Point Block Diagonal 1232 //------------------------------------------------------------------------------ 1233 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip(CeedOperator op, 1234 CeedVector assembled, CeedRequest *request) { 1235 int ierr; 1236 bool isComposite; 1237 ierr = CeedOperatorIsComposite(op, &isComposite); CeedChkBackend(ierr); 1238 if (isComposite) { 1239 return CeedOperatorLinearAssembleAddDiagonalCompositeCore_Hip(op, assembled, 1240 request, true); 1241 } else { 1242 return CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, true); 1243 } 1244 } 1245 1246 1247 //------------------------------------------------------------------------------ 1248 // Create operator 1249 //------------------------------------------------------------------------------ 1250 int CeedOperatorCreate_Hip(CeedOperator op) { 1251 int ierr; 1252 Ceed ceed; 1253 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1254 CeedOperator_Hip *impl; 1255 1256 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1257 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1258 1259 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1260 CeedOperatorLinearAssembleQFunction_Hip); 1261 CeedChkBackend(ierr); 1262 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1263 "LinearAssembleQFunctionUpdate", 1264 CeedOperatorLinearAssembleQFunctionUpdate_Hip); 1265 CeedChkBackend(ierr); 1266 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1267 CeedOperatorLinearAssembleAddDiagonal_Hip); 1268 CeedChkBackend(ierr); 1269 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1270 "LinearAssembleAddPointBlockDiagonal", 1271 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1272 CeedChkBackend(ierr); 1273 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1274 CeedOperatorApplyAdd_Hip); CeedChkBackend(ierr); 1275 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1276 CeedOperatorDestroy_Hip); CeedChkBackend(ierr); 1277 return CEED_ERROR_SUCCESS; 1278 } 1279 1280 //------------------------------------------------------------------------------ 1281 // Composite Operator Create 1282 //------------------------------------------------------------------------------ 1283 int CeedCompositeOperatorCreate_Hip(CeedOperator op) { 1284 int ierr; 1285 Ceed ceed; 1286 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1287 1288 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1289 CeedOperatorLinearAssembleAddDiagonal_Hip); 1290 CeedChkBackend(ierr); 1291 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1292 "LinearAssembleAddPointBlockDiagonal", 1293 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1294 CeedChkBackend(ierr); 1295 return CEED_ERROR_SUCCESS; 1296 } 1297 //------------------------------------------------------------------------------ 1298