<T>LAPACK 0.1.2
C++ Template Linear Algebra PACKage
Loading...
Searching...
No Matches
generalized_schur_swap.hpp
Go to the documentation of this file.
1
5//
6// Copyright (c) 2025, University of Colorado Denver. All rights reserved.
7//
8// This file is part of <T>LAPACK.
9// <T>LAPACK is free software: you can redistribute it and/or modify it under
10// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
11
12#ifndef TLAPACK_GENERALIZED_SCHUR_SWAP_HH
13#define TLAPACK_GENERALIZED_SCHUR_SWAP_HH
14
15#include "tlapack/base/utils.hpp"
16#include "tlapack/blas/gemm.hpp"
17#include "tlapack/blas/rot.hpp"
18#include "tlapack/blas/rotg.hpp"
19#include "tlapack/blas/swap.hpp"
20#include "tlapack/blas/trsv.hpp"
21#include "tlapack/lapack/getrf.hpp"
22#include "tlapack/lapack/inv_house3.hpp"
23#include "tlapack/lapack/lahqz_eig22.hpp"
24#include "tlapack/lapack/lahqz_schur22.hpp"
25#include "tlapack/lapack/lange.hpp"
26#include "tlapack/lapack/larfg.hpp"
27#include "tlapack/lapack/svd22.hpp"
28
29namespace tlapack {
30
58template <TLAPACK_CSMATRIX matrix_t,
59 enable_if_t<is_real<type_t<matrix_t>>, bool> = true>
60int generalized_schur_swap(bool want_q,
61 bool want_z,
62 matrix_t& A,
63 matrix_t& B,
64 matrix_t& Q,
65 matrix_t& Z,
66 const size_type<matrix_t>& j0,
67 const size_type<matrix_t>& n1,
68 const size_type<matrix_t>& n2)
69{
70 using idx_t = size_type<matrix_t>;
71 using T = type_t<matrix_t>;
72 using range = pair<idx_t, idx_t>;
73
74 // Functor for creating new matrices
75 Create<matrix_t> new_matrix;
76
77 const idx_t n = ncols(A);
78 const T zero(0);
79
80 tlapack_check(nrows(A) == n);
81 tlapack_check(nrows(Q) == n);
82 tlapack_check(ncols(Q) == n);
83 tlapack_check(0 <= j0);
84 tlapack_check(j0 + n1 + n2 <= n);
85 tlapack_check(n1 == 1 or n1 == 2);
86 tlapack_check(n2 == 1 or n2 == 2);
87
88 const idx_t j1 = j0 + 1;
89 const idx_t j2 = j0 + 2;
90 const idx_t j3 = j0 + 3;
91
92 // Check if the 2x2 eigenvalue blocks consist of 2 1x1 blocks
93 // If so, treat them separately
94 if (n1 == 2)
95 if (A(j1, j0) == zero) {
96 int info;
97 info = generalized_schur_swap(want_q, want_z, A, B, Q, Z, j1,
98 (idx_t)1, n2);
99 if (info != 0) return info;
100 info = generalized_schur_swap(want_q, want_z, A, B, Q, Z, j0,
101 (idx_t)1, n2);
102 if (info != 0) return info;
103 return 0;
104 }
105 if (n2 == 2)
106 if (A(j0 + n1 + 1, j0 + n1) == zero) {
107 int info;
108 info = generalized_schur_swap(want_q, want_z, A, B, Q, Z, j0, n1,
109 (idx_t)1);
110 if (info != 0) return info;
111 info = generalized_schur_swap(want_q, want_z, A, B, Q, Z, j1, n1,
112 (idx_t)1);
113 if (info != 0) return info;
114 return 0;
115 }
116
117 if (n1 == 1 and n2 == 1) {
118 //
119 // Swap two 1-by-1 blocks.
120 //
121 const T a00 = A(j0, j0);
122 const T a01 = A(j0, j1);
123 const T a11 = A(j1, j1);
124 const T b00 = B(j0, j0);
125 const T b01 = B(j0, j1);
126 const T b11 = B(j1, j1);
127
128 const bool use_b = abs(b11 * a00) > abs(b00 * a11);
129
130 //
131 // Determine the transformation to perform the interchange
132 //
133 T cl, sl, cr, sr;
134 T temp = b11 * a00 - a11 * b00;
135 T temp2 = b11 * a01 - a11 * b01;
136 rotg(temp2, temp, cr, sr);
137
138 // Apply transformation from the right
139 {
140 auto a1 = slice(A, range{0, j1 + 1}, j0);
141 auto a2 = slice(A, range{0, j1 + 1}, j1);
142 rot(a2, a1, cr, sr);
143 auto b1 = slice(B, range{0, j1 + 1}, j0);
144 auto b2 = slice(B, range{0, j1 + 1}, j1);
145 rot(b2, b1, cr, sr);
146 if (want_z) {
147 auto z1 = col(Z, j0);
148 auto z2 = col(Z, j1);
149 rot(z2, z1, cr, sr);
150 }
151 }
152
153 if (use_b) {
154 temp = B(j0, j0);
155 temp2 = B(j1, j0);
156 }
157 else {
158 temp = A(j0, j0);
159 temp2 = A(j1, j0);
160 }
161 rotg(temp, temp2, cl, sl);
162
163 // Apply transformation from the left
164 {
165 auto a1 = slice(A, j0, range{j0, n});
166 auto a2 = slice(A, j1, range{j0, n});
167 rot(a1, a2, cl, sl);
168 auto b1 = slice(B, j0, range{j0, n});
169 auto b2 = slice(B, j1, range{j0, n});
170 rot(b1, b2, cl, sl);
171 if (want_q) {
172 auto q1 = col(Q, j0);
173 auto q2 = col(Q, j1);
174 rot(q1, q2, cl, sl);
175 }
176 }
177
178 A(j1, j0) = (T)0;
179 B(j1, j0) = (T)0;
180 }
181 if (n1 == 1 and n2 == 2) {
182 //
183 // Swap 1-by-1 block with 2-by-2 block
184 //
185
186 std::vector<T> H_;
187 auto H = new_matrix(H_, 2, 3);
188 std::vector<T> vl(3);
189 std::vector<T> vrA(3);
190 std::vector<T> vrB(3);
191 T taul, taurA, taurB;
192
193 std::vector<T> AA_;
194 auto AA = new_matrix(AA_, 3, 3);
195 std::vector<T> BB_;
196 auto BB = new_matrix(BB_, 3, 3);
197 lacpy(GENERAL, slice(A, range(j0, j3), range(j0, j3)), AA);
198 lacpy(GENERAL, slice(B, range(j0, j3), range(j0, j3)), BB);
199
200 auto a00 = AA(0, 0);
201 auto b00 = BB(0, 0);
202
203 T norma = lange(FROB_NORM, AA);
204 T normb = lange(FROB_NORM, BB);
205
206 H(0, 0) = b00 * AA(2, 1) - a00 * BB(2, 1);
207 H(0, 1) = b00 * AA(1, 1) - a00 * BB(1, 1);
208 H(0, 2) = b00 * AA(0, 1) - a00 * BB(0, 1);
209 H(1, 0) = b00 * AA(2, 2) - a00 * BB(2, 2);
210 H(1, 1) = b00 * AA(1, 2) - a00 * BB(1, 2);
211 H(1, 2) = b00 * AA(0, 2) - a00 * BB(0, 2);
212
213 inv_house3(H, vl, taul);
214
215 // Tentatively apply update from the left to local matrices
216 for (idx_t j = 0; j < 3; ++j) {
217 T sum = AA(2, j) + vl[1] * AA(1, j) + vl[2] * AA(0, j);
218 AA(2, j) = AA(2, j) - sum * taul;
219 AA(1, j) = AA(1, j) - sum * taul * vl[1];
220 AA(0, j) = AA(0, j) - sum * taul * vl[2];
221 }
222 for (idx_t j = 0; j < 3; ++j) {
223 T sum = BB(2, j) + vl[1] * BB(1, j) + vl[2] * BB(0, j);
224 BB(2, j) = BB(2, j) - sum * taul;
225 BB(1, j) = BB(1, j) - sum * taul * vl[1];
226 BB(0, j) = BB(0, j) - sum * taul * vl[2];
227 }
228
229 // Determine two sets of right transformations
230 // one that zeroes out the last row of A and one that zeroes out the
231 // last row of B
232 // We will later choose the one that gives the smaller error
233 vrA[0] = AA(2, 2);
234 vrA[1] = AA(2, 1);
235 vrA[2] = AA(2, 0);
236 larfg(FORWARD, COLUMNWISE_STORAGE, vrA, taurA);
237
238 vrB[0] = BB(2, 2);
239 vrB[1] = BB(2, 1);
240 vrB[2] = BB(2, 0);
241 larfg(FORWARD, COLUMNWISE_STORAGE, vrB, taurB);
242
243 // Apply the update calculated using BB to AA
244 // and vice versa. We choose the one that gives the smaller error
245 for (idx_t j = 2; j < 3; ++j) {
246 T sum = AA(j, 2) + vrB[1] * AA(j, 1) + vrB[2] * AA(j, 0);
247 AA(j, 2) = AA(j, 2) - sum * taurB;
248 AA(j, 1) = AA(j, 1) - sum * taurB * vrB[1];
249 AA(j, 0) = AA(j, 0) - sum * taurB * vrB[2];
250 }
251 for (idx_t j = 2; j < 3; ++j) {
252 T sum = BB(j, 2) + vrA[1] * BB(j, 1) + vrA[2] * BB(j, 0);
253 BB(j, 2) = BB(j, 2) - sum * taurA;
254 BB(j, 1) = BB(j, 1) - sum * taurA * vrA[1];
255 BB(j, 0) = BB(j, 0) - sum * taurA * vrA[2];
256 }
257
258 //
259 // Determine if the swap was successful
260 //
261 T errA = lapy2(AA(2, 0), AA(2, 1));
262 T errB = lapy2(BB(2, 0), BB(2, 1));
263 const T eps = ulp<T>();
264 const T small_num = safe_min<T>();
265
266 if (errA > max((T)20 * norma * eps, small_num) and
267 errB > max((T)20 * normb * eps, small_num)) {
268 // The swap failed, return with error
269 // Note, though we don't have a proof that this will always be the
270 // case, there are currently no known cases where this swap can
271 // fail.
272 return 1;
273 }
274
275 //
276 // Swap is accepted, apply the updates to the original matrices
277 //
278 for (idx_t j = j0; j < n; ++j) {
279 T sum = A(j2, j) + vl[1] * A(j1, j) + vl[2] * A(j0, j);
280 A(j2, j) = A(j2, j) - sum * taul;
281 A(j1, j) = A(j1, j) - sum * taul * vl[1];
282 A(j0, j) = A(j0, j) - sum * taul * vl[2];
283 }
284 for (idx_t j = j0; j < n; ++j) {
285 T sum = B(j2, j) + vl[1] * B(j1, j) + vl[2] * B(j0, j);
286 B(j2, j) = B(j2, j) - sum * taul;
287 B(j1, j) = B(j1, j) - sum * taul * vl[1];
288 B(j0, j) = B(j0, j) - sum * taul * vl[2];
289 }
290 if (want_q) {
291 for (idx_t j = 0; j < n; ++j) {
292 T sum = Q(j, j2) + vl[1] * Q(j, j1) + vl[2] * Q(j, j0);
293 Q(j, j2) = Q(j, j2) - sum * taul;
294 Q(j, j1) = Q(j, j1) - sum * taul * vl[1];
295 Q(j, j0) = Q(j, j0) - sum * taul * vl[2];
296 }
297 }
298 if (errB * norma < errA * normb) {
299 // The error is smallest when using vrA
300 for (idx_t j = 0; j < j3; ++j) {
301 T sum = A(j, j2) + vrA[1] * A(j, j1) + vrA[2] * A(j, j0);
302 A(j, j2) = A(j, j2) - sum * taurA;
303 A(j, j1) = A(j, j1) - sum * taurA * vrA[1];
304 A(j, j0) = A(j, j0) - sum * taurA * vrA[2];
305 }
306 for (idx_t j = 0; j < j3; ++j) {
307 T sum = B(j, j2) + vrA[1] * B(j, j1) + vrA[2] * B(j, j0);
308 B(j, j2) = B(j, j2) - sum * taurA;
309 B(j, j1) = B(j, j1) - sum * taurA * vrA[1];
310 B(j, j0) = B(j, j0) - sum * taurA * vrA[2];
311 }
312 if (want_z) {
313 for (idx_t j = 0; j < n; ++j) {
314 T sum = Z(j, j2) + vrA[1] * Z(j, j1) + vrA[2] * Z(j, j0);
315 Z(j, j2) = Z(j, j2) - sum * taurA;
316 Z(j, j1) = Z(j, j1) - sum * taurA * vrA[1];
317 Z(j, j0) = Z(j, j0) - sum * taurA * vrA[2];
318 }
319 }
320 }
321 else {
322 // The error is smallest when using vrB
323 for (idx_t j = 0; j < j3; ++j) {
324 T sum = A(j, j2) + vrB[1] * A(j, j1) + vrB[2] * A(j, j0);
325 A(j, j2) = A(j, j2) - sum * taurB;
326 A(j, j1) = A(j, j1) - sum * taurB * vrB[1];
327 A(j, j0) = A(j, j0) - sum * taurB * vrB[2];
328 }
329 for (idx_t j = 0; j < j3; ++j) {
330 T sum = B(j, j2) + vrB[1] * B(j, j1) + vrB[2] * B(j, j0);
331 B(j, j2) = B(j, j2) - sum * taurB;
332 B(j, j1) = B(j, j1) - sum * taurB * vrB[1];
333 B(j, j0) = B(j, j0) - sum * taurB * vrB[2];
334 }
335 if (want_z) {
336 for (idx_t j = 0; j < n; ++j) {
337 T sum = Z(j, j2) + vrB[1] * Z(j, j1) + vrB[2] * Z(j, j0);
338 Z(j, j2) = Z(j, j2) - sum * taurB;
339 Z(j, j1) = Z(j, j1) - sum * taurB * vrB[1];
340 Z(j, j0) = Z(j, j0) - sum * taurB * vrB[2];
341 }
342 }
343 }
344
345 A(j2, j0) = (T)0;
346 A(j2, j1) = (T)0;
347 B(j2, j0) = (T)0;
348 B(j2, j1) = (T)0;
349 }
350 if (n1 == 2 and n2 == 1) {
351 //
352 // Swap 2-by-2 block with 1-by-1 block
353 //
354 std::vector<T> H_;
355 auto H = new_matrix(H_, 2, 3);
356 T taur, taulA, taulB;
357 std::vector<T> vr(3);
358 std::vector<T> vlA(3);
359 std::vector<T> vlB(3);
360
361 std::vector<T> AA_;
362 auto AA = new_matrix(AA_, 3, 3);
363 std::vector<T> BB_;
364 auto BB = new_matrix(BB_, 3, 3);
365 lacpy(GENERAL, slice(A, range(j0, j3), range(j0, j3)), AA);
366 lacpy(GENERAL, slice(B, range(j0, j3), range(j0, j3)), BB);
367
368 auto a22 = A(j2, j2);
369 auto b22 = B(j2, j2);
370
371 T norma = lange(FROB_NORM, AA);
372 T normb = lange(FROB_NORM, BB);
373
374 H(0, 0) = b22 * A(j0, j0) - a22 * B(j0, j0);
375 H(0, 1) = b22 * A(j0, j1) - a22 * B(j0, j1);
376 H(0, 2) = b22 * A(j0, j2) - a22 * B(j0, j2);
377 H(1, 0) = b22 * A(j1, j0) - a22 * B(j1, j0);
378 H(1, 1) = b22 * A(j1, j1) - a22 * B(j1, j1);
379 H(1, 2) = b22 * A(j1, j2) - a22 * B(j1, j2);
380
381 inv_house3(H, vr, taur);
382
383 // Apply update from the right to the local matrices
384 for (idx_t j = 0; j < 3; ++j) {
385 T sum = AA(j, 0) + vr[1] * AA(j, 1) + vr[2] * AA(j, 2);
386 AA(j, 0) = AA(j, 0) - sum * taur;
387 AA(j, 1) = AA(j, 1) - sum * taur * vr[1];
388 AA(j, 2) = AA(j, 2) - sum * taur * vr[2];
389 }
390 for (idx_t j = 0; j < 3; ++j) {
391 T sum = BB(j, 0) + vr[1] * BB(j, 1) + vr[2] * BB(j, 2);
392 BB(j, 0) = BB(j, 0) - sum * taur;
393 BB(j, 1) = BB(j, 1) - sum * taur * vr[1];
394 BB(j, 2) = BB(j, 2) - sum * taur * vr[2];
395 }
396
397 // Determine two sets of left transformations
398 // one that zeroes out the first column of A and one that zeroes out the
399 // first column of B
400 // We will later choose the one that gives the smaller error
401 vlA[0] = AA(0, 0);
402 vlA[1] = AA(1, 0);
403 vlA[2] = AA(2, 0);
404 larfg(FORWARD, COLUMNWISE_STORAGE, vlA, taulA);
405
406 vlB[0] = BB(0, 0);
407 vlB[1] = BB(1, 0);
408 vlB[2] = BB(2, 0);
409 larfg(FORWARD, COLUMNWISE_STORAGE, vlB, taulB);
410
411 // Apply update from the left
412 for (idx_t j = 0; j < 3; ++j) {
413 T sum = AA(0, j) + vlB[1] * AA(1, j) + vlB[2] * AA(2, j);
414 AA(0, j) = AA(0, j) - sum * taulB;
415 AA(1, j) = AA(1, j) - sum * taulB * vlB[1];
416 AA(2, j) = AA(2, j) - sum * taulB * vlB[2];
417 }
418 for (idx_t j = 0; j < 3; ++j) {
419 T sum = BB(0, j) + vlA[1] * BB(1, j) + vlA[2] * BB(2, j);
420 BB(0, j) = BB(0, j) - sum * taulA;
421 BB(1, j) = BB(1, j) - sum * taulA * vlA[1];
422 BB(2, j) = BB(2, j) - sum * taulA * vlA[2];
423 }
424
425 //
426 // Determine if the swap was successful
427 //
428 T errA = lapy2(AA(1, 0), AA(2, 0));
429 T errB = lapy2(BB(1, 0), BB(2, 0));
430 const T eps = ulp<T>();
431 const T small_num = safe_min<T>();
432
433 if (errA > max((T)20 * norma * eps, small_num) and
434 errB > max((T)20 * normb * eps, small_num)) {
435 // The swap failed, return with error
436 // Note, though we don't have a proof that this will always be the
437 // case, there are currently no known cases where this swap can
438 // fail.
439 return 1;
440 }
441
442 //
443 // Swap is accepted, apply the updates to the original matrices
444 //
445 for (idx_t j = 0; j < j3; ++j) {
446 T sum = A(j, j0) + vr[1] * A(j, j1) + vr[2] * A(j, j2);
447 A(j, j0) = A(j, j0) - sum * taur;
448 A(j, j1) = A(j, j1) - sum * taur * vr[1];
449 A(j, j2) = A(j, j2) - sum * taur * vr[2];
450 }
451 for (idx_t j = 0; j < j3; ++j) {
452 T sum = B(j, j0) + vr[1] * B(j, j1) + vr[2] * B(j, j2);
453 B(j, j0) = B(j, j0) - sum * taur;
454 B(j, j1) = B(j, j1) - sum * taur * vr[1];
455 B(j, j2) = B(j, j2) - sum * taur * vr[2];
456 }
457 if (want_z) {
458 for (idx_t j = 0; j < n; ++j) {
459 T sum = Z(j, j0) + vr[1] * Z(j, j1) + vr[2] * Z(j, j2);
460 Z(j, j0) = Z(j, j0) - sum * taur;
461 Z(j, j1) = Z(j, j1) - sum * taur * vr[1];
462 Z(j, j2) = Z(j, j2) - sum * taur * vr[2];
463 }
464 }
465 if (errB * norma < errA * normb) {
466 // The error is smallest when using vlA
467 for (idx_t j = j0; j < n; ++j) {
468 T sum = A(j0, j) + vlA[1] * A(j1, j) + vlA[2] * A(j2, j);
469 A(j0, j) = A(j0, j) - sum * taulA;
470 A(j1, j) = A(j1, j) - sum * taulA * vlA[1];
471 A(j2, j) = A(j2, j) - sum * taulA * vlA[2];
472 }
473 for (idx_t j = j0; j < n; ++j) {
474 T sum = B(j0, j) + vlA[1] * B(j1, j) + vlA[2] * B(j2, j);
475 B(j0, j) = B(j0, j) - sum * taulA;
476 B(j1, j) = B(j1, j) - sum * taulA * vlA[1];
477 B(j2, j) = B(j2, j) - sum * taulA * vlA[2];
478 }
479 if (want_q) {
480 for (idx_t j = 0; j < n; ++j) {
481 T sum = Q(j, j0) + vlA[1] * Q(j, j1) + vlA[2] * Q(j, j2);
482 Q(j, j0) = Q(j, j0) - sum * taulA;
483 Q(j, j1) = Q(j, j1) - sum * taulA * vlA[1];
484 Q(j, j2) = Q(j, j2) - sum * taulA * vlA[2];
485 }
486 }
487 }
488 else {
489 // The error is smallest when using vlB
490 for (idx_t j = j0; j < n; ++j) {
491 T sum = A(j0, j) + vlB[1] * A(j1, j) + vlB[2] * A(j2, j);
492 A(j0, j) = A(j0, j) - sum * taulB;
493 A(j1, j) = A(j1, j) - sum * taulB * vlB[1];
494 A(j2, j) = A(j2, j) - sum * taulB * vlB[2];
495 }
496 for (idx_t j = j0; j < n; ++j) {
497 T sum = B(j0, j) + vlB[1] * B(j1, j) + vlB[2] * B(j2, j);
498 B(j0, j) = B(j0, j) - sum * taulB;
499 B(j1, j) = B(j1, j) - sum * taulB * vlB[1];
500 B(j2, j) = B(j2, j) - sum * taulB * vlB[2];
501 }
502 if (want_q) {
503 for (idx_t j = 0; j < n; ++j) {
504 T sum = Q(j, j0) + vlB[1] * Q(j, j1) + vlB[2] * Q(j, j2);
505 Q(j, j0) = Q(j, j0) - sum * taulB;
506 Q(j, j1) = Q(j, j1) - sum * taulB * vlB[1];
507 Q(j, j2) = Q(j, j2) - sum * taulB * vlB[2];
508 }
509 }
510 }
511
512 A(j1, j0) = (T)0;
513 A(j2, j0) = (T)0;
514 B(j1, j0) = (T)0;
515 B(j2, j0) = (T)0;
516 }
517 if (n1 == 2 and n2 == 2) {
518 //
519 // Swap 2-by-2 block with 2-by-2 block
520 //
521 std::vector<T> M_;
522 auto M = new_matrix(M_, 8, 8);
523 std::vector<T> x(8);
524 std::vector<idx_t> piv(8);
525
526 for (idx_t j = 0; j < 8; ++j)
527 for (idx_t i = 0; i < 8; ++i)
528 M(i, j) = (T)0;
529
530 // Construct matrix with kronecker structure
531 // I (x) A00
532 M(0, 0) = A(j0, j0);
533 M(0, 1) = A(j0, j1);
534 M(1, 0) = A(j1, j0);
535 M(1, 1) = A(j1, j1);
536 M(2, 2) = A(j0, j0);
537 M(2, 3) = A(j0, j1);
538 M(3, 2) = A(j1, j0);
539 M(3, 3) = A(j1, j1);
540 // I (x) B00
541 M(4, 0) = B(j0, j0);
542 M(4, 1) = B(j0, j1);
543 M(5, 0) = B(j1, j0);
544 M(5, 1) = B(j1, j1);
545 M(6, 2) = B(j0, j0);
546 M(6, 3) = B(j0, j1);
547 M(7, 2) = B(j1, j0);
548 M(7, 3) = B(j1, j1);
549 // A11T (x) I
550 M(0, 4) = -A(j2, j2);
551 M(0, 5) = -A(j3, j2);
552 M(1, 6) = -A(j2, j2);
553 M(1, 7) = -A(j3, j2);
554 M(2, 4) = -A(j2, j3);
555 M(2, 5) = -A(j3, j3);
556 M(3, 6) = -A(j2, j3);
557 M(3, 7) = -A(j3, j3);
558 // B11T (x) I
559 M(4, 4) = -B(j2, j2);
560 M(4, 5) = -B(j3, j2);
561 M(5, 6) = -B(j2, j2);
562 M(5, 7) = -B(j3, j2);
563 M(6, 4) = -B(j2, j3);
564 M(6, 5) = -B(j3, j3);
565 M(7, 6) = -B(j2, j3);
566 M(7, 7) = -B(j3, j3);
567 // RHS
568 x[0] = A(j0, j2);
569 x[1] = A(j1, j2);
570 x[2] = A(j0, j3);
571 x[3] = A(j1, j3);
572 x[4] = B(j0, j2);
573 x[5] = B(j1, j2);
574 x[6] = B(j0, j3);
575 x[7] = B(j1, j3);
576 // LU of M
577 int ierr = getrf(M, piv);
578 if (ierr != 0) {
579 return 1;
580 }
581 // Apply pivot to rhs
582 for (idx_t i = 0; i < 8; ++i) {
583 if (i != piv[i]) std::swap(x[i], x[piv[i]]);
584 }
585 // Solve Ly = rhs
586 trsv(LOWER_TRIANGLE, NO_TRANS, UNIT_DIAG, M, x);
587 // Solve Ux = y
588 trsv(UPPER_TRIANGLE, NO_TRANS, NON_UNIT_DIAG, M, x);
589
590 // Find Zc so that
591 // [ -x[0] -x[2] ] [ * * ]
592 // Zc^T [ -x[1] -x[3] ] = [ * * ]
593 // [ 1 0 ] [ 0 0 ]
594 // [ 0 1 ] [ 0 0 ]
595
596 // Rotation to make X upper triangular
597 T cxl1, sxl1;
598 rotg(x[0], x[1], cxl1, sxl1);
599 x[1] = (T)0;
600 T rottemp = cxl1 * x[2] + sxl1 * x[3];
601 x[3] = -sxl1 * x[2] + cxl1 * x[3];
602 x[2] = rottemp;
603 // SVD of (upper triangular) X
604 T cxl2, sxl2, cxr, sxr, ssx1, ssx2;
605 svd22(x[0], x[2], x[3], ssx2, ssx1, cxl2, sxl2, cxr, sxr);
606 // Fuse left rotations
607 T cxl, sxl;
608 cxl = cxl1 * cxl2 - sxl1 * sxl2;
609 sxl = cxl2 * sxl1 + sxl2 * cxl1;
610 // Rotations based on the singular values
611 ssx1 = -ssx1;
612 ssx2 = -ssx2;
613 T temp = (T)1;
614 T cx1, sx1, cx2, sx2;
615 rotg(ssx1, temp, cx1, sx1);
616 temp = (T)1;
617 rotg(ssx2, temp, cx2, sx2);
618
619 // Find Qc so that
620 // [ 1 0 ] [ 0 0 ]
621 // Qc^T [ 0 1 ] = [ 0 0 ]
622 // [ x[4] x[6] ] [ * * ]
623 // [ x[5] x[7] ] [ * * ]
624
625 // Rotation to make Y^T upper triangular
626 T cyl1, syl1;
627 rotg(x[4], x[5], cyl1, syl1);
628 x[5] = (T)0;
629 rottemp = cyl1 * x[6] + syl1 * x[7];
630 x[7] = -syl1 * x[6] + cyl1 * x[7];
631 x[6] = rottemp;
632 // SVD of (upper triangular) Y
633 T cyl2, syl2, cyr, syr, ssy1, ssy2;
634 svd22(x[4], x[6], x[7], ssy2, ssy1, cyl2, syl2, cyr, syr);
635 // Fuse left rotations
636 T cyl, syl;
637 cyl = cyl1 * cyl2 - syl1 * syl2;
638 syl = cyl2 * syl1 + syl2 * cyl1;
639 // Rotations based on the singular values
640 temp = (T)1;
641 T cy1, sy1, cy2, sy2;
642 rotg(ssy1, temp, cy1, sy1);
643 temp = (T)1;
644 rotg(ssy2, temp, cy2, sy2);
645
646 // Perform the swap on a local matrix and check the error
647 std::vector<T> AA_;
648 auto AA = new_matrix(AA_, 4, 4);
649 std::vector<T> BB_;
650 auto BB = new_matrix(BB_, 4, 4);
651 std::vector<T> QQ_;
652 auto QQ = new_matrix(QQ_, 4, 4);
653 std::vector<T> ZZ_;
654 auto ZZ = new_matrix(ZZ_, 4, 4);
655
656 lacpy(GENERAL, slice(A, range(j0, j3 + 1), range(j0, j3 + 1)), AA);
657 lacpy(GENERAL, slice(B, range(j0, j3 + 1), range(j0, j3 + 1)), BB);
658 laset(GENERAL, T(0), T(1), QQ);
659 laset(GENERAL, T(0), T(1), ZZ);
660
661 auto norma = lange(FROB_NORM, AA);
662 auto normb = lange(FROB_NORM, BB);
663
664 // Apply rotations from the left to local matrices
665 {
666 auto a0 = row(AA, 0);
667 auto a1 = row(AA, 1);
668 auto a2 = row(AA, 2);
669 auto a3 = row(AA, 3);
670 rot(a0, a1, cyr, syr);
671 rot(a2, a3, cyl, syl);
672 rot(a2, a0, cy1, sy1);
673 rot(a3, a1, cy2, sy2);
674
675 auto b0 = row(BB, 0);
676 auto b1 = row(BB, 1);
677 auto b2 = row(BB, 2);
678 auto b3 = row(BB, 3);
679 rot(b0, b1, cyr, syr);
680 rot(b2, b3, cyl, syl);
681 rot(b2, b0, cy1, sy1);
682 rot(b3, b1, cy2, sy2);
683
684 auto q0 = col(QQ, 0);
685 auto q1 = col(QQ, 1);
686 auto q2 = col(QQ, 2);
687 auto q3 = col(QQ, 3);
688
689 rot(q0, q1, cyr, syr);
690 rot(q2, q3, cyl, syl);
691 rot(q2, q0, cy1, sy1);
692 rot(q3, q1, cy2, sy2);
693 }
694 // Apply rotations from the right to local matrices
695 {
696 auto a0 = col(AA, 0);
697 auto a1 = col(AA, 1);
698 auto a2 = col(AA, 2);
699 auto a3 = col(AA, 3);
700 rot(a0, a1, cxl, sxl);
701 rot(a2, a3, cxr, sxr);
702 rot(a0, a2, cx1, sx1);
703 rot(a1, a3, cx2, sx2);
704
705 auto b0 = col(BB, 0);
706 auto b1 = col(BB, 1);
707 auto b2 = col(BB, 2);
708 auto b3 = col(BB, 3);
709 rot(b0, b1, cxl, sxl);
710 rot(b2, b3, cxr, sxr);
711 rot(b0, b2, cx1, sx1);
712 rot(b1, b3, cx2, sx2);
713
714 auto z0 = col(ZZ, 0);
715 auto z1 = col(ZZ, 1);
716 auto z2 = col(ZZ, 2);
717 auto z3 = col(ZZ, 3);
718 rot(z0, z1, cxl, sxl);
719 rot(z2, z3, cxr, sxr);
720 rot(z0, z2, cx1, sx1);
721 rot(z1, z3, cx2, sx2);
722 }
723
724 // Weak stability test
725 auto enorma = lange(FROB_NORM, slice(AA, range(2, 4), range(0, 2)));
726 auto enormb = lange(FROB_NORM, slice(BB, range(2, 4), range(0, 2)));
727 const T eps = ulp<T>();
728 const T small_num = safe_min<T>();
729
730 idx_t iter = 0;
731 T tolA = max((T)20 * norma * eps, small_num);
732 T tolB = max((T)20 * normb * eps, small_num);
733 const idx_t max_iter = 6;
734 while (iter < max_iter) {
735 if (enorma <= tolA and enormb <= tolB) break;
736 if (iter == max_iter - 1) {
737 return 1;
738 }
739 // The swap is not (yet) accepted, apply iterative refinement to
740 // try to improve the swap.
741
742 for (idx_t j = 0; j < 8; ++j)
743 for (idx_t i = 0; i < 8; ++i)
744 M(i, j) = (T)0;
745
747 // I (x) AA(2:3,2:3)
748 M(0, 0) = AA(2, 2);
749 M(0, 1) = AA(2, 3);
750 M(1, 0) = AA(3, 2);
751 M(1, 1) = AA(3, 3);
752 M(2, 2) = AA(2, 2);
753 M(2, 3) = AA(2, 3);
754 M(3, 2) = AA(3, 2);
755 M(3, 3) = AA(3, 3);
756 // -AA(0:1,0:1)^T (x) I
757 M(0, 4) = -AA(0, 0);
758 M(0, 5) = -AA(1, 0);
759 M(1, 6) = -AA(0, 0);
760 M(1, 7) = -AA(1, 0);
761 M(2, 4) = -AA(0, 1);
762 M(2, 5) = -AA(1, 1);
763 M(3, 6) = -AA(0, 1);
764 M(3, 7) = -AA(1, 1);
765
766 // I (x) BB(2:3,2:3)
767 M(4, 0) = BB(2, 2);
768 M(4, 1) = BB(2, 3);
769 M(5, 0) = BB(3, 2);
770 M(5, 1) = BB(3, 3);
771 M(6, 2) = BB(2, 2);
772 M(6, 3) = BB(2, 3);
773 M(7, 2) = BB(3, 2);
774 M(7, 3) = BB(3, 3);
775 // -BB(0:1,0:1)^T (x) I
776 M(4, 4) = -BB(0, 0);
777 M(4, 5) = -BB(1, 0);
778 M(5, 6) = -BB(0, 0);
779 M(5, 7) = -BB(1, 0);
780 M(6, 4) = -BB(0, 1);
781 M(6, 5) = -BB(1, 1);
782 M(7, 6) = -BB(0, 1);
783 M(7, 7) = -BB(1, 1);
784
785 // RHS
786 x[0] = AA(2, 0);
787 x[1] = AA(3, 0);
788 x[2] = AA(2, 1);
789 x[3] = AA(3, 1);
790 x[4] = BB(2, 0);
791 x[5] = BB(3, 0);
792 x[6] = BB(2, 1);
793 x[7] = BB(3, 1);
794
795 // LU of M
796 int ierr = getrf(M, piv);
797 if (ierr != 0) return 1;
798 // Apply pivot to rhs
799 for (idx_t i = 0; i < 8; ++i) {
800 if (i != piv[i]) std::swap(x[i], x[piv[i]]);
801 }
802 // Solve Ly = rhs
803 trsv(LOWER_TRIANGLE, NO_TRANS, UNIT_DIAG, M, x);
804 // Solve Ux = y
805 trsv(UPPER_TRIANGLE, NO_TRANS, NON_UNIT_DIAG, M, x);
806
807 // Find Zc so that
808 // [ 1 0 ] [ * * ]
809 // Zc^T [ 0 1 ] = [ * * ]
810 // [ -x[0] -x[2] ] [ 0 0 ]
811 // [ -x[1] -x[3] ] [ 0 0 ]
812
813 // Rotation to make X upper triangular
814 T cxl1, sxl1;
815 rotg(x[0], x[1], cxl1, sxl1);
816 x[1] = (T)0;
817 T rottemp = cxl1 * x[2] + sxl1 * x[3];
818 x[3] = -sxl1 * x[2] + cxl1 * x[3];
819 x[2] = rottemp;
820 // SVD of (upper triangular) X
821 T cxl2, sxl2, cxr, sxr, ssx1, ssx2;
822 svd22(x[0], x[2], x[3], ssx2, ssx1, cxl2, sxl2, cxr, sxr);
823 // Fuse left rotations
824 T cxl, sxl;
825 cxl = cxl1 * cxl2 - sxl1 * sxl2;
826 sxl = cxl2 * sxl1 + sxl2 * cxl1;
827 // Rotations based on the singular values
828 ssx1 = -ssx1;
829 ssx2 = -ssx2;
830 T temp = (T)1;
831 T cx1, sx1, cx2, sx2;
832 rotg(temp, ssx1, cx1, sx1);
833 temp = (T)1;
834 rotg(temp, ssx2, cx2, sx2);
835
836 // Find Qc so that
837 // [ 1 0 ] [ 0 0 ]
838 // Qc^T [ 0 1 ] = [ 0 0 ]
839 // [ x[4] x[6] ] [ * * ]
840 // [ x[5] x[7] ] [ * * ]
841
842 std::swap(x[5], x[6]);
843 T cyl1, syl1;
844 rotg(x[4], x[5], cyl1, syl1);
845 x[5] = (T)0;
846 rottemp = cyl1 * x[6] + syl1 * x[7];
847 x[7] = -syl1 * x[6] + cyl1 * x[7];
848 x[6] = rottemp;
849 // SVD of (upper triangular) X
850 T cyl2, syl2, cyr, syr, ssy1, ssy2;
851 svd22(x[4], x[6], x[7], ssy2, ssy1, cyl2, syl2, cyr, syr);
852 // Fuse left rotations
853 T cyl, syl;
854 cyl = cyl1 * cyl2 - syl1 * syl2;
855 syl = cyl2 * syl1 + syl2 * cyl1;
856 // Rotations based on the singular values
857 ssy1 = -ssy1;
858 ssy2 = -ssy2;
859 temp = (T)1;
860 T cy1, sy1, cy2, sy2;
861 rotg(temp, ssy1, cy1, sy1);
862 temp = (T)1;
863 rotg(temp, ssy2, cy2, sy2);
864
865 // Apply rotations from the left to local matrices
866 {
867 auto a0 = row(AA, 0);
868 auto a1 = row(AA, 1);
869 auto a2 = row(AA, 2);
870 auto a3 = row(AA, 3);
871 rot(a0, a1, cyr, syr);
872 rot(a2, a3, cyl, syl);
873 rot(a0, a2, cy1, sy1);
874 rot(a1, a3, cy2, sy2);
875
876 auto b0 = row(BB, 0);
877 auto b1 = row(BB, 1);
878 auto b2 = row(BB, 2);
879 auto b3 = row(BB, 3);
880 rot(b0, b1, cyr, syr);
881 rot(b2, b3, cyl, syl);
882 rot(b0, b2, cy1, sy1);
883 rot(b1, b3, cy2, sy2);
884
885 auto q0 = col(QQ, 0);
886 auto q1 = col(QQ, 1);
887 auto q2 = col(QQ, 2);
888 auto q3 = col(QQ, 3);
889
890 rot(q0, q1, cyr, syr);
891 rot(q2, q3, cyl, syl);
892 rot(q0, q2, cy1, sy1);
893 rot(q1, q3, cy2, sy2);
894 }
895 // Apply rotations from the right to local matrices
896 {
897 auto a0 = col(AA, 0);
898 auto a1 = col(AA, 1);
899 auto a2 = col(AA, 2);
900 auto a3 = col(AA, 3);
901 rot(a0, a1, cxr, sxr);
902 rot(a2, a3, cxl, sxl);
903 rot(a0, a2, cx1, sx1);
904 rot(a1, a3, cx2, sx2);
905
906 auto b0 = col(BB, 0);
907 auto b1 = col(BB, 1);
908 auto b2 = col(BB, 2);
909 auto b3 = col(BB, 3);
910 rot(b0, b1, cxr, sxr);
911 rot(b2, b3, cxl, sxl);
912 rot(b0, b2, cx1, sx1);
913 rot(b1, b3, cx2, sx2);
914
915 auto z0 = col(ZZ, 0);
916 auto z1 = col(ZZ, 1);
917 auto z2 = col(ZZ, 2);
918 auto z3 = col(ZZ, 3);
919 rot(z0, z1, cxr, sxr);
920 rot(z2, z3, cxl, sxl);
921 rot(z0, z2, cx1, sx1);
922 rot(z1, z3, cx2, sx2);
923 }
924
925 enorma = lange(FROB_NORM, slice(AA, range(2, 4), range(0, 2)));
926 enormb = lange(FROB_NORM, slice(BB, range(2, 4), range(0, 2)));
927 iter++;
928 }
929
930 // TODO: this is a large workspace, add it to the interface
931 std::vector<T> workl_;
932 auto workl = new_matrix(workl_, 4, n);
933 std::vector<T> workr_;
934 auto workr = new_matrix(workr_, n, 4);
935
936 // Apply QQ
937 {
938 auto A_slice = slice(A, range(j0, j3 + 1), range(j0, n));
939 auto work_slice = slice(workl, range(0, 4), range(j0, n));
940 gemm(TRANSPOSE, NO_TRANS, (T)1, QQ, A_slice, StrongZero(T(0)),
941 work_slice);
942 lacpy(GENERAL, work_slice, A_slice);
943 }
944 {
945 auto B_slice = slice(B, range(j0, j3 + 1), range(j0, n));
946 auto work_slice = slice(workl, range(0, 4), range(j0, n));
947 gemm(TRANSPOSE, NO_TRANS, (T)1, QQ, B_slice, StrongZero(T(0)),
948 work_slice);
949 lacpy(GENERAL, work_slice, B_slice);
950 }
951 {
952 auto Q_slice = slice(Q, range(0, n), range(j0, j3 + 1));
953 auto work_slice = slice(workr, range(0, n), range(0, 4));
954 gemm(NO_TRANS, NO_TRANS, (T)1, Q_slice, QQ, StrongZero(T(0)),
955 work_slice);
956 lacpy(GENERAL, work_slice, Q_slice);
957 }
958 // Apply ZZ
959 {
960 auto A_slice =
961 slice(A, range(0, min(n, j3 + 1)), range(j0, j3 + 1));
962 auto work_slice =
963 slice(workr, range(0, min(n, j3 + 1)), range(0, 4));
964 gemm(NO_TRANS, NO_TRANS, (T)1, A_slice, ZZ, StrongZero(T(0)),
965 work_slice);
966 lacpy(GENERAL, work_slice, A_slice);
967 }
968 {
969 auto B_slice =
970 slice(B, range(0, min(n, j3 + 1)), range(j0, j3 + 1));
971 auto work_slice =
972 slice(workr, range(0, min(n, j3 + 1)), range(0, 4));
973 gemm(NO_TRANS, NO_TRANS, (T)1, B_slice, ZZ, StrongZero(T(0)),
974 work_slice);
975 lacpy(GENERAL, work_slice, B_slice);
976 }
977 {
978 auto Z_slice = slice(Z, range(0, n), range(j0, j3 + 1));
979 auto work_slice = slice(workr, range(0, n), range(0, 4));
980 gemm(NO_TRANS, NO_TRANS, (T)1, Z_slice, ZZ, StrongZero(T(0)),
981 work_slice);
982 lacpy(GENERAL, work_slice, Z_slice);
983 }
984
985 // Set relevant parts to zero
986 A(j2, j0) = (T)0;
987 A(j3, j0) = (T)0;
988 A(j2, j1) = (T)0;
989 A(j3, j1) = (T)0;
990 B(j2, j0) = (T)0;
991 B(j3, j0) = (T)0;
992 B(j2, j1) = (T)0;
993 B(j3, j1) = (T)0;
994 }
995
996 // Standardize the 2x2 Schur blocks (if any)
997 if (n2 == 2) {
998 auto A22 = slice(A, range(j0, j2), range(j0, j2));
999 auto B22 = slice(B, range(j0, j2), range(j0, j2));
1000
1001 complex_type<T> alpha1, alpha2;
1002 real_type<T> beta1, beta2;
1003
1004 // Make B upper triangular
1005 T cl1, sl1;
1006 rotg(B(j0, j0), B(j1, j0), cl1, sl1);
1007 B(j1, j0) = (T)0;
1008 {
1009 auto b1 = slice(B, j0, range(j1, j2));
1010 auto b2 = slice(B, j1, range(j1, j2));
1011 rot(b1, b2, cl1, sl1);
1012 auto a1 = slice(A, j0, range(j0, j2));
1013 auto a2 = slice(A, j1, range(j0, j2));
1014 rot(a1, a2, cl1, sl1);
1015 }
1016
1017 T cl2, sl2, cr, sr, scal0, scal1;
1018 lahqz_schur22(A22, B22, alpha1, alpha2, beta1, beta2, cl2, sl2, cr, sr,
1019 scal0, scal1);
1020 // Fuse left rotations
1021 T cl, sl;
1022 cl = cl1 * cl2 - sl1 * sl2;
1023 sl = cl2 * sl1 + sl2 * cl1;
1024 // Apply left rotation
1025 {
1026 auto a1 = slice(A, j0, range(j2, n));
1027 auto a2 = slice(A, j1, range(j2, n));
1028 rot(a1, a2, cl, sl);
1029 auto b1 = slice(B, j0, range(j2, n));
1030 auto b2 = slice(B, j1, range(j2, n));
1031 rot(b1, b2, cl, sl);
1032 auto q0 = col(Q, j0);
1033 auto q1 = col(Q, j1);
1034 rot(q0, q1, cl, sl);
1035 }
1036 // Apply right rotation and scaling
1037 {
1038 auto a1 = slice(A, range(0, j0), j0);
1039 auto a2 = slice(A, range(0, j0), j1);
1040 rot(a1, a2, cr, sr);
1041 auto b1 = slice(B, range(0, j0), j0);
1042 auto b2 = slice(B, range(0, j0), j1);
1043 rot(b1, b2, cr, sr);
1044 auto z0 = col(Z, j0);
1045 auto z1 = col(Z, j1);
1046 rot(z0, z1, cr, sr);
1047
1048 if (scal0 != (T)1) {
1049 scal(scal0, a1);
1050 scal(scal0, b1);
1051 scal(scal0, z0);
1052 }
1053
1054 if (scal1 != (T)1) {
1055 scal(scal1, a2);
1056 scal(scal1, b2);
1057 scal(scal1, z1);
1058 }
1059 }
1060 }
1061 if (n1 == 2) {
1062 // Make B upper triangular
1063 T cl1, sl1;
1064 rotg(B(j0 + n2, j0 + n2), B(j1 + n2, j0 + n2), cl1, sl1);
1065 B(j1 + n2, j0 + n2) = (T)0;
1066 {
1067 auto b1 = slice(B, j0 + n2, range(j1 + n2, j2 + n2));
1068 auto b2 = slice(B, j1 + n2, range(j1 + n2, j2 + n2));
1069 rot(b1, b2, cl1, sl1);
1070 auto a1 = slice(A, j0 + n2, range(j0 + n2, j2 + n2));
1071 auto a2 = slice(A, j1 + n2, range(j0 + n2, j2 + n2));
1072 rot(a1, a2, cl1, sl1);
1073 }
1074
1075 complex_type<T> alpha1, alpha2;
1076 real_type<T> beta1, beta2;
1077 auto A22 = slice(A, range(j0 + n2, j2 + n2), range(j0 + n2, j2 + n2));
1078 auto B22 = slice(B, range(j0 + n2, j2 + n2), range(j0 + n2, j2 + n2));
1079 T cl2, sl2, cr, sr, scal0, scal1;
1080 lahqz_schur22(A22, B22, alpha1, alpha2, beta1, beta2, cl2, sl2, cr, sr,
1081 scal0, scal1);
1082 // Fuse left rotations
1083 T cl, sl;
1084 cl = cl1 * cl2 - sl1 * sl2;
1085 sl = cl2 * sl1 + sl2 * cl1;
1086 // Apply left rotation
1087 {
1088 auto a1 = slice(A, j0 + n2, range(j2 + n2, n));
1089 auto a2 = slice(A, j1 + n2, range(j2 + n2, n));
1090 rot(a1, a2, cl, sl);
1091 auto b1 = slice(B, j0 + n2, range(j2 + n2, n));
1092 auto b2 = slice(B, j1 + n2, range(j2 + n2, n));
1093 rot(b1, b2, cl, sl);
1094 auto q0 = col(Q, j0 + n2);
1095 auto q1 = col(Q, j1 + n2);
1096 rot(q0, q1, cl, sl);
1097 }
1098 // Apply right rotation
1099 {
1100 auto a1 = slice(A, range(0, j0 + n2), j0 + n2);
1101 auto a2 = slice(A, range(0, j0 + n2), j1 + n2);
1102 rot(a1, a2, cr, sr);
1103 auto b1 = slice(B, range(0, j0 + n2), j0 + n2);
1104 auto b2 = slice(B, range(0, j0 + n2), j1 + n2);
1105 rot(b1, b2, cr, sr);
1106 auto z0 = col(Z, j0 + n2);
1107 auto z1 = col(Z, j1 + n2);
1108 rot(z0, z1, cr, sr);
1109
1110 if (scal0 != (T)1) {
1111 scal(scal0, a1);
1112 scal(scal0, b1);
1113 scal(scal0, z0);
1114 }
1115
1116 if (scal1 != (T)1) {
1117 scal(scal1, a2);
1118 scal(scal1, b2);
1119 scal(scal1, z1);
1120 }
1121 }
1122 }
1123
1124 return 0;
1125}
1126
1133template <TLAPACK_CSMATRIX matrix_t,
1134 enable_if_t<is_complex<type_t<matrix_t>>, bool> = true>
1135int generalized_schur_swap(bool want_q,
1136 bool want_z,
1137 matrix_t& A,
1138 matrix_t& B,
1139 matrix_t& Q,
1140 matrix_t& Z,
1141 const size_type<matrix_t>& j0,
1142 const size_type<matrix_t>& n1,
1143 const size_type<matrix_t>& n2)
1144{
1145 using idx_t = size_type<matrix_t>;
1146 using T = type_t<matrix_t>;
1147 using real_t = real_type<T>;
1148 using range = pair<idx_t, idx_t>;
1149
1150 const idx_t n = ncols(A);
1151
1152 tlapack_check(nrows(A) == n);
1153 tlapack_check(nrows(Q) == n);
1154 tlapack_check(ncols(Q) == n);
1155 tlapack_check(0 <= j0 and j0 < n);
1156 tlapack_check(n1 == 1);
1157 tlapack_check(n2 == 1);
1158
1159 const idx_t j1 = j0 + 1;
1160
1161 //
1162 // In the complex case, there can only be 1x1 blocks to swap
1163 //
1164 const T a00 = A(j0, j0);
1165 const T a01 = A(j0, j1);
1166 const T a11 = A(j1, j1);
1167 const T b00 = B(j0, j0);
1168 const T b01 = B(j0, j1);
1169 const T b11 = B(j1, j1);
1170
1171 const bool use_b = abs(b11 * a00) > abs(b00 * a11);
1172
1173 //
1174 // Determine the transformation to perform the interchange
1175 //
1176 real_t cl, cr;
1177 T sl, sr;
1178 T temp = b11 * a00 - a11 * b00;
1179 T temp2 = b11 * a01 - a11 * b01;
1180 rotg(temp2, temp, cr, sr);
1181
1182 // Apply transformation from the right
1183 {
1184 auto a1 = slice(A, range{0, j1 + 1}, j0);
1185 auto a2 = slice(A, range{0, j1 + 1}, j1);
1186 rot(a2, a1, cr, sr);
1187 auto b1 = slice(B, range{0, j1 + 1}, j0);
1188 auto b2 = slice(B, range{0, j1 + 1}, j1);
1189 rot(b2, b1, cr, sr);
1190 if (want_z) {
1191 auto z1 = col(Z, j0);
1192 auto z2 = col(Z, j1);
1193 rot(z2, z1, cr, sr);
1194 }
1195 }
1196
1197 if (use_b) {
1198 temp = B(j0, j0);
1199 temp2 = B(j1, j0);
1200 }
1201 else {
1202 temp = A(j0, j0);
1203 temp2 = A(j1, j0);
1204 }
1205 rotg(temp, temp2, cl, sl);
1206
1207 // Apply transformation from the left
1208 {
1209 auto a1 = slice(A, j0, range{j0, n});
1210 auto a2 = slice(A, j1, range{j0, n});
1211 rot(a1, a2, cl, sl);
1212 auto b1 = slice(B, j0, range{j0, n});
1213 auto b2 = slice(B, j1, range{j0, n});
1214 rot(b1, b2, cl, sl);
1215 if (want_q) {
1216 auto q1 = col(Q, j0);
1217 auto q2 = col(Q, j1);
1218 rot(q1, q2, cl, conj(sl));
1219 }
1220 }
1221
1222 A(j1, j0) = (T)0;
1223 B(j1, j0) = (T)0;
1224
1225 return 0;
1226}
1227
1228} // namespace tlapack
1229
1230#endif // TLAPACK_GENERALIZED_SCHUR_SWAP_HH
TLAPACK_CSMATRIX
#define TLAPACK_CSMATRIX
Macro for tlapack::concepts::ConstructableAndSliceableMatrix compatible with C++17.
Definition concepts.hpp:961
tlapack::lapy2
real_type< TX, TY > lapy2(const TX &x, const TY &y)
Finds , taking care not to cause unnecessary overflow.
Definition lapy2.hpp:34
tlapack::generalized_schur_swap
int generalized_schur_swap(bool want_q, bool want_z, matrix_t &A, matrix_t &B, matrix_t &Q, matrix_t &Z, const size_type< matrix_t > &j0, const size_type< matrix_t > &n1, const size_type< matrix_t > &n2)
schur_swap, swaps 2 eigenvalues of A.
Definition generalized_schur_swap.hpp:60
tlapack::lange
auto lange(norm_t normType, const matrix_t &A)
Calculates the norm of a matrix.
Definition lange.hpp:38
tlapack::svd22
void svd22(const T &f, const T &g, const T &h, T &ssmin, T &ssmax, T &csl, T &snl, T &csr, T &snr)
Computes the singular value decomposition of a 2-by-2 real triangular matrix.
Definition svd22.hpp:55
tlapack::laset
void laset(uplo_t uplo, const type_t< matrix_t > &alpha, const type_t< matrix_t > &beta, matrix_t &A)
Initializes a matrix to diagonal and off-diagonal values.
Definition laset.hpp:38
tlapack::larfg
void larfg(storage_t storeMode, type_t< vector_t > &alpha, vector_t &x, type_t< vector_t > &tau)
Generates a elementary Householder reflection.
Definition larfg.hpp:73
tlapack::lacpy
void lacpy(uplo_t uplo, const matrixA_t &A, matrixB_t &B)
Copies a matrix from A to B.
Definition lacpy.hpp:38
tlapack::inv_house3
void inv_house3(const matrix_t &A, vector_t &v, type_t< vector_t > &tau)
Inv_house calculates a reflector to reduce the first column in a 2x3 matrix A from the right to zero.
Definition inv_house3.hpp:44
tlapack::rotg
void rotg(T &a, T &b, T &c, T &s)
Construct plane rotation that eliminates b, such that:
Definition rotg.hpp:39
tlapack::rot
void rot(vectorX_t &x, vectorY_t &y, const c_type &c, const s_type &s)
Apply plane rotation:
Definition rot.hpp:44
tlapack::scal
void scal(const alpha_t &alpha, vector_t &x)
Scale vector by constant, .
Definition scal.hpp:30
tlapack::trsv
void trsv(Uplo uplo, Op trans, Diag diag, const matrixA_t &A, vectorX_t &x)
Solve the triangular matrix-vector equation.
Definition trsv.hpp:64
tlapack::syr
void syr(Uplo uplo, const alpha_t &alpha, const vectorX_t &x, matrixA_t &A)
Symmetric matrix rank-1 update:
Definition syr.hpp:45
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_check
#define tlapack_check(cond)
Throw an error if cond is false.
Definition exceptionHandling.hpp:98
tlapack::getrf
int getrf(matrix_t &A, piv_t &piv, const GetrfOpts &opts={})
getrf computes an LU factorization of a general m-by-n matrix A.
Definition getrf.hpp:64
inv_house3.hpp
lahqz_eig22.hpp
lahqz_schur22.hpp
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::FROB_NORM
constexpr internal::FrobNorm FROB_NORM
Frobenius norm of matrices.
Definition types.hpp:354
tlapack::LOWER_TRIANGLE
constexpr internal::LowerTriangle LOWER_TRIANGLE
Lower Triangle access.
Definition types.hpp:188
tlapack::UPPER_TRIANGLE
constexpr internal::UpperTriangle UPPER_TRIANGLE
Upper Triangle access.
Definition types.hpp:186
tlapack::conj
constexpr T conj(const T &x) noexcept
Extends std::conj() to real datatypes.
Definition utils.hpp:100
tlapack::FORWARD
constexpr internal::Forward FORWARD
Forward direction.
Definition types.hpp:388
tlapack::TRANSPOSE
constexpr internal::Transpose TRANSPOSE
transpose
Definition types.hpp:262
tlapack::lahqz_schur22
void lahqz_schur22(A_t &A, B_t &B, complex_type< type_t< A_t > > &alpha1, complex_type< type_t< A_t > > &alpha2, real_type< type_t< A_t > > &beta1, real_type< type_t< A_t > > &beta2, real_type< type_t< A_t > > &cl, type_t< A_t > &sl, real_type< type_t< A_t > > &cr, type_t< A_t > &sr, type_t< A_t > &scal0, type_t< A_t > &scal1)
Computes the generalized Schur factorization of a 2x2 pencil (A,B) with B upper triangular.
Definition lahqz_schur22.hpp:63
tlapack::UNIT_DIAG
constexpr internal::UnitDiagonal UNIT_DIAG
The main diagonal is assumed to consist of 1's.
Definition types.hpp:222
tlapack::GENERAL
constexpr internal::GeneralAccess GENERAL
General access.
Definition types.hpp:180
tlapack::NON_UNIT_DIAG
constexpr internal::NonUnitDiagonal NON_UNIT_DIAG
The main diagonal is not assumed to consist of 1's.
Definition types.hpp:220
tlapack::COLUMNWISE_STORAGE
constexpr internal::ColumnwiseStorage COLUMNWISE_STORAGE
Columnwise storage.
Definition types.hpp:421
tlapack::NO_TRANS
constexpr internal::NoTranspose NO_TRANS
no transpose
Definition types.hpp:260
tlapack::StrongZero
Strong zero type.
Definition StrongZero.hpp:43