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 DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
36 $ LDVR, WORK, LWORK, INFO )
37 *
38 * -- LAPACK driver routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_EIG_EVX_TCC
45 #define FLENS_LAPACK_EIG_EVX_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA>
54 Pair<typename GeMatrix<MA>::IndexType>
55 evx_generic_wsq(bool computeVL,
56 bool computeVR,
57 SENSE::Sense sense,
58 GeMatrix<MA> &A)
59 {
60 using std::max;
61
62 typedef typename GeMatrix<MA>::ElementType T;
63 typedef typename GeMatrix<MA>::IndexType IndexType;
64
65 const IndexType n = A.numRows();
66 const bool wantSN = (sense==SENSE::None);
67 const bool wantSE = (sense==SENSE::EigenvaluesOnly);
68
69 IndexType minWork, maxWork;
70
71 if (n==0) {
72 minWork = 1;
73 maxWork = 1;
74 } else {
75 maxWork = n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
76
77 IndexType hseqrWork;
78 if (computeVL) {
79 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
80 IndexType(1), n, A);
81 } else if (computeVL) {
82 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
83 IndexType(1), n, A);
84 } else {
85 if (wantSN) {
86 hseqrWork = hseqr_wsq(HSEQR::Eigenvalues, HSEQR::No,
87 IndexType(1), n, A);
88 } else {
89 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::No,
90 IndexType(1), n, A);
91 }
92 }
93
94 if ((!computeVL) && (!computeVR)) {
95 minWork = 2*n;
96 if (!wantSN) {
97 minWork = max(minWork, n*n+6*n);
98 }
99 maxWork = max(maxWork, hseqrWork);
100 if (!wantSN) {
101 maxWork = max(maxWork, n*n+6*n);
102 }
103 } else {
104 minWork = 3*n;
105 if ((!wantSN) && (!wantSE)) {
106 minWork = max(minWork, n*n + 6*n);
107 }
108 maxWork = max(maxWork, hseqrWork);
109 maxWork = max(maxWork, n+(n-1*ilaenv<T>(1, "ORGHR", "", n, 1, n)));
110 if ((!wantSN) && (!wantSE)) {
111 maxWork = max(maxWork, n*n + 6*n);
112 }
113 maxWork = max(maxWork, 3*n);
114 }
115 maxWork = max(maxWork, minWork);
116 }
117
118 return Pair<typename GeMatrix<MA>::IndexType>(minWork, maxWork);
119 }
120
121 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
122 typename IndexType, typename VSCALE, typename ABNORM,
123 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
124 typename GeMatrix<MA>::IndexType
125 evx_generic(BALANCE::Balance balance,
126 bool computeVL,
127 bool computeVR,
128 SENSE::Sense sense,
129 GeMatrix<MA> &A,
130 DenseVector<VWR> &wr,
131 DenseVector<VWI> &wi,
132 GeMatrix<MVL> &VL,
133 GeMatrix<MVR> &VR,
134 IndexType &iLo,
135 IndexType &iHi,
136 DenseVector<VSCALE> &scale,
137 ABNORM &ABNorm,
138 DenseVector<RCONDE> &rCondE,
139 DenseVector<RCONDV> &rCondV,
140 DenseVector<VWORK> &work,
141 DenseVector<VIWORK> &iWork)
142 {
143 typedef typename GeMatrix<MA>::ElementType T;
144
145 const IndexType n = A.numRows();
146 const Underscore<IndexType> _;
147
148 const T Zero(0), One(1);
149 //
150 // .. Local Arrays ..
151 //
152 bool _selectData[1];
153 DenseVectorView<bool>
154 select = typename DenseVectorView<bool>::Engine(1, _selectData);
155 //
156 // Test the input arguments
157 //
158 IndexType info = 0;
159
160 const bool wantSN = (sense==SENSE::None);
161 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
162 const bool wantSB = (sense==SENSE::Both);
163 //
164 // Compute workspace
165 // (Note: Comments in the code beginning "Workspace:" describe the
166 // minimal amount of workspace needed at that point in the code,
167 // as well as the preferred amount for good performance.
168 // NB refers to the optimal block size for the immediately
169 // following subroutine, as returned by ILAENV.
170 // HSWORK refers to the workspace preferred by DHSEQR, as
171 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
172 // the worst case.)
173 //
174 Pair<IndexType> wsQuery = evx_wsq(computeVL, computeVR, sense, A);
175 IndexType minWork = wsQuery.first;
176 IndexType maxWork = wsQuery.second;
177
178 if (work.length()!=0 && work.length()<minWork) {
179 ASSERT(0);
180 } else if (work.length()==0) {
181 work.resize(maxWork);
182 }
183 work(1) = maxWork;
184
185 //
186 // Quick return if possible
187 //
188 if (n==0) {
189 return info;
190 }
191 //
192 // Get machine constants
193 //
194 const T eps = lamch<T>(Precision);
195 T smallNum = lamch<T>(SafeMin);
196 T bigNum = One / smallNum;
197 labad(smallNum, bigNum);
198 smallNum = sqrt(smallNum) / eps;
199 bigNum = One / smallNum;
200 //
201 // Scale A if max element outside range [SMLNUM,BIGNUM]
202 //
203 IndexType iCond = 0;
204 const T ANorm = lan(MaximumNorm, A);
205 bool scaleA = false;
206 T cScale;
207 if (ANorm>Zero && ANorm<smallNum) {
208 scaleA = true;
209 cScale = smallNum;
210 } else if (ANorm>bigNum) {
211 scaleA = true;
212 cScale = bigNum;
213 }
214 if (scaleA) {
215 lascl(LASCL::FullMatrix, 0, 0, ANorm, cScale, A);
216 }
217 //
218 // Balance the matrix and compute ABNRM
219 //
220 bal(balance, A, iLo, iHi, scale);
221 ABNorm = lan(OneNorm, A);
222 if (scaleA) {
223 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, ABNorm);
224 }
225 //
226 // Reduce to upper Hessenberg form
227 // (Workspace: need 2*N, prefer N+N*NB)
228 //
229 IndexType iTau = 1;
230 IndexType iWrk = iTau + n;
231 IndexType lWork = work.length();
232
233 auto tau = work(_(iTau, iTau+n-2));
234 auto hrdWork = work(_(iWrk, lWork));
235
236 hrd(iLo, iHi, A, tau, hrdWork);
237
238 if (computeVL) {
239 //
240 // Want left eigenvectors
241 // Copy Householder vectors to VL
242 //
243 VL.lower() = A.lower();
244 //
245 // Generate orthogonal matrix in VL
246 // (Workspace: need 2*N-1, prefer N+(N-1)*NB)
247 //
248 orghr(iLo, iHi, VL, tau, hrdWork);
249 //
250 // Perform QR iteration, accumulating Schur vectors in VL
251 // (Workspace: need 1, prefer HSWORK (see comments) )
252 //
253 iWrk = iTau;
254 auto hseqrWork = work(_(iWrk, lWork));
255 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
256 wr, wi, VL, hseqrWork);
257
258 if (computeVR) {
259 //
260 // Want left and right eigenvectors
261 // Copy Schur vectors to VR
262 //
263 VR = VL;
264 }
265
266 } else if (computeVR) {
267 //
268 // Want right eigenvectors
269 // Copy Householder vectors to VR
270 //
271 VR.lower() = A.lower();
272 //
273 // Generate orthogonal matrix in VR
274 // (Workspace: need 2*N-1, prefer N+(N-1)*NB)
275 //
276 orghr(iLo, iHi, VR, tau, hrdWork);
277 //
278 // Perform QR iteration, accumulating Schur vectors in VR
279 // (Workspace: need 1, prefer HSWORK (see comments) )
280 //
281 iWrk = iTau;
282 auto hseqrWork = work(_(iWrk, lWork));
283 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
284 wr, wi, VR, hseqrWork);
285
286 } else {
287 //
288 // Compute eigenvalues only
289 // If condition numbers desired, compute Schur form
290 //
291 HSEQR::Job job = (wantSN) ? HSEQR::Eigenvalues
292 : HSEQR::Schur;
293 //
294 // (Workspace: need 1, prefer HSWORK (see comments) )
295 //
296 iWrk = iTau;
297 auto hseqrWork = work(_(iWrk, lWork));
298 info = hseqr(job, HSEQR::No, iLo, iHi, A, wr, wi, VR, hseqrWork);
299 }
300 //
301 // If INFO > 0 from DHSEQR, then quit
302 //
303 if (info==0) {
304
305 if (computeVL || computeVR) {
306 //
307 // Compute left and/or right eigenvectors
308 // (Workspace: need 3*N)
309 //
310 IndexType nOut;
311 auto trevcWork = work(_(iWrk, lWork));
312
313 trevc(computeVL, computeVR, TREVC::Backtransform, select,
314 A, VL, VR, n, nOut, trevcWork);
315 }
316 //
317 // Compute condition numbers if desired
318 // (Workspace: need N*N+6*N unless SENSE = 'E')
319 //
320 if (!wantSN) {
321 IndexType nOut;
322
323 IndexType _n = (sense!=SENSE::EigenvaluesOnly) ? n : 0;
324 GeMatrixView<T> Work(_n, n+6, work(_(iWrk, lWork)), n);
325
326 trsna(TRSNA::Job(sense), TRSNA::All, select, A, VL, VR,
327 rCondE, rCondV, n, nOut, Work, iWork);
328 }
329
330 if (computeVL) {
331 //
332 // Undo balancing of left eigenvectors
333 //
334 bak(balance, Left, iLo, iHi, scale, VL);
335 //
336 // Normalize left eigenvectors and make largest component real
337 //
338 for (IndexType i=1; i<=n; ++i) {
339 if (wi(i)==Zero) {
340 VL(_,i) *= One / blas::nrm2(VL(_,i));
341 } else if (wi(i)>Zero) {
342 const T scl = One / lapy2(blas::nrm2(VL(_,i)),
343 blas::nrm2(VL(_,i+1)));
344 VL(_,i) *= scl;
345 VL(_,i+1) *= scl;
346 for (IndexType k=1; k<=n; ++k) {
347 work(k) = pow(VL(k,i),2) + pow(VL(k,i+1),2);
348 }
349 IndexType k = blas::iamax(work(_(1,n)));
350 T cs, sn, r;
351 lartg(VL(k,i), VL(k,i+1), cs, sn, r);
352 blas::rot(VL(_,i), VL(_,i+1), cs, sn);
353 VL(k,i+1) = Zero;
354 }
355 }
356 }
357
358 if (computeVR) {
359 //
360 // Undo balancing of right eigenvectors
361 //
362 bak(BALANCE::Both, Right, iLo, iHi, scale, VR);
363 //
364 // Normalize right eigenvectors and make largest component real
365 //
366 for (IndexType i=1; i<=n; ++i) {
367 if (wi(i)==Zero) {
368 VR(_,i) *= One / blas::nrm2(VR(_,i));
369 } else if (wi(i)>Zero) {
370 const T scl = One / lapy2(blas::nrm2(VR(_,i)),
371 blas::nrm2(VR(_,i+1)));
372 VR(_,i) *= scl;
373 VR(_,i+1) *= scl;
374 for (IndexType k=1; k<=n; ++k) {
375 work(k) = pow(VR(k,i), 2) + pow(VR(k,i+1), 2);
376 }
377 IndexType k = blas::iamax(work(_(1,n)));
378 T cs, sn, r;
379 lartg(VR(k,i), VR(k,i+1), cs, sn, r);
380 blas::rot(VR(_,i), VR(_,i+1), cs, sn);
381 VR(k,i+1) = Zero;
382 }
383 }
384 }
385 }
386
387 //
388 // Undo scaling if necessary
389 //
390 if (scaleA) {
391 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(info+1,n)));
392 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(info+1,n)));
393
394 if (info==0) {
395 if ((wantSV || wantSB) && iCond==0) {
396 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, rCondV);
397 }
398 } else {
399 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(1,iLo-1)));
400 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(1,iLo-1)));
401 }
402 }
403
404 work(1) = maxWork;
405 return info;
406 }
407
408 //== interface for native lapack ===============================================
409
410 #ifdef CHECK_CXXLAPACK
411
412 template <typename MA>
413 Pair<typename GeMatrix<MA>::IndexType>
414 evx_native_wsq(bool computeVL,
415 bool computeVR,
416 SENSE::Sense sense,
417 GeMatrix<MA> &A)
418 {
419 typedef typename GeMatrix<MA>::ElementType T;
420
421 const char BALANC = 'N';
422 const char JOBVL = (computeVL) ? 'V' : 'N';
423 const char JOBVR = (computeVR) ? 'V' : 'N';
424 const char SENSE = getF77LapackChar(sense);
425 const INTEGER N = A.numRows();
426 const INTEGER LDA = A.leadingDimension();
427 const INTEGER LDVL = N;
428 const INTEGER LDVR = N;
429 INTEGER IDUMMY;
430 T DUMMY;
431 const INTEGER LWORK = -1;
432 T WORK;
433 INTEGER INFO;
434
435 if (IsSame<T,DOUBLE>::value) {
436 LAPACK_IMPL(dgeevx)(&BALANC,
437 &JOBVL,
438 &JOBVR,
439 &SENSE,
440 &N,
441 A.data(),
442 &LDA,
443 &DUMMY,
444 &DUMMY,
445 &DUMMY,
446 &LDVL,
447 &DUMMY,
448 &LDVR,
449 &IDUMMY,
450 &IDUMMY,
451 &DUMMY,
452 &DUMMY,
453 &DUMMY,
454 &DUMMY,
455 &WORK,
456 &LWORK,
457 &IDUMMY,
458 &INFO);
459 } else {
460 ASSERT(0);
461 }
462 ASSERT(INFO>=0);
463
464 return Pair<typename GeMatrix<MA>::IndexType>(WORK, WORK);
465 }
466
467 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
468 typename IndexType, typename VSCALE, typename ABNORM,
469 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
470 typename GeMatrix<MA>::IndexType
471 evx_native(BALANCE::Balance balance,
472 bool computeVL,
473 bool computeVR,
474 SENSE::Sense sense,
475 GeMatrix<MA> &A,
476 DenseVector<VWR> &wr,
477 DenseVector<VWI> &wi,
478 GeMatrix<MVL> &VL,
479 GeMatrix<MVR> &VR,
480 IndexType &iLo,
481 IndexType &iHi,
482 DenseVector<VSCALE> &scale,
483 ABNORM &abNorm,
484 DenseVector<RCONDE> &rCondE,
485 DenseVector<RCONDV> &rCondV,
486 DenseVector<VWORK> &work,
487 DenseVector<VIWORK> &iWork)
488 {
489 typedef typename GeMatrix<MA>::ElementType T;
490
491 const char BALANC = char(balance);
492 const char JOBVL = (computeVL) ? 'V' : 'N';
493 const char JOBVR = (computeVR) ? 'V' : 'N';
494 const char SENSE = char(sense);
495 const INTEGER N = A.numRows();
496 const INTEGER LDA = A.leadingDimension();
497 const INTEGER LDVL = VL.leadingDimension();
498 const INTEGER LDVR = VR.leadingDimension();
499 INTEGER ILO = iLo;
500 INTEGER IHI = iHi;
501 T _ABNRM = abNorm;
502 const INTEGER LWORK = work.length();
503 INTEGER INFO;
504
505 DenseVector<Array<INTEGER> > _iWork(iWork.length());
506 for (INTEGER i=1; i<=_iWork.length(); ++i) {
507 _iWork(i) = iWork(i);
508 }
509
510 ASSERT(BALANC==char(balance));
511
512 if (IsSame<T,DOUBLE>::value) {
513 LAPACK_IMPL(dgeevx)(&BALANC,
514 &JOBVL,
515 &JOBVR,
516 &SENSE,
517 &N,
518 A.data(),
519 &LDA,
520 wr.data(),
521 wi.data(),
522 VL.data(),
523 &LDVL,
524 VR.data(),
525 &LDVR,
526 &ILO,
527 &IHI,
528 scale.data(),
529 &_ABNRM,
530 rCondE.data(),
531 rCondV.data(),
532 work.data(),
533 &LWORK,
534 _iWork.data(),
535 &INFO);
536 } else {
537 ASSERT(0);
538 }
539 ASSERT(INFO>=0);
540
541 for (INTEGER i=1; i<=_iWork.length(); ++i) {
542 iWork(i) = _iWork(i);
543 }
544 iLo = ILO;
545 iHi = IHI;
546 abNorm = _ABNRM;
547 return INFO;
548 }
549
550 #endif // CHECK_CXXLAPACK
551
552 //== public interface ==========================================================
553 template <typename MA>
554 Pair<typename GeMatrix<MA>::IndexType>
555 evx_wsq(bool computeVL,
556 bool computeVR,
557 SENSE::Sense sense,
558 GeMatrix<MA> &A)
559 {
560 return evx_generic_wsq(computeVL, computeVR, sense, A);
561 }
562
563 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
564 typename IndexType, typename VSCALE, typename ABNORM,
565 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
566 typename GeMatrix<MA>::IndexType
567 evx(BALANCE::Balance balance,
568 bool computeVL,
569 bool computeVR,
570 SENSE::Sense sense,
571 GeMatrix<MA> &A,
572 DenseVector<VWR> &wr,
573 DenseVector<VWI> &wi,
574 GeMatrix<MVL> &VL,
575 GeMatrix<MVR> &VR,
576 IndexType &iLo,
577 IndexType &iHi,
578 DenseVector<VSCALE> &scale,
579 ABNORM &abNorm,
580 DenseVector<RCONDE> &rCondE,
581 DenseVector<RCONDV> &rCondV,
582 DenseVector<VWORK> &work,
583 DenseVector<VIWORK> &iWork)
584 {
585 LAPACK_DEBUG_OUT("evx");
586
587 //
588 // Test the input parameters
589 //
590 # ifndef NDEBUG
591 ASSERT(A.numRows()==A.numCols());
592 ASSERT(A.firstRow()==1);
593 ASSERT(A.firstCol()==1);
594 ASSERT(work.firstIndex()==1);
595
596 typename GeMatrix<MA>::IndexType n = A.numRows();
597
598 ASSERT(wr.firstIndex()==1);
599 ASSERT(wr.length()==n);
600
601 ASSERT(wi.firstIndex()==1);
602 ASSERT(wi.length()==n);
603
604 if (computeVL) {
605 ASSERT(VL.numRows()==VL.numCols());
606 ASSERT(VL.numRows()==n);
607 ASSERT(VL.firstRow()==1);
608 ASSERT(VL.firstCol()==1);
609 }
610
611 if (computeVR) {
612 ASSERT(VR.numRows()==VR.numCols());
613 ASSERT(VR.numRows()==n);
614 ASSERT(VR.firstRow()==1);
615 ASSERT(VR.firstCol()==1);
616 }
617
618 ASSERT(scale.firstIndex()==1);
619 ASSERT(scale.length()==n);
620
621 ASSERT(rCondE.firstIndex()==1);
622 ASSERT(rCondE.length()==n);
623
624 ASSERT(rCondV.firstIndex()==1);
625 ASSERT(rCondV.length()==n);
626
627 ASSERT(iWork.length()==2*(n-1));
628 # endif
629
630 //
631 // Make copies of output arguments
632 //
633 # ifdef CHECK_CXXLAPACK
634 typename GeMatrix<MA>::NoView A_org = A;
635 typename DenseVector<VWR>::NoView wr_org = wr;
636 typename DenseVector<VWI>::NoView wi_org = wi;
637 typename GeMatrix<MVL>::NoView VL_org = VL;
638 typename GeMatrix<MVR>::NoView VR_org = VR;
639 IndexType iLo_org = iLo;
640 IndexType iHi_org = iHi;
641 typename DenseVector<VSCALE>::NoView scale_org = scale;
642 ABNORM abNorm_org = abNorm;
643 typename DenseVector<RCONDE>::NoView rCondE_org = rCondE;
644 typename DenseVector<RCONDV>::NoView rCondV_org = rCondV;
645 typename DenseVector<VWORK>::NoView work_org = work;
646 typename DenseVector<VIWORK>::NoView iWork_org = iWork;
647 # endif
648
649 //
650 // Call implementation
651 //
652 IndexType info = evx_generic(balance, computeVL, computeVR, sense,
653 A, wr, wi, VL, VR, iLo, iHi, scale, abNorm,
654 rCondE, rCondV, work, iWork);
655 # ifdef CHECK_CXXLAPACK
656 //
657 // Make copies of results computed by the generic implementation
658 //
659 typename GeMatrix<MA>::NoView A_generic = A;
660 typename DenseVector<VWR>::NoView wr_generic = wr;
661 typename DenseVector<VWI>::NoView wi_generic = wi;
662 typename GeMatrix<MVL>::NoView VL_generic = VL;
663 typename GeMatrix<MVR>::NoView VR_generic = VR;
664 IndexType iLo_generic = iLo;
665 IndexType iHi_generic = iHi;
666 typename DenseVector<VSCALE>::NoView scale_generic = scale;
667 ABNORM abNorm_generic = abNorm;
668 typename DenseVector<RCONDE>::NoView rCondE_generic = rCondE;
669 typename DenseVector<RCONDV>::NoView rCondV_generic = rCondV;
670 typename DenseVector<VWORK>::NoView work_generic = work;
671 typename DenseVector<VIWORK>::NoView iWork_generic = iWork;
672 //
673 // restore output arguments
674 //
675 A = A_org;
676 wr = wr_org;
677 wi = wi_org;
678 VL = VL_org;
679 VR = VR_org;
680 iLo = iLo_org;
681 iHi = iHi_org;
682 scale = scale_org;
683 abNorm = abNorm_org;
684 rCondE = rCondE_org;
685 rCondV = rCondV_org;
686 work = work_org;
687 iWork = iWork_org;
688 //
689 // Compare generic results with results from the native implementation
690 //
691 IndexType _info = evx_native(balance, computeVL, computeVR, sense,
692 A, wr, wi, VL, VR, iLo, iHi, scale, abNorm,
693 rCondE, rCondV, work, iWork);
694
695 bool failed = false;
696 if (! isIdentical(A_generic, A, "A_generic", "A")) {
697 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
698 std::cerr << "F77LAPACK: A = " << A << std::endl;
699 failed = true;
700 }
701
702 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
703 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
704 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
705 failed = true;
706 }
707
708 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
709 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
710 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
711 failed = true;
712 }
713
714 if (! isIdentical(VL_generic, VL, "VL_generic", "VL")) {
715 std::cerr << "CXXLAPACK: VL_generic = " << VL_generic << std::endl;
716 std::cerr << "F77LAPACK: VL = " << VL << std::endl;
717 failed = true;
718 }
719
720 if (! isIdentical(VR_generic, VR, "VR_generic", "VR")) {
721 std::cerr << "CXXLAPACK: VR_generic = " << VR_generic << std::endl;
722 std::cerr << "F77LAPACK: VR = " << VR << std::endl;
723 failed = true;
724 }
725
726 if (! isIdentical(iLo_generic, iLo, "iLo_generic", "iLo")) {
727 std::cerr << "CXXLAPACK: iLo_generic = " << iLo_generic << std::endl;
728 std::cerr << "F77LAPACK: iLo = " << iLo << std::endl;
729 failed = true;
730 }
731
732 if (! isIdentical(iHi_generic, iHi, "iHi_generic", "iHi")) {
733 std::cerr << "CXXLAPACK: iHi_generic = " << iHi_generic << std::endl;
734 std::cerr << "F77LAPACK: iHi = " << iHi << std::endl;
735 failed = true;
736 }
737
738 if (! isIdentical(scale_generic, scale, "scale_generic", "scale")) {
739 std::cerr << "CXXLAPACK: scale_generic = "
740 << scale_generic << std::endl;
741 std::cerr << "F77LAPACK: scale = " << scale << std::endl;
742 failed = true;
743 }
744
745 if (! isIdentical(abNorm_generic, abNorm, "abNorm_generic", "abNorm")) {
746 std::cerr << "CXXLAPACK: abNorm_generic = "
747 << abNorm_generic << std::endl;
748 std::cerr << "F77LAPACK: abNorm = " << abNorm << std::endl;
749 failed = true;
750 }
751
752 if (! isIdentical(rCondE_generic, rCondE, "rCondE_generic", "rCondE")) {
753 std::cerr << "CXXLAPACK: rCondE_generic = "
754 << rCondE_generic << std::endl;
755 std::cerr << "F77LAPACK: rCondE = " << rCondE << std::endl;
756 failed = true;
757 }
758
759 if (! isIdentical(rCondV_generic, rCondV, "rCondV_generic", "rCondV")) {
760 std::cerr << "CXXLAPACK: rCondV_generic = "
761 << rCondV_generic << std::endl;
762 std::cerr << "F77LAPACK: rCondV = " << rCondV << std::endl;
763 failed = true;
764 }
765
766 if (! isIdentical(work_generic, work, "work_generic", "work")) {
767 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
768 std::cerr << "F77LAPACK: work = " << work << std::endl;
769 failed = true;
770 }
771
772 if (! isIdentical(iWork_generic, iWork, "iWork_generic", "iWork")) {
773 std::cerr << "CXXLAPACK: iWork_generic = "
774 << iWork_generic << std::endl;
775 std::cerr << "F77LAPACK: iWork = " << iWork << std::endl;
776 failed = true;
777 }
778
779 if (! isIdentical(info, _info, " info", "_info")) {
780 std::cerr << "CXXLAPACK: info = " << info << std::endl;
781 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
782 failed = true;
783 }
784
785 if (failed) {
786 ASSERT(0);
787 } else {
788 // std::cerr << "passed: evx.tcc" << std::endl;
789 }
790 # endif
791
792 return info;
793 }
794
795 //-- forwarding ----------------------------------------------------------------
796
797 } } // namespace lapack, flens
798
799 #endif // FLENS_LAPACK_EIG_EVX_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 DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
36 $ LDVR, WORK, LWORK, INFO )
37 *
38 * -- LAPACK driver routine (version 3.3.1) --
39 * -- LAPACK is a software package provided by Univ. of Tennessee, --
40 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
41 * -- April 2011 --
42 */
43
44 #ifndef FLENS_LAPACK_EIG_EVX_TCC
45 #define FLENS_LAPACK_EIG_EVX_TCC 1
46
47 #include <flens/blas/blas.h>
48 #include <flens/lapack/lapack.h>
49
50 namespace flens { namespace lapack {
51
52 //== generic lapack implementation =============================================
53 template <typename MA>
54 Pair<typename GeMatrix<MA>::IndexType>
55 evx_generic_wsq(bool computeVL,
56 bool computeVR,
57 SENSE::Sense sense,
58 GeMatrix<MA> &A)
59 {
60 using std::max;
61
62 typedef typename GeMatrix<MA>::ElementType T;
63 typedef typename GeMatrix<MA>::IndexType IndexType;
64
65 const IndexType n = A.numRows();
66 const bool wantSN = (sense==SENSE::None);
67 const bool wantSE = (sense==SENSE::EigenvaluesOnly);
68
69 IndexType minWork, maxWork;
70
71 if (n==0) {
72 minWork = 1;
73 maxWork = 1;
74 } else {
75 maxWork = n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
76
77 IndexType hseqrWork;
78 if (computeVL) {
79 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
80 IndexType(1), n, A);
81 } else if (computeVL) {
82 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
83 IndexType(1), n, A);
84 } else {
85 if (wantSN) {
86 hseqrWork = hseqr_wsq(HSEQR::Eigenvalues, HSEQR::No,
87 IndexType(1), n, A);
88 } else {
89 hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::No,
90 IndexType(1), n, A);
91 }
92 }
93
94 if ((!computeVL) && (!computeVR)) {
95 minWork = 2*n;
96 if (!wantSN) {
97 minWork = max(minWork, n*n+6*n);
98 }
99 maxWork = max(maxWork, hseqrWork);
100 if (!wantSN) {
101 maxWork = max(maxWork, n*n+6*n);
102 }
103 } else {
104 minWork = 3*n;
105 if ((!wantSN) && (!wantSE)) {
106 minWork = max(minWork, n*n + 6*n);
107 }
108 maxWork = max(maxWork, hseqrWork);
109 maxWork = max(maxWork, n+(n-1*ilaenv<T>(1, "ORGHR", "", n, 1, n)));
110 if ((!wantSN) && (!wantSE)) {
111 maxWork = max(maxWork, n*n + 6*n);
112 }
113 maxWork = max(maxWork, 3*n);
114 }
115 maxWork = max(maxWork, minWork);
116 }
117
118 return Pair<typename GeMatrix<MA>::IndexType>(minWork, maxWork);
119 }
120
121 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
122 typename IndexType, typename VSCALE, typename ABNORM,
123 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
124 typename GeMatrix<MA>::IndexType
125 evx_generic(BALANCE::Balance balance,
126 bool computeVL,
127 bool computeVR,
128 SENSE::Sense sense,
129 GeMatrix<MA> &A,
130 DenseVector<VWR> &wr,
131 DenseVector<VWI> &wi,
132 GeMatrix<MVL> &VL,
133 GeMatrix<MVR> &VR,
134 IndexType &iLo,
135 IndexType &iHi,
136 DenseVector<VSCALE> &scale,
137 ABNORM &ABNorm,
138 DenseVector<RCONDE> &rCondE,
139 DenseVector<RCONDV> &rCondV,
140 DenseVector<VWORK> &work,
141 DenseVector<VIWORK> &iWork)
142 {
143 typedef typename GeMatrix<MA>::ElementType T;
144
145 const IndexType n = A.numRows();
146 const Underscore<IndexType> _;
147
148 const T Zero(0), One(1);
149 //
150 // .. Local Arrays ..
151 //
152 bool _selectData[1];
153 DenseVectorView<bool>
154 select = typename DenseVectorView<bool>::Engine(1, _selectData);
155 //
156 // Test the input arguments
157 //
158 IndexType info = 0;
159
160 const bool wantSN = (sense==SENSE::None);
161 const bool wantSV = (sense==SENSE::InvariantSubspaceOnly);
162 const bool wantSB = (sense==SENSE::Both);
163 //
164 // Compute workspace
165 // (Note: Comments in the code beginning "Workspace:" describe the
166 // minimal amount of workspace needed at that point in the code,
167 // as well as the preferred amount for good performance.
168 // NB refers to the optimal block size for the immediately
169 // following subroutine, as returned by ILAENV.
170 // HSWORK refers to the workspace preferred by DHSEQR, as
171 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
172 // the worst case.)
173 //
174 Pair<IndexType> wsQuery = evx_wsq(computeVL, computeVR, sense, A);
175 IndexType minWork = wsQuery.first;
176 IndexType maxWork = wsQuery.second;
177
178 if (work.length()!=0 && work.length()<minWork) {
179 ASSERT(0);
180 } else if (work.length()==0) {
181 work.resize(maxWork);
182 }
183 work(1) = maxWork;
184
185 //
186 // Quick return if possible
187 //
188 if (n==0) {
189 return info;
190 }
191 //
192 // Get machine constants
193 //
194 const T eps = lamch<T>(Precision);
195 T smallNum = lamch<T>(SafeMin);
196 T bigNum = One / smallNum;
197 labad(smallNum, bigNum);
198 smallNum = sqrt(smallNum) / eps;
199 bigNum = One / smallNum;
200 //
201 // Scale A if max element outside range [SMLNUM,BIGNUM]
202 //
203 IndexType iCond = 0;
204 const T ANorm = lan(MaximumNorm, A);
205 bool scaleA = false;
206 T cScale;
207 if (ANorm>Zero && ANorm<smallNum) {
208 scaleA = true;
209 cScale = smallNum;
210 } else if (ANorm>bigNum) {
211 scaleA = true;
212 cScale = bigNum;
213 }
214 if (scaleA) {
215 lascl(LASCL::FullMatrix, 0, 0, ANorm, cScale, A);
216 }
217 //
218 // Balance the matrix and compute ABNRM
219 //
220 bal(balance, A, iLo, iHi, scale);
221 ABNorm = lan(OneNorm, A);
222 if (scaleA) {
223 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, ABNorm);
224 }
225 //
226 // Reduce to upper Hessenberg form
227 // (Workspace: need 2*N, prefer N+N*NB)
228 //
229 IndexType iTau = 1;
230 IndexType iWrk = iTau + n;
231 IndexType lWork = work.length();
232
233 auto tau = work(_(iTau, iTau+n-2));
234 auto hrdWork = work(_(iWrk, lWork));
235
236 hrd(iLo, iHi, A, tau, hrdWork);
237
238 if (computeVL) {
239 //
240 // Want left eigenvectors
241 // Copy Householder vectors to VL
242 //
243 VL.lower() = A.lower();
244 //
245 // Generate orthogonal matrix in VL
246 // (Workspace: need 2*N-1, prefer N+(N-1)*NB)
247 //
248 orghr(iLo, iHi, VL, tau, hrdWork);
249 //
250 // Perform QR iteration, accumulating Schur vectors in VL
251 // (Workspace: need 1, prefer HSWORK (see comments) )
252 //
253 iWrk = iTau;
254 auto hseqrWork = work(_(iWrk, lWork));
255 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
256 wr, wi, VL, hseqrWork);
257
258 if (computeVR) {
259 //
260 // Want left and right eigenvectors
261 // Copy Schur vectors to VR
262 //
263 VR = VL;
264 }
265
266 } else if (computeVR) {
267 //
268 // Want right eigenvectors
269 // Copy Householder vectors to VR
270 //
271 VR.lower() = A.lower();
272 //
273 // Generate orthogonal matrix in VR
274 // (Workspace: need 2*N-1, prefer N+(N-1)*NB)
275 //
276 orghr(iLo, iHi, VR, tau, hrdWork);
277 //
278 // Perform QR iteration, accumulating Schur vectors in VR
279 // (Workspace: need 1, prefer HSWORK (see comments) )
280 //
281 iWrk = iTau;
282 auto hseqrWork = work(_(iWrk, lWork));
283 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
284 wr, wi, VR, hseqrWork);
285
286 } else {
287 //
288 // Compute eigenvalues only
289 // If condition numbers desired, compute Schur form
290 //
291 HSEQR::Job job = (wantSN) ? HSEQR::Eigenvalues
292 : HSEQR::Schur;
293 //
294 // (Workspace: need 1, prefer HSWORK (see comments) )
295 //
296 iWrk = iTau;
297 auto hseqrWork = work(_(iWrk, lWork));
298 info = hseqr(job, HSEQR::No, iLo, iHi, A, wr, wi, VR, hseqrWork);
299 }
300 //
301 // If INFO > 0 from DHSEQR, then quit
302 //
303 if (info==0) {
304
305 if (computeVL || computeVR) {
306 //
307 // Compute left and/or right eigenvectors
308 // (Workspace: need 3*N)
309 //
310 IndexType nOut;
311 auto trevcWork = work(_(iWrk, lWork));
312
313 trevc(computeVL, computeVR, TREVC::Backtransform, select,
314 A, VL, VR, n, nOut, trevcWork);
315 }
316 //
317 // Compute condition numbers if desired
318 // (Workspace: need N*N+6*N unless SENSE = 'E')
319 //
320 if (!wantSN) {
321 IndexType nOut;
322
323 IndexType _n = (sense!=SENSE::EigenvaluesOnly) ? n : 0;
324 GeMatrixView<T> Work(_n, n+6, work(_(iWrk, lWork)), n);
325
326 trsna(TRSNA::Job(sense), TRSNA::All, select, A, VL, VR,
327 rCondE, rCondV, n, nOut, Work, iWork);
328 }
329
330 if (computeVL) {
331 //
332 // Undo balancing of left eigenvectors
333 //
334 bak(balance, Left, iLo, iHi, scale, VL);
335 //
336 // Normalize left eigenvectors and make largest component real
337 //
338 for (IndexType i=1; i<=n; ++i) {
339 if (wi(i)==Zero) {
340 VL(_,i) *= One / blas::nrm2(VL(_,i));
341 } else if (wi(i)>Zero) {
342 const T scl = One / lapy2(blas::nrm2(VL(_,i)),
343 blas::nrm2(VL(_,i+1)));
344 VL(_,i) *= scl;
345 VL(_,i+1) *= scl;
346 for (IndexType k=1; k<=n; ++k) {
347 work(k) = pow(VL(k,i),2) + pow(VL(k,i+1),2);
348 }
349 IndexType k = blas::iamax(work(_(1,n)));
350 T cs, sn, r;
351 lartg(VL(k,i), VL(k,i+1), cs, sn, r);
352 blas::rot(VL(_,i), VL(_,i+1), cs, sn);
353 VL(k,i+1) = Zero;
354 }
355 }
356 }
357
358 if (computeVR) {
359 //
360 // Undo balancing of right eigenvectors
361 //
362 bak(BALANCE::Both, Right, iLo, iHi, scale, VR);
363 //
364 // Normalize right eigenvectors and make largest component real
365 //
366 for (IndexType i=1; i<=n; ++i) {
367 if (wi(i)==Zero) {
368 VR(_,i) *= One / blas::nrm2(VR(_,i));
369 } else if (wi(i)>Zero) {
370 const T scl = One / lapy2(blas::nrm2(VR(_,i)),
371 blas::nrm2(VR(_,i+1)));
372 VR(_,i) *= scl;
373 VR(_,i+1) *= scl;
374 for (IndexType k=1; k<=n; ++k) {
375 work(k) = pow(VR(k,i), 2) + pow(VR(k,i+1), 2);
376 }
377 IndexType k = blas::iamax(work(_(1,n)));
378 T cs, sn, r;
379 lartg(VR(k,i), VR(k,i+1), cs, sn, r);
380 blas::rot(VR(_,i), VR(_,i+1), cs, sn);
381 VR(k,i+1) = Zero;
382 }
383 }
384 }
385 }
386
387 //
388 // Undo scaling if necessary
389 //
390 if (scaleA) {
391 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(info+1,n)));
392 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(info+1,n)));
393
394 if (info==0) {
395 if ((wantSV || wantSB) && iCond==0) {
396 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, rCondV);
397 }
398 } else {
399 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(1,iLo-1)));
400 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(1,iLo-1)));
401 }
402 }
403
404 work(1) = maxWork;
405 return info;
406 }
407
408 //== interface for native lapack ===============================================
409
410 #ifdef CHECK_CXXLAPACK
411
412 template <typename MA>
413 Pair<typename GeMatrix<MA>::IndexType>
414 evx_native_wsq(bool computeVL,
415 bool computeVR,
416 SENSE::Sense sense,
417 GeMatrix<MA> &A)
418 {
419 typedef typename GeMatrix<MA>::ElementType T;
420
421 const char BALANC = 'N';
422 const char JOBVL = (computeVL) ? 'V' : 'N';
423 const char JOBVR = (computeVR) ? 'V' : 'N';
424 const char SENSE = getF77LapackChar(sense);
425 const INTEGER N = A.numRows();
426 const INTEGER LDA = A.leadingDimension();
427 const INTEGER LDVL = N;
428 const INTEGER LDVR = N;
429 INTEGER IDUMMY;
430 T DUMMY;
431 const INTEGER LWORK = -1;
432 T WORK;
433 INTEGER INFO;
434
435 if (IsSame<T,DOUBLE>::value) {
436 LAPACK_IMPL(dgeevx)(&BALANC,
437 &JOBVL,
438 &JOBVR,
439 &SENSE,
440 &N,
441 A.data(),
442 &LDA,
443 &DUMMY,
444 &DUMMY,
445 &DUMMY,
446 &LDVL,
447 &DUMMY,
448 &LDVR,
449 &IDUMMY,
450 &IDUMMY,
451 &DUMMY,
452 &DUMMY,
453 &DUMMY,
454 &DUMMY,
455 &WORK,
456 &LWORK,
457 &IDUMMY,
458 &INFO);
459 } else {
460 ASSERT(0);
461 }
462 ASSERT(INFO>=0);
463
464 return Pair<typename GeMatrix<MA>::IndexType>(WORK, WORK);
465 }
466
467 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
468 typename IndexType, typename VSCALE, typename ABNORM,
469 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
470 typename GeMatrix<MA>::IndexType
471 evx_native(BALANCE::Balance balance,
472 bool computeVL,
473 bool computeVR,
474 SENSE::Sense sense,
475 GeMatrix<MA> &A,
476 DenseVector<VWR> &wr,
477 DenseVector<VWI> &wi,
478 GeMatrix<MVL> &VL,
479 GeMatrix<MVR> &VR,
480 IndexType &iLo,
481 IndexType &iHi,
482 DenseVector<VSCALE> &scale,
483 ABNORM &abNorm,
484 DenseVector<RCONDE> &rCondE,
485 DenseVector<RCONDV> &rCondV,
486 DenseVector<VWORK> &work,
487 DenseVector<VIWORK> &iWork)
488 {
489 typedef typename GeMatrix<MA>::ElementType T;
490
491 const char BALANC = char(balance);
492 const char JOBVL = (computeVL) ? 'V' : 'N';
493 const char JOBVR = (computeVR) ? 'V' : 'N';
494 const char SENSE = char(sense);
495 const INTEGER N = A.numRows();
496 const INTEGER LDA = A.leadingDimension();
497 const INTEGER LDVL = VL.leadingDimension();
498 const INTEGER LDVR = VR.leadingDimension();
499 INTEGER ILO = iLo;
500 INTEGER IHI = iHi;
501 T _ABNRM = abNorm;
502 const INTEGER LWORK = work.length();
503 INTEGER INFO;
504
505 DenseVector<Array<INTEGER> > _iWork(iWork.length());
506 for (INTEGER i=1; i<=_iWork.length(); ++i) {
507 _iWork(i) = iWork(i);
508 }
509
510 ASSERT(BALANC==char(balance));
511
512 if (IsSame<T,DOUBLE>::value) {
513 LAPACK_IMPL(dgeevx)(&BALANC,
514 &JOBVL,
515 &JOBVR,
516 &SENSE,
517 &N,
518 A.data(),
519 &LDA,
520 wr.data(),
521 wi.data(),
522 VL.data(),
523 &LDVL,
524 VR.data(),
525 &LDVR,
526 &ILO,
527 &IHI,
528 scale.data(),
529 &_ABNRM,
530 rCondE.data(),
531 rCondV.data(),
532 work.data(),
533 &LWORK,
534 _iWork.data(),
535 &INFO);
536 } else {
537 ASSERT(0);
538 }
539 ASSERT(INFO>=0);
540
541 for (INTEGER i=1; i<=_iWork.length(); ++i) {
542 iWork(i) = _iWork(i);
543 }
544 iLo = ILO;
545 iHi = IHI;
546 abNorm = _ABNRM;
547 return INFO;
548 }
549
550 #endif // CHECK_CXXLAPACK
551
552 //== public interface ==========================================================
553 template <typename MA>
554 Pair<typename GeMatrix<MA>::IndexType>
555 evx_wsq(bool computeVL,
556 bool computeVR,
557 SENSE::Sense sense,
558 GeMatrix<MA> &A)
559 {
560 return evx_generic_wsq(computeVL, computeVR, sense, A);
561 }
562
563 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
564 typename IndexType, typename VSCALE, typename ABNORM,
565 typename RCONDE, typename RCONDV, typename VWORK, typename VIWORK>
566 typename GeMatrix<MA>::IndexType
567 evx(BALANCE::Balance balance,
568 bool computeVL,
569 bool computeVR,
570 SENSE::Sense sense,
571 GeMatrix<MA> &A,
572 DenseVector<VWR> &wr,
573 DenseVector<VWI> &wi,
574 GeMatrix<MVL> &VL,
575 GeMatrix<MVR> &VR,
576 IndexType &iLo,
577 IndexType &iHi,
578 DenseVector<VSCALE> &scale,
579 ABNORM &abNorm,
580 DenseVector<RCONDE> &rCondE,
581 DenseVector<RCONDV> &rCondV,
582 DenseVector<VWORK> &work,
583 DenseVector<VIWORK> &iWork)
584 {
585 LAPACK_DEBUG_OUT("evx");
586
587 //
588 // Test the input parameters
589 //
590 # ifndef NDEBUG
591 ASSERT(A.numRows()==A.numCols());
592 ASSERT(A.firstRow()==1);
593 ASSERT(A.firstCol()==1);
594 ASSERT(work.firstIndex()==1);
595
596 typename GeMatrix<MA>::IndexType n = A.numRows();
597
598 ASSERT(wr.firstIndex()==1);
599 ASSERT(wr.length()==n);
600
601 ASSERT(wi.firstIndex()==1);
602 ASSERT(wi.length()==n);
603
604 if (computeVL) {
605 ASSERT(VL.numRows()==VL.numCols());
606 ASSERT(VL.numRows()==n);
607 ASSERT(VL.firstRow()==1);
608 ASSERT(VL.firstCol()==1);
609 }
610
611 if (computeVR) {
612 ASSERT(VR.numRows()==VR.numCols());
613 ASSERT(VR.numRows()==n);
614 ASSERT(VR.firstRow()==1);
615 ASSERT(VR.firstCol()==1);
616 }
617
618 ASSERT(scale.firstIndex()==1);
619 ASSERT(scale.length()==n);
620
621 ASSERT(rCondE.firstIndex()==1);
622 ASSERT(rCondE.length()==n);
623
624 ASSERT(rCondV.firstIndex()==1);
625 ASSERT(rCondV.length()==n);
626
627 ASSERT(iWork.length()==2*(n-1));
628 # endif
629
630 //
631 // Make copies of output arguments
632 //
633 # ifdef CHECK_CXXLAPACK
634 typename GeMatrix<MA>::NoView A_org = A;
635 typename DenseVector<VWR>::NoView wr_org = wr;
636 typename DenseVector<VWI>::NoView wi_org = wi;
637 typename GeMatrix<MVL>::NoView VL_org = VL;
638 typename GeMatrix<MVR>::NoView VR_org = VR;
639 IndexType iLo_org = iLo;
640 IndexType iHi_org = iHi;
641 typename DenseVector<VSCALE>::NoView scale_org = scale;
642 ABNORM abNorm_org = abNorm;
643 typename DenseVector<RCONDE>::NoView rCondE_org = rCondE;
644 typename DenseVector<RCONDV>::NoView rCondV_org = rCondV;
645 typename DenseVector<VWORK>::NoView work_org = work;
646 typename DenseVector<VIWORK>::NoView iWork_org = iWork;
647 # endif
648
649 //
650 // Call implementation
651 //
652 IndexType info = evx_generic(balance, computeVL, computeVR, sense,
653 A, wr, wi, VL, VR, iLo, iHi, scale, abNorm,
654 rCondE, rCondV, work, iWork);
655 # ifdef CHECK_CXXLAPACK
656 //
657 // Make copies of results computed by the generic implementation
658 //
659 typename GeMatrix<MA>::NoView A_generic = A;
660 typename DenseVector<VWR>::NoView wr_generic = wr;
661 typename DenseVector<VWI>::NoView wi_generic = wi;
662 typename GeMatrix<MVL>::NoView VL_generic = VL;
663 typename GeMatrix<MVR>::NoView VR_generic = VR;
664 IndexType iLo_generic = iLo;
665 IndexType iHi_generic = iHi;
666 typename DenseVector<VSCALE>::NoView scale_generic = scale;
667 ABNORM abNorm_generic = abNorm;
668 typename DenseVector<RCONDE>::NoView rCondE_generic = rCondE;
669 typename DenseVector<RCONDV>::NoView rCondV_generic = rCondV;
670 typename DenseVector<VWORK>::NoView work_generic = work;
671 typename DenseVector<VIWORK>::NoView iWork_generic = iWork;
672 //
673 // restore output arguments
674 //
675 A = A_org;
676 wr = wr_org;
677 wi = wi_org;
678 VL = VL_org;
679 VR = VR_org;
680 iLo = iLo_org;
681 iHi = iHi_org;
682 scale = scale_org;
683 abNorm = abNorm_org;
684 rCondE = rCondE_org;
685 rCondV = rCondV_org;
686 work = work_org;
687 iWork = iWork_org;
688 //
689 // Compare generic results with results from the native implementation
690 //
691 IndexType _info = evx_native(balance, computeVL, computeVR, sense,
692 A, wr, wi, VL, VR, iLo, iHi, scale, abNorm,
693 rCondE, rCondV, work, iWork);
694
695 bool failed = false;
696 if (! isIdentical(A_generic, A, "A_generic", "A")) {
697 std::cerr << "CXXLAPACK: A_generic = " << A_generic << std::endl;
698 std::cerr << "F77LAPACK: A = " << A << std::endl;
699 failed = true;
700 }
701
702 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
703 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
704 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
705 failed = true;
706 }
707
708 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
709 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
710 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
711 failed = true;
712 }
713
714 if (! isIdentical(VL_generic, VL, "VL_generic", "VL")) {
715 std::cerr << "CXXLAPACK: VL_generic = " << VL_generic << std::endl;
716 std::cerr << "F77LAPACK: VL = " << VL << std::endl;
717 failed = true;
718 }
719
720 if (! isIdentical(VR_generic, VR, "VR_generic", "VR")) {
721 std::cerr << "CXXLAPACK: VR_generic = " << VR_generic << std::endl;
722 std::cerr << "F77LAPACK: VR = " << VR << std::endl;
723 failed = true;
724 }
725
726 if (! isIdentical(iLo_generic, iLo, "iLo_generic", "iLo")) {
727 std::cerr << "CXXLAPACK: iLo_generic = " << iLo_generic << std::endl;
728 std::cerr << "F77LAPACK: iLo = " << iLo << std::endl;
729 failed = true;
730 }
731
732 if (! isIdentical(iHi_generic, iHi, "iHi_generic", "iHi")) {
733 std::cerr << "CXXLAPACK: iHi_generic = " << iHi_generic << std::endl;
734 std::cerr << "F77LAPACK: iHi = " << iHi << std::endl;
735 failed = true;
736 }
737
738 if (! isIdentical(scale_generic, scale, "scale_generic", "scale")) {
739 std::cerr << "CXXLAPACK: scale_generic = "
740 << scale_generic << std::endl;
741 std::cerr << "F77LAPACK: scale = " << scale << std::endl;
742 failed = true;
743 }
744
745 if (! isIdentical(abNorm_generic, abNorm, "abNorm_generic", "abNorm")) {
746 std::cerr << "CXXLAPACK: abNorm_generic = "
747 << abNorm_generic << std::endl;
748 std::cerr << "F77LAPACK: abNorm = " << abNorm << std::endl;
749 failed = true;
750 }
751
752 if (! isIdentical(rCondE_generic, rCondE, "rCondE_generic", "rCondE")) {
753 std::cerr << "CXXLAPACK: rCondE_generic = "
754 << rCondE_generic << std::endl;
755 std::cerr << "F77LAPACK: rCondE = " << rCondE << std::endl;
756 failed = true;
757 }
758
759 if (! isIdentical(rCondV_generic, rCondV, "rCondV_generic", "rCondV")) {
760 std::cerr << "CXXLAPACK: rCondV_generic = "
761 << rCondV_generic << std::endl;
762 std::cerr << "F77LAPACK: rCondV = " << rCondV << std::endl;
763 failed = true;
764 }
765
766 if (! isIdentical(work_generic, work, "work_generic", "work")) {
767 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
768 std::cerr << "F77LAPACK: work = " << work << std::endl;
769 failed = true;
770 }
771
772 if (! isIdentical(iWork_generic, iWork, "iWork_generic", "iWork")) {
773 std::cerr << "CXXLAPACK: iWork_generic = "
774 << iWork_generic << std::endl;
775 std::cerr << "F77LAPACK: iWork = " << iWork << std::endl;
776 failed = true;
777 }
778
779 if (! isIdentical(info, _info, " info", "_info")) {
780 std::cerr << "CXXLAPACK: info = " << info << std::endl;
781 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
782 failed = true;
783 }
784
785 if (failed) {
786 ASSERT(0);
787 } else {
788 // std::cerr << "passed: evx.tcc" << std::endl;
789 }
790 # endif
791
792 return info;
793 }
794
795 //-- forwarding ----------------------------------------------------------------
796
797 } } // namespace lapack, flens
798
799 #endif // FLENS_LAPACK_EIG_EVX_TCC