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 *
35 SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
36 $ WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
37 $ IWORK, LIWORK, BWORK, INFO )
38 *
39 * -- LAPACK driver routine (version 3.2.2) --
40 * -- LAPACK is a software package provided by Univ. of Tennessee, --
41 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
42 * June 2010
43 */
44
45 #ifndef FLENS_LAPACK_EIG_ESX_TCC
46 #define FLENS_LAPACK_EIG_ESX_TCC 1
47
48 #include <cmath>
49 #include <flens/blas/blas.h>
50 #include <flens/lapack/lapack.h>
51
52 namespace flens { namespace lapack {
53
54 //== generic lapack implementation =============================================
55
56 template <typename MA>
57 Quadruple<typename GeMatrix<MA>::IndexType>
58 esx_generic_wsq(bool computeSchurVectors,
59 SENSE::Sense sense,
60 const GeMatrix<MA> &A)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MA>::ElementType T;
65 typedef typename GeMatrix<MA>::IndexType IndexType;
66
67 const IndexType n = A.numRows();
68 const bool wantSN = (sense==SENSE::None);
69 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
70 const bool wantSB = (sense==SENSE::Both);
71
72 IndexType minWork, maxWork, liWork, lWork;
73
74 liWork = 1;
75 if (n==0) {
76 minWork = 1;
77 lWork = 1;
78 } else {
79 maxWork = 2*n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
80 minWork = 3*n;
81
82 HSEQR::ComputeZ computeZ = (computeSchurVectors) ? HSEQR::NoInit
83 : HSEQR::No;
84 IndexType hsWork = hseqr_wsq(HSEQR::Schur, computeZ,
85 IndexType(1), n, A);
86
87 if (!computeSchurVectors) {
88 maxWork = max(maxWork, n + hsWork);
89 } else {
90 maxWork = max(maxWork,
91 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
92 maxWork = max(maxWork, n +hsWork);
93 }
94 lWork = maxWork;
95 if (!wantSN) {
96 lWork = max(lWork, n + (n*n)/2);
97 }
98 if (wantSV || wantSB) {
99 liWork = (n*n)/4;
100 }
101
102 }
103 return Quadruple<IndexType>(minWork, maxWork, lWork, liWork);
104 }
105
106 template <typename SelectFunction, typename MA, typename IndexType,
107 typename VWR, typename VWI, typename MVS,
108 typename RCONDE, typename RCONDV,
109 typename VWORK, typename VIWORK, typename BWORK>
110 IndexType
111 esx_generic(bool computeSchurVectors,
112 bool sortEigenvalues,
113 SelectFunction selectFunction,
114 SENSE::Sense sense,
115 GeMatrix<MA> &A,
116 IndexType &sDim,
117 DenseVector<VWR> &wr,
118 DenseVector<VWI> &wi,
119 GeMatrix<MVS> &VS,
120 RCONDE &rCondE,
121 RCONDV &rCondV,
122 DenseVector<VWORK> &work,
123 DenseVector<VIWORK> &iWork,
124 DenseVector<BWORK> &bWork)
125 {
126 using std::sqrt;
127
128 typedef typename GeMatrix<MA>::ElementType T;
129
130 const Underscore<IndexType> _;
131 const IndexType n = A.numRows();
132 const T Zero(0), One(1);
133
134 //
135 // .. Local Arrays ..
136 // this array is used to save variables between calls to lacn2
137 //
138 T _dum[1];
139 DenseVectorView<T> dum = typename DenseVectorView<T>::Engine(1, _dum);
140 //
141 // Test the input arguments
142 //
143 IndexType info = 0;
144 const bool wantSN = (sense==SENSE::None);
145 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
146 const bool wantSB = (sense==SENSE::Both);
147 //
148 // Compute workspace
149 // (Note: Comments in the code beginning "RWorkspace:" describe the
150 // minimal amount of real workspace needed at that point in the
151 // code, as well as the preferred amount for good performance.
152 // IWorkspace refers to integer workspace.
153 // NB refers to the optimal block size for the immediately
154 // following subroutine, as returned by ILAENV.
155 // HSWORK refers to the workspace preferred by DHSEQR, as
156 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
157 // the worst case.
158 // If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
159 // depends on SDIM, which is computed by the routine DTRSEN later
160 // in the code.)
161 //
162 Quadruple<IndexType> wsQuery = esx_wsq(computeSchurVectors, sense, A);
163 IndexType minWork = wsQuery.first;
164 IndexType maxWork = wsQuery.second;
165 IndexType lWork = wsQuery.third;
166 IndexType liWork = wsQuery.forth;
167
168 if (work.length()!=0 && work.length()<minWork) {
169 ASSERT(0);
170 } else if (work.length()==0) {
171 work.resize(lWork);
172 }
173 work(1) = lWork;
174
175 if (iWork.length()!=0 && iWork.length()<liWork) {
176 ASSERT(0);
177 } else if (work.length()==0) {
178 iWork.resize(liWork);
179 }
180 iWork(1) = liWork;
181
182 //
183 // Quick return if possible
184 //
185 if (n==0) {
186 sDim = 0;
187 return info;
188 }
189 //
190 // Get machine constants
191 //
192 const T eps = lamch<T>(Precision);
193 T smallNum = lamch<T>(SafeMin);
194 T bigNum = One / smallNum;
195 labad(smallNum, bigNum);
196
197 smallNum = sqrt(smallNum) / eps;
198 bigNum = One / smallNum;
199 //
200 // Scale A if max element outside range [SMLNUM,BIGNUM]
201 //
202 const T normA = lan(MaximumNorm, A);
203
204 bool scaleA = false;
205 T cScale;
206
207 if (normA>Zero && normA<smallNum) {
208 scaleA = true;
209 cScale = smallNum;
210 } else if (normA>bigNum) {
211 scaleA = true;
212 cScale = bigNum;
213 }
214 if (scaleA) {
215 lascl(LASCL::FullMatrix, 0, 0, normA, cScale, A);
216 }
217
218 //
219 // Permute the matrix to make it more nearly triangular
220 // (RWorkspace: need N)
221 //
222 IndexType iBal = 1;
223 auto balWork = work(_(iBal,iBal+n-1));
224
225 IndexType iLo, iHi;
226 bal(BALANCE::PermuteOnly, A, iLo, iHi, balWork);
227 //
228 // Reduce to upper Hessenberg form
229 // (RWorkspace: need 3*N, prefer 2*N+N*NB)
230 //
231 IndexType iTau = iBal + n;
232 IndexType iWrk = iTau + n;
233 lWork = work.length();
234
235 auto tau = work(_(iTau, iTau+n-2));
236 auto hrdWork = work(_(iWrk, lWork));
237
238 hrd(iLo, iHi, A, tau, hrdWork);
239
240 if (computeSchurVectors) {
241 //
242 // Copy Householder vectors to VS
243 //
244 VS.lower() = A.lower();
245 //
246 // Generate orthogonal matrix in VS
247 // (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
248 //
249 orghr(iLo, iHi, VS, tau, hrdWork);
250 }
251
252 sDim = 0;
253 //
254 // Perform QR iteration, accumulating Schur vectors in VS if desired
255 // (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
256 //
257 iWrk = iTau;
258 auto hseqrWork = work(_(iWrk, lWork));
259 HSEQR::ComputeZ computeZ = (computeSchurVectors) ? HSEQR::NoInit
260 : HSEQR::No;
261 IndexType iEval = hseqr(HSEQR::Schur, computeZ, iLo, iHi,
262 A, wr, wi, VS, hseqrWork);
263 if (iEval>0) {
264 info = iEval;
265 }
266 //
267 // Sort eigenvalues if desired
268 //
269 if (sortEigenvalues && info==0) {
270 if (scaleA) {
271 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wr);
272 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi);
273 }
274 for (IndexType i=1; i<=n; ++i) {
275 bWork(i) = selectFunction(wr(i), wi(i));
276 }
277 //
278 // Reorder eigenvalues, transform Schur vectors, and compute
279 // reciprocal condition numbers
280 // (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
281 // otherwise, need N )
282 // (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
283 // otherwise, need 0 )
284 //
285 IndexType iCond = trsen(TRSEN::Job(sense), computeSchurVectors, bWork,
286 A, VS, wr, wi, sDim, rCondE, rCondV,
287 hseqrWork, iWork);
288 if (!wantSN) {
289 maxWork = max(maxWork, n+2*sDim*(n-sDim));
290 }
291 if (iCond>0) {
292 //
293 // DTRSEN failed to reorder or to restore standard Schur form
294 //
295 info = iCond + n;
296 }
297 }
298
299 if (computeSchurVectors) {
300 //
301 // Undo balancing
302 // (Workspace: need N)
303 //
304 bak(BALANCE::PermuteOnly, Right, iLo, iHi, balWork, VS);
305 }
306
307 if (scaleA) {
308 //
309 // Undo scaling for the Schur form of A
310 //
311 lascl(LASCL::UpperHessenberg, 0, 0, cScale, normA, A);
312 wr = A.diag(0);
313 if ((wantSV || wantSB) && info==0) {
314 dum(1) = rCondV;
315 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, dum);
316 rCondV = dum(1);
317 }
318 if (cScale==smallNum) {
319 //
320 // If scaling back towards underflow, adjust WI if an
321 // offdiagonal element of a 2-by-2 block in the Schur form
322 // underflows.
323 //
324 IndexType i1, i2;
325 if (iEval>0) {
326 i1 = iEval + 1;
327 i2 = iHi - 1;
328 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi(_(1,iLo-1)));
329 } else if (sortEigenvalues) {
330 i1 = 1;
331 i2 = n - 1;
332 } else {
333 i1 = iLo;
334 i2 = iHi - 1;
335 }
336 IndexType iNext = i1-1;
337 for (IndexType i=i1; i<=i2; ++i) {
338 if (i<iNext) {
339 continue;
340 }
341 if (wi(i)==Zero) {
342 iNext = i + 1;
343 } else {
344 if (A(i+1,i)==Zero) {
345 wi(i) = Zero;
346 wi(i+1) = Zero;
347 } else if (A(i+1,i)!=Zero && A(i,i+1)==Zero) {
348 wi(i) = Zero;
349 wi(i+1) = Zero;
350 if (i>1) {
351 blas::swap(A(_(1,i-1),i), A(_(1,i-1),i+1));
352 }
353 if (n>i+1) {
354 blas::swap(A(i,_(i+2,n)), A(i+1,_(i+2,n)));
355 }
356 blas::swap(VS(_,i), VS(_,i+1));
357 A(i,i+1) = A(i+1,i);
358 A(i+1,i) = Zero;
359 }
360 iNext = i + 2;
361 }
362 }
363 }
364 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi(_(iEval+1,n)));
365 }
366
367 if (sortEigenvalues && info==0) {
368 //
369 // Check if reordering successful
370 //
371 bool lastSelect = true;
372 bool last2Select = true;
373 sDim = 0;
374 IndexType ip = 0;
375 for (IndexType i=1; i<=n; ++i) {
376 bool currentSelect = selectFunction(wr(i), wi(i));
377 if (wi(i)==Zero) {
378 if (currentSelect) {
379 ++sDim;
380 }
381 ip = 0;
382 if (currentSelect && !lastSelect) {
383 info = n + 2;
384 }
385 } else {
386 if (ip==1) {
387 //
388 // Last eigenvalue of conjugate pair
389 //
390 currentSelect = currentSelect || lastSelect;
391 lastSelect = currentSelect;
392 if (currentSelect) {
393 sDim += 2;
394 }
395 ip = -1;
396 if (currentSelect && !last2Select) {
397 info = n + 2;
398 }
399 } else {
400 //
401 // First eigenvalue of conjugate pair
402 //
403 ip = 1;
404 }
405 }
406 last2Select = lastSelect;
407 lastSelect = currentSelect;
408 }
409 }
410
411 work(1) = maxWork;
412 if (wantSV || wantSB) {
413 iWork(1) = max(IndexType(1), sDim*(n-sDim));
414 } else {
415 iWork(1) = 1;
416 }
417
418 return info;
419 }
420
421 //== interface for native lapack ===============================================
422
423 #ifdef CHECK_CXXLAPACK
424
425 template <typename SelectFunction, typename MA, typename IndexType,
426 typename VWR, typename VWI, typename MVS,
427 typename RCONDE, typename RCONDV,
428 typename VWORK, typename VIWORK, typename BWORK>
429 IndexType
430 esx_native(bool computeSchurVectors,
431 bool sortEigenvalues,
432 SelectFunction selectFunction,
433 SENSE::Sense sense,
434 GeMatrix<MA> &A,
435 IndexType &sDim,
436 DenseVector<VWR> &wr,
437 DenseVector<VWI> &wi,
438 GeMatrix<MVS> &VS,
439 RCONDE &rCondE,
440 RCONDV &rCondV,
441 DenseVector<VWORK> &work,
442 DenseVector<VIWORK> &iWork,
443 DenseVector<BWORK> &bWork)
444 {
445 typedef typename GeMatrix<MA>::ElementType T;
446
447 const char JOBVS = (computeSchurVectors) ? 'V' : 'N';
448 const char SORT = (sortEigenvalues) ? 'S' : 'N';
449 const char SENSE = sense;
450 LOGICAL (*SELECT)(const T*, const T*) = selectFunction.select;
451 const INTEGER N = A.numRows();
452 const INTEGER LDA = A.leadingDimension();
453 INTEGER SDIM = sDim;
454 const INTEGER LDVS = VS.leadingDimension();
455 T _RCONDE = rCondE;
456 T _RCONDV = rCondV;
457 const INTEGER LWORK = work.length();
458 const INTEGER LIWORK = iWork.length();
459 INTEGER INFO;
460
461 DenseVector<Array<INTEGER> > _iWork(LIWORK);
462 for (IndexType i=1; i<=LIWORK; ++i) {
463 _iWork(i) = iWork(i);
464 }
465 DenseVector<Array<LOGICAL> > _bWork(N);
466 for (IndexType i=1; i<=N; ++i) {
467 _bWork(i) = bWork(i);
468 }
469
470 if (IsSame<T,DOUBLE>::value) {
471 LAPACK_IMPL(dgeesx)(&JOBVS,
472 &SORT,
473 SELECT,
474 &SENSE,
475 &N,
476 A.data(),
477 &LDA,
478 &SDIM,
479 wr.data(),
480 wi.data(),
481 VS.data(),
482 &LDVS,
483 &_RCONDE,
484 &_RCONDV,
485 work.data(),
486 &LWORK,
487 _iWork.data(),
488 &LIWORK,
489 _bWork.data(),
490 &INFO);
491 } else {
492 ASSERT(0);
493 }
494 if (INFO<0) {
495 std::cerr << "dgeesx: INFO = " << INFO << std::endl;
496 }
497 ASSERT(INFO>=0);
498
499 sDim = SDIM;
500 rCondE = _RCONDE;
501 rCondV = _RCONDV;
502
503 for (IndexType i=1; i<=LIWORK; ++i) {
504 iWork(i) = _iWork(i);
505 }
506 for (IndexType i=1; i<=N; ++i) {
507 bWork(i) = _bWork(i);
508 }
509
510 return INFO;
511 }
512
513 #endif // CHECK_CXXLAPACK
514
515 //== public interface ==========================================================
516 template <typename MA>
517 Quadruple<typename GeMatrix<MA>::IndexType>
518 esx_wsq(bool computeSchurVectors,
519 SENSE::Sense sense,
520 const GeMatrix<MA> &A)
521 {
522 return esx_generic_wsq(computeSchurVectors, sense, A);
523 }
524
525 template <typename SelectFunction, typename MA, typename IndexType,
526 typename VWR, typename VWI, typename MVS,
527 typename RCONDE, typename RCONDV,
528 typename VWORK, typename VIWORK, typename BWORK>
529 IndexType
530 esx(bool computeSchurVectors,
531 bool sortEigenvalues,
532 SelectFunction selectFunction,
533 SENSE::Sense sense,
534 GeMatrix<MA> &A,
535 IndexType &sDim,
536 DenseVector<VWR> &wr,
537 DenseVector<VWI> &wi,
538 GeMatrix<MVS> &VS,
539 RCONDE &rCondE,
540 RCONDV &rCondV,
541 DenseVector<VWORK> &work,
542 DenseVector<VIWORK> &iWork,
543 DenseVector<BWORK> &bWork)
544 {
545 LAPACK_DEBUG_OUT("esx");
546
547 //
548 // Test the input parameters
549 //
550 # ifndef NDEBUG
551 ASSERT(A.numRows()==A.numCols());
552 ASSERT(A.firstRow()==1);
553 ASSERT(A.firstCol()==1);
554 ASSERT(work.firstIndex()==1);
555 ASSERT(iWork.firstIndex()==1);
556 ASSERT(bWork.firstIndex()==1);
557
558 typename GeMatrix<MA>::IndexType n = A.numRows();
559
560 ASSERT(wr.firstIndex()==1);
561 ASSERT(wr.length()==n);
562
563 ASSERT(wi.firstIndex()==1);
564 ASSERT(wi.length()==n);
565
566 if (computeSchurVectors) {
567 ASSERT(VS.numRows()==VS.numCols());
568 ASSERT(VS.numRows()==n);
569 ASSERT(VS.firstRow()==1);
570 ASSERT(VS.firstCol()==1);
571 }
572
573 if (sortEigenvalues) {
574 ASSERT(bWork.firstIndex()==1);
575 ASSERT(bWork.length()>=n);
576 }
577 # endif
578
579 //
580 // Make copies of output arguments
581 //
582 # ifdef CHECK_CXXLAPACK
583
584 typename GeMatrix<MA>::NoView A_org = A;
585 IndexType sDim_org = sDim;
586 typename DenseVector<VWR>::NoView wr_org = wr;
587 typename DenseVector<VWI>::NoView wi_org = wi;
588 typename GeMatrix<MVS>::NoView VS_org = VS;
589 RCONDE rCondE_org = rCondE;
590 RCONDV rCondV_org = rCondV;
591 typename DenseVector<VWORK>::NoView work_org = work;
592 typename DenseVector<VIWORK>::NoView iWork_org = iWork;
593 typename DenseVector<BWORK>::NoView bWork_org = bWork;
594
595 # endif
596
597 //
598 // Call implementation
599 //
600 IndexType result = esx_generic(computeSchurVectors,
601 sortEigenvalues,
602 selectFunction,
603 sense,
604 A,
605 sDim,
606 wr,
607 wi,
608 VS,
609 rCondE,
610 rCondV,
611 work,
612 iWork,
613 bWork);
614
615 # ifdef CHECK_CXXLAPACK
616 //
617 // Make copies of results computed by the generic implementation
618 //
619 typename GeMatrix<MA>::NoView A_generic = A;
620 IndexType sDim_generic = sDim;
621 typename DenseVector<VWR>::NoView wr_generic = wr;
622 typename DenseVector<VWI>::NoView wi_generic = wi;
623 typename GeMatrix<MVS>::NoView VS_generic = VS;
624 RCONDE rCondE_generic = rCondE;
625 RCONDV rCondV_generic = rCondV;
626 typename DenseVector<VWORK>::NoView work_generic = work;
627 typename DenseVector<VIWORK>::NoView iWork_generic = iWork;
628 typename DenseVector<BWORK>::NoView bWork_generic = bWork;
629
630 //
631 // restore output arguments
632 //
633 A = A_org;
634 sDim = sDim_org;
635 wr = wr_org;
636 wi = wi_org;
637 VS = VS_org;
638 rCondE = rCondE_org;
639 rCondV = rCondV_org;
640 work = work_org;
641 iWork = iWork_org;
642 bWork = bWork_org;
643
644 //
645 // Compare generic results with results from the native implementation
646 //
647 IndexType _result = esx_native(computeSchurVectors,
648 sortEigenvalues,
649 selectFunction,
650 sense,
651 A,
652 sDim,
653 wr,
654 wi,
655 VS,
656 rCondE,
657 rCondV,
658 work,
659 iWork,
660 bWork);
661
662 bool failed = false;
663 if (! isIdentical(A_generic, A, "A_generic", "A")) {
664 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
665 std::cerr << "F77LAPACK: A = " << A << std::endl;
666 failed = true;
667 }
668
669 if (! isIdentical(sDim_generic, sDim, "sDim_generic", "sDim")) {
670 std::cerr << "CXXLAPACK: sDim_generic = " << sDim_generic << std::endl;
671 std::cerr << "F77LAPACK: sDim = " << sDim << std::endl;
672 failed = true;
673 }
674
675 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
676 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
677 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
678 failed = true;
679 }
680
681 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
682 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
683 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
684 failed = true;
685 }
686
687 if (! isIdentical(VS_generic, VS, "VS_generic", "VS")) {
688 std::cerr << "CXXLAPACK: VS_generic = " << VS_generic << std::endl;
689 std::cerr << "F77LAPACK: VS = " << VS << std::endl;
690 failed = true;
691 }
692
693 if (! isIdentical(rCondE_generic, rCondE, "rCondE_generic", "rCondE")) {
694 std::cerr << "CXXLAPACK: rCondE_generic = "
695 << rCondE_generic << std::endl;
696 std::cerr << "F77LAPACK: rCondE = " << rCondE << std::endl;
697 failed = true;
698 }
699
700 if (! isIdentical(rCondV_generic, rCondV, "rCondV_generic", "rCondV")) {
701 std::cerr << "CXXLAPACK: rCondV_generic = "
702 << rCondV_generic << std::endl;
703 std::cerr << "F77LAPACK: rCondV = " << rCondV << std::endl;
704 failed = true;
705 }
706
707 if (! isIdentical(work_generic, work, "work_generic", "work")) {
708 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
709 std::cerr << "F77LAPACK: work = " << work << std::endl;
710 failed = true;
711 }
712
713 if (! isIdentical(iWork_generic, iWork, "iWork_generic", "iWork")) {
714 std::cerr << "CXXLAPACK: iWork_generic = "
715 << iWork_generic << std::endl;
716 std::cerr << "F77LAPACK: iWork = " << iWork << std::endl;
717 failed = true;
718 }
719
720 if (! isIdentical(bWork_generic, bWork, "bWork_generic", "bWork")) {
721 std::cerr << "CXXLAPACK: bWork_generic = "
722 << bWork_generic << std::endl;
723 std::cerr << "F77LAPACK: bWork = " << bWork << std::endl;
724 failed = true;
725 }
726
727 if (! isIdentical(result, _result, " result", "_result")) {
728 std::cerr << "CXXLAPACK: result = " << result << std::endl;
729 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
730 failed = true;
731 }
732
733 if (failed) {
734 ASSERT(0);
735 } else {
736 // std::cerr << "passed: esx.tcc" << std::endl;
737 }
738
739 # endif
740
741 return result;
742 }
743
744 //-- forwarding ----------------------------------------------------------------
745
746 } } // namespace lapack, flens
747
748 #endif // FLENS_LAPACK_EIG_ESX_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 *
35 SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
36 $ WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
37 $ IWORK, LIWORK, BWORK, INFO )
38 *
39 * -- LAPACK driver routine (version 3.2.2) --
40 * -- LAPACK is a software package provided by Univ. of Tennessee, --
41 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
42 * June 2010
43 */
44
45 #ifndef FLENS_LAPACK_EIG_ESX_TCC
46 #define FLENS_LAPACK_EIG_ESX_TCC 1
47
48 #include <cmath>
49 #include <flens/blas/blas.h>
50 #include <flens/lapack/lapack.h>
51
52 namespace flens { namespace lapack {
53
54 //== generic lapack implementation =============================================
55
56 template <typename MA>
57 Quadruple<typename GeMatrix<MA>::IndexType>
58 esx_generic_wsq(bool computeSchurVectors,
59 SENSE::Sense sense,
60 const GeMatrix<MA> &A)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MA>::ElementType T;
65 typedef typename GeMatrix<MA>::IndexType IndexType;
66
67 const IndexType n = A.numRows();
68 const bool wantSN = (sense==SENSE::None);
69 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
70 const bool wantSB = (sense==SENSE::Both);
71
72 IndexType minWork, maxWork, liWork, lWork;
73
74 liWork = 1;
75 if (n==0) {
76 minWork = 1;
77 lWork = 1;
78 } else {
79 maxWork = 2*n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
80 minWork = 3*n;
81
82 HSEQR::ComputeZ computeZ = (computeSchurVectors) ? HSEQR::NoInit
83 : HSEQR::No;
84 IndexType hsWork = hseqr_wsq(HSEQR::Schur, computeZ,
85 IndexType(1), n, A);
86
87 if (!computeSchurVectors) {
88 maxWork = max(maxWork, n + hsWork);
89 } else {
90 maxWork = max(maxWork,
91 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
92 maxWork = max(maxWork, n +hsWork);
93 }
94 lWork = maxWork;
95 if (!wantSN) {
96 lWork = max(lWork, n + (n*n)/2);
97 }
98 if (wantSV || wantSB) {
99 liWork = (n*n)/4;
100 }
101
102 }
103 return Quadruple<IndexType>(minWork, maxWork, lWork, liWork);
104 }
105
106 template <typename SelectFunction, typename MA, typename IndexType,
107 typename VWR, typename VWI, typename MVS,
108 typename RCONDE, typename RCONDV,
109 typename VWORK, typename VIWORK, typename BWORK>
110 IndexType
111 esx_generic(bool computeSchurVectors,
112 bool sortEigenvalues,
113 SelectFunction selectFunction,
114 SENSE::Sense sense,
115 GeMatrix<MA> &A,
116 IndexType &sDim,
117 DenseVector<VWR> &wr,
118 DenseVector<VWI> &wi,
119 GeMatrix<MVS> &VS,
120 RCONDE &rCondE,
121 RCONDV &rCondV,
122 DenseVector<VWORK> &work,
123 DenseVector<VIWORK> &iWork,
124 DenseVector<BWORK> &bWork)
125 {
126 using std::sqrt;
127
128 typedef typename GeMatrix<MA>::ElementType T;
129
130 const Underscore<IndexType> _;
131 const IndexType n = A.numRows();
132 const T Zero(0), One(1);
133
134 //
135 // .. Local Arrays ..
136 // this array is used to save variables between calls to lacn2
137 //
138 T _dum[1];
139 DenseVectorView<T> dum = typename DenseVectorView<T>::Engine(1, _dum);
140 //
141 // Test the input arguments
142 //
143 IndexType info = 0;
144 const bool wantSN = (sense==SENSE::None);
145 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
146 const bool wantSB = (sense==SENSE::Both);
147 //
148 // Compute workspace
149 // (Note: Comments in the code beginning "RWorkspace:" describe the
150 // minimal amount of real workspace needed at that point in the
151 // code, as well as the preferred amount for good performance.
152 // IWorkspace refers to integer workspace.
153 // NB refers to the optimal block size for the immediately
154 // following subroutine, as returned by ILAENV.
155 // HSWORK refers to the workspace preferred by DHSEQR, as
156 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
157 // the worst case.
158 // If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
159 // depends on SDIM, which is computed by the routine DTRSEN later
160 // in the code.)
161 //
162 Quadruple<IndexType> wsQuery = esx_wsq(computeSchurVectors, sense, A);
163 IndexType minWork = wsQuery.first;
164 IndexType maxWork = wsQuery.second;
165 IndexType lWork = wsQuery.third;
166 IndexType liWork = wsQuery.forth;
167
168 if (work.length()!=0 && work.length()<minWork) {
169 ASSERT(0);
170 } else if (work.length()==0) {
171 work.resize(lWork);
172 }
173 work(1) = lWork;
174
175 if (iWork.length()!=0 && iWork.length()<liWork) {
176 ASSERT(0);
177 } else if (work.length()==0) {
178 iWork.resize(liWork);
179 }
180 iWork(1) = liWork;
181
182 //
183 // Quick return if possible
184 //
185 if (n==0) {
186 sDim = 0;
187 return info;
188 }
189 //
190 // Get machine constants
191 //
192 const T eps = lamch<T>(Precision);
193 T smallNum = lamch<T>(SafeMin);
194 T bigNum = One / smallNum;
195 labad(smallNum, bigNum);
196
197 smallNum = sqrt(smallNum) / eps;
198 bigNum = One / smallNum;
199 //
200 // Scale A if max element outside range [SMLNUM,BIGNUM]
201 //
202 const T normA = lan(MaximumNorm, A);
203
204 bool scaleA = false;
205 T cScale;
206
207 if (normA>Zero && normA<smallNum) {
208 scaleA = true;
209 cScale = smallNum;
210 } else if (normA>bigNum) {
211 scaleA = true;
212 cScale = bigNum;
213 }
214 if (scaleA) {
215 lascl(LASCL::FullMatrix, 0, 0, normA, cScale, A);
216 }
217
218 //
219 // Permute the matrix to make it more nearly triangular
220 // (RWorkspace: need N)
221 //
222 IndexType iBal = 1;
223 auto balWork = work(_(iBal,iBal+n-1));
224
225 IndexType iLo, iHi;
226 bal(BALANCE::PermuteOnly, A, iLo, iHi, balWork);
227 //
228 // Reduce to upper Hessenberg form
229 // (RWorkspace: need 3*N, prefer 2*N+N*NB)
230 //
231 IndexType iTau = iBal + n;
232 IndexType iWrk = iTau + n;
233 lWork = work.length();
234
235 auto tau = work(_(iTau, iTau+n-2));
236 auto hrdWork = work(_(iWrk, lWork));
237
238 hrd(iLo, iHi, A, tau, hrdWork);
239
240 if (computeSchurVectors) {
241 //
242 // Copy Householder vectors to VS
243 //
244 VS.lower() = A.lower();
245 //
246 // Generate orthogonal matrix in VS
247 // (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
248 //
249 orghr(iLo, iHi, VS, tau, hrdWork);
250 }
251
252 sDim = 0;
253 //
254 // Perform QR iteration, accumulating Schur vectors in VS if desired
255 // (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
256 //
257 iWrk = iTau;
258 auto hseqrWork = work(_(iWrk, lWork));
259 HSEQR::ComputeZ computeZ = (computeSchurVectors) ? HSEQR::NoInit
260 : HSEQR::No;
261 IndexType iEval = hseqr(HSEQR::Schur, computeZ, iLo, iHi,
262 A, wr, wi, VS, hseqrWork);
263 if (iEval>0) {
264 info = iEval;
265 }
266 //
267 // Sort eigenvalues if desired
268 //
269 if (sortEigenvalues && info==0) {
270 if (scaleA) {
271 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wr);
272 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi);
273 }
274 for (IndexType i=1; i<=n; ++i) {
275 bWork(i) = selectFunction(wr(i), wi(i));
276 }
277 //
278 // Reorder eigenvalues, transform Schur vectors, and compute
279 // reciprocal condition numbers
280 // (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
281 // otherwise, need N )
282 // (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
283 // otherwise, need 0 )
284 //
285 IndexType iCond = trsen(TRSEN::Job(sense), computeSchurVectors, bWork,
286 A, VS, wr, wi, sDim, rCondE, rCondV,
287 hseqrWork, iWork);
288 if (!wantSN) {
289 maxWork = max(maxWork, n+2*sDim*(n-sDim));
290 }
291 if (iCond>0) {
292 //
293 // DTRSEN failed to reorder or to restore standard Schur form
294 //
295 info = iCond + n;
296 }
297 }
298
299 if (computeSchurVectors) {
300 //
301 // Undo balancing
302 // (Workspace: need N)
303 //
304 bak(BALANCE::PermuteOnly, Right, iLo, iHi, balWork, VS);
305 }
306
307 if (scaleA) {
308 //
309 // Undo scaling for the Schur form of A
310 //
311 lascl(LASCL::UpperHessenberg, 0, 0, cScale, normA, A);
312 wr = A.diag(0);
313 if ((wantSV || wantSB) && info==0) {
314 dum(1) = rCondV;
315 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, dum);
316 rCondV = dum(1);
317 }
318 if (cScale==smallNum) {
319 //
320 // If scaling back towards underflow, adjust WI if an
321 // offdiagonal element of a 2-by-2 block in the Schur form
322 // underflows.
323 //
324 IndexType i1, i2;
325 if (iEval>0) {
326 i1 = iEval + 1;
327 i2 = iHi - 1;
328 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi(_(1,iLo-1)));
329 } else if (sortEigenvalues) {
330 i1 = 1;
331 i2 = n - 1;
332 } else {
333 i1 = iLo;
334 i2 = iHi - 1;
335 }
336 IndexType iNext = i1-1;
337 for (IndexType i=i1; i<=i2; ++i) {
338 if (i<iNext) {
339 continue;
340 }
341 if (wi(i)==Zero) {
342 iNext = i + 1;
343 } else {
344 if (A(i+1,i)==Zero) {
345 wi(i) = Zero;
346 wi(i+1) = Zero;
347 } else if (A(i+1,i)!=Zero && A(i,i+1)==Zero) {
348 wi(i) = Zero;
349 wi(i+1) = Zero;
350 if (i>1) {
351 blas::swap(A(_(1,i-1),i), A(_(1,i-1),i+1));
352 }
353 if (n>i+1) {
354 blas::swap(A(i,_(i+2,n)), A(i+1,_(i+2,n)));
355 }
356 blas::swap(VS(_,i), VS(_,i+1));
357 A(i,i+1) = A(i+1,i);
358 A(i+1,i) = Zero;
359 }
360 iNext = i + 2;
361 }
362 }
363 }
364 lascl(LASCL::FullMatrix, 0, 0, cScale, normA, wi(_(iEval+1,n)));
365 }
366
367 if (sortEigenvalues && info==0) {
368 //
369 // Check if reordering successful
370 //
371 bool lastSelect = true;
372 bool last2Select = true;
373 sDim = 0;
374 IndexType ip = 0;
375 for (IndexType i=1; i<=n; ++i) {
376 bool currentSelect = selectFunction(wr(i), wi(i));
377 if (wi(i)==Zero) {
378 if (currentSelect) {
379 ++sDim;
380 }
381 ip = 0;
382 if (currentSelect && !lastSelect) {
383 info = n + 2;
384 }
385 } else {
386 if (ip==1) {
387 //
388 // Last eigenvalue of conjugate pair
389 //
390 currentSelect = currentSelect || lastSelect;
391 lastSelect = currentSelect;
392 if (currentSelect) {
393 sDim += 2;
394 }
395 ip = -1;
396 if (currentSelect && !last2Select) {
397 info = n + 2;
398 }
399 } else {
400 //
401 // First eigenvalue of conjugate pair
402 //
403 ip = 1;
404 }
405 }
406 last2Select = lastSelect;
407 lastSelect = currentSelect;
408 }
409 }
410
411 work(1) = maxWork;
412 if (wantSV || wantSB) {
413 iWork(1) = max(IndexType(1), sDim*(n-sDim));
414 } else {
415 iWork(1) = 1;
416 }
417
418 return info;
419 }
420
421 //== interface for native lapack ===============================================
422
423 #ifdef CHECK_CXXLAPACK
424
425 template <typename SelectFunction, typename MA, typename IndexType,
426 typename VWR, typename VWI, typename MVS,
427 typename RCONDE, typename RCONDV,
428 typename VWORK, typename VIWORK, typename BWORK>
429 IndexType
430 esx_native(bool computeSchurVectors,
431 bool sortEigenvalues,
432 SelectFunction selectFunction,
433 SENSE::Sense sense,
434 GeMatrix<MA> &A,
435 IndexType &sDim,
436 DenseVector<VWR> &wr,
437 DenseVector<VWI> &wi,
438 GeMatrix<MVS> &VS,
439 RCONDE &rCondE,
440 RCONDV &rCondV,
441 DenseVector<VWORK> &work,
442 DenseVector<VIWORK> &iWork,
443 DenseVector<BWORK> &bWork)
444 {
445 typedef typename GeMatrix<MA>::ElementType T;
446
447 const char JOBVS = (computeSchurVectors) ? 'V' : 'N';
448 const char SORT = (sortEigenvalues) ? 'S' : 'N';
449 const char SENSE = sense;
450 LOGICAL (*SELECT)(const T*, const T*) = selectFunction.select;
451 const INTEGER N = A.numRows();
452 const INTEGER LDA = A.leadingDimension();
453 INTEGER SDIM = sDim;
454 const INTEGER LDVS = VS.leadingDimension();
455 T _RCONDE = rCondE;
456 T _RCONDV = rCondV;
457 const INTEGER LWORK = work.length();
458 const INTEGER LIWORK = iWork.length();
459 INTEGER INFO;
460
461 DenseVector<Array<INTEGER> > _iWork(LIWORK);
462 for (IndexType i=1; i<=LIWORK; ++i) {
463 _iWork(i) = iWork(i);
464 }
465 DenseVector<Array<LOGICAL> > _bWork(N);
466 for (IndexType i=1; i<=N; ++i) {
467 _bWork(i) = bWork(i);
468 }
469
470 if (IsSame<T,DOUBLE>::value) {
471 LAPACK_IMPL(dgeesx)(&JOBVS,
472 &SORT,
473 SELECT,
474 &SENSE,
475 &N,
476 A.data(),
477 &LDA,
478 &SDIM,
479 wr.data(),
480 wi.data(),
481 VS.data(),
482 &LDVS,
483 &_RCONDE,
484 &_RCONDV,
485 work.data(),
486 &LWORK,
487 _iWork.data(),
488 &LIWORK,
489 _bWork.data(),
490 &INFO);
491 } else {
492 ASSERT(0);
493 }
494 if (INFO<0) {
495 std::cerr << "dgeesx: INFO = " << INFO << std::endl;
496 }
497 ASSERT(INFO>=0);
498
499 sDim = SDIM;
500 rCondE = _RCONDE;
501 rCondV = _RCONDV;
502
503 for (IndexType i=1; i<=LIWORK; ++i) {
504 iWork(i) = _iWork(i);
505 }
506 for (IndexType i=1; i<=N; ++i) {
507 bWork(i) = _bWork(i);
508 }
509
510 return INFO;
511 }
512
513 #endif // CHECK_CXXLAPACK
514
515 //== public interface ==========================================================
516 template <typename MA>
517 Quadruple<typename GeMatrix<MA>::IndexType>
518 esx_wsq(bool computeSchurVectors,
519 SENSE::Sense sense,
520 const GeMatrix<MA> &A)
521 {
522 return esx_generic_wsq(computeSchurVectors, sense, A);
523 }
524
525 template <typename SelectFunction, typename MA, typename IndexType,
526 typename VWR, typename VWI, typename MVS,
527 typename RCONDE, typename RCONDV,
528 typename VWORK, typename VIWORK, typename BWORK>
529 IndexType
530 esx(bool computeSchurVectors,
531 bool sortEigenvalues,
532 SelectFunction selectFunction,
533 SENSE::Sense sense,
534 GeMatrix<MA> &A,
535 IndexType &sDim,
536 DenseVector<VWR> &wr,
537 DenseVector<VWI> &wi,
538 GeMatrix<MVS> &VS,
539 RCONDE &rCondE,
540 RCONDV &rCondV,
541 DenseVector<VWORK> &work,
542 DenseVector<VIWORK> &iWork,
543 DenseVector<BWORK> &bWork)
544 {
545 LAPACK_DEBUG_OUT("esx");
546
547 //
548 // Test the input parameters
549 //
550 # ifndef NDEBUG
551 ASSERT(A.numRows()==A.numCols());
552 ASSERT(A.firstRow()==1);
553 ASSERT(A.firstCol()==1);
554 ASSERT(work.firstIndex()==1);
555 ASSERT(iWork.firstIndex()==1);
556 ASSERT(bWork.firstIndex()==1);
557
558 typename GeMatrix<MA>::IndexType n = A.numRows();
559
560 ASSERT(wr.firstIndex()==1);
561 ASSERT(wr.length()==n);
562
563 ASSERT(wi.firstIndex()==1);
564 ASSERT(wi.length()==n);
565
566 if (computeSchurVectors) {
567 ASSERT(VS.numRows()==VS.numCols());
568 ASSERT(VS.numRows()==n);
569 ASSERT(VS.firstRow()==1);
570 ASSERT(VS.firstCol()==1);
571 }
572
573 if (sortEigenvalues) {
574 ASSERT(bWork.firstIndex()==1);
575 ASSERT(bWork.length()>=n);
576 }
577 # endif
578
579 //
580 // Make copies of output arguments
581 //
582 # ifdef CHECK_CXXLAPACK
583
584 typename GeMatrix<MA>::NoView A_org = A;
585 IndexType sDim_org = sDim;
586 typename DenseVector<VWR>::NoView wr_org = wr;
587 typename DenseVector<VWI>::NoView wi_org = wi;
588 typename GeMatrix<MVS>::NoView VS_org = VS;
589 RCONDE rCondE_org = rCondE;
590 RCONDV rCondV_org = rCondV;
591 typename DenseVector<VWORK>::NoView work_org = work;
592 typename DenseVector<VIWORK>::NoView iWork_org = iWork;
593 typename DenseVector<BWORK>::NoView bWork_org = bWork;
594
595 # endif
596
597 //
598 // Call implementation
599 //
600 IndexType result = esx_generic(computeSchurVectors,
601 sortEigenvalues,
602 selectFunction,
603 sense,
604 A,
605 sDim,
606 wr,
607 wi,
608 VS,
609 rCondE,
610 rCondV,
611 work,
612 iWork,
613 bWork);
614
615 # ifdef CHECK_CXXLAPACK
616 //
617 // Make copies of results computed by the generic implementation
618 //
619 typename GeMatrix<MA>::NoView A_generic = A;
620 IndexType sDim_generic = sDim;
621 typename DenseVector<VWR>::NoView wr_generic = wr;
622 typename DenseVector<VWI>::NoView wi_generic = wi;
623 typename GeMatrix<MVS>::NoView VS_generic = VS;
624 RCONDE rCondE_generic = rCondE;
625 RCONDV rCondV_generic = rCondV;
626 typename DenseVector<VWORK>::NoView work_generic = work;
627 typename DenseVector<VIWORK>::NoView iWork_generic = iWork;
628 typename DenseVector<BWORK>::NoView bWork_generic = bWork;
629
630 //
631 // restore output arguments
632 //
633 A = A_org;
634 sDim = sDim_org;
635 wr = wr_org;
636 wi = wi_org;
637 VS = VS_org;
638 rCondE = rCondE_org;
639 rCondV = rCondV_org;
640 work = work_org;
641 iWork = iWork_org;
642 bWork = bWork_org;
643
644 //
645 // Compare generic results with results from the native implementation
646 //
647 IndexType _result = esx_native(computeSchurVectors,
648 sortEigenvalues,
649 selectFunction,
650 sense,
651 A,
652 sDim,
653 wr,
654 wi,
655 VS,
656 rCondE,
657 rCondV,
658 work,
659 iWork,
660 bWork);
661
662 bool failed = false;
663 if (! isIdentical(A_generic, A, "A_generic", "A")) {
664 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
665 std::cerr << "F77LAPACK: A = " << A << std::endl;
666 failed = true;
667 }
668
669 if (! isIdentical(sDim_generic, sDim, "sDim_generic", "sDim")) {
670 std::cerr << "CXXLAPACK: sDim_generic = " << sDim_generic << std::endl;
671 std::cerr << "F77LAPACK: sDim = " << sDim << std::endl;
672 failed = true;
673 }
674
675 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
676 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
677 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
678 failed = true;
679 }
680
681 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
682 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
683 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
684 failed = true;
685 }
686
687 if (! isIdentical(VS_generic, VS, "VS_generic", "VS")) {
688 std::cerr << "CXXLAPACK: VS_generic = " << VS_generic << std::endl;
689 std::cerr << "F77LAPACK: VS = " << VS << std::endl;
690 failed = true;
691 }
692
693 if (! isIdentical(rCondE_generic, rCondE, "rCondE_generic", "rCondE")) {
694 std::cerr << "CXXLAPACK: rCondE_generic = "
695 << rCondE_generic << std::endl;
696 std::cerr << "F77LAPACK: rCondE = " << rCondE << std::endl;
697 failed = true;
698 }
699
700 if (! isIdentical(rCondV_generic, rCondV, "rCondV_generic", "rCondV")) {
701 std::cerr << "CXXLAPACK: rCondV_generic = "
702 << rCondV_generic << std::endl;
703 std::cerr << "F77LAPACK: rCondV = " << rCondV << std::endl;
704 failed = true;
705 }
706
707 if (! isIdentical(work_generic, work, "work_generic", "work")) {
708 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
709 std::cerr << "F77LAPACK: work = " << work << std::endl;
710 failed = true;
711 }
712
713 if (! isIdentical(iWork_generic, iWork, "iWork_generic", "iWork")) {
714 std::cerr << "CXXLAPACK: iWork_generic = "
715 << iWork_generic << std::endl;
716 std::cerr << "F77LAPACK: iWork = " << iWork << std::endl;
717 failed = true;
718 }
719
720 if (! isIdentical(bWork_generic, bWork, "bWork_generic", "bWork")) {
721 std::cerr << "CXXLAPACK: bWork_generic = "
722 << bWork_generic << std::endl;
723 std::cerr << "F77LAPACK: bWork = " << bWork << std::endl;
724 failed = true;
725 }
726
727 if (! isIdentical(result, _result, " result", "_result")) {
728 std::cerr << "CXXLAPACK: result = " << result << std::endl;
729 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
730 failed = true;
731 }
732
733 if (failed) {
734 ASSERT(0);
735 } else {
736 // std::cerr << "passed: esx.tcc" << std::endl;
737 }
738
739 # endif
740
741 return result;
742 }
743
744 //-- forwarding ----------------------------------------------------------------
745
746 } } // namespace lapack, flens
747
748 #endif // FLENS_LAPACK_EIG_ESX_TCC