<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
71template <typename matrixA_t, typename matrixB_t, TLAPACK_SCALAR alpha_t>
72void trsm_tri(Side sideA,
73 Uplo uploB,
74 Op transA,
75 Diag diagA,
76 const alpha_t& alpha,
77 const matrixA_t& A,
78 matrixB_t& B)
79{
80 using T = tlapack::type_t<matrixB_t>;
81 using idx_t = tlapack::size_type<matrixB_t>;
82 using range = std::pair<idx_t, idx_t>;
83 using real_t = tlapack::real_type<T>;
84
85 tlapack::Uplo uploA;
86 if (transA != tlapack::Op::NoTrans)
87 uploA = (uploB == tlapack::Uplo::Upper) ? tlapack::Uplo::Lower
88 : tlapack::Uplo::Upper;
89 else
90 uploA = uploB;
91
92 idx_t n = nrows(A);
93
94 // lascl(uploB, T(1.), alpha, B); this does not work
95
96 if (uploB == Uplo::Upper) {
97 for (idx_t j = 0; j < n; j++) {
98 for (idx_t i = 0; i < j + 1; i++) {
99 B(i, j) *= alpha;
100 }
101 }
102 }
103 else {
104 for (idx_t j = 0; j < n; j++) {
105 for (idx_t i = j; i < n; i++) {
106 B(i, j) *= alpha;
107 }
108 }
109 }
110
111 if (n == 1) {
112 if (diagA == Diag::NonUnit) {
113 if (transA == tlapack::Op::ConjTrans)
114 B(0, 0) /= conj(A(0, 0));
115 else
116 B(0, 0) /= A(0, 0);
117 return;
118 }
119 else {
120 // if (transA == tlapack::Op::ConjTrans)
121 // B(0, 0) /= conj(A(0, 0));
122 // else
123 // B(0, 0) /= A(0, 0);
124 return;
125 }
126 }
127
128 idx_t nd = n / 2;
129
130 auto A00 = slice(A, range(0, nd), range(0, nd));
131 auto A01 = slice(A, range(0, nd), range(nd, n));
132 auto A10 = slice(A, range(nd, n), range(0, nd));
133 auto A11 = slice(A, range(nd, n), range(nd, n));
134
135 auto B00 = slice(B, range(0, nd), range(0, nd));
136 auto B01 = slice(B, range(0, nd), range(nd, n));
137 auto B10 = slice(B, range(nd, n), range(0, nd));
138 auto B11 = slice(B, range(nd, n), range(nd, n));
139
140 if (sideA == tlapack::Side::Left) {
141 if (transA == tlapack::Op::NoTrans) {
142 // Form: B := alpha*inv(A)*B
143 if (uploB == tlapack::Uplo::Upper) {
144 // Left, NoTrans, Upper, Diag
145 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
146
147 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
148
149 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
150 real_t(-1), B11, A01, real_t(1), B01);
151
152 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B01);
153
154 return;
155 }
156 else {
157 // Left, NoTrans, Lower, Diag
158 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
159
160 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
161
162 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
163 real_t(-1), B00, A10, real_t(1), B10);
164
165 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B10);
166
167 return;
168 }
169 }
170 else {
171 // Form: B := alpha*inv(A**T)*B
172 // Form: B := alpha*inv(A**H)*B
173 if (uploB == tlapack::Uplo::Upper) {
174 // Left, Trans or ConjTrans, Upper, Diag
175
176 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
177
178 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
179
180 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
181 real_t(-1), B11, A10, real_t(1), B01);
182
183 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B01);
184
185 return;
186 }
187 else {
188 // Left, Trans or ConjTrans, Lower, Diag
189
190 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
191
192 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
193
194 trmm_out(Side::Right, uploB, Op::NoTrans, Diag::NonUnit, transA,
195 real_t(-1), B00, A01, real_t(1), B10);
196
197 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B10);
198 return;
199 }
200 }
201 }
202 else {
203 if (transA == tlapack::Op::NoTrans) {
204 // Form: B := alpha*B*inv(A)
205 if (uploB == tlapack::Uplo::Upper) {
206 // Right, NoTrans, Upper, Diag
207
208 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
209
210 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
211
212 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
213 real_t(-1), B00, A01, real_t(1), B01);
214
215 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B01);
216
217 return;
218 }
219 else {
220 // Right, NoTrans, Lower, Diag
221 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
222
223 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
224
225 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
226 real_t(-1), B11, A10, real_t(1), B10);
227
228 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B10);
229
230 return;
231 }
232 }
233 else {
234 // Form: B := alpha*B*inv(A**T)
235 // Form: B := alpha*B*inv(A**H)
236
237 if (uploB == tlapack::Uplo::Upper) {
238 // Right, Trans or ConjTrans, Upper, Diag
239 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
240
241 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
242
243 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
244 real_t(-1), B00, A10, real_t(1), B01);
245
246 trsm(sideA, uploA, transA, diagA, real_t(1), A11, B01);
247
248 return;
249 }
250 else {
251 // Right, Trans or ConjTrans, Lower, Diag
252
253 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A00, B00);
254
255 trsm_tri(sideA, uploB, transA, diagA, real_t(1), A11, B11);
256
257 trmm_out(Side::Left, uploB, Op::NoTrans, Diag::NonUnit, transA,
258 real_t(-1), B11, A01, real_t(1), B10);
259
260 trsm(sideA, uploA, transA, diagA, real_t(1), A00, B10);
261 return;
262 }
263 }
264 }
265}
266} // namespace tlapack
267
268#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-vector equation.
Definition trsm_tri.hpp:72
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
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::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