1 // Copyright (c) 2017-2024, 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.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-tools.h> 11 #include <assert.h> 12 #include <stdbool.h> 13 #include <string.h> 14 #include <hip/hip_runtime.h> 15 16 #include "../hip/ceed-hip-common.h" 17 #include "../hip/ceed-hip-compile.h" 18 #include "ceed-hip-ref.h" 19 20 //------------------------------------------------------------------------------ 21 // Destroy operator 22 //------------------------------------------------------------------------------ 23 static int CeedOperatorDestroy_Hip(CeedOperator op) { 24 CeedOperator_Hip *impl; 25 26 CeedCallBackend(CeedOperatorGetData(op, &impl)); 27 28 // Apply data 29 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 30 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs[i])); 31 } 32 CeedCallBackend(CeedFree(&impl->e_vecs)); 33 CeedCallBackend(CeedFree(&impl->input_states)); 34 35 for (CeedInt i = 0; i < impl->num_inputs; i++) { 36 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 37 } 38 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 39 40 for (CeedInt i = 0; i < impl->num_outputs; i++) { 41 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 42 } 43 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 44 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 45 46 // QFunction assembly data 47 for (CeedInt i = 0; i < impl->num_active_in; i++) { 48 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 49 } 50 CeedCallBackend(CeedFree(&impl->qf_active_in)); 51 52 // Diag data 53 if (impl->diag) { 54 Ceed ceed; 55 56 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 57 if (impl->diag->module) { 58 CeedCallHip(ceed, hipModuleUnload(impl->diag->module)); 59 } 60 if (impl->diag->module_point_block) { 61 CeedCallHip(ceed, hipModuleUnload(impl->diag->module_point_block)); 62 } 63 CeedCallHip(ceed, hipFree(impl->diag->d_eval_modes_in)); 64 CeedCallHip(ceed, hipFree(impl->diag->d_eval_modes_out)); 65 CeedCallHip(ceed, hipFree(impl->diag->d_identity)); 66 CeedCallHip(ceed, hipFree(impl->diag->d_interp_in)); 67 CeedCallHip(ceed, hipFree(impl->diag->d_interp_out)); 68 CeedCallHip(ceed, hipFree(impl->diag->d_grad_in)); 69 CeedCallHip(ceed, hipFree(impl->diag->d_grad_out)); 70 CeedCallHip(ceed, hipFree(impl->diag->d_div_in)); 71 CeedCallHip(ceed, hipFree(impl->diag->d_div_out)); 72 CeedCallHip(ceed, hipFree(impl->diag->d_curl_in)); 73 CeedCallHip(ceed, hipFree(impl->diag->d_curl_out)); 74 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->diag_rstr)); 75 CeedCallBackend(CeedElemRestrictionDestroy(&impl->diag->point_block_diag_rstr)); 76 CeedCallBackend(CeedVectorDestroy(&impl->diag->elem_diag)); 77 CeedCallBackend(CeedVectorDestroy(&impl->diag->point_block_elem_diag)); 78 } 79 CeedCallBackend(CeedFree(&impl->diag)); 80 81 if (impl->asmb) { 82 Ceed ceed; 83 84 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 85 CeedCallHip(ceed, hipModuleUnload(impl->asmb->module)); 86 CeedCallHip(ceed, hipFree(impl->asmb->d_B_in)); 87 CeedCallHip(ceed, hipFree(impl->asmb->d_B_out)); 88 } 89 CeedCallBackend(CeedFree(&impl->asmb)); 90 91 CeedCallBackend(CeedFree(&impl)); 92 return CEED_ERROR_SUCCESS; 93 } 94 95 //------------------------------------------------------------------------------ 96 // Setup infields or outfields 97 //------------------------------------------------------------------------------ 98 static int CeedOperatorSetupFields_Hip(CeedQFunction qf, CeedOperator op, bool is_input, bool is_at_points, CeedVector *e_vecs, CeedVector *q_vecs, 99 CeedInt start_e, CeedInt num_fields, CeedInt Q, CeedInt num_elem) { 100 Ceed ceed; 101 CeedQFunctionField *qf_fields; 102 CeedOperatorField *op_fields; 103 104 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 105 if (is_input) { 106 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 107 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 108 } else { 109 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 110 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 111 } 112 113 // Loop over fields 114 for (CeedInt i = 0; i < num_fields; i++) { 115 bool is_strided = false, skip_restriction = false; 116 CeedSize q_size; 117 CeedInt size; 118 CeedEvalMode eval_mode; 119 CeedBasis basis; 120 121 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 122 if (eval_mode != CEED_EVAL_WEIGHT) { 123 CeedElemRestriction elem_rstr; 124 125 // Check whether this field can skip the element restriction: 126 // Must be passive input, with eval_mode NONE, and have a strided restriction with CEED_STRIDES_BACKEND. 127 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 128 129 // First, check whether the field is input or output: 130 if (is_input) { 131 CeedVector vec; 132 133 // Check for passive input 134 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 135 if (vec != CEED_VECTOR_ACTIVE) { 136 // Check eval_mode 137 if (eval_mode == CEED_EVAL_NONE) { 138 // Check for strided restriction 139 CeedCallBackend(CeedElemRestrictionIsStrided(elem_rstr, &is_strided)); 140 if (is_strided) { 141 // Check if vector is already in preferred backend ordering 142 CeedCallBackend(CeedElemRestrictionHasBackendStrides(elem_rstr, &skip_restriction)); 143 } 144 } 145 } 146 } 147 if (skip_restriction) { 148 // We do not need an E-Vector, but will use the input field vector's data directly in the operator application. 149 e_vecs[i + start_e] = NULL; 150 } else { 151 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs[i + start_e])); 152 } 153 } 154 155 switch (eval_mode) { 156 case CEED_EVAL_NONE: 157 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 158 q_size = (CeedSize)num_elem * Q * size; 159 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 160 break; 161 case CEED_EVAL_INTERP: 162 case CEED_EVAL_GRAD: 163 case CEED_EVAL_DIV: 164 case CEED_EVAL_CURL: 165 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 166 q_size = (CeedSize)num_elem * Q * size; 167 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 168 break; 169 case CEED_EVAL_WEIGHT: // Only on input fields 170 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 171 q_size = (CeedSize)num_elem * Q; 172 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 173 if (is_at_points) { 174 CeedInt num_points[num_elem]; 175 176 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = Q; 177 CeedCallBackend( 178 CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 179 } else { 180 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 181 } 182 break; 183 } 184 } 185 return CEED_ERROR_SUCCESS; 186 } 187 188 //------------------------------------------------------------------------------ 189 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction. 190 //------------------------------------------------------------------------------ 191 static int CeedOperatorSetup_Hip(CeedOperator op) { 192 Ceed ceed; 193 bool is_setup_done; 194 CeedInt Q, num_elem, num_input_fields, num_output_fields; 195 CeedQFunctionField *qf_input_fields, *qf_output_fields; 196 CeedQFunction qf; 197 CeedOperatorField *op_input_fields, *op_output_fields; 198 CeedOperator_Hip *impl; 199 200 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 201 if (is_setup_done) return CEED_ERROR_SUCCESS; 202 203 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 204 CeedCallBackend(CeedOperatorGetData(op, &impl)); 205 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 206 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 207 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 208 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 209 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 210 211 // Allocate 212 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs)); 213 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 214 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 215 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 216 impl->num_inputs = num_input_fields; 217 impl->num_outputs = num_output_fields; 218 219 // Set up infield and outfield e_vecs and q_vecs 220 // Infields 221 CeedCallBackend(CeedOperatorSetupFields_Hip(qf, op, true, false, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, Q, num_elem)); 222 // Outfields 223 CeedCallBackend( 224 CeedOperatorSetupFields_Hip(qf, op, false, false, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, Q, num_elem)); 225 226 CeedCallBackend(CeedOperatorSetSetupDone(op)); 227 return CEED_ERROR_SUCCESS; 228 } 229 230 //------------------------------------------------------------------------------ 231 // Setup Operator Inputs 232 //------------------------------------------------------------------------------ 233 static inline int CeedOperatorSetupInputs_Hip(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 234 CeedVector in_vec, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 235 CeedOperator_Hip *impl, CeedRequest *request) { 236 for (CeedInt i = 0; i < num_input_fields; i++) { 237 CeedEvalMode eval_mode; 238 CeedVector vec; 239 CeedElemRestriction elem_rstr; 240 241 // Get input vector 242 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 243 if (vec == CEED_VECTOR_ACTIVE) { 244 if (skip_active) continue; 245 else vec = in_vec; 246 } 247 248 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 249 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 250 } else { 251 // Get input vector 252 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 253 // Get input element restriction 254 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 255 if (vec == CEED_VECTOR_ACTIVE) vec = in_vec; 256 // Restrict, if necessary 257 if (!impl->e_vecs[i]) { 258 // No restriction for this field; read data directly from vec. 259 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 260 } else { 261 uint64_t state; 262 263 CeedCallBackend(CeedVectorGetState(vec, &state)); 264 if (state != impl->input_states[i]) { 265 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs[i], request)); 266 impl->input_states[i] = state; 267 } 268 // Get evec 269 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs[i], CEED_MEM_DEVICE, (const CeedScalar **)&e_data[i])); 270 } 271 } 272 } 273 return CEED_ERROR_SUCCESS; 274 } 275 276 //------------------------------------------------------------------------------ 277 // Input Basis Action 278 //------------------------------------------------------------------------------ 279 static inline int CeedOperatorInputBasis_Hip(CeedInt num_elem, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 280 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 281 CeedOperator_Hip *impl) { 282 for (CeedInt i = 0; i < num_input_fields; i++) { 283 CeedInt elem_size, size; 284 CeedEvalMode eval_mode; 285 CeedElemRestriction elem_rstr; 286 CeedBasis basis; 287 288 // Skip active input 289 if (skip_active) { 290 CeedVector vec; 291 292 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 293 if (vec == CEED_VECTOR_ACTIVE) continue; 294 } 295 // Get elem_size, eval_mode, size 296 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 297 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 298 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 299 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 300 // Basis action 301 switch (eval_mode) { 302 case CEED_EVAL_NONE: 303 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 304 break; 305 case CEED_EVAL_INTERP: 306 case CEED_EVAL_GRAD: 307 case CEED_EVAL_DIV: 308 case CEED_EVAL_CURL: 309 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 310 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs[i], impl->q_vecs_in[i])); 311 break; 312 case CEED_EVAL_WEIGHT: 313 break; // No action 314 } 315 } 316 return CEED_ERROR_SUCCESS; 317 } 318 319 //------------------------------------------------------------------------------ 320 // Restore Input Vectors 321 //------------------------------------------------------------------------------ 322 static inline int CeedOperatorRestoreInputs_Hip(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 323 const bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Hip *impl) { 324 for (CeedInt i = 0; i < num_input_fields; i++) { 325 CeedEvalMode eval_mode; 326 CeedVector vec; 327 328 // Skip active input 329 if (skip_active) { 330 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 331 if (vec == CEED_VECTOR_ACTIVE) continue; 332 } 333 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 334 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 335 } else { 336 if (!impl->e_vecs[i]) { // This was a skip_restriction case 337 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 338 CeedCallBackend(CeedVectorRestoreArrayRead(vec, (const CeedScalar **)&e_data[i])); 339 } else { 340 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs[i], (const CeedScalar **)&e_data[i])); 341 } 342 } 343 } 344 return CEED_ERROR_SUCCESS; 345 } 346 347 //------------------------------------------------------------------------------ 348 // Apply and add to output 349 //------------------------------------------------------------------------------ 350 static int CeedOperatorApplyAdd_Hip(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 351 CeedInt Q, num_elem, elem_size, num_input_fields, num_output_fields, size; 352 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 353 CeedQFunctionField *qf_input_fields, *qf_output_fields; 354 CeedQFunction qf; 355 CeedOperatorField *op_input_fields, *op_output_fields; 356 CeedOperator_Hip *impl; 357 358 CeedCallBackend(CeedOperatorGetData(op, &impl)); 359 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 360 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 361 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 362 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 363 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 364 365 // Setup 366 CeedCallBackend(CeedOperatorSetup_Hip(op)); 367 368 // Input Evecs and Restriction 369 CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request)); 370 371 // Input basis apply if needed 372 CeedCallBackend(CeedOperatorInputBasis_Hip(num_elem, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl)); 373 374 // Output pointers, as necessary 375 for (CeedInt i = 0; i < num_output_fields; i++) { 376 CeedEvalMode eval_mode; 377 378 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 379 if (eval_mode == CEED_EVAL_NONE) { 380 // Set the output Q-Vector to use the E-Vector data directly. 381 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 382 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 383 } 384 } 385 386 // Q function 387 CeedCallBackend(CeedQFunctionApply(qf, num_elem * Q, impl->q_vecs_in, impl->q_vecs_out)); 388 389 // Output basis apply if needed 390 for (CeedInt i = 0; i < num_output_fields; i++) { 391 CeedEvalMode eval_mode; 392 CeedElemRestriction elem_rstr; 393 CeedBasis basis; 394 395 // Get elem_size, eval_mode, size 396 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 397 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 398 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 399 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 400 // Basis action 401 switch (eval_mode) { 402 case CEED_EVAL_NONE: 403 break; // No action 404 case CEED_EVAL_INTERP: 405 case CEED_EVAL_GRAD: 406 case CEED_EVAL_DIV: 407 case CEED_EVAL_CURL: 408 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 409 CeedCallBackend(CeedBasisApply(basis, num_elem, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs[i + impl->num_inputs])); 410 break; 411 // LCOV_EXCL_START 412 case CEED_EVAL_WEIGHT: { 413 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 414 // LCOV_EXCL_STOP 415 } 416 } 417 } 418 419 // Output restriction 420 for (CeedInt i = 0; i < num_output_fields; i++) { 421 CeedEvalMode eval_mode; 422 CeedVector vec; 423 CeedElemRestriction elem_rstr; 424 425 // Restore evec 426 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 427 if (eval_mode == CEED_EVAL_NONE) { 428 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 429 } 430 // Get output vector 431 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 432 // Restrict 433 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 434 // Active 435 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 436 437 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request)); 438 } 439 440 // Restore input arrays 441 CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl)); 442 return CEED_ERROR_SUCCESS; 443 } 444 445 //------------------------------------------------------------------------------ 446 // CeedOperator needs to connect all the named fields (be they active or passive) to the named inputs and outputs of its CeedQFunction. 447 //------------------------------------------------------------------------------ 448 static int CeedOperatorSetupAtPoints_Hip(CeedOperator op) { 449 Ceed ceed; 450 bool is_setup_done; 451 CeedInt max_num_points = -1, num_elem, num_input_fields, num_output_fields; 452 CeedQFunctionField *qf_input_fields, *qf_output_fields; 453 CeedQFunction qf; 454 CeedOperatorField *op_input_fields, *op_output_fields; 455 CeedOperator_Hip *impl; 456 457 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 458 if (is_setup_done) return CEED_ERROR_SUCCESS; 459 460 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 461 CeedCallBackend(CeedOperatorGetData(op, &impl)); 462 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 463 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 464 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 465 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 466 { 467 CeedElemRestriction elem_rstr = NULL; 468 469 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &elem_rstr, NULL)); 470 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(elem_rstr, &max_num_points)); 471 } 472 impl->max_num_points = max_num_points; 473 474 // Allocate 475 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs)); 476 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 477 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 478 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 479 impl->num_inputs = num_input_fields; 480 impl->num_outputs = num_output_fields; 481 482 // Set up infield and outfield e_vecs and q_vecs 483 // Infields 484 CeedCallBackend(CeedOperatorSetupFields_Hip(qf, op, true, true, impl->e_vecs, impl->q_vecs_in, 0, num_input_fields, max_num_points, num_elem)); 485 // Outfields 486 CeedCallBackend(CeedOperatorSetupFields_Hip(qf, op, false, true, impl->e_vecs, impl->q_vecs_out, num_input_fields, num_output_fields, 487 max_num_points, num_elem)); 488 489 CeedCallBackend(CeedOperatorSetSetupDone(op)); 490 return CEED_ERROR_SUCCESS; 491 } 492 493 //------------------------------------------------------------------------------ 494 // Input Basis Action AtPoints 495 //------------------------------------------------------------------------------ 496 static inline int CeedOperatorInputBasisAtPoints_Hip(CeedInt num_elem, const CeedInt *num_points, CeedQFunctionField *qf_input_fields, 497 CeedOperatorField *op_input_fields, CeedInt num_input_fields, const bool skip_active, 498 CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Hip *impl) { 499 for (CeedInt i = 0; i < num_input_fields; i++) { 500 CeedInt elem_size, size; 501 CeedEvalMode eval_mode; 502 CeedElemRestriction elem_rstr; 503 CeedBasis basis; 504 505 // Skip active input 506 if (skip_active) { 507 CeedVector vec; 508 509 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 510 if (vec == CEED_VECTOR_ACTIVE) continue; 511 } 512 // Get elem_size, eval_mode, size 513 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 514 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 515 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 516 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 517 // Basis action 518 switch (eval_mode) { 519 case CEED_EVAL_NONE: 520 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 521 break; 522 case CEED_EVAL_INTERP: 523 case CEED_EVAL_GRAD: 524 case CEED_EVAL_DIV: 525 case CEED_EVAL_CURL: 526 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 527 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs[i], 528 impl->q_vecs_in[i])); 529 break; 530 case CEED_EVAL_WEIGHT: 531 break; // No action 532 } 533 } 534 return CEED_ERROR_SUCCESS; 535 } 536 537 //------------------------------------------------------------------------------ 538 // Apply and add to output AtPoints 539 //------------------------------------------------------------------------------ 540 static int CeedOperatorApplyAddAtPoints_Hip(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 541 CeedInt max_num_points, num_elem, elem_size, num_input_fields, num_output_fields, size; 542 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 543 CeedQFunctionField *qf_input_fields, *qf_output_fields; 544 CeedQFunction qf; 545 CeedOperatorField *op_input_fields, *op_output_fields; 546 CeedOperator_Hip *impl; 547 548 CeedCallBackend(CeedOperatorGetData(op, &impl)); 549 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 550 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 551 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 552 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 553 CeedInt num_points[num_elem]; 554 555 // Setup 556 CeedCallBackend(CeedOperatorSetupAtPoints_Hip(op)); 557 max_num_points = impl->max_num_points; 558 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = max_num_points; 559 560 // Input Evecs and Restriction 561 CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data, impl, request)); 562 563 // Get point coordinates 564 if (!impl->point_coords_elem) { 565 CeedVector point_coords = NULL; 566 CeedElemRestriction rstr_points = NULL; 567 568 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 569 CeedCallBackend(CeedElemRestrictionCreateVector(rstr_points, NULL, &impl->point_coords_elem)); 570 CeedCallBackend(CeedElemRestrictionApply(rstr_points, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 571 } 572 573 // Input basis apply if needed 574 CeedCallBackend(CeedOperatorInputBasisAtPoints_Hip(num_elem, num_points, qf_input_fields, op_input_fields, num_input_fields, false, e_data, impl)); 575 576 // Output pointers, as necessary 577 for (CeedInt i = 0; i < num_output_fields; i++) { 578 CeedEvalMode eval_mode; 579 580 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 581 if (eval_mode == CEED_EVAL_NONE) { 582 // Set the output Q-Vector to use the E-Vector data directly. 583 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 584 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 585 } 586 } 587 588 // Q function 589 CeedCallBackend(CeedQFunctionApply(qf, num_elem * max_num_points, impl->q_vecs_in, impl->q_vecs_out)); 590 591 // Output basis apply if needed 592 for (CeedInt i = 0; i < num_output_fields; i++) { 593 CeedEvalMode eval_mode; 594 CeedElemRestriction elem_rstr; 595 CeedBasis basis; 596 597 // Get elem_size, eval_mode, size 598 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 599 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 600 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 601 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 602 // Basis action 603 switch (eval_mode) { 604 case CEED_EVAL_NONE: 605 break; // No action 606 case CEED_EVAL_INTERP: 607 case CEED_EVAL_GRAD: 608 case CEED_EVAL_DIV: 609 case CEED_EVAL_CURL: 610 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 611 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], 612 impl->e_vecs[i + impl->num_inputs])); 613 break; 614 // LCOV_EXCL_START 615 case CEED_EVAL_WEIGHT: { 616 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 617 // LCOV_EXCL_STOP 618 } 619 } 620 } 621 622 // Output restriction 623 for (CeedInt i = 0; i < num_output_fields; i++) { 624 CeedEvalMode eval_mode; 625 CeedVector vec; 626 CeedElemRestriction elem_rstr; 627 628 // Restore evec 629 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 630 if (eval_mode == CEED_EVAL_NONE) { 631 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 632 } 633 // Get output vector 634 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 635 // Restrict 636 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 637 // Active 638 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 639 640 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], vec, request)); 641 } 642 643 // Restore input arrays 644 CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, false, e_data, impl)); 645 return CEED_ERROR_SUCCESS; 646 } 647 648 //------------------------------------------------------------------------------ 649 // Linear QFunction Assembly Core 650 //------------------------------------------------------------------------------ 651 static inline int CeedOperatorLinearAssembleQFunctionCore_Hip(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 652 CeedRequest *request) { 653 Ceed ceed, ceed_parent; 654 CeedInt num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size; 655 CeedScalar *assembled_array, *e_data[2 * CEED_FIELD_MAX] = {NULL}; 656 CeedVector *active_inputs; 657 CeedQFunctionField *qf_input_fields, *qf_output_fields; 658 CeedQFunction qf; 659 CeedOperatorField *op_input_fields, *op_output_fields; 660 CeedOperator_Hip *impl; 661 662 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 663 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 664 CeedCallBackend(CeedOperatorGetData(op, &impl)); 665 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 666 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 667 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 668 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 669 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 670 active_inputs = impl->qf_active_in; 671 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 672 673 // Setup 674 CeedCallBackend(CeedOperatorSetup_Hip(op)); 675 676 // Input Evecs and Restriction 677 CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 678 679 // Count number of active input fields 680 if (!num_active_in) { 681 for (CeedInt i = 0; i < num_input_fields; i++) { 682 CeedScalar *q_vec_array; 683 CeedVector vec; 684 685 // Get input vector 686 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 687 // Check if active input 688 if (vec == CEED_VECTOR_ACTIVE) { 689 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 690 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 691 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, &q_vec_array)); 692 CeedCallBackend(CeedRealloc(num_active_in + size, &active_inputs)); 693 for (CeedInt field = 0; field < size; field++) { 694 CeedSize q_size = (CeedSize)Q * num_elem; 695 696 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_inputs[num_active_in + field])); 697 CeedCallBackend( 698 CeedVectorSetArray(active_inputs[num_active_in + field], CEED_MEM_DEVICE, CEED_USE_POINTER, &q_vec_array[field * Q * num_elem])); 699 } 700 num_active_in += size; 701 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 702 } 703 } 704 impl->num_active_in = num_active_in; 705 impl->qf_active_in = active_inputs; 706 } 707 708 // Count number of active output fields 709 if (!num_active_out) { 710 for (CeedInt i = 0; i < num_output_fields; i++) { 711 CeedVector vec; 712 713 // Get output vector 714 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 715 // Check if active output 716 if (vec == CEED_VECTOR_ACTIVE) { 717 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 718 num_active_out += size; 719 } 720 } 721 impl->num_active_out = num_active_out; 722 } 723 724 // Check sizes 725 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 726 727 // Build objects if needed 728 if (build_objects) { 729 CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 730 CeedInt strides[3] = {1, num_elem * Q, Q}; /* *NOPAD* */ 731 732 // Create output restriction 733 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out, 734 (CeedSize)num_active_in * (CeedSize)num_active_out * (CeedSize)num_elem * (CeedSize)Q, strides, 735 rstr)); 736 // Create assembled vector 737 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 738 } 739 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 740 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_DEVICE, &assembled_array)); 741 742 // Input basis apply 743 CeedCallBackend(CeedOperatorInputBasis_Hip(num_elem, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 744 745 // Assemble QFunction 746 for (CeedInt in = 0; in < num_active_in; in++) { 747 // Set Inputs 748 CeedCallBackend(CeedVectorSetValue(active_inputs[in], 1.0)); 749 if (num_active_in > 1) { 750 CeedCallBackend(CeedVectorSetValue(active_inputs[(in + num_active_in - 1) % num_active_in], 0.0)); 751 } 752 // Set Outputs 753 for (CeedInt out = 0; out < num_output_fields; out++) { 754 CeedVector vec; 755 756 // Get output vector 757 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 758 // Check if active output 759 if (vec == CEED_VECTOR_ACTIVE) { 760 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, CEED_USE_POINTER, assembled_array)); 761 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 762 assembled_array += size * Q * num_elem; // Advance the pointer by the size of the output 763 } 764 } 765 // Apply QFunction 766 CeedCallBackend(CeedQFunctionApply(qf, Q * num_elem, impl->q_vecs_in, impl->q_vecs_out)); 767 } 768 769 // Un-set output q_vecs to prevent accidental overwrite of Assembled 770 for (CeedInt out = 0; out < num_output_fields; out++) { 771 CeedVector vec; 772 773 // Get output vector 774 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 775 // Check if active output 776 if (vec == CEED_VECTOR_ACTIVE) { 777 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_DEVICE, NULL)); 778 } 779 } 780 781 // Restore input arrays 782 CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 783 784 // Restore output 785 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 786 return CEED_ERROR_SUCCESS; 787 } 788 789 //------------------------------------------------------------------------------ 790 // Assemble Linear QFunction 791 //------------------------------------------------------------------------------ 792 static int CeedOperatorLinearAssembleQFunction_Hip(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 793 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, true, assembled, rstr, request); 794 } 795 796 //------------------------------------------------------------------------------ 797 // Update Assembled Linear QFunction 798 //------------------------------------------------------------------------------ 799 static int CeedOperatorLinearAssembleQFunctionUpdate_Hip(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 800 return CeedOperatorLinearAssembleQFunctionCore_Hip(op, false, &assembled, &rstr, request); 801 } 802 803 //------------------------------------------------------------------------------ 804 // Assemble Diagonal Setup 805 //------------------------------------------------------------------------------ 806 static inline int CeedOperatorAssembleDiagonalSetup_Hip(CeedOperator op) { 807 Ceed ceed; 808 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 809 CeedInt q_comp, num_nodes, num_qpts; 810 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 811 CeedBasis basis_in = NULL, basis_out = NULL; 812 CeedQFunctionField *qf_fields; 813 CeedQFunction qf; 814 CeedOperatorField *op_fields; 815 CeedOperator_Hip *impl; 816 817 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 818 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 819 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 820 821 // Determine active input basis 822 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 823 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 824 for (CeedInt i = 0; i < num_input_fields; i++) { 825 CeedVector vec; 826 827 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 828 if (vec == CEED_VECTOR_ACTIVE) { 829 CeedBasis basis; 830 CeedEvalMode eval_mode; 831 832 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 833 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, 834 "Backend does not implement operator diagonal assembly with multiple active bases"); 835 basis_in = basis; 836 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 837 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 838 if (eval_mode != CEED_EVAL_WEIGHT) { 839 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 840 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 841 for (CeedInt d = 0; d < q_comp; d++) eval_modes_in[num_eval_modes_in + d] = eval_mode; 842 num_eval_modes_in += q_comp; 843 } 844 } 845 } 846 847 // Determine active output basis 848 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 849 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 850 for (CeedInt i = 0; i < num_output_fields; i++) { 851 CeedVector vec; 852 853 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 854 if (vec == CEED_VECTOR_ACTIVE) { 855 CeedBasis basis; 856 CeedEvalMode eval_mode; 857 858 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 859 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 860 "Backend does not implement operator diagonal assembly with multiple active bases"); 861 basis_out = basis; 862 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 863 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 864 if (eval_mode != CEED_EVAL_WEIGHT) { 865 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF assembly 866 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 867 for (CeedInt d = 0; d < q_comp; d++) eval_modes_out[num_eval_modes_out + d] = eval_mode; 868 num_eval_modes_out += q_comp; 869 } 870 } 871 } 872 873 // Operator data struct 874 CeedCallBackend(CeedOperatorGetData(op, &impl)); 875 CeedCallBackend(CeedCalloc(1, &impl->diag)); 876 CeedOperatorDiag_Hip *diag = impl->diag; 877 878 // Basis matrices 879 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 880 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 881 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 882 const CeedInt interp_bytes = num_nodes * num_qpts * sizeof(CeedScalar); 883 const CeedInt eval_modes_bytes = sizeof(CeedEvalMode); 884 bool has_eval_none = false; 885 886 // CEED_EVAL_NONE 887 for (CeedInt i = 0; i < num_eval_modes_in; i++) has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 888 for (CeedInt i = 0; i < num_eval_modes_out; i++) has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 889 if (has_eval_none) { 890 CeedScalar *identity = NULL; 891 892 CeedCallBackend(CeedCalloc(num_nodes * num_qpts, &identity)); 893 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 894 CeedCallHip(ceed, hipMalloc((void **)&diag->d_identity, interp_bytes)); 895 CeedCallHip(ceed, hipMemcpy(diag->d_identity, identity, interp_bytes, hipMemcpyHostToDevice)); 896 CeedCallBackend(CeedFree(&identity)); 897 } 898 899 // CEED_EVAL_INTERP, CEED_EVAL_GRAD, CEED_EVAL_DIV, and CEED_EVAL_CURL 900 for (CeedInt in = 0; in < 2; in++) { 901 CeedFESpace fespace; 902 CeedBasis basis = in ? basis_in : basis_out; 903 904 CeedCallBackend(CeedBasisGetFESpace(basis, &fespace)); 905 switch (fespace) { 906 case CEED_FE_SPACE_H1: { 907 CeedInt q_comp_interp, q_comp_grad; 908 const CeedScalar *interp, *grad; 909 CeedScalar *d_interp, *d_grad; 910 911 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 912 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad)); 913 914 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 915 CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 916 CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice)); 917 CeedCallBackend(CeedBasisGetGrad(basis, &grad)); 918 CeedCallHip(ceed, hipMalloc((void **)&d_grad, interp_bytes * q_comp_grad)); 919 CeedCallHip(ceed, hipMemcpy(d_grad, grad, interp_bytes * q_comp_grad, hipMemcpyHostToDevice)); 920 if (in) { 921 diag->d_interp_in = d_interp; 922 diag->d_grad_in = d_grad; 923 } else { 924 diag->d_interp_out = d_interp; 925 diag->d_grad_out = d_grad; 926 } 927 } break; 928 case CEED_FE_SPACE_HDIV: { 929 CeedInt q_comp_interp, q_comp_div; 930 const CeedScalar *interp, *div; 931 CeedScalar *d_interp, *d_div; 932 933 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 934 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp_div)); 935 936 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 937 CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 938 CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice)); 939 CeedCallBackend(CeedBasisGetDiv(basis, &div)); 940 CeedCallHip(ceed, hipMalloc((void **)&d_div, interp_bytes * q_comp_div)); 941 CeedCallHip(ceed, hipMemcpy(d_div, div, interp_bytes * q_comp_div, hipMemcpyHostToDevice)); 942 if (in) { 943 diag->d_interp_in = d_interp; 944 diag->d_div_in = d_div; 945 } else { 946 diag->d_interp_out = d_interp; 947 diag->d_div_out = d_div; 948 } 949 } break; 950 case CEED_FE_SPACE_HCURL: { 951 CeedInt q_comp_interp, q_comp_curl; 952 const CeedScalar *interp, *curl; 953 CeedScalar *d_interp, *d_curl; 954 955 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp)); 956 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp_curl)); 957 958 CeedCallBackend(CeedBasisGetInterp(basis, &interp)); 959 CeedCallHip(ceed, hipMalloc((void **)&d_interp, interp_bytes * q_comp_interp)); 960 CeedCallHip(ceed, hipMemcpy(d_interp, interp, interp_bytes * q_comp_interp, hipMemcpyHostToDevice)); 961 CeedCallBackend(CeedBasisGetCurl(basis, &curl)); 962 CeedCallHip(ceed, hipMalloc((void **)&d_curl, interp_bytes * q_comp_curl)); 963 CeedCallHip(ceed, hipMemcpy(d_curl, curl, interp_bytes * q_comp_curl, hipMemcpyHostToDevice)); 964 if (in) { 965 diag->d_interp_in = d_interp; 966 diag->d_curl_in = d_curl; 967 } else { 968 diag->d_interp_out = d_interp; 969 diag->d_curl_out = d_curl; 970 } 971 } break; 972 } 973 } 974 975 // Arrays of eval_modes 976 CeedCallHip(ceed, hipMalloc((void **)&diag->d_eval_modes_in, num_eval_modes_in * eval_modes_bytes)); 977 CeedCallHip(ceed, hipMemcpy(diag->d_eval_modes_in, eval_modes_in, num_eval_modes_in * eval_modes_bytes, hipMemcpyHostToDevice)); 978 CeedCallHip(ceed, hipMalloc((void **)&diag->d_eval_modes_out, num_eval_modes_out * eval_modes_bytes)); 979 CeedCallHip(ceed, hipMemcpy(diag->d_eval_modes_out, eval_modes_out, num_eval_modes_out * eval_modes_bytes, hipMemcpyHostToDevice)); 980 CeedCallBackend(CeedFree(&eval_modes_in)); 981 CeedCallBackend(CeedFree(&eval_modes_out)); 982 return CEED_ERROR_SUCCESS; 983 } 984 985 //------------------------------------------------------------------------------ 986 // Assemble Diagonal Setup (Compilation) 987 //------------------------------------------------------------------------------ 988 static inline int CeedOperatorAssembleDiagonalSetupCompile_Hip(CeedOperator op, CeedInt use_ceedsize_idx, const bool is_point_block) { 989 Ceed ceed; 990 char *diagonal_kernel_source; 991 const char *diagonal_kernel_path; 992 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 993 CeedInt num_comp, q_comp, num_nodes, num_qpts; 994 CeedBasis basis_in = NULL, basis_out = NULL; 995 CeedQFunctionField *qf_fields; 996 CeedQFunction qf; 997 CeedOperatorField *op_fields; 998 CeedOperator_Hip *impl; 999 1000 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1001 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1002 CeedCallBackend(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 1003 1004 // Determine active input basis 1005 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1006 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1007 for (CeedInt i = 0; i < num_input_fields; i++) { 1008 CeedVector vec; 1009 1010 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1011 if (vec == CEED_VECTOR_ACTIVE) { 1012 CeedEvalMode eval_mode; 1013 1014 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1015 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1016 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1017 if (eval_mode != CEED_EVAL_WEIGHT) { 1018 num_eval_modes_in += q_comp; 1019 } 1020 } 1021 } 1022 1023 // Determine active output basis 1024 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1025 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1026 for (CeedInt i = 0; i < num_output_fields; i++) { 1027 CeedVector vec; 1028 1029 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1030 if (vec == CEED_VECTOR_ACTIVE) { 1031 CeedEvalMode eval_mode; 1032 1033 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1034 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1035 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1036 if (eval_mode != CEED_EVAL_WEIGHT) { 1037 num_eval_modes_out += q_comp; 1038 } 1039 } 1040 } 1041 1042 // Operator data struct 1043 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1044 CeedOperatorDiag_Hip *diag = impl->diag; 1045 1046 // Assemble kernel 1047 hipModule_t *module = is_point_block ? &diag->module_point_block : &diag->module; 1048 CeedInt elems_per_block = 1; 1049 CeedCallBackend(CeedBasisGetNumNodes(basis_in, &num_nodes)); 1050 CeedCallBackend(CeedBasisGetNumComponents(basis_in, &num_comp)); 1051 if (basis_in == CEED_BASIS_NONE) num_qpts = num_nodes; 1052 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 1053 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-ref-operator-assemble-diagonal.h", &diagonal_kernel_path)); 1054 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Kernel Source -----\n"); 1055 CeedCallBackend(CeedLoadSourceToBuffer(ceed, diagonal_kernel_path, &diagonal_kernel_source)); 1056 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Diagonal Assembly Source Complete! -----\n"); 1057 CeedCallHip(ceed, CeedCompile_Hip(ceed, diagonal_kernel_source, module, 8, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1058 num_eval_modes_out, "NUM_COMP", num_comp, "NUM_NODES", num_nodes, "NUM_QPTS", num_qpts, "USE_CEEDSIZE", 1059 use_ceedsize_idx, "USE_POINT_BLOCK", is_point_block ? 1 : 0, "BLOCK_SIZE", num_nodes * elems_per_block)); 1060 CeedCallHip(ceed, CeedGetKernel_Hip(ceed, *module, "LinearDiagonal", is_point_block ? &diag->LinearPointBlock : &diag->LinearDiagonal)); 1061 CeedCallBackend(CeedFree(&diagonal_kernel_path)); 1062 CeedCallBackend(CeedFree(&diagonal_kernel_source)); 1063 return CEED_ERROR_SUCCESS; 1064 } 1065 1066 //------------------------------------------------------------------------------ 1067 // Assemble Diagonal Core 1068 //------------------------------------------------------------------------------ 1069 static inline int CeedOperatorAssembleDiagonalCore_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request, const bool is_point_block) { 1070 Ceed ceed; 1071 CeedInt num_elem, num_nodes; 1072 CeedScalar *elem_diag_array; 1073 const CeedScalar *assembled_qf_array; 1074 CeedVector assembled_qf = NULL, elem_diag; 1075 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out, diag_rstr; 1076 CeedOperator_Hip *impl; 1077 1078 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1079 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1080 1081 // Assemble QFunction 1082 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, request)); 1083 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1084 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1085 1086 // Setup 1087 if (!impl->diag) CeedCallBackend(CeedOperatorAssembleDiagonalSetup_Hip(op)); 1088 CeedOperatorDiag_Hip *diag = impl->diag; 1089 1090 assert(diag != NULL); 1091 1092 // Assemble kernel if needed 1093 if ((!is_point_block && !diag->LinearDiagonal) || (is_point_block && !diag->LinearPointBlock)) { 1094 CeedSize assembled_length, assembled_qf_length; 1095 CeedInt use_ceedsize_idx = 0; 1096 CeedCallBackend(CeedVectorGetLength(assembled, &assembled_length)); 1097 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1098 if ((assembled_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1099 1100 CeedCallBackend(CeedOperatorAssembleDiagonalSetupCompile_Hip(op, use_ceedsize_idx, is_point_block)); 1101 } 1102 1103 // Restriction and diagonal vector 1104 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1105 CeedCheck(rstr_in == rstr_out, ceed, CEED_ERROR_BACKEND, 1106 "Cannot assemble operator diagonal with different input and output active element restrictions"); 1107 if (!is_point_block && !diag->diag_rstr) { 1108 CeedCallBackend(CeedElemRestrictionCreateUnsignedCopy(rstr_out, &diag->diag_rstr)); 1109 CeedCallBackend(CeedElemRestrictionCreateVector(diag->diag_rstr, NULL, &diag->elem_diag)); 1110 } else if (is_point_block && !diag->point_block_diag_rstr) { 1111 CeedCallBackend(CeedOperatorCreateActivePointBlockRestriction(rstr_out, &diag->point_block_diag_rstr)); 1112 CeedCallBackend(CeedElemRestrictionCreateVector(diag->point_block_diag_rstr, NULL, &diag->point_block_elem_diag)); 1113 } 1114 diag_rstr = is_point_block ? diag->point_block_diag_rstr : diag->diag_rstr; 1115 elem_diag = is_point_block ? diag->point_block_elem_diag : diag->elem_diag; 1116 CeedCallBackend(CeedVectorSetValue(elem_diag, 0.0)); 1117 1118 // Only assemble diagonal if the basis has nodes, otherwise inputs are null pointers 1119 CeedCallBackend(CeedElemRestrictionGetElementSize(diag_rstr, &num_nodes)); 1120 if (num_nodes > 0) { 1121 // Assemble element operator diagonals 1122 CeedCallBackend(CeedVectorGetArray(elem_diag, CEED_MEM_DEVICE, &elem_diag_array)); 1123 CeedCallBackend(CeedElemRestrictionGetNumElements(diag_rstr, &num_elem)); 1124 1125 // Compute the diagonal of B^T D B 1126 CeedInt elems_per_block = 1; 1127 CeedInt grid = CeedDivUpInt(num_elem, elems_per_block); 1128 void *args[] = {(void *)&num_elem, &diag->d_identity, &diag->d_interp_in, &diag->d_grad_in, &diag->d_div_in, 1129 &diag->d_curl_in, &diag->d_interp_out, &diag->d_grad_out, &diag->d_div_out, &diag->d_curl_out, 1130 &diag->d_eval_modes_in, &diag->d_eval_modes_out, &assembled_qf_array, &elem_diag_array}; 1131 1132 if (is_point_block) { 1133 CeedCallBackend(CeedRunKernelDim_Hip(ceed, diag->LinearPointBlock, grid, num_nodes, 1, elems_per_block, args)); 1134 } else { 1135 CeedCallBackend(CeedRunKernelDim_Hip(ceed, diag->LinearDiagonal, grid, num_nodes, 1, elems_per_block, args)); 1136 } 1137 1138 // Restore arrays 1139 CeedCallBackend(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 1140 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1141 } 1142 1143 // Assemble local operator diagonal 1144 CeedCallBackend(CeedElemRestrictionApply(diag_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 1145 1146 // Cleanup 1147 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1148 return CEED_ERROR_SUCCESS; 1149 } 1150 1151 //------------------------------------------------------------------------------ 1152 // Assemble Linear Diagonal 1153 //------------------------------------------------------------------------------ 1154 static int CeedOperatorLinearAssembleAddDiagonal_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1155 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, false)); 1156 return CEED_ERROR_SUCCESS; 1157 } 1158 1159 //------------------------------------------------------------------------------ 1160 // Assemble Linear Point Block Diagonal 1161 //------------------------------------------------------------------------------ 1162 static int CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1163 CeedCallBackend(CeedOperatorAssembleDiagonalCore_Hip(op, assembled, request, true)); 1164 return CEED_ERROR_SUCCESS; 1165 } 1166 1167 //------------------------------------------------------------------------------ 1168 // Single Operator Assembly Setup 1169 //------------------------------------------------------------------------------ 1170 static int CeedSingleOperatorAssembleSetup_Hip(CeedOperator op, CeedInt use_ceedsize_idx) { 1171 Ceed ceed; 1172 char *assembly_kernel_source; 1173 const char *assembly_kernel_path; 1174 CeedInt num_input_fields, num_output_fields, num_eval_modes_in = 0, num_eval_modes_out = 0; 1175 CeedInt elem_size_in, num_qpts_in = 0, num_comp_in, elem_size_out, num_qpts_out, num_comp_out, q_comp; 1176 CeedEvalMode *eval_modes_in = NULL, *eval_modes_out = NULL; 1177 CeedElemRestriction rstr_in = NULL, rstr_out = NULL; 1178 CeedBasis basis_in = NULL, basis_out = NULL; 1179 CeedQFunctionField *qf_fields; 1180 CeedQFunction qf; 1181 CeedOperatorField *input_fields, *output_fields; 1182 CeedOperator_Hip *impl; 1183 1184 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1185 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1186 1187 // Get intput and output fields 1188 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1189 1190 // Determine active input basis eval mode 1191 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1192 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 1193 for (CeedInt i = 0; i < num_input_fields; i++) { 1194 CeedVector vec; 1195 1196 CeedCallBackend(CeedOperatorFieldGetVector(input_fields[i], &vec)); 1197 if (vec == CEED_VECTOR_ACTIVE) { 1198 CeedBasis basis; 1199 CeedEvalMode eval_mode; 1200 1201 CeedCallBackend(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 1202 CeedCheck(!basis_in || basis_in == basis, ceed, CEED_ERROR_BACKEND, "Backend does not implement operator assembly with multiple active bases"); 1203 basis_in = basis; 1204 CeedCallBackend(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_in)); 1205 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1206 if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; 1207 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); 1208 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1209 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1210 if (eval_mode != CEED_EVAL_WEIGHT) { 1211 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1212 CeedCallBackend(CeedRealloc(num_eval_modes_in + q_comp, &eval_modes_in)); 1213 for (CeedInt d = 0; d < q_comp; d++) { 1214 eval_modes_in[num_eval_modes_in + d] = eval_mode; 1215 } 1216 num_eval_modes_in += q_comp; 1217 } 1218 } 1219 } 1220 1221 // Determine active output basis; basis_out and rstr_out only used if same as input, TODO 1222 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 1223 for (CeedInt i = 0; i < num_output_fields; i++) { 1224 CeedVector vec; 1225 1226 CeedCallBackend(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1227 if (vec == CEED_VECTOR_ACTIVE) { 1228 CeedBasis basis; 1229 CeedEvalMode eval_mode; 1230 1231 CeedCallBackend(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1232 CeedCheck(!basis_out || basis_out == basis, ceed, CEED_ERROR_BACKEND, 1233 "Backend does not implement operator assembly with multiple active bases"); 1234 basis_out = basis; 1235 CeedCallBackend(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_out)); 1236 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1237 if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; 1238 else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); 1239 CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED, 1240 "Active input and output bases must have the same number of quadrature points"); 1241 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1242 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1243 if (eval_mode != CEED_EVAL_WEIGHT) { 1244 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1245 CeedCallBackend(CeedRealloc(num_eval_modes_out + q_comp, &eval_modes_out)); 1246 for (CeedInt d = 0; d < q_comp; d++) { 1247 eval_modes_out[num_eval_modes_out + d] = eval_mode; 1248 } 1249 num_eval_modes_out += q_comp; 1250 } 1251 } 1252 } 1253 CeedCheck(num_eval_modes_in > 0 && num_eval_modes_out > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs"); 1254 1255 CeedCallBackend(CeedCalloc(1, &impl->asmb)); 1256 CeedOperatorAssemble_Hip *asmb = impl->asmb; 1257 asmb->elems_per_block = 1; 1258 asmb->block_size_x = elem_size_in; 1259 asmb->block_size_y = elem_size_out; 1260 1261 bool fallback = asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block > 1024; 1262 1263 if (fallback) { 1264 // Use fallback kernel with 1D threadblock 1265 asmb->block_size_y = 1; 1266 } 1267 1268 // Compile kernels 1269 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); 1270 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); 1271 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/hip/hip-ref-operator-assemble.h", &assembly_kernel_path)); 1272 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Kernel Source -----\n"); 1273 CeedCallBackend(CeedLoadSourceToBuffer(ceed, assembly_kernel_path, &assembly_kernel_source)); 1274 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Assembly Source Complete! -----\n"); 1275 CeedCallBackend(CeedCompile_Hip(ceed, assembly_kernel_source, &asmb->module, 10, "NUM_EVAL_MODES_IN", num_eval_modes_in, "NUM_EVAL_MODES_OUT", 1276 num_eval_modes_out, "NUM_COMP_IN", num_comp_in, "NUM_COMP_OUT", num_comp_out, "NUM_NODES_IN", elem_size_in, 1277 "NUM_NODES_OUT", elem_size_out, "NUM_QPTS", num_qpts_in, "BLOCK_SIZE", 1278 asmb->block_size_x * asmb->block_size_y * asmb->elems_per_block, "BLOCK_SIZE_Y", asmb->block_size_y, "USE_CEEDSIZE", 1279 use_ceedsize_idx)); 1280 CeedCallBackend(CeedGetKernel_Hip(ceed, asmb->module, "LinearAssemble", &asmb->LinearAssemble)); 1281 CeedCallBackend(CeedFree(&assembly_kernel_path)); 1282 CeedCallBackend(CeedFree(&assembly_kernel_source)); 1283 1284 // Load into B_in, in order that they will be used in eval_modes_in 1285 { 1286 const CeedInt in_bytes = elem_size_in * num_qpts_in * num_eval_modes_in * sizeof(CeedScalar); 1287 CeedInt d_in = 0; 1288 CeedEvalMode eval_modes_in_prev = CEED_EVAL_NONE; 1289 bool has_eval_none = false; 1290 CeedScalar *identity = NULL; 1291 1292 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1293 has_eval_none = has_eval_none || (eval_modes_in[i] == CEED_EVAL_NONE); 1294 } 1295 if (has_eval_none) { 1296 CeedCallBackend(CeedCalloc(elem_size_in * num_qpts_in, &identity)); 1297 for (CeedInt i = 0; i < (elem_size_in < num_qpts_in ? elem_size_in : num_qpts_in); i++) identity[i * elem_size_in + i] = 1.0; 1298 } 1299 1300 CeedCallHip(ceed, hipMalloc((void **)&asmb->d_B_in, in_bytes)); 1301 for (CeedInt i = 0; i < num_eval_modes_in; i++) { 1302 const CeedScalar *h_B_in; 1303 1304 CeedCallBackend(CeedOperatorGetBasisPointer(basis_in, eval_modes_in[i], identity, &h_B_in)); 1305 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_in, eval_modes_in[i], &q_comp)); 1306 if (q_comp > 1) { 1307 if (i == 0 || eval_modes_in[i] != eval_modes_in_prev) d_in = 0; 1308 else h_B_in = &h_B_in[(++d_in) * elem_size_in * num_qpts_in]; 1309 } 1310 eval_modes_in_prev = eval_modes_in[i]; 1311 1312 CeedCallHip(ceed, hipMemcpy(&asmb->d_B_in[i * elem_size_in * num_qpts_in], h_B_in, elem_size_in * num_qpts_in * sizeof(CeedScalar), 1313 hipMemcpyHostToDevice)); 1314 } 1315 1316 if (identity) { 1317 CeedCallBackend(CeedFree(&identity)); 1318 } 1319 } 1320 1321 // Load into B_out, in order that they will be used in eval_modes_out 1322 { 1323 const CeedInt out_bytes = elem_size_out * num_qpts_out * num_eval_modes_out * sizeof(CeedScalar); 1324 CeedInt d_out = 0; 1325 CeedEvalMode eval_modes_out_prev = CEED_EVAL_NONE; 1326 bool has_eval_none = false; 1327 CeedScalar *identity = NULL; 1328 1329 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1330 has_eval_none = has_eval_none || (eval_modes_out[i] == CEED_EVAL_NONE); 1331 } 1332 if (has_eval_none) { 1333 CeedCallBackend(CeedCalloc(elem_size_out * num_qpts_out, &identity)); 1334 for (CeedInt i = 0; i < (elem_size_out < num_qpts_out ? elem_size_out : num_qpts_out); i++) identity[i * elem_size_out + i] = 1.0; 1335 } 1336 1337 CeedCallHip(ceed, hipMalloc((void **)&asmb->d_B_out, out_bytes)); 1338 for (CeedInt i = 0; i < num_eval_modes_out; i++) { 1339 const CeedScalar *h_B_out; 1340 1341 CeedCallBackend(CeedOperatorGetBasisPointer(basis_out, eval_modes_out[i], identity, &h_B_out)); 1342 CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis_out, eval_modes_out[i], &q_comp)); 1343 if (q_comp > 1) { 1344 if (i == 0 || eval_modes_out[i] != eval_modes_out_prev) d_out = 0; 1345 else h_B_out = &h_B_out[(++d_out) * elem_size_out * num_qpts_out]; 1346 } 1347 eval_modes_out_prev = eval_modes_out[i]; 1348 1349 CeedCallHip(ceed, hipMemcpy(&asmb->d_B_out[i * elem_size_out * num_qpts_out], h_B_out, elem_size_out * num_qpts_out * sizeof(CeedScalar), 1350 hipMemcpyHostToDevice)); 1351 } 1352 1353 if (identity) { 1354 CeedCallBackend(CeedFree(&identity)); 1355 } 1356 } 1357 return CEED_ERROR_SUCCESS; 1358 } 1359 1360 //------------------------------------------------------------------------------ 1361 // Assemble matrix data for COO matrix of assembled operator. 1362 // The sparsity pattern is set by CeedOperatorLinearAssembleSymbolic. 1363 // 1364 // Note that this (and other assembly routines) currently assume only one active input restriction/basis per operator (could have multiple basis eval 1365 // modes). 1366 // TODO: allow multiple active input restrictions/basis objects 1367 //------------------------------------------------------------------------------ 1368 static int CeedSingleOperatorAssemble_Hip(CeedOperator op, CeedInt offset, CeedVector values) { 1369 Ceed ceed; 1370 CeedSize values_length = 0, assembled_qf_length = 0; 1371 CeedInt use_ceedsize_idx = 0, num_elem_in, num_elem_out, elem_size_in, elem_size_out; 1372 CeedScalar *values_array; 1373 const CeedScalar *assembled_qf_array; 1374 CeedVector assembled_qf = NULL; 1375 CeedElemRestriction assembled_rstr = NULL, rstr_in, rstr_out; 1376 CeedRestrictionType rstr_type_in, rstr_type_out; 1377 const bool *orients_in = NULL, *orients_out = NULL; 1378 const CeedInt8 *curl_orients_in = NULL, *curl_orients_out = NULL; 1379 CeedOperator_Hip *impl; 1380 1381 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1382 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1383 1384 // Assemble QFunction 1385 CeedCallBackend(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_rstr, CEED_REQUEST_IMMEDIATE)); 1386 CeedCallBackend(CeedElemRestrictionDestroy(&assembled_rstr)); 1387 CeedCallBackend(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_DEVICE, &assembled_qf_array)); 1388 1389 CeedCallBackend(CeedVectorGetLength(values, &values_length)); 1390 CeedCallBackend(CeedVectorGetLength(assembled_qf, &assembled_qf_length)); 1391 if ((values_length > INT_MAX) || (assembled_qf_length > INT_MAX)) use_ceedsize_idx = 1; 1392 1393 // Setup 1394 if (!impl->asmb) CeedCallBackend(CeedSingleOperatorAssembleSetup_Hip(op, use_ceedsize_idx)); 1395 CeedOperatorAssemble_Hip *asmb = impl->asmb; 1396 1397 assert(asmb != NULL); 1398 1399 // Assemble element operator 1400 CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_DEVICE, &values_array)); 1401 values_array += offset; 1402 1403 CeedCallBackend(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 1404 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); 1405 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 1406 1407 CeedCallBackend(CeedElemRestrictionGetType(rstr_in, &rstr_type_in)); 1408 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1409 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_in, CEED_MEM_DEVICE, &orients_in)); 1410 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1411 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_in, CEED_MEM_DEVICE, &curl_orients_in)); 1412 } 1413 1414 if (rstr_in != rstr_out) { 1415 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); 1416 CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, 1417 "Active input and output operator restrictions must have the same number of elements"); 1418 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 1419 1420 CeedCallBackend(CeedElemRestrictionGetType(rstr_out, &rstr_type_out)); 1421 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1422 CeedCallBackend(CeedElemRestrictionGetOrientations(rstr_out, CEED_MEM_DEVICE, &orients_out)); 1423 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1424 CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr_out, CEED_MEM_DEVICE, &curl_orients_out)); 1425 } 1426 } else { 1427 elem_size_out = elem_size_in; 1428 orients_out = orients_in; 1429 curl_orients_out = curl_orients_in; 1430 } 1431 1432 // Compute B^T D B 1433 CeedInt shared_mem = 1434 ((curl_orients_in || curl_orients_out ? elem_size_in * elem_size_out : 0) + (curl_orients_in ? elem_size_in * asmb->block_size_y : 0)) * 1435 sizeof(CeedScalar); 1436 CeedInt grid = CeedDivUpInt(num_elem_in, asmb->elems_per_block); 1437 void *args[] = {(void *)&num_elem_in, &asmb->d_B_in, &asmb->d_B_out, &orients_in, &curl_orients_in, 1438 &orients_out, &curl_orients_out, &assembled_qf_array, &values_array}; 1439 1440 CeedCallBackend( 1441 CeedRunKernelDimShared_Hip(ceed, asmb->LinearAssemble, grid, asmb->block_size_x, asmb->block_size_y, asmb->elems_per_block, shared_mem, args)); 1442 1443 // Restore arrays 1444 CeedCallBackend(CeedVectorRestoreArray(values, &values_array)); 1445 CeedCallBackend(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 1446 1447 // Cleanup 1448 CeedCallBackend(CeedVectorDestroy(&assembled_qf)); 1449 if (rstr_type_in == CEED_RESTRICTION_ORIENTED) { 1450 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_in, &orients_in)); 1451 } else if (rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 1452 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_in, &curl_orients_in)); 1453 } 1454 if (rstr_in != rstr_out) { 1455 if (rstr_type_out == CEED_RESTRICTION_ORIENTED) { 1456 CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr_out, &orients_out)); 1457 } else if (rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 1458 CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr_out, &curl_orients_out)); 1459 } 1460 } 1461 return CEED_ERROR_SUCCESS; 1462 } 1463 1464 //------------------------------------------------------------------------------ 1465 // Assemble Linear QFunction AtPoints 1466 //------------------------------------------------------------------------------ 1467 static int CeedOperatorLinearAssembleQFunctionAtPoints_Hip(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1468 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Backend does not implement CeedOperatorLinearAssembleQFunction"); 1469 } 1470 1471 //------------------------------------------------------------------------------ 1472 // Assemble Linear Diagonal AtPoints 1473 //------------------------------------------------------------------------------ 1474 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Hip(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1475 bool is_active_at_points = true; 1476 CeedSize e_vec_size = 0; 1477 CeedInt max_num_points, num_elem, num_input_fields, num_output_fields, elem_size_active = 1, num_comp_active = 1; 1478 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {NULL}; 1479 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1480 CeedQFunction qf; 1481 CeedOperatorField *op_input_fields, *op_output_fields; 1482 CeedOperator_Hip *impl; 1483 1484 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1485 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1486 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1487 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1488 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1489 CeedInt num_points[num_elem]; 1490 1491 // Setup 1492 CeedCallBackend(CeedOperatorSetupAtPoints_Hip(op)); 1493 max_num_points = impl->max_num_points; 1494 for (CeedInt i = 0; i < num_elem; i++) num_points[i] = max_num_points; 1495 1496 // Input Evecs and Restriction 1497 CeedCallBackend(CeedOperatorSetupInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1498 1499 // Check if active field is at points 1500 for (CeedInt i = 0; i < num_input_fields; i++) { 1501 CeedRestrictionType rstr_type; 1502 CeedVector vec; 1503 CeedElemRestriction elem_rstr; 1504 1505 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1506 // Skip non-active input 1507 if (vec != CEED_VECTOR_ACTIVE) continue; 1508 1509 // Get active restriction type 1510 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1511 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1512 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1513 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1514 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); 1515 } 1516 1517 // Get point coordinates 1518 if (!impl->point_coords_elem) { 1519 CeedVector point_coords = NULL; 1520 CeedElemRestriction rstr_points = NULL; 1521 1522 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1523 CeedCallBackend(CeedElemRestrictionCreateVector(rstr_points, NULL, &impl->point_coords_elem)); 1524 CeedCallBackend(CeedElemRestrictionApply(rstr_points, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1525 } 1526 1527 // Input basis apply if needed 1528 CeedCallBackend(CeedOperatorInputBasisAtPoints_Hip(num_elem, num_points, qf_input_fields, op_input_fields, num_input_fields, true, e_data, impl)); 1529 1530 // Output pointers, as necessary 1531 for (CeedInt i = 0; i < num_output_fields; i++) { 1532 CeedEvalMode eval_mode; 1533 1534 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1535 if (eval_mode == CEED_EVAL_NONE) { 1536 // Set the output Q-Vector to use the E-Vector data directly. 1537 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 1538 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 1539 } 1540 } 1541 1542 // Loop over active fields 1543 e_vec_size = (is_active_at_points ? max_num_points : elem_size_active) * num_comp_active; 1544 for (CeedInt s = 0; s < e_vec_size; s++) { 1545 for (CeedInt i = 0; i < num_input_fields; i++) { 1546 bool is_active_input = false; 1547 CeedEvalMode eval_mode; 1548 CeedVector vec; 1549 CeedBasis basis; 1550 1551 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1552 // Skip non-active input 1553 is_active_input = vec == CEED_VECTOR_ACTIVE; 1554 if (!is_active_input) continue; 1555 1556 // Update unit vector 1557 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs[i], 0.0)); 1558 else CeedCallBackend(CeedVectorSetValueStrided(impl->e_vecs[i], s - 1, e_vec_size, 0.0)); 1559 CeedCallBackend(CeedVectorSetValueStrided(impl->e_vecs[i], s, e_vec_size, 1.0)); 1560 1561 // Basis action 1562 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1563 switch (eval_mode) { 1564 case CEED_EVAL_NONE: 1565 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i])); 1566 break; 1567 case CEED_EVAL_INTERP: 1568 case CEED_EVAL_GRAD: 1569 case CEED_EVAL_DIV: 1570 case CEED_EVAL_CURL: 1571 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1572 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs[i], 1573 impl->q_vecs_in[i])); 1574 break; 1575 case CEED_EVAL_WEIGHT: 1576 break; // No action 1577 } 1578 } 1579 1580 // Q function 1581 CeedCallBackend(CeedQFunctionApply(qf, num_elem * max_num_points, impl->q_vecs_in, impl->q_vecs_out)); 1582 1583 // Output basis apply if needed 1584 for (CeedInt i = 0; i < num_output_fields; i++) { 1585 bool is_active_output = false; 1586 CeedEvalMode eval_mode; 1587 CeedVector vec; 1588 CeedElemRestriction elem_rstr; 1589 CeedBasis basis; 1590 1591 // Get output vector 1592 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1593 is_active_output = vec == CEED_VECTOR_ACTIVE; 1594 if (!is_active_output) continue; 1595 1596 // Basis action 1597 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1598 switch (eval_mode) { 1599 case CEED_EVAL_NONE: 1600 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 1601 break; 1602 case CEED_EVAL_INTERP: 1603 case CEED_EVAL_GRAD: 1604 case CEED_EVAL_DIV: 1605 case CEED_EVAL_CURL: 1606 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1607 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_elem, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], 1608 impl->e_vecs[i + impl->num_inputs])); 1609 break; 1610 // LCOV_EXCL_START 1611 case CEED_EVAL_WEIGHT: { 1612 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1613 // LCOV_EXCL_STOP 1614 } 1615 } 1616 1617 // Mask output e-vec 1618 { 1619 CeedInt j = num_input_fields; 1620 CeedSize out_size; 1621 1622 CeedCallBackend(CeedVectorGetLength(impl->e_vecs[i + impl->num_inputs], &out_size)); 1623 for (j = 0; j < num_input_fields; j++) { 1624 bool is_active_input = false; 1625 CeedSize in_size; 1626 CeedVector vec; 1627 1628 // Skip non-active input 1629 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[j], &vec)); 1630 is_active_input = vec == CEED_VECTOR_ACTIVE; 1631 if (!is_active_input) continue; 1632 CeedCallBackend(CeedVectorGetLength(impl->e_vecs[j], &in_size)); 1633 if (in_size == out_size) break; 1634 } 1635 CeedCheck(j < num_input_fields, CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Matching input field not found"); 1636 CeedCallBackend(CeedVectorPointwiseMult(impl->e_vecs[i + impl->num_inputs], impl->e_vecs[j], impl->e_vecs[i + impl->num_inputs])); 1637 } 1638 1639 // Restrict 1640 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1641 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs[i + impl->num_inputs], assembled, request)); 1642 1643 // Reset q_vec for 1644 if (eval_mode == CEED_EVAL_NONE) { 1645 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs[i + impl->num_inputs], CEED_MEM_DEVICE, &e_data[i + num_input_fields])); 1646 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_DEVICE, CEED_USE_POINTER, e_data[i + num_input_fields])); 1647 } 1648 } 1649 } 1650 1651 // Restore CEED_EVAL_NONE 1652 for (CeedInt i = 0; i < num_output_fields; i++) { 1653 CeedEvalMode eval_mode; 1654 CeedElemRestriction elem_rstr; 1655 1656 // Get eval_mode 1657 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1658 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1659 1660 // Restore evec 1661 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1662 if (eval_mode == CEED_EVAL_NONE) { 1663 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs[i + impl->num_inputs], &e_data[i + num_input_fields])); 1664 } 1665 } 1666 1667 // Restore input arrays 1668 CeedCallBackend(CeedOperatorRestoreInputs_Hip(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1669 return CEED_ERROR_SUCCESS; 1670 } 1671 1672 //------------------------------------------------------------------------------ 1673 // Create operator 1674 //------------------------------------------------------------------------------ 1675 int CeedOperatorCreate_Hip(CeedOperator op) { 1676 Ceed ceed; 1677 CeedOperator_Hip *impl; 1678 1679 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1680 CeedCallBackend(CeedCalloc(1, &impl)); 1681 CeedCallBackend(CeedOperatorSetData(op, impl)); 1682 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Hip)); 1683 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Hip)); 1684 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonal_Hip)); 1685 CeedCallBackend( 1686 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddPointBlockDiagonal", CeedOperatorLinearAssembleAddPointBlockDiagonal_Hip)); 1687 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedSingleOperatorAssemble_Hip)); 1688 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Hip)); 1689 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Hip)); 1690 return CEED_ERROR_SUCCESS; 1691 } 1692 1693 //------------------------------------------------------------------------------ 1694 // Create operator AtPoints 1695 //------------------------------------------------------------------------------ 1696 int CeedOperatorCreateAtPoints_Hip(CeedOperator op) { 1697 Ceed ceed; 1698 CeedOperator_Hip *impl; 1699 1700 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1701 CeedCallBackend(CeedCalloc(1, &impl)); 1702 CeedCallBackend(CeedOperatorSetData(op, impl)); 1703 1704 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Hip)); 1705 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Hip)); 1706 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Hip)); 1707 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Hip)); 1708 return CEED_ERROR_SUCCESS; 1709 } 1710 1711 //------------------------------------------------------------------------------ 1712