1! ------------------------------------------------------------------ 2! Programmer(s): Daniel R. Reynolds @ SMU 3! ------------------------------------------------------------------ 4! SUNDIALS Copyright Start 5! Copyright (c) 2002-2021, Lawrence Livermore National Security 6! and Southern Methodist University. 7! All rights reserved. 8! 9! See the top-level LICENSE and NOTICE files for details. 10! 11! SPDX-License-Identifier: BSD-3-Clause 12! SUNDIALS Copyright End 13! ------------------------------------------------------------------ 14! This is an example custom NVECTOR containing complex data. 15! ------------------------------------------------------------------ 16 17module fnvector_complex_mod 18 19 use, intrinsic :: iso_c_binding 20 use fsundials_nvector_mod 21 22 implicit none 23 24 ! ---------------------------------------------------------------- 25 type, public :: FVec 26 logical :: own_data 27 integer(c_long) :: len 28 complex(c_double_complex), pointer :: data(:) 29 end type FVec 30 ! ---------------------------------------------------------------- 31 32contains 33 34 ! ---------------------------------------------------------------- 35 function FN_VNew_Complex(n) result(sunvec_y) 36 37 implicit none 38 integer(c_long), value :: n 39 type(N_Vector), pointer :: sunvec_y 40 type(N_Vector_Ops), pointer :: ops 41 type(FVec), pointer :: content 42 43 ! allocate output N_Vector structure 44 sunvec_y => FN_VNewEmpty() 45 46 ! allocate and fill content structure 47 allocate(content) 48 allocate(content%data(n)) 49 content%own_data = .true. 50 content%len = n 51 52 ! attach the content structure to the output N_Vector 53 sunvec_y%content = c_loc(content) 54 55 ! access the N_Vector ops structure, and set internal function pointers 56 call c_f_pointer(sunvec_y%ops, ops) 57 ops%nvgetvectorid = c_funloc(FN_VGetVectorID_Complex) 58 ops%nvdestroy = c_funloc(FN_VDestroy_Complex) 59 ops%nvgetlength = c_funloc(FN_VGetLength_Complex) 60 ops%nvconst = c_funloc(FN_VConst_Complex) 61 ops%nvclone = c_funloc(FN_VClone_Complex) 62 ops%nvspace = c_funloc(FN_VSpace_Complex) 63 ops%nvlinearsum = c_funloc(FN_VLinearSum_Complex) 64 ops%nvprod = c_funloc(FN_VProd_Complex) 65 ops%nvdiv = c_funloc(FN_VDiv_Complex) 66 ops%nvscale = c_funloc(FN_VScale_Complex) 67 ops%nvabs = c_funloc(FN_VAbs_Complex) 68 ops%nvinv = c_funloc(FN_VInv_Complex) 69 ops%nvaddconst = c_funloc(FN_VAddConst_Complex) 70 ops%nvmaxnorm = c_funloc(FN_VMaxNorm_Complex) 71 ops%nvwrmsnorm = c_funloc(FN_VWRMSNorm_Complex) 72 ops%nvwrmsnormmask = c_funloc(FN_VWRMSNormMask_Complex) 73 ops%nvmin = c_funloc(FN_VMin_Complex) 74 ops%nvwl2norm = c_funloc(FN_VWL2Norm_Complex) 75 ops%nvl1norm = c_funloc(FN_VL1Norm_Complex) 76 ops%nvinvtest = c_funloc(FN_VInvTest_Complex) 77 ops%nvmaxnormlocal = c_funloc(FN_VMaxNorm_Complex) 78 ops%nvminlocal = c_funloc(FN_VMin_Complex) 79 ops%nvl1normlocal = c_funloc(FN_VL1Norm_Complex) 80 ops%nvinvtestlocal = c_funloc(FN_VInvTest_Complex) 81 ops%nvwsqrsumlocal = c_funloc(FN_VWSqrSum_Complex) 82 ops%nvwsqrsummasklocal = c_funloc(FN_VWSqrSumMask_Complex) 83 84 end function FN_VNew_Complex 85 86 ! ---------------------------------------------------------------- 87 function FN_VMake_Complex(n, data) result(sunvec_y) 88 89 implicit none 90 integer(c_long), value :: n 91 type(N_Vector), pointer :: sunvec_y 92 type(N_Vector_Ops), pointer :: ops 93 type(FVec), pointer :: content 94 complex(c_double_complex), target :: data(:) 95 96 ! allocate output N_Vector structure 97 sunvec_y => FN_VNewEmpty() 98 99 ! allocate and fill content structure 100 allocate(content) 101 content%own_data = .false. 102 content%len = n 103 content%data => data 104 105 ! attach the content structure to the output N_Vector 106 sunvec_y%content = c_loc(content) 107 108 ! access the N_Vector ops structure, and set internal function pointers 109 call c_f_pointer(sunvec_y%ops, ops) 110 ops%nvgetvectorid = c_funloc(FN_VGetVectorID_Complex) 111 ops%nvdestroy = c_funloc(FN_VDestroy_Complex) 112 ops%nvgetlength = c_funloc(FN_VGetLength_Complex) 113 ops%nvconst = c_funloc(FN_VConst_Complex) 114 ops%nvclone = c_funloc(FN_VClone_Complex) 115 ops%nvspace = c_funloc(FN_VSpace_Complex) 116 ops%nvlinearsum = c_funloc(FN_VLinearSum_Complex) 117 ops%nvprod = c_funloc(FN_VProd_Complex) 118 ops%nvdiv = c_funloc(FN_VDiv_Complex) 119 ops%nvscale = c_funloc(FN_VScale_Complex) 120 ops%nvabs = c_funloc(FN_VAbs_Complex) 121 ops%nvinv = c_funloc(FN_VInv_Complex) 122 ops%nvaddconst = c_funloc(FN_VAddConst_Complex) 123 ops%nvmaxnorm = c_funloc(FN_VMaxNorm_Complex) 124 ops%nvwrmsnorm = c_funloc(FN_VWRMSNorm_Complex) 125 ops%nvwrmsnormmask = c_funloc(FN_VWRMSNormMask_Complex) 126 ops%nvmin = c_funloc(FN_VMin_Complex) 127 ops%nvwl2norm = c_funloc(FN_VWL2Norm_Complex) 128 ops%nvl1norm = c_funloc(FN_VL1Norm_Complex) 129 ops%nvinvtest = c_funloc(FN_VInvTest_Complex) 130 ops%nvmaxnormlocal = c_funloc(FN_VMaxNorm_Complex) 131 ops%nvminlocal = c_funloc(FN_VMin_Complex) 132 ops%nvl1normlocal = c_funloc(FN_VL1Norm_Complex) 133 ops%nvinvtestlocal = c_funloc(FN_VInvTest_Complex) 134 ops%nvwsqrsumlocal = c_funloc(FN_VWSqrSum_Complex) 135 ops%nvwsqrsummasklocal = c_funloc(FN_VWSqrSumMask_Complex) 136 137 end function FN_VMake_Complex 138 139 ! ---------------------------------------------------------------- 140 function FN_VGetFVec(sunvec_x) result(x) 141 142 implicit none 143 type(N_Vector) :: sunvec_x 144 type(FVec), pointer :: x 145 146 ! extract Fortran matrix structure to output 147 call c_f_pointer(sunvec_x%content, x) 148 149 return 150 151 end function FN_VGetFVec 152 153 ! ---------------------------------------------------------------- 154 integer(N_Vector_ID) function FN_VGetVectorID_Complex(sunvec_y) & 155 result(id) bind(C) 156 157 implicit none 158 type(N_Vector) :: sunvec_y 159 160 id = SUNDIALS_NVEC_CUSTOM 161 return 162 163 end function FN_VGetVectorID_Complex 164 165 ! ---------------------------------------------------------------- 166 subroutine FN_VDestroy_Complex(sunvec_y) bind(C) 167 168 implicit none 169 type(N_Vector), target :: sunvec_y 170 type(FVec), pointer :: y 171 172 ! access FVec structure 173 y => FN_VGetFVec(sunvec_y) 174 175 ! if vector owns the data, then deallocate 176 if (y%own_data) deallocate(y%data) 177 178 ! deallocate the underlying Fortran object (the content) 179 deallocate(y) 180 181 ! set N_Vector structure members to NULL and return 182 sunvec_y%content = C_NULL_PTR 183 184 ! deallocate overall N_Vector structure 185 call FN_VFreeEmpty(sunvec_y) 186 187 return 188 189 end subroutine FN_VDestroy_Complex 190 191 ! ---------------------------------------------------------------- 192 integer(c_long) function FN_VGetLength_Complex(sunvec_y) & 193 bind(C) result(length) 194 195 implicit none 196 type(N_Vector) :: sunvec_y 197 type(FVec), pointer :: y 198 199 y => FN_VGetFVec(sunvec_y) 200 length = y%len 201 return 202 203 end function FN_VGetLength_Complex 204 205 ! ---------------------------------------------------------------- 206 subroutine FN_VConst_Complex(const, sunvec_y) bind(C) 207 208 implicit none 209 type(N_Vector) :: sunvec_y 210 real(c_double), value :: const 211 type(FVec), pointer :: y 212 213 ! extract Fortran vector structure to work with 214 y => FN_VGetFVec(sunvec_y) 215 216 ! set all entries to const (whole array operation) 217 y%data = dcmplx(const, 0.d0) 218 return 219 220 end subroutine FN_VConst_Complex 221 222 ! ---------------------------------------------------------------- 223 function FN_VClone_Complex(sunvec_x) result(y_ptr) bind(C) 224 225 implicit none 226 type(N_Vector) :: sunvec_x 227 type(N_Vector), pointer :: sunvec_y 228 type(c_ptr) :: y_ptr 229 integer(c_int) :: retval 230 type(FVec), pointer :: x, y 231 232 ! extract Fortran vector structure to work with 233 x => FN_VGetFVec(sunvec_x) 234 235 ! allocate output N_Vector structure 236 sunvec_y => FN_VNewEmpty() 237 238 ! copy operations from x into y 239 retval = FN_VCopyOps(sunvec_x, sunvec_y) 240 241 ! allocate and clone content structure 242 allocate(y) 243 allocate(y%data(x%len)) 244 y%own_data = .true. 245 y%len = x%len 246 247 ! attach the content structure to the output N_Vector 248 sunvec_y%content = c_loc(y) 249 250 ! set the c_ptr output 251 y_ptr = c_loc(sunvec_y) 252 return 253 254 end function FN_VClone_Complex 255 256 ! ---------------------------------------------------------------- 257 subroutine FN_VSpace_Complex(sunvec_x, lrw, liw) bind(C) 258 259 implicit none 260 type(N_Vector) :: sunvec_x 261 integer(c_int64_t) :: lrw(1) 262 integer(c_int64_t) :: liw(1) 263 type(FVec), pointer :: x 264 265 ! extract Fortran vector structure to work with 266 x => FN_VGetFVec(sunvec_x) 267 268 ! set output arguments and return (multiply lrw by 2 since complex) 269 lrw(1) = 2*x%len 270 liw(1) = 3 271 return 272 273 end subroutine FN_VSpace_Complex 274 275 ! ---------------------------------------------------------------- 276 subroutine FN_VLinearSum_Complex(a, sunvec_x, b, sunvec_y, sunvec_z) & 277 bind(C) 278 279 implicit none 280 type(N_Vector) :: sunvec_x 281 type(N_Vector) :: sunvec_y 282 type(N_Vector) :: sunvec_z 283 real(c_double), value :: a 284 real(c_double), value :: b 285 type(FVec), pointer :: x, y, z 286 287 ! extract Fortran vector structures to work with 288 x => FN_VGetFVec(sunvec_x) 289 y => FN_VGetFVec(sunvec_y) 290 z => FN_VGetFVec(sunvec_z) 291 292 ! perform computation (via whole array ops) and return 293 z%data = a*x%data + b*y%data 294 return 295 296 end subroutine FN_VLinearSum_Complex 297 298 ! ---------------------------------------------------------------- 299 subroutine FN_VProd_Complex(sunvec_x, sunvec_y, sunvec_z) bind(C) 300 301 implicit none 302 type(N_Vector) :: sunvec_x 303 type(N_Vector) :: sunvec_y 304 type(N_Vector) :: sunvec_z 305 type(FVec), pointer :: x, y, z 306 307 ! extract Fortran vector structures to work with 308 x => FN_VGetFVec(sunvec_x) 309 y => FN_VGetFVec(sunvec_y) 310 z => FN_VGetFVec(sunvec_z) 311 312 ! perform computation (via whole array ops) and return 313 z%data = x%data * y%data 314 return 315 316 end subroutine FN_VProd_Complex 317 318 ! ---------------------------------------------------------------- 319 subroutine FN_VDiv_Complex(sunvec_x, sunvec_y, sunvec_z) bind(C) 320 321 implicit none 322 type(N_Vector) :: sunvec_x 323 type(N_Vector) :: sunvec_y 324 type(N_Vector) :: sunvec_z 325 type(FVec), pointer :: x, y, z 326 327 ! extract Fortran vector structures to work with 328 x => FN_VGetFVec(sunvec_x) 329 y => FN_VGetFVec(sunvec_y) 330 z => FN_VGetFVec(sunvec_z) 331 332 ! perform computation (via whole array ops) and return 333 z%data = x%data / y%data 334 return 335 336 end subroutine FN_VDiv_Complex 337 338 ! ---------------------------------------------------------------- 339 subroutine FN_VScale_Complex(c, sunvec_x, sunvec_z) bind(C) 340 341 implicit none 342 real(c_double), value :: c 343 type(N_Vector) :: sunvec_x 344 type(N_Vector) :: sunvec_z 345 type(FVec), pointer :: x, z 346 347 ! extract Fortran vector structures to work with 348 x => FN_VGetFVec(sunvec_x) 349 z => FN_VGetFVec(sunvec_z) 350 351 ! perform computation (via whole array ops) and return 352 z%data = c * x%data 353 return 354 355 end subroutine FN_VScale_Complex 356 357 ! ---------------------------------------------------------------- 358 subroutine FN_VAbs_Complex(sunvec_x, sunvec_z) bind(C) 359 360 implicit none 361 type(N_Vector) :: sunvec_x 362 type(N_Vector) :: sunvec_z 363 type(FVec), pointer :: x, z 364 365 ! extract Fortran vector structures to work with 366 x => FN_VGetFVec(sunvec_x) 367 z => FN_VGetFVec(sunvec_z) 368 369 ! perform computation (via whole array ops) and return 370 z%data = abs(x%data) 371 return 372 373 end subroutine FN_VAbs_Complex 374 375 ! ---------------------------------------------------------------- 376 subroutine FN_VInv_Complex(sunvec_x, sunvec_z) bind(C) 377 378 implicit none 379 type(N_Vector) :: sunvec_x 380 type(N_Vector) :: sunvec_z 381 type(FVec), pointer :: x, z 382 383 ! extract Fortran vector structures to work with 384 x => FN_VGetFVec(sunvec_x) 385 z => FN_VGetFVec(sunvec_z) 386 387 ! perform computation (via whole array ops) and return 388 z%data = 1.d0 / x%data 389 return 390 391 end subroutine FN_VInv_Complex 392 393 ! ---------------------------------------------------------------- 394 subroutine FN_VAddConst_Complex(sunvec_x, b, sunvec_z) bind(C) 395 396 implicit none 397 type(N_Vector) :: sunvec_x 398 real(c_double), value :: b 399 type(N_Vector) :: sunvec_z 400 type(FVec), pointer :: x, z 401 402 ! extract Fortran vector structures to work with 403 x => FN_VGetFVec(sunvec_x) 404 z => FN_VGetFVec(sunvec_z) 405 406 ! perform computation (via whole array ops) and return 407 z%data = x%data + dcmplx(b, 0.d0) 408 return 409 410 end subroutine FN_VAddConst_Complex 411 412 ! ---------------------------------------------------------------- 413 real(c_double) function FN_VMaxNorm_Complex(sunvec_x) & 414 result(maxnorm) bind(C) 415 416 implicit none 417 type(N_Vector) :: sunvec_x 418 type(FVec), pointer :: x 419 420 ! extract Fortran vector structure to work with 421 x => FN_VGetFVec(sunvec_x) 422 423 ! perform computation (via whole array ops) and return 424 maxnorm = maxval(abs(x%data)) 425 return 426 427 end function FN_VMaxNorm_Complex 428 429 ! ---------------------------------------------------------------- 430 real(c_double) function FN_VWSqrSum_Complex(sunvec_x, sunvec_w) & 431 result(sqrsum) bind(C) 432 433 implicit none 434 type(N_Vector) :: sunvec_x 435 type(N_Vector) :: sunvec_w 436 type(FVec), pointer :: x, w 437 438 ! extract Fortran vector structures to work with 439 x => FN_VGetFVec(sunvec_x) 440 w => FN_VGetFVec(sunvec_w) 441 442 ! perform computation (via whole array ops) and return 443 sqrsum = sum(abs(x%data)**2 * abs(w%data)**2) 444 return 445 446 end function FN_VWSqrSum_Complex 447 448 ! ---------------------------------------------------------------- 449 real(c_double) function FN_VWSqrSumMask_Complex(sunvec_x, sunvec_w, sunvec_id) & 450 result(sqrsum) bind(C) 451 452 implicit none 453 type(N_Vector) :: sunvec_x 454 type(N_Vector) :: sunvec_w 455 type(N_Vector) :: sunvec_id 456 type(FVec), pointer :: x, w, id 457 integer(c_long) :: i 458 459 ! extract Fortran vector structures to work with 460 x => FN_VGetFVec(sunvec_x) 461 w => FN_VGetFVec(sunvec_w) 462 id => FN_VGetFVec(sunvec_id) 463 464 ! perform computation and return 465 sqrsum = 0.d0 466 do i = 1,x%len 467 if (real(id%data(i)) > 0.d0) then 468 sqrsum = sqrsum + (abs(x%data(i)) * abs(w%data(i)))**2 469 end if 470 end do 471 return 472 473 end function FN_VWSqrSumMask_Complex 474 475 ! ---------------------------------------------------------------- 476 real(c_double) function FN_VWRMSNorm_Complex(sunvec_x, sunvec_w) & 477 result(wrmsnorm) bind(C) 478 479 implicit none 480 type(N_Vector) :: sunvec_x 481 type(N_Vector) :: sunvec_w 482 type(FVec), pointer :: x 483 484 ! extract Fortran vector structure to work with 485 x => FN_VGetFVec(sunvec_x) 486 487 ! postprocess result from FN_VWSqrSum for result 488 wrmsnorm = dsqrt( FN_VWSqrSum_Complex(sunvec_x, sunvec_w) / x%len ) 489 return 490 491 end function FN_VWRMSNorm_Complex 492 493 ! ---------------------------------------------------------------- 494 real(c_double) function FN_VWRMSNormMask_Complex(sunvec_x, sunvec_w, sunvec_id) & 495 result(wrmsnorm) bind(C) 496 497 implicit none 498 type(N_Vector) :: sunvec_x 499 type(N_Vector) :: sunvec_w 500 type(N_Vector) :: sunvec_id 501 type(FVec), pointer :: x 502 503 ! extract Fortran vector structure to work with 504 x => FN_VGetFVec(sunvec_x) 505 506 ! postprocess result from FN_VWSqrSumMask for result 507 wrmsnorm = dsqrt( FN_VWSqrSumMask_Complex(sunvec_x, sunvec_w, sunvec_id) / x%len ) 508 return 509 510 end function FN_VWRMSNormMask_Complex 511 512 ! ---------------------------------------------------------------- 513 real(c_double) function FN_VMin_Complex(sunvec_x) & 514 result(mnval) bind(C) 515 516 implicit none 517 type(N_Vector) :: sunvec_x 518 type(FVec), pointer :: x 519 520 ! extract Fortran vector structure to work with 521 x => FN_VGetFVec(sunvec_x) 522 523 ! perform computation (via whole array ops) and return 524 ! Note: SUNDIALS only wants the minimum of all real components 525 mnval = minval(real(x%data)) 526 return 527 528 end function FN_VMin_Complex 529 530 ! ---------------------------------------------------------------- 531 real(c_double) function FN_VWL2Norm_Complex(sunvec_x, sunvec_w) & 532 result(wl2norm) bind(C) 533 534 implicit none 535 type(N_Vector) :: sunvec_x 536 type(N_Vector) :: sunvec_w 537 type(FVec), pointer :: x 538 539 ! extract Fortran vector structure to work with 540 x => FN_VGetFVec(sunvec_x) 541 542 ! postprocess result from FN_VWSqrSum for result 543 wl2norm = dsqrt(FN_VWSqrSum_Complex(sunvec_x, sunvec_w)) 544 return 545 546 end function FN_VWL2Norm_Complex 547 548 ! ---------------------------------------------------------------- 549 real(c_double) function FN_VL1Norm_Complex(sunvec_x) & 550 result(l1norm) bind(C) 551 552 implicit none 553 type(N_Vector) :: sunvec_x 554 type(FVec), pointer :: x 555 556 ! extract Fortran vector structure to work with 557 x => FN_VGetFVec(sunvec_x) 558 559 ! perform computation (via whole array ops) and return 560 l1norm = sum(abs(x%data)) 561 return 562 563 end function FN_VL1Norm_Complex 564 565 ! ---------------------------------------------------------------- 566 integer(c_int) function FN_VInvTest_Complex(sunvec_x, sunvec_z) & 567 result(no_zero_found) bind(C) 568 569 implicit none 570 type(N_Vector) :: sunvec_x 571 type(N_Vector) :: sunvec_z 572 type(FVec), pointer :: x, z 573 integer(c_long) :: i 574 575 ! extract Fortran vector structures to work with 576 x => FN_VGetFVec(sunvec_x) 577 z => FN_VGetFVec(sunvec_z) 578 579 ! perform operation and return 580 no_zero_found = 1 581 do i = 1,x%len 582 if (x%data(i) == dcmplx(0.d0, 0.d0)) then 583 no_zero_found = 0 584 else 585 z%data(i) = 1.d0 / x%data(i) 586 end if 587 end do 588 return 589 590 end function FN_VInvTest_Complex 591 592end module fnvector_complex_mod 593! ------------------------------------------------------------------ 594