GCC Code Coverage Report


Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 96.5% 274 / 0 / 284
Functions: 77.8% 7 / 0 / 9
Branches: 67.7% 593 / 0 / 876

tensorprod/m_tensorprod_matrix.f90
Line Branch Exec Source
1 !> @brief Module for the mass matrix of a 3D tensor product B-spline space
2 !>
3 !> The mass matrix is stored in a sparse dedicated format, which can easily be converted to a PETSc matrix
4 module m_tensorprod_matrix
5 use m_common, only: wp, user_function_3d_interface, user_function_1d_interface
6 use m_tensorprod_basis, only: TensorProdSpace
7 use m_tensorprod_indices, only: TensorProdIndices
8 #include "petsc.fi"
9
10 implicit none
11
12 private
13 public :: TensorProdMat
14 public :: get_petsc, mass_matrix
15
16 !> @brief Type for the mass matrix of a 3D tensor product B-spline basis
17 type TensorProdMat
18 !> The tensor product B-spline space for the row index
19 type(TensorProdSpace) :: space_row
20 !> The tensor product B-spline space for the column index
21 type(TensorProdSpace) :: space_col
22
23 !> The values of the mass matrix (:, :, :, local_row)
24 real(wp), allocatable :: vals(:, :, :, :, :, :)
25
26 !> The column indices of the mass matrix (:, :, :, local_row)
27 integer, allocatable :: cols(:, :, :, :, :, :)
28
29 !> The row indices of the mass matrix (local_row -> global_row)
30 integer, allocatable :: rows(:, :, :)
31
32 !>
33 type(TensorProdIndices) :: rank_actv_bsplines
34
35 !> The number of non-zero entries per row
36 integer :: total_nz_per_row
37
38 !> The numer of rows
39 integer :: nr_rows
40 !> The number of rows for the current rank
41 integer :: my_nr_rows
42
43 !> The number of columns
44 integer :: nr_cols
45 !> The number of columns for the current rank
46 integer :: my_nr_cols
47
48 !> Whether the matrix is symmetric
49 logical :: symmetric
50
51 !> Whether the matrix is symmetric positive definite
52 logical :: spd
53 contains
54 procedure :: init => init_mat
55 procedure :: destroy => destroy_mat
56 end type TensorProdMat
57
58 !> @brief Construct a mass matrix for two 3D tensor product B-spline spaces
59 interface mass_matrix
60 procedure mass_matrix_3d
61 end interface mass_matrix
62
63 !> @brief Get the PETSc matrix from the (array of) TensorProdMat(s)
64 interface get_petsc
65 module procedure get_petsc_mat_single, get_petsc_mat_blocks, get_petsc_mat_diagonal
66 end interface
67 contains
68
69 !> @brief Initialize the TensorProdMat
70 !>
71 !> @param[inout] this The TensorProdMat object
72 !> @param[in] space1 The first tensor product B-spline basis
73 !> @param[in] space2 The second tensor product B-spline basis
74 !> @param[in] symmetric _(optional)_ Whether the matrix is symmetric
75 !> @param[in] spd _(optional)_ Whether the matrix is symmetric positive definite
76 1957 subroutine init_mat(this, space1, space2, symmetric, spd)
77 use m_tensorprod_basis, only: size
78 use m_tensorprod_indices, only: size
79 implicit none
80 class(TensorProdMat), intent(inout) :: this
81 type(TensorProdSpace), intent(in) :: space1, space2
82 logical, intent(in), optional :: symmetric, spd
83
84 integer :: sze1, sze2, sze3
85
86 1957 this%space_row = space1
87 1957 this%space_col = space2
88
89 1957 sze1 = 1 + space1%spaces(1)%degree + space2%spaces(1)%degree
90 1957 sze2 = 1 + space1%spaces(2)%degree + space2%spaces(2)%degree
91 1957 sze3 = 1 + space1%spaces(3)%degree + space2%spaces(3)%degree
92
93 ! The number of local rows is equal to the active B-splines on this rank
94 1957 this%rank_actv_bsplines = space1%rank_resp_intervals
95 1957 this%rank_actv_bsplines%i1 = this%rank_actv_bsplines%i1 + space1%spaces(1)%degree
96 1957 this%rank_actv_bsplines%j1 = this%rank_actv_bsplines%j1 + space1%spaces(2)%degree
97 1957 this%rank_actv_bsplines%k1 = this%rank_actv_bsplines%k1 + space1%spaces(3)%degree
98
99 allocate (this%vals(0:sze1 - 1, 0:sze2 - 1, 0:sze3 - 1, &
100 this%rank_actv_bsplines%i0:this%rank_actv_bsplines%i1, &
101 this%rank_actv_bsplines%j0:this%rank_actv_bsplines%j1, &
102
22/34
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 1957 times.
✓ Branch 4 → 3 taken 1957 times.
✗ Branch 4 → 5 not taken.
✓ Branch 5 → 6 taken 47 times.
✓ Branch 5 → 7 taken 1910 times.
✓ Branch 7 → 6 taken 1910 times.
✗ Branch 7 → 8 not taken.
✓ Branch 8 → 9 taken 336 times.
✓ Branch 8 → 10 taken 1621 times.
✓ Branch 10 → 9 taken 1621 times.
✗ Branch 10 → 11 not taken.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 13 taken 1957 times.
✓ Branch 13 → 12 taken 1957 times.
✗ Branch 13 → 14 not taken.
✓ Branch 14 → 15 taken 47 times.
✓ Branch 14 → 16 taken 1910 times.
✓ Branch 16 → 15 taken 1910 times.
✗ Branch 16 → 17 not taken.
✓ Branch 17 → 18 taken 336 times.
✓ Branch 17 → 19 taken 1621 times.
✓ Branch 19 → 18 taken 1621 times.
✗ Branch 19 → 20 not taken.
✓ Branch 20 → 21 taken 1957 times.
✗ Branch 20 → 22 not taken.
✓ Branch 22 → 23 taken 1574 times.
✓ Branch 22 → 24 taken 383 times.
✗ Branch 24 → 25 not taken.
✓ Branch 24 → 26 taken 1957 times.
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 1957 times.
✗ Branch 28 → 29 not taken.
✓ Branch 28 → 30 taken 1957 times.
15656 this%rank_actv_bsplines%k0:this%rank_actv_bsplines%k1))
103 allocate (this%cols(0:sze1 - 1, 0:sze2 - 1, 0:sze3 - 1, &
104 this%rank_actv_bsplines%i0:this%rank_actv_bsplines%i1, &
105 this%rank_actv_bsplines%j0:this%rank_actv_bsplines%j1, &
106
22/34
✗ Branch 30 → 31 not taken.
✓ Branch 30 → 32 taken 1957 times.
✓ Branch 32 → 31 taken 1957 times.
✗ Branch 32 → 33 not taken.
✓ Branch 33 → 34 taken 47 times.
✓ Branch 33 → 35 taken 1910 times.
✓ Branch 35 → 34 taken 1910 times.
✗ Branch 35 → 36 not taken.
✓ Branch 36 → 37 taken 336 times.
✓ Branch 36 → 38 taken 1621 times.
✓ Branch 38 → 37 taken 1621 times.
✗ Branch 38 → 39 not taken.
✗ Branch 39 → 40 not taken.
✓ Branch 39 → 41 taken 1957 times.
✓ Branch 41 → 40 taken 1957 times.
✗ Branch 41 → 42 not taken.
✓ Branch 42 → 43 taken 47 times.
✓ Branch 42 → 44 taken 1910 times.
✓ Branch 44 → 43 taken 1910 times.
✗ Branch 44 → 45 not taken.
✓ Branch 45 → 46 taken 336 times.
✓ Branch 45 → 47 taken 1621 times.
✓ Branch 47 → 46 taken 1621 times.
✗ Branch 47 → 48 not taken.
✓ Branch 48 → 49 taken 1957 times.
✗ Branch 48 → 50 not taken.
✓ Branch 50 → 51 taken 1574 times.
✓ Branch 50 → 52 taken 383 times.
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 1957 times.
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 1957 times.
✗ Branch 56 → 57 not taken.
✓ Branch 56 → 58 taken 1957 times.
15656 this%rank_actv_bsplines%k0:this%rank_actv_bsplines%k1))
107 allocate (this%rows(this%rank_actv_bsplines%i0:this%rank_actv_bsplines%i1, &
108 this%rank_actv_bsplines%j0:this%rank_actv_bsplines%j1, &
109
14/22
✗ Branch 58 → 59 not taken.
✓ Branch 58 → 60 taken 1957 times.
✓ Branch 60 → 59 taken 1957 times.
✗ Branch 60 → 61 not taken.
✓ Branch 61 → 62 taken 47 times.
✓ Branch 61 → 63 taken 1910 times.
✓ Branch 63 → 62 taken 1910 times.
✗ Branch 63 → 64 not taken.
✓ Branch 64 → 65 taken 336 times.
✓ Branch 64 → 66 taken 1621 times.
✓ Branch 66 → 65 taken 1621 times.
✗ Branch 66 → 67 not taken.
✓ Branch 67 → 68 taken 1957 times.
✗ Branch 67 → 69 not taken.
✓ Branch 69 → 70 taken 1574 times.
✓ Branch 69 → 71 taken 383 times.
✗ Branch 71 → 72 not taken.
✓ Branch 71 → 73 taken 1957 times.
✗ Branch 73 → 74 not taken.
✓ Branch 73 → 75 taken 1957 times.
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 1957 times.
9785 this%rank_actv_bsplines%k0:this%rank_actv_bsplines%k1))
110
111
12/12
✓ Branch 78 → 79 taken 11838 times.
✓ Branch 78 → 95 taken 1957 times.
✓ Branch 80 → 81 taken 171765 times.
✓ Branch 80 → 94 taken 11838 times.
✓ Branch 82 → 83 taken 2740871 times.
✓ Branch 82 → 93 taken 171765 times.
✓ Branch 84 → 85 taken 13596635 times.
✓ Branch 84 → 92 taken 2740871 times.
✓ Branch 86 → 87 taken 87471679 times.
✓ Branch 86 → 91 taken 13596635 times.
✓ Branch 88 → 89 taken 612037867 times.
✓ Branch 88 → 90 taken 87471679 times.
716032612 this%vals(:, :, :, :, :, :) = 0._wp
112
12/12
✓ Branch 95 → 96 taken 11838 times.
✓ Branch 95 → 112 taken 1957 times.
✓ Branch 97 → 98 taken 171765 times.
✓ Branch 97 → 111 taken 11838 times.
✓ Branch 99 → 100 taken 2740871 times.
✓ Branch 99 → 110 taken 171765 times.
✓ Branch 101 → 102 taken 13596635 times.
✓ Branch 101 → 109 taken 2740871 times.
✓ Branch 103 → 104 taken 87471679 times.
✓ Branch 103 → 108 taken 13596635 times.
✓ Branch 105 → 106 taken 612037867 times.
✓ Branch 105 → 107 taken 87471679 times.
716032612 this%cols(:, :, :, :, :, :) = -1
113
6/6
✓ Branch 112 → 113 taken 11838 times.
✓ Branch 112 → 120 taken 1957 times.
✓ Branch 114 → 115 taken 171765 times.
✓ Branch 114 → 119 taken 11838 times.
✓ Branch 116 → 117 taken 2740871 times.
✓ Branch 116 → 118 taken 171765 times.
2926431 this%rows(:, :, :) = -1
114
115 1957 this%total_nz_per_row = sze1 * sze2 * sze3
116 1957 this%nr_rows = size(space1)
117 1957 this%nr_cols = size(space2)
118
119 1957 this%my_nr_rows = size(space1, my_rank=.true.)
120 1957 this%my_nr_cols = size(space2, my_rank=.true.)
121
122 1957 this%symmetric = .false.
123 1957 this%spd = .false.
124
125
1/2
✓ Branch 120 → 121 taken 1957 times.
✗ Branch 120 → 122 not taken.
1957 if (present(symmetric)) this%symmetric = symmetric
126
1/2
✓ Branch 122 → 123 taken 1957 times.
✗ Branch 122 → 124 not taken.
1957 if (present(spd)) this%spd = spd
127 1957 end subroutine init_mat
128
129 !> @brief destroy for the TensorProdMat
130 !>
131 !> @param[inout] this The TensorProdMat object
132 4939 subroutine destroy_mat(this)
133 implicit none
134 class(TensorProdMat), intent(inout) :: this
135
136
2/2
✓ Branch 2 → 3 taken 1957 times.
✓ Branch 2 → 4 taken 2982 times.
4939 if (allocated(this%vals)) deallocate (this%vals)
137
2/2
✓ Branch 4 → 5 taken 1957 times.
✓ Branch 4 → 6 taken 2982 times.
4939 if (allocated(this%cols)) deallocate (this%cols)
138 4939 end subroutine destroy_mat
139
140 !> @brief Get the PETSc matrix from the TensorProdMat
141 !>
142 !> @param[out] M The PETSc matrix
143 !> @param[in] this The TensorProdMat object
144 232 subroutine get_petsc_mat_single(M, this)
145 implicit none
146 Mat, intent(out) :: M
147 type(TensorProdMat), intent(in) :: this
148
149 integer :: ierr, nonlocal_nz_per_row
150 integer :: i, j, k
151
152 232 nonlocal_nz_per_row = this%total_nz_per_row ! TODO a better upper bound of the local/non-local nz per row
153 call MatCreateAIJ(this%space_row%domain%comm, this%my_nr_rows, this%my_nr_cols, this%nr_rows, &
154 & this%nr_cols, this%total_nz_per_row, PETSC_NULL_INTEGER_ARRAY, nonlocal_nz_per_row, &
155
1/2
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 232 times.
232 & PETSC_NULL_INTEGER_ARRAY, M, ierr); CHKERRQ(ierr)
156
157
1/2
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 232 times.
232 PetscCall(MatSetFromOptions(M, ierr))
158
1/2
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 232 times.
232 PetscCall(MatSetOption(M, MAT_SYMMETRIC, merge(PETSC_TRUE, PETSC_FALSE, this%symmetric), ierr))
159
1/2
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 232 times.
232 PetscCall(MatSetOption(M, MAT_SPD, merge(PETSC_TRUE, PETSC_FALSE, this%spd), ierr))
160
1/2
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 17 taken 232 times.
232 PetscCall(MatSetUp(M, ierr))
161
162 ! Copy the values to the PETSC matrix object
163
2/2
✓ Branch 18 → 19 taken 1577 times.
✓ Branch 18 → 29 taken 232 times.
1809 do k = this%rank_actv_bsplines%k0, this%rank_actv_bsplines%k1
164
2/2
✓ Branch 20 → 21 taken 16807 times.
✓ Branch 20 → 27 taken 1577 times.
18616 do j = this%rank_actv_bsplines%j0, this%rank_actv_bsplines%j1
165
2/2
✓ Branch 22 → 23 taken 290069 times.
✓ Branch 22 → 25 taken 16807 times.
308453 do i = this%rank_actv_bsplines%i0, this%rank_actv_bsplines%i1
166 call petsc_insert_single_row(M, this%rows(i, j, k), this%total_nz_per_row, this%cols(:, :, :, i, j, k), &
167 306876 this%vals(:, :, :, i, j, k))
168 end do
169 end do
170 end do
171
172
1/2
✗ Branch 31 → 32 not taken.
✓ Branch 31 → 33 taken 232 times.
232 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY, ierr))
173
1/2
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 232 times.
232 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY, ierr))
174
175 end subroutine get_petsc_mat_single
176
177 !> @brief Insert a single row into the PETSc matrix
178 !>
179 !> @param[inout] M The PETSc matrix
180 !> @param[in] row The row index to insert
181 !> @param[in] nr_elements The number of elements to insert
182 !> @param[in] cols The column indices of the elements to insert
183 !> @param[in] vals The values of the elements to insert
184 2740871 subroutine petsc_insert_single_row(M, row, nr_elements, cols, vals)
185 implicit none
186
187 Mat, intent(inout) :: M
188 integer, intent(in) :: row
189 integer, intent(in) :: nr_elements
190 integer, contiguous, target, intent(in) :: cols(:, :, :)
191 real(wp), contiguous, target, intent(in) :: vals(:, :, :)
192
193 2740871 integer, pointer :: col_pointer(:)
194 real(wp), pointer :: val_pointer(:)
195 integer :: ierr
196
197 ! Avoid allocating for reshape
198
199 2740871 col_pointer(1:nr_elements) => cols(:, :, :)
200 2740871 val_pointer(1:nr_elements) => vals(:, :, :)
201
5/8
✓ Branch 3 → 4 taken 2740871 times.
✓ Branch 3 → 5 taken 2740871 times.
✗ Branch 8 → 9 not taken.
✓ Branch 8 → 11 taken 2740871 times.
✗ Branch 11 → 12 not taken.
✓ Branch 11 → 14 taken 2740871 times.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 17 taken 2740871 times.
5481742 PetscCall(MatSetValues(M, 1, [row], nr_elements, col_pointer, val_pointer, ADD_VALUES, ierr))
202 ! NOTE: use ADD_VALUES because a cols might not be unique
203 end subroutine
204
205 !> @brief Get the PETSc matrix from the block-diagonal TensorProdMat blocks
206 !>
207 !> @param[out] M The PETSc matrix
208 !> @param[in] this The block-diagonal TensorProdMat objects
209 !>
210 !> @note It is assumed thaat the column indices of the blocks of 'this' are global indices
211
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 497 times.
497 subroutine get_petsc_mat_diagonal(M, this)
212 implicit none
213 Mat, intent(out) :: M
214 type(TensorProdMat), intent(in) :: this(:)
215
216 integer :: ierr, nonlocal_nz_per_row, local_nz_per_row
217
218 integer :: blk, i, j, k
219 integer :: nr_rows, nr_cols, my_nr_rows, my_nr_cols
220
221
2/2
✓ Branch 5 → 6 taken 1491 times.
✓ Branch 5 → 7 taken 497 times.
1988 nr_rows = sum(this%nr_rows)
222
2/2
✓ Branch 8 → 9 taken 1491 times.
✓ Branch 8 → 10 taken 497 times.
1988 nr_cols = sum(this%nr_cols)
223
2/2
✓ Branch 11 → 12 taken 1491 times.
✓ Branch 11 → 13 taken 497 times.
1988 my_nr_rows = sum(this%my_nr_rows)
224
2/2
✓ Branch 14 → 15 taken 1491 times.
✓ Branch 14 → 16 taken 497 times.
1988 my_nr_cols = sum(this%my_nr_cols)
225
226 ! TODO a better upper bound of the local/non-local nz per row
227
2/2
✓ Branch 17 → 18 taken 1491 times.
✓ Branch 17 → 19 taken 497 times.
1988 local_nz_per_row = maxval(this%total_nz_per_row)
228 497 nonlocal_nz_per_row = local_nz_per_row
229
230 call MatCreateAIJ(this(1)%space_row%domain%comm, &
231 & my_nr_rows, &
232 & my_nr_cols, &
233 & nr_rows, &
234 & nr_cols, &
235 & local_nz_per_row, &
236
1/2
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 22 taken 497 times.
497 & PETSC_NULL_INTEGER_ARRAY, nonlocal_nz_per_row, PETSC_NULL_INTEGER_ARRAY, M, ierr); CHKERRQ(ierr)
237
238
1/2
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 497 times.
497 PetscCall(MatSetFromOptions(M, ierr))
239
1/2
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 497 times.
497 PetscCall(MatSetOption(M, MAT_SYMMETRIC, PETSC_TRUE, ierr))
240
1/2
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 497 times.
497 PetscCall(MatSetOption(M, MAT_SPD, PETSC_TRUE, ierr))
241
1/2
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 497 times.
497 PetscCall(MatSetUp(M, ierr))
242
243 ! Copy the values to the PETSC matrix object
244
2/2
✓ Branch 35 → 36 taken 1491 times.
✓ Branch 35 → 50 taken 497 times.
1988 do blk = 1, size(this)
245 ! The indices are already global
246
2/2
✓ Branch 37 → 38 taken 6505 times.
✓ Branch 37 → 48 taken 1491 times.
8493 do k = this(blk)%rank_actv_bsplines%k0, this(blk)%rank_actv_bsplines%k1
247
2/2
✓ Branch 39 → 40 taken 77552 times.
✓ Branch 39 → 46 taken 6505 times.
85548 do j = this(blk)%rank_actv_bsplines%j0, this(blk)%rank_actv_bsplines%j1
248
2/2
✓ Branch 41 → 42 taken 1464486 times.
✓ Branch 41 → 44 taken 77552 times.
1548543 do i = this(blk)%rank_actv_bsplines%i0, this(blk)%rank_actv_bsplines%i1
249 call petsc_insert_single_row(M, this(blk)%rows(i, j, k), this(blk)%total_nz_per_row, this(blk)%cols(:, :, :, i, j, k), &
250 1542038 this(blk)%vals(:, :, :, i, j, k))
251 end do
252 end do
253 end do
254 end do
255
256
1/2
✗ Branch 52 → 53 not taken.
✓ Branch 52 → 54 taken 497 times.
497 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY, ierr))
257
1/2
✗ Branch 55 → 56 not taken.
✓ Branch 55 → 57 taken 497 times.
497 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY, ierr))
258
259 end subroutine get_petsc_mat_diagonal
260
261 !> @brief Get the PETSc matrix from the TensorProdMat blocks
262 !>
263 !> @param[out] M The PETSc matrix
264 !> @param[in] this The TensorProdMat objects
265 !>
266 !> @note It is assumed thaat the column indices of the blocks of 'this' are global indices
267
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 26 times.
26 subroutine get_petsc_mat_blocks(M, this)
268 implicit none
269 Mat, intent(out) :: M
270 type(TensorProdMat), intent(in) :: this(:, :)
271
272 integer :: ierr, nonlocal_nz_per_row, local_nz_per_row
273 integer :: blk_row, blk_col, i, j, k
274
275 integer :: nr_rows, nr_cols, my_nr_rows, my_nr_cols
276
277
2/2
✓ Branch 5 → 6 taken 78 times.
✓ Branch 5 → 7 taken 26 times.
104 nr_rows = sum(this(:, 1)%nr_rows)
278
2/2
✓ Branch 8 → 9 taken 78 times.
✓ Branch 8 → 10 taken 26 times.
104 nr_cols = sum(this(1, :)%nr_cols)
279
2/2
✓ Branch 11 → 12 taken 78 times.
✓ Branch 11 → 13 taken 26 times.
104 my_nr_rows = sum(this(:, 1)%my_nr_rows)
280
2/2
✓ Branch 14 → 15 taken 78 times.
✓ Branch 14 → 16 taken 26 times.
104 my_nr_cols = sum(this(1, :)%my_nr_cols)
281
282
4/4
✓ Branch 17 → 18 taken 78 times.
✓ Branch 17 → 22 taken 26 times.
✓ Branch 19 → 20 taken 234 times.
✓ Branch 19 → 21 taken 78 times.
338 local_nz_per_row = maxval(sum(this%total_nz_per_row, dim=2))
283 26 nonlocal_nz_per_row = local_nz_per_row ! TODO a better upper bound of the local/non-local nz per row
284 call MatCreateAIJ(this(1, 1)%space_row%domain%comm, &
285 & my_nr_rows, &
286 & my_nr_cols, &
287 & nr_rows, &
288 & nr_cols, &
289 & local_nz_per_row, &
290
1/2
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 26 times.
26 & PETSC_NULL_INTEGER_ARRAY, nonlocal_nz_per_row, PETSC_NULL_INTEGER_ARRAY, M, ierr); CHKERRQ(ierr)
291
292
1/2
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 26 times.
26 PetscCall(MatSetFromOptions(M, ierr))
293
1/2
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 31 taken 26 times.
26 PetscCall(MatSetOption(M, MAT_SYMMETRIC, PETSC_TRUE, ierr))
294
1/2
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 26 times.
26 PetscCall(MatSetOption(M, MAT_SPD, PETSC_TRUE, ierr))
295
1/2
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 26 times.
26 PetscCall(MatSetUp(M, ierr))
296
297 ! Copy the values to the PETSC matrix object
298
2/2
✓ Branch 38 → 39 taken 78 times.
✓ Branch 38 → 57 taken 26 times.
104 do blk_row = 1, size(this, dim=1)
299
2/2
✓ Branch 40 → 41 taken 234 times.
✓ Branch 40 → 55 taken 78 times.
338 do blk_col = 1, size(this, dim=2)
300
2/2
✓ Branch 42 → 43 taken 3756 times.
✓ Branch 42 → 53 taken 234 times.
4068 do k = this(blk_row, blk_col)%rank_actv_bsplines%k0, this(blk_row, blk_col)%rank_actv_bsplines%k1
301
2/2
✓ Branch 44 → 45 taken 77406 times.
✓ Branch 44 → 51 taken 3756 times.
81396 do j = this(blk_row, blk_col)%rank_actv_bsplines%j0, this(blk_row, blk_col)%rank_actv_bsplines%j1
302
2/2
✓ Branch 46 → 47 taken 986316 times.
✓ Branch 46 → 49 taken 77406 times.
1067478 do i = this(blk_row, blk_col)%rank_actv_bsplines%i0, this(blk_row, blk_col)%rank_actv_bsplines%i1
303 call petsc_insert_single_row(M, this(blk_row, blk_col)%rows(i, j, k), this(blk_row, blk_col)%total_nz_per_row, &
304 1063722 & this(blk_row, blk_col)%cols(:, :, :, i, j, k), this(blk_row, blk_col)%vals(:, :, :, i, j, k))
305 end do
306 end do
307 end do
308 end do
309 end do
310
311
1/2
✗ Branch 59 → 60 not taken.
✓ Branch 59 → 61 taken 26 times.
26 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY, ierr))
312
1/2
✗ Branch 62 → 63 not taken.
✓ Branch 62 → 64 taken 26 times.
26 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY, ierr))
313
314 end subroutine get_petsc_mat_blocks
315
316 !> @brief Compute the mass matrix for a 3D tensor product B-spline basis
317 !>
318 !> The resulting matrix is such that the (row, col) element correspnds to the L2 inner product of the 'row' B-spline basis
319 !> function from the first basis and the 'col' B-spline basis function from the second basis, optionally weighted with coordfun
320 !>
321 !> @param M The mass matrix of type TensorProdMat
322 !> @param tp_space1 The first tensor product B-spline basis
323 !> @param tp_space2 The second tensor product B-spline basis
324 !> @param coordfun _(optional)_ The coordinate function
325 !> @param n_quad_extra _(optional)_ The extra number of quadrature points on top of those needed to exactly integrate the B-spline
326 !> spaces
327
2/4
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 1957 times.
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 14 taken 1957 times.
1957 subroutine mass_matrix_3d(M, tp_space1, tp_space2, coordfun, n_quad_extra)
328 use m_bspline_basis, only: MAX_BSPLINE_DEGREE, get_internal_and_actual_inds
329 use m_bspline_quadrature, only: BSplineQuadrature, n_quad_exact
330 use m_common, only: set_union
331 use m_tensorprod_basis, only: TensorProdSpace, cartesian_to_linear, size, actv_interval_inds
332 use m_tensorprod_quadrature, only: TensorProdQuadrature, evaluate_at_nodes, DEFAULT_EXTRA_QUADRATURE_POINTS
333
334 implicit none
335
336 type(TensorProdMat), intent(out) :: M
337 type(TensorProdSpace), intent(in) :: tp_space1, tp_space2
338 procedure(user_function_3d_interface), optional :: coordfun
339 integer, intent(in), optional :: n_quad_extra
340
341 integer :: row, col ! Linearly indexed B-spline indices: row = first basis, col = second basis
342 integer :: degree_x1, degree_y1, degree_z1 ! B-spline degree per direction of the first B-spline basis
343 integer :: degree_x2, degree_y2, degree_z2 ! B-spline degree per direction of the second B-spline basis
344 integer :: n_quad1, n_quad2, n_quad3 ! Number of quadrature points in each direction
345 integer :: nr_intervals_x, nr_intervals_y, nr_intervals_z ! Number of intervals per direction (must be equal for both B-spline
346 ! spaces)
347 logical :: equal_spaces_x, equal_spaces_y, equal_spaces_z ! Whether the B-spline spaces are equal in the respective direction
348 logical :: equal_spaces_all ! Whether the B-spline spaces are equal in each direction
349 integer :: i1, j1, k1, i2, j2, k2 ! The B-spline indices (for the first and second B-spline basis)
350 integer :: ii, jj, kk ! The interval indices
351 integer :: qdx, qdy, qdz ! The quadrature point indices
352 integer :: i1_minus_ii, j1_minus_jj, k1_minus_kk ! Equal to i1 - ii, j1 - jj, k1 - kk
353 integer :: i2_minus_ii, j2_minus_jj, k2_minus_kk ! Equal to i2 - ii, j2 - jj, k2 - kk
354 integer, dimension(0:MAX_BSPLINE_DEGREE) :: i1_internals, j1_internals, k1_internals
355 integer, dimension(0:MAX_BSPLINE_DEGREE) :: i2_internals, j2_internals, k2_internals
356 integer :: i1_internal, j1_internal, k1_internal ! The internal B-spline indices (there are at most degree + 1 different
357 ! B-splines per direction)
358 integer :: i2_internal, j2_internal, k2_internal ! Same for the second B-spline basis
359 integer :: n_quad_extra_ ! The number of extra quadrature points
360
361 integer :: local_cdx_i, local_cdx_j, local_cdx_k
362 logical, dimension(0:MAX_BSPLINE_DEGREE) :: bspline_mirrored_x1s, bspline_mirrored_y1s, bspline_mirrored_z1s
363 logical, dimension(0:MAX_BSPLINE_DEGREE) :: bspline_mirrored_x2s, bspline_mirrored_y2s, bspline_mirrored_z2s
364 logical :: bspline_mirrored_x1, bspline_mirrored_y1, bspline_mirrored_z1
365 logical :: bspline_mirrored_x2, bspline_mirrored_y2, bspline_mirrored_z2
366 logical :: is_periodic_x, is_periodic_y, is_periodic_z
367 integer :: i2_minus_ii_max, j2_minus_jj_max, k2_minus_kk_max
368
369 1957 real(wp), allocatable :: coordfun_vals(:, :, :) !> Values of coordfun at the quadrature points of the current interval
370 real(wp) :: val !> Value of the mass matrix
371 real(wp), allocatable :: vals(:, :, :, :, :, :) !> Values of the mass matrix per interval (i2-ii, i1-ii, j2-jj, j1-jj, k2-kk,
372 !> k1-kk)
373 integer, allocatable :: cols(:, :, :) !> Column indices corresponding to vals (i1-ii, j1-jj, k1-kk)
374 integer, allocatable :: rows(:, :, :) !> Row indices corresponding to vals (i2-ii, j2-jj, k2-kk)
375
376 real(wp), allocatable :: bspline1_vals_x(:), bspline1_vals_y(:), bspline1_vals_z(:)
377 real(wp), allocatable :: bspline_vals_x(:), bspline_vals_y(:), bspline_vals_z(:)
378
379 real(wp), allocatable :: partial_sum_z(:, :), partial_sum_yz(:)
380
381
5/6
✓ Branch 3 → 4 taken 5871 times.
✓ Branch 3 → 5 taken 1957 times.
✓ Branch 6 → 7 taken 5871 times.
✓ Branch 6 → 8 taken 1957 times.
✗ Branch 8 → 9 not taken.
✓ Branch 8 → 10 taken 1957 times.
13699 type(TensorProdQuadrature) :: tp_space1_quad, tp_space2_quad
382
383 ! The number of quadrature points in each diriection is determined such that
384 ! the quadrature is exact for the given B-spline spaces
385 1957 n_quad1 = n_quad_exact(tp_space1%spaces(1), tp_space2%spaces(1))
386 1957 n_quad2 = n_quad_exact(tp_space1%spaces(2), tp_space2%spaces(2))
387 1957 n_quad3 = n_quad_exact(tp_space1%spaces(3), tp_space2%spaces(3))
388
2/2
✓ Branch 14 → 15 taken 1110 times.
✓ Branch 14 → 18 taken 847 times.
1957 if (present(coordfun)) then
389
2/2
✓ Branch 15 → 16 taken 289 times.
✓ Branch 15 → 17 taken 821 times.
1110 if (present(n_quad_extra)) then
390 289 n_quad_extra_ = n_quad_extra
391 else
392 n_quad_extra_ = DEFAULT_EXTRA_QUADRATURE_POINTS
393 end if
394
395 ! The number of quadrature points in each direction is increased
396 1110 n_quad1 = n_quad1 + n_quad_extra_
397 1110 n_quad2 = n_quad2 + n_quad_extra_
398 1110 n_quad3 = n_quad3 + n_quad_extra_
399 end if
400
401
2/2
✓ Branch 19 → 20 taken 5871 times.
✓ Branch 19 → 26 taken 1957 times.
7828 do i1 = 1, size(tp_space1%spaces)
402
4/6
✓ Branch 20 → 21 taken 2712 times.
✓ Branch 20 → 22 taken 3159 times.
✗ Branch 21 → 22 not taken.
✓ Branch 21 → 23 taken 2712 times.
✗ Branch 23 → 24 not taken.
✓ Branch 23 → 25 taken 5871 times.
9030 if (.not. (tp_space1%spaces(i1)%is_periodic .eqv. tp_space2%spaces(i1)%is_periodic .and. &
403 1957 & tp_space1%spaces(i1)%nr_intervals == tp_space2%spaces(i1)%nr_intervals)) then
404 error stop "B-spline spaces must have the same periodicity and number of intervals for the computation of the mass matrix"
405 end if
406 end do
407
408 ! We precompute the value of the B-spline basis functions at the quadrature points
409 tp_space1_quad = TensorProdQuadrature(BSplineQuadrature(tp_space1%spaces(1), n_quad1), &
410 BSplineQuadrature(tp_space1%spaces(2), n_quad2), &
411
29/50
✓ Branch 33 → 34 taken 5871 times.
✓ Branch 33 → 41 taken 1957 times.
✗ Branch 34 → 35 not taken.
✓ Branch 34 → 36 taken 5871 times.
✗ Branch 36 → 37 not taken.
✓ Branch 36 → 38 taken 5871 times.
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 40 taken 5871 times.
✓ Branch 41 → 42 taken 1957 times.
✗ Branch 41 → 43 not taken.
✓ Branch 43 → 44 taken 1957 times.
✗ Branch 43 → 45 not taken.
✓ Branch 45 → 46 taken 1957 times.
✗ Branch 45 → 47 not taken.
✓ Branch 47 → 48 taken 1957 times.
✗ Branch 47 → 49 not taken.
✓ Branch 49 → 50 taken 1957 times.
✗ Branch 49 → 51 not taken.
✓ Branch 51 → 52 taken 1957 times.
✗ Branch 51 → 53 not taken.
✓ Branch 53 → 54 taken 1957 times.
✗ Branch 53 → 55 not taken.
✓ Branch 55 → 56 taken 1957 times.
✗ Branch 55 → 57 not taken.
✓ Branch 57 → 58 taken 1957 times.
✗ Branch 57 → 59 not taken.
✓ Branch 60 → 61 taken 5871 times.
✓ Branch 60 → 71 taken 1957 times.
✓ Branch 61 → 62 taken 5871 times.
✗ Branch 61 → 63 not taken.
✓ Branch 64 → 65 taken 5871 times.
✗ Branch 64 → 66 not taken.
✓ Branch 67 → 68 taken 5871 times.
✗ Branch 67 → 69 not taken.
✓ Branch 71 → 72 taken 5871 times.
✓ Branch 71 → 79 taken 1957 times.
✗ Branch 72 → 73 not taken.
✓ Branch 72 → 74 taken 5871 times.
✗ Branch 74 → 75 not taken.
✓ Branch 74 → 76 taken 5871 times.
✗ Branch 76 → 77 not taken.
✓ Branch 76 → 78 taken 5871 times.
✓ Branch 80 → 81 taken 5871 times.
✓ Branch 80 → 88 taken 1957 times.
✓ Branch 81 → 82 taken 5871 times.
✗ Branch 81 → 83 not taken.
✓ Branch 83 → 84 taken 5871 times.
✗ Branch 83 → 85 not taken.
✓ Branch 85 → 86 taken 5871 times.
✗ Branch 85 → 87 not taken.
27398 BSplineQuadrature(tp_space1%spaces(3), n_quad3))
412 tp_space2_quad = TensorProdQuadrature(BSplineQuadrature(tp_space2%spaces(1), n_quad1), &
413 BSplineQuadrature(tp_space2%spaces(2), n_quad2), &
414
29/50
✓ Branch 94 → 95 taken 5871 times.
✓ Branch 94 → 102 taken 1957 times.
✗ Branch 95 → 96 not taken.
✓ Branch 95 → 97 taken 5871 times.
✗ Branch 97 → 98 not taken.
✓ Branch 97 → 99 taken 5871 times.
✗ Branch 99 → 100 not taken.
✓ Branch 99 → 101 taken 5871 times.
✓ Branch 102 → 103 taken 1957 times.
✗ Branch 102 → 104 not taken.
✓ Branch 104 → 105 taken 1957 times.
✗ Branch 104 → 106 not taken.
✓ Branch 106 → 107 taken 1957 times.
✗ Branch 106 → 108 not taken.
✓ Branch 108 → 109 taken 1957 times.
✗ Branch 108 → 110 not taken.
✓ Branch 110 → 111 taken 1957 times.
✗ Branch 110 → 112 not taken.
✓ Branch 112 → 113 taken 1957 times.
✗ Branch 112 → 114 not taken.
✓ Branch 114 → 115 taken 1957 times.
✗ Branch 114 → 116 not taken.
✓ Branch 116 → 117 taken 1957 times.
✗ Branch 116 → 118 not taken.
✓ Branch 118 → 119 taken 1957 times.
✗ Branch 118 → 120 not taken.
✓ Branch 121 → 122 taken 5871 times.
✓ Branch 121 → 132 taken 1957 times.
✓ Branch 122 → 123 taken 5871 times.
✗ Branch 122 → 124 not taken.
✓ Branch 125 → 126 taken 5871 times.
✗ Branch 125 → 127 not taken.
✓ Branch 128 → 129 taken 5871 times.
✗ Branch 128 → 130 not taken.
✓ Branch 132 → 133 taken 5871 times.
✓ Branch 132 → 140 taken 1957 times.
✗ Branch 133 → 134 not taken.
✓ Branch 133 → 135 taken 5871 times.
✗ Branch 135 → 136 not taken.
✓ Branch 135 → 137 taken 5871 times.
✗ Branch 137 → 138 not taken.
✓ Branch 137 → 139 taken 5871 times.
✓ Branch 141 → 142 taken 5871 times.
✓ Branch 141 → 149 taken 1957 times.
✓ Branch 142 → 143 taken 5871 times.
✗ Branch 142 → 144 not taken.
✓ Branch 144 → 145 taken 5871 times.
✗ Branch 144 → 146 not taken.
✓ Branch 146 → 147 taken 5871 times.
✗ Branch 146 → 148 not taken.
27398 BSplineQuadrature(tp_space2%spaces(3), n_quad3))
415
416 1957 degree_x1 = tp_space1%spaces(1)%degree
417 1957 degree_y1 = tp_space1%spaces(2)%degree
418 1957 degree_z1 = tp_space1%spaces(3)%degree
419 1957 degree_x2 = tp_space2%spaces(1)%degree
420 1957 degree_y2 = tp_space2%spaces(2)%degree
421 1957 degree_z2 = tp_space2%spaces(3)%degree
422
423 1957 is_periodic_x = tp_space1%spaces(1)%is_periodic
424 1957 is_periodic_y = tp_space1%spaces(2)%is_periodic
425 1957 is_periodic_z = tp_space1%spaces(3)%is_periodic
426
427 ! Per interval the symmetry can be exploited and computations can be re-used
428 ! TODO even if the spaces are not equal, we can still exploit symmetry: the only difference is scaling, which is a constant
429 ! TODO factor per B-spline
430 equal_spaces_x = size(tp_space1%spaces(1)) == size(tp_space2%spaces(1)) .and. &
431 (tp_space1%spaces(1)%is_scaled .eqv. tp_space2%spaces(1)%is_scaled) .and. &
432 (tp_space1%spaces(1)%is_periodic .eqv. tp_space2%spaces(1)%is_periodic) .and. &
433
4/6
✓ Branch 149 → 150 taken 1853 times.
✓ Branch 149 → 152 taken 104 times.
✓ Branch 150 → 151 taken 1853 times.
✗ Branch 150 → 152 not taken.
✗ Branch 151 → 152 not taken.
✓ Branch 151 → 153 taken 1853 times.
1957 degree_x1 == degree_x2
434
435 equal_spaces_y = size(tp_space1%spaces(2)) == size(tp_space2%spaces(2)) .and. &
436 (tp_space1%spaces(2)%is_scaled .eqv. tp_space2%spaces(2)%is_scaled) .and. &
437 (tp_space1%spaces(2)%is_periodic .eqv. tp_space2%spaces(2)%is_periodic) .and. &
438
4/6
✓ Branch 153 → 154 taken 1957 times.
✗ Branch 153 → 156 not taken.
✓ Branch 154 → 155 taken 1853 times.
✓ Branch 154 → 156 taken 104 times.
✗ Branch 155 → 156 not taken.
✓ Branch 155 → 157 taken 1853 times.
1957 degree_y1 == degree_y2
439
440 equal_spaces_z = size(tp_space1%spaces(3)) == size(tp_space2%spaces(3)) .and. &
441 (tp_space1%spaces(3)%is_scaled .eqv. tp_space2%spaces(3)%is_scaled) .and. &
442 (tp_space1%spaces(3)%is_periodic .eqv. tp_space2%spaces(3)%is_periodic) .and. &
443
4/6
✓ Branch 157 → 158 taken 1957 times.
✗ Branch 157 → 160 not taken.
✓ Branch 158 → 159 taken 1853 times.
✓ Branch 158 → 160 taken 104 times.
✗ Branch 159 → 160 not taken.
✓ Branch 159 → 161 taken 1853 times.
1957 degree_z1 == degree_z2
444
445
3/4
✓ Branch 161 → 162 taken 1801 times.
✓ Branch 161 → 163 taken 156 times.
✗ Branch 162 → 163 not taken.
✓ Branch 162 → 164 taken 1801 times.
1957 equal_spaces_all = equal_spaces_x .and. equal_spaces_y .and. equal_spaces_z
446
447 i2_minus_ii_max = degree_x2
448 j2_minus_jj_max = degree_y2
449 k2_minus_kk_max = degree_z2
450
451
13/20
✓ Branch 164 → 165 taken 119 times.
✓ Branch 164 → 166 taken 1838 times.
✓ Branch 166 → 165 taken 1838 times.
✗ Branch 166 → 167 not taken.
✓ Branch 167 → 168 taken 47 times.
✓ Branch 167 → 169 taken 1910 times.
✓ Branch 169 → 168 taken 1910 times.
✗ Branch 169 → 170 not taken.
✗ Branch 170 → 171 not taken.
✓ Branch 170 → 172 taken 1957 times.
✓ Branch 172 → 171 taken 1957 times.
✗ Branch 172 → 173 not taken.
✓ Branch 173 → 174 taken 1957 times.
✗ Branch 173 → 175 not taken.
✓ Branch 175 → 176 taken 1791 times.
✓ Branch 175 → 177 taken 166 times.
✗ Branch 177 → 178 not taken.
✓ Branch 177 → 179 taken 1957 times.
✗ Branch 179 → 180 not taken.
✓ Branch 179 → 181 taken 1957 times.
9785 allocate (coordfun_vals(n_quad3, n_quad2, n_quad1))
452
21/32
✗ Branch 181 → 182 not taken.
✓ Branch 181 → 183 taken 1957 times.
✓ Branch 183 → 182 taken 1957 times.
✗ Branch 183 → 184 not taken.
✗ Branch 184 → 185 not taken.
✓ Branch 184 → 186 taken 1957 times.
✓ Branch 186 → 185 taken 1957 times.
✗ Branch 186 → 187 not taken.
✓ Branch 187 → 188 taken 47 times.
✓ Branch 187 → 189 taken 1910 times.
✓ Branch 189 → 188 taken 1910 times.
✗ Branch 189 → 190 not taken.
✓ Branch 190 → 191 taken 47 times.
✓ Branch 190 → 192 taken 1910 times.
✓ Branch 192 → 191 taken 1910 times.
✗ Branch 192 → 193 not taken.
✓ Branch 193 → 194 taken 336 times.
✓ Branch 193 → 195 taken 1621 times.
✓ Branch 195 → 194 taken 1621 times.
✗ Branch 195 → 196 not taken.
✓ Branch 196 → 197 taken 336 times.
✓ Branch 196 → 198 taken 1621 times.
✓ Branch 198 → 197 taken 1621 times.
✗ Branch 198 → 199 not taken.
✓ Branch 199 → 200 taken 1957 times.
✗ Branch 199 → 201 not taken.
✓ Branch 201 → 202 taken 1574 times.
✓ Branch 201 → 203 taken 383 times.
✗ Branch 203 → 204 not taken.
✓ Branch 203 → 205 taken 1957 times.
✗ Branch 205 → 206 not taken.
✓ Branch 205 → 207 taken 1957 times.
15656 allocate (vals(0:degree_x2, 0:degree_x1, 0:degree_y2, 0:degree_y1, 0:degree_z2, 0:degree_z1))
453
13/20
✗ Branch 207 → 208 not taken.
✓ Branch 207 → 209 taken 1957 times.
✓ Branch 209 → 208 taken 1957 times.
✗ Branch 209 → 210 not taken.
✓ Branch 210 → 211 taken 47 times.
✓ Branch 210 → 212 taken 1910 times.
✓ Branch 212 → 211 taken 1910 times.
✗ Branch 212 → 213 not taken.
✓ Branch 213 → 214 taken 336 times.
✓ Branch 213 → 215 taken 1621 times.
✓ Branch 215 → 214 taken 1621 times.
✗ Branch 215 → 216 not taken.
✓ Branch 216 → 217 taken 1957 times.
✗ Branch 216 → 218 not taken.
✓ Branch 218 → 219 taken 1574 times.
✓ Branch 218 → 220 taken 383 times.
✗ Branch 220 → 221 not taken.
✓ Branch 220 → 222 taken 1957 times.
✗ Branch 222 → 223 not taken.
✓ Branch 222 → 224 taken 1957 times.
9785 allocate (cols(0:degree_x2, 0:degree_y2, 0:degree_z2))
454
13/20
✗ Branch 224 → 225 not taken.
✓ Branch 224 → 226 taken 1957 times.
✓ Branch 226 → 225 taken 1957 times.
✗ Branch 226 → 227 not taken.
✓ Branch 227 → 228 taken 47 times.
✓ Branch 227 → 229 taken 1910 times.
✓ Branch 229 → 228 taken 1910 times.
✗ Branch 229 → 230 not taken.
✓ Branch 230 → 231 taken 336 times.
✓ Branch 230 → 232 taken 1621 times.
✓ Branch 232 → 231 taken 1621 times.
✗ Branch 232 → 233 not taken.
✓ Branch 233 → 234 taken 1957 times.
✗ Branch 233 → 235 not taken.
✓ Branch 235 → 236 taken 1574 times.
✓ Branch 235 → 237 taken 383 times.
✗ Branch 237 → 238 not taken.
✓ Branch 237 → 239 taken 1957 times.
✗ Branch 239 → 240 not taken.
✓ Branch 239 → 241 taken 1957 times.
9785 allocate (rows(0:degree_x1, 0:degree_y1, 0:degree_z1))
455
456
6/12
✗ Branch 241 → 242 not taken.
✓ Branch 241 → 243 taken 1957 times.
✓ Branch 243 → 242 taken 1957 times.
✗ Branch 243 → 244 not taken.
✓ Branch 244 → 245 taken 1957 times.
✗ Branch 244 → 246 not taken.
✓ Branch 246 → 247 taken 1957 times.
✗ Branch 246 → 248 not taken.
✗ Branch 248 → 249 not taken.
✓ Branch 248 → 250 taken 1957 times.
✗ Branch 250 → 251 not taken.
✓ Branch 250 → 252 taken 1957 times.
5871 allocate (bspline1_vals_x(n_quad1))
457
8/12
✓ Branch 252 → 253 taken 47 times.
✓ Branch 252 → 254 taken 1910 times.
✓ Branch 254 → 253 taken 1910 times.
✗ Branch 254 → 255 not taken.
✓ Branch 255 → 256 taken 1957 times.
✗ Branch 255 → 257 not taken.
✓ Branch 257 → 258 taken 1910 times.
✓ Branch 257 → 259 taken 47 times.
✗ Branch 259 → 260 not taken.
✓ Branch 259 → 261 taken 1957 times.
✗ Branch 261 → 262 not taken.
✓ Branch 261 → 263 taken 1957 times.
5871 allocate (bspline1_vals_y(n_quad2))
458
8/12
✓ Branch 263 → 264 taken 119 times.
✓ Branch 263 → 265 taken 1838 times.
✓ Branch 265 → 264 taken 1838 times.
✗ Branch 265 → 266 not taken.
✓ Branch 266 → 267 taken 1957 times.
✗ Branch 266 → 268 not taken.
✓ Branch 268 → 269 taken 1838 times.
✓ Branch 268 → 270 taken 119 times.
✗ Branch 270 → 271 not taken.
✓ Branch 270 → 272 taken 1957 times.
✗ Branch 272 → 273 not taken.
✓ Branch 272 → 274 taken 1957 times.
5871 allocate (bspline1_vals_z(n_quad3))
459
460
4/8
✗ Branch 274 → 275 not taken.
✓ Branch 274 → 276 taken 1957 times.
✓ Branch 276 → 275 taken 1957 times.
✗ Branch 276 → 277 not taken.
✗ Branch 277 → 278 not taken.
✓ Branch 277 → 279 taken 1957 times.
✗ Branch 279 → 280 not taken.
✓ Branch 279 → 281 taken 1957 times.
3914 allocate (bspline_vals_x(n_quad1))
461
5/8
✓ Branch 281 → 282 taken 47 times.
✓ Branch 281 → 283 taken 1910 times.
✓ Branch 283 → 282 taken 1910 times.
✗ Branch 283 → 284 not taken.
✗ Branch 284 → 285 not taken.
✓ Branch 284 → 286 taken 1957 times.
✗ Branch 286 → 287 not taken.
✓ Branch 286 → 288 taken 1957 times.
3914 allocate (bspline_vals_y(n_quad2))
462
5/8
✓ Branch 288 → 289 taken 119 times.
✓ Branch 288 → 290 taken 1838 times.
✓ Branch 290 → 289 taken 1838 times.
✗ Branch 290 → 291 not taken.
✗ Branch 291 → 292 not taken.
✓ Branch 291 → 293 taken 1957 times.
✗ Branch 293 → 294 not taken.
✓ Branch 293 → 295 taken 1957 times.
3914 allocate (bspline_vals_z(n_quad3))
463
464
10/16
✓ Branch 295 → 296 taken 47 times.
✓ Branch 295 → 297 taken 1910 times.
✓ Branch 297 → 296 taken 1910 times.
✗ Branch 297 → 298 not taken.
✗ Branch 298 → 299 not taken.
✓ Branch 298 → 300 taken 1957 times.
✓ Branch 300 → 299 taken 1957 times.
✗ Branch 300 → 301 not taken.
✓ Branch 301 → 302 taken 1957 times.
✗ Branch 301 → 303 not taken.
✓ Branch 303 → 304 taken 1910 times.
✓ Branch 303 → 305 taken 47 times.
✗ Branch 305 → 306 not taken.
✓ Branch 305 → 307 taken 1957 times.
✗ Branch 307 → 308 not taken.
✓ Branch 307 → 309 taken 1957 times.
7828 allocate (partial_sum_z(n_quad2, n_quad1))
465
4/8
✗ Branch 309 → 310 not taken.
✓ Branch 309 → 311 taken 1957 times.
✓ Branch 311 → 310 taken 1957 times.
✗ Branch 311 → 312 not taken.
✗ Branch 312 → 313 not taken.
✓ Branch 312 → 314 taken 1957 times.
✗ Branch 314 → 315 not taken.
✓ Branch 314 → 316 taken 1957 times.
3914 allocate (partial_sum_yz(n_quad1))
466
467 1957 call M%init(tp_space1, tp_space2, symmetric=equal_spaces_all, spd=equal_spaces_all)
468
469 ! Loop over all intervals first such that the coordfun is evaluated only once per quadrature point
470
2/2
✓ Branch 318 → 319 taken 10843 times.
✓ Branch 318 → 666 taken 1957 times.
12800 do kk = tp_space1%rank_resp_intervals%k0, tp_space1%rank_resp_intervals%k1
471
2/2
✓ Branch 320 → 321 taken 138108 times.
✓ Branch 320 → 664 taken 10843 times.
150908 do jj = tp_space1%rank_resp_intervals%j0, tp_space1%rank_resp_intervals%j1
472
2/2
✓ Branch 322 → 323 taken 1880819 times.
✓ Branch 322 → 662 taken 138108 times.
2029770 do ii = tp_space1%rank_resp_intervals%i0, tp_space1%rank_resp_intervals%i1
473
2/2
✓ Branch 323 → 324 taken 1492907 times.
✓ Branch 323 → 325 taken 387912 times.
1880819 if (present(coordfun)) then
474 1492907 call evaluate_at_nodes(tp_space1_quad, coordfun, ii, jj, kk, coordfun_vals)
475 end if
476
477
2/2
✓ Branch 326 → 327 taken 4821682 times.
✓ Branch 326 → 329 taken 1880819 times.
6702501 do k1_minus_kk = 0, degree_z1
478 call get_internal_and_actual_inds(tp_space1%spaces(3), kk, k1_minus_kk, k1, k1_internals(k1_minus_kk), &
479 6702501 & bspline_mirrored_z1s(k1_minus_kk))
480 ! TODO simplify get_internal_and_actual_inds: don't need k1 from there
481 end do
482
2/2
✓ Branch 331 → 332 taken 4821682 times.
✓ Branch 331 → 334 taken 1880819 times.
6702501 do k2_minus_kk = 0, degree_z2
483 call get_internal_and_actual_inds(tp_space2%spaces(3), kk, k2_minus_kk, k2, k2_internals(k2_minus_kk), &
484 6702501 & bspline_mirrored_z2s(k2_minus_kk))
485 end do
486
2/2
✓ Branch 336 → 337 taken 6750069 times.
✓ Branch 336 → 339 taken 1880819 times.
8630888 do j1_minus_jj = 0, degree_y1
487 call get_internal_and_actual_inds(tp_space1%spaces(2), jj, j1_minus_jj, j1, j1_internals(j1_minus_jj), &
488 8630888 & bspline_mirrored_y1s(j1_minus_jj))
489 end do
490
2/2
✓ Branch 341 → 342 taken 6750069 times.
✓ Branch 341 → 344 taken 1880819 times.
8630888 do j2_minus_jj = 0, degree_y2
491 call get_internal_and_actual_inds(tp_space2%spaces(2), jj, j2_minus_jj, j2, j2_internals(j2_minus_jj), &
492 8630888 & bspline_mirrored_y2s(j2_minus_jj))
493 end do
494
2/2
✓ Branch 346 → 347 taken 6955325 times.
✓ Branch 346 → 349 taken 1880819 times.
8836144 do i1_minus_ii = 0, degree_x1
495 call get_internal_and_actual_inds(tp_space1%spaces(1), ii, i1_minus_ii, i1, i1_internals(i1_minus_ii), &
496 8836144 & bspline_mirrored_x1s(i1_minus_ii))
497 end do
498
2/2
✓ Branch 351 → 352 taken 6955325 times.
✓ Branch 351 → 354 taken 1880819 times.
8836144 do i2_minus_ii = 0, degree_x2
499 call get_internal_and_actual_inds(tp_space2%spaces(1), ii, i2_minus_ii, i2, i2_internals(i2_minus_ii), &
500 8836144 & bspline_mirrored_x2s(i2_minus_ii))
501 end do
502
503 ! Loop over the third index (k or z-direction) of the B-spline basis function of the first basis
504 ! The B-splines have compact support: they contribute only if k1 is in the interval [kk, kk+degree_z1]
505
2/2
✓ Branch 356 → 357 taken 4821682 times.
✓ Branch 356 → 548 taken 1880819 times.
6702501 do k1_minus_kk = 0, degree_z1
506 4821682 k1_internal = k1_internals(k1_minus_kk)
507 4821682 bspline_mirrored_z1 = bspline_mirrored_z1s(k1_minus_kk)
508
509
2/2
✓ Branch 357 → 358 taken 29928 times.
✓ Branch 357 → 366 taken 4791754 times.
4821682 if (bspline_mirrored_z1) then
510
4/8
✓ Branch 358 → 359 taken 29928 times.
✗ Branch 358 → 360 not taken.
✗ Branch 359 → 360 not taken.
✓ Branch 359 → 363 taken 29928 times.
✗ Branch 360 → 361 not taken.
✗ Branch 360 → 362 not taken.
✓ Branch 364 → 365 taken 176612 times.
✓ Branch 364 → 374 taken 29928 times.
236468 bspline1_vals_z = tp_space1_quad%bspline_quad(3)%bspline_at_nodes(k1_internal, degree_z1 - k1_minus_kk)%data(n_quad3:1:-1)
511 else
512
4/8
✓ Branch 366 → 367 taken 4791754 times.
✗ Branch 366 → 368 not taken.
✗ Branch 367 → 368 not taken.
✓ Branch 367 → 371 taken 4791754 times.
✗ Branch 368 → 369 not taken.
✗ Branch 368 → 370 not taken.
✓ Branch 372 → 373 taken 24051962 times.
✓ Branch 372 → 374 taken 4791754 times.
33635470 bspline1_vals_z = tp_space1_quad%bspline_quad(3)%bspline_at_nodes(k1_internal, k1_minus_kk)%data(1:n_quad3)
513 end if
514
3/4
✗ Branch 374 → 375 not taken.
✓ Branch 374 → 376 taken 4821682 times.
✓ Branch 377 → 378 taken 24228574 times.
✓ Branch 377 → 379 taken 4821682 times.
33871938 bspline1_vals_z = bspline1_vals_z * tp_space1_quad%bspline_quad(3)%weights
515
516
2/2
✓ Branch 379 → 380 taken 4148434 times.
✓ Branch 379 → 381 taken 673248 times.
4821682 if (equal_spaces_z) then
517 ! If the spaces are equal in the z direction, we can exploit symmetry
518 ! by computing "the lower triangular part" of the contribution to the
519 ! mass matrix. This yields a factor 1/2 of computation per equal direction
520 k2_minus_kk_max = k1_minus_kk
521 end if
522
523 ! Loop over the third index (k or z-direction) of the B-spline basis function of the second basis
524
2/2
✓ Branch 382 → 383 taken 11311754 times.
✓ Branch 382 → 546 taken 4821682 times.
18014255 do k2_minus_kk = 0, k2_minus_kk_max
525 11311754 k2_internal = k2_internals(k2_minus_kk)
526 11311754 bspline_mirrored_z2 = bspline_mirrored_z2s(k2_minus_kk)
527 ! Given kk (interval in the z-direction) and k2_minus_kk (B-spline index in the z-direction of the second basis relative
528 ! to kk), determine the actual B-spline index k2 (in the z-direction), the internally used B-spline index k2_internal,
529 ! and determine if the B-spline is mirrored
530
2/2
✓ Branch 383 → 384 taken 55736 times.
✓ Branch 383 → 392 taken 11256018 times.
11311754 if (bspline_mirrored_z2) then
531 ! This means that the clamped B-spline is at the right-hand side boundary, and is
532 ! therefore defined as the mirror image of a B-spline at the left-hand side boundary
533 bspline_vals_z = bspline1_vals_z * &
534
4/8
✓ Branch 384 → 385 taken 55736 times.
✗ Branch 384 → 386 not taken.
✗ Branch 385 → 386 not taken.
✓ Branch 385 → 389 taken 55736 times.
✗ Branch 386 → 387 not taken.
✗ Branch 386 → 388 not taken.
✓ Branch 390 → 391 taken 339088 times.
✓ Branch 390 → 400 taken 55736 times.
450560 & tp_space2_quad%bspline_quad(3)%bspline_at_nodes(k2_internal, degree_z2 - k2_minus_kk)%data(n_quad3:1:-1)
535 else
536 bspline_vals_z = bspline1_vals_z * &
537
4/8
✓ Branch 392 → 393 taken 11256018 times.
✗ Branch 392 → 394 not taken.
✗ Branch 393 → 394 not taken.
✓ Branch 393 → 397 taken 11256018 times.
✗ Branch 394 → 395 not taken.
✗ Branch 394 → 396 not taken.
✓ Branch 398 → 399 taken 59810098 times.
✓ Branch 398 → 400 taken 11256018 times.
82322134 & tp_space2_quad%bspline_quad(3)%bspline_at_nodes(k2_internal, k2_minus_kk)%data(1:n_quad3)
538 end if
539
540
2/2
✓ Branch 400 → 401 taken 7808092 times.
✓ Branch 400 → 413 taken 3503662 times.
11311754 if (present(coordfun)) then
541 ! We precompute the partial sum over the z-direction: this reduces the
542 ! complexity from (assuming n_quad = degree for simplicity)
543 ! nr_intervals^3 * degree^9
544 ! to
545 ! nr_intervals^3 * degree^7
546
2/2
✓ Branch 402 → 403 taken 45398980 times.
✓ Branch 402 → 412 taken 7808092 times.
53207072 do qdx = 1, n_quad1
547
2/2
✓ Branch 404 → 405 taken 268460413 times.
✓ Branch 404 → 410 taken 45398980 times.
321667485 do qdy = 1, n_quad2
548 val = 0._wp
549
2/2
✓ Branch 406 → 407 taken 1552040359 times.
✓ Branch 406 → 408 taken 268460413 times.
1820500772 do qdz = 1, n_quad3
550 1820500772 val = val + bspline_vals_z(qdz) * coordfun_vals(qdz, qdy, qdx)
551 end do
552 313859393 partial_sum_z(qdy, qdx) = val
553 end do
554 end do
555 else
556 ! partial_sum_z(qdx, qdy) = sum(bspline_vals_z)
557 end if
558
559 ! Loop over the second index (j or y-direction) of the B-spline basis function of the first basis
560
2/2
✓ Branch 414 → 415 taken 41117636 times.
✓ Branch 414 → 527 taken 11311754 times.
52429390 do j1_minus_jj = 0, degree_y1
561 41117636 j1_internal = j1_internals(j1_minus_jj)
562 41117636 bspline_mirrored_y1 = bspline_mirrored_y1s(j1_minus_jj)
563
2/2
✓ Branch 415 → 416 taken 4211212 times.
✓ Branch 415 → 424 taken 36906424 times.
41117636 if (bspline_mirrored_y1) then
564 bspline1_vals_y = &
565
4/8
✓ Branch 416 → 417 taken 4211212 times.
✗ Branch 416 → 418 not taken.
✗ Branch 417 → 418 not taken.
✓ Branch 417 → 421 taken 4211212 times.
✗ Branch 418 → 419 not taken.
✗ Branch 418 → 420 not taken.
✓ Branch 422 → 423 taken 17605408 times.
✓ Branch 422 → 432 taken 4211212 times.
26027832 & tp_space1_quad%bspline_quad(2)%bspline_at_nodes(j1_internal, degree_y1 - j1_minus_jj)%data(n_quad2:1:-1)
566 else
567
4/8
✓ Branch 424 → 425 taken 36906424 times.
✗ Branch 424 → 426 not taken.
✗ Branch 425 → 426 not taken.
✓ Branch 425 → 429 taken 36906424 times.
✗ Branch 426 → 427 not taken.
✗ Branch 426 → 428 not taken.
✓ Branch 430 → 431 taken 202279240 times.
✓ Branch 430 → 432 taken 36906424 times.
276092088 bspline1_vals_y = tp_space1_quad%bspline_quad(2)%bspline_at_nodes(j1_internal, j1_minus_jj)%data(1:n_quad2)
568 end if
569
3/4
✗ Branch 432 → 433 not taken.
✓ Branch 432 → 434 taken 41117636 times.
✓ Branch 435 → 436 taken 219884648 times.
✓ Branch 435 → 437 taken 41117636 times.
302119920 bspline1_vals_y = bspline1_vals_y * tp_space1_quad%bspline_quad(2)%weights
570
571
2/2
✓ Branch 437 → 438 taken 35976308 times.
✓ Branch 437 → 439 taken 5141328 times.
41117636 if (equal_spaces_y) then
572 j2_minus_jj_max = j1_minus_jj
573 end if
574
575 ! Loop over the second index (j or y-direction) of the B-spline basis function of the second basis
576
2/2
✓ Branch 440 → 441 taken 105861448 times.
✓ Branch 440 → 525 taken 41117636 times.
158290838 do j2_minus_jj = 0, j2_minus_jj_max
577 105861448 j2_internal = j2_internals(j2_minus_jj)
578 105861448 bspline_mirrored_y2 = bspline_mirrored_y2s(j2_minus_jj)
579
2/2
✓ Branch 441 → 442 taken 6962900 times.
✓ Branch 441 → 450 taken 98898548 times.
105861448 if (bspline_mirrored_y2) then
580 bspline_vals_y = bspline1_vals_y * &
581
4/8
✓ Branch 442 → 443 taken 6962900 times.
✗ Branch 442 → 444 not taken.
✗ Branch 443 → 444 not taken.
✓ Branch 443 → 447 taken 6962900 times.
✗ Branch 444 → 445 not taken.
✗ Branch 444 → 446 not taken.
✓ Branch 448 → 449 taken 31103866 times.
✓ Branch 448 → 458 taken 6962900 times.
45029666 & tp_space2_quad%bspline_quad(2)%bspline_at_nodes(j2_internal, degree_y2 - j2_minus_jj)%data(n_quad2:1:-1)
582 else
583 bspline_vals_y = bspline1_vals_y * &
584
4/8
✓ Branch 450 → 451 taken 98898548 times.
✗ Branch 450 → 452 not taken.
✗ Branch 451 → 452 not taken.
✓ Branch 451 → 455 taken 98898548 times.
✗ Branch 452 → 453 not taken.
✗ Branch 452 → 454 not taken.
✓ Branch 456 → 457 taken 558613872 times.
✓ Branch 456 → 458 taken 98898548 times.
756410968 & tp_space2_quad%bspline_quad(2)%bspline_at_nodes(j2_internal, j2_minus_jj)%data(1:n_quad2)
585 end if
586
587
2/2
✓ Branch 458 → 459 taken 77038934 times.
✓ Branch 458 → 466 taken 28822514 times.
105861448 if (present(coordfun)) then
588 ! We precompute the partial sum over the yz-directions
589
2/2
✓ Branch 460 → 461 taken 468059234 times.
✓ Branch 460 → 473 taken 77038934 times.
545098168 do qdx = 1, n_quad1
590 val = 0._wp
591
2/2
✓ Branch 462 → 463 taken 2954472438 times.
✓ Branch 462 → 464 taken 468059234 times.
3422531672 do qdy = 1, n_quad2
592 3422531672 val = val + bspline_vals_y(qdy) * partial_sum_z(qdy, qdx)
593 end do
594 545098168 partial_sum_yz(qdx) = val
595 end do
596 else
597
6/6
✓ Branch 466 → 467 taken 110741624 times.
✓ Branch 466 → 468 taken 28822514 times.
✓ Branch 468 → 469 taken 133751366 times.
✓ Branch 468 → 470 taken 28822514 times.
✓ Branch 471 → 472 taken 118126266 times.
✓ Branch 471 → 474 taken 28822514 times.
391441770 partial_sum_yz = sum(bspline_vals_y) * sum(bspline_vals_z)
598 end if
599
600 ! Loop over the first index (i or x-direction) of the B-spline basis function of the first basis
601
2/2
✓ Branch 475 → 476 taken 422628348 times.
✓ Branch 475 → 514 taken 105861448 times.
528489796 do i1_minus_ii = 0, degree_x1
602 422628348 i1_internal = i1_internals(i1_minus_ii)
603 422628348 bspline_mirrored_x1 = bspline_mirrored_x1s(i1_minus_ii)
604
2/2
✓ Branch 476 → 477 taken 44540094 times.
✓ Branch 476 → 485 taken 378088254 times.
422628348 if (bspline_mirrored_x1) then
605 bspline1_vals_x = &
606
4/8
✓ Branch 477 → 478 taken 44540094 times.
✗ Branch 477 → 479 not taken.
✗ Branch 478 → 479 not taken.
✓ Branch 478 → 482 taken 44540094 times.
✗ Branch 479 → 480 not taken.
✗ Branch 479 → 481 not taken.
✓ Branch 483 → 484 taken 285299450 times.
✓ Branch 483 → 493 taken 44540094 times.
374379638 & tp_space1_quad%bspline_quad(1)%bspline_at_nodes(i1_internal, degree_x1 - i1_minus_ii)%data(n_quad1:1:-1)
607 else
608 bspline1_vals_x = &
609
4/8
✓ Branch 485 → 486 taken 378088254 times.
✗ Branch 485 → 487 not taken.
✗ Branch 486 → 487 not taken.
✓ Branch 486 → 490 taken 378088254 times.
✗ Branch 487 → 488 not taken.
✗ Branch 487 → 489 not taken.
✓ Branch 491 → 492 taken 2129254674 times.
✓ Branch 491 → 493 taken 378088254 times.
2885431182 & tp_space1_quad%bspline_quad(1)%bspline_at_nodes(i1_internal, i1_minus_ii)%data(1:n_quad1)
610 end if
611
3/4
✗ Branch 493 → 494 not taken.
✓ Branch 493 → 495 taken 422628348 times.
✓ Branch 496 → 497 taken 2414554124 times.
✓ Branch 496 → 498 taken 422628348 times.
3259810820 bspline1_vals_x = bspline1_vals_x * tp_space1_quad%bspline_quad(1)%weights
612
613
2/2
✓ Branch 498 → 499 taken 370869468 times.
✓ Branch 498 → 500 taken 51758880 times.
422628348 if (equal_spaces_x) then
614 i2_minus_ii_max = i1_minus_ii
615 end if
616
617 ! Loop over the first index (i or x-direction) of the B-spline basis function of the second basis
618
2/2
✓ Branch 501 → 502 taken 1159948880 times.
✓ Branch 501 → 512 taken 422628348 times.
1688438676 do i2_minus_ii = 0, i2_minus_ii_max
619 1159948880 i2_internal = i2_internals(i2_minus_ii)
620 1159948880 bspline_mirrored_x2 = bspline_mirrored_x2s(i2_minus_ii)
621
622 val = 0._wp
623
2/2
✓ Branch 502 → 503 taken 94656044 times.
✓ Branch 502 → 506 taken 1065292836 times.
1159948880 if (bspline_mirrored_x2) then
624
2/2
✓ Branch 504 → 505 taken 615998044 times.
✓ Branch 504 → 510 taken 94656044 times.
710654088 do qdx = 1, n_quad1
625 val = val + bspline1_vals_x(qdx) * partial_sum_yz(qdx) * &
626 710654088 & tp_space2_quad%bspline_quad(1)%bspline_at_nodes(i2_internal, degree_x2 - i2_minus_ii)%data(n_quad1 + 1 - qdx)
627 end do
628 else
629 do qdx = 1, n_quad1
630 val = val + bspline1_vals_x(qdx) * partial_sum_yz(qdx) * &
631 & tp_space2_quad%bspline_quad(1)%bspline_at_nodes(i2_internal, i2_minus_ii)%data(qdx)
632 end do
633 end if
634
635 1582577228 vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk) = val
636
637 end do ! i2_minus_ii
638 end do ! i1_minus_ii
639
640
2/2
✓ Branch 515 → 516 taken 91220488 times.
✓ Branch 515 → 524 taken 14640960 times.
146979084 if (equal_spaces_x) then
641 ! Impose symmetry of the local values (set remainder of i2_minus_ii > i1_minus_ii)
642
2/2
✓ Branch 517 → 518 taken 370869468 times.
✓ Branch 517 → 523 taken 91220488 times.
462089956 do i1_minus_ii = 0, degree_x1
643
2/2
✓ Branch 519 → 520 taken 604909172 times.
✓ Branch 519 → 521 taken 370869468 times.
1066999128 do i2_minus_ii = i1_minus_ii + 1, degree_x2
644 ! SWAPPED SWAPPED
645 vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk) = &
646 975778640 vals(i1_minus_ii, i2_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk)
647 end do
648 end do
649 end if
650 end do ! j2_minus_jj
651 end do ! j1_minus_jj
652
653
2/2
✓ Branch 528 → 529 taken 9728490 times.
✓ Branch 528 → 545 taken 1583264 times.
16133436 if (equal_spaces_y) then
654 ! Impose symmetry of the local values (set remainder of j2_minus_jj > j1_minus_jj)
655
2/2
✓ Branch 530 → 531 taken 35976308 times.
✓ Branch 530 → 544 taken 9728490 times.
45704798 do j1_minus_jj = 0, degree_y1
656
2/2
✓ Branch 532 → 533 taken 52813780 times.
✓ Branch 532 → 542 taken 35976308 times.
98518578 do j2_minus_jj = j1_minus_jj + 1, degree_y2
657
2/2
✓ Branch 534 → 535 taken 217284475 times.
✓ Branch 534 → 540 taken 52813780 times.
306074563 do i1_minus_ii = 0, degree_x1
658
2/2
✓ Branch 536 → 537 taken 933658933 times.
✓ Branch 536 → 538 taken 217284475 times.
1203757188 do i2_minus_ii = 0, degree_x2
659 ! SWAPPED SWAPPED
660 vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk) = &
661 1150943408 vals(i2_minus_ii, i1_minus_ii, j1_minus_jj, j2_minus_jj, k2_minus_kk, k1_minus_kk)
662 end do
663 end do
664 end do
665 end do
666 end if
667 end do ! k2_minus_kk
668 end do ! k1_minus_kk
669
670
2/2
✓ Branch 549 → 550 taken 1589427 times.
✓ Branch 549 → 574 taken 291392 times.
1880819 if (equal_spaces_z) then
671 ! Impose symmetry of the local values (set remainder of k2_minus_kk > k1_minus_kk)
672
2/2
✓ Branch 551 → 552 taken 4148434 times.
✓ Branch 551 → 573 taken 1589427 times.
5737861 do k1_minus_kk = 0, degree_z1
673
2/2
✓ Branch 553 → 554 taken 5455032 times.
✓ Branch 553 → 571 taken 4148434 times.
11192893 do k2_minus_kk = k1_minus_kk + 1, degree_z2
674
2/2
✓ Branch 555 → 556 taken 20387740 times.
✓ Branch 555 → 569 taken 5455032 times.
29991206 do j1_minus_jj = 0, degree_y1
675
2/2
✓ Branch 557 → 558 taken 80521628 times.
✓ Branch 557 → 567 taken 20387740 times.
106364400 do j2_minus_jj = 0, degree_y2
676
2/2
✓ Branch 559 → 560 taken 333047232 times.
✓ Branch 559 → 565 taken 80521628 times.
433956600 do i1_minus_ii = 0, degree_x1
677
2/2
✓ Branch 561 → 562 taken 1433793152 times.
✓ Branch 561 → 563 taken 333047232 times.
1847362012 do i2_minus_ii = 0, degree_x2
678 ! SWAPPED SWAPPED
679 vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk) = &
680 1766840384 vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k1_minus_kk, k2_minus_kk)
681 end do
682 end do
683 end do
684 end do
685 end do
686 end do
687 end if
688
689 ! Determine row indices
690
2/2
✓ Branch 575 → 576 taken 4821682 times.
✓ Branch 575 → 592 taken 1880819 times.
6702501 do k1_minus_kk = 0, degree_z1
691 4821682 k1 = k1_minus_kk + kk
692
2/2
✓ Branch 576 → 577 taken 4419168 times.
✓ Branch 576 → 578 taken 402514 times.
4821682 if (is_periodic_z) k1 = modulo(k1, tp_space1%spaces(3)%nr_bsplines)
693
2/2
✓ Branch 579 → 580 taken 17215112 times.
✓ Branch 579 → 590 taken 4821682 times.
23917613 do j1_minus_jj = 0, degree_y1
694 17215112 j1 = j1_minus_jj + jj
695
2/2
✓ Branch 580 → 581 taken 12405552 times.
✓ Branch 580 → 582 taken 4809560 times.
17215112 if (is_periodic_y) j1 = modulo(j1, tp_space1%spaces(2)%nr_bsplines)
696
2/2
✓ Branch 583 → 584 taken 66455345 times.
✓ Branch 583 → 588 taken 17215112 times.
88492139 do i1_minus_ii = 0, degree_x1
697 66455345 i1 = i1_minus_ii + ii
698
2/2
✓ Branch 584 → 585 taken 16615880 times.
✓ Branch 584 → 586 taken 49839465 times.
66455345 if (is_periodic_x) i1 = modulo(i1, tp_space1%spaces(1)%nr_bsplines)
699
700 66455345 call cartesian_to_linear(tp_space1, row, i1, j1, k1)
701 66455345 rows(i1_minus_ii, j1_minus_jj, k1_minus_kk) = row
702 83670457 M%rows(i1_minus_ii + ii, j1_minus_jj + jj, k1_minus_kk + kk) = row
703 end do
704 end do
705 end do
706
707 ! Determine column indices
708
2/2
✓ Branch 593 → 594 taken 1443731 times.
✓ Branch 593 → 616 taken 437088 times.
1880819 if (equal_spaces_all) then
709
10/22
✓ Branch 594 → 595 taken 1443731 times.
✗ Branch 594 → 598 not taken.
✓ Branch 595 → 596 taken 1443731 times.
✗ Branch 595 → 598 not taken.
✓ Branch 596 → 597 taken 1443731 times.
✗ Branch 596 → 598 not taken.
✗ Branch 597 → 598 not taken.
✓ Branch 597 → 607 taken 1443731 times.
✗ Branch 598 → 599 not taken.
✗ Branch 598 → 600 not taken.
✗ Branch 600 → 601 not taken.
✗ Branch 600 → 602 not taken.
✗ Branch 602 → 603 not taken.
✗ Branch 602 → 604 not taken.
✗ Branch 604 → 605 not taken.
✗ Branch 604 → 606 not taken.
✓ Branch 608 → 609 taken 3770066 times.
✓ Branch 608 → 635 taken 1443731 times.
✓ Branch 610 → 611 taken 13810152 times.
✓ Branch 610 → 615 taken 3770066 times.
✓ Branch 612 → 613 taken 54656561 times.
✓ Branch 612 → 614 taken 13810152 times.
75124241 cols = rows
710 else
711
2/2
✓ Branch 617 → 618 taken 1051616 times.
✓ Branch 617 → 634 taken 437088 times.
1488704 do k2_minus_kk = 0, degree_z2
712 1051616 k2 = k2_minus_kk + kk
713
1/2
✓ Branch 618 → 619 taken 1051616 times.
✗ Branch 618 → 620 not taken.
1051616 if (is_periodic_z) k2 = modulo(k2, tp_space2%spaces(3)%nr_bsplines)
714
2/2
✓ Branch 621 → 622 taken 3404960 times.
✓ Branch 621 → 632 taken 1051616 times.
4893664 do j2_minus_jj = 0, degree_y2
715 3404960 j2 = j2_minus_jj + jj
716
1/2
✓ Branch 622 → 623 taken 3404960 times.
✗ Branch 622 → 624 not taken.
3404960 if (is_periodic_y) j2 = modulo(j2, tp_space2%spaces(2)%nr_bsplines)
717
2/2
✓ Branch 625 → 626 taken 11798784 times.
✓ Branch 625 → 630 taken 3404960 times.
16255360 do i2_minus_ii = 0, degree_x2
718 11798784 i2 = i2_minus_ii + ii
719
1/2
✗ Branch 626 → 627 not taken.
✓ Branch 626 → 628 taken 11798784 times.
11798784 if (is_periodic_x) i2 = modulo(i2, tp_space2%spaces(1)%nr_bsplines)
720
721 11798784 call cartesian_to_linear(tp_space2, col, i2, j2, k2)
722 15203744 cols(i2_minus_ii, j2_minus_jj, k2_minus_kk) = col
723 end do
724 end do
725 end do
726 end if
727
728 ! Accumulate entries of the mass matrix
729
2/2
✓ Branch 636 → 637 taken 4821682 times.
✓ Branch 636 → 660 taken 1880819 times.
6840609 do k1_minus_kk = 0, degree_z1
730
2/2
✓ Branch 638 → 639 taken 16766786 times.
✓ Branch 638 → 658 taken 4821682 times.
23469287 do k2_minus_kk = 0, degree_z2
731 16766786 local_cdx_k = degree_z1 + k2_minus_kk - k1_minus_kk
732 16766786 k1 = k1_minus_kk + kk
733
2/2
✓ Branch 640 → 641 taken 61505376 times.
✓ Branch 640 → 656 taken 16766786 times.
83093844 do j1_minus_jj = 0, degree_y1
734
2/2
✓ Branch 642 → 643 taken 239196856 times.
✓ Branch 642 → 654 taken 61505376 times.
317469018 do j2_minus_jj = 0, degree_y2
735 239196856 local_cdx_j = degree_y1 + j2_minus_jj - j1_minus_jj
736 239196856 j1 = j1_minus_jj + jj
737
2/2
✓ Branch 644 → 645 taken 972960055 times.
✓ Branch 644 → 652 taken 239196856 times.
1273662287 do i1_minus_ii = 0, degree_x1
738 do i2_minus_ii = 0, degree_x2
739 4132310137 i1 = i1_minus_ii + ii
740 4132310137 local_cdx_i = degree_x1 + i2_minus_ii - i1_minus_ii
741
742
2/2
✓ Branch 647 → 648 taken 437809477 times.
✓ Branch 647 → 649 taken 3694500660 times.
4132310137 if (M%cols(local_cdx_i, local_cdx_j, local_cdx_k, i1, j1, k1) == -1) then
743 437809477 M%cols(local_cdx_i, local_cdx_j, local_cdx_k, i1, j1, k1) = cols(i2_minus_ii, j2_minus_jj, k2_minus_kk)
744 end if
745
746 M%vals(local_cdx_i, local_cdx_j, local_cdx_k, i1, j1, k1) = M%vals(local_cdx_i, local_cdx_j, local_cdx_k, i1, j1, k1) + &
747 & vals(i2_minus_ii, i1_minus_ii, j2_minus_jj, j1_minus_jj, k2_minus_kk, k1_minus_kk)
748 end do
749 end do
750 end do
751 end do
752 end do
753 end do
754 end do ! ii
755 end do ! jj
756 end do ! kk
757
758
1/2
✗ Branch 667 → 668 not taken.
✓ Branch 667 → 669 taken 1957 times.
1957 deallocate (coordfun_vals)
759
1/2
✗ Branch 669 → 670 not taken.
✓ Branch 669 → 671 taken 1957 times.
1957 deallocate (vals)
760
1/2
✗ Branch 671 → 672 not taken.
✓ Branch 671 → 673 taken 1957 times.
1957 deallocate (cols)
761
1/2
✗ Branch 673 → 674 not taken.
✓ Branch 673 → 675 taken 1957 times.
1957 deallocate (rows)
762
763
1/2
✗ Branch 675 → 676 not taken.
✓ Branch 675 → 677 taken 1957 times.
1957 deallocate (bspline1_vals_x)
764
1/2
✗ Branch 677 → 678 not taken.
✓ Branch 677 → 679 taken 1957 times.
1957 deallocate (bspline1_vals_y)
765
1/2
✗ Branch 679 → 680 not taken.
✓ Branch 679 → 681 taken 1957 times.
1957 deallocate (bspline1_vals_z)
766
767
1/2
✗ Branch 681 → 682 not taken.
✓ Branch 681 → 683 taken 1957 times.
1957 deallocate (bspline_vals_x)
768
1/2
✗ Branch 683 → 684 not taken.
✓ Branch 683 → 685 taken 1957 times.
1957 deallocate (bspline_vals_y)
769
1/2
✗ Branch 685 → 686 not taken.
✓ Branch 685 → 687 taken 1957 times.
1957 deallocate (bspline_vals_z)
770
771
1/2
✗ Branch 687 → 688 not taken.
✓ Branch 687 → 689 taken 1957 times.
1957 deallocate (partial_sum_z)
772
1/2
✗ Branch 689 → 690 not taken.
✓ Branch 689 → 691 taken 1957 times.
1957 deallocate (partial_sum_yz)
773
774 1957 call tp_space1_quad%destroy()
775 1957 call tp_space2_quad%destroy()
776
777
10/16
✓ Branch 694 → 695 taken 5871 times.
✓ Branch 694 → 702 taken 1957 times.
✗ Branch 695 → 696 not taken.
✓ Branch 695 → 697 taken 5871 times.
✗ Branch 697 → 698 not taken.
✓ Branch 697 → 699 taken 5871 times.
✗ Branch 699 → 700 not taken.
✓ Branch 699 → 701 taken 5871 times.
✓ Branch 702 → 703 taken 5871 times.
✓ Branch 702 → 710 taken 1957 times.
✗ Branch 703 → 704 not taken.
✓ Branch 703 → 705 taken 5871 times.
✗ Branch 705 → 706 not taken.
✓ Branch 705 → 707 taken 5871 times.
✗ Branch 707 → 708 not taken.
✓ Branch 707 → 709 taken 5871 times.
13699 end subroutine mass_matrix_3d
778
779
4/30
None:
✓ Branch 27 → 28 taken 5871 times.
✓ Branch 27 → 29 taken 1957 times.
✓ Branch 88 → 89 taken 5871 times.
✓ Branch 88 → 90 taken 1957 times.
__m_tensorprod_matrix_MOD___copy_m_tensorprod_matrix_Tensorprodmat:
✗ Branch 2 → 3 not taken.
✗ Branch 2 → 12 not taken.
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 5 not taken.
✗ Branch 6 → 7 not taken.
✗ Branch 6 → 8 not taken.
✗ Branch 9 → 10 not taken.
✗ Branch 9 → 11 not taken.
__m_tensorprod_matrix_MOD___final_m_tensorprod_matrix_Tensorprodmat:
✗ Branch 3 → 4 not taken.
✗ Branch 3 → 7 not taken.
✗ Branch 4 → 5 not taken.
✗ Branch 4 → 6 not taken.
✗ Branch 8 → 9 not taken.
✗ Branch 8 → 23 not taken.
✗ Branch 10 → 11 not taken.
✗ Branch 10 → 12 not taken.
✗ Branch 12 → 13 not taken.
✗ Branch 12 → 18 not taken.
✗ Branch 13 → 14 not taken.
✗ Branch 13 → 15 not taken.
✗ Branch 15 → 16 not taken.
✗ Branch 15 → 17 not taken.
✗ Branch 18 → 19 not taken.
✗ Branch 18 → 22 not taken.
✗ Branch 19 → 20 not taken.
✗ Branch 19 → 21 not taken.
15656 end module m_tensorprod_matrix
780