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