<T>LAPACK 0.1.2
C++ Template Linear Algebra PACKage
Loading...
Searching...
No Matches
trmm_out.hpp
Go to the documentation of this file.
1
3//
4// Copyright (c) 2025, University of Colorado Denver. All rights reserved.
5//
6// This file is part of <T>LAPACK.
7// <T>LAPACK is free software: you can redistribute it and/or modify it under
8// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
9
10#ifndef TLAPACK_TRMM_OUT
11#define TLAPACK_TRMM_OUT
12
13#include "../../../test/include/MatrixMarket.hpp"
14#include "tlapack/base/utils.hpp"
15#include "tlapack/blas/gemm.hpp"
16
17namespace tlapack {
76template <TLAPACK_MATRIX matrixA_t,
77 TLAPACK_MATRIX matrixB_t,
78 TLAPACK_MATRIX matrixC_t,
79 TLAPACK_SCALAR alpha_t,
80 TLAPACK_SCALAR beta_t,
81 class T = type_t<matrixB_t>>
82void trmm_out(Side side,
83 Uplo uplo,
84 Op transA,
85 Diag diag,
86 Op transB,
87 const alpha_t& alpha,
88 const matrixA_t& A,
89 const matrixB_t& B,
90 const beta_t& beta,
91 matrixC_t& C)
92{
93 using idx_t = tlapack::size_type<matrixA_t>;
94 using range = pair<idx_t, idx_t>;
95 using real_t = real_type<T>;
96 // only working with transA = notrans right now
97
98 idx_t m, n;
99 m = (transB == Op::NoTrans) ? nrows(B) : ncols(B);
100 n = (transB == Op::NoTrans) ? ncols(B) : nrows(B);
101
102 if (transB == Op::NoTrans) {
103 if (uplo == Uplo::Upper) {
104 if (diag == Diag::NonUnit) {
105 if (side == Side::Right) { // notransA, notransB, upper,
106 // nonunit, right works
107 idx_t n0 = n / 2;
108 if (n == 1) {
109 for (idx_t i = 0; i < m; ++i) {
110 C(i, 0) =
111 alpha * B(i, 0) * A(0, 0) + beta * C(i, 0);
112 }
113 }
114 else {
115 auto C0 = slice(C, range(0, m), range(0, n0));
116 auto C1 = slice(C, range(0, m), range(n0, n));
117
118 auto A00 = slice(A, range(0, n0), range(0, n0));
119 auto A01 = slice(A, range(0, n0), range(n0, n));
120 auto A11 = slice(A, range(n0, n), range(n0, n));
121
122 auto B0 = slice(B, range(0, m), range(0, n0));
123 auto B1 = slice(B, range(0, m), range(n0, n));
124
125 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
126 B0, beta, C0);
127
128 gemm(transB, transA, alpha, B0, A01, beta, C1);
129
130 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
131 B1, real_t(1), C1);
132 }
133 }
134 else { // notransA, notransB, upper, nonunit, left works
135
136 idx_t m0 = m / 2;
137
138 if (m == 1) {
139 for (idx_t j = 0; j < n; ++j)
140 C(0, j) =
141 alpha * B(0, j) * A(0, 0) + beta * C(0, j);
142 }
143 else {
144 auto C0 = slice(C, range(0, m0), range(0, n));
145 auto C1 = slice(C, range(m0, m), range(0, n));
146
147 auto A00 = slice(A, range(0, m0), range(0, m0));
148 auto A01 = slice(A, range(0, m0), range(m0, m));
149 auto A11 = slice(A, range(m0, m), range(m0, m));
150
151 auto B0 = slice(B, range(0, m0), range(0, n));
152 auto B1 = slice(B, range(m0, m), range(0, n));
153
154 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
155 B1, beta, C1);
156
157 gemm(transA, transB, alpha, A01, B1, beta, C0);
158
159 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
160 B0, real_t(1), C0);
161 }
162 }
163 }
164 else {
165 if (side == Side::Right) { // notransA, notransB, upper, unit,
166 // right works
167 idx_t n0 = n / 2;
168
169 if (n == 1) {
170 for (idx_t i = 0; i < m; ++i) {
171 C(i, 0) = alpha * B(i, 0) + beta * C(i, 0);
172 }
173 }
174 else {
175 auto C0 = slice(C, range(0, m), range(0, n0));
176 auto C1 = slice(C, range(0, m), range(n0, n));
177
178 auto A00 = slice(A, range(0, n0), range(0, n0));
179 auto A01 = slice(A, range(0, n0), range(n0, n));
180 auto A11 = slice(A, range(n0, n), range(n0, n));
181
182 auto B0 = slice(B, range(0, m), range(0, n0));
183 auto B1 = slice(B, range(0, m), range(n0, n));
184
185 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
186 B0, beta, C0);
187
188 gemm(transB, transA, alpha, B0, A01, beta, C1);
189
190 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
191 B1, real_t(1), C1);
192 }
193 }
194 else { // notransA, notransB, upper, unit, left works
195 idx_t m0 = m / 2;
196
197 if (m == 1) {
198 for (idx_t j = 0; j < n; ++j)
199 C(0, j) = alpha * B(0, j) + beta * C(0, j);
200 }
201 else {
202 auto C0 = slice(C, range(0, m0), range(0, n));
203 auto C1 = slice(C, range(m0, m), range(0, n));
204
205 auto A00 = slice(A, range(0, m0), range(0, m0));
206 auto A01 = slice(A, range(0, m0), range(m0, m));
207 auto A11 = slice(A, range(m0, m), range(m0, m));
208
209 auto B0 = slice(B, range(0, m0), range(0, n));
210 auto B1 = slice(B, range(m0, m), range(0, n));
211
212 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
213 B1, beta, C1);
214
215 gemm(transA, transB, alpha, A01, B1, beta, C0);
216
217 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
218 B0, real_t(1), C0);
219 }
220 }
221 }
222 }
223 else {
224 if (diag == Diag::NonUnit) {
225 if (side == Side::Right) { // notransA, notransB, right,
226 // nonunit and lower works
227
228 idx_t n0 = n / 2;
229 if (n == 1) {
230 for (idx_t i = 0; i < m; ++i) {
231 C(i, 0) =
232 alpha * B(i, 0) * A(0, 0) + beta * C(i, 0);
233 }
234 }
235 else {
236 auto C0 = slice(C, range(0, m), range(0, n0));
237 auto C1 = slice(C, range(0, m), range(n0, n));
238
239 auto A00 = slice(A, range(0, n0), range(0, n0));
240 auto A10 = slice(A, range(n0, n), range(0, n0));
241 auto A11 = slice(A, range(n0, n), range(n0, n));
242
243 auto B0 = slice(B, range(0, m), range(0, n0));
244 auto B1 = slice(B, range(0, m), range(n0, n));
245
246 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
247 B1, beta, C1);
248
249 gemm(transB, transA, alpha, B1, A10, beta, C0);
250
251 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
252 B0, real_t(1), C0);
253 }
254 }
255 else { // notransA, notransB, left, nonunit, lower works
256 idx_t m0 = m / 2;
257
258 if (m == 1) {
259 for (idx_t j = 0; j < n; ++j)
260 C(0, j) =
261 alpha * B(0, j) * A(0, 0) + beta * C(0, j);
262 }
263 else {
264 auto C0 = slice(C, range(0, m0), range(0, n));
265 auto C1 = slice(C, range(m0, m), range(0, n));
266
267 auto A00 = slice(A, range(0, m0), range(0, m0));
268 auto A10 = slice(A, range(m0, m), range(0, m0));
269 auto A11 = slice(A, range(m0, m), range(m0, m));
270
271 auto B0 = slice(B, range(0, m0), range(0, n));
272 auto B1 = slice(B, range(m0, m), range(0, n));
273
274 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
275 B0, beta, C0);
276
277 gemm(transA, transB, alpha, A10, B0, beta, C1);
278
279 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
280 B1, real_t(1), C1);
281 }
282 }
283 }
284 else {
285 if (side == Side::Right) { // notransA, notransB, right, unit
286 // and lower works
287
288 idx_t n0 = n / 2;
289 if (n == 1) {
290 for (idx_t i = 0; i < m; ++i) {
291 C(i, 0) = alpha * B(i, 0) + beta * C(i, 0);
292 }
293 }
294 else {
295 auto C0 = slice(C, range(0, m), range(0, n0));
296 auto C1 = slice(C, range(0, m), range(n0, n));
297
298 auto A00 = slice(A, range(0, n0), range(0, n0));
299 auto A10 = slice(A, range(n0, n), range(0, n0));
300 auto A11 = slice(A, range(n0, n), range(n0, n));
301
302 auto B0 = slice(B, range(0, m), range(0, n0));
303 auto B1 = slice(B, range(0, m), range(n0, n));
304
305 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
306 B1, beta, C1);
307
308 gemm(transB, transA, alpha, B1, A10, beta, C0);
309
310 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
311 B0, real_t(1), C0);
312 }
313 }
314 else { // notransA, notransB, left, unit and lower works
315 idx_t m0 = m / 2;
316
317 if (m == 1) {
318 for (idx_t j = 0; j < n; ++j)
319 C(0, j) = alpha * B(0, j) + beta * C(0, j);
320 }
321 else {
322 auto C0 = slice(C, range(0, m0), range(0, n));
323 auto C1 = slice(C, range(m0, m), range(0, n));
324
325 auto A00 = slice(A, range(0, m0), range(0, m0));
326 auto A10 = slice(A, range(m0, m), range(0, m0));
327 auto A11 = slice(A, range(m0, m), range(m0, m));
328
329 auto B0 = slice(B, range(0, m0), range(0, n));
330 auto B1 = slice(B, range(m0, m), range(0, n));
331
332 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
333 B0, beta, C0);
334
335 gemm(transA, transB, alpha, A10, B0, beta, C1);
336
337 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
338 B1, real_t(1), C1);
339 }
340 }
341 }
342 }
343 }
344 else if (transB == Op::ConjTrans) {
345 if (uplo == Uplo::Upper) {
346 if (diag == Diag::NonUnit) {
347 if (side == Side::Right) { // notransA, conjtransB, nonunit,
348 // right, upper works
349
350 idx_t n0 = n / 2;
351 if (n == 1) {
352 for (idx_t j = 0; j < m; ++j) {
353 C(j, 0) = alpha * conj(B(0, j)) * A(0, 0) +
354 beta * C(j, 0);
355 }
356 }
357 else {
358 auto C0 = slice(C, range(0, m), range(0, n0));
359 auto C1 = slice(C, range(0, m), range(n0, n));
360
361 auto A00 = slice(A, range(0, n0), range(0, n0));
362 auto A01 = slice(A, range(0, n0), range(n0, n));
363 auto A11 = slice(A, range(n0, n), range(n0, n));
364
365 auto B0 = slice(B, range(0, n0), range(0, m));
366 auto B1 = slice(B, range(n0, n), range(0, m));
367
368 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
369 B0, beta, C0);
370
371 gemm(transB, transA, alpha, B0, A01, beta, C1);
372
373 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
374 B1, real_t(1), C1);
375 }
376 }
377 else { // notransA, conjtransB, nonunit, left, upper works
378 idx_t m0 = m / 2;
379
380 if (m == 1) {
381 for (idx_t j = 0; j < n; ++j)
382 C(0, j) = alpha * conj(B(j, 0)) * A(0, 0) +
383 beta * C(0, j);
384 }
385 else {
386 auto C0 = slice(C, range(0, m0), range(0, n));
387 auto C1 = slice(C, range(m0, m), range(0, n));
388
389 auto A00 = slice(A, range(0, m0), range(0, m0));
390 auto A01 = slice(A, range(0, m0), range(m0, m));
391 auto A11 = slice(A, range(m0, m), range(m0, m));
392
393 auto B0 = slice(B, range(0, n), range(0, m0));
394 auto B1 = slice(B, range(0, n), range(m0, m));
395
396 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
397 B1, beta, C1);
398
399 gemm(transA, transB, alpha, A01, B1, beta, C0);
400
401 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
402 B0, real_t(1), C0);
403 }
404 }
405 }
406 else {
407 if (side == Side::Right) { // notransA, conjtransB, unit,
408 // right, upper works
409
410 idx_t n0 = n / 2;
411 if (n == 1) {
412 for (idx_t j = 0; j < m; ++j) {
413 C(j, 0) = alpha * conj(B(0, j)) + beta * C(j, 0);
414 }
415 }
416 else {
417 auto C0 = slice(C, range(0, m), range(0, n0));
418 auto C1 = slice(C, range(0, m), range(n0, n));
419
420 auto A00 = slice(A, range(0, n0), range(0, n0));
421 auto A01 = slice(A, range(0, n0), range(n0, n));
422 auto A11 = slice(A, range(n0, n), range(n0, n));
423
424 auto B0 = slice(B, range(0, n0), range(0, m));
425 auto B1 = slice(B, range(n0, n), range(0, m));
426
427 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
428 B0, beta, C0);
429
430 gemm(transB, transA, alpha, B0, A01, beta, C1);
431
432 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
433 B1, real_t(1), C1);
434 }
435 }
436 else { // notransA, conjtransB, unit, left, upper works
437 idx_t m0 = m / 2;
438
439 if (m == 1) {
440 for (idx_t j = 0; j < n; ++j)
441 C(0, j) = alpha * conj(B(j, 0)) + beta * C(0, j);
442 }
443 else {
444 auto C0 = slice(C, range(0, m0), range(0, n));
445 auto C1 = slice(C, range(m0, m), range(0, n));
446
447 auto A00 = slice(A, range(0, m0), range(0, m0));
448 auto A01 = slice(A, range(0, m0), range(m0, m));
449 auto A11 = slice(A, range(m0, m), range(m0, m));
450
451 auto B0 = slice(B, range(0, n), range(0, m0));
452 auto B1 = slice(B, range(0, n), range(m0, m));
453
454 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
455 B1, beta, C1);
456
457 gemm(transA, transB, alpha, A01, B1, beta, C0);
458
459 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
460 B0, real_t(1), C0);
461 }
462 }
463 }
464 }
465 else {
466 if (diag == Diag::NonUnit) {
467 if (side == Side::Right) { // notransA, conjtransB, nonunit,
468 // right, lower works
469
470 idx_t n0 = n / 2;
471 if (n == 1) {
472 for (idx_t j = 0; j < m; ++j) {
473 C(j, 0) = alpha * conj(B(0, j)) * A(0, 0) +
474 beta * C(j, 0);
475 }
476 }
477 else {
478 auto C0 = slice(C, range(0, m), range(0, n0));
479 auto C1 = slice(C, range(0, m), range(n0, n));
480
481 auto A00 = slice(A, range(0, n0), range(0, n0));
482 auto A10 = slice(A, range(n0, n), range(0, n0));
483 auto A11 = slice(A, range(n0, n), range(n0, n));
484
485 auto B0 = slice(B, range(0, n0), range(0, m));
486 auto B1 = slice(B, range(n0, n), range(0, m));
487
488 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
489 B1, beta, C1);
490
491 gemm(transB, transA, alpha, B1, A10, beta, C0);
492
493 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
494 B0, real_t(1), C0);
495 }
496 }
497 else { // notransA, conjtransB, nonunit, left, lower works
498 idx_t m0 = m / 2;
499
500 if (m == 1) {
501 for (idx_t j = 0; j < n; ++j)
502 C(0, j) = alpha * conj(B(j, 0)) * A(0, 0) +
503 beta * C(0, j);
504 }
505 else {
506 auto C0 = slice(C, range(0, m0), range(0, n));
507 auto C1 = slice(C, range(m0, m), range(0, n));
508
509 auto A00 = slice(A, range(0, m0), range(0, m0));
510 auto A10 = slice(A, range(m0, m), range(0, m0));
511 auto A11 = slice(A, range(m0, m), range(m0, m));
512
513 auto B0 = slice(B, range(0, n), range(0, m0));
514 auto B1 = slice(B, range(0, n), range(m0, m));
515
516 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
517 B0, beta, C0);
518
519 gemm(transA, transB, alpha, A10, B0, beta, C1);
520
521 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
522 B1, real_t(1), C1);
523 }
524 }
525 }
526 else {
527 if (side == Side::Right) { // notransA, conjtransB, unit,
528 // right, lower works
529
530 idx_t n0 = n / 2;
531 if (n == 1) {
532 for (idx_t j = 0; j < m; ++j) {
533 C(j, 0) = alpha * conj(B(0, j)) + beta * C(j, 0);
534 }
535 }
536 else {
537 auto C0 = slice(C, range(0, m), range(0, n0));
538 auto C1 = slice(C, range(0, m), range(n0, n));
539
540 auto A00 = slice(A, range(0, n0), range(0, n0));
541 auto A10 = slice(A, range(n0, n), range(0, n0));
542 auto A11 = slice(A, range(n0, n), range(n0, n));
543
544 auto B0 = slice(B, range(0, n0), range(0, m));
545 auto B1 = slice(B, range(n0, n), range(0, m));
546
547 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
548 B1, beta, C1);
549
550 gemm(transB, transA, alpha, B1, A10, beta, C0);
551
552 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
553 B0, real_t(1), C0);
554 }
555 }
556 else { // notransA, conjtransB, unit, left, lower
557 idx_t m0 = m / 2;
558
559 if (m == 1) {
560 for (idx_t j = 0; j < n; ++j)
561 C(0, j) = alpha * conj(B(j, 0)) + beta * C(0, j);
562 }
563 else {
564 auto C0 = slice(C, range(0, m0), range(0, n));
565 auto C1 = slice(C, range(m0, m), range(0, n));
566
567 auto A00 = slice(A, range(0, m0), range(0, m0));
568 auto A10 = slice(A, range(m0, m), range(0, m0));
569 auto A11 = slice(A, range(m0, m), range(m0, m));
570
571 auto B0 = slice(B, range(0, n), range(0, m0));
572 auto B1 = slice(B, range(0, n), range(m0, m));
573
574 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
575 B0, beta, C0);
576
577 gemm(transA, transB, alpha, A10, B0, beta, C1);
578
579 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
580 B1, real_t(1), C1);
581 }
582 }
583 }
584 }
585 }
586 else {
587 if (uplo == Uplo::Upper) {
588 if (diag == Diag::NonUnit) {
589 if (side == Side::Right) { // notransA, transB, nonunit, right,
590 // upper works
591
592 idx_t n0 = n / 2;
593 if (n == 1) {
594 for (idx_t j = 0; j < m; ++j) {
595 C(j, 0) =
596 alpha * B(0, j) * A(0, 0) + beta * C(j, 0);
597 }
598 }
599 else {
600 auto C0 = slice(C, range(0, m), range(0, n0));
601 auto C1 = slice(C, range(0, m), range(n0, n));
602
603 auto A00 = slice(A, range(0, n0), range(0, n0));
604 auto A01 = slice(A, range(0, n0), range(n0, n));
605 auto A11 = slice(A, range(n0, n), range(n0, n));
606
607 auto B0 = slice(B, range(0, n0), range(0, m));
608 auto B1 = slice(B, range(n0, n), range(0, m));
609
610 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
611 B0, beta, C0);
612
613 gemm(transB, transA, alpha, B0, A01, beta, C1);
614
615 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
616 B1, real_t(1), C1);
617 }
618 }
619 else { // notransA, transB, nonunit, left, upper
620 idx_t m0 = m / 2;
621
622 if (m == 1) {
623 for (idx_t j = 0; j < n; ++j)
624 C(0, j) =
625 alpha * B(j, 0) * A(0, 0) + beta * C(0, j);
626 }
627 else {
628 auto C0 = slice(C, range(0, m0), range(0, n));
629 auto C1 = slice(C, range(m0, m), range(0, n));
630
631 auto A00 = slice(A, range(0, m0), range(0, m0));
632 auto A01 = slice(A, range(0, m0), range(m0, m));
633 auto A11 = slice(A, range(m0, m), range(m0, m));
634
635 auto B0 = slice(B, range(0, n), range(0, m0));
636 auto B1 = slice(B, range(0, n), range(m0, m));
637
638 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
639 B1, beta, C1);
640
641 gemm(transA, transB, alpha, A01, B1, beta, C0);
642
643 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
644 B0, real_t(1), C0);
645 }
646 }
647 }
648 else {
649 if (side == Side::Right) { // notransA, transB, unit, right,
650 // upper works
651
652 idx_t n0 = n / 2;
653 if (n == 1) {
654 for (idx_t j = 0; j < m; ++j) {
655 C(j, 0) = alpha * B(0, j) + beta * C(j, 0);
656 }
657 }
658 else {
659 auto C0 = slice(C, range(0, m), range(0, n0));
660 auto C1 = slice(C, range(0, m), range(n0, n));
661
662 auto A00 = slice(A, range(0, n0), range(0, n0));
663 auto A01 = slice(A, range(0, n0), range(n0, n));
664 auto A11 = slice(A, range(n0, n), range(n0, n));
665
666 auto B0 = slice(B, range(0, n0), range(0, m));
667 auto B1 = slice(B, range(n0, n), range(0, m));
668
669 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
670 B0, beta, C0);
671
672 gemm(transB, transA, alpha, B0, A01, beta, C1);
673
674 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
675 B1, real_t(1), C1);
676 }
677 }
678 else { // notransA, transB, unit, left, upper
679 idx_t m0 = m / 2;
680
681 if (m == 1) {
682 for (idx_t j = 0; j < n; ++j)
683 C(0, j) = alpha * B(j, 0) + beta * C(0, j);
684 }
685 else {
686 auto C0 = slice(C, range(0, m0), range(0, n));
687 auto C1 = slice(C, range(m0, m), range(0, n));
688
689 auto A00 = slice(A, range(0, m0), range(0, m0));
690 auto A01 = slice(A, range(0, m0), range(m0, m));
691 auto A11 = slice(A, range(m0, m), range(m0, m));
692
693 auto B0 = slice(B, range(0, n), range(0, m0));
694 auto B1 = slice(B, range(0, n), range(m0, m));
695
696 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
697 B1, beta, C1);
698
699 gemm(transA, transB, alpha, A01, B1, beta, C0);
700
701 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
702 B0, real_t(1), C0);
703 }
704 }
705 }
706 }
707 else {
708 if (diag == Diag::NonUnit) {
709 if (side == Side::Right) { // notransA, transB, nonunit, right,
710 // lower works
711
712 idx_t n0 = n / 2;
713 if (n == 1) {
714 for (idx_t j = 0; j < m; ++j) {
715 C(j, 0) =
716 alpha * B(0, j) * A(0, 0) + beta * C(j, 0);
717 }
718 }
719 else {
720 auto C0 = slice(C, range(0, m), range(0, n0));
721 auto C1 = slice(C, range(0, m), range(n0, n));
722
723 auto A00 = slice(A, range(0, n0), range(0, n0));
724 auto A10 = slice(A, range(n0, n), range(0, n0));
725 auto A11 = slice(A, range(n0, n), range(n0, n));
726
727 auto B0 = slice(B, range(0, n0), range(0, m));
728 auto B1 = slice(B, range(n0, n), range(0, m));
729
730 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
731 B1, beta, C1);
732
733 gemm(transB, transA, alpha, B1, A10, beta, C0);
734
735 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
736 B0, real_t(1), C0);
737 }
738 }
739 else { // notransA, transB, nonunit, left, lower
740 idx_t m0 = m / 2;
741
742 if (m == 1) {
743 for (idx_t j = 0; j < n; ++j)
744 C(0, j) =
745 alpha * B(j, 0) * A(0, 0) + beta * C(0, j);
746 }
747 else {
748 auto C0 = slice(C, range(0, m0), range(0, n));
749 auto C1 = slice(C, range(m0, m), range(0, n));
750
751 auto A00 = slice(A, range(0, m0), range(0, m0));
752 auto A10 = slice(A, range(m0, m), range(0, m0));
753 auto A11 = slice(A, range(m0, m), range(m0, m));
754
755 auto B0 = slice(B, range(0, n), range(0, m0));
756 auto B1 = slice(B, range(0, n), range(m0, m));
757
758 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
759 B0, beta, C0);
760
761 gemm(transA, transB, alpha, A10, B0, beta, C1);
762
763 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
764 B1, real_t(1), C1);
765 }
766 }
767 }
768 else {
769 if (side == Side::Right) { // notransA, transB, unit, right,
770 // lower works
771
772 idx_t n0 = n / 2;
773 if (n == 1) {
774 for (idx_t j = 0; j < m; ++j) {
775 C(j, 0) = alpha * B(0, j) + beta * C(j, 0);
776 }
777 }
778 else {
779 auto C0 = slice(C, range(0, m), range(0, n0));
780 auto C1 = slice(C, range(0, m), range(n0, n));
781
782 auto A00 = slice(A, range(0, n0), range(0, n0));
783 auto A10 = slice(A, range(n0, n), range(0, n0));
784 auto A11 = slice(A, range(n0, n), range(n0, n));
785
786 auto B0 = slice(B, range(0, n0), range(0, m));
787 auto B1 = slice(B, range(n0, n), range(0, m));
788
789 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
790 B1, beta, C1);
791
792 gemm(transB, transA, alpha, B1, A10, beta, C0);
793
794 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
795 B0, real_t(1), C0);
796 }
797 }
798 else { // notransA, transB, unit, left, lower
799 idx_t m0 = m / 2;
800
801 if (m == 1) {
802 for (idx_t j = 0; j < n; ++j)
803 C(0, j) = alpha * B(j, 0) + beta * C(0, j);
804 }
805 else {
806 auto C0 = slice(C, range(0, m0), range(0, n));
807 auto C1 = slice(C, range(m0, m), range(0, n));
808
809 auto A00 = slice(A, range(0, m0), range(0, m0));
810 auto A10 = slice(A, range(m0, m), range(0, m0));
811 auto A11 = slice(A, range(m0, m), range(m0, m));
812
813 auto B0 = slice(B, range(0, n), range(0, m0));
814 auto B1 = slice(B, range(0, n), range(m0, m));
815
816 trmm_out(side, uplo, transA, diag, transB, alpha, A00,
817 B0, beta, C0);
818
819 gemm(transA, transB, alpha, A10, B0, beta, C1);
820
821 trmm_out(side, uplo, transA, diag, transB, alpha, A11,
822 B1, real_t(1), C1);
823 }
824 }
825 }
826 }
827 }
828}
829
830} // namespace tlapack
831
832#endif // TLAPACK_TRMM_OUT
TLAPACK_SCALAR
#define TLAPACK_SCALAR
Macro for tlapack::concepts::Scalar compatible with C++17.
Definition concepts.hpp:915
TLAPACK_MATRIX
#define TLAPACK_MATRIX
Macro for tlapack::concepts::Matrix compatible with C++17.
Definition concepts.hpp:896
tlapack::trmm_out
void trmm_out(Side side, Uplo uplo, Op transA, Diag diag, Op transB, const alpha_t &alpha, const matrixA_t &A, const matrixB_t &B, const beta_t &beta, matrixC_t &C)
Triangular matrix-matrix multiply:
Definition trmm_out.hpp:82
tlapack::gemm
void gemm(Op transA, Op transB, const alpha_t &alpha, const matrixA_t &A, const matrixB_t &B, const beta_t &beta, matrixC_t &C)
General matrix-matrix multiply:
Definition gemm.hpp:61
tlapack
Sort the numbers in D in increasing order (if ID = 'I') or in decreasing order (if ID = 'D' ).
Definition arrayTraits.hpp:15
tlapack::real_type
typename traits::real_type_traits< Types..., int >::type real_type
The common real type of the list of types.
Definition scalar_type_traits.hpp:113
tlapack::conj
constexpr T conj(const T &x) noexcept
Extends std::conj() to real datatypes.
Definition utils.hpp:100
tlapack::Diag
Diag
Definition types.hpp:197
tlapack::Diag::NonUnit
@ NonUnit
The main diagonal is not assumed to consist of 1's.
tlapack::Side
Side
Definition types.hpp:271
tlapack::Side::Right
@ Right
right side
tlapack::diag
constexpr auto diag(T &A, int diagIdx=0) noexcept
Get the Diagonal of an Eigen Matrix.
Definition eigen.hpp:576
tlapack::Op
Op
Definition types.hpp:227
tlapack::Op::NoTrans
@ NoTrans
no transpose
tlapack::Op::ConjTrans
@ ConjTrans
conjugate transpose
tlapack::Uplo
Uplo
Definition types.hpp:50
tlapack::Uplo::Upper
@ Upper
0 <= i <= j, 0 <= j <= n.