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 file contains procedure implementations of [pm_matrixInv](@ref pm_matrixInv).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Set the type and kind.
28 : #if CK_ENABLED
29 : #define GET_CONJG(X) conjg(X)
30 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
31 : #define TYPE_KIND complex(TKC)
32 : #elif RK_ENABLED
33 : #define GET_CONJG(X) X
34 : real(TKC), parameter :: ZERO = 0._TKC, ONE = 1._TKC
35 : #define TYPE_KIND real(TKC)
36 : #else
37 : #error "Unrecognized interface."
38 : #endif
39 : ! Define the runtime checks.
40 : #define CHECK_SHAPE_MAT \
41 : CHECK_ASSERTION(__LINE__, size(mat, 1, IK) == size(mat, 2, IK), SK_"@setMatInv(): The condition `size(mat, 1) == size(mat, 2)` must hold. shape(mat) = "//getStr(shape(mat, IK)))
42 : #define CHECK_SHAPE_INV \
43 : CHECK_ASSERTION(__LINE__, all(shape(inv, IK) == shape(mat, IK)), SK_"@setMatInv(): The condition `all(shape(inv) == shape(mat))` must hold. shape(inv), shape(mat) = "//getStr([shape(inv, IK), shape(mat, IK)]))
44 :
45 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
46 : #if getMatInv_ENABLED && (Def_ENABLED || Det_ENABLED || Inf_ENABLED)
47 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 :
49 316 : integer(IK) :: rperm(size(mat, 1, IK))
50 160 : TYPE_KIND :: lup(size(mat, 1, IK), size(mat, 2, IK))
51 : #if Def_ENABLED || Det_ENABLED
52 : #define AUXIL info
53 : integer(IK) :: info
54 : #elif !Inf_ENABLED
55 : #error "Unrecognized interface."
56 : #endif
57 2490 : lup = mat
58 158 : call setMatLUP(lup, rperm, AUXIL)
59 158 : if (AUXIL /= 0_IK) then
60 : #if Inf_ENABLED
61 : return
62 : #else
63 0 : error stop MODULE_NAME//SK_"@getMatInv(): LUP factorization of the input `mat` failed."
64 : #endif
65 : end if
66 158 : call setMatInv(inv, lup, rperm)
67 : #if Det_ENABLED
68 : block
69 : integer(IK) :: idim
70 2 : auxil = ONE
71 14 : do idim = 1, size(rperm, 1, IK)
72 12 : auxil = auxil * lup(idim, idim)
73 14 : if (rperm(idim) /= idim) auxil = -auxil
74 : end do
75 : end block
76 : #endif
77 : #undef AUXIL
78 :
79 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80 : #elif getMatInv_ENABLED && (CUD_ENABLED || CUU_ENABLED || CLD_ENABLED || CLU_ENABLED || LUP_ENABLED || CCU_ENABLED || CCL_ENABLED)
81 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82 :
83 133 : call setMatInv(inv, mat, auxil)
84 : !#if CCU_ENABLED
85 : ! call setMatCopy(inv(1:,2:), rdpack, inv(2:,1:), rdpack, lowDia, transHerm)
86 : !#elif CCL_ENABLED
87 : ! call setMatCopy(inv(2:,1:), rdpack, inv(1:,2:), rdpack, uppDia, transHerm)
88 : !#endif
89 :
90 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91 : #elif setMatInv_ENABLED && (CUD_ENABLED || CLD_ENABLED || CUU_ENABLED || CLU_ENABLED)
92 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
93 :
94 : ! Set the indexing rules based on the subset type.
95 : #if CUD_ENABLED || CUU_ENABLED
96 : #define GET_INDEX(I,J) J,I
97 : #elif CLD_ENABLED || CLU_ENABLED
98 : #define GET_INDEX(I,J) I,J
99 : #else
100 : #error "Unrecognized interface."
101 : #endif
102 : TYPE_KIND :: summ
103 : integer(IK) :: irow, icol, idim, ndim
104 180 : ndim = size(mat, 1, IK)
105 180 : CHECK_SHAPE_MAT
106 1980 : CHECK_SHAPE_INV
107 835 : do irow = 1_IK, ndim
108 : #if CLD_ENABLED || CLU_ENABLED
109 809 : inv(1 : irow - 1, irow) = ZERO
110 : #endif
111 : #if CUD_ENABLED || CLD_ENABLED
112 358 : CHECK_ASSERTION(__LINE__, mat(irow, irow) /= ZERO, SK_"@setMatInv(): The condition `mat(irow, irow) /= 0` must hold. irow = "//getStr(irow))
113 163 : inv(irow, irow) = ONE / mat(irow, irow)
114 : #elif CUU_ENABLED || CLU_ENABLED
115 297 : CHECK_ASSERTION(__LINE__, mat(irow, irow) == ONE, SK_"@setMatInv(): The condition `mat(irow, irow) == 1` must hold. irow = "//getStr(irow))
116 137 : inv(irow, irow) = ONE
117 : #endif
118 : #if CUD_ENABLED || CUU_ENABLED
119 918 : inv(irow + 1 : ndim, irow) = ZERO
120 : #endif
121 1907 : do icol = 1_IK, irow - 1_IK
122 0 : summ = ZERO
123 3153 : do idim = icol, irow - 1_IK
124 3153 : summ = summ + mat(GET_INDEX(irow, idim)) * inv(GET_INDEX(idim, icol))
125 : end do
126 : #if CUD_ENABLED || CLD_ENABLED
127 935 : inv(GET_INDEX(irow, icol)) = -summ * inv(irow, irow)
128 : #elif CUU_ENABLED || CLU_ENABLED
129 792 : inv(GET_INDEX(irow, icol)) = -summ
130 : #endif
131 : end do
132 : end do
133 : #undef GET_INDEX
134 :
135 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 : #elif setMatInv_ENABLED && (CCU_ENABLED || CCL_ENABLED)
137 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138 :
139 : ! Set the indexing rules based on the subset type.
140 : #if CCL_ENABLED && (FUL_ENABLED || UXD_ENABLED)
141 : #define INDMAT(I,J)I,J
142 : #define INDINV(I,J)I,J
143 : #elif CCU_ENABLED && (FUL_ENABLED || XLD_ENABLED)
144 : #define INDMAT(I,J)J,I
145 : #define INDINV(I,J)J,I
146 : #elif CCL_ENABLED && XLD_ENABLED
147 : #define INDMAT(I,J)I,J
148 : #define INDINV(I,J)J,I
149 : #elif CCU_ENABLED && UXD_ENABLED
150 : #define INDMAT(I,J)J,I
151 : #define INDINV(I,J)I,J
152 : #else
153 : #error "Unrecognized interface."
154 : #endif
155 : !TYPE_KIND :: summ
156 : integer(IK) :: irow, icol, ndim !, idim
157 107 : ndim = size(mat, 1, IK)
158 107 : CHECK_SHAPE_MAT
159 1177 : CHECK_SHAPE_INV
160 107 : if (1_IK < ndim) then
161 461 : CHECK_ASSERTION(__LINE__, all(real(getMatCopy(lfpack, mat, rdpack, dia), TKC) /= 0._TKC), \
162 : SK_"@setMatInv(): The The diagonal elements of the Cholesky factorization must be non-zero. getMatCopy(lfpack, mat, rdpack, dia) = "//getStr(getMatCopy(lfpack, mat, rdpack, dia)))
163 373 : do irow = 1_IK, ndim - 1_IK
164 285 : inv(irow, irow) = ONE / mat(irow, irow)
165 1062 : do icol = irow + 1_IK, ndim
166 : #if FUL_ENABLED || (CCL_ENABLED && UXD_ENABLED) || (CCU_ENABLED && XLD_ENABLED)
167 1568 : inv(INDINV(irow, icol)) = -dot_product(mat(INDMAT(icol, irow : icol - 1)), inv(INDINV(irow, irow : icol - 1))) / mat(icol, icol)
168 : #elif (CCU_ENABLED && UXD_ENABLED) || (CCL_ENABLED && XLD_ENABLED)
169 815 : inv(INDINV(irow, icol)) = -GET_CONJG(dot_product(mat(INDMAT(icol, irow : icol - 1)), GET_CONJG(inv(INDINV(irow, irow : icol - 1))))) / mat(icol, icol)
170 : #else
171 : #error "Unrecognized interface."
172 : #endif
173 : end do
174 : end do
175 88 : inv(irow, irow) = ONE / mat(irow, irow)
176 461 : do irow = 1_IK, ndim
177 1523 : do icol = irow, ndim
178 3783 : inv(INDINV(irow, icol)) = dot_product(inv(INDINV(icol, icol : ndim)), inv(INDINV(irow, icol : ndim)))
179 : #if FUL_ENABLED
180 465 : inv(INDINV(icol, irow)) = GET_CONJG(inv(INDINV(irow, icol)))
181 : #endif
182 : end do
183 : end do
184 19 : elseif (ndim == 1_IK) then
185 19 : inv(1,1) = ONE / mat(1,1)**2
186 : end if
187 : #undef INDINV
188 : #undef INDMAT
189 :
190 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
191 : #elif setMatInv_ENABLED && LUP_ENABLED
192 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
193 :
194 : TYPE_KIND :: summ
195 : integer(IK) :: irow, icol, idim, ndim
196 180 : ndim = size(mat, 1, IK)
197 180 : CHECK_SHAPE_MAT
198 1980 : CHECK_SHAPE_INV
199 1080 : CHECK_ASSERTION(__LINE__, size(auxil, 1, IK) == ndim, SK_"@setMatInv(): The condition `size(auxil) == size(inv, 1)` must hold. size(auxil), shape(inv) = "//getStr([size(auxil, 1, IK), shape(inv, IK)]))
200 713 : do icol = 1, ndim
201 1376 : inv(1 : icol - 1, icol) = ZERO
202 533 : inv(icol, icol) = ONE
203 1376 : inv(icol + 1 : ndim, icol) = ZERO
204 : idim = 0_IK
205 2752 : do irow = 1_IK, ndim
206 2219 : summ = inv(auxil(irow), icol)
207 2219 : inv(auxil(irow), icol) = inv(irow, icol)
208 2219 : if (idim /= 0_IK) then
209 2816 : summ = summ - dot_product(GET_CONJG(mat(irow, idim : irow - 1)), inv(idim : irow - 1, icol))
210 1376 : elseif (summ /= ZERO) then
211 : idim = irow
212 : end if
213 2752 : inv(irow, icol) = summ
214 : end do
215 2932 : do irow = ndim, 1_IK, -1_IK
216 2219 : CHECK_ASSERTION(__LINE__, mat(irow, irow) /= ZERO, SK_"@setMatInv(): The condition `mat(irow, irow) /= ZERO` must hold. irow = "//getStr(irow))
217 7828 : inv(irow, icol) = (inv(irow, icol) - dot_product(GET_CONJG(mat(irow, irow + 1 : ndim)), inv(irow + 1 : ndim, icol))) / mat(irow, irow)
218 : end do
219 : end do
220 :
221 : #else
222 : !%%%%%%%%%%%%%%%%%%%%%%%%
223 : #error "Unrecognized interface."
224 : !%%%%%%%%%%%%%%%%%%%%%%%%
225 : #endif
226 : #undef CHECK_SHAPE_INV
227 : #undef CHECK_SHAPE_MAT
228 : #undef TYPE_KIND
229 : #undef GET_CONJG
|
