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, *pbOffsets; 733 CeedSize l_size; 734 ierr = CeedElemRestrictionGetNumElements(rstr, &nelem); CeedChkBackend(ierr); 735 ierr = CeedElemRestrictionGetNumComponents(rstr, &ncomp); CeedChkBackend(ierr); 736 ierr = CeedElemRestrictionGetElementSize(rstr, &elemsize); CeedChkBackend(ierr); 737 ierr = CeedElemRestrictionGetCompStride(rstr, &compstride); 738 CeedChkBackend(ierr); 739 ierr = CeedElemRestrictionGetLVectorSize(rstr, &l_size); CeedChkBackend(ierr); 740 CeedInt shift = ncomp; 741 if (compstride != 1) 742 shift *= ncomp; 743 ierr = CeedCalloc(nelem*elemsize, &pbOffsets); CeedChkBackend(ierr); 744 for (CeedInt i = 0; i < nelem*elemsize; i++) { 745 pbOffsets[i] = offsets[i]*shift; 746 } 747 748 // Create new restriction 749 ierr = CeedElemRestrictionCreate(ceed, nelem, elemsize, ncomp*ncomp, 1, 750 l_size * 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 991 // Setup 992 if (!impl->diag) { 993 ierr = CeedOperatorAssembleDiagonalSetup_Hip(op, pointBlock); 994 CeedChkBackend(ierr); 995 } 996 CeedOperatorDiag_Hip *diag = impl->diag; 997 assert(diag != NULL); 998 999 // Restriction 1000 if (pointBlock && !diag->pbdiagrstr) { 1001 CeedElemRestriction pbdiagrstr; 1002 ierr = CreatePBRestriction(diag->diagrstr, &pbdiagrstr); CeedChkBackend(ierr); 1003 diag->pbdiagrstr = pbdiagrstr; 1004 } 1005 CeedElemRestriction diagrstr = pointBlock ? diag->pbdiagrstr : diag->diagrstr; 1006 1007 // Create diagonal vector 1008 CeedVector elemdiag = pointBlock ? diag->pbelemdiag : diag->elemdiag; 1009 if (!elemdiag) { 1010 // Element diagonal vector 1011 ierr = CeedElemRestrictionCreateVector(diagrstr, NULL, &elemdiag); 1012 CeedChkBackend(ierr); 1013 if (pointBlock) 1014 diag->pbelemdiag = elemdiag; 1015 else 1016 diag->elemdiag = elemdiag; 1017 } 1018 ierr = CeedVectorSetValue(elemdiag, 0.0); CeedChkBackend(ierr); 1019 1020 // Assemble element operator diagonals 1021 CeedScalar *elemdiagarray; 1022 const CeedScalar *assembledqfarray; 1023 ierr = CeedVectorGetArray(elemdiag, CEED_MEM_DEVICE, &elemdiagarray); 1024 CeedChkBackend(ierr); 1025 ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_DEVICE, &assembledqfarray); 1026 CeedChkBackend(ierr); 1027 CeedInt nelem; 1028 ierr = CeedElemRestrictionGetNumElements(diagrstr, &nelem); 1029 CeedChkBackend(ierr); 1030 1031 // Compute the diagonal of B^T D B 1032 int elemsPerBlock = 1; 1033 int grid = nelem/elemsPerBlock+((nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1034 void *args[] = {(void *) &nelem, &diag->d_identity, 1035 &diag->d_interpin, &diag->d_gradin, &diag->d_interpout, 1036 &diag->d_gradout, &diag->d_emodein, &diag->d_emodeout, 1037 &assembledqfarray, &elemdiagarray 1038 }; 1039 if (pointBlock) { 1040 ierr = CeedRunKernelDimHip(ceed, diag->linearPointBlock, grid, 1041 diag->nnodes, 1, elemsPerBlock, args); 1042 CeedChkBackend(ierr); 1043 } else { 1044 ierr = CeedRunKernelDimHip(ceed, diag->linearDiagonal, grid, 1045 diag->nnodes, 1, elemsPerBlock, args); 1046 CeedChkBackend(ierr); 1047 } 1048 1049 // Restore arrays 1050 ierr = CeedVectorRestoreArray(elemdiag, &elemdiagarray); CeedChkBackend(ierr); 1051 ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray); 1052 CeedChkBackend(ierr); 1053 1054 // Assemble local operator diagonal 1055 ierr = CeedElemRestrictionApply(diagrstr, CEED_TRANSPOSE, elemdiag, 1056 assembled, request); CeedChkBackend(ierr); 1057 1058 // Cleanup 1059 ierr = CeedVectorDestroy(&assembledqf); CeedChkBackend(ierr); 1060 1061 return CEED_ERROR_SUCCESS; 1062 } 1063 1064 //------------------------------------------------------------------------------ 1065 // Assemble Linear Diagonal 1066 //------------------------------------------------------------------------------ 1067 static int CeedOperatorLinearAssembleAddDiagonal_Hip(CeedOperator op, 1068 CeedVector assembled, CeedRequest *request) { 1069 int ierr = CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, false); 1070 CeedChkBackend(ierr); 1071 return CEED_ERROR_SUCCESS; 1072 } 1073 1074 //------------------------------------------------------------------------------ 1075 // Assemble Linear Point Block Diagonal 1076 //------------------------------------------------------------------------------ 1077 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip(CeedOperator op, 1078 CeedVector assembled, CeedRequest *request) { 1079 int ierr = CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, true); 1080 CeedChkBackend(ierr); 1081 return CEED_ERROR_SUCCESS; 1082 } 1083 1084 //------------------------------------------------------------------------------ 1085 // Single operator assembly setup 1086 //------------------------------------------------------------------------------ 1087 static int CeedSingleOperatorAssembleSetup_Hip(CeedOperator op) { 1088 int ierr; 1089 Ceed ceed; 1090 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1091 CeedOperator_Hip *impl; 1092 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1093 1094 // Get intput and output fields 1095 CeedInt num_input_fields, num_output_fields; 1096 CeedOperatorField *input_fields; 1097 CeedOperatorField *output_fields; 1098 ierr = CeedOperatorGetFields(op, &num_input_fields, &input_fields, 1099 &num_output_fields, &output_fields); CeedChkBackend(ierr); 1100 1101 // Determine active input basis eval mode 1102 CeedQFunction qf; 1103 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 1104 CeedQFunctionField *qf_fields; 1105 ierr = CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL); 1106 CeedChkBackend(ierr); 1107 // Note that the kernel will treat each dimension of a gradient action separately; 1108 // i.e., when an active input has a CEED_EVAL_GRAD mode, num_emode_in will increment 1109 // by dim. However, for the purposes of loading the B matrices, it will be treated 1110 // as one mode, and we will load/copy the entire gradient matrix at once, so 1111 // num_B_in_mats_to_load will be incremented by 1. 1112 CeedInt num_emode_in = 0, dim = 1, num_B_in_mats_to_load = 0, size_B_in = 0; 1113 CeedEvalMode *eval_mode_in = NULL; //will be of size num_B_in_mats_load 1114 CeedBasis basis_in = NULL; 1115 CeedInt nqpts = 0, esize = 0; 1116 CeedElemRestriction rstr_in = NULL; 1117 for (CeedInt i=0; i<num_input_fields; i++) { 1118 CeedVector vec; 1119 ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChkBackend(ierr); 1120 if (vec == CEED_VECTOR_ACTIVE) { 1121 ierr = CeedOperatorFieldGetBasis(input_fields[i], &basis_in); 1122 CeedChkBackend(ierr); 1123 ierr = CeedBasisGetDimension(basis_in, &dim); CeedChkBackend(ierr); 1124 ierr = CeedBasisGetNumQuadraturePoints(basis_in, &nqpts); CeedChkBackend(ierr); 1125 ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in); 1126 CeedChkBackend(ierr); 1127 ierr = CeedElemRestrictionGetElementSize(rstr_in, &esize); CeedChkBackend(ierr); 1128 CeedEvalMode eval_mode; 1129 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1130 CeedChkBackend(ierr); 1131 if (eval_mode != CEED_EVAL_NONE) { 1132 ierr = CeedRealloc(num_B_in_mats_to_load + 1, &eval_mode_in); 1133 CeedChkBackend(ierr); 1134 eval_mode_in[num_B_in_mats_to_load] = eval_mode; 1135 num_B_in_mats_to_load += 1; 1136 if (eval_mode == CEED_EVAL_GRAD) { 1137 num_emode_in += dim; 1138 size_B_in += dim * esize * nqpts; 1139 } else { 1140 num_emode_in +=1; 1141 size_B_in += esize * nqpts; 1142 } 1143 } 1144 } 1145 } 1146 1147 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1148 ierr = CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields); 1149 CeedChkBackend(ierr); 1150 CeedInt num_emode_out = 0, num_B_out_mats_to_load = 0, size_B_out = 0; 1151 CeedEvalMode *eval_mode_out = NULL; 1152 CeedBasis basis_out = NULL; 1153 CeedElemRestriction rstr_out = NULL; 1154 for (CeedInt i=0; i<num_output_fields; i++) { 1155 CeedVector vec; 1156 ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChkBackend(ierr); 1157 if (vec == CEED_VECTOR_ACTIVE) { 1158 ierr = CeedOperatorFieldGetBasis(output_fields[i], &basis_out); 1159 CeedChkBackend(ierr); 1160 ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out); 1161 CeedChkBackend(ierr); 1162 if (rstr_out && rstr_out != rstr_in) 1163 // LCOV_EXCL_START 1164 return CeedError(ceed, CEED_ERROR_BACKEND, 1165 "Multi-field non-composite operator assembly not supported"); 1166 // LCOV_EXCL_STOP 1167 CeedEvalMode eval_mode; 1168 ierr = CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode); 1169 CeedChkBackend(ierr); 1170 if (eval_mode != CEED_EVAL_NONE) { 1171 ierr = CeedRealloc(num_B_out_mats_to_load + 1, &eval_mode_out); 1172 CeedChkBackend(ierr); 1173 eval_mode_out[num_B_out_mats_to_load] = eval_mode; 1174 num_B_out_mats_to_load += 1; 1175 if (eval_mode == CEED_EVAL_GRAD) { 1176 num_emode_out += dim; 1177 size_B_out += dim * esize * nqpts; 1178 } else { 1179 num_emode_out +=1; 1180 size_B_out += esize * nqpts; 1181 } 1182 } 1183 } 1184 } 1185 1186 if (num_emode_in == 0 || num_emode_out == 0) 1187 // LCOV_EXCL_START 1188 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1189 "Cannot assemble operator without inputs/outputs"); 1190 // LCOV_EXCL_STOP 1191 1192 CeedInt nelem, ncomp; 1193 ierr = CeedElemRestrictionGetNumElements(rstr_in, &nelem); CeedChkBackend(ierr); 1194 ierr = CeedElemRestrictionGetNumComponents(rstr_in, &ncomp); 1195 CeedChkBackend(ierr); 1196 1197 ierr = CeedCalloc(1, &impl->asmb); CeedChkBackend(ierr); 1198 CeedOperatorAssemble_Hip *asmb = impl->asmb; 1199 asmb->nelem = nelem; 1200 1201 // Compile kernels 1202 int elemsPerBlock = 1; 1203 asmb->elemsPerBlock = elemsPerBlock; 1204 CeedInt block_size = esize * esize * elemsPerBlock; 1205 char *assembly_kernel_path, *assembly_kernel_source; 1206 ierr = CeedGetJitAbsolutePath(ceed, 1207 "ceed/jit-source/hip/hip-ref-operator-assemble.h", 1208 &assembly_kernel_path); CeedChkBackend(ierr); 1209 CeedDebug256(ceed, 2, "----- Loading Assembly Kernel Source -----\n"); 1210 ierr = CeedLoadSourceToBuffer(ceed, assembly_kernel_path, 1211 &assembly_kernel_source); 1212 CeedChkBackend(ierr); 1213 CeedDebug256(ceed, 2, "----- Loading Assembly Source Complete! -----\n"); 1214 bool fallback = block_size > 1024; 1215 if (fallback) { // Use fallback kernel with 1D threadblock 1216 block_size = esize * elemsPerBlock; 1217 asmb->block_size_x = esize; 1218 asmb->block_size_y = 1; 1219 } else { // Use kernel with 2D threadblock 1220 asmb->block_size_x = esize; 1221 asmb->block_size_y = esize; 1222 } 1223 ierr = CeedCompileHip(ceed, assembly_kernel_source, &asmb->module, 7, 1224 "NELEM", nelem, 1225 "NUMEMODEIN", num_emode_in, 1226 "NUMEMODEOUT", num_emode_out, 1227 "NQPTS", nqpts, 1228 "NNODES", esize, 1229 "BLOCK_SIZE", block_size, 1230 "NCOMP", ncomp 1231 ); CeedChk_Hip(ceed, ierr); 1232 ierr = CeedGetKernelHip(ceed, asmb->module, 1233 fallback ? "linearAssembleFallback" : "linearAssemble", 1234 &asmb->linearAssemble); CeedChk_Hip(ceed, ierr); 1235 ierr = CeedFree(&assembly_kernel_path); CeedChkBackend(ierr); 1236 ierr = CeedFree(&assembly_kernel_source); CeedChkBackend(ierr); 1237 1238 // Build 'full' B matrices (not 1D arrays used for tensor-product matrices) 1239 const CeedScalar *interp_in, *grad_in; 1240 ierr = CeedBasisGetInterp(basis_in, &interp_in); CeedChkBackend(ierr); 1241 ierr = CeedBasisGetGrad(basis_in, &grad_in); CeedChkBackend(ierr); 1242 1243 // Load into B_in, in order that they will be used in eval_mode 1244 const CeedInt inBytes = size_B_in * sizeof(CeedScalar); 1245 CeedInt mat_start = 0; 1246 ierr = hipMalloc((void **) &asmb->d_B_in, inBytes); CeedChk_Hip(ceed, ierr); 1247 for (int i = 0; i < num_B_in_mats_to_load; i++) { 1248 CeedEvalMode eval_mode = eval_mode_in[i]; 1249 if (eval_mode == CEED_EVAL_INTERP) { 1250 ierr = hipMemcpy(&asmb->d_B_in[mat_start], interp_in, 1251 esize * nqpts * sizeof(CeedScalar), 1252 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1253 mat_start += esize * nqpts; 1254 } else if (eval_mode == CEED_EVAL_GRAD) { 1255 ierr = hipMemcpy(&asmb->d_B_in[mat_start], grad_in, 1256 dim * esize * nqpts * sizeof(CeedScalar), 1257 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1258 mat_start += dim * esize * nqpts; 1259 } 1260 } 1261 1262 const CeedScalar *interp_out, *grad_out; 1263 // Note that this function currently assumes 1 basis, so this should always be true 1264 // for now 1265 if (basis_out == basis_in) { 1266 interp_out = interp_in; 1267 grad_out = grad_in; 1268 } else { 1269 ierr = CeedBasisGetInterp(basis_out, &interp_out); CeedChkBackend(ierr); 1270 ierr = CeedBasisGetGrad(basis_out, &grad_out); CeedChkBackend(ierr); 1271 } 1272 1273 // Load into B_out, in order that they will be used in eval_mode 1274 const CeedInt outBytes = size_B_out * sizeof(CeedScalar); 1275 mat_start = 0; 1276 ierr = hipMalloc((void **) &asmb->d_B_out, outBytes); CeedChk_Hip(ceed, ierr); 1277 for (int i = 0; i < num_B_out_mats_to_load; i++) { 1278 CeedEvalMode eval_mode = eval_mode_out[i]; 1279 if (eval_mode == CEED_EVAL_INTERP) { 1280 ierr = hipMemcpy(&asmb->d_B_out[mat_start], interp_out, 1281 esize * nqpts * sizeof(CeedScalar), 1282 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1283 mat_start += esize * nqpts; 1284 } else if (eval_mode == CEED_EVAL_GRAD) { 1285 ierr = hipMemcpy(&asmb->d_B_out[mat_start], grad_out, 1286 dim * esize * nqpts * sizeof(CeedScalar), 1287 hipMemcpyHostToDevice); CeedChk_Hip(ceed, ierr); 1288 mat_start += dim * esize * nqpts; 1289 } 1290 } 1291 return CEED_ERROR_SUCCESS; 1292 } 1293 1294 //------------------------------------------------------------------------------ 1295 // Assemble matrix data for COO matrix of assembled operator. 1296 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1297 // 1298 // Note that this (and other assembly routines) currently assume only one 1299 // active input restriction/basis per operator (could have multiple basis eval 1300 // modes). 1301 // TODO: allow multiple active input restrictions/basis objects 1302 //------------------------------------------------------------------------------ 1303 static int CeedSingleOperatorAssemble_Hip(CeedOperator op, CeedInt offset, 1304 CeedVector values) { 1305 1306 int ierr; 1307 Ceed ceed; 1308 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1309 CeedOperator_Hip *impl; 1310 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 1311 1312 // Setup 1313 if (!impl->asmb) { 1314 ierr = CeedSingleOperatorAssembleSetup_Hip(op); 1315 CeedChkBackend(ierr); 1316 assert(impl->asmb != NULL); 1317 } 1318 1319 // Assemble QFunction 1320 CeedVector assembled_qf; 1321 CeedElemRestriction rstr_q; 1322 ierr = CeedOperatorLinearAssembleQFunctionBuildOrUpdate( 1323 op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChkBackend(ierr); 1324 ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChkBackend(ierr); 1325 CeedScalar *values_array; 1326 ierr = CeedVectorGetArrayWrite(values, CEED_MEM_DEVICE, &values_array); 1327 CeedChkBackend(ierr); 1328 values_array += offset; 1329 const CeedScalar *qf_array; 1330 ierr = CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &qf_array); 1331 CeedChkBackend(ierr); 1332 1333 // Compute B^T D B 1334 const CeedInt nelem = impl->asmb->nelem; // to satisfy clang-tidy 1335 const CeedInt elemsPerBlock = impl->asmb->elemsPerBlock; 1336 const CeedInt grid = nelem/elemsPerBlock+(( 1337 nelem/elemsPerBlock*elemsPerBlock<nelem)?1:0); 1338 void *args[] = {&impl->asmb->d_B_in, &impl->asmb->d_B_out, 1339 &qf_array, &values_array 1340 }; 1341 ierr = CeedRunKernelDimHip(ceed, impl->asmb->linearAssemble, grid, 1342 impl->asmb->block_size_x, impl->asmb->block_size_y, 1343 elemsPerBlock, args); 1344 CeedChkBackend(ierr); 1345 1346 1347 // Restore arrays 1348 ierr = CeedVectorRestoreArray(values, &values_array); CeedChkBackend(ierr); 1349 ierr = CeedVectorRestoreArrayRead(assembled_qf, &qf_array); 1350 CeedChkBackend(ierr); 1351 1352 // Cleanup 1353 ierr = CeedVectorDestroy(&assembled_qf); CeedChkBackend(ierr); 1354 1355 return CEED_ERROR_SUCCESS; 1356 } 1357 1358 //------------------------------------------------------------------------------ 1359 // Create operator 1360 //------------------------------------------------------------------------------ 1361 int CeedOperatorCreate_Hip(CeedOperator op) { 1362 int ierr; 1363 Ceed ceed; 1364 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 1365 CeedOperator_Hip *impl; 1366 1367 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 1368 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 1369 1370 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", 1371 CeedOperatorLinearAssembleQFunction_Hip); 1372 CeedChkBackend(ierr); 1373 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1374 "LinearAssembleQFunctionUpdate", 1375 CeedOperatorLinearAssembleQFunctionUpdate_Hip); 1376 CeedChkBackend(ierr); 1377 ierr = CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", 1378 CeedOperatorLinearAssembleAddDiagonal_Hip); 1379 CeedChkBackend(ierr); 1380 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1381 "LinearAssembleAddPointBlockDiagonal", 1382 CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip); 1383 CeedChkBackend(ierr); 1384 ierr = CeedSetBackendFunction(ceed, "Operator", op, 1385 "LinearAssembleSingle", CeedSingleOperatorAssemble_Hip); 1386 CeedChkBackend(ierr); 1387 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 1388 CeedOperatorApplyAdd_Hip); CeedChkBackend(ierr); 1389 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 1390 CeedOperatorDestroy_Hip); CeedChkBackend(ierr); 1391 return CEED_ERROR_SUCCESS; 1392 } 1393 1394 //------------------------------------------------------------------------------ 1395