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