1 /*
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
35 $ INFO )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LAQTR_TCC
44 #define FLENS_LAPACK_EIG_LAQTR_TCC 1
45
46 #include <flens/aux/aux.h>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
55 typename VWORK>
56 typename DenseVector<VX>::IndexType
57 laqtr_generic(bool trans,
58 bool real,
59 const GeMatrix<MT> &T,
60 const DenseVector<VB> &b,
61 const W &w,
62 SCALE &scale,
63 DenseVector<VX> &x,
64 DenseVector<VWORK> &work)
65 {
66 using std::abs;
67 using flens::max;
68
69 typedef typename DenseVector<VX>::ElementType ElementType;
70 typedef typename DenseVector<VX>::IndexType IndexType;
71
72 const Underscore<IndexType> _;
73 const IndexType n = T.numRows();
74
75 const ElementType Zero(0), One(1);
76 //
77 // .. Local Arrays ..
78 //
79 ElementType _dataD[4], _dataV[4];
80 GeMatrixView<ElementType>
81 D = typename GeMatrixView<ElementType>::Engine(2, 2, _dataD, 2),
82 V = typename GeMatrixView<ElementType>::Engine(2, 2, _dataV, 2);
83
84 IndexType info = 0;
85 //
86 // Quick return if possible
87 //
88 if (n==0) {
89 return info;
90 }
91 //
92 // Set constants to control overflow
93 //
94 const ElementType eps = lamch<ElementType>(Precision);
95 const ElementType smallNum = lamch<ElementType>(SafeMin)/eps;
96 const ElementType bigNum = One/smallNum;
97
98 ElementType xNorm = lan(MaximumNorm, T);
99 if (!real) {
100 typedef typename DenseVector<VB>::ElementType TB;
101 const GeMatrixConstView<TB> B(n, 1, b, n);
102
103 xNorm = max(xNorm, abs(w), lan(MaximumNorm, B));
104 }
105 const ElementType sMin = max(smallNum, eps*xNorm);
106 //
107 // Compute 1-norm of each column of strictly upper triangular
108 // part of T to control overflow in triangular solver.
109 //
110 work(1) = Zero;
111 for (IndexType j=2; j<=n; ++j) {
112 work(j) = blas::asum(T(_(1,j-1),j));
113 }
114
115 if (!real) {
116 for (IndexType i=2; i<=n; ++i) {
117 work(i) += abs(b(i));
118 }
119 }
120
121 const IndexType n2 = 2*n;
122 const IndexType n1 = (real) ? n : n2;
123
124 const IndexType k = blas::iamax(x(_(1,n1)));
125 ElementType xMax = abs(x(k));
126 scale = One;
127
128 if (xMax>bigNum) {
129 scale = bigNum / xMax;
130 x(_(1,n1)) *= scale;
131 xMax = bigNum;
132 }
133
134 if (real) {
135
136 if (!trans) {
137 //
138 // Solve T*p = scale*c
139 //
140 IndexType jNext = n;
141 for (IndexType j=n; j>=1; --j) {
142 if (j>jNext) {
143 continue;
144 }
145 IndexType j1 = j;
146 IndexType j2 = j;
147 jNext = j - 1;
148 if (j>1) {
149 if (T(j,j-1)!=Zero) {
150 j1 = j - 1;
151 jNext = j - 2;
152 }
153 }
154
155 if (j1==j2) {
156 //
157 // Meet 1 by 1 diagonal block
158 //
159 // Scale to avoid overflow when computing
160 // x(j) = b(j)/T(j,j)
161 //
162 ElementType xj = abs(x(j1));
163 ElementType Tjj = abs(T(j1,j1));
164 ElementType tmp = T(j1,j1);
165 if (Tjj<sMin) {
166 tmp = sMin;
167 Tjj = sMin;
168 info = 1;
169 }
170
171 if (xj==Zero) {
172 continue;
173 }
174
175 if (Tjj<One) {
176 if (xj>bigNum*Tjj) {
177 const ElementType rec = One / xj;
178 x(_(1,n)) *= rec;
179 scale *= rec;
180 xMax *= rec;
181 }
182 }
183 x(j1) /= tmp;
184 xj = abs(x(j1));
185 //
186 // Scale x if necessary to avoid overflow when adding a
187 // multiple of column j1 of T.
188 //
189 if (xj>One) {
190 const ElementType rec = One / xj;
191 if (work(j1)>(bigNum-xMax)*rec) {
192 x(_(1,n)) *= rec;
193 scale *= rec;
194 }
195 }
196 if (j1>1) {
197 auto _x = x(_(1,j1-1));
198
199 _x -= x(j1)*T(_(1,j1-1),j1);
200 const IndexType k = blas::iamax(_x);
201 xMax = abs(x(k));
202 }
203
204 } else {
205 //
206 // Meet 2 by 2 diagonal block
207 //
208 // Call 2 by 2 linear system solve, to take
209 // care of possible overflow by scaling factor.
210 //
211 D(1,1) = x(j1);
212 D(2,1) = x(j2);
213
214 ElementType scaleC;
215 IndexType iErr = laln2(false, IndexType(1), sMin, One,
216 T(_(j1,j1+1),_(j1,j1+1)), One, One,
217 D(_(1,2),_(1,1)), Zero, Zero,
218 V(_(1,2),_(1,1)), scaleC, xNorm);
219 if (iErr!=0) {
220 info = 2;
221 }
222
223 if (scaleC!=One) {
224 x(_(1,n)) *= scaleC;
225 scale *= scaleC;
226 }
227 x(j1) = V(1,1);
228 x(j2) = V(2,1);
229 //
230 // Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
231 // to avoid overflow in updating right-hand side.
232 //
233 const ElementType xj = max(abs(V(1,1)), abs(V(2,1)));
234 if (xj>One) {
235 const ElementType rec = One / xj;
236 if (max(work(j1),work(j2))>(bigNum-xMax)*rec) {
237 x(_(1,n)) *= rec;
238 scale *= rec;
239 }
240 }
241 //
242 // Update right-hand side
243 //
244 if (j1>1) {
245 auto _x = x(_(1,j1-1));
246 _x -= x(j1) * T(_(1,j1-1),j1);
247 _x -= x(j2) * T(_(1,j1-1),j2);
248 const IndexType k = blas::iamax(_x);
249 xMax = abs(x(k));
250 }
251
252 }
253
254 }
255
256 } else {
257 //
258 // Solve T**T*p = scale*c
259 //
260 IndexType jNext = 1;
261 for (IndexType j=1; j<=n; ++j) {
262 if (j<jNext) {
263 continue;
264 }
265 IndexType j1 = j;
266 IndexType j2 = j;
267 jNext = j + 1;
268 if (j<n) {
269 if (T(j+1,j)!=Zero) {
270 j2 = j + 1;
271 jNext = j + 2;
272 }
273 }
274
275 if (j1==j2) {
276 //
277 // 1 by 1 diagonal block
278 //
279 // Scale if necessary to avoid overflow in forming the
280 // right-hand side element by inner product.
281 //
282 ElementType xj = abs(x(j1));
283 if (xMax>One) {
284 const ElementType rec = One / xMax;
285 if (work(j1)>(bigNum-xj)*rec) {
286 x(_(1,n)) *= rec;
287 scale *= rec;
288 xMax *= rec;
289 }
290 }
291
292 x(j1) -= T(_(1,j1-1),j1) * x(_(1,j1-1));
293
294 xj = abs(x(j1));
295 ElementType Tjj = abs(T(j1,j1));
296 ElementType tmp = T(j1,j1);
297 if (Tjj<sMin) {
298 tmp = sMin;
299 Tjj = sMin;
300 info = 1;
301 }
302
303 if (Tjj<One) {
304 if (xj>bigNum*Tjj) {
305 const ElementType rec = One / xj;
306 x(_(1,n)) *= rec;
307 scale *= rec;
308 xMax *= rec;
309 }
310 }
311 x(j1) /= tmp;
312 xMax = max(xMax, abs(x(j1)));
313
314 } else {
315 //
316 // 2 by 2 diagonal block
317 //
318 // Scale if necessary to avoid overflow in forming the
319 // right-hand side elements by inner product.
320 //
321 const ElementType xj = max(abs(x(j1)), abs(x(j2)));
322 if (xMax>One) {
323 const ElementType rec = One / xMax;
324 if (max(work(j2),work(j1))>(bigNum-xj)*rec) {
325 x(_(1,n)) *= rec;
326 scale *= rec;
327 xMax *= rec;
328 }
329 }
330
331 D(1,1) = x(j1) - T(_(1,j1-1),j1) * x(_(1,j1-1));
332 D(2,1) = x(j2) - T(_(1,j1-1),j2) * x(_(1,j1-1));
333
334 ElementType scaleC;
335 IndexType iErr = laln2(true, IndexType(1), sMin, One,
336 T(_(j1,j1+1),_(j1,j1+1)), One, One,
337 D(_(1,2),_(1,1)), Zero, Zero,
338 V(_(1,2),_(1,1)), scaleC, xNorm);
339 if (iErr!=0) {
340 info = 2;
341 }
342
343 if (scaleC!=One) {
344 x(_(1,n)) *= scaleC;
345 scale *= scaleC;
346 }
347 x(j1) = V(1,1);
348 x(j2) = V(2,1);
349 xMax = max(abs(x(j1)), abs(x(j2)), xMax);
350
351 }
352 }
353 }
354
355 } else {
356
357 const ElementType sMinW = max(eps*abs(w), sMin);
358 if (!trans) {
359 //
360 // Solve (T + iB)*(p+iq) = c+id
361 //
362 IndexType jNext = n;
363 for (IndexType j=n; j>=1; --j) {
364 if (j>jNext) {
365 continue;
366 }
367 IndexType j1 = j;
368 IndexType j2 = j;
369 jNext = j - 1;
370 if (j>1) {
371 if (T(j,j-1)!=Zero) {
372 j1 = j - 1;
373 jNext = j - 2;
374 }
375 }
376
377 if (j1==j2) {
378 //
379 // 1 by 1 diagonal block
380 //
381 // Scale if necessary to avoid overflow in division
382 //
383 ElementType z = w;
384 if (j1==1) {
385 z = b(1);
386 }
387 ElementType xj = abs(x(j1)) + abs(x(n+j1));
388 ElementType Tjj = abs(T(j1,j1)) + abs(z);
389 ElementType tmp = T(j1,j1);
390 if (Tjj<sMinW) {
391 tmp = sMinW;
392 Tjj = sMinW;
393 info = 1;
394 }
395
396 if (xj==Zero) {
397 continue;
398 }
399
400 if (Tjj<One) {
401 if (xj>bigNum*Tjj) {
402 const ElementType rec = One / xj;
403 x *= rec;
404 scale *= rec;
405 xMax *= rec;
406 }
407 }
408 ElementType sr, si;
409 ladiv(x(j1), x(n+j1), tmp, z, sr, si);
410 x(j1) = sr;
411 x(n+j1) = si;
412 xj = abs(x(j1)) + abs(x(n+j1));
413 //
414 // Scale x if necessary to avoid overflow when adding a
415 // multiple of column j1 of T.
416 //
417 if (xj>One) {
418 const ElementType rec = One / xj;
419 if (work(j1)>(bigNum-xMax)*rec) {
420 x *= rec;
421 scale *= rec;
422 }
423 }
424
425 if (j1>1) {
426 x(_(1,j1-1)) -= x(j1) * T(_(1,j1-1),j1);
427 x(_(n+1,n+j1-1)) -= x(n+j1) * T(_(1,j1-1),j1);
428
429 x(1) += b(j1) * x(n+j1);
430 x(n+1) -= b(j1) * x(j1);
431
432 xMax = Zero;
433 for (IndexType k=1; k<=j1-1; ++k) {
434 xMax = max(xMax, abs(x(k)) + abs(x(k+n)));
435 }
436 }
437
438 } else {
439 //
440 // Meet 2 by 2 diagonal block
441 //
442 D(1,1) = x(j1);
443 D(2,1) = x(j2);
444 D(1,2) = x(n+j1);
445 D(2,2) = x(n+j2);
446
447 ElementType scaleC;
448 IndexType iErr = laln2(false, IndexType(2), sMinW, One,
449 T(_(j1,j1+1),_(j1,j1+1)), One, One,
450 D, Zero, -w, V, scaleC, xNorm);
451 if (iErr!=0) {
452 info = 2;
453 }
454
455 if (scaleC!=One) {
456 x *= scaleC;
457 scale *= scaleC;
458 }
459 x(j1) = V(1,1);
460 x(j2) = V(2,1);
461 x(n+j1) = V(1,2);
462 x(n+j2) = V(2,2);
463 //
464 // Scale X(J1), .... to avoid overflow in
465 // updating right hand side.
466 //
467 const ElementType xj = max(abs(V(1,1))+abs(V(1,2)),
468 abs(V(2,1))+abs(V(2,2)));
469 if (xj>One) {
470 const ElementType rec = One / xj;
471 if (max(work(j1), work(j2))>(bigNum-xMax)*rec) {
472 x *= rec;
473 scale *= rec;
474 }
475 }
476 //
477 // Update the right-hand side.
478 //
479 if (j1>1) {
480 x(_(1,j1-1)) -= x(j1) * T(_(1,j1-1),j1);
481 x(_(1,j1-1)) -= x(j2) * T(_(1,j1-1),j2);
482
483 x(_(n+1,n+j1-1)) -= x(n+j1) * T(_(1,j1-1),j1);
484 x(_(n+1,n+j1-1)) -= x(n+j2) * T(_(1,j1-1),j2);
485
486 x(1) += b(j1)*x(n+j1);
487 x(1) += b(j2)*x(n+j2);
488 x(n+1) -= b(j1)*x(j1);
489 x(n+1) -= b(j2)*x(j2);
490
491 xMax = Zero;
492 for (IndexType k=1; k<=j1-1; ++k) {
493 xMax = max(abs(x(k)) + abs(x(k+n)), xMax);
494 }
495 }
496
497 }
498 }
499
500 } else {
501 //
502 // Solve (T + iB)**T*(p+iq) = c+id
503 //
504 IndexType jNext = 1;
505 for (IndexType j=1; j<=n; ++j) {
506 if (j<jNext) {
507 continue;
508 }
509 IndexType j1 = j;
510 IndexType j2 = j;
511 jNext = j + 1;
512 if (j<n) {
513 if (T(j+1,j)!=Zero) {
514 j2 = j + 1;
515 jNext = j + 2;
516 }
517 }
518
519 if (j1==j2) {
520 //
521 // 1 by 1 diagonal block
522 //
523 // Scale if necessary to avoid overflow in forming the
524 // right-hand side element by inner product.
525 //
526 ElementType xj = abs(x(j1)) + abs(x(j1+n));
527 if (xMax>One) {
528 const ElementType rec = One / xMax;
529 if (work(j1)>(bigNum-xj)*rec) {
530 x *= rec;
531 scale *= rec;
532 xMax *= rec;
533 }
534 }
535
536 x(j1) -= T(_(1,j1-1),j1) * x(_(1,j1-1));
537 x(n+j1) -= T(_(1,j1-1),j1) * x(_(n+1,n+j1-1));
538 if (j1>1) {
539 x(j1) -= b(j1) * x(n+1);
540 x(n+j1) += b(j1) * x(1);
541 }
542 xj = abs(x(j1)) + abs(x(j1+n));
543
544 ElementType z = w;
545 if (j1==1) {
546 z = b(1);
547 }
548 //
549 // Scale if necessary to avoid overflow in
550 // complex division
551 //
552 ElementType Tjj = abs(T(j1,j1)) + abs(z);
553 ElementType tmp = T(j1,j1);
554 if (Tjj<sMinW) {
555 tmp = sMinW;
556 Tjj = sMinW;
557 info = 1;
558 }
559 if (Tjj<One) {
560 if (xj>bigNum*Tjj) {
561 const ElementType rec = One / xj;
562 x *= rec;
563 scale *= rec;
564 xMax *= rec;
565 }
566 }
567 ElementType sr, si;
568 ladiv(x(j1), x(n+j1), tmp, -z, sr, si);
569 x(j1) = sr;
570 x(j1+n) = si;
571 xMax = max(abs(x(j1)) + abs(x(j1+n)), xMax);
572
573 } else {
574 //
575 // 2 by 2 diagonal block
576 //
577 // Scale if necessary to avoid overflow in forming the
578 // right-hand side element by inner product.
579 //
580 const ElementType xj = max(abs(x(j1)) + abs(x(n+j1)),
581 abs(x(j2)) + abs(x(n+j2)));
582 if (xMax>One) {
583 const ElementType rec = One / xMax;
584 if (max(work(j1), work(j2))>(bigNum-xj)/xj) {
585 x *= rec;
586 scale *= rec;
587 xMax *= rec;
588 }
589 }
590
591 D(1,1) = x(j1) - T(_(1,j1-1),j1) * x(_(1,j1-1));
592 D(2,1) = x(j2) - T(_(1,j1-1),j2) * x(_(1,j1-1));
593 D(1,2) = x(n+j1) - T(_(1,j1-1),j1) * x(_(n+1,n+j1-1));
594 D(2,2) = x(n+j2) - T(_(1,j1-1),j2) * x(_(n+1,n+j1-1));
595 D(1,1) -= b(j1) * x(n+1);
596 D(2,1) -= b(j2) * x(n+1);
597 D(1,2) += b(j1) * x(1);
598 D(2,2) += b(j2) * x(1);
599
600 ElementType scaleC;
601 IndexType iErr = laln2(true, IndexType(2), sMinW, One,
602 T(_(j1,j1+1),_(j1,j1+1)), One, One,
603 D, Zero, w, V, scaleC, xNorm);
604 if (iErr!=0) {
605 info = 2;
606 }
607
608 if (scaleC!=One) {
609 x *= scaleC;
610 scale *= scaleC;
611 }
612 x(j1) = V(1,1);
613 x(j2) = V(2,1);
614 x(n+j1) = V(1,2);
615 x(n+j2) = V(2,2);
616 xMax = max(abs(x(j1))+abs(x(n+j1)),
617 abs(x(j2))+abs(x(n+j2)),
618 xMax);
619
620 }
621
622 }
623
624 }
625
626 }
627
628 return info;
629 }
630
631 //== interface for native lapack ===============================================
632
633 #ifdef CHECK_CXXLAPACK
634
635 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
636 typename VWORK>
637 typename DenseVector<VX>::IndexType
638 laqtr_native(bool trans,
639 bool real,
640 const GeMatrix<MT> &T,
641 const DenseVector<VB> &b,
642 const W &w,
643 SCALE &scale,
644 DenseVector<VX> &x,
645 DenseVector<VWORK> &work)
646 {
647 typedef typename DenseVector<VX>::ElementType ElementType;
648
649 const LOGICAL LTRAN = trans;
650 const LOGICAL LREAL = real;
651 const INTEGER N = T.numRows();
652 const INTEGER LDT = T.leadingDimension();
653 const DOUBLE _W = w;
654 DOUBLE _SCALE = scale;
655 INTEGER INFO;
656
657 if (IsSame<ElementType, DOUBLE>::value) {
658 LAPACK_IMPL(dlaqtr)(<RAN,
659 &LREAL,
660 &N,
661 T.data(),
662 &LDT,
663 b.data(),
664 &_W,
665 &_SCALE,
666 x.data(),
667 work.data(),
668 &INFO);
669 } else {
670 ASSERT(0);
671 }
672 ASSERT(INFO>=0);
673
674 scale = _SCALE;
675
676 return INFO;
677 }
678
679 #endif // CHECK_CXXLAPACK
680
681 //== public interface ==========================================================
682 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
683 typename VWORK>
684 typename DenseVector<VX>::IndexType
685 laqtr(bool trans,
686 bool real,
687 const GeMatrix<MT> &T,
688 const DenseVector<VB> &b,
689 const W &w,
690 SCALE &scale,
691 DenseVector<VX> &x,
692 DenseVector<VWORK> &work)
693 {
694 typedef typename DenseVector<VX>::IndexType IndexType;
695
696 # ifndef NDEBUG
697 //
698 // Test the input parameters
699 //
700 ASSERT(T.firstRow()==1);
701 ASSERT(T.firstCol()==1);
702 ASSERT(T.numRows()==T.numCols());
703
704 const IndexType n = T.numRows();
705
706 if (!real) {
707 ASSERT(b.firstIndex()==1);
708 ASSERT(b.length()==n);
709 }
710
711 ASSERT(x.length()==2*n);
712 ASSERT(work.length()==n);
713 # endif
714
715 # ifdef CHECK_CXXLAPACK
716 //
717 // Make copies of output arguments
718 //
719 SCALE scale_org = scale;
720 typename DenseVector<VX>::NoView x_org = x;
721 typename DenseVector<VWORK>::NoView work_org = work;
722 # endif
723
724 IndexType info = laqtr_generic(trans, real, T, b, w, scale, x, work);
725
726 # ifdef CHECK_CXXLAPACK
727 //
728 // Make copies of results computed by the generic implementation
729 //
730 SCALE scale_generic = scale;
731 typename DenseVector<VX>::NoView x_generic = x;
732 typename DenseVector<VWORK>::NoView work_generic = work;
733 //
734 // restore output arguments
735 //
736 scale = scale_org;
737 x = x_org;
738 work = work_org;
739 //
740 // Compare results
741 //
742 IndexType _info = laqtr_native(trans, real, T, b, w, scale, x, work);
743
744 bool failed = false;
745 if (! isIdentical(scale_generic, scale, "scale_generic", "scale")) {
746 std::cerr << "CXXLAPACK: scale_generic = "
747 << scale_generic << std::endl;
748 std::cerr << "F77LAPACK: scale = " << scale << std::endl;
749 failed = true;
750 }
751
752 if (! isIdentical(x_generic, x, "x_generic", "x")) {
753 std::cerr << "CXXLAPACK: x_generic = " << x_generic << std::endl;
754 std::cerr << "F77LAPACK: x = " << x << std::endl;
755 failed = true;
756 }
757
758 if (! isIdentical(work_generic, work, "work_generic", "work")) {
759 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
760 std::cerr << "F77LAPACK: work = " << work << std::endl;
761 failed = true;
762 }
763
764 if (! isIdentical(info, _info, "info", "_info")) {
765 std::cerr << "CXXLAPACK: info= " << info<< std::endl;
766 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
767 failed = true;
768 }
769
770 if (failed) {
771 ASSERT(0);
772 }
773 # endif
774
775 return info;
776 }
777
778 //-- forwarding ----------------------------------------------------------------
779 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
780 typename VWORK>
781 typename VX::IndexType
782 laqtr(bool trans,
783 bool real,
784 const MT &T,
785 const VB &b,
786 const W &w,
787 SCALE &scale,
788 VX &x,
789 VWORK &work)
790 {
791 typedef typename VX::IndexType IndexType;
792
793 CHECKPOINT_ENTER;
794 const IndexType info = laqtr(trans, real, T, b, w, scale, x, work);
795 CHECKPOINT_LEAVE;
796
797 return info;
798 }
799
800 } } // namespace lapack, flens
801
802 #endif // FLENS_LAPACK_EIG_LAQTR_TCC
2 * Copyright (c) 2011, Michael Lehn
3 *
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * 1) Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2) Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in
14 * the documentation and/or other materials provided with the
15 * distribution.
16 * 3) Neither the name of the FLENS development group nor the names of
17 * its contributors may be used to endorse or promote products derived
18 * from this software without specific prior written permission.
19 *
20 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
23 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
24 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
25 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
26 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
27 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
28 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
29 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 */
32
33 /* Based on
34 SUBROUTINE DLAQTR( LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK,
35 $ INFO )
36 *
37 * -- LAPACK auxiliary routine (version 3.3.1) --
38 * -- LAPACK is a software package provided by Univ. of Tennessee, --
39 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
40 * -- April 2011 --
41 */
42
43 #ifndef FLENS_LAPACK_EIG_LAQTR_TCC
44 #define FLENS_LAPACK_EIG_LAQTR_TCC 1
45
46 #include <flens/aux/aux.h>
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53
54 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
55 typename VWORK>
56 typename DenseVector<VX>::IndexType
57 laqtr_generic(bool trans,
58 bool real,
59 const GeMatrix<MT> &T,
60 const DenseVector<VB> &b,
61 const W &w,
62 SCALE &scale,
63 DenseVector<VX> &x,
64 DenseVector<VWORK> &work)
65 {
66 using std::abs;
67 using flens::max;
68
69 typedef typename DenseVector<VX>::ElementType ElementType;
70 typedef typename DenseVector<VX>::IndexType IndexType;
71
72 const Underscore<IndexType> _;
73 const IndexType n = T.numRows();
74
75 const ElementType Zero(0), One(1);
76 //
77 // .. Local Arrays ..
78 //
79 ElementType _dataD[4], _dataV[4];
80 GeMatrixView<ElementType>
81 D = typename GeMatrixView<ElementType>::Engine(2, 2, _dataD, 2),
82 V = typename GeMatrixView<ElementType>::Engine(2, 2, _dataV, 2);
83
84 IndexType info = 0;
85 //
86 // Quick return if possible
87 //
88 if (n==0) {
89 return info;
90 }
91 //
92 // Set constants to control overflow
93 //
94 const ElementType eps = lamch<ElementType>(Precision);
95 const ElementType smallNum = lamch<ElementType>(SafeMin)/eps;
96 const ElementType bigNum = One/smallNum;
97
98 ElementType xNorm = lan(MaximumNorm, T);
99 if (!real) {
100 typedef typename DenseVector<VB>::ElementType TB;
101 const GeMatrixConstView<TB> B(n, 1, b, n);
102
103 xNorm = max(xNorm, abs(w), lan(MaximumNorm, B));
104 }
105 const ElementType sMin = max(smallNum, eps*xNorm);
106 //
107 // Compute 1-norm of each column of strictly upper triangular
108 // part of T to control overflow in triangular solver.
109 //
110 work(1) = Zero;
111 for (IndexType j=2; j<=n; ++j) {
112 work(j) = blas::asum(T(_(1,j-1),j));
113 }
114
115 if (!real) {
116 for (IndexType i=2; i<=n; ++i) {
117 work(i) += abs(b(i));
118 }
119 }
120
121 const IndexType n2 = 2*n;
122 const IndexType n1 = (real) ? n : n2;
123
124 const IndexType k = blas::iamax(x(_(1,n1)));
125 ElementType xMax = abs(x(k));
126 scale = One;
127
128 if (xMax>bigNum) {
129 scale = bigNum / xMax;
130 x(_(1,n1)) *= scale;
131 xMax = bigNum;
132 }
133
134 if (real) {
135
136 if (!trans) {
137 //
138 // Solve T*p = scale*c
139 //
140 IndexType jNext = n;
141 for (IndexType j=n; j>=1; --j) {
142 if (j>jNext) {
143 continue;
144 }
145 IndexType j1 = j;
146 IndexType j2 = j;
147 jNext = j - 1;
148 if (j>1) {
149 if (T(j,j-1)!=Zero) {
150 j1 = j - 1;
151 jNext = j - 2;
152 }
153 }
154
155 if (j1==j2) {
156 //
157 // Meet 1 by 1 diagonal block
158 //
159 // Scale to avoid overflow when computing
160 // x(j) = b(j)/T(j,j)
161 //
162 ElementType xj = abs(x(j1));
163 ElementType Tjj = abs(T(j1,j1));
164 ElementType tmp = T(j1,j1);
165 if (Tjj<sMin) {
166 tmp = sMin;
167 Tjj = sMin;
168 info = 1;
169 }
170
171 if (xj==Zero) {
172 continue;
173 }
174
175 if (Tjj<One) {
176 if (xj>bigNum*Tjj) {
177 const ElementType rec = One / xj;
178 x(_(1,n)) *= rec;
179 scale *= rec;
180 xMax *= rec;
181 }
182 }
183 x(j1) /= tmp;
184 xj = abs(x(j1));
185 //
186 // Scale x if necessary to avoid overflow when adding a
187 // multiple of column j1 of T.
188 //
189 if (xj>One) {
190 const ElementType rec = One / xj;
191 if (work(j1)>(bigNum-xMax)*rec) {
192 x(_(1,n)) *= rec;
193 scale *= rec;
194 }
195 }
196 if (j1>1) {
197 auto _x = x(_(1,j1-1));
198
199 _x -= x(j1)*T(_(1,j1-1),j1);
200 const IndexType k = blas::iamax(_x);
201 xMax = abs(x(k));
202 }
203
204 } else {
205 //
206 // Meet 2 by 2 diagonal block
207 //
208 // Call 2 by 2 linear system solve, to take
209 // care of possible overflow by scaling factor.
210 //
211 D(1,1) = x(j1);
212 D(2,1) = x(j2);
213
214 ElementType scaleC;
215 IndexType iErr = laln2(false, IndexType(1), sMin, One,
216 T(_(j1,j1+1),_(j1,j1+1)), One, One,
217 D(_(1,2),_(1,1)), Zero, Zero,
218 V(_(1,2),_(1,1)), scaleC, xNorm);
219 if (iErr!=0) {
220 info = 2;
221 }
222
223 if (scaleC!=One) {
224 x(_(1,n)) *= scaleC;
225 scale *= scaleC;
226 }
227 x(j1) = V(1,1);
228 x(j2) = V(2,1);
229 //
230 // Scale V(1,1) (= X(J1)) and/or V(2,1) (=X(J2))
231 // to avoid overflow in updating right-hand side.
232 //
233 const ElementType xj = max(abs(V(1,1)), abs(V(2,1)));
234 if (xj>One) {
235 const ElementType rec = One / xj;
236 if (max(work(j1),work(j2))>(bigNum-xMax)*rec) {
237 x(_(1,n)) *= rec;
238 scale *= rec;
239 }
240 }
241 //
242 // Update right-hand side
243 //
244 if (j1>1) {
245 auto _x = x(_(1,j1-1));
246 _x -= x(j1) * T(_(1,j1-1),j1);
247 _x -= x(j2) * T(_(1,j1-1),j2);
248 const IndexType k = blas::iamax(_x);
249 xMax = abs(x(k));
250 }
251
252 }
253
254 }
255
256 } else {
257 //
258 // Solve T**T*p = scale*c
259 //
260 IndexType jNext = 1;
261 for (IndexType j=1; j<=n; ++j) {
262 if (j<jNext) {
263 continue;
264 }
265 IndexType j1 = j;
266 IndexType j2 = j;
267 jNext = j + 1;
268 if (j<n) {
269 if (T(j+1,j)!=Zero) {
270 j2 = j + 1;
271 jNext = j + 2;
272 }
273 }
274
275 if (j1==j2) {
276 //
277 // 1 by 1 diagonal block
278 //
279 // Scale if necessary to avoid overflow in forming the
280 // right-hand side element by inner product.
281 //
282 ElementType xj = abs(x(j1));
283 if (xMax>One) {
284 const ElementType rec = One / xMax;
285 if (work(j1)>(bigNum-xj)*rec) {
286 x(_(1,n)) *= rec;
287 scale *= rec;
288 xMax *= rec;
289 }
290 }
291
292 x(j1) -= T(_(1,j1-1),j1) * x(_(1,j1-1));
293
294 xj = abs(x(j1));
295 ElementType Tjj = abs(T(j1,j1));
296 ElementType tmp = T(j1,j1);
297 if (Tjj<sMin) {
298 tmp = sMin;
299 Tjj = sMin;
300 info = 1;
301 }
302
303 if (Tjj<One) {
304 if (xj>bigNum*Tjj) {
305 const ElementType rec = One / xj;
306 x(_(1,n)) *= rec;
307 scale *= rec;
308 xMax *= rec;
309 }
310 }
311 x(j1) /= tmp;
312 xMax = max(xMax, abs(x(j1)));
313
314 } else {
315 //
316 // 2 by 2 diagonal block
317 //
318 // Scale if necessary to avoid overflow in forming the
319 // right-hand side elements by inner product.
320 //
321 const ElementType xj = max(abs(x(j1)), abs(x(j2)));
322 if (xMax>One) {
323 const ElementType rec = One / xMax;
324 if (max(work(j2),work(j1))>(bigNum-xj)*rec) {
325 x(_(1,n)) *= rec;
326 scale *= rec;
327 xMax *= rec;
328 }
329 }
330
331 D(1,1) = x(j1) - T(_(1,j1-1),j1) * x(_(1,j1-1));
332 D(2,1) = x(j2) - T(_(1,j1-1),j2) * x(_(1,j1-1));
333
334 ElementType scaleC;
335 IndexType iErr = laln2(true, IndexType(1), sMin, One,
336 T(_(j1,j1+1),_(j1,j1+1)), One, One,
337 D(_(1,2),_(1,1)), Zero, Zero,
338 V(_(1,2),_(1,1)), scaleC, xNorm);
339 if (iErr!=0) {
340 info = 2;
341 }
342
343 if (scaleC!=One) {
344 x(_(1,n)) *= scaleC;
345 scale *= scaleC;
346 }
347 x(j1) = V(1,1);
348 x(j2) = V(2,1);
349 xMax = max(abs(x(j1)), abs(x(j2)), xMax);
350
351 }
352 }
353 }
354
355 } else {
356
357 const ElementType sMinW = max(eps*abs(w), sMin);
358 if (!trans) {
359 //
360 // Solve (T + iB)*(p+iq) = c+id
361 //
362 IndexType jNext = n;
363 for (IndexType j=n; j>=1; --j) {
364 if (j>jNext) {
365 continue;
366 }
367 IndexType j1 = j;
368 IndexType j2 = j;
369 jNext = j - 1;
370 if (j>1) {
371 if (T(j,j-1)!=Zero) {
372 j1 = j - 1;
373 jNext = j - 2;
374 }
375 }
376
377 if (j1==j2) {
378 //
379 // 1 by 1 diagonal block
380 //
381 // Scale if necessary to avoid overflow in division
382 //
383 ElementType z = w;
384 if (j1==1) {
385 z = b(1);
386 }
387 ElementType xj = abs(x(j1)) + abs(x(n+j1));
388 ElementType Tjj = abs(T(j1,j1)) + abs(z);
389 ElementType tmp = T(j1,j1);
390 if (Tjj<sMinW) {
391 tmp = sMinW;
392 Tjj = sMinW;
393 info = 1;
394 }
395
396 if (xj==Zero) {
397 continue;
398 }
399
400 if (Tjj<One) {
401 if (xj>bigNum*Tjj) {
402 const ElementType rec = One / xj;
403 x *= rec;
404 scale *= rec;
405 xMax *= rec;
406 }
407 }
408 ElementType sr, si;
409 ladiv(x(j1), x(n+j1), tmp, z, sr, si);
410 x(j1) = sr;
411 x(n+j1) = si;
412 xj = abs(x(j1)) + abs(x(n+j1));
413 //
414 // Scale x if necessary to avoid overflow when adding a
415 // multiple of column j1 of T.
416 //
417 if (xj>One) {
418 const ElementType rec = One / xj;
419 if (work(j1)>(bigNum-xMax)*rec) {
420 x *= rec;
421 scale *= rec;
422 }
423 }
424
425 if (j1>1) {
426 x(_(1,j1-1)) -= x(j1) * T(_(1,j1-1),j1);
427 x(_(n+1,n+j1-1)) -= x(n+j1) * T(_(1,j1-1),j1);
428
429 x(1) += b(j1) * x(n+j1);
430 x(n+1) -= b(j1) * x(j1);
431
432 xMax = Zero;
433 for (IndexType k=1; k<=j1-1; ++k) {
434 xMax = max(xMax, abs(x(k)) + abs(x(k+n)));
435 }
436 }
437
438 } else {
439 //
440 // Meet 2 by 2 diagonal block
441 //
442 D(1,1) = x(j1);
443 D(2,1) = x(j2);
444 D(1,2) = x(n+j1);
445 D(2,2) = x(n+j2);
446
447 ElementType scaleC;
448 IndexType iErr = laln2(false, IndexType(2), sMinW, One,
449 T(_(j1,j1+1),_(j1,j1+1)), One, One,
450 D, Zero, -w, V, scaleC, xNorm);
451 if (iErr!=0) {
452 info = 2;
453 }
454
455 if (scaleC!=One) {
456 x *= scaleC;
457 scale *= scaleC;
458 }
459 x(j1) = V(1,1);
460 x(j2) = V(2,1);
461 x(n+j1) = V(1,2);
462 x(n+j2) = V(2,2);
463 //
464 // Scale X(J1), .... to avoid overflow in
465 // updating right hand side.
466 //
467 const ElementType xj = max(abs(V(1,1))+abs(V(1,2)),
468 abs(V(2,1))+abs(V(2,2)));
469 if (xj>One) {
470 const ElementType rec = One / xj;
471 if (max(work(j1), work(j2))>(bigNum-xMax)*rec) {
472 x *= rec;
473 scale *= rec;
474 }
475 }
476 //
477 // Update the right-hand side.
478 //
479 if (j1>1) {
480 x(_(1,j1-1)) -= x(j1) * T(_(1,j1-1),j1);
481 x(_(1,j1-1)) -= x(j2) * T(_(1,j1-1),j2);
482
483 x(_(n+1,n+j1-1)) -= x(n+j1) * T(_(1,j1-1),j1);
484 x(_(n+1,n+j1-1)) -= x(n+j2) * T(_(1,j1-1),j2);
485
486 x(1) += b(j1)*x(n+j1);
487 x(1) += b(j2)*x(n+j2);
488 x(n+1) -= b(j1)*x(j1);
489 x(n+1) -= b(j2)*x(j2);
490
491 xMax = Zero;
492 for (IndexType k=1; k<=j1-1; ++k) {
493 xMax = max(abs(x(k)) + abs(x(k+n)), xMax);
494 }
495 }
496
497 }
498 }
499
500 } else {
501 //
502 // Solve (T + iB)**T*(p+iq) = c+id
503 //
504 IndexType jNext = 1;
505 for (IndexType j=1; j<=n; ++j) {
506 if (j<jNext) {
507 continue;
508 }
509 IndexType j1 = j;
510 IndexType j2 = j;
511 jNext = j + 1;
512 if (j<n) {
513 if (T(j+1,j)!=Zero) {
514 j2 = j + 1;
515 jNext = j + 2;
516 }
517 }
518
519 if (j1==j2) {
520 //
521 // 1 by 1 diagonal block
522 //
523 // Scale if necessary to avoid overflow in forming the
524 // right-hand side element by inner product.
525 //
526 ElementType xj = abs(x(j1)) + abs(x(j1+n));
527 if (xMax>One) {
528 const ElementType rec = One / xMax;
529 if (work(j1)>(bigNum-xj)*rec) {
530 x *= rec;
531 scale *= rec;
532 xMax *= rec;
533 }
534 }
535
536 x(j1) -= T(_(1,j1-1),j1) * x(_(1,j1-1));
537 x(n+j1) -= T(_(1,j1-1),j1) * x(_(n+1,n+j1-1));
538 if (j1>1) {
539 x(j1) -= b(j1) * x(n+1);
540 x(n+j1) += b(j1) * x(1);
541 }
542 xj = abs(x(j1)) + abs(x(j1+n));
543
544 ElementType z = w;
545 if (j1==1) {
546 z = b(1);
547 }
548 //
549 // Scale if necessary to avoid overflow in
550 // complex division
551 //
552 ElementType Tjj = abs(T(j1,j1)) + abs(z);
553 ElementType tmp = T(j1,j1);
554 if (Tjj<sMinW) {
555 tmp = sMinW;
556 Tjj = sMinW;
557 info = 1;
558 }
559 if (Tjj<One) {
560 if (xj>bigNum*Tjj) {
561 const ElementType rec = One / xj;
562 x *= rec;
563 scale *= rec;
564 xMax *= rec;
565 }
566 }
567 ElementType sr, si;
568 ladiv(x(j1), x(n+j1), tmp, -z, sr, si);
569 x(j1) = sr;
570 x(j1+n) = si;
571 xMax = max(abs(x(j1)) + abs(x(j1+n)), xMax);
572
573 } else {
574 //
575 // 2 by 2 diagonal block
576 //
577 // Scale if necessary to avoid overflow in forming the
578 // right-hand side element by inner product.
579 //
580 const ElementType xj = max(abs(x(j1)) + abs(x(n+j1)),
581 abs(x(j2)) + abs(x(n+j2)));
582 if (xMax>One) {
583 const ElementType rec = One / xMax;
584 if (max(work(j1), work(j2))>(bigNum-xj)/xj) {
585 x *= rec;
586 scale *= rec;
587 xMax *= rec;
588 }
589 }
590
591 D(1,1) = x(j1) - T(_(1,j1-1),j1) * x(_(1,j1-1));
592 D(2,1) = x(j2) - T(_(1,j1-1),j2) * x(_(1,j1-1));
593 D(1,2) = x(n+j1) - T(_(1,j1-1),j1) * x(_(n+1,n+j1-1));
594 D(2,2) = x(n+j2) - T(_(1,j1-1),j2) * x(_(n+1,n+j1-1));
595 D(1,1) -= b(j1) * x(n+1);
596 D(2,1) -= b(j2) * x(n+1);
597 D(1,2) += b(j1) * x(1);
598 D(2,2) += b(j2) * x(1);
599
600 ElementType scaleC;
601 IndexType iErr = laln2(true, IndexType(2), sMinW, One,
602 T(_(j1,j1+1),_(j1,j1+1)), One, One,
603 D, Zero, w, V, scaleC, xNorm);
604 if (iErr!=0) {
605 info = 2;
606 }
607
608 if (scaleC!=One) {
609 x *= scaleC;
610 scale *= scaleC;
611 }
612 x(j1) = V(1,1);
613 x(j2) = V(2,1);
614 x(n+j1) = V(1,2);
615 x(n+j2) = V(2,2);
616 xMax = max(abs(x(j1))+abs(x(n+j1)),
617 abs(x(j2))+abs(x(n+j2)),
618 xMax);
619
620 }
621
622 }
623
624 }
625
626 }
627
628 return info;
629 }
630
631 //== interface for native lapack ===============================================
632
633 #ifdef CHECK_CXXLAPACK
634
635 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
636 typename VWORK>
637 typename DenseVector<VX>::IndexType
638 laqtr_native(bool trans,
639 bool real,
640 const GeMatrix<MT> &T,
641 const DenseVector<VB> &b,
642 const W &w,
643 SCALE &scale,
644 DenseVector<VX> &x,
645 DenseVector<VWORK> &work)
646 {
647 typedef typename DenseVector<VX>::ElementType ElementType;
648
649 const LOGICAL LTRAN = trans;
650 const LOGICAL LREAL = real;
651 const INTEGER N = T.numRows();
652 const INTEGER LDT = T.leadingDimension();
653 const DOUBLE _W = w;
654 DOUBLE _SCALE = scale;
655 INTEGER INFO;
656
657 if (IsSame<ElementType, DOUBLE>::value) {
658 LAPACK_IMPL(dlaqtr)(<RAN,
659 &LREAL,
660 &N,
661 T.data(),
662 &LDT,
663 b.data(),
664 &_W,
665 &_SCALE,
666 x.data(),
667 work.data(),
668 &INFO);
669 } else {
670 ASSERT(0);
671 }
672 ASSERT(INFO>=0);
673
674 scale = _SCALE;
675
676 return INFO;
677 }
678
679 #endif // CHECK_CXXLAPACK
680
681 //== public interface ==========================================================
682 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
683 typename VWORK>
684 typename DenseVector<VX>::IndexType
685 laqtr(bool trans,
686 bool real,
687 const GeMatrix<MT> &T,
688 const DenseVector<VB> &b,
689 const W &w,
690 SCALE &scale,
691 DenseVector<VX> &x,
692 DenseVector<VWORK> &work)
693 {
694 typedef typename DenseVector<VX>::IndexType IndexType;
695
696 # ifndef NDEBUG
697 //
698 // Test the input parameters
699 //
700 ASSERT(T.firstRow()==1);
701 ASSERT(T.firstCol()==1);
702 ASSERT(T.numRows()==T.numCols());
703
704 const IndexType n = T.numRows();
705
706 if (!real) {
707 ASSERT(b.firstIndex()==1);
708 ASSERT(b.length()==n);
709 }
710
711 ASSERT(x.length()==2*n);
712 ASSERT(work.length()==n);
713 # endif
714
715 # ifdef CHECK_CXXLAPACK
716 //
717 // Make copies of output arguments
718 //
719 SCALE scale_org = scale;
720 typename DenseVector<VX>::NoView x_org = x;
721 typename DenseVector<VWORK>::NoView work_org = work;
722 # endif
723
724 IndexType info = laqtr_generic(trans, real, T, b, w, scale, x, work);
725
726 # ifdef CHECK_CXXLAPACK
727 //
728 // Make copies of results computed by the generic implementation
729 //
730 SCALE scale_generic = scale;
731 typename DenseVector<VX>::NoView x_generic = x;
732 typename DenseVector<VWORK>::NoView work_generic = work;
733 //
734 // restore output arguments
735 //
736 scale = scale_org;
737 x = x_org;
738 work = work_org;
739 //
740 // Compare results
741 //
742 IndexType _info = laqtr_native(trans, real, T, b, w, scale, x, work);
743
744 bool failed = false;
745 if (! isIdentical(scale_generic, scale, "scale_generic", "scale")) {
746 std::cerr << "CXXLAPACK: scale_generic = "
747 << scale_generic << std::endl;
748 std::cerr << "F77LAPACK: scale = " << scale << std::endl;
749 failed = true;
750 }
751
752 if (! isIdentical(x_generic, x, "x_generic", "x")) {
753 std::cerr << "CXXLAPACK: x_generic = " << x_generic << std::endl;
754 std::cerr << "F77LAPACK: x = " << x << std::endl;
755 failed = true;
756 }
757
758 if (! isIdentical(work_generic, work, "work_generic", "work")) {
759 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
760 std::cerr << "F77LAPACK: work = " << work << std::endl;
761 failed = true;
762 }
763
764 if (! isIdentical(info, _info, "info", "_info")) {
765 std::cerr << "CXXLAPACK: info= " << info<< std::endl;
766 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
767 failed = true;
768 }
769
770 if (failed) {
771 ASSERT(0);
772 }
773 # endif
774
775 return info;
776 }
777
778 //-- forwarding ----------------------------------------------------------------
779 template <typename MT, typename VB, typename W, typename SCALE, typename VX,
780 typename VWORK>
781 typename VX::IndexType
782 laqtr(bool trans,
783 bool real,
784 const MT &T,
785 const VB &b,
786 const W &w,
787 SCALE &scale,
788 VX &x,
789 VWORK &work)
790 {
791 typedef typename VX::IndexType IndexType;
792
793 CHECKPOINT_ENTER;
794 const IndexType info = laqtr(trans, real, T, b, w, scale, x, work);
795 CHECKPOINT_LEAVE;
796
797 return info;
798 }
799
800 } } // namespace lapack, flens
801
802 #endif // FLENS_LAPACK_EIG_LAQTR_TCC