GCC Code Coverage Report

Directory: src/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 39.5% 77 / 0 / 195
Functions: 41.2% 7 / 0 / 17
Branches: 34.2% 128 / 0 / 374
other/s_qr_factorisation.f90
Line Branch Exec Source
1 submodule(m_qr_factorisation) s_qr_factorisation
2 implicit none
3 contains
4
5 !> @brief Computes the orthogonal matrix Q from QR factorization using Householder reflectors (out-of-place)
6
2/4
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 3 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 3 times.
 3 module subroutine orthogonalize_out_of_place(Q, A, success)
7 real(wp), intent(out) :: Q(1:, 1:)
8 real(wp), intent(in) :: A(1:, 1:)
9 logical, intent(out), optional :: success
10
11 3 real(wp), allocatable :: A_(:, :)
12
13
8/16
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 3 times.
✓ Branch 8 → 7 taken 3 times.
✗ Branch 8 → 9 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 3 times.
✓ Branch 11 → 10 taken 3 times.
✗ Branch 11 → 12 not taken.
✓ Branch 12 → 13 taken 3 times.
✗ Branch 12 → 14 not taken.
✓ Branch 14 → 15 taken 3 times.
✗ Branch 14 → 16 not taken.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 3 times.
✗ Branch 18 → 19 not taken.
✓ Branch 18 → 20 taken 3 times.
 12 allocate (A_(size(A, 1), size(A, 2)))
14
7/12
✓ Branch 20 → 21 taken 3 times.
✗ Branch 20 → 23 not taken.
✓ Branch 21 → 22 taken 3 times.
✗ Branch 21 → 23 not taken.
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 27 taken 3 times.
✗ Branch 23 → 24 not taken.
✗ Branch 23 → 25 not taken.
✓ Branch 28 → 29 taken 9 times.
✓ Branch 28 → 33 taken 3 times.
✓ Branch 30 → 31 taken 33 times.
✓ Branch 30 → 32 taken 9 times.
 48 A_ = A
15 3 call orthogonalize_inplace(A_, success)
16
4/4
✓ Branch 35 → 36 taken 9 times.
✓ Branch 35 → 40 taken 3 times.
✓ Branch 37 → 38 taken 33 times.
✓ Branch 37 → 39 taken 9 times.
 45 Q = A_
17
1/2
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 42 taken 3 times.
 3 deallocate (A_)
18 3 end subroutine orthogonalize_out_of_place
19
20 !> @brief Computes the orthogonal matrix Q from QR factorization using Householder reflectors (in-place)
21
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 1660 times.
 1660 module subroutine orthogonalize_inplace(A, success)
22 real(wp), intent(inout) :: A(1:, 1:)
23 logical, intent(out), optional :: success
24
25 real(wp), allocatable :: Q(:, :)
26 3320 real(wp) :: alpha(size(A, 2))
27 ✗ integer :: p(size(A, 2))
28 logical :: success_
29 integer :: rows, cols, i, j, qCols
30 character(len=256) :: msg
31
32 1660 rows = size(A, 1)
33 1660 cols = size(A, 2)
34 1660 qCols = min(rows, cols)
35
36
1/2
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 22 taken 1660 times.
 1660 if (size(A, 2) /= qCols) then
37 ✗ if (qCols == 0) then
38 ✗ if (present(success)) success = .true.
39 ✗ return
40 end if
41 ✗ if (present(success)) then
42 ✗ success = .false.
43 ✗ return
44 else
45 write (msg, '(A,I0,A,I0,A,I0,A,I0)') "ERROR in orthogonalize_inplace: A has wrong dimensions. " &
46 ✗ //"Expected ", rows, " x ", qCols, " but got ", size(A, 1), " x ", size(A, 2)
47 ✗ error stop msg
48 end if
49 end if
50
51
10/16
✗ Branch 22 → 23 not taken.
✓ Branch 22 → 25 taken 1660 times.
✓ Branch 25 → 23 taken 1660 times.
✗ Branch 25 → 26 not taken.
✓ Branch 26 → 27 taken 588 times.
✓ Branch 26 → 28 taken 1072 times.
✓ Branch 28 → 27 taken 1072 times.
✗ Branch 28 → 29 not taken.
✓ Branch 29 → 30 taken 1660 times.
✗ Branch 29 → 31 not taken.
✓ Branch 31 → 32 taken 1072 times.
✓ Branch 31 → 33 taken 588 times.
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 1660 times.
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 1660 times.
 6640 allocate (Q(rows, qCols))
52 1660 call qr_factorise_inplace(A, alpha, p, success_)
53
54
1/2
✗ Branch 38 → 39 not taken.
✓ Branch 38 → 53 taken 1660 times.
 1660 if (.not. success_) then
55 ✗ if (present(success)) then
56 ✗ success = .false.
57 ✗ deallocate (Q)
58 ✗ return
59 else
60 ✗ deallocate (Q)
61 ✗ write (msg, '(A,I0,A,I0,A)') "orthogonalize_inplace: QR factorization failed because matrix (", &
62 ✗ rows, " x ", cols, ") not full rank."
63 ✗ error stop msg
64 end if
65 end if
66
67
4/4
✓ Branch 53 → 54 taken 7541 times.
✓ Branch 53 → 58 taken 1660 times.
✓ Branch 55 → 56 taken 1964676 times.
✓ Branch 55 → 57 taken 7541 times.
 1973877 Q = 0.0_wp
68
2/2
✓ Branch 59 → 60 taken 7541 times.
✓ Branch 59 → 61 taken 1660 times.
 9201 do i = 1, qCols
69 9201 Q(i, i) = 1.0_wp
70 end do
71
72
2/2
✓ Branch 62 → 63 taken 7541 times.
✓ Branch 62 → 69 taken 1660 times.
 9201 do j = 1, qCols
73
2/2
✓ Branch 64 → 65 taken 113941 times.
✓ Branch 64 → 67 taken 7541 times.
 123142 do i = min(cols, rows), 1, -1
74 121482 call apply_householder_transformation(Q(:, j), i, A, alpha, p)
75 end do
76 end do
77
78
4/4
✓ Branch 70 → 71 taken 7541 times.
✓ Branch 70 → 75 taken 1660 times.
✓ Branch 72 → 73 taken 1964676 times.
✓ Branch 72 → 74 taken 7541 times.
 1973877 A = Q
79
1/2
✗ Branch 75 → 76 not taken.
✓ Branch 75 → 77 taken 1660 times.
 1660 if (present(success)) success = .true.
80
1/2
✗ Branch 77 → 78 not taken.
✓ Branch 77 → 79 taken 1660 times.
 1660 deallocate (Q)
81 end subroutine orthogonalize_inplace
82
83 !> @brief Computes the orthogonal matrix Q from QR factorization (in-place, 4D array)
84 1656 module subroutine orthogonalize_inplace_4d(A, success)
85 real(wp), contiguous, target, intent(inout) :: A(:, :, :, :)
86 logical, intent(out), optional :: success
87
88 real(wp), pointer :: A_flat(:, :)
89 integer :: dim1, dim2, dim3, dim4, rows, cols
90
91 1656 dim1 = size(A, 1)
92 1656 dim2 = size(A, 2)
93 1656 dim3 = size(A, 3)
94 1656 dim4 = size(A, 4)
95 1656 rows = dim1 * dim2 * dim3
96 cols = dim4
97
98 1656 A_flat(1:rows, 1:cols) => A(:, :, :, :)
99 1656 call orthogonalize(A_flat, success)
100 1656 end subroutine orthogonalize_inplace_4d
101
102 !> @brief Computes the column space of a dense matrix A
103 ✗ module subroutine dense_column_space(column_space, rank, A, tolerance, sparsify)
104 real(wp), intent(in) :: A(1:, 1:)
105 real(wp), intent(out) :: column_space(:, :)
106 integer, intent(out) :: rank
107 real(wp), intent(in), optional :: tolerance
108 logical, intent(in), optional :: sparsify
109
110 integer :: col, pCol, j
111 ✗ real(wp) :: tolerance_, norms(size(A, 2)), normX, v(size(A, 1))
112 real(wp), parameter :: COLUMN_HAS_BEEN_VISITED = -42.0
113 logical :: sparsify_
114
115 ✗ if (present(tolerance)) then
116 ✗ tolerance_ = tolerance
117 else
118 tolerance_ = 1.0e-10_wp
119 end if
120
121 ✗ if (present(sparsify)) then
122 ✗ sparsify_ = sparsify
123 else
124 sparsify_ = .false.
125 end if
126
127 ✗ do col = 1, size(A, 2)
128 ✗ if (sparsify_) then
129 ✗ norms(col) = size(A, 1) - count(abs(A(:, col)) > tolerance_)
130 else
131 ✗ norms(col) = dot_product(A(:, col), A(:, col))
132 end if
133 end do
134
135 ✗ rank = 0
136 ✗ column_space = 0.0_wp
137 ✗ do col = 1, size(A, 2)
138 ✗ pCol = my_maxloc(norms)
139 ✗ norms(pCol) = COLUMN_HAS_BEEN_VISITED
140
141 ✗ v = A(:, pCol)
142 ✗ do j = 1, col - 1
143 ✗ v = v - column_space(:, j) * dot_product(column_space(:, j), v)
144 end do
145 ✗ normX = norm2(v)
146
147 ✗ if (normX >= tolerance_) then
148 ✗ v = v / normX
149 ✗ rank = rank + 1
150 ✗ column_space(:, rank) = v
151 end if
152 end do
153 ✗ end subroutine dense_column_space
154
155 !> @brief Computes the rank of a dense matrix A
156
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 471 times.
 471 module integer function dense_matrix_rank(A, tolerance)
157 real(wp), intent(in) :: A(1:, 1:)
158 real(wp), intent(in), optional :: tolerance
159
160 942 integer :: p(size(A, 2))
161 942 real(wp) :: alpha(size(A, 2))
162 logical :: success
163 real(wp), allocatable :: A_(:, :)
164 real(wp) :: tolerance_
165
166
10/16
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 471 times.
✓ Branch 6 → 5 taken 471 times.
✗ Branch 6 → 7 not taken.
✓ Branch 7 → 8 taken 162 times.
✓ Branch 7 → 9 taken 309 times.
✓ Branch 9 → 8 taken 309 times.
✗ Branch 9 → 10 not taken.
✓ Branch 10 → 11 taken 471 times.
✗ Branch 10 → 12 not taken.
✓ Branch 12 → 13 taken 309 times.
✓ Branch 12 → 14 taken 162 times.
✗ Branch 14 → 15 not taken.
✓ Branch 14 → 16 taken 471 times.
✗ Branch 16 → 17 not taken.
✓ Branch 16 → 18 taken 471 times.
 1884 allocate (A_(size(A, 1), size(A, 2)))
167
7/12
✓ Branch 18 → 19 taken 471 times.
✗ Branch 18 → 21 not taken.
✓ Branch 19 → 20 taken 471 times.
✗ Branch 19 → 21 not taken.
✗ Branch 20 → 21 not taken.
✓ Branch 20 → 25 taken 471 times.
✗ Branch 21 → 22 not taken.
✗ Branch 21 → 23 not taken.
✓ Branch 26 → 27 taken 8386 times.
✓ Branch 26 → 31 taken 471 times.
✓ Branch 28 → 29 taken 11346394 times.
✓ Branch 28 → 30 taken 8386 times.
 11355722 A_ = A
168 471 call qr_factorise_inplace(A_, alpha, p, success)
169
170
1/2
✗ Branch 32 → 33 not taken.
✓ Branch 32 → 34 taken 471 times.
 471 if (present(tolerance)) then
171 ✗ tolerance_ = tolerance
172 else
173 tolerance_ = 1.0e-10_wp
174 end if
175
176
4/4
✓ Branch 35 → 36 taken 8386 times.
✓ Branch 35 → 39 taken 471 times.
✓ Branch 36 → 37 taken 8382 times.
✓ Branch 36 → 38 taken 4 times.
 8857 dense_matrix_rank = count(abs(alpha) > tolerance_)
177
1/2
✓ Branch 39 → 40 taken 471 times.
✗ Branch 39 → 41 not taken.
 471 end function dense_matrix_rank
178
179 !> @brief Computes the inverse of a dense square matrix A
180 ✗ module subroutine dense_matrix_inverse(mat)
181 real(wp), intent(inout) :: mat(1:, 1:)
182 integer :: n, m
183 logical :: success
184 ✗ real(wp) :: qrOper(size(mat, 1), size(mat, 2))
185
186 ✗ n = size(mat, 1)
187 ✗ m = size(mat, 2)
188
189 ✗ if (n /= m) then
190 write (*, '(A,I0,A,I0)') "ERROR in dense_matrix_inverse: matrix must be square. " &
191 ✗ //"Got dimensions ", n, " x ", m
192 ✗ error stop 1
193 end if
194
195 ✗ call qr_get_solution_operator_inplace(qrOper, mat, success)
196 ✗ if (.not. success) then
197 write (*, '(A,I0)') "ERROR in dense_matrix_inverse: matrix is singular (cannot be inverted). " &
198 ✗ //"Matrix size: ", n, " x ", m
199 ✗ error stop 1
200 end if
201
202 ✗ mat = qrOper
203 ✗ end subroutine dense_matrix_inverse
204
205 !> @brief Computes the in-place solution operator R^{-1} * Q^T
206 ✗ subroutine qr_get_solution_operator_inplace(qrOper, A, success)
207 real(wp), intent(out) :: qrOper(1:, 1:)
208 real(wp), intent(inout) :: A(1:, 1:)
209 logical, intent(out) :: success
210
211 integer :: ii, jj
212 ✗ real(wp) :: B(size(A, 1), size(A, 1))
213
214 ✗ do ii = 1, size(B, 1)
215 ✗ do jj = 1, size(B, 2)
216 ✗ if (ii == jj) then
217 ✗ B(ii, jj) = 1.0_wp
218 else
219 ✗ B(ii, jj) = 0.0_wp
220 end if
221 end do
222 end do
223
224 ✗ call qr_solve_inplace_multiple_rhs(qrOper, A, B, success)
225 ✗ end subroutine qr_get_solution_operator_inplace
226
227 !> @brief Computes the solution operator R^{-1} * Q^T (out-of-place)
228 ✗ subroutine qr_get_solution_operator(qrOper, A, success)
229 real(wp), intent(out) :: qrOper(1:, 1:)
230 real(wp), intent(in) :: A(1:, 1:)
231 logical, intent(out) :: success
232
233 ✗ real(wp) :: A_(size(A, 1), size(A, 2))
234
235 ✗ A_ = A
236 ✗ call qr_get_solution_operator_inplace(qrOper, A_, success)
237 ✗ end subroutine qr_get_solution_operator
238
239 !> @brief Solves a linear system Ax = b using QR factorisation (out-of-place)
240 ✗ module subroutine qr_solve_single_rhs(x, A, b, success, resNorm)
241 real(wp), intent(out) :: x(1:)
242 real(wp), intent(in) :: A(1:, 1:), b(1:)
243 logical, intent(out) :: success
244 real(wp), intent(out), optional :: resNorm
245
246 ✗ real(wp) :: A_(size(A, 1), size(A, 2)), b_(size(b, 1))
247
248 ✗ A_ = A
249 ✗ b_ = b
250 ✗ call qr_solve_inplace_single_rhs(x, A_, b_, success, resNorm)
251 ✗ end subroutine qr_solve_single_rhs
252
253 !> @brief Solves a linear system Ax = b using QR factorisation (in-place)
254 ✗ module subroutine qr_solve_inplace_single_rhs(x, A, b, success, resNorm)
255 real(wp), intent(out) :: x(1:)
256 real(wp), intent(inout) :: A(1:, 1:), b(1:)
257 logical, intent(out) :: success
258 real(wp), intent(out), optional :: resNorm
259
260 ✗ real(wp) :: alpha(size(A, 2))
261 ✗ integer :: p(size(A, 2))
262
263 ✗ success = .false.
264 ✗ if (size(A, 1) /= size(b, 1) .or. size(A, 2) /= size(x, 1)) return
265
266 ✗ call qr_factorise_inplace(A, alpha, p, success)
267 ✗ if (.not. success) return
268
269 ✗ call qr_solve_inplace_given_fractorisation(x, b, A, alpha, p, resNorm, present(resNorm))
270 end subroutine qr_solve_inplace_single_rhs
271
272 !> @brief Solves a linear system Ax = b for multiple right-hand sides (out-of-place)
273 ✗ module subroutine qr_solve_multiple_rhs(x, A, b, success, resNorm)
274 real(wp), intent(out) :: x(1:, 1:)
275 real(wp), intent(in) :: A(1:, 1:), b(1:, 1:)
276 logical, intent(out) :: success
277 real(wp), intent(out), optional :: resNorm(:)
278
279 ✗ real(wp) :: A_(size(A, 1), size(A, 2)), b_(size(b, 1), size(b, 2))
280
281 ✗ A_ = A
282 ✗ b_ = b
283 ✗ call qr_solve_inplace_multiple_rhs(x, A_, b_, success, resNorm)
284 ✗ end subroutine qr_solve_multiple_rhs
285
286 !> @brief Solves a linear system Ax = b for multiple right-hand sides (in-place)
287 ✗ module subroutine qr_solve_inplace_multiple_rhs(x, A, b, success, resNorm)
288 real(wp), intent(out) :: x(1:, 1:)
289 real(wp), intent(inout) :: A(1:, 1:), b(1:, 1:)
290 logical, intent(out) :: success
291 real(wp), intent(out), optional :: resNorm(:)
292
293 ✗ real(wp) :: alpha(size(A, 2))
294 ✗ integer :: p(size(A, 2))
295 integer :: sol, sols
296
297 ✗ success = .false.
298 ✗ if (size(A, 1) /= size(b, 1) .or. size(A, 2) /= size(x, 1) .or. size(x, 2) /= size(b, 2)) return
299
300 ✗ call qr_factorise_inplace(A, alpha, p, success)
301 ✗ if (.not. success) return
302
303 sols = size(x, 2)
304 ✗ do sol = 1, sols
305 ✗ call qr_solve_inplace_given_fractorisation(x(:, sol), b(:, sol), A, alpha, p, resNorm(sol), present(resNorm))
306 end do
307 end subroutine qr_solve_inplace_multiple_rhs
308
309 !> @brief In-place QR factorisation of a dense matrix A
310
3/6
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 2131 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 2131 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 2131 times.
 2131 subroutine qr_factorise_inplace(A, alpha, p, success)
311 real(wp), intent(inout) :: A(1:, 1:)
312 real(wp), intent(out) :: alpha(:)
313 integer, intent(out) :: p(:)
314 logical, intent(out) :: success
315
316 real(wp), parameter :: COLUMN_HAS_BEEN_VISITED = -42.0
317 2130 real(wp) :: norms(size(A, 2)), normX, beta
318 integer :: row, col, pCol, rows, cols
319
320 2131 rows = size(A, 1)
321 cols = size(A, 2)
322
323 ! Check array dimensions
324
2/4
✓ Branch 8 → 9 taken 2131 times.
✗ Branch 8 → 10 not taken.
✗ Branch 9 → 10 not taken.
✓ Branch 9 → 11 taken 2131 times.
 2131 if (size(alpha) /= cols .or. size(p) /= cols) then
325 ✗ error stop "qr_factorise_inplace: alpha and p must have size equal to number of columns of A"
326 end if
327
328
2/2
✓ Branch 11 → 12 taken 15927 times.
✓ Branch 11 → 20 taken 2131 times.
 18058 do col = 1, cols
329
6/6
✓ Branch 13 → 14 taken 13311070 times.
✓ Branch 13 → 19 taken 15927 times.
✓ Branch 14 → 15 taken 1653749 times.
✓ Branch 14 → 18 taken 11657321 times.
✓ Branch 15 → 16 taken 10210 times.
✓ Branch 15 → 17 taken 1643539 times.
 13329128 norms(col) = norm2(A(:, col), 1)**2
330 end do
331
2/2
✓ Branch 21 → 22 taken 15927 times.
✓ Branch 21 → 23 taken 2131 times.
 18058 alpha = 0.0_wp
332
2/2
✓ Branch 23 → 24 taken 15927 times.
✓ Branch 23 → 25 taken 2131 times.
 18058 p = 0
333
334 2131 success = .true.
335
2/2
✓ Branch 26 → 27 taken 15925 times.
✓ Branch 26 → 53 taken 2130 times.
 18055 do col = 1, cols
336 15925 pCol = my_maxloc(norms)
337
338
6/6
✓ Branch 28 → 29 taken 12907368 times.
✓ Branch 28 → 34 taken 15925 times.
✓ Branch 29 → 30 taken 8589679 times.
✓ Branch 29 → 33 taken 4317689 times.
✓ Branch 30 → 31 taken 8182 times.
✓ Branch 30 → 32 taken 8581497 times.
 12923293 normX = norm2(A(col:rows, pCol), 1)
339
3/4
✓ Branch 34 → 35 taken 15924 times.
✓ Branch 34 → 36 taken 1 time.
✗ Branch 35 → 36 not taken.
✓ Branch 35 → 37 taken 15924 times.
 15925 if (normX == 0 .or. norms(pCol) == COLUMN_HAS_BEEN_VISITED) then
340 1 success = .false.
341 return
342 end if
343 15924 p(col) = pCol
344
345 15924 alpha(col) = -sign(normX, A(col, pCol))
346 15924 beta = sqrt(normX * (normX + abs(A(col, pCol))))
347 15924 A(col, pCol) = A(col, pCol) - alpha(col)
348
2/2
✓ Branch 38 → 39 taken 12907365 times.
✓ Branch 38 → 40 taken 15924 times.
 12923289 A(col:rows, pCol) = A(col:rows, pCol) / beta
349
350
2/2
✓ Branch 40 → 41 taken 14544 times.
✓ Branch 40 → 52 taken 1380 times.
 18054 if (col < cols) then
351 14544 norms(pCol) = COLUMN_HAS_BEEN_VISITED
352
2/2
✓ Branch 42 → 43 taken 807392 times.
✓ Branch 42 → 51 taken 14544 times.
 821936 do row = 1, cols
353
2/2
✓ Branch 43 → 44 taken 403696 times.
✓ Branch 43 → 50 taken 403696 times.
 821936 if (norms(row) > 0) then
354
4/4
✓ Branch 45 → 46 taken 655332170 times.
✓ Branch 45 → 47 taken 403696 times.
✓ Branch 47 → 48 taken 655332170 times.
✓ Branch 47 → 49 taken 403696 times.
 1311068036 A(col:rows, row) = A(col:rows, row) - A(col:rows, pCol) * dot_product(A(col:rows, pCol), A(col:rows, row))
355 403696 norms(row) = norms(row) - A(col, row)**2
356 end if
357 end do
358 end if
359 end do
360 end subroutine qr_factorise_inplace
361
362 !> @brief Solves a linear system given the QR factorisation
363 ✗ subroutine qr_solve_inplace_given_fractorisation(x, b, A, alpha, p, resNorm, computeResNorm)
364 real(wp), intent(out) :: x(1:)
365 real(wp), intent(in) :: A(1:, 1:), alpha(1:)
366 real(wp), intent(inout) :: b(1:)
367 integer, intent(in) :: p(1:)
368 real(wp), intent(out) :: resNorm
369 logical, intent(in) :: computeResNorm
370
371 ✗ real(wp) :: backupB(size(b))
372 integer :: col, cols, rows
373
374 ✗ rows = size(A, 1)
375 ✗ cols = size(A, 2)
376 ✗ if (computeResNorm) backupB = b
377
378 ✗ do col = 1, cols
379 ✗ call apply_householder_transformation(b, col, A, alpha, p)
380 end do
381
382 ✗ call solve_upper_triangular(x, b, A, alpha, p)
383
384 ✗ if (computeResNorm) then
385 ✗ b(cols + 1:rows) = 0.0
386 ✗ do col = cols, 1, -1
387 ✗ call apply_householder_transformation(b, col, A, alpha, p)
388 end do
389 ✗ resNorm = norm2(b - backupB)
390 end if
391 ✗ end subroutine qr_solve_inplace_given_fractorisation
392
393 !> @brief Applies the Householder transformation to the vector b
394
3/6
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 113941 times.
✗ Branch 4 → 5 not taken.
✓ Branch 4 → 6 taken 113941 times.
✗ Branch 6 → 7 not taken.
✓ Branch 6 → 8 taken 113941 times.
 113941 subroutine apply_householder_transformation(b, idx, A, alpha, p)
395 real(wp), intent(inout) :: b(1:)
396 real(wp), intent(in) :: A(1:, 1:), alpha(1:)
397 integer, intent(in) :: p(1:), idx
398
399 integer :: rows
400
401 113941 rows = size(A, 1)
402
4/4
✓ Branch 9 → 10 taken 113941 times.
✓ Branch 9 → 11 taken 33801058 times.
✓ Branch 12 → 13 taken 33801058 times.
✓ Branch 12 → 14 taken 113941 times.
 67829998 b(idx:rows) = b(idx:rows) - A(idx:rows, p(idx)) * dot_product(A(idx:rows, p(idx)), b(idx:rows))
403 113941 end subroutine apply_householder_transformation
404
405 !> @brief Solves the upper triangular system R * x = b
406 ✗ subroutine solve_upper_triangular(x, b, A, alpha, p)
407 real(wp), intent(out) :: x(1:)
408 real(wp), intent(in) :: b(1:)
409 real(wp), intent(in) :: A(1:, 1:), alpha(1:)
410 integer, intent(in) :: p(1:)
411
412 integer :: col, cols
413
414 ✗ cols = size(A, 2)
415 ✗ do col = cols, 1, -1
416 ✗ x(p(col)) = (b(col) - dot_product(A(col, p(col + 1:cols)), x(p(col + 1:cols)))) / alpha(col)
417 end do
418 ✗ end subroutine solve_upper_triangular
419
420 !> @brief Custom implementation of maxloc for 1D real arrays (returns the index of the maximum element)
421 !>
422 !> @param[in] arr The input array
423 !>
424 !> @return The index of the maximum element in the array
425 !> @note: nvfortran seems to have a bug with maxloc, hence this custom implementation
426
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 15925 times.
 15925 pure integer function my_maxloc(arr)
427 real(wp), intent(in) :: arr(:)
428 integer :: i, maxIndex
429 real(wp) :: maxValue
430
431 maxValue = -huge(1.0_wp)
432 maxIndex = lbound(arr, 1)
433
3/4
✓ Branch 4 → 5 taken 15925 times.
✗ Branch 4 → 6 not taken.
✓ Branch 7 → 8 taken 823319 times.
✓ Branch 7 → 11 taken 15925 times.
 855169 do i = lbound(arr, 1), ubound(arr, 1)
434
2/2
✓ Branch 8 → 9 taken 43805 times.
✓ Branch 8 → 10 taken 779514 times.
 839244 if (arr(i) > maxValue) then
435 maxValue = arr(i)
436 maxIndex = i
437 end if
438 end do
439 my_maxloc = maxIndex
440 15925 end function my_maxloc
441
442 end submodule s_qr_factorisation
443