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 interface [pm_matrixDet](@ref pm_matrixDet).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Apr 21, 2017, 1:54 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : #if CK_ENABLED
28 : #define TYPE_KIND complex(TKC)
29 : #define GET_CONJG(X) conjg(X)
30 : #elif RK_ENABLED
31 : #define GET_CONJG(X) X
32 : #define TYPE_KIND real(TKC)
33 : #else
34 : #error "Unrecognized interface."
35 : #endif
36 : ! Define the output variable.
37 : #if getMatDetSqrt_ENABLED || setMatDetSqrt_ENABLED
38 : #define GETLOG(X) X
39 : #define OUTPUT detSqrt
40 : #define INCREMENT(X,Y) X = X * Y
41 : #elif getMatDetSqrtLog_ENABLED || setMatDetSqrtLog_ENABLED
42 : #define GETLOG(X) log(X)
43 : #define OUTPUT detSqrtLog
44 : #define INCREMENT(X,Y) X = X + log(Y)
45 : #elif !(getMatDet_ENABLED || setMatDet_ENABLED)
46 : #error "Unrecognized interface."
47 : #endif
48 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49 : #if getMatDet_ENABLED || setMatDet_ENABLED
50 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51 :
52 : #if getMatDet_ENABLED
53 : integer(IK) :: info
54 260 : TYPE_KIND :: matLUP(size(mat, 1, IK), size(mat, 2, IK))
55 : #elif setMatDet_ENABLED
56 : #define matLUP mat
57 : #else
58 : #error "Unrecognized interface."
59 : #endif
60 : integer(IK) :: idim, ndim
61 20 : integer(IK) :: rperm(size(mat,1))
62 : ndim = size(mat, 1, IK)
63 150 : CHECK_ASSERTION(__LINE__, size(mat, 2, IK) == ndim, SK_"The input `mat` must be square. shape(mat) = "//getStr(shape(mat)))
64 : #if getMatDet_ENABLED
65 4046 : matLUP = mat
66 : #endif
67 150 : call setMatLUP(matLUP, rperm, info) ! This returns det as +-1.
68 150 : if (info /= 0_IK) then
69 : #if getMatDet_ENABLED
70 : error stop "LU factorization failed." ! LCOV_EXCL_LINE
71 : #elif setMatDet_ENABLED
72 : return
73 : #endif
74 : end if
75 95 : det = 1 ! parity
76 834 : do idim = 1_IK, ndim
77 684 : det = det * matLUP(idim, idim)
78 834 : if (rperm(idim) /= idim) det = -det
79 : end do
80 :
81 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 : #elif getMatDetSqrt_ENABLED || getMatDetSqrtLog_ENABLED
83 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84 :
85 : #if getMatDetSqrt_ENABLED
86 : #define SETDET(SUBSET) \
87 : call setMatDetSqrt(mat, SUBSET, detSqrt, info, chol, nothing)
88 : #elif getMatDetSqrtLog_ENABLED
89 : #define SETDET(SUBSET) \
90 : call setMatDetSqrtLog(mat, SUBSET, detSqrtLog, info, chol, nothing)
91 : #else
92 : #error "Unrecognized interface."
93 : #endif
94 : integer(IK) :: info
95 500204 : TYPE_KIND :: chol(size(mat, 1, IK), size(mat, 2, IK))
96 250102 : if (size(mat, 1, IK) < 64_IK) then
97 : ! Use the unblocked algorithm from this module.
98 250102 : if (present(subset)) then
99 : select type (subset)
100 : type is (uppDia_type)
101 250027 : SETDET(uppDia)
102 : type is (lowDia_type)
103 24 : SETDET(lowDia)
104 : class default
105 : error stop "Unrecognized `subset`. Only objects of `uppDia_type()` and `lowDia_type()` are recognized." ! LCOV_EXCL_LINE
106 : end select
107 : else
108 51 : SETDET(uppDia)
109 : end if
110 250102 : if (info /= 0_IK) error stop "The Cholesky factorization of small matrix failed."
111 : else
112 : ! Use the blocked algorithm potentially dispatched to LAPACK if available.
113 0 : if (present(subset)) then
114 : select type (subset)
115 : type is (uppDia_type)
116 0 : call setMatCopy(chol, rdpack, mat, rdpack, uppDia)
117 0 : call setMatChol(chol, uppDia, info, iteration)
118 : type is (lowDia_type)
119 0 : call setMatCopy(chol, rdpack, mat, rdpack, lowDia)
120 0 : call setMatChol(chol, lowDia, info, iteration)
121 : class default
122 : error stop "Unrecognized `subset`. Only objects of `uppDia_type()` and `lowDia_type()` are recognized." ! LCOV_EXCL_LINE
123 : end select
124 : else
125 0 : call setMatCopy(chol, rdpack, mat, rdpack, uppDia)
126 0 : call setMatChol(chol, uppDia, info, iteration)
127 : end if
128 0 : if (info /= 0_IK) error stop "The Cholesky factorization of large matrix failed."
129 0 : OUTPUT = GETLOG(real(chol(1, 1), TKC))
130 0 : do info = 2, size(chol, 1, IK)
131 0 : INCREMENT(OUTPUT,real(chol(info, info), TKC))
132 : end do
133 : info = 0_IK
134 : end if
135 : #undef SETDET
136 :
137 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138 : #elif setMatDetSqrt_ENABLED || setMatDetSqrtLog_ENABLED
139 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 :
141 : real(TKC) :: summ
142 : integer(IK) :: irow
143 : !#if IMP_ENABLED
144 : integer(IK) :: ndim
145 250153 : ndim = size(mat, 1, IK)
146 750459 : 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)]))
147 2751683 : 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)]))
148 : !#elif EXP_ENABLED
149 : ! ! Ensure offsets and subset bounds are logical.
150 : ! check_assertion(__LINE__, 0_IK <= ndim, SK_"@setMatChol(): The condition `0 <= ndim` must hold. ndim = "//getStr(ndim))
151 : ! 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)))
152 : ! 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)]))
153 : ! 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)]))
154 : !#else
155 : !#error "Unrecognized interface."
156 : !#endif
157 3001836 : CHECK_ASSERTION(__LINE__, all(ndim <= ubound(chol, kind = IK)), SK_"@setMatChol(): The condition `all(shape(mat) == shape(chol))` must hold. shape(mat), ubound(chol) = "//getStr([shape(mat, IK), ubound(chol, kind = IK)]))
158 : ! Define operations based on the specified triangular side.
159 : #if UXD_ENABLED
160 : #define GET_MATMUL(I,J) matmul(J,I)
161 : #define GETMAT(MAT,I,J) MAT(J,I)
162 : #elif XLD_ENABLED
163 : #define GET_MATMUL(I,J) matmul(I,J)
164 : #define GETMAT(MAT,I,J) MAT(I,J)
165 : #else
166 : #error "Unrecognized interface."
167 : #endif
168 250153 : if (1_IK < ndim) then
169 250118 : info = 1_IK
170 250118 : summ = real(mat(info, info), TKC)
171 250118 : if (0._TKC < summ) then
172 250118 : summ = sqrt(summ)
173 250118 : chol(info, info) = summ
174 250118 : OUTPUT = GETLOG(summ)
175 250118 : summ = 1._TKC / summ
176 : #if ONO_ENABLED
177 500284 : do irow = info + 1, ndim
178 500284 : GETMAT(chol, irow, info) = (GETMAT(mat, irow, info) - dot_product(GETMAT(chol, info, 1 : info - 1), GETMAT(chol, irow, 1 : info - 1))) * summ
179 : end do
180 : #elif OTH_ENABLED
181 122 : do irow = info + 1, ndim
182 122 : 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
183 : end do
184 : #endif
185 500406 : do info = 2_IK, ndim
186 : #if ONO_ENABLED
187 500628 : summ = real(mat(info, info), TKC) - real(dot_product(GETMAT(chol, info, 1 : info - 1), GETMAT(chol, info, 1 : info - 1)), TKC)
188 : #elif OTH_ENABLED
189 222 : summ = real(mat(info, info), TKC) - real(dot_product(GETMAT(chol, 1 : info - 1, info), GETMAT(chol, 1 : info - 1, info)), TKC)
190 : #else
191 : #error "Unrecognized interface."
192 : #endif
193 250288 : if (0._TKC < summ) then
194 250288 : summ = sqrt(summ)
195 250288 : chol(info, info) = summ
196 250288 : INCREMENT(OUTPUT,summ)
197 250288 : summ = 1._TKC / summ
198 : #if ONO_ENABLED
199 : !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
200 250422 : do irow = info + 1, ndim
201 250750 : GETMAT(chol, irow, info) = (GETMAT(mat, irow, info) - dot_product(GETMAT(chol, info, 1 : info - 1), GETMAT(chol, irow, 1 : info - 1))) * summ
202 : end do
203 : #elif OTH_ENABLED
204 140 : do irow = info + 1, ndim
205 218 : 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
206 : end do
207 : #endif
208 : cycle
209 : end if
210 250118 : return
211 : end do
212 250118 : info = 0_IK
213 : end if
214 35 : elseif (1_IK == ndim) then
215 35 : summ = real(mat(1, 1), TKC)
216 35 : if (0._TKC < real(summ, TKC)) then
217 35 : summ = sqrt(summ)
218 35 : OUTPUT = GETLOG(summ)
219 35 : chol(1,1) = summ
220 35 : info = 0_IK
221 : else
222 0 : info = 1_IK
223 : end if
224 : else
225 0 : info = 0_IK
226 : end if
227 : #undef GET_MATMUL
228 : #undef GETMAT
229 :
230 : #else
231 : !%%%%%%%%%%%%%%%%%%%%%%%%
232 : #error "Unrecognized interface."
233 : !%%%%%%%%%%%%%%%%%%%%%%%%
234 : #endif
235 : #undef INCREMENT
236 : #undef TYPE_KIND
237 : #undef GET_CONJG
238 : #undef OUTPUT
239 : #undef GETLOG
240 : #undef matLUP
|
