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_EV_TCC
45 #define FLENS_LAPACK_EIG_EV_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
54 //-- ev: workspace query
55 template <typename MA>
56 Pair<typename GeMatrix<MA>::IndexType>
57 ev_generic_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
58 {
59 using std::max;
60
61 typedef typename GeMatrix<MA>::ElementType T;
62 typedef typename GeMatrix<MA>::IndexType IndexType;
63
64 const IndexType n = A.numRows();
65
66 IndexType minWork, maxWork;
67
68 if (n==0) {
69 minWork = 1;
70 maxWork = 1;
71 } else {
72 maxWork = 2*n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
73 if (computeVL) {
74 minWork = 4*n;
75 maxWork = max(maxWork,
76 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
77 IndexType hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
78 IndexType(1), n, A);
79 maxWork = max(max(maxWork, n+1), n+hseqrWork);
80 maxWork = max(maxWork, 4*n);
81 } else if (computeVR) {
82 minWork = 4*n;
83 maxWork = max(maxWork,
84 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
85 IndexType hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
86 IndexType(1), n, A);
87 maxWork = max(max(maxWork, n+1), n+hseqrWork);
88 maxWork = max(maxWork, 4*n);
89 } else {
90 minWork = 3*n;
91 IndexType hseqrWork = hseqr_wsq(HSEQR::Eigenvalues, HSEQR::No,
92 IndexType(1), n, A);
93 maxWork = max(max(maxWork, n+1), n+hseqrWork);
94 }
95 maxWork = max(maxWork, minWork);
96 }
97 return Pair<IndexType>(minWork, maxWork);
98 }
99
100 //-- ev: computation
101 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
102 typename VWORK>
103 typename GeMatrix<MA>::IndexType
104 ev_generic(bool computeVL, bool computeVR,
105 GeMatrix<MA> &A,
106 DenseVector<VWR> &wr, DenseVector<VWI> &wi,
107 GeMatrix<MVL> &VL, GeMatrix<MVR> &VR,
108 DenseVector<VWORK> &work)
109 {
110 using std::pow;
111 using std::sqrt;
112
113 typedef typename GeMatrix<MA>::ElementType T;
114 typedef typename GeMatrix<MA>::IndexType IndexType;
115
116 const Underscore<IndexType> _;
117 const IndexType n = A.numRows();
118 const T Zero(0), One(1);
119 //
120 // Compute workspace
121 // (Note: Comments in the code beginning "Workspace:" describe the
122 // minimal amount of workspace needed at that point in the code,
123 // as well as the preferred amount for good performance.
124 // NB refers to the optimal block size for the immediately
125 // following subroutine, as returned by ILAENV.
126 // HSWORK refers to the workspace preferred by DHSEQR, as
127 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
128 // the worst case.)
129 //
130 Pair<IndexType> wsQuery = ev_wsq(computeVL, computeVR, A);
131 IndexType minWork = wsQuery.first;
132 IndexType maxWork = wsQuery.second;
133
134 if (work.length()!=0 && work.length()<minWork) {
135 ASSERT(0);
136 } else if (work.length()==0) {
137 work.resize(maxWork);
138 }
139
140 work(1) = maxWork;
141
142 //
143 // Quick return if possible
144 //
145 if (n==0) {
146 return 0;
147 }
148 //
149 // Get machine constants
150 //
151 const T Eps = lamch<T>(Precision);
152 T smallNum = lamch<T>(SafeMin);
153 T bigNum = One / smallNum;
154 labad(smallNum, bigNum);
155 smallNum = sqrt(smallNum) / Eps;
156 bigNum = One / smallNum;
157 //
158 // Scale A if max element outside range [SMLNUM,BIGNUM]
159 //
160 const T ANorm = lan(MaximumNorm, A);
161 bool scaleA = false;
162 T cScale;
163 if (ANorm>Zero && ANorm<smallNum) {
164 scaleA = true;
165 cScale = smallNum;
166 } else if (ANorm>bigNum) {
167 scaleA = true;
168 cScale = bigNum;
169 }
170 if (scaleA) {
171 lascl(LASCL::FullMatrix, 0, 0, ANorm, cScale, A);
172 }
173 //
174 // Balance the matrix
175 // (Workspace: need N)
176 //
177 IndexType iBal = 1;
178 IndexType iLo, iHi;
179
180 auto balWork = work(_(iBal,iBal+n-1));
181 bal(BALANCE::Both, A, iLo, iHi, balWork);
182
183 //
184 // Reduce to upper Hessenberg form
185 // (Workspace: need 3*N, prefer 2*N+N*NB)
186 //
187 IndexType iTau = iBal + n;
188 IndexType iWork = iTau + n;
189 IndexType lWork = work.length();
190
191 auto tau = work(_(iTau, iTau+n-2));
192 auto hrdWork = work(_(iWork, lWork));
193
194 hrd(iLo, iHi, A, tau, hrdWork);
195
196 IndexType info = 0;
197 if (computeVL) {
198 //
199 // Want left eigenvectors
200 // Copy Householder vectors to VL
201 //
202 VL.lower() = A.lower();
203 //
204 // Generate orthogonal matrix in VL
205 // (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
206 //
207 orghr(iLo, iHi, VL, tau, hrdWork);
208 //
209 // Perform QR iteration, accumulating Schur vectors in VL
210 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
211 //
212 iWork = iTau;
213 auto hseqrWork = work(_(iWork, lWork));
214 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
215 wr, wi, VL, hseqrWork);
216
217 if (computeVR) {
218 //
219 // Want left and right eigenvectors
220 // Copy Schur vectors to VR
221 //
222 VR = VL;
223 }
224
225 } else if (computeVR) {
226 //
227 // Want right eigenvectors
228 // Copy Householder vectors to VR
229 //
230 VR.lower() = A.lower();
231
232 //
233 // Generate orthogonal matrix in VR
234 // (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
235 //
236 orghr(iLo, iHi, VR, tau, hrdWork);
237 //
238 // Perform QR iteration, accumulating Schur vectors in VR
239 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
240 //
241 iWork = iTau;
242 auto hseqrWork = work(_(iWork, lWork));
243 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
244 wr, wi, VR, hseqrWork);
245
246 } else {
247 //
248 // Compute eigenvalues only
249 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
250 //
251 iWork = iTau;
252 auto hseqrWork = work(_(iWork, lWork));
253 info = hseqr(HSEQR::Eigenvalues, HSEQR::No, iLo, iHi, A,
254 wr, wi, VR, hseqrWork);
255 }
256
257 //
258 // If INFO > 0 from DHSEQR, then quit
259 //
260 if (info==0) {
261
262 if (computeVL || computeVR) {
263 //
264 // Compute left and/or right eigenvectors
265 // (Workspace: need 4*N)
266 //
267 // TODO: shouldn't that be "Workspace: need 3*N" ?? instead of 4*N?
268 IndexType nOut;
269 auto trevcWork = work(_(iWork, lWork));
270
271 // TODO: find a way that we don't have to pass an empty "select"
272 DenseVector<Array<bool> > select;
273 trevc(computeVL, computeVR, TREVC::Backtransform, select,
274 A, VL, VR, n, nOut, trevcWork);
275 }
276
277 if (computeVL) {
278 //
279 // Undo balancing of left eigenvectors
280 // (Workspace: need N)
281 //
282 bak(BALANCE::Both, Left, iLo, iHi, balWork, VL);
283 //
284 // Normalize left eigenvectors and make largest component real
285 //
286 auto _work = work(_(iWork, iWork+n-1));
287 for (IndexType i=1; i<=n; ++i) {
288 if (wi(i)==Zero) {
289 VL(_,i) *= One / blas::nrm2(VL(_,i));
290 } else if (wi(i)>Zero) {
291 T scale = One / lapy2(blas::nrm2(VL(_,i )),
292 blas::nrm2(VL(_,i+1)));
293 VL(_,i) *= scale;
294 VL(_,i+1) *= scale;
295 for (IndexType k=1; k<=n; ++k) {
296 _work(k) = pow(VL(k,i), 2) + pow(VL(k,i+1), 2);
297 }
298 IndexType k = blas::iamax(_work);
299 T cs, sn, r;
300 lartg(VL(k,i), VL(k,i+1), cs, sn, r);
301 blas::rot(VL(_,i), VL(_,i+1), cs, sn);
302 VL(k,i+1) = Zero;
303 }
304 }
305 }
306
307 if (computeVR) {
308 //
309 // Undo balancing of right eigenvectors
310 // (Workspace: need N)
311 //
312 bak(BALANCE::Both, Right, iLo, iHi, balWork, VR);
313 //
314 // Normalize right eigenvectors and make largest component real
315 //
316 auto _work = work(_(iWork, iWork+n-1));
317 for (IndexType i=1; i<=n; ++i) {
318 if (wi(i)==Zero) {
319 VR(_,i) *= One / blas::nrm2(VR(_,i));
320 } else if (wi(i)>Zero) {
321 T scale = One / lapy2(blas::nrm2(VR(_,i )),
322 blas::nrm2(VR(_,i+1)));
323 VR(_,i) *= scale;
324 VR(_,i+1) *= scale;
325 for (IndexType k=1; k<=n; ++k) {
326 _work(k) = pow(VR(k,i), 2) + pow(VR(k,i+1), 2);
327 }
328 IndexType k = blas::iamax(_work);
329 T cs, sn, r;
330 lartg(VR(k,i), VR(k,i+1), cs, sn, r);
331 blas::rot(VR(_,i), VR(_,i+1), cs, sn);
332 VR(k,i+1) = Zero;
333 }
334 }
335 }
336 }
337 //
338 // Undo scaling if necessary
339 //
340 if (scaleA) {
341 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(info+1,n)));
342 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(info+1,n)));
343
344 if (info>0) {
345 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(1,iLo-1)));
346 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(1,iLo-1)));
347 }
348 }
349
350 work(1) = maxWork;
351 return info;
352 }
353
354 //== interface for native lapack ===============================================
355
356 #ifdef CHECK_CXXLAPACK
357
358 //-- ev: workspace query
359 template <typename MA>
360 typename GeMatrix<MA>::IndexType
361 ev_native_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
362 {
363 typedef typename GeMatrix<MA>::ElementType T;
364
365 const char JOBVL = (computeVL) ? 'V' : 'N';
366 const char JOBVR = (computeVR) ? 'V' : 'N';
367 const INTEGER N = A.numRows();
368 const INTEGER LDA = A.leadingDimension();
369 const INTEGER LDVL = (JOBVL=='V') ? LDA : 1;
370 const INTEGER LDVR = (JOBVR=='V') ? LDA : 1;
371 T WORK = 0;
372 T DUMMY = 0;
373 const INTEGER LWORK = -1;
374 INTEGER INFO;
375
376 if (IsSame<T,DOUBLE>::value) {
377 LAPACK_IMPL(dgeev)(&JOBVL,
378 &JOBVR,
379 &N,
380 A.data(),
381 &LDA,
382 &DUMMY,
383 &DUMMY,
384 &DUMMY,
385 &LDVL,
386 &DUMMY,
387 &LDVR,
388 &WORK,
389 &LWORK,
390 &INFO);
391 } else {
392 ASSERT(0);
393 }
394 ASSERT(INFO>=0);
395 return WORK;
396 }
397
398 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
399 typename VWORK>
400 typename GeMatrix<MA>::IndexType
401 ev_native(bool computeVL,
402 bool computeVR,
403 GeMatrix<MA> &A,
404 DenseVector<VWR> &wr,
405 DenseVector<VWI> &wi,
406 GeMatrix<MVL> &VL,
407 GeMatrix<MVR> &VR,
408 DenseVector<VWORK> &work)
409 {
410 typedef typename GeMatrix<MA>::ElementType T;
411
412 const char JOBVL = (computeVL) ? 'V' : 'N';
413 const char JOBVR = (computeVR) ? 'V' : 'N';
414 const INTEGER N = A.numRows();
415 const INTEGER LDA = A.leadingDimension();
416 const INTEGER LDVL = VL.leadingDimension();
417 const INTEGER LDVR = VR.leadingDimension();
418 const INTEGER LWORK = work.length();
419 INTEGER INFO;
420
421 if (IsSame<T,DOUBLE>::value) {
422 LAPACK_IMPL(dgeev)(&JOBVL,
423 &JOBVR,
424 &N,
425 A.data(),
426 &LDA,
427 wr.data(),
428 wi.data(),
429 VL.data(),
430 &LDVL,
431 VR.data(),
432 &LDVR,
433 work.data(),
434 &LWORK,
435 &INFO);
436 } else {
437 ASSERT(0);
438 }
439
440 ASSERT(INFO>=0);
441 return INFO;
442 }
443
444 #endif // CHECK_CXXLAPACK
445
446 //== public interface ==========================================================
447
448 template <typename MA>
449 Pair<typename GeMatrix<MA>::IndexType>
450 ev_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
451 {
452 LAPACK_DEBUG_OUT("ev_wsq");
453
454 //
455 // Test the input parameters
456 //
457 # ifndef NDEBUG
458 ASSERT(A.numRows()==A.numCols());
459 ASSERT(A.firstRow()==1);
460 ASSERT(A.firstCol()==1);
461 # endif
462
463 //
464 // Call implementation
465 //
466 const auto ws = ev_generic_wsq(computeVL, computeVR, A);
467
468 # ifdef CHECK_CXXLAPACK
469 //
470 // Compare results
471 //
472 auto optWorkSize = ev_native_wsq(computeVL, computeVR, A);
473 if (! isIdentical(optWorkSize, ws.second, "optWorkSize", "ws.second")) {
474 ASSERT(0);
475 }
476 # endif
477
478 return ws;
479 }
480
481 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
482 typename VWORK>
483 typename GeMatrix<MA>::IndexType
484 ev(bool computeVL,
485 bool computeVR,
486 GeMatrix<MA> &A,
487 DenseVector<VWR> &wr,
488 DenseVector<VWI> &wi,
489 GeMatrix<MVL> &VL,
490 GeMatrix<MVR> &VR,
491 DenseVector<VWORK> &work)
492 {
493 LAPACK_DEBUG_OUT("ev");
494
495 typedef typename GeMatrix<MA>::IndexType IndexType;
496
497 const IndexType n = A.numRows();
498
499 //
500 // Test the input parameters
501 //
502 # ifndef NDEBUG
503 ASSERT(A.numRows()==A.numCols());
504 ASSERT(A.firstRow()==1);
505 ASSERT(A.firstCol()==1);
506 ASSERT(work.firstIndex()==1);
507
508
509 ASSERT(wr.firstIndex()==1);
510 ASSERT(wr.length()==0 || wr.length()==n);
511
512 ASSERT(wi.firstIndex()==1);
513 ASSERT(wi.length()==0 || wi.length()==n);
514
515 if (computeVL) {
516 ASSERT(VL.numRows()==VL.numCols());
517 ASSERT(VL.numRows()==0 || VL.numRows()==n);
518 ASSERT(VL.firstRow()==1);
519 ASSERT(VL.firstCol()==1);
520 }
521
522 if (computeVR) {
523 ASSERT(VR.numRows()==VR.numCols());
524 ASSERT(VR.numRows()==0 || VR.numRows()==n);
525 ASSERT(VR.firstRow()==1);
526 ASSERT(VR.firstCol()==1);
527 }
528 # endif
529
530 //
531 // Resize output arguments if they are empty and needed
532 //
533 if (wr.length()==0) {
534 wr.resize(n);
535 }
536 if (wi.length()==0) {
537 wi.resize(n);
538 }
539 if (computeVL && VL.numRows()==0) {
540 VL.resize(n, n);
541 }
542 if (computeVR && VR.numRows()==0) {
543 VR.resize(n, n);
544 }
545
546 //
547 // Make copies of output arguments
548 //
549 # ifdef CHECK_CXXLAPACK
550 typename GeMatrix<MA>::NoView A_org = A;
551
552 typename GeMatrix<MA>::NoView _A = A;
553 typename DenseVector<VWR>::NoView _wr = wr;
554 typename DenseVector<VWI>::NoView _wi = wi;
555 typename GeMatrix<MVL>::NoView _VL = VL;
556 typename GeMatrix<MVR>::NoView _VR = VR;
557 typename DenseVector<VWORK>::NoView _work = work;
558 # endif
559
560 //
561 // Call implementation
562 //
563 IndexType result = ev_generic(computeVL, computeVR,
564 A, wr, wi, VL, VR,
565 work);
566 # ifdef CHECK_CXXLAPACK
567 //
568 // Compare results
569 //
570 if (_work.length()==0) {
571 _work.resize(work.length());
572 }
573 IndexType _result = ev_native(computeVL, computeVR,
574 _A, _wr, _wi, _VL, _VR,
575 _work);
576
577 bool failed = false;
578 if (! isIdentical(A, _A, " A", "_A")) {
579 std::cerr << "CXXLAPACK: A = " << A << std::endl;
580 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
581 failed = true;
582 }
583
584 if (! isIdentical(wr, _wr, " wr", "_wr")) {
585 std::cerr << "CXXLAPACK: wr = " << wr << std::endl;
586 std::cerr << "F77LAPACK: _wr = " << _wr << std::endl;
587 failed = true;
588 }
589
590 if (! isIdentical(wi, _wi, " wi", "_wi")) {
591 std::cerr << "CXXLAPACK: wi = " << wi << std::endl;
592 std::cerr << "F77LAPACK: _wi = " << _wi << std::endl;
593 failed = true;
594 }
595
596 if (! isIdentical(VL, _VL, " VL", "_VL")) {
597 std::cerr << "CXXLAPACK: VL = " << VL << std::endl;
598 std::cerr << "F77LAPACK: _VL = " << _VL << std::endl;
599 failed = true;
600 }
601
602 if (! isIdentical(VR, _VR, " VR", "_VR")) {
603 std::cerr << "CXXLAPACK: VR = " << VR << std::endl;
604 std::cerr << "F77LAPACK: _VR = " << _VR << std::endl;
605 failed = true;
606 }
607
608 if (! isIdentical(work, _work, " work", "_work")) {
609 std::cerr << "CXXLAPACK: work = " << work << std::endl;
610 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
611 failed = true;
612 }
613
614 if (! isIdentical(result, _result, " result", "_result")) {
615 std::cerr << "CXXLAPACK: result = " << result << std::endl;
616 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
617 failed = true;
618 }
619
620 if (failed) {
621 std::cerr << "-- List of all matrices/vectors --------------------"
622 << std::endl;
623
624 std::cerr << "ORIGINAL A = " << A_org << std::endl;
625
626 std::cerr << "CXXLAPACK: A = " << A << std::endl;
627 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
628
629 std::cerr << "CXXLAPACK: wr = " << wr << std::endl;
630 std::cerr << "F77LAPACK: _wr = " << _wr << std::endl;
631
632 std::cerr << "CXXLAPACK: wi = " << wi << std::endl;
633 std::cerr << "F77LAPACK: _wi = " << _wi << std::endl;
634
635 std::cerr << "CXXLAPACK: VL = " << VL << std::endl;
636 std::cerr << "F77LAPACK: _VL = " << _VL << std::endl;
637
638 std::cerr << "CXXLAPACK: VR = " << VR << std::endl;
639 std::cerr << "F77LAPACK: _VR = " << _VR << std::endl;
640
641 std::cerr << "CXXLAPACK: work = " << work << std::endl;
642 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
643
644 std::cerr << "CXXLAPACK: result = " << result << std::endl;
645 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
646
647 std::cerr << "---------------------------------------" << std::endl;
648 ASSERT(0);
649 } else {
650 // std::cerr << "passed: ev.tcc" << std::endl;
651 }
652
653 # endif
654
655 return result;
656 }
657
658 //-- forwarding ----------------------------------------------------------------
659 template <typename MA>
660 Pair<typename MA::IndexType>
661 ev_wsq(bool computeVL, bool computeVR, MA &&A)
662 {
663 CHECKPOINT_ENTER;
664 const Pair<typename MA::IndexType> result = ev_wsq(computeVL, computeVR, A);
665 CHECKPOINT_LEAVE;
666
667 return result;
668 }
669
670 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
671 typename VWORK>
672 typename MA::IndexType
673 ev(bool computeVL, bool computeVR,
674 MA &&A, VWR &&wr, VWI &&wi,
675 MVL &&VL, MVR &&VR,
676 VWORK &&work)
677 {
678 typedef typename MA::IndexType IndexType;
679
680 CHECKPOINT_ENTER;
681 const IndexType info = ev(computeVL, computeVR, A, wr, wi, VL, VR, work);
682 CHECKPOINT_LEAVE;
683
684 return info;
685 }
686
687 } } // namespace lapack, flens
688
689 #endif // FLENS_LAPACK_EIG_EV_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_EV_TCC
45 #define FLENS_LAPACK_EIG_EV_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
54 //-- ev: workspace query
55 template <typename MA>
56 Pair<typename GeMatrix<MA>::IndexType>
57 ev_generic_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
58 {
59 using std::max;
60
61 typedef typename GeMatrix<MA>::ElementType T;
62 typedef typename GeMatrix<MA>::IndexType IndexType;
63
64 const IndexType n = A.numRows();
65
66 IndexType minWork, maxWork;
67
68 if (n==0) {
69 minWork = 1;
70 maxWork = 1;
71 } else {
72 maxWork = 2*n + n*ilaenv<T>(1, "GEHRD", "", n, 1, n, 0);
73 if (computeVL) {
74 minWork = 4*n;
75 maxWork = max(maxWork,
76 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
77 IndexType hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
78 IndexType(1), n, A);
79 maxWork = max(max(maxWork, n+1), n+hseqrWork);
80 maxWork = max(maxWork, 4*n);
81 } else if (computeVR) {
82 minWork = 4*n;
83 maxWork = max(maxWork,
84 2*n + (n-1)*ilaenv<T>(1, "ORGHR", "", n, 1, n));
85 IndexType hseqrWork = hseqr_wsq(HSEQR::Schur, HSEQR::NoInit,
86 IndexType(1), n, A);
87 maxWork = max(max(maxWork, n+1), n+hseqrWork);
88 maxWork = max(maxWork, 4*n);
89 } else {
90 minWork = 3*n;
91 IndexType hseqrWork = hseqr_wsq(HSEQR::Eigenvalues, HSEQR::No,
92 IndexType(1), n, A);
93 maxWork = max(max(maxWork, n+1), n+hseqrWork);
94 }
95 maxWork = max(maxWork, minWork);
96 }
97 return Pair<IndexType>(minWork, maxWork);
98 }
99
100 //-- ev: computation
101 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
102 typename VWORK>
103 typename GeMatrix<MA>::IndexType
104 ev_generic(bool computeVL, bool computeVR,
105 GeMatrix<MA> &A,
106 DenseVector<VWR> &wr, DenseVector<VWI> &wi,
107 GeMatrix<MVL> &VL, GeMatrix<MVR> &VR,
108 DenseVector<VWORK> &work)
109 {
110 using std::pow;
111 using std::sqrt;
112
113 typedef typename GeMatrix<MA>::ElementType T;
114 typedef typename GeMatrix<MA>::IndexType IndexType;
115
116 const Underscore<IndexType> _;
117 const IndexType n = A.numRows();
118 const T Zero(0), One(1);
119 //
120 // Compute workspace
121 // (Note: Comments in the code beginning "Workspace:" describe the
122 // minimal amount of workspace needed at that point in the code,
123 // as well as the preferred amount for good performance.
124 // NB refers to the optimal block size for the immediately
125 // following subroutine, as returned by ILAENV.
126 // HSWORK refers to the workspace preferred by DHSEQR, as
127 // calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
128 // the worst case.)
129 //
130 Pair<IndexType> wsQuery = ev_wsq(computeVL, computeVR, A);
131 IndexType minWork = wsQuery.first;
132 IndexType maxWork = wsQuery.second;
133
134 if (work.length()!=0 && work.length()<minWork) {
135 ASSERT(0);
136 } else if (work.length()==0) {
137 work.resize(maxWork);
138 }
139
140 work(1) = maxWork;
141
142 //
143 // Quick return if possible
144 //
145 if (n==0) {
146 return 0;
147 }
148 //
149 // Get machine constants
150 //
151 const T Eps = lamch<T>(Precision);
152 T smallNum = lamch<T>(SafeMin);
153 T bigNum = One / smallNum;
154 labad(smallNum, bigNum);
155 smallNum = sqrt(smallNum) / Eps;
156 bigNum = One / smallNum;
157 //
158 // Scale A if max element outside range [SMLNUM,BIGNUM]
159 //
160 const T ANorm = lan(MaximumNorm, A);
161 bool scaleA = false;
162 T cScale;
163 if (ANorm>Zero && ANorm<smallNum) {
164 scaleA = true;
165 cScale = smallNum;
166 } else if (ANorm>bigNum) {
167 scaleA = true;
168 cScale = bigNum;
169 }
170 if (scaleA) {
171 lascl(LASCL::FullMatrix, 0, 0, ANorm, cScale, A);
172 }
173 //
174 // Balance the matrix
175 // (Workspace: need N)
176 //
177 IndexType iBal = 1;
178 IndexType iLo, iHi;
179
180 auto balWork = work(_(iBal,iBal+n-1));
181 bal(BALANCE::Both, A, iLo, iHi, balWork);
182
183 //
184 // Reduce to upper Hessenberg form
185 // (Workspace: need 3*N, prefer 2*N+N*NB)
186 //
187 IndexType iTau = iBal + n;
188 IndexType iWork = iTau + n;
189 IndexType lWork = work.length();
190
191 auto tau = work(_(iTau, iTau+n-2));
192 auto hrdWork = work(_(iWork, lWork));
193
194 hrd(iLo, iHi, A, tau, hrdWork);
195
196 IndexType info = 0;
197 if (computeVL) {
198 //
199 // Want left eigenvectors
200 // Copy Householder vectors to VL
201 //
202 VL.lower() = A.lower();
203 //
204 // Generate orthogonal matrix in VL
205 // (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
206 //
207 orghr(iLo, iHi, VL, tau, hrdWork);
208 //
209 // Perform QR iteration, accumulating Schur vectors in VL
210 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
211 //
212 iWork = iTau;
213 auto hseqrWork = work(_(iWork, lWork));
214 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
215 wr, wi, VL, hseqrWork);
216
217 if (computeVR) {
218 //
219 // Want left and right eigenvectors
220 // Copy Schur vectors to VR
221 //
222 VR = VL;
223 }
224
225 } else if (computeVR) {
226 //
227 // Want right eigenvectors
228 // Copy Householder vectors to VR
229 //
230 VR.lower() = A.lower();
231
232 //
233 // Generate orthogonal matrix in VR
234 // (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
235 //
236 orghr(iLo, iHi, VR, tau, hrdWork);
237 //
238 // Perform QR iteration, accumulating Schur vectors in VR
239 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
240 //
241 iWork = iTau;
242 auto hseqrWork = work(_(iWork, lWork));
243 info = hseqr(HSEQR::Schur, HSEQR::NoInit, iLo, iHi, A,
244 wr, wi, VR, hseqrWork);
245
246 } else {
247 //
248 // Compute eigenvalues only
249 // (Workspace: need N+1, prefer N+HSWORK (see comments) )
250 //
251 iWork = iTau;
252 auto hseqrWork = work(_(iWork, lWork));
253 info = hseqr(HSEQR::Eigenvalues, HSEQR::No, iLo, iHi, A,
254 wr, wi, VR, hseqrWork);
255 }
256
257 //
258 // If INFO > 0 from DHSEQR, then quit
259 //
260 if (info==0) {
261
262 if (computeVL || computeVR) {
263 //
264 // Compute left and/or right eigenvectors
265 // (Workspace: need 4*N)
266 //
267 // TODO: shouldn't that be "Workspace: need 3*N" ?? instead of 4*N?
268 IndexType nOut;
269 auto trevcWork = work(_(iWork, lWork));
270
271 // TODO: find a way that we don't have to pass an empty "select"
272 DenseVector<Array<bool> > select;
273 trevc(computeVL, computeVR, TREVC::Backtransform, select,
274 A, VL, VR, n, nOut, trevcWork);
275 }
276
277 if (computeVL) {
278 //
279 // Undo balancing of left eigenvectors
280 // (Workspace: need N)
281 //
282 bak(BALANCE::Both, Left, iLo, iHi, balWork, VL);
283 //
284 // Normalize left eigenvectors and make largest component real
285 //
286 auto _work = work(_(iWork, iWork+n-1));
287 for (IndexType i=1; i<=n; ++i) {
288 if (wi(i)==Zero) {
289 VL(_,i) *= One / blas::nrm2(VL(_,i));
290 } else if (wi(i)>Zero) {
291 T scale = One / lapy2(blas::nrm2(VL(_,i )),
292 blas::nrm2(VL(_,i+1)));
293 VL(_,i) *= scale;
294 VL(_,i+1) *= scale;
295 for (IndexType k=1; k<=n; ++k) {
296 _work(k) = pow(VL(k,i), 2) + pow(VL(k,i+1), 2);
297 }
298 IndexType k = blas::iamax(_work);
299 T cs, sn, r;
300 lartg(VL(k,i), VL(k,i+1), cs, sn, r);
301 blas::rot(VL(_,i), VL(_,i+1), cs, sn);
302 VL(k,i+1) = Zero;
303 }
304 }
305 }
306
307 if (computeVR) {
308 //
309 // Undo balancing of right eigenvectors
310 // (Workspace: need N)
311 //
312 bak(BALANCE::Both, Right, iLo, iHi, balWork, VR);
313 //
314 // Normalize right eigenvectors and make largest component real
315 //
316 auto _work = work(_(iWork, iWork+n-1));
317 for (IndexType i=1; i<=n; ++i) {
318 if (wi(i)==Zero) {
319 VR(_,i) *= One / blas::nrm2(VR(_,i));
320 } else if (wi(i)>Zero) {
321 T scale = One / lapy2(blas::nrm2(VR(_,i )),
322 blas::nrm2(VR(_,i+1)));
323 VR(_,i) *= scale;
324 VR(_,i+1) *= scale;
325 for (IndexType k=1; k<=n; ++k) {
326 _work(k) = pow(VR(k,i), 2) + pow(VR(k,i+1), 2);
327 }
328 IndexType k = blas::iamax(_work);
329 T cs, sn, r;
330 lartg(VR(k,i), VR(k,i+1), cs, sn, r);
331 blas::rot(VR(_,i), VR(_,i+1), cs, sn);
332 VR(k,i+1) = Zero;
333 }
334 }
335 }
336 }
337 //
338 // Undo scaling if necessary
339 //
340 if (scaleA) {
341 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(info+1,n)));
342 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(info+1,n)));
343
344 if (info>0) {
345 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wr(_(1,iLo-1)));
346 lascl(LASCL::FullMatrix, 0, 0, cScale, ANorm, wi(_(1,iLo-1)));
347 }
348 }
349
350 work(1) = maxWork;
351 return info;
352 }
353
354 //== interface for native lapack ===============================================
355
356 #ifdef CHECK_CXXLAPACK
357
358 //-- ev: workspace query
359 template <typename MA>
360 typename GeMatrix<MA>::IndexType
361 ev_native_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
362 {
363 typedef typename GeMatrix<MA>::ElementType T;
364
365 const char JOBVL = (computeVL) ? 'V' : 'N';
366 const char JOBVR = (computeVR) ? 'V' : 'N';
367 const INTEGER N = A.numRows();
368 const INTEGER LDA = A.leadingDimension();
369 const INTEGER LDVL = (JOBVL=='V') ? LDA : 1;
370 const INTEGER LDVR = (JOBVR=='V') ? LDA : 1;
371 T WORK = 0;
372 T DUMMY = 0;
373 const INTEGER LWORK = -1;
374 INTEGER INFO;
375
376 if (IsSame<T,DOUBLE>::value) {
377 LAPACK_IMPL(dgeev)(&JOBVL,
378 &JOBVR,
379 &N,
380 A.data(),
381 &LDA,
382 &DUMMY,
383 &DUMMY,
384 &DUMMY,
385 &LDVL,
386 &DUMMY,
387 &LDVR,
388 &WORK,
389 &LWORK,
390 &INFO);
391 } else {
392 ASSERT(0);
393 }
394 ASSERT(INFO>=0);
395 return WORK;
396 }
397
398 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
399 typename VWORK>
400 typename GeMatrix<MA>::IndexType
401 ev_native(bool computeVL,
402 bool computeVR,
403 GeMatrix<MA> &A,
404 DenseVector<VWR> &wr,
405 DenseVector<VWI> &wi,
406 GeMatrix<MVL> &VL,
407 GeMatrix<MVR> &VR,
408 DenseVector<VWORK> &work)
409 {
410 typedef typename GeMatrix<MA>::ElementType T;
411
412 const char JOBVL = (computeVL) ? 'V' : 'N';
413 const char JOBVR = (computeVR) ? 'V' : 'N';
414 const INTEGER N = A.numRows();
415 const INTEGER LDA = A.leadingDimension();
416 const INTEGER LDVL = VL.leadingDimension();
417 const INTEGER LDVR = VR.leadingDimension();
418 const INTEGER LWORK = work.length();
419 INTEGER INFO;
420
421 if (IsSame<T,DOUBLE>::value) {
422 LAPACK_IMPL(dgeev)(&JOBVL,
423 &JOBVR,
424 &N,
425 A.data(),
426 &LDA,
427 wr.data(),
428 wi.data(),
429 VL.data(),
430 &LDVL,
431 VR.data(),
432 &LDVR,
433 work.data(),
434 &LWORK,
435 &INFO);
436 } else {
437 ASSERT(0);
438 }
439
440 ASSERT(INFO>=0);
441 return INFO;
442 }
443
444 #endif // CHECK_CXXLAPACK
445
446 //== public interface ==========================================================
447
448 template <typename MA>
449 Pair<typename GeMatrix<MA>::IndexType>
450 ev_wsq(bool computeVL, bool computeVR, GeMatrix<MA> &A)
451 {
452 LAPACK_DEBUG_OUT("ev_wsq");
453
454 //
455 // Test the input parameters
456 //
457 # ifndef NDEBUG
458 ASSERT(A.numRows()==A.numCols());
459 ASSERT(A.firstRow()==1);
460 ASSERT(A.firstCol()==1);
461 # endif
462
463 //
464 // Call implementation
465 //
466 const auto ws = ev_generic_wsq(computeVL, computeVR, A);
467
468 # ifdef CHECK_CXXLAPACK
469 //
470 // Compare results
471 //
472 auto optWorkSize = ev_native_wsq(computeVL, computeVR, A);
473 if (! isIdentical(optWorkSize, ws.second, "optWorkSize", "ws.second")) {
474 ASSERT(0);
475 }
476 # endif
477
478 return ws;
479 }
480
481 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
482 typename VWORK>
483 typename GeMatrix<MA>::IndexType
484 ev(bool computeVL,
485 bool computeVR,
486 GeMatrix<MA> &A,
487 DenseVector<VWR> &wr,
488 DenseVector<VWI> &wi,
489 GeMatrix<MVL> &VL,
490 GeMatrix<MVR> &VR,
491 DenseVector<VWORK> &work)
492 {
493 LAPACK_DEBUG_OUT("ev");
494
495 typedef typename GeMatrix<MA>::IndexType IndexType;
496
497 const IndexType n = A.numRows();
498
499 //
500 // Test the input parameters
501 //
502 # ifndef NDEBUG
503 ASSERT(A.numRows()==A.numCols());
504 ASSERT(A.firstRow()==1);
505 ASSERT(A.firstCol()==1);
506 ASSERT(work.firstIndex()==1);
507
508
509 ASSERT(wr.firstIndex()==1);
510 ASSERT(wr.length()==0 || wr.length()==n);
511
512 ASSERT(wi.firstIndex()==1);
513 ASSERT(wi.length()==0 || wi.length()==n);
514
515 if (computeVL) {
516 ASSERT(VL.numRows()==VL.numCols());
517 ASSERT(VL.numRows()==0 || VL.numRows()==n);
518 ASSERT(VL.firstRow()==1);
519 ASSERT(VL.firstCol()==1);
520 }
521
522 if (computeVR) {
523 ASSERT(VR.numRows()==VR.numCols());
524 ASSERT(VR.numRows()==0 || VR.numRows()==n);
525 ASSERT(VR.firstRow()==1);
526 ASSERT(VR.firstCol()==1);
527 }
528 # endif
529
530 //
531 // Resize output arguments if they are empty and needed
532 //
533 if (wr.length()==0) {
534 wr.resize(n);
535 }
536 if (wi.length()==0) {
537 wi.resize(n);
538 }
539 if (computeVL && VL.numRows()==0) {
540 VL.resize(n, n);
541 }
542 if (computeVR && VR.numRows()==0) {
543 VR.resize(n, n);
544 }
545
546 //
547 // Make copies of output arguments
548 //
549 # ifdef CHECK_CXXLAPACK
550 typename GeMatrix<MA>::NoView A_org = A;
551
552 typename GeMatrix<MA>::NoView _A = A;
553 typename DenseVector<VWR>::NoView _wr = wr;
554 typename DenseVector<VWI>::NoView _wi = wi;
555 typename GeMatrix<MVL>::NoView _VL = VL;
556 typename GeMatrix<MVR>::NoView _VR = VR;
557 typename DenseVector<VWORK>::NoView _work = work;
558 # endif
559
560 //
561 // Call implementation
562 //
563 IndexType result = ev_generic(computeVL, computeVR,
564 A, wr, wi, VL, VR,
565 work);
566 # ifdef CHECK_CXXLAPACK
567 //
568 // Compare results
569 //
570 if (_work.length()==0) {
571 _work.resize(work.length());
572 }
573 IndexType _result = ev_native(computeVL, computeVR,
574 _A, _wr, _wi, _VL, _VR,
575 _work);
576
577 bool failed = false;
578 if (! isIdentical(A, _A, " A", "_A")) {
579 std::cerr << "CXXLAPACK: A = " << A << std::endl;
580 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
581 failed = true;
582 }
583
584 if (! isIdentical(wr, _wr, " wr", "_wr")) {
585 std::cerr << "CXXLAPACK: wr = " << wr << std::endl;
586 std::cerr << "F77LAPACK: _wr = " << _wr << std::endl;
587 failed = true;
588 }
589
590 if (! isIdentical(wi, _wi, " wi", "_wi")) {
591 std::cerr << "CXXLAPACK: wi = " << wi << std::endl;
592 std::cerr << "F77LAPACK: _wi = " << _wi << std::endl;
593 failed = true;
594 }
595
596 if (! isIdentical(VL, _VL, " VL", "_VL")) {
597 std::cerr << "CXXLAPACK: VL = " << VL << std::endl;
598 std::cerr << "F77LAPACK: _VL = " << _VL << std::endl;
599 failed = true;
600 }
601
602 if (! isIdentical(VR, _VR, " VR", "_VR")) {
603 std::cerr << "CXXLAPACK: VR = " << VR << std::endl;
604 std::cerr << "F77LAPACK: _VR = " << _VR << std::endl;
605 failed = true;
606 }
607
608 if (! isIdentical(work, _work, " work", "_work")) {
609 std::cerr << "CXXLAPACK: work = " << work << std::endl;
610 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
611 failed = true;
612 }
613
614 if (! isIdentical(result, _result, " result", "_result")) {
615 std::cerr << "CXXLAPACK: result = " << result << std::endl;
616 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
617 failed = true;
618 }
619
620 if (failed) {
621 std::cerr << "-- List of all matrices/vectors --------------------"
622 << std::endl;
623
624 std::cerr << "ORIGINAL A = " << A_org << std::endl;
625
626 std::cerr << "CXXLAPACK: A = " << A << std::endl;
627 std::cerr << "F77LAPACK: _A = " << _A << std::endl;
628
629 std::cerr << "CXXLAPACK: wr = " << wr << std::endl;
630 std::cerr << "F77LAPACK: _wr = " << _wr << std::endl;
631
632 std::cerr << "CXXLAPACK: wi = " << wi << std::endl;
633 std::cerr << "F77LAPACK: _wi = " << _wi << std::endl;
634
635 std::cerr << "CXXLAPACK: VL = " << VL << std::endl;
636 std::cerr << "F77LAPACK: _VL = " << _VL << std::endl;
637
638 std::cerr << "CXXLAPACK: VR = " << VR << std::endl;
639 std::cerr << "F77LAPACK: _VR = " << _VR << std::endl;
640
641 std::cerr << "CXXLAPACK: work = " << work << std::endl;
642 std::cerr << "F77LAPACK: _work = " << _work << std::endl;
643
644 std::cerr << "CXXLAPACK: result = " << result << std::endl;
645 std::cerr << "F77LAPACK: _result = " << _result << std::endl;
646
647 std::cerr << "---------------------------------------" << std::endl;
648 ASSERT(0);
649 } else {
650 // std::cerr << "passed: ev.tcc" << std::endl;
651 }
652
653 # endif
654
655 return result;
656 }
657
658 //-- forwarding ----------------------------------------------------------------
659 template <typename MA>
660 Pair<typename MA::IndexType>
661 ev_wsq(bool computeVL, bool computeVR, MA &&A)
662 {
663 CHECKPOINT_ENTER;
664 const Pair<typename MA::IndexType> result = ev_wsq(computeVL, computeVR, A);
665 CHECKPOINT_LEAVE;
666
667 return result;
668 }
669
670 template <typename MA, typename VWR, typename VWI, typename MVL, typename MVR,
671 typename VWORK>
672 typename MA::IndexType
673 ev(bool computeVL, bool computeVR,
674 MA &&A, VWR &&wr, VWI &&wi,
675 MVL &&VL, MVR &&VR,
676 VWORK &&work)
677 {
678 typedef typename MA::IndexType IndexType;
679
680 CHECKPOINT_ENTER;
681 const IndexType info = ev(computeVL, computeVR, A, wr, wi, VL, VR, work);
682 CHECKPOINT_LEAVE;
683
684 return info;
685 }
686
687 } } // namespace lapack, flens
688
689 #endif // FLENS_LAPACK_EIG_EV_TCC