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_matrixMulTri](@ref pm_matrixMulTri).
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 dispatch function.
28 : #if trmv_ENABLED
29 : #define blasTRXX blasTRMV
30 : #elif trsv_ENABLED
31 : #define blasTRXX blasTRSV
32 : #elif trmm_ENABLED
33 : #define blasTRXX blasTRMM
34 : #elif trsm_ENABLED
35 : #define blasTRXX blasTRSM
36 : #else
37 : #error "Unrecognized interface."
38 : #endif
39 : ! Define the blas dispatch transposition.
40 : #if (CGMB_ENABLED && ONOA_ENABLED) || (CGMA_ENABLED && ONOB_ENABLED)
41 : #define BLAS_TRANSA "N"
42 : #elif (CGMB_ENABLED && INVA_ENABLED) || (CGMA_ENABLED && INVB_ENABLED)
43 : #define BLAS_TRANSA "N"
44 : #elif (CGMB_ENABLED && (OTSA_ENABLED || OTOA_ENABLED)) || (CGMA_ENABLED && (OTSB_ENABLED || OTOB_ENABLED))
45 : #define BLAS_TRANSA "T"
46 : #elif (CGMB_ENABLED && (OTHA_ENABLED || OTUA_ENABLED)) || (CGMA_ENABLED && (OTHB_ENABLED || OTUB_ENABLED))
47 : #define BLAS_TRANSA "C"
48 : #else
49 : #error "Unrecognized interface."
50 : #endif
51 : ! Define the blas dispatch matrix class.
52 : #if (trmv_ENABLED || trsv_ENABLED || trmm_ENABLED || trsm_ENABLED) && ((CGMB_ENABLED && (CLDA_ENABLED || CLUA_ENABLED)) || (CGMA_ENABLED && (CLDB_ENABLED || CLUB_ENABLED)))
53 : #define BLAS_UPLO "L"
54 : #elif (trmv_ENABLED || trsv_ENABLED || trmm_ENABLED || trsm_ENABLED) && ((CGMB_ENABLED && (CUDA_ENABLED || CUUA_ENABLED)) || (CGMA_ENABLED && (CUDB_ENABLED || CUUB_ENABLED)))
55 : #define BLAS_UPLO "U"
56 : #else
57 : #error "Unrecognized interface."
58 : #endif
59 : ! Define the constants.
60 : #if CK_ENABLED
61 : #define TYPE_KIND complex(CKC)
62 : complex(CKC), parameter :: ZERO = (0._CKC, 0._CKC), ONE = (1._CKC, 0._CKC)
63 : #elif RK_ENABLED
64 : #define TYPE_KIND real(RKC)
65 : real(RKC), parameter :: ZERO = 0._RKC, ONE = 1._RKC
66 : #else
67 : #error "Unrecognized interface."
68 : #endif
69 : ! Define the input and output matrices.
70 : #if CGMA_ENABLED
71 : #define BLAS_SIDE "R"
72 : #define SOLMAT matA
73 : #define TRIMAT matB
74 : #elif CGMB_ENABLED
75 : #define BLAS_SIDE "L"
76 : #define SOLMAT matB
77 : #define TRIMAT matA
78 : #else
79 : #error "Unrecognized interface."
80 : #endif
81 : ! Define the conjugation of the triangular matrix.
82 : #if CK_ENABLED && ((OTUA_ENABLED || OTUB_ENABLED) || (OTHA_ENABLED || OTHB_ENABLED))
83 : #define GET_CONJG(X)conjg(X)
84 : #elif (ONOA_ENABLED || ONOB_ENABLED) || (INVA_ENABLED || INVB_ENABLED) || (OTOA_ENABLED || OTOB_ENABLED) || (OTUA_ENABLED || OTUB_ENABLED)
85 : #define GET_CONJG(X)X
86 : #else
87 : #error "Unrecognized interface."
88 : #endif
89 : ! Define the diag/unit triangular matrices.
90 : #if (CGMB_ENABLED && (CLDA_ENABLED || CUDA_ENABLED)) || (CGMA_ENABLED && (CLDB_ENABLED || CUDB_ENABLED))
91 : #define BLAS_DIAG "N"
92 : #define CXD_ENABLED 1
93 : #elif (CGMB_ENABLED && (CLUA_ENABLED || CUUA_ENABLED)) || (CGMA_ENABLED && (CLUB_ENABLED || CUUB_ENABLED))
94 : #define BLAS_DIAG "U"
95 : #define CXD_ENABLED 0
96 : #else
97 : #error "Unrecognized interface."
98 : #endif
99 :
100 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
101 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
102 : #if setMatMulTri_ENABLED && (trmv_ENABLED || trsv_ENABLED)
103 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105 :
106 : #if !(BLAS_ENABLED && DISPATCH_ENABLED)
107 : integer(IK) :: irow, icol, bstart, ib, jb
108 : TYPE_KIND :: temp
109 : #endif
110 : ! Define the assumed shape of the matrices.
111 : #if ASS_ENABLED
112 : integer(IK) :: nrow
113 : integer(IK) , parameter :: incB = 1
114 19206 : nrow = size(TRIMAT, 1, IK)
115 19206 : CHECK_ASSERTION(__LINE__, size(matA, 1, IK) == size(matA, 2, IK), SK_"@setMatMulTri(): The condition `size(matA, 1) == size(matA, 2)` must hold. shape(matA) = "//getStr(shape(matA, IK))) ! fpp
116 134442 : CHECK_ASSERTION(__LINE__, size(matA, 1, IK) == size(matB, 1, IK), SK_"@setMatMulTri(): The condition `size(matA, 1) == size(matB, 1)` must hold. shape(matA), shape(matB) = "//getStr([shape(matA), shape(matB)])) ! fpp
117 : #elif EXP_ENABLED
118 : ! Ensure offsets and subset bounds are logical.
119 19206 : CHECK_ASSERTION(__LINE__, 0_IK <= nrow, SK_"@setMatMulTri(): The condition `0 <= nrow` must hold. nrow = "//getStr(nrow)) ! fpp
120 19206 : CHECK_ASSERTION(__LINE__, 0_IK /= incB, SK_"@setMatMulTri(): The condition `0 /= incB` must hold. nrow = "//getStr(incB)) ! fpp
121 96030 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(matA, kind = IK)), SK_"@setMatMulTri(): The condition `all(0 <= [roffA, coffA])` must hold. roffA, coffA = "//getStr(1_IK - lbound(matA, kind = IK))) ! fpp
122 76824 : CHECK_ASSERTION(__LINE__, abs(incB) * (nrow - 1) < size(matB, 1, IK), SK_"@setMatMulTri(): The condition `size(matB(1::abs(incB))) <= nrow` must hold. incB, nrow, size(matB) = "//getStr([incB, nrow, size(matB, kind = IK)])) ! fpp
123 230472 : CHECK_ASSERTION(__LINE__, all(nrow <= ubound(matA, kind = IK)), SK_"@setMatMulTri(): The condition `all(nrow == shape(matA) - [roffA, coffA])` must hold. nrow, roffA, coffA, shape(matA) = "//getStr([nrow, 1_IK - lbound(matA, kind = IK), shape(matA, IK)])) ! fpp
124 : ! Set up the start point in matB if the increment is not unity.
125 : ! This will be (nrow - 1) * incB too small for descending loops.
126 : #else
127 : #error "Unrecognized interface."
128 : #endif
129 : ! Quick return if possible.
130 38412 : if (nrow == 0_IK) return
131 : #if BLAS_ENABLED && DISPATCH_ENABLED
132 : call blasTRXX(BLAS_UPLO, BLAS_TRANSA, BLAS_DIAG, nrow, TRIMAT(1,1), size(TRIMAT, 1, IK), SOLMAT(1), incB)
133 : #else
134 : ! Start the operations.
135 : ! In this version the elements of matA are accessed sequentially with one pass through matA.
136 16350 : if (incB < 0_IK) then
137 7944 : bstart = 1 - (nrow - 1) * incB
138 8406 : elseif (incB /= 1) then
139 : bstart = 1
140 : end if
141 :
142 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143 : #if CGMB_ENABLED && ONOB_ENABLED && ONOA_ENABLED
144 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 :
146 : ! BLAS 2 TRMV - Form matB := matA * matB.
147 : #if CUDA_ENABLED || CUUA_ENABLED
148 1362 : if (incB == 1) then
149 7938 : do icol = 1, nrow
150 7938 : if (matB(icol) /= ZERO) then
151 4080 : temp = matB(icol)
152 19228 : do irow = 1, icol - 1
153 19228 : matB(irow) = matB(irow) + temp * matA(irow, icol)
154 : end do
155 : #if CXD_ENABLED
156 3174 : matB(icol) = matB(icol) * matA(icol, icol)
157 : #endif
158 : end if
159 : end do
160 : else
161 : jb = bstart
162 5746 : do icol = 1, nrow
163 4612 : if (matB(jb) /= ZERO) then
164 2912 : temp = matB(jb)
165 : ib = bstart
166 14152 : do irow = 1, icol - 1
167 9540 : matB(ib) = matB(ib) + temp * matA(irow, icol)
168 14152 : ib = ib + incB
169 : end do
170 : #if CXD_ENABLED
171 2306 : matB(jb) = matB(jb) * matA(icol, icol)
172 : #endif
173 : end if
174 5746 : jb = jb + incB
175 : end do
176 : end if
177 : #elif CLDA_ENABLED || CLUA_ENABLED
178 1363 : if (incB == 1) then
179 7943 : do icol = nrow, 1, -1
180 7943 : if (matB(icol) /= ZERO) then
181 4080 : temp = matB(icol)
182 19238 : do irow = nrow, icol + 1, -1
183 19238 : matB(irow) = matB(irow) + temp * matA(irow, icol)
184 : end do
185 : #if CXD_ENABLED
186 3174 : matB(icol) = matB(icol) * matA(icol, icol)
187 : #endif
188 : end if
189 : end do
190 : else
191 1135 : bstart = bstart + (nrow - 1) * incB
192 : jb = bstart
193 5751 : do icol = nrow, 1, -1
194 4616 : if (matB(jb) /= ZERO) then
195 2912 : temp = matB(jb)
196 : ib = bstart
197 14162 : do irow = nrow, icol + 1, -1
198 9546 : matB(ib) = matB(ib) + temp * matA(irow, icol)
199 14162 : ib = ib - incB
200 : end do
201 : #if CXD_ENABLED
202 2306 : matB(jb) = matB(jb) * matA(icol, icol)
203 : #endif
204 : end if
205 5751 : jb = jb - incB
206 : end do
207 : end if
208 : #else
209 : #error "Unrecognized interface."
210 : #endif
211 :
212 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213 : #elif CGMB_ENABLED && ONOB_ENABLED && INVA_ENABLED
214 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215 :
216 : ! BLAS TRSV - Form matB := inv(matA) * matB.
217 :
218 : #if CUDA_ENABLED || CUUA_ENABLED
219 1362 : if (incB == 1) then
220 7938 : do icol = nrow, 1, -1
221 7938 : if (matB(icol) /= ZERO) then
222 : #if CXD_ENABLED
223 3174 : matB(icol) = matB(icol) / matA(icol, icol)
224 : #endif
225 4080 : temp = matB(icol)
226 19228 : do irow = icol - 1, 1, -1
227 19228 : matB(irow) = matB(irow) - temp * matA(irow,icol)
228 : end do
229 : end if
230 : end do
231 : else
232 1134 : jb = bstart + (nrow - 1) * incB
233 5746 : do icol = nrow, 1, -1
234 4612 : if (matB(jb) /= ZERO) then
235 : #if CXD_ENABLED
236 2306 : matB(jb) = matB(jb) / matA(icol, icol)
237 : #endif
238 2912 : temp = matB(jb)
239 : ib = jb
240 14152 : do irow = icol - 1, 1, -1
241 9540 : ib = ib - incB
242 14152 : matB(ib) = matB(ib) - temp * matA(irow, icol)
243 : end do
244 : end if
245 5746 : jb = jb - incB
246 : end do
247 : end if
248 : #elif CLDA_ENABLED || CLUA_ENABLED
249 1363 : if (incB == 1) then
250 7943 : do icol = 1, nrow
251 7943 : if (matB(icol) /= ZERO) then
252 : #if CXD_ENABLED
253 3174 : matB(icol) = matB(icol) / matA(icol, icol)
254 : #endif
255 4080 : temp = matB(icol)
256 19238 : do irow = icol + 1, nrow
257 19238 : matB(irow) = matB(irow) - temp * matA(irow, icol)
258 : end do
259 : end if
260 : end do
261 : else
262 : jb = bstart
263 5751 : do icol = 1, nrow
264 4616 : if (matB(jb) /= ZERO) then
265 : #if CXD_ENABLED
266 2306 : matB(jb) = matB(jb) / matA(icol, icol)
267 : #endif
268 2912 : temp = matB(jb)
269 : ib = jb
270 14162 : do irow = icol + 1, nrow
271 9546 : ib = ib + incB
272 14162 : matB(ib) = matB(ib) - temp * matA(irow, icol)
273 : end do
274 : end if
275 5751 : jb = jb + incB
276 : end do
277 : end if
278 : #else
279 : #error "Unrecognized interface."
280 : #endif
281 :
282 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283 : #elif CGMB_ENABLED && ONOB_ENABLED && (OTSA_ENABLED || OTHA_ENABLED)
284 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
285 :
286 : ! BLAS TRMV - Form matB := matA^T * matB or matB := matA^H * matB.
287 :
288 : #if CUDA_ENABLED || CUUA_ENABLED
289 2726 : if (incB == 1) then
290 15896 : do icol = nrow, 1, -1
291 12712 : temp = matB(icol)
292 : #if CXD_ENABLED
293 6356 : temp = temp * GET_CONJG(matA(icol, icol))
294 : #endif
295 38496 : do irow = icol - 1, 1, -1
296 38496 : temp = temp + GET_CONJG(matA(irow, icol)) * matB(irow)
297 : end do
298 15896 : matB(icol) = temp
299 : end do
300 : else
301 2268 : jb = bstart + (nrow - 1) * incB
302 11492 : do icol = nrow, 1, -1
303 9224 : temp = matB(jb)
304 : ib = jb
305 : #if CXD_ENABLED
306 4612 : temp = temp*GET_CONJG(matA(icol,icol))
307 : #endif
308 28304 : do irow = icol - 1, 1, -1
309 19080 : ib = ib - incB
310 28304 : temp = temp + GET_CONJG(matA(irow, icol)) * matB(ib)
311 : end do
312 9224 : matB(jb) = temp
313 11492 : jb = jb - incB
314 : end do
315 : end if
316 : #elif CLDA_ENABLED || CLUA_ENABLED
317 2724 : if (incB == 1) then
318 15876 : do icol = 1, nrow
319 12696 : temp = matB(icol)
320 : #if CXD_ENABLED
321 6348 : temp = temp * GET_CONJG(matA(icol, icol))
322 : #endif
323 38456 : do irow = icol + 1, nrow
324 38456 : temp = temp + GET_CONJG(matA(irow, icol)) * matB(irow)
325 : end do
326 15876 : matB(icol) = temp
327 : end do
328 : else
329 : jb = bstart
330 11492 : do icol = 1, nrow
331 9224 : temp = matB(jb)
332 : ib = jb
333 : #if CXD_ENABLED
334 4612 : temp = temp * GET_CONJG(matA(icol, icol))
335 : #endif
336 28304 : do irow = icol + 1, nrow
337 19080 : ib = ib + incB
338 28304 : temp = temp + GET_CONJG(matA(irow,icol))*matB(ib)
339 : end do
340 9224 : matB(jb) = temp
341 11492 : jb = jb + incB
342 : end do
343 : end if
344 : #else
345 : #error "Unrecognized interface."
346 : #endif
347 :
348 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
349 : #elif CGMB_ENABLED && ONOB_ENABLED && (OTOA_ENABLED || OTUA_ENABLED)
350 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
351 :
352 : ! BLAS 2 TRSV - Form matB := inv(matA^T) * matB or matB := inv(matA^H) * matB.
353 : #if CUDA_ENABLED || CUUA_ENABLED
354 2726 : if (incB == 1) then
355 15896 : do icol = 1, nrow
356 12712 : temp = matB(icol)
357 38496 : do irow = 1, icol - 1
358 38496 : temp = temp - GET_CONJG(matA(irow, icol)) * matB(irow)
359 : end do
360 : #if CXD_ENABLED
361 6356 : temp = temp / GET_CONJG(matA(icol, icol))
362 : #endif
363 15896 : matB(icol) = temp
364 : end do
365 : else
366 : jb = bstart
367 11492 : do icol = 1, nrow
368 : ib = bstart
369 9224 : temp = matB(jb)
370 28304 : do irow = 1, icol - 1
371 19080 : temp = temp - GET_CONJG(matA(irow, icol)) * matB(ib)
372 28304 : ib = ib + incB
373 : end do
374 : #if CXD_ENABLED
375 4612 : temp = temp / GET_CONJG(matA(icol, icol))
376 : #endif
377 9224 : matB(jb) = temp
378 11492 : jb = jb + incB
379 : end do
380 : end if
381 : #elif CLDA_ENABLED || CLUA_ENABLED
382 2724 : if (incB == 1) then
383 15876 : do icol = nrow, 1, -1
384 12696 : temp = matB(icol)
385 38456 : do irow = nrow,icol + 1, -1
386 38456 : temp = temp - GET_CONJG(matA(irow, icol)) * matB(irow)
387 : end do
388 : #if CXD_ENABLED
389 6348 : temp = temp / GET_CONJG(matA(icol, icol))
390 : #endif
391 15876 : matB(icol) = temp
392 : end do
393 : else
394 2268 : bstart = bstart + (nrow-1)*incB
395 : jb = bstart
396 11492 : do icol = nrow, 1, -1
397 : ib = bstart
398 9224 : temp = matB(jb)
399 28304 : do irow = nrow, icol + 1, -1
400 19080 : temp = temp - GET_CONJG(matA(irow, icol)) * matB(ib)
401 28304 : ib = ib - incB
402 : end do
403 : #if CXD_ENABLED
404 4612 : temp = temp / GET_CONJG(matA(icol, icol))
405 : #endif
406 9224 : matB(jb) = temp
407 11492 : jb = jb - incB
408 : end do
409 : end if
410 : #else
411 : #error "Unrecognized interface."
412 : #endif
413 :
414 : #else
415 : !%%%%%%%%%%%%%%%%%%%%%%%%
416 : #error "Unrecognized interface."
417 : !%%%%%%%%%%%%%%%%%%%%%%%%
418 : #endif
419 : #endif
420 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
421 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
422 : #elif setMatMulTri_ENABLED && (trmm_ENABLED || trsm_ENABLED)
423 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
424 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425 :
426 : #if !(BLAS_ENABLED && DISPATCH_ENABLED)
427 : integer(IK) :: irow, icol, idum
428 : #endif
429 : ! Define the assumed shape of the matrices.
430 : #if ASS_ENABLED
431 : #define ALPHA alpha_def
432 : integer(IK) :: nrow, ncol
433 : TYPE_KIND :: alpha_def
434 76812 : if (present(alpha)) then
435 38402 : alpha_def = alpha
436 : else
437 24000 : alpha_def = ONE
438 : end if
439 76812 : nrow = size(SOLMAT, 1, IK)
440 38404 : ncol = size(SOLMAT, 2, IK)
441 : #elif EXP_ENABLED
442 : ! Ensure offsets and subset bounds are logical.
443 188523 : CHECK_ASSERTION(__LINE__, all(0_IK <= [nrow, ncol]), SK_"@setMatMulTri(): The condition `all(0_IK <= [nrow, ncol])` must hold. nrow, ncol = "//getStr([nrow, ncol])) ! fpp
444 314205 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(matA, kind = IK)), SK_"@setMatMulTri(): The condition `all(0 <= [roffA, coffA])` must hold. roffA, coffA = "//getStr(1_IK - lbound(matA, kind = IK))) ! fpp
445 314205 : CHECK_ASSERTION(__LINE__, all(0_IK <= 1_IK - lbound(matB, kind = IK)), SK_"@setMatMulTri(): The condition `all(0 <= [roffB, coffB])` must hold. roffB, coffB = "//getStr(1_IK - lbound(matB, kind = IK))) ! fpp
446 : #else
447 : #error "Unrecognized interface."
448 : #endif
449 : ! Ensure subset of `TRIMAT` matches the corresponding subset of `SOLMAT`.
450 : #if CGMA_ENABLED
451 830974 : CHECK_ASSERTION(__LINE__, all([nrow, ncol] <= ubound(matA, kind = IK)), SK_"@setMatMulTri(): The condition `all([nrow + roffA, ncol + coffA] <= shape(matA, kind = IK))` must hold. nrow, ncol, roffA, coffA, shape(matA) = "//getStr([nrow, ncol, 1_IK - lbound(matA, kind = IK), shape(matA, IK)])) ! fpp
452 761144 : CHECK_ASSERTION(__LINE__, all(ncol <= ubound(matB, kind = IK)), SK_"@setMatMulTri(): The condition `all([roffB, coffB] + ncol <= shape(matB))` must hold. roffB, coffB, ncol, shape(matB) = "//getStr([1_IK - lbound(matB, kind = IK), ncol, shape(matB, IK)])) ! fpp
453 : #elif CGMB_ENABLED
454 761068 : CHECK_ASSERTION(__LINE__, all(nrow <= ubound(matA, kind = IK)), \
455 : SK_"@setMatMulTri(): The condition `all([roffA, coffA] + nrow <= shape(matA))` must hold. roffA, coffA, nrow, shape(matA) = "//getStr([1_IK - lbound(matA, kind = IK), nrow, shape(matA, IK)])) ! fpp
456 1180006 : CHECK_ASSERTION(__LINE__, all([nrow, ncol] <= ubound(matB, kind = IK)), SK_"@setMatMulTri(): The condition `all([nrow + roffB, ncol + coffB] <= shape(matB, kind = IK))` must hold. nrow, ncol, roffB, coffB, shape(matB) = "//getStr([nrow, ncol, 1_IK - lbound(matB, kind = IK), shape(matB, IK), ubound(matB, kind = IK)])) ! fpp
457 : #endif
458 : ! Set the iteration limits based on the storage format.
459 : #if (CUDA_ENABLED || CUUA_ENABLED) || (CUDB_ENABLED || CUUB_ENABLED)
460 : #define GET_HALF_RANGE(i,j,k) i, j - 1
461 : #elif (CLDA_ENABLED || CLUA_ENABLED) || (CLDB_ENABLED || CLUB_ENABLED)
462 : #define GET_HALF_RANGE(i,j,k) j + 1, k
463 : #else
464 : #error "Unrecognized interface."
465 : #endif
466 : ! Quick return if possible.
467 139653 : if (nrow == 0_IK .or. ncol == 0_IK) return
468 110820 : if (ALPHA == ZERO) then
469 379200 : SOLMAT(1 : nrow, 1 : ncol) = ZERO
470 : return
471 : end if
472 : #if BLAS_ENABLED && DISPATCH_ENABLED
473 : call blasTRXX(BLAS_SIDE, BLAS_UPLO, BLAS_TRANSA, BLAS_DIAG, nrow, ncol, ALPHA, TRIMAT(1,1), size(TRIMAT, 1, IK), SOLMAT(1,1), size(SOLMAT, 1, IK))
474 : #else
475 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
476 : #if trmm_ENABLED && CGMB_ENABLED && ONOB_ENABLED && ONOA_ENABLED
477 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
478 :
479 : ! BLAS 3 TRMM - Form SOLMAT := alpha * TRIMAT * SOLMAT.
480 : block
481 : TYPE_KIND :: temp
482 : #if CUDA_ENABLED || CUUA_ENABLED
483 13990 : do icol = 1, ncol
484 59232 : do idum = 1, nrow
485 56410 : if (SOLMAT(idum, icol) /= ZERO) then
486 45240 : temp = ALPHA * SOLMAT(idum, icol)
487 138738 : do irow = 1, idum - 1
488 138738 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * TRIMAT(irow, idum)
489 : end do
490 : #if CXD_ENABLED
491 22634 : temp = temp * TRIMAT(idum, idum)
492 : #endif
493 45240 : SOLMAT(idum, icol) = temp
494 : end if
495 : end do
496 : end do
497 : #elif CLDA_ENABLED || CLUA_ENABLED
498 13982 : do icol = 1, ncol
499 59194 : do idum = nrow, 1, -1
500 56374 : if (SOLMAT(idum, icol) /= ZERO) then
501 45212 : temp = ALPHA * SOLMAT(idum, icol)
502 22606 : SOLMAT(idum, icol) = temp
503 : #if CXD_ENABLED
504 22606 : SOLMAT(idum, icol) = SOLMAT(idum, icol) * TRIMAT(idum, idum)
505 : #endif
506 138654 : do irow = idum + 1, nrow
507 138654 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * TRIMAT(irow, idum)
508 : end do
509 : end if
510 : end do
511 : end do
512 : #else
513 : #error "Unrecognized interface."
514 : #endif
515 : end block
516 :
517 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
518 : #elif trsm_ENABLED && CGMB_ENABLED && ONOB_ENABLED && INVA_ENABLED
519 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
520 :
521 : ! BLAS 3 TRSM - Form SOLMAT := alpha * inv(TRIMAT) * SOLMAT.
522 : #if CUDA_ENABLED || CUUA_ENABLED
523 : #define GET_FULL_RANGE(i,j) j, i, -1
524 : #elif CLDA_ENABLED || CLUA_ENABLED
525 : #define GET_FULL_RANGE(i,j) i, j
526 : #else
527 : #error "Unrecognized interface."
528 : #endif
529 27972 : do icol = 1, ncol
530 22330 : if (ALPHA /= ONE) then
531 32728 : do irow = 1, nrow
532 32728 : SOLMAT(irow, icol) = ALPHA * SOLMAT(irow, icol)
533 : end do
534 : end if
535 118426 : do idum = GET_FULL_RANGE(1, nrow)
536 112784 : if (SOLMAT(idum, icol) /= ZERO) then
537 : #if CXD_ENABLED
538 22634 : SOLMAT(idum, icol) = SOLMAT(idum, icol) / TRIMAT(idum, idum)
539 : #endif
540 277392 : do irow = GET_HALF_RANGE(1, idum, nrow)
541 277392 : SOLMAT(irow, icol) = SOLMAT(irow, icol) - SOLMAT(idum, icol) * TRIMAT(irow, idum)
542 : end do
543 : end if
544 : end do
545 : end do
546 :
547 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
548 : #elif trmm_ENABLED && CGMB_ENABLED && ONOB_ENABLED && (OTSA_ENABLED || OTHA_ENABLED)
549 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
550 :
551 : ! BLAS 3 TRMM - Form SOLMAT := alpha * TRIMAT^T * SOLMAT or SOLMAT := alpha * TRIMAT^H * SOLMAT
552 :
553 : block
554 : TYPE_KIND :: temp
555 55938 : do icol = 1, ncol
556 : #if CUDA_ENABLED || CUUA_ENABLED
557 118438 : do irow = nrow, 1, -1
558 45212 : temp = SOLMAT(irow, icol)
559 : #if CXD_ENABLED
560 45252 : temp = temp * GET_CONJG(TRIMAT(irow, irow))
561 : #endif
562 277428 : do idum = 1, irow - 1
563 277428 : temp = temp + GET_CONJG(TRIMAT(idum, irow)) * SOLMAT(idum, icol)
564 : end do
565 112796 : SOLMAT(irow, icol) = ALPHA * temp
566 : end do
567 : #elif CLDA_ENABLED || CLUA_ENABLED
568 118388 : do irow = 1, nrow
569 90424 : temp = SOLMAT(irow, icol)
570 : #if CXD_ENABLED
571 45212 : temp = temp * GET_CONJG(TRIMAT(irow, irow))
572 : #endif
573 277308 : do idum = irow + 1, nrow
574 277308 : temp = temp + GET_CONJG(TRIMAT(idum,irow))*SOLMAT(idum, icol)
575 : end do
576 112748 : SOLMAT(irow, icol) = ALPHA*temp
577 : end do
578 : #else
579 : #error "Unrecognized interface."
580 : #endif
581 : end do
582 : end block
583 :
584 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
585 : #elif trsm_ENABLED && CGMB_ENABLED && ONOB_ENABLED && (OTOA_ENABLED || OTUA_ENABLED)
586 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
587 :
588 : ! BLAS 3 TRSM - Form SOLMAT := alpha * inv(transpose(TRIMAT)) * SOLMAT
589 : ! BLAS 3 TRSM - Form SOLMAT := alpha * inv(conjg(transpose(TRIMAT))) * SOLMAT
590 : #if CUDA_ENABLED || CUUA_ENABLED
591 : #define GET_FULL_RANGE(i,j) i, j
592 : #elif CLDA_ENABLED || CLUA_ENABLED
593 : #define GET_FULL_RANGE(i,j) j, i, -1
594 : #else
595 : #error "Unrecognized interface."
596 : #endif
597 : !error stop "SOLMAT := alpha * inv(transpose(TRIMAT)) * SOLMAT"
598 : block
599 : TYPE_KIND :: temp
600 88329 : do icol = 1, ncol
601 301945 : do irow = GET_FULL_RANGE(1, nrow)
602 123192 : temp = ALPHA * SOLMAT(irow, icol)
603 604146 : do idum = GET_HALF_RANGE(1, irow, nrow)
604 604146 : temp = temp - GET_CONJG(TRIMAT(idum, irow)) * SOLMAT(idum, icol)
605 : end do
606 : #if CXD_ENABLED
607 123192 : temp = temp / GET_CONJG(TRIMAT(irow, irow))
608 : #endif
609 278609 : SOLMAT(irow, icol) = temp
610 : end do
611 : end do
612 : end block
613 :
614 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
615 : #elif trmm_ENABLED && CGMA_ENABLED && ONOA_ENABLED && ONOB_ENABLED
616 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
617 :
618 : ! BLAS 3 TRMM - Form SOLMAT := alpha * SOLMAT * TRIMAT.
619 :
620 : block
621 : TYPE_KIND :: temp
622 : #if CUDB_ENABLED || CUUB_ENABLED
623 14008 : do icol = ncol, 1, -1
624 3413 : temp = ALPHA
625 : #if CXD_ENABLED
626 5591 : temp = temp * TRIMAT(icol, icol)
627 : #endif
628 56438 : do irow = 1, nrow
629 56438 : SOLMAT(irow, icol) = temp * SOLMAT(irow, icol)
630 : end do
631 36314 : do idum = 1, icol - 1
632 33490 : if (TRIMAT(idum, icol) /= ZERO) then
633 22300 : temp = ALPHA * TRIMAT(idum, icol)
634 113468 : do irow = 1, nrow
635 113468 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * SOLMAT(irow, idum)
636 : end do
637 : end if
638 : end do
639 : end do
640 : #elif CLDB_ENABLED || CLUB_ENABLED
641 14006 : do icol = 1, ncol
642 3413 : temp = ALPHA
643 : #if CXD_ENABLED
644 5601 : temp = temp * TRIMAT(icol, icol)
645 : #endif
646 56454 : do irow = 1, nrow
647 56454 : SOLMAT(irow, icol) = temp * SOLMAT(irow, icol)
648 : end do
649 36302 : do idum = icol + 1, ncol
650 33478 : if (TRIMAT(idum, icol) /= ZERO) then
651 22290 : temp = ALPHA * TRIMAT(idum, icol)
652 113480 : do irow = 1, nrow
653 113480 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * SOLMAT(irow, idum)
654 : end do
655 : end if
656 : end do
657 : end do
658 : #else
659 : #error "Unrecognized interface."
660 : #endif
661 : end block
662 :
663 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
664 : #elif trsm_ENABLED && CGMA_ENABLED && ONOA_ENABLED && INVB_ENABLED
665 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
666 :
667 : ! BLAS 3 TRSM - Form SOLMAT := alpha * SOLMAT * inv(TRIMAT).
668 : #if CUDB_ENABLED || CUUB_ENABLED
669 : #define GET_FULL_RANGE(i,j) i, j
670 : #elif CLDB_ENABLED || CLUB_ENABLED
671 : #define GET_FULL_RANGE(i,j) j, i, -1
672 : #else
673 : #error "Unrecognized interface."
674 : #endif
675 : block
676 : #if CXD_ENABLED
677 : TYPE_KIND :: temp
678 : #endif
679 28014 : do icol = GET_FULL_RANGE(1, ncol)
680 22366 : if (ALPHA /= ONE) then
681 32728 : do irow = 1, nrow
682 32728 : SOLMAT(irow,icol) = ALPHA * SOLMAT(irow, icol)
683 : end do
684 : end if
685 69790 : do idum = GET_HALF_RANGE(1, icol, ncol)
686 66968 : if (TRIMAT(idum, icol) /= ZERO) then
687 226948 : do irow = 1, nrow
688 226948 : SOLMAT(irow, icol) = SOLMAT(irow, icol) - TRIMAT(idum, icol) * SOLMAT(irow, idum)
689 : end do
690 : end if
691 : end do
692 : #if CXD_ENABLED
693 11192 : temp = ONE / TRIMAT(icol, icol)
694 59320 : do irow = 1, nrow
695 56494 : SOLMAT(irow, icol) = temp * SOLMAT(irow, icol)
696 : end do
697 : #endif
698 : end do
699 : end block
700 :
701 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
702 : #elif trmm_ENABLED && CGMA_ENABLED && ONOA_ENABLED && (OTSB_ENABLED || OTHB_ENABLED)
703 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
704 :
705 : ! BLAS 3 TRMM - Form SOLMAT := alpha * SOLMAT * TRIMAT**T or SOLMAT := alpha * SOLMAT * TRIMAT^H.
706 :
707 : block
708 : TYPE_KIND :: temp
709 : #if CUDB_ENABLED || CUUB_ENABLED
710 27964 : do idum = 1, ncol
711 66836 : do icol = 1, idum - 1
712 66836 : if (TRIMAT(icol, idum) /= ZERO) then
713 44512 : temp = ALPHA * GET_CONJG(TRIMAT(icol, idum))
714 226688 : do irow = 1, nrow
715 226688 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * SOLMAT(irow, idum)
716 : end do
717 : end if
718 : end do
719 6826 : temp = ALPHA
720 : #if CXD_ENABLED
721 11162 : temp = temp * GET_CONJG(TRIMAT(idum, idum))
722 : #endif
723 27964 : if (temp /= ONE) then
724 72738 : do irow = 1, nrow
725 72738 : SOLMAT(irow, idum) = temp * SOLMAT(irow, idum)
726 : end do
727 : end if
728 : end do
729 : #elif CLDB_ENABLED || CLUB_ENABLED
730 27964 : do idum = ncol, 1, -1
731 66836 : do icol = idum + 1, ncol
732 66836 : if (TRIMAT(icol, idum) /= ZERO) then
733 44512 : temp = ALPHA * GET_CONJG(TRIMAT(icol, idum))
734 226688 : do irow = 1, nrow
735 226688 : SOLMAT(irow, icol) = SOLMAT(irow, icol) + temp * SOLMAT(irow, idum)
736 : end do
737 : end if
738 : end do
739 6826 : temp = ALPHA
740 : #if CXD_ENABLED
741 11162 : temp = temp * GET_CONJG(TRIMAT(idum, idum))
742 : #endif
743 27964 : if (temp /= ONE) then
744 72738 : do irow = 1, nrow
745 72738 : SOLMAT(irow, idum) = temp * SOLMAT(irow, idum)
746 : end do
747 : end if
748 : end do
749 : #else
750 : #error "Unrecognized interface."
751 : #endif
752 : end block
753 :
754 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
755 : #elif trsm_ENABLED && CGMA_ENABLED && ONOA_ENABLED && (OTOB_ENABLED || OTUB_ENABLED)
756 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
757 :
758 : ! BLAS 3 TRSM - Form SOLMAT := alpha * SOLMAT * inv(TRIMAT**T).
759 : ! BLAS 3 TRSM - Form SOLMAT := alpha * SOLMAT * inv(TRIMAT**H).
760 : #if CUDB_ENABLED || CUUB_ENABLED
761 : #define GET_FULL_RANGE(i,j) j, i, -1
762 : #elif CLDB_ENABLED || CLUB_ENABLED
763 : #define GET_FULL_RANGE(i,j) i, j
764 : #else
765 : #error "Unrecognized interface."
766 : #endif
767 : block
768 : TYPE_KIND :: temp
769 84268 : do idum = GET_FULL_RANGE(1, ncol)
770 : #if CXD_ENABLED
771 38610 : temp = ONE / GET_CONJG(TRIMAT(idum, idum))
772 161762 : do irow = 1, nrow
773 161762 : SOLMAT(irow, idum) = temp * SOLMAT(irow, idum)
774 : end do
775 : #endif
776 155549 : do icol = GET_HALF_RANGE(1, idum, ncol)
777 155549 : if (TRIMAT(icol, idum) /= ZERO) then
778 61803 : temp = GET_CONJG(TRIMAT(icol, idum))
779 466728 : do irow = 1, nrow
780 466728 : SOLMAT(irow, icol) = SOLMAT(irow, icol) - temp * SOLMAT(irow, idum)
781 : end do
782 : end if
783 : end do
784 84268 : if (ALPHA /= ONE) then
785 65456 : do irow = 1, nrow
786 65456 : SOLMAT(irow, idum) = ALPHA * SOLMAT(irow, idum)
787 : end do
788 : end if
789 : end do
790 : end block
791 : #else
792 : !%%%%%%%%%%%%%%%%%%%%%%%%
793 : #error "Unrecognized interface."
794 : !%%%%%%%%%%%%%%%%%%%%%%%%
795 : #endif
796 : #endif
797 :
798 : #else
799 : !%%%%%%%%%%%%%%%%%%%%%%%%
800 : !%%%%%%%%%%%%%%%%%%%%%%%%
801 : #error "Unrecognized interface."
802 : !%%%%%%%%%%%%%%%%%%%%%%%%
803 : !%%%%%%%%%%%%%%%%%%%%%%%%
804 : #endif
805 : #undef GET_FULL_RANGE
806 : #undef GET_HALF_RANGE
807 : #undef CXD_ENABLED
808 : #undef BLAS_TRANSA
809 : #undef BLAS_SIDE
810 : #undef BLAS_UPLO
811 : #undef BLAS_DIAG
812 : #undef TYPE_KIND
813 : #undef GET_CONJG
814 : #undef blasTRXX
815 : #undef SOLMAT
816 : #undef TRIMAT
817 : #undef ALPHA
|
