Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains procedure implementation of the generic interfaces of [pm_matrixChol](@ref pm_matrixChol).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Thursday 01:00 AM, September 23, 2021, Dallas, TX
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if CK_ENABLED
28 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
29 : #define TYPE_KIND complex(TKC)
30 : #define GET_CONJG(X) conjg(X)
31 : #define GET_RE(x) x%re
32 : #elif RK_ENABLED
33 : real(TKC), parameter :: ZERO = 0._TKC, ONE = 1._TKC
34 : #define TYPE_KIND real(TKC)
35 : #define GET_CONJG(X) X
36 : #define GET_RE(x) x
37 : #else
38 : #error "Unrecognized interface."
39 : #endif
40 : !%%%%%%%%%%%%%%%%
41 : #if setChoLow_ENABLED
42 : !%%%%%%%%%%%%%%%%
43 :
44 : real(TKC) :: summ
45 : integer(IK) :: idim
46 3083 : do idim = 1_IK, ndim
47 7338 : summ = mat(idim,idim) - dot_product(mat(idim,1:idim-1), mat(idim,1:idim-1))
48 3083 : if (0._TKC < summ) then
49 2280 : dia(idim) = sqrt(summ)
50 23100 : mat(idim+1:ndim,idim) = (mat(idim,idim+1:ndim) - matmul(mat(idim+1:ndim,1:idim-1), mat(idim,1:idim-1))) / dia(idim)
51 : else
52 400 : dia(1) = -idim
53 400 : return
54 : end if
55 : end do
56 :
57 : !%%%%%%%%%%%%%%%%%
58 : #elif getMatChol_ENABLED
59 : !%%%%%%%%%%%%%%%%%
60 :
61 : integer(IK) :: info
62 : #if ULD_ENABLED
63 : type(uppDia_type), parameter :: subset = uppDia
64 : #elif !(UXD_ENABLED || XLD_ENABLED)
65 : #error "Unrecognized interface."
66 : #endif
67 123391 : chol = ZERO
68 : optionalBlock: block
69 5141 : if (present(operation)) then
70 3200 : if (same_type_as(operation, transHerm)) then
71 39360 : call setMatChol(mat, subset, info, chol, transHerm)
72 1600 : exit optionalBlock
73 1600 : elseif (.not. same_type_as(operation, nothing)) then
74 : error stop MODULE_NAME//SK_"@getMatChol(): Unrecognized `operation` other than `nothing` or `transHerm` specified." ! LCOV_EXCL_LINE
75 : end if
76 : end if
77 : ! default operation.
78 84031 : call setMatChol(mat, subset, info, chol, nothing)
79 : end block optionalBlock
80 5141 : if (info /= 0_IK) error stop MODULE_NAME//SK_"@getMatChol(): Cholesky factorization failed."
81 :
82 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
83 : #elif setMatChol_ENABLED && AXX_ENABLED && ONO_ENABLED
84 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
85 :
86 1600 : call setMatChol(mat, subset, info, iteration)
87 :
88 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89 : #elif setMatChol_ENABLED && ANI_ENABLED
90 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 :
92 : real(TKC) :: summ
93 : integer(IK) :: irow
94 : #if IMP_ENABLED
95 : integer(IK) :: ndim
96 1356301 : ndim = size(mat, 1, IK)
97 4068903 : CHECK_ASSERTION(__LINE__, size(mat, 1, IK) == size(mat, 2, IK), SK_"@setMatChol(): The condition `size(mat, 1) == size(mat, 2)` must hold. ndim = "//getStr([size(mat, 1, IK), size(mat, 2, IK)]))
98 14919311 : CHECK_ASSERTION(__LINE__, all(shape(mat, IK) == shape(chol, IK)), SK_"@setMatChol(): The condition `all(shape(mat) == shape(chol))` must hold. shape(mat), shape(chol) = "//getStr([shape(mat, IK), shape(chol, IK)]))
99 : #elif EXP_ENABLED
100 : ! Ensure offsets and subset bounds are logical.
101 41 : CHECK_ASSERTION(__LINE__, 0_IK <= ndim, SK_"@setMatChol(): The condition `0 <= ndim` must hold. ndim = "//getStr(ndim))
102 205 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(mat, kind = IK)), SK_"@setMatChol(): The condition `all(0 <= [roff, coff])` must hold. roff, coff = "//getStr(1_IK - lbound(mat)))
103 492 : CHECK_ASSERTION(__LINE__, all(ndim <= ubound(mat, kind = IK)), SK_"@setMatChol(): The condition `ndim + [roff, coff] <= shape(mat)` must hold. ndim, roff, coff, shape(mat) = "//getStr([ndim, 1 - lbound(mat, kind = IK), shape(mat, kind = IK)]))
104 492 : CHECK_ASSERTION(__LINE__, all(ndim <= ubound(chol, kind = IK)), SK_"@setMatChol(): The condition `ndim + [roffc, coffc] <= shape(chol)` must hold. ndim, roffc, coffc, shape(chol) = "//getStr([ndim, 1 - lbound(chol, kind = IK), shape(chol, kind = IK)]))
105 : #else
106 : #error "Unrecognized interface."
107 : #endif
108 : ! Define operations based on the specified triangular side.
109 : #if UXD_ENABLED
110 : #define GET_MATMUL(I,J) matmul(J,I)
111 : #define GETMAT(MAT,I,J) MAT(J,I)
112 : #elif XLD_ENABLED
113 : #define GET_MATMUL(I,J) matmul(I,J)
114 : #define GETMAT(MAT,I,J) MAT(I,J)
115 : #else
116 : #error "Unrecognized interface."
117 : #endif
118 4214052 : do info = 1_IK, ndim
119 : #if ONO_ENABLED
120 3822111 : summ = real(mat(info, info), TKC) - real(dot_product(GETMAT(chol, info, 1 : info - 1), GETMAT(chol, info, 1 : info - 1)), TKC)
121 : #elif OTH_ENABLED
122 850163 : summ = real(mat(info, info), TKC) - real(dot_product(GETMAT(chol, 1 : info - 1, info), GETMAT(chol, 1 : info - 1, info)), TKC)
123 : #else
124 : #error "Unrecognized interface."
125 : #endif
126 2861323 : if (0._TKC < summ) then
127 2857710 : summ = sqrt(summ)
128 2857710 : chol(info, info) = summ
129 2857710 : summ = 1._TKC / summ
130 : #if ONO_ENABLED
131 : !GETMAT(chol, info + 1 : ndim, info) = (GETMAT(mat, info + 1 : ndim, info) - GET_MATMUL(GET_CONJG(GETMAT(chol, info + 1 : ndim, 1 : info - 1)),GETMAT(chol, info, 1 : info - 1))) * summ
132 3823215 : do irow = info + 1, ndim
133 3927454 : GETMAT(chol, irow, info) = (GETMAT(mat, irow, info) - dot_product(GETMAT(chol, info, 1 : info - 1), GETMAT(chol, irow, 1 : info - 1))) * summ
134 : end do
135 : #elif OTH_ENABLED
136 851645 : do irow = info + 1, ndim
137 1217383 : GETMAT(chol, info, irow) = (GET_CONJG(GETMAT(mat, irow, info) - dot_product(GETMAT(chol, 1 : info - 1, irow), GETMAT(chol, 1 : info - 1, info)))) * summ
138 : end do
139 : #endif
140 : cycle
141 : end if
142 1352729 : return
143 : end do
144 1352729 : info = 0_IK
145 : #undef GET_MATMUL
146 : #undef GETMAT
147 :
148 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
149 : #elif setMatChol_ENABLED && (ABI_ENABLED || ABR_ENABLED) && ONO_ENABLED
150 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
151 :
152 : #if CK_ENABLED
153 : #define OPERATION transHerm
154 : #elif RK_ENABLED
155 : #define OPERATION transSymm
156 : #else
157 : #error "Unrecognized interface."
158 : #endif
159 : #if ABI_ENABLED
160 : ! Set the block size for this environment.
161 : ! Per LAPACK documentation, the current value (ILAENV( 1, 'dpotrf', uplo, ndim, -1, -1, -1 )) may not be optimal.
162 : ! Further runtime experimentation might be needed.
163 : integer(IK) , parameter :: NBLOCK = 64_IK
164 : integer(IK) :: icol, jcol, bdim_def
165 : #elif ABR_ENABLED
166 : integer(IK) :: ndimHalf1, ndimHalf2
167 : #endif
168 : ! Set/check bounds.
169 : #if IMP_ENABLED
170 : integer(IK) , parameter :: roff = 0, coff = 0
171 : integer(IK) :: ndim
172 11200 : ndim = size(mat, 1, IK)
173 33600 : CHECK_ASSERTION(__LINE__, size(mat, 1, IK) == size(mat, 2, IK), SK_"@setMatChol(): The condition `size(mat, 1) == size(mat, 2)` must hold. ndim = "//getStr([size(mat, 1, IK), size(mat, 2, IK)]))
174 : #elif EXP_ENABLED
175 : ! Ensure offsets and subset bounds are logical.
176 60910 : CHECK_ASSERTION(__LINE__, 0_IK <= ndim, SK_"@setMatChol(): The condition `0 <= ndim` must hold. ndim = "//getStr(ndim))
177 304550 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(mat, kind = IK)), SK_"@setMatChol(): The condition `all(0 <= [roff, coff])` must hold. roff, coff = "//getStr(1_IK - lbound(mat)))
178 730920 : CHECK_ASSERTION(__LINE__, all(ndim <= ubound(mat, kind = IK)), SK_"@setMatChol(): The condition `ndim + [roff, coff] <= shape(mat)` must hold. ndim, roff, coff, shape(mat) = "//getStr([ndim, 1 - lbound(mat, kind = IK), shape(mat, kind = IK)]))
179 : #else
180 : #error "Unrecognized interface."
181 : #endif
182 72110 : info = 0_IK
183 90172 : if (ndim == 0_IK) return ! Quick return if possible.
184 :
185 : !%%%%%%%%%%
186 : #if ABI_ENABLED
187 : !%%%%%%%%%%
188 :
189 8005 : if (present(bdim)) then
190 3203 : CHECK_ASSERTION(__LINE__, 0_IK < bdim, SK_"@setMatChol(): The condition `0 < bdim` must hold. bdim = "//getStr(bdim))
191 : bdim_def = bdim
192 : else
193 : bdim_def = NBLOCK
194 : end if
195 8005 : if (bdim_def <= 1_IK .or. ndim <= bdim_def) then ! Use unblocked code.
196 7388 : call setMatChol(mat, subset, info, recursion, ndim, roff, coff)
197 : else ! Use blocked code.
198 : ! Compute the Cholesky factorization mat = U**T*U.
199 1505 : do icol = 0_IK, ndim - 1_IK, bdim_def
200 : ! Update and factorize the current diagonal block and test for non-positive-definiteness.
201 1184 : jcol = min(bdim_def, ndim - icol)
202 : #if UXD_ENABLED
203 : call setMatUpdate ( mat = mat & ! LCOV_EXCL_LINE
204 : , class = hermitian & ! LCOV_EXCL_LINE
205 : , subset = uppDia & ! LCOV_EXCL_LINE
206 : , matA = mat & ! LCOV_EXCL_LINE
207 : , operationA = trans & ! LCOV_EXCL_LINE
208 : , alpha = -ONE & ! LCOV_EXCL_LINE
209 : , beta = ONE & ! LCOV_EXCL_LINE
210 : , ndim = jcol & ! LCOV_EXCL_LINE
211 : , ndum = icol & ! LCOV_EXCL_LINE
212 : , roff = roff + icol & ! LCOV_EXCL_LINE
213 : , coff = coff + icol & ! LCOV_EXCL_LINE
214 : , roffA = roff & ! LCOV_EXCL_LINE
215 : , coffA = coff + icol & ! LCOV_EXCL_LINE
216 593 : ) ! Update and factor MatPosDef22.
217 593 : call setMatChol(mat, subset, info, recursion, jcol, roff + icol, coff + icol)
218 593 : if (info /= 0_IK) then
219 148 : info = info + icol
220 148 : return
221 : end if
222 606 : if (icol + jcol <= ndim) then ! Compute the current block row.
223 : call setMatMulAdd ( matA = mat & ! LCOV_EXCL_LINE
224 : , operationA = transHerm & ! LCOV_EXCL_LINE
225 : , matB = mat & ! LCOV_EXCL_LINE
226 : , matC = mat & ! LCOV_EXCL_LINE
227 : , alpha = -ONE & ! LCOV_EXCL_LINE
228 : , beta = ONE & ! LCOV_EXCL_LINE
229 : , nrow = jcol & ! LCOV_EXCL_LINE
230 : , ncol = ndim - icol - jcol & ! LCOV_EXCL_LINE
231 : , ndum = icol & ! LCOV_EXCL_LINE
232 : , roffA = roff & ! LCOV_EXCL_LINE
233 : , coffA = coff + icol & ! LCOV_EXCL_LINE
234 : , roffB = roff & ! LCOV_EXCL_LINE
235 : , coffB = coff + icol + jcol & ! LCOV_EXCL_LINE
236 : , roffC = roff + icol & ! LCOV_EXCL_LINE
237 : , coffC = coff + icol + jcol & ! LCOV_EXCL_LINE
238 445 : )
239 : ! syslin = AXB
240 : call setMatMulTri ( matA = mat & ! LCOV_EXCL_LINE
241 : , classA = upperDiag & ! LCOV_EXCL_LINE
242 : , operationA = transUnit & ! LCOV_EXCL_LINE
243 : , matB = mat & ! LCOV_EXCL_LINE
244 : , alpha = ONE & ! LCOV_EXCL_LINE
245 : , nrow = jcol & ! LCOV_EXCL_LINE
246 : , ncol = ndim - icol - jcol & ! LCOV_EXCL_LINE
247 : , roffA = roff + icol & ! LCOV_EXCL_LINE
248 : , coffA = coff + icol & ! LCOV_EXCL_LINE
249 : , roffB = roff + icol & ! LCOV_EXCL_LINE
250 : , coffB = coff + icol + jcol & ! LCOV_EXCL_LINE
251 445 : )
252 : end if
253 : #elif XLD_ENABLED
254 : ! Compute the Cholesky factorization mat = l*l**t.
255 : call setMatUpdate ( mat = mat & ! LCOV_EXCL_LINE
256 : , class = hermitian & ! LCOV_EXCL_LINE
257 : , subset = lowDia & ! LCOV_EXCL_LINE
258 : , matA = mat & ! LCOV_EXCL_LINE
259 : , operationA = nothing & ! LCOV_EXCL_LINE
260 : , alpha = -ONE & ! LCOV_EXCL_LINE
261 : , beta = ONE & ! LCOV_EXCL_LINE
262 : , ndim = jcol & ! LCOV_EXCL_LINE
263 : , ndum = icol & ! LCOV_EXCL_LINE
264 : , roff = roff + icol & ! LCOV_EXCL_LINE
265 : , coff = coff + icol & ! LCOV_EXCL_LINE
266 : , roffA = roff + icol & ! LCOV_EXCL_LINE
267 : , coffA = coff & ! LCOV_EXCL_LINE
268 591 : ) ! Update and factor MatPosDef22.
269 591 : call setMatChol(mat, subset, info, recursion, jcol, roff + icol, coff + icol)
270 591 : if (info /= 0_IK) then
271 148 : info = info + icol
272 148 : return
273 : end if
274 603 : if (icol + jcol <= ndim) then ! Compute the current block column.
275 : call setMatMulAdd ( matA = mat & ! LCOV_EXCL_LINE
276 : , matB = mat & ! LCOV_EXCL_LINE
277 : , matC = mat & ! LCOV_EXCL_LINE
278 : , alpha = -ONE & ! LCOV_EXCL_LINE
279 : , beta = ONE & ! LCOV_EXCL_LINE
280 : , nrow = ndim - icol - jcol & ! LCOV_EXCL_LINE
281 : , ncol = jcol & ! LCOV_EXCL_LINE
282 : , ndum = icol & ! LCOV_EXCL_LINE
283 : , roffA = roff + icol + jcol & ! LCOV_EXCL_LINE
284 : , coffA = coff & ! LCOV_EXCL_LINE
285 : , roffB = roff + icol & ! LCOV_EXCL_LINE
286 : , coffB = coff & ! LCOV_EXCL_LINE
287 : , roffC = roff + icol + jcol & ! LCOV_EXCL_LINE
288 : , coffC = coff + icol & ! LCOV_EXCL_LINE
289 : , operationB = transHerm & ! LCOV_EXCL_LINE
290 443 : )
291 : ! syslin = XAB
292 : call setMatMulTri ( matA = mat & ! LCOV_EXCL_LINE
293 : , matB = mat & ! LCOV_EXCL_LINE
294 : , classB = lowerDiag & ! LCOV_EXCL_LINE
295 : , operationB = transUnit & ! LCOV_EXCL_LINE
296 : , alpha = ONE & ! LCOV_EXCL_LINE
297 : , nrow = ndim - icol - jcol & ! LCOV_EXCL_LINE
298 : , ncol = jcol & ! LCOV_EXCL_LINE
299 : , roffA = roff + icol + jcol & ! LCOV_EXCL_LINE
300 : , coffA = coff + icol & ! LCOV_EXCL_LINE
301 : , roffB = roff + icol & ! LCOV_EXCL_LINE
302 : , coffB = coff + icol & ! LCOV_EXCL_LINE
303 443 : )
304 : end if
305 : #else
306 : #error "Unrecognized interface."
307 : #endif
308 : end do
309 : end if
310 :
311 : !%%%%%%%%%%
312 : #elif ABR_ENABLED
313 : !%%%%%%%%%%
314 :
315 64105 : if(1_IK < ndim) then ! use recursive code
316 28789 : ndimHalf1 = ndim / 2_IK
317 28789 : ndimHalf2 = ndim - ndimHalf1
318 28789 : call setMatChol(mat, subset, info, control, ndimHalf1, roff, coff) ! Factor MatPosDef11.
319 28789 : if (info /= 0_IK) return
320 : #if UXD_ENABLED
321 : ! Compute the Cholesky factorization mat = U**T*U :: syslin = AXB
322 : call setMatMulTri ( matA = mat & ! LCOV_EXCL_LINE
323 : , classA = upperDiag & ! LCOV_EXCL_LINE
324 : , operationA = transUnit & ! LCOV_EXCL_LINE
325 : , matB = mat & ! LCOV_EXCL_LINE
326 : , alpha = ONE & ! LCOV_EXCL_LINE
327 : , nrow = ndimHalf1 & ! LCOV_EXCL_LINE
328 : , ncol = ndimHalf2 & ! LCOV_EXCL_LINE
329 : , roffA = roff & ! LCOV_EXCL_LINE
330 : , coffA = coff & ! LCOV_EXCL_LINE
331 : , roffB = roff & ! LCOV_EXCL_LINE
332 : , coffB = coff + ndimHalf1 & ! LCOV_EXCL_LINE
333 11770 : ) ! Update and scale MatPosDef12.
334 : !block
335 : !use pm_blas, only: blasTRSM
336 : !CALL blasTRSM('L', 'U', 'C', 'N', ndimHalf1, ndimHalf2, ONE, mat(1, 1), size(mat, 1), mat(1, ndimHalf1 + 1), size(mat, 1))
337 : !end block
338 : call setMatUpdate ( mat = mat & ! LCOV_EXCL_LINE
339 : , class = hermitian & ! LCOV_EXCL_LINE
340 : , subset = uppDia & ! LCOV_EXCL_LINE
341 : , matA = mat & ! LCOV_EXCL_LINE
342 : , operationA = trans & ! LCOV_EXCL_LINE
343 : , alpha = -ONE & ! LCOV_EXCL_LINE
344 : , beta = ONE & ! LCOV_EXCL_LINE
345 : , ndim = ndimHalf2 & ! LCOV_EXCL_LINE
346 : , ndum = ndimHalf1 & ! LCOV_EXCL_LINE
347 : , roff = roff + ndimHalf1 & ! LCOV_EXCL_LINE
348 : , coff = coff + ndimHalf1 & ! LCOV_EXCL_LINE
349 : , roffA = roff & ! LCOV_EXCL_LINE
350 : , coffA = coff + ndimHalf1 & ! LCOV_EXCL_LINE
351 11770 : ) ! Update and factor MatPosDef22.
352 : #elif XLD_ENABLED
353 : ! Compute the Cholesky factorization mat = L*L**T :: syslin = XAB
354 : call setMatMulTri ( matA = mat & ! LCOV_EXCL_LINE
355 : , matB = mat & ! LCOV_EXCL_LINE
356 : , classB = lowerDiag & ! LCOV_EXCL_LINE
357 : , operationB = transUnit & ! LCOV_EXCL_LINE
358 : , alpha = ONE & ! LCOV_EXCL_LINE
359 : , nrow = ndimHalf2 & ! LCOV_EXCL_LINE
360 : , ncol = ndimHalf1 & ! LCOV_EXCL_LINE
361 : , roffA = roff + ndimHalf1 & ! LCOV_EXCL_LINE
362 : , coffA = coff & ! LCOV_EXCL_LINE
363 : , roffB = roff & ! LCOV_EXCL_LINE
364 : , coffB = coff & ! LCOV_EXCL_LINE
365 11771 : ) ! Update and scale MatPosDef21.
366 : call setMatUpdate ( mat = mat & ! LCOV_EXCL_LINE
367 : , class = hermitian & ! LCOV_EXCL_LINE
368 : , subset = lowDia & ! LCOV_EXCL_LINE
369 : , matA = mat & ! LCOV_EXCL_LINE
370 : , operationA = nothing & ! LCOV_EXCL_LINE
371 : , alpha = -ONE & ! LCOV_EXCL_LINE
372 : , beta = ONE & ! LCOV_EXCL_LINE
373 : , ndim = ndimHalf2 & ! LCOV_EXCL_LINE
374 : , ndum = ndimHalf1 & ! LCOV_EXCL_LINE
375 : , roff = roff + ndimHalf1 & ! LCOV_EXCL_LINE
376 : , coff = coff + ndimHalf1 & ! LCOV_EXCL_LINE
377 : , roffA = roff + ndimHalf1 & ! LCOV_EXCL_LINE
378 : , coffA = coff & ! LCOV_EXCL_LINE
379 11771 : ) ! Update and factor MatPosDef22.
380 : #else
381 : #error "Unrecognized interface."
382 : #endif
383 23541 : call setMatChol(mat, subset, info, control, ndimHalf2, roff + ndimHalf1, coff + ndimHalf1)
384 23541 : if (info /= 0_IK) then
385 6118 : info = info + ndimHalf1
386 6118 : return
387 : end if
388 : else
389 : ! Test for non-positive-definiteness.
390 35316 : if (GET_RE(ZERO) < GET_RE(mat(1, 1))) then
391 28916 : mat(1, 1) = sqrt(GET_RE(mat(1, 1)))
392 : else
393 : !write(*,*) GET_RE(mat(1, 1))
394 6400 : info = 1_IK
395 6400 : return
396 : end if
397 : end if
398 : #else
399 : !%%%%%%%%%%%%%%%%%%%%%%%%
400 : #error "Unrecognized interface."
401 : !%%%%%%%%%%%%%%%%%%%%%%%%
402 : #endif
403 : #else
404 : !%%%%%%%%%%%%%%%%%%%%%%%%
405 : #error "Unrecognized interface."
406 : !%%%%%%%%%%%%%%%%%%%%%%%%
407 : #endif
408 : #undef TYPE_KIND
409 : #undef GET_CONJG
410 : #undef OPERATION
411 : #undef GETMAT
412 : #undef GET_RE
413 : #undef CHOMAT
|
