<T>LAPACK 0.1.2
C++ Template Linear Algebra PACKage
Loading...
Searching...
No Matches
trsm_tri.hpp
Go to the documentation of this file.
1
4//
5// Copyright (c) 2025, University of Colorado Denver. All rights reserved.
6//
7// This file is part of <T>LAPACK.
8// <T>LAPACK is free software: you can redistribute it and/or modify it under
9// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
10
11#ifndef TLAPACK_TRSM_TRI
12#define TLAPACK_TRSM_TRI
13
14#include "tlapack/base/utils.hpp"
15#include "tlapack/blas/trmm.hpp"
16#include "tlapack/lapack/lascl.hpp"
17#include "tlapack/lapack/trmm_out.hpp"
18
19namespace tlapack {
20
98template <typename matrixA_t, typename matrixB_t, TLAPACK_SCALAR alpha_t>
99void trsm_tri(Side sideA,
100 Uplo uploB,
101 Op transA,
102 Diag diagA,
103 const alpha_t& alpha,
104 const matrixA_t& A,
105 matrixB_t& B)
106{
107 using T = tlapack::type_t<matrixB_t>;
108 using idx_t = tlapack::size_type<matrixB_t>;
109 using range = std::pair<idx_t, idx_t>;
110 using real_t = tlapack::real_type<T>;
111
112 tlapack::Uplo uploA;
113 if (transA != tlapack::Op::NoTrans)
114 uploA = (uploB == tlapack::Uplo::Upper) ? tlapack::Uplo::Lower
115 : tlapack::Uplo::Upper;
116 else
117 uploA = uploB;
118
119 idx_t n = nrows(A);
120
121 // lascl(uploB, T(1.), alpha, B); this does not work
122
123 if (uploB == Uplo::Upper) {
124 for (idx_t j = 0; j < n; j++) {
125 for (idx_t i = 0; i < j + 1; i++) {
126 B(i, j) *= alpha;
127 }
128 }
129 }
130 else {
131 for (idx_t j = 0; j < n; j++) {
132 for (idx_t i = j; i < n; i++) {
133 B(i, j) *= alpha;
134 }
135 }
136 }
137
138 if (n == 1) {
139 if (diagA == Diag::NonUnit) {
140 if (transA == tlapack::Op::ConjTrans)
141 B(0, 0) /= conj(A(0, 0));
142 else
143 B(0, 0) /= A(0, 0);
144 return;
145 }
146 else {
147 return;
148 }
149 }
150
151 idx_t nd = n / 2;
152
153 auto A00 = slice(A, range(0, nd), range(0, nd));
154 auto A01 = slice(A, range(0, nd), range(nd, n));
155 auto A10 = slice(A, range(nd, n), range(0, nd));
156 auto A11 = slice(A, range(nd, n), range(nd, n));
157
158 auto B00 = slice(B, range(0, nd), range(0, nd));
159 auto B01 = slice(B, range(0, nd), range(nd, n));
160 auto B10 = slice(B, range(nd, n), range(0, nd));
161 auto B11 = slice(B, range(nd, n), range(nd, n));
162
163 if (sideA == tlapack::Side::Left) {
164 if (transA == tlapack::Op::NoTrans) {
165 // Form: B := alpha*inv(A)*B
166 if (uploB == tlapack::Uplo::Upper) {
167 // Left, NoTrans, Upper, Diag
168 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
169
170 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
171
172 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
173 real_t(-1), B11, A01, real_t(1), B01);
174
175 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B01);
176
177 return;
178 }
179 else {
180 // Left, NoTrans, Lower, Diag
181 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
182
183 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
184
185 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
186 real_t(-1), B00, A10, real_t(1), B10);
187
188 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B10);
189
190 return;
191 }
192 }
193 else {
194 // Form: B := alpha*inv(A**T)*B
195 // Form: B := alpha*inv(A**H)*B
196 if (uploB == tlapack::Uplo::Upper) {
197 // Left, Trans or ConjTrans, Upper, Diag
198
199 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
200
201 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
202
203 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
204 real_t(-1), B11, A10, real_t(1), B01);
205
206 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B01);
207
208 return;
209 }
210 else {
211 // Left, Trans or ConjTrans, Lower, Diag
212
213 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
214
215 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
216
217 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
218 real_t(-1), B00, A01, real_t(1), B10);
219
220 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B10);
221 return;
222 }
223 }
224 }
225 else {
226 if (transA == tlapack::Op::NoTrans) {
227 // Form: B := alpha*B*inv(A)
228 if (uploB == tlapack::Uplo::Upper) {
229 // Right, NoTrans, Upper, Diag
230
231 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
232
233 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
234
235 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
236 real_t(-1), B00, A01, real_t(1), B01);
237
238 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B01);
239
240 return;
241 }
242 else {
243 // Right, NoTrans, Lower, Diag
244 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
245
246 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
247
248 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
249 real_t(-1), B11, A10, real_t(1), B10);
250
251 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B10);
252
253 return;
254 }
255 }
256 else {
257 // Form: B := alpha*B*inv(A**T)
258 // Form: B := alpha*B*inv(A**H)
259
260 if (uploB == tlapack::Uplo::Upper) {
261 // Right, Trans or ConjTrans, Upper, Diag
262 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
263
264 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
265
266 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
267 real_t(-1), B00, A10, real_t(1), B01);
268
269 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B01);
270
271 return;
272 }
273 else {
274 // Right, Trans or ConjTrans, Lower, Diag
275
276 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
277
278 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
279
280 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
281 real_t(-1), B11, A01, real_t(1), B10);
282
283 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B10);
284 return;
285 }
286 }
287 }
288}
289} // namespace tlapack
290
291#endif // TLAPACK_TRSM_TRI
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::trsm_tri
void trsm_tri(Side sideA, Uplo uploB, Op transA, Diag diagA, const alpha_t &alpha, const matrixA_t &A, matrixB_t &B)
Solve the triangular matrix-matrix equation.
Definition trsm_tri.hpp:99
tlapack::trsm
void trsm(Side side, Uplo uplo, Op trans, Diag diag, const alpha_t &alpha, const matrixA_t &A, matrixB_t &B)
Solve the triangular matrix-vector equation.
Definition trsm.hpp:76
lascl.hpp
Multiplies a matrix by a scalar.
tlapack
Recursive QR factorization using compact WY Householder representation.
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::Side::Left
@ Left
left side
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.
tlapack::Uplo::Lower
@ Lower
0 <= i <= m, 0 <= j <= i.
trmm_out.hpp