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_matrixLUP](@ref pm_matrixLUP).
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 : !real(TKC), parameter :: IHUGE = 1._TKC / huge(0._TKC)
28 : !real(TKC), parameter :: SMALL = merge(tiny(0._TKC), IHUGE * (1._TKC + epsilon(0._TKC)), IHUGE < tiny(0._TKC))
29 : real(TKC), parameter :: SMALL = tiny(0._TKC)**.999
30 : #if CK_ENABLED
31 : #define GET_CONJG(X) conjg(X)
32 : #define TYPE_KIND complex(TKC)
33 : #define GET_ABSL1(X) abs(X%re) + abs(X%im)
34 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
35 : #elif RK_ENABLED
36 : #define GET_CONJG(X) X
37 : #define GET_ABSL1(X) abs(X)
38 : #define TYPE_KIND real(TKC)
39 : real(TKC), parameter :: ZERO = 0._TKC, ONE = 1._TKC
40 : #else
41 : #error "Unrecognized interface."
42 : #endif
43 :
44 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45 : #if setMatLUP_ENABLED && SQM_ENABLED && LAPACK_ENABLED && DISPATCH_ENABLED
46 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47 :
48 : !check_assertion(__LINE__, all(size(rperm, 1, IK) == shape(mat, IK)), \
49 : !SK_"@setMatLUP(): The condition `all(size(rperm) == shape(mat))` must hold. size(rperm), shape(mat) = "//\
50 : !getStr([size(rperm, 1, IK), shape(mat, IK)])) ! fpp
51 : call lapackGETRF(size(mat, 1, IK), size(mat, 2, IK), mat(1,1), size(mat, 1, IK), rperm, info)
52 : !if (present(parity)) then
53 : ! block
54 : ! integer(IK) :: idim
55 : ! parity = ONE
56 : ! do idim = 1, size(rperm, 1, IK)
57 : ! if (rperm(idim) /= idim) parity = -parity
58 : ! end do
59 : ! end block
60 : !end if
61 :
62 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 : #elif setMatLUP_ENABLED && SQM_ENABLED
64 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65 :
66 : TYPE_KIND :: dumm, summ !, parity_def
67 336 : real(TKC) :: rowMax, temp, rowMaxInv(size(rperm, 1, IK))
68 : integer(IK) :: irow, imax, icol, idim, ndim
69 2688 : CHECK_ASSERTION(__LINE__, all(size(rperm, 1, IK) == shape(mat, IK)), \
70 : SK_"@setMatLUP(): The condition `all(size(rperm) == shape(mat))` must hold. size(rperm), shape(mat) = "//\
71 : getStr([size(rperm, 1, IK), shape(mat, IK)])) ! fpp
72 336 : info = 0_IK
73 : !parity_def = 1_IK
74 : ndim = size(rperm, 1, IK)
75 336 : if (ndim == 0_IK) return ! warning: the value of parity is undefined on return here.
76 1607 : do info = 1, ndim
77 : ! either this or the following after works, with the L1 norm perhaps more robust.
78 8969 : rowMaxInv(info) = abs(mat(maxloc(abs(mat(1 : ndim, info)), 1, kind = IK), info))
79 : !rowMaxInv(info) = abs(mat(maxloc(GET_ABSL1(mat(1 : ndim, info)), 1, kind = IK), info))
80 1271 : if (rowMaxInv(info) == ZERO) return
81 1607 : rowMaxInv(info) = 1._TKC / rowMaxInv(info)
82 : end do
83 336 : info = 0_IK
84 1607 : do icol = 1, ndim
85 3849 : do irow = 1, icol - 1
86 2578 : summ = mat(irow, icol)
87 6382 : do idim = 1, irow - 1
88 6382 : summ = summ - mat(idim, icol) * mat(irow, idim)
89 : end do
90 3849 : mat(irow, icol) = summ
91 : end do
92 525 : rowMax = 0._TKC
93 5120 : do irow = icol, ndim
94 3849 : summ = mat(irow, icol)
95 10231 : do idim = 1, icol - 1
96 10231 : summ = summ - mat(idim, icol) * mat(irow, idim)
97 : end do
98 3849 : mat(irow, icol) = summ
99 3849 : temp = rowMaxInv(irow) * abs(summ)
100 5120 : if (rowMax <= temp) then
101 915 : rowMax = temp
102 : imax = irow
103 : end if
104 : end do
105 1271 : if (icol /= imax)then
106 2876 : do idim = 1, ndim
107 2492 : dumm = mat(imax, idim)
108 2492 : mat(imax, idim) = mat(icol, idim)
109 2876 : mat(icol, idim) = dumm
110 : end do
111 : !parity_def = -parity_def
112 384 : rowMaxInv(imax) = rowMaxInv(icol)
113 : end if
114 1271 : rperm(icol) = imax
115 1271 : if (mat(icol, icol) == ZERO) then
116 0 : mat(icol, icol) = SMALL
117 0 : info = icol
118 : end if
119 1607 : if (icol /= ndim) then
120 935 : dumm = ONE / mat(icol, icol)
121 3513 : do irow = icol + 1, ndim
122 3513 : mat(irow, icol) = mat(irow, icol) * dumm
123 : end do
124 : endif
125 : end do
126 : !if (present(parity)) parity = parity_def
127 :
128 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 : #elif setMatLUP_ENABLED && REC_ENABLED
130 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 :
132 : TYPE_KIND :: scaler
133 : integer(IK) :: irow, nmin, nminHalf
134 : info = 0_IK
135 : ! Quick return if possible
136 : if (nrow == 0_IK .or. ncol == 0_IK) return
137 : if (nrow == 1_IK) then
138 : ! Use unblocked code for ONE row case just need to handle rperm and info.
139 : rperm(1) = 1_IK
140 : if (mat(1, 1) == ZERO) info = 1_IK
141 : elseif (ncol == 1_IK) then
142 : ! Use unblocked code for ONE column case and find pivot and test for singularity.
143 : irow = maxloc(GET_ABSL1(mat(1 : nrow)), 1)
144 : rperm(1) = irow
145 : if (mat(irow, 1) /= ZERO) then
146 : ! apply the interchange.
147 : if (irow /= 1_IK) then
148 : temp = mat(1, 1)
149 : mat(1, 1) = mat(irow, 1)
150 : mat(irow, 1) = temp
151 : end if
152 : ! compute elements 2:nrow of the column
153 : if (SMALL <= abs(mat(1, 1))) then
154 : scaler = ONE / mat(1, 1)
155 : do irow = 2_IK, nrow
156 : mat(irow, 1) = mat(irow, 1) * scaler
157 : end do
158 : else
159 : do irow = 2_IK, nrow
160 : mat(irow, 1) = mat(irow, 1) / mat(1, 1)
161 : end do
162 : end if
163 : else
164 : info = 1_IK
165 : end if
166 : else
167 : ! use recursive code.
168 : nmin = min(nrow, ncol)
169 : nminHalf = nmin / 2_IK
170 : ndif = ncol - nminHalf
171 : ! [ a11 ]
172 : ! factor [ --- ]
173 : ! [ a21 ]
174 : call zgetrf2(nrow, nminHalf, mat, lda, rperm, iinfo)
175 : if (info == 0_IK .and. 0_IK < iinfo) info = iinfo
176 : ! [ a12 ]
177 : ! apply interchanges to [ --- ]
178 : ! [ a22 ]
179 : call setMatPerm(mat(:, nminHalf + 1 : ndimHalf + ndif), rperm(:), roff, incr)
180 : !call ZLASWP(ndif, mat(1, nminHalf + 1), lda, 1, nminHalf, rperm, 1)
181 : ! solve a12
182 : call ztrsm('l', 'l', 'ncol', 'u', nminHalf, ndif, ONE, mat, lda, mat(1, nminHalf + 1), lda)
183 : ! update a22
184 : call zgemm('ncol', 'ncol', nrow - nminHalf, ndif, nminHalf, -ONE, mat(nminHalf + 1, 1), lda, mat(1, nminHalf + 1), lda, ONE, mat(nminHalf + 1, nminHalf + 1), lda)
185 : ! factor a22
186 : call zgetrf2(nrow - nminHalf, ndif, mat(nminHalf + 1, nminHalf + 1), lda, rperm(nminHalf + 1), iinfo)
187 : ! adjust info and the pivot indices
188 : if (info == 0_IK .and. 0_IK < iinfo) info = iinfo + nminHalf
189 : do irow = nminHalf + 1_IK, nminHalf
190 : rperm(irow) = rperm(irow) + nmin
191 : end do
192 : ! apply interchanges to a21
193 : call zlaswp(nminHalf, mat(1, 1), lda, nminHalf + 1_IK, nminHalf, rperm, 1_IK)
194 : end if
195 :
196 : !%%%%%%%%%%%%%%%
197 : #elif Blocking_ENABLED
198 : !%%%%%%%%%%%%%%%
199 :
200 : integer(IK) :: nmin, bdim_def, iinfo, irow, icol, jcol
201 : #if IMP_ENABLED
202 : integer(IK) :: nrow, ncol
203 : nrow = size(mat, 1, IK)
204 : ncol = size(mat, 2, IK)
205 : #elif EXP_ENABLED
206 : CHECK_ASSERTION(__LINE__, all([nrow, ncol] <= ubound(mat, kind = IK)), \
207 : SK_"@setMatLUP(): The condition `all([nrow + roff, ncol + coff] <= shape(mat, kind = IK))` must hold. nrow, roff, ncol, coff, shape(mat) = "//\
208 : getStr([nrow, roff, ncol, coff, shape(mat, IK)])) ! fpp
209 : #else
210 : #error "Unrecognized interface."
211 : #endif
212 : nmin = min(nrow, ncol)
213 : CHECK_ASSERTION(__LINE__, size(rperm, 1, IK) == nmin, \
214 : SK_"@setMatLUP(): The condition `size(rperm) == min(shape(mat))` must hold. size(rperm), shape(mat) = "//\
215 : getStr([size(rperm, 1, IK), shape(mat, IK)])) ! fpp
216 : info = 0_IK
217 : ! Quick return if possible.
218 : if (nmin == 0_IK) return
219 : if (present(bdim)) then
220 : CHECK_ASSERTION(__LINE__, 0_IK < bdim, SK_"@setMatChol(): The condition `0 < bdim` must hold. bdim = "//getStr(bdim))
221 : bdim_def = bdim
222 : else
223 : bdim_def = 64_IK
224 : end if
225 : if (bdim_def <= 1_IK .or. nmin <= bdim_def) ! use unblocked code.
226 : #if IMP_ENABLED
227 : call setMatLUP(mat, rperm, info, recursion)
228 : #elif EXP_ENABLED
229 : call setMatLUP(mat, rperm, info, recursion, nrow, ncol, roff, coff)
230 : #endif
231 : else ! use blocked code.
232 : do icol = 1_IK, nmin, bdim_def
233 : jcol = min(nmin - icol + 1_IK, bdim_def)
234 : ! Factor diagonal and subdiagonal blocks and test for exact singularity.
235 : call setMatLUP(mat, rperm(icol), iinfo, recursion, nrow - icol + 1_IK, jcol, icol - 1_IK, icol - 1_IK)
236 : ! Adjust info and the pivot indices.
237 : if (info == 0_IK .and. 0_IK < iinfo) info = iinfo + icol - 1_IK
238 : do irow = icol, min(nrow, icol + jcol - 1_IK)
239 : rperm(irow) = icol - 1 + rperm(irow)
240 : end do
241 : ! Apply interchanges to columns 1 : icol - 1.
242 : call zlaswp(icol - 1_IK, mat, lda, icol, icol + jcol - 1_IK, rperm, 1_IK)
243 : if (icol + jcol <= ncol) then
244 : ! Apply interchanges to columns icol + jcol:ncol.
245 : call zlaswp(ncol - icol - jcol + 1_IK, mat(1, icol + jcol), lda, icol, icol + jcol - 1_IK, rperm, 1_IK)
246 : ! Compute block row of `u`.
247 : call ztrsm('left', 'lower', 'no transpose', 'unit', jcol, ncol - icol - jcol + 1_IK, ONE, mat(icol, icol), lda, mat(icol, icol + jcol), lda)
248 : if (icol + jcol <= nrow) then
249 : ! Update trailing submatrix.
250 : call zgemm ( 'no transpose', 'no transpose' &
251 : , nrow - icol - jcol + 1_IK, ncol - icol - jcol + 1_IK, jcol, -ONE &
252 : , mat(icol + jcol, icol), lda, mat(icol, icol + jcol) &
253 : , lda, ONE, mat(icol + jcol, icol + jcol), lda &
254 : )
255 : end if
256 : end if
257 : end do
258 : end if
259 :
260 : #else
261 : !%%%%%%%%%%%%%%%%%%%%%%%%
262 : #error "Unrecognized interface."
263 : !%%%%%%%%%%%%%%%%%%%%%%%%
264 : #endif
265 : #undef TYPE_KIND
266 : #undef GET_ABSL1
267 : #undef GET_CONJG
|
