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_matrixUpdate](@ref pm_matrixUpdate).
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 : ! Define the BLAS storage triangle.
28 : #if shrk_ENABLED && SLD_ENABLED
29 : #define BLAS_UPLO "L"
30 : #elif shrk_ENABLED && SUD_ENABLED
31 : #define BLAS_UPLO "U"
32 : #elif !(setMatUpdateR1_ENABLED || setMatUpdateTriang_ENABLED)
33 : #error "Unrecognized interface."
34 : #endif
35 : ! Define the BLAS storage triangle.
36 : #if shrk_ENABLED && OTP_ENABLED && CHM_ENABLED && CK_ENABLED
37 : #define blasXXRK blasHERK
38 : #define BLAS_TRANS "C"
39 : #elif shrk_ENABLED && ONO_ENABLED && CHM_ENABLED && CK_ENABLED
40 : #define blasXXRK blasHERK
41 : #define BLAS_TRANS "N"
42 : #elif shrk_ENABLED && OTP_ENABLED && (CHM_ENABLED || CSM_ENABLED) && (CK_ENABLED || RK_ENABLED)
43 : #define blasXXRK blasSYRK
44 : #define BLAS_TRANS "T"
45 : #elif shrk_ENABLED && ONO_ENABLED && (CHM_ENABLED || CSM_ENABLED) && (CK_ENABLED || RK_ENABLED)
46 : #define blasXXRK blasSYRK
47 : #define BLAS_TRANS "N"
48 : #elif !(IK_ENABLED || setMatUpdateR1_ENABLED || setMatUpdateTriang_ENABLED)
49 : #error "Unrecognized interface."
50 : #endif
51 : ! Define the constants.
52 : #if IK_ENABLED
53 : #define TYPE_KIND integer(TKC)
54 : integer(TKC), parameter :: ZERO = 0_TKC, ONE = 1_TKC
55 : #elif CK_ENABLED
56 : #define TYPE_KIND complex(TKC)
57 : complex(TKC), parameter :: ZERO = (0._TKC, 0._TKC), ONE = (1._TKC, 0._TKC)
58 : #elif RK_ENABLED
59 : #define TYPE_KIND real(TKC)
60 : real(TKC) , parameter :: ZERO = 0._TKC, ONE = 1._TKC
61 : #else
62 : #error "Unrecognized interface."
63 : #endif
64 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
65 : #if setMatUpdate_ENABLED && shrk_ENABLED
66 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
67 :
68 : #if !(BLAS_ENABLED && DISPATCH_ENABLED)
69 : TYPE_KIND :: temp
70 : integer(IK) :: irow, icol, idum
71 : #endif
72 : #if ASS_ENABLED
73 : #define BETA beta_def
74 : #define ALPHA alpha_def
75 : integer(IK) :: ndim, ndum
76 : TYPE_KIND :: alpha_def, beta_def
77 6 : if (present(beta)) then; beta_def = beta; else; beta_def = ONE; end if
78 6 : if (present(alpha)) then; alpha_def = alpha; else; alpha_def = ONE; end if
79 6 : ndim = size(mat, 1, IK)
80 6 : CHECK_ASSERTION(__LINE__, size(mat, 1, IK) == size(mat, 2, IK), SK_"@setMatUpdate(): The condition `size(mat, 1) == size(mat, 2)` must hold. shape(mat) = "//getStr(shape(mat, IK)))
81 : #if ONO_ENABLED
82 4 : ndum = size(matA, 2, IK)
83 36 : CHECK_ASSERTION(__LINE__, size(mat, 2, IK) == size(matA, 1, IK), SK_"@setMatUpdate(): The condition `size(mat, 2) == size(matA, 1)` must hold. shape(mat), shape(matA) = "//getStr([shape(mat, IK), shape(matA, IK)]))
84 : #elif OTP_ENABLED
85 2 : ndum = size(matA, 1, IK)
86 18 : CHECK_ASSERTION(__LINE__, size(mat, 2, IK) == size(matA, 2, IK), SK_"@setMatUpdate(): The condition `size(mat, 2) == size(matA, 2)` must hold. shape(mat), shape(matA) = "//getStr([shape(mat, IK), shape(matA, IK)]))
87 : #else
88 : #error "Unrecognized interface."
89 : #endif
90 : #elif EXP_ENABLED
91 : ! Ensure roffs and subset bounds are logical.
92 74193 : CHECK_ASSERTION(__LINE__, all(0_IK <= [ndim, ndum]), SK_"@setMatUpdate(): The condition `all(0_IK <= [ndim, ndum])` must hold. ndim, ndum = "//getStr([ndim, ndum]))
93 173117 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(mat, kind = IK)), SK_"@setMatUpdate(): The condition `all(0 <= [roff, coff])` must hold. roff, coff = "//getStr([1_IK - lbound(mat)]))
94 173117 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(matA, kind = IK)), SK_"@setMatUpdate(): The condition `all(0 <= [roffA, coffA])` must hold. roffA, coffA = "//getStr([1_IK - lbound(matA)]))
95 : #if ONO_ENABLED
96 160758 : CHECK_ASSERTION(__LINE__, all([ndim, ndum] <= shape(matA, IK)), SK_"@setMatUpdate(): The condition `all([roffA + ndim, coffA + ndum] <= shape(matA))` must hold. roffA, coffA, ndim, ndum, shape(matA) = "//getStr([1_IK - lbound(matA, kind = IK), ndim, ndum, shape(matA, IK)]))
97 : #elif OTP_ENABLED
98 160745 : CHECK_ASSERTION(__LINE__, all([ndum, ndim] <= shape(matA, IK)), SK_"@setMatUpdate(): The condition `all([roffA + ndum, coffA + ndim] <= shape(matA))` must hold. roffA, coffA, ndim, ndum, shape(matA) = "//getStr([1_IK - lbound(matA, kind = IK), ndim, ndum, shape(matA, IK)]))
99 : #else
100 : #error "Unrecognized interface."
101 : #endif
102 : #else
103 : #error "Unrecognized interface."
104 : #endif
105 : ! Set the updating components.
106 : #if CHM_ENABLED && CK_ENABLED
107 : #define GET_CONJG(X) conjg(X)
108 : #define GET_RE(X) X%re
109 : #define GET_REAL(X)real(X, TKC)
110 : #elif CHM_ENABLED || (CSM_ENABLED && (ONO_ENABLED || OTP_ENABLED))
111 : #define GET_CONJG(X) X
112 : #define GET_REAL(X)X
113 : #define GET_RE(X) X
114 : #else
115 : #error "Unrecognized interface."
116 : #endif
117 : ! Set loop bounds.
118 : #if SLD_ENABLED
119 : #define POFFSET + 1
120 : #define GET_SLICE(i,j,k) j, k
121 : #define SLD_SET(X,Y) X = Y
122 : #define SUD_SET(X,Y)
123 : #elif SUD_ENABLED
124 : #define POFFSET - 1
125 : #define SLD_SET(X,Y)
126 : #define SUD_SET(X,Y) X = Y
127 : #define GET_SLICE(i,j,k) i, j
128 : #else
129 : #error "Unrecognized interface."
130 : #endif
131 : #if BLAS_ENABLED && DISPATCH_ENABLED
132 : call blasXXRK(BLAS_UPLO, BLAS_TRANS, ndim, ndum, GET_RE(ALPHA), matA(1,1), size(matA, 1, IK), GET_RE(BETA), mat(1,1), size(mat, 1, IK))
133 : #else
134 : ! Quick return if possible.
135 24737 : if (ndim == 0_IK .or. ((GET_RE(ALPHA) == GET_RE(ZERO) .or. ndum == 0_IK) .and. GET_RE(BETA) == GET_RE(ONE))) return
136 24120 : if (GET_RE(ALPHA) == GET_RE(ZERO)) then
137 0 : if (GET_RE(BETA) == GET_RE(ZERO)) then
138 0 : do icol = 1, ndim
139 0 : do irow = GET_SLICE(1, icol, ndim)
140 0 : mat(irow, icol) = ZERO
141 : end do
142 : end do
143 : else
144 0 : do icol = 1, ndim
145 0 : SLD_SET(mat(icol, icol), GET_RE(BETA) * GET_RE(mat(icol, icol)))
146 0 : do irow = GET_SLICE(1, icol POFFSET, ndim)
147 0 : mat(irow, icol) = GET_RE(BETA) * mat(irow, icol)
148 : end do
149 0 : SUD_SET(mat(icol, icol), GET_RE(BETA) * GET_RE(mat(icol, icol)))
150 : end do
151 : end if
152 0 : return
153 : end if
154 : ! Perform Update.
155 : #if ONO_ENABLED
156 : ! Form mat := GET_RE(ALPHA) * matA * matA ** T + GET_RE(BETA) * mat.
157 : ! Form mat := GET_RE(ALPHA) * matA * matA ** H + GET_RE(BETA) * mat.
158 32287 : do icol = 1, ndim
159 20225 : if (GET_RE(BETA) == GET_RE(ZERO)) then
160 0 : do irow = GET_SLICE(1, icol, ndim)
161 0 : mat(irow, icol) = ZERO
162 : end do
163 20225 : else if (GET_RE(BETA) /= GET_RE(ONE)) then
164 0 : SLD_SET(mat(icol, icol), GET_RE(BETA) * GET_RE(mat(icol, icol)))
165 12 : do irow = GET_SLICE(1, icol POFFSET, ndim)
166 12 : mat(irow, icol) = GET_RE(BETA) * mat(irow, icol)
167 : end do
168 6 : SUD_SET(mat(icol, icol), GET_RE(BETA) * GET_RE(mat(icol, icol)))
169 : #if CHM_ENABLED && CK_ENABLED
170 : else
171 9728 : mat(icol, icol)%im = ZERO%im
172 : #endif
173 : end if
174 65064 : do idum = 1, ndum
175 53002 : if (matA(icol, idum) /= ZERO) then
176 88 : temp = GET_RE(ALPHA) * GET_CONJG(matA(icol, idum))
177 19467 : SLD_SET(mat(icol, icol), GET_RE(mat(icol, icol)) + GET_REAL(temp * matA(icol, idum)))
178 33672 : do irow = GET_SLICE(1, icol POFFSET, ndim)
179 33672 : mat(irow, icol) = mat(irow, icol) + temp * matA(irow, idum)
180 : end do
181 88 : SUD_SET(mat(icol, icol), GET_RE(mat(icol, icol)) + GET_REAL(temp * matA(icol, idum)))
182 : end if
183 : end do
184 : end do
185 : #elif OTP_ENABLED && CSM_ENABLED
186 : ! Form mat := GET_RE(ALPHA) * matA ** T * matA + GET_RE(BETA) * mat.
187 : ! This block is indeed the same as the block for HA_ENABLED.
188 : ! However, it is kept separately for now for its cleaner implementation.
189 18 : do icol = 1, ndim
190 90 : do irow = GET_SLICE(1, icol, ndim)
191 0 : temp = ZERO
192 288 : do idum = 1, ndum
193 288 : temp = temp + matA(idum, irow) * matA(idum, icol)
194 : end do
195 88 : if (GET_RE(BETA) == ZERO) then
196 0 : mat(irow, icol) = GET_RE(ALPHA) * temp
197 : else
198 72 : mat(irow, icol) = GET_RE(ALPHA) * temp + GET_RE(BETA) * mat(irow, icol)
199 : end if
200 : end do
201 : end do
202 : #elif OTP_ENABLED && CHM_ENABLED
203 : ! Form mat := GET_RE(ALPHA) * matA ** H * matA + GET_RE(BETA) * mat.
204 : #define SET_DIAGONAL \
205 : GET_RE(temp) = GET_RE(ZERO); \
206 : do idum = 1, ndum; \
207 : GET_RE(temp) = GET_RE(temp) + GET_REAL(GET_CONJG(matA(idum, icol)) * matA(idum, icol)); \
208 : end do; \
209 : if (GET_RE(BETA) == GET_RE(ZERO)) then; \
210 : mat(icol, icol) = GET_RE(ALPHA) * GET_RE(temp); \
211 : else; \
212 : mat(icol, icol) = GET_RE(ALPHA) * GET_RE(temp) + GET_RE(BETA) * GET_RE(mat(icol, icol)); \
213 : end if;
214 : #define SET_OFFDIAG \
215 : do irow = GET_SLICE(1, icol POFFSET, ndim); \
216 : temp = ZERO; \
217 : do idum = 1, ndum; \
218 : temp = temp + GET_CONJG(matA(idum, irow)) * matA(idum, icol); \
219 : end do; \
220 : if (GET_RE(BETA) == GET_RE(ZERO)) then; \
221 : mat(irow, icol) = GET_RE(ALPHA) * temp; \
222 : else; \
223 : mat(irow, icol) = GET_RE(ALPHA) * temp + GET_RE(BETA) * mat(irow, icol); \
224 : end if; \
225 : end do;
226 32252 : do icol = 1, ndim
227 : #if SLD_ENABLED
228 36 : SET_DIAGONAL
229 44 : SET_OFFDIAG
230 : #elif SUD_ENABLED
231 55389 : SET_OFFDIAG
232 64899 : SET_DIAGONAL
233 : #endif
234 : end do
235 : #else
236 : #error "Unrecognized interface."
237 : #endif
238 : #endif
239 :
240 : !%%%%%%%%%%%%%%%%%%%%%
241 : #elif setMatUpdateR1_ENABLED
242 : !%%%%%%%%%%%%%%%%%%%%%
243 :
244 : integer(IK) :: i, j, jy, ix, kx, effLenX, srow, erow
245 : TYPE_KIND :: temp
246 : ! Define transpose type for complex arguments.
247 : #if setMatUpdateR1H_CK_ENABLED || setMatUpdateR1AH_CK_ENABLED
248 : #define GET_CONJG(X) conjg(X)
249 : #else
250 : #define GET_CONJG(X) X
251 : #endif
252 : ! Define the `alpha` factor of the transformation.
253 : #if Fixed_ENABLED
254 : #define ALPHA_TIMES
255 : #elif Alpha_ENABLED
256 : #define ALPHA_TIMES alpha *
257 1 : if (alpha == ZERO) return ! Quick return if possible.
258 : #else
259 : #error "Unrecognized interface."
260 : #endif
261 : kx = 1_IK
262 6 : if (incA /= 1_IK) then
263 1 : CHECK_ASSERTION(__LINE__, incA /= 0_IK, SK_"@setMatUpdateR1(): The condition `incA /= 0_IK` must hold. incA = "//getStr(incA)) ! fpp
264 1 : effLenX = (size(vecA, 1, IK) - 1_IK) / abs(incA) + 1_IK ! the effective length of `vecA`.
265 3 : CHECK_ASSERTION(__LINE__, incA * effLenX /= size(vecA, 1, IK), \
266 : SK_"@setMatUpdateR1(): The condition `size(vecA) == incA * ((size(vecA, 1) - 1) / abs(incA) + 1)` must hold. incA, size(vecA) = "//\
267 : getStr([incA, size(vecA, 1, IK)])) ! fpp
268 1 : if (incA < 0_IK) kx = size(vecA, 1, IK)
269 : else
270 5 : effLenX = size(vecA, 1, IK)
271 : end if
272 :
273 : jy = 1_IK
274 6 : if (incB /= 1_IK) then
275 3 : CHECK_ASSERTION(__LINE__, incB /= 0_IK, SK_"@setMatUpdateR1(): The condition `incB /= 0_IK` must hold. incB = "//getStr(incB)) ! fpp
276 3 : if (incB < 0_IK) jy = size(vecB, 1, IK)
277 : end if
278 :
279 6 : CHECK_ASSERTION(__LINE__, 0_IK <= roff, SK_"@setMatUpdateR1(): The conditions `0_IK <= roff` must hold. roff = "//getStr(roff)) ! fpp
280 :
281 18 : CHECK_ASSERTION(__LINE__, size(mat, 2, IK) == (size(vecB, 1, IK) - 1_IK) / abs(incB) + 1_IK, \
282 : SK_"@setMatUpdateR1(): The condition `size(mat, 2) == (size(vecB) - 1) / abs(incB) + 1` must hold. size(mat, 2), (size(vecB) - 1) / abs(incB) + 1 = "//\
283 : getStr([size(mat, 2, IK), (size(vecB, 1, IK) - 1_IK) / abs(incB) + 1_IK])) ! fpp
284 30 : CHECK_ASSERTION(__LINE__, size(mat, 1, IK) >= effLenX + roff, \
285 : SK_"@setMatUpdateR1(): The conditions `size(mat, 1) >= (size(vecA, 1) - 1) / abs(incA) + 1 + roff` must hold. size(mat, 1), size(vecA), incA, roff = "//\
286 : getStr([size(mat, 1, IK), size(vecA, 1, IK), incA, roff])) ! fpp
287 :
288 : ! Start the operations. In this version the elements of `mat` are accessed sequentially with one pass through `mat`.
289 :
290 6 : srow = roff + 1_IK
291 : erow = roff + effLenX
292 6 : if (incA == 1_IK) then
293 20 : do j = 1_IK, size(mat, 2, IK)
294 15 : if (vecB(jy) /= ZERO) then
295 3 : temp = ALPHA_TIMES GET_CONJG(vecB(jy))
296 81 : do i = srow, erow
297 81 : mat(i,j) = mat(i,j) + vecA(i) * temp
298 : end do
299 : end if
300 20 : jy = jy + incB
301 : end do
302 : else
303 4 : do j = 1, size(mat, 2, IK)
304 3 : if (vecB(jy) /= ZERO) then
305 0 : temp = ALPHA_TIMES GET_CONJG(vecB(jy))
306 : ix = kx
307 15 : do i = srow, erow
308 12 : mat(i,j) = mat(i,j) + vecA(ix) * temp
309 15 : ix = ix + incA
310 : end do
311 : end if
312 4 : jy = jy + incB
313 : end do
314 : end if
315 : #undef ALPHA_TIMES
316 :
317 : !%%%%%%%%%%%%%%%%%%%%%%%%%
318 : #elif setMatUpdateTriang_ENABLED
319 : !%%%%%%%%%%%%%%%%%%%%%%%%%
320 :
321 : integer(IK) :: irow, icol, idum
322 : TYPE_KIND :: temp
323 :
324 : ! Ensure roffs and subset bounds are logical.
325 18 : CHECK_ASSERTION(__LINE__, all(0_IK <= [ndim, ndum]), \
326 : SK_"@setMatUpdateTriang(): The condition `all(0_IK <= [ndim, ndum])` must hold. ndim, ndum = "//getStr([ndim, ndum])) ! fpp
327 42 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(matA, kind = IK)), \
328 : SK_"@setMatUpdateTriang(): The condition `all(0 <= [roffA, coffA])` must hold. roffA, coffA = "//getStr([1_IK - lbound(matA)])) ! fpp
329 42 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(mat, kind = IK)), \
330 : SK_"@setMatUpdateTriang(): The condition `all(0 <= [roff, coff])` must hold. roff, coff = "//getStr([1_IK - lbound(mat)])) ! fpp
331 : ! Ensure subset of `matA` does not overflow the matrix bounds.
332 : #if AS_ENABLED || AH_ENABLED
333 52 : CHECK_ASSERTION(__LINE__, all([ndim, ndum] <= shape(matA, kind = IK)), \
334 : SK_"@setMatUpdateTriang(): The condition `all([roffA + ndim, coffA + ndum] <= shape(matA))` must hold. roffA, coffA, ndim, ndum, shape(matA) = "\
335 : //getStr([1_IK - lbound(matA, kind = IK), ndim, ndum, shape(matA, IK)])) ! fpp
336 : #elif SA_ENABLED || HA_ENABLED
337 26 : CHECK_ASSERTION(__LINE__, all([ndum, ndim] <= shape(matA,IK)), \
338 : SK_"@setMatUpdateTriang(): The condition `all([roffA + ndum, coffA + ndim] <= shape(matA))` must hold. roffA, coffA, ndim, ndum, shape(matA) = "\
339 : //getStr([1_IK - lbound(matA, kind = IK), ndim, ndum, shape(matA, IK)])) ! fpp
340 : #else
341 : #error "Unrecognized interface."
342 : #endif
343 : ! Set the updating components.
344 : #if AS_ENABLED || SA_ENABLED
345 : #define GET_CONJG(X) X
346 : #define GET_REAL(X)X
347 : #define GET_RE(X) X
348 : #elif AH_ENABLED || HA_ENABLED
349 : #define GET_CONJG(X) conjg(X)
350 : #define GET_RE(X) X%re
351 : #define GET_REAL(X)real(X, TKC)
352 : #else
353 : #error "Unrecognized interface."
354 : #endif
355 : ! Set loop bounds.
356 : #if lowDiaC_ENABLED
357 : #define POFFSET + 1
358 : #define GET_SLICE(i,j,k) j, k
359 : #define lowDiaC_SET(X,Y) X = Y
360 : #define uppDiaC_SET(X,Y)
361 : #elif uppDiaC_ENABLED
362 : #define POFFSET - 1
363 : #define lowDiaC_SET(X,Y)
364 : #define uppDiaC_SET(X,Y) X = Y
365 : #define GET_SLICE(i,j,k) i, j
366 : #else
367 : #error "Unrecognized interface."
368 : #endif
369 : ! Quick return if possible.
370 6 : if (ndim == 0_IK .or. ((alpha == ZERO .or. ndum == 0_IK) .and. beta == ONE)) return
371 6 : if (alpha == GET_RE(ZERO)) then
372 0 : if (beta == GET_RE(ZERO)) then
373 0 : do icol = 1, ndim
374 0 : do irow = GET_SLICE(1, icol, ndim)
375 0 : mat(irow, icol) = ZERO
376 : end do
377 : end do
378 : else
379 0 : do icol = 1, ndim
380 0 : lowDiaC_SET(mat(icol, icol), beta * GET_RE(mat(icol, icol)))
381 0 : do irow = GET_SLICE(1, icol POFFSET, ndim)
382 0 : mat(irow, icol) = beta * mat(irow, icol)
383 : end do
384 0 : uppDiaC_SET(mat(icol, icol), beta * GET_RE(mat(icol, icol)))
385 : end do
386 : end if
387 0 : return
388 : end if
389 : ! Perform Update.
390 : #if AS_ENABLED || AH_ENABLED
391 : ! Form mat := alpha * matA * matA ** T + beta * mat.
392 : ! Form mat := alpha * matA * matA ** H + beta * mat.
393 21 : do icol = 1, ndim
394 17 : if (beta == GET_RE(ZERO)) then
395 0 : do irow = GET_SLICE(1, icol, ndim)
396 0 : mat(irow, icol) = ZERO
397 : end do
398 17 : else if (beta /= GET_RE(ONE)) then
399 0 : lowDiaC_SET(mat(icol, icol), beta * GET_RE(mat(icol, icol)))
400 6 : do irow = GET_SLICE(1, icol POFFSET, ndim)
401 6 : mat(irow, icol) = beta * mat(irow, icol)
402 : end do
403 3 : uppDiaC_SET(mat(icol, icol), beta * GET_RE(mat(icol, icol)))
404 : #if AH_ENABLED
405 : else
406 6 : mat(icol, icol)%im = ZERO%im
407 : #endif
408 : end if
409 82 : do idum = 1, ndum
410 78 : if (matA(icol, idum) /= ZERO) then
411 44 : temp = alpha * GET_CONJG(matA(icol, idum))
412 15 : lowDiaC_SET(mat(icol, icol), GET_RE(mat(icol, icol)) + GET_REAL(temp * matA(icol, idum)))
413 160 : do irow = GET_SLICE(1, icol POFFSET, ndim)
414 160 : mat(irow, icol) = mat(irow, icol) + temp * matA(irow, idum)
415 : end do
416 44 : uppDiaC_SET(mat(icol, icol), GET_RE(mat(icol, icol)) + GET_REAL(temp * matA(icol, idum)))
417 : end if
418 : end do
419 : end do
420 : #elif SA_ENABLED
421 : ! Form mat := alpha * matA ** T * matA + beta * mat.
422 : ! This block is indeed the same as the block for HA_ENABLED.
423 : ! However, it is kept separately for now for its cleaner implementation.
424 9 : do icol = 1, ndim
425 45 : do irow = GET_SLICE(1, icol, ndim)
426 0 : temp = ZERO
427 144 : do idum = 1, ndum
428 144 : temp = temp + matA(idum, irow) * matA(idum, icol)
429 : end do
430 44 : if (beta == ZERO) then
431 0 : mat(irow, icol) = alpha * temp
432 : else
433 36 : mat(irow, icol) = alpha * temp + beta * mat(irow, icol)
434 : end if
435 : end do
436 : end do
437 : #elif HA_ENABLED
438 : ! Form mat := alpha * matA ** H * matA + beta * mat.
439 : #define SET_DIAGONAL \
440 : GET_RE(temp) = GET_RE(ZERO); \
441 : do idum = 1, ndum; \
442 : GET_RE(temp) = GET_RE(temp) + GET_REAL(GET_CONJG(matA(idum, icol)) * matA(idum, icol)); \
443 : end do; \
444 : if (beta == GET_RE(ZERO)) then; \
445 : mat(icol, icol) = alpha * GET_RE(temp); \
446 : else; \
447 : mat(icol, icol) = alpha * GET_RE(temp) + beta * GET_RE(mat(icol, icol)); \
448 : end if;
449 : #define SET_OFFDIAG \
450 : do irow = GET_SLICE(1, icol POFFSET, ndim); \
451 : temp = ZERO; \
452 : do idum = 1, ndum; \
453 : temp = temp + GET_CONJG(matA(idum, irow)) * matA(idum, icol); \
454 : end do; \
455 : if (beta == GET_RE(ZERO)) then; \
456 : mat(irow, icol) = alpha * temp; \
457 : else; \
458 : mat(irow, icol) = alpha * temp + beta * mat(irow, icol); \
459 : end if; \
460 : end do;
461 4 : do icol = 1, ndim
462 : #if lowDiaC_ENABLED
463 18 : SET_DIAGONAL
464 22 : SET_OFFDIAG
465 : #elif uppDiaC_ENABLED
466 0 : SET_OFFDIAG
467 0 : SET_DIAGONAL
468 : #endif
469 : end do
470 : #else
471 : #error "Unrecognized interface."
472 : #endif
473 :
474 : #else
475 : !%%%%%%%%%%%%%%%%%%%%%%%%
476 : #error "Unrecognized interface."
477 : !%%%%%%%%%%%%%%%%%%%%%%%%
478 : #endif
479 : #undef SET_DIAGONAL
480 : #undef SET_OFFDIAG
481 : #undef lowDiaC_SET
482 : #undef uppDiaC_SET
483 : #undef BLAS_TRANS
484 : #undef BLAS_UPLO
485 : #undef TYPE_KIND
486 : #undef GET_CONJG
487 : #undef GET_SLICE
488 : #undef blasXXRK
489 : #undef GET_REAL
490 : #undef SLD_SET
491 : #undef SUD_SET
492 : #undef POFFSET
493 : #undef GET_RE
494 : #undef ALPHA
495 : #undef BETA
|
