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 DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
36 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
37 *
38 * -- LAPACK 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
45 #ifndef FLENS_LAPACK_EIG_TRSEN_TCC
46 #define FLENS_LAPACK_EIG_TRSEN_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 template <typename SELECT, typename MT, typename IndexType>
56 Pair<IndexType>
57 trsen_generic_wsq(TRSEN::Job job,
58 const DenseVector<SELECT> &select,
59 GeMatrix<MT> &T,
60 IndexType &m)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MT>::ElementType ElementType;
65 const ElementType Zero(0);
66
67 const IndexType n = T.numRows();
68
69 //
70 // Decode and test the input parameters
71 //
72 const bool wantBH = (job==TRSEN::Both);
73 const bool wantSP = (job==TRSEN::InvariantSubspaceOnly) || wantBH;
74 //
75 // Set M to the dimension of the specified invariant subspace,
76 // and test LWORK and LIWORK.
77 //
78 m = 0;
79 bool pair = false;
80 for (IndexType k=1; k<=n; ++k) {
81 if (pair) {
82 pair = false;
83 } else {
84 if (k<n) {
85 if (T(k+1,k)==Zero) {
86 if (select(k)) {
87 ++m;
88 }
89 } else {
90 pair = true;
91 if (select(k) || select(k+1)) {
92 m += 2;
93 }
94 }
95 } else {
96 if (select(n)) {
97 ++m;
98 }
99 }
100 }
101 }
102
103 IndexType n1 = m;
104 IndexType n2 = n - m;
105 IndexType nn = n1*n2;
106
107 IndexType workMin = 0,
108 iWorkMin = 0;
109
110 if (wantSP) {
111 workMin = max(IndexType(1), 2*nn);
112 iWorkMin = max(IndexType(1), nn);
113 } else if (job==TRSEN::None) {
114 workMin = max(IndexType(1), n);
115 iWorkMin = 1;
116 } else if (job==TRSEN::EigenvaluesOnly) {
117 workMin = max(IndexType(1), nn);
118 iWorkMin = 1;
119 }
120
121 return Pair<IndexType>(workMin, iWorkMin);
122 }
123
124 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
125 typename IndexType, typename S, typename SEP,
126 typename WORK, typename IWORK>
127 IndexType
128 trsen_generic(TRSEN::Job job,
129 bool computeQ,
130 const DenseVector<SELECT> &select,
131 GeMatrix<MT> &T,
132 GeMatrix<MQ> &Q,
133 DenseVector<WR> &wr,
134 DenseVector<WI> &wi,
135 IndexType &m,
136 S &s,
137 SEP &sep,
138 DenseVector<WORK> &work,
139 DenseVector<IWORK> &iwork)
140 {
141 using std::abs;
142
143 typedef typename GeMatrix<MT>::ElementType ElementType;
144 const ElementType Zero(0), One(1);
145
146 const IndexType n = T.numRows();
147 const Underscore<IndexType> _;
148
149 IndexType info = 0;
150
151 //
152 // .. Local Arrays ..
153 // this array is used to save variables between calls to lacn2
154 //
155 IndexType _isaveData[3] = {0, 0, 0};
156 DenseVectorView<IndexType>
157 isave = typename DenseVectorView<IndexType>::Engine(3, _isaveData);
158 //
159 // Decode and test the input parameters
160 //
161 const bool wantBH = (job==TRSEN::Both);
162 const bool wantS = (job==TRSEN::EigenvaluesOnly) || wantBH;
163 const bool wantSP = (job==TRSEN::InvariantSubspaceOnly) || wantBH;
164 //
165 // Set M to the dimension of the specified invariant subspace,
166 // and compute minimal size of WORK and IWORK
167 //
168 auto wsq = trsen_generic_wsq(job, select, T, m);
169 const IndexType workMin = wsq.first;
170 const IndexType iWorkMin = wsq.second;
171
172 IndexType lWork = work.length();
173 IndexType liWork = iwork.length();
174
175 if (lWork<workMin && lWork>0) {
176 ASSERT(0);
177 return -1;
178 } else if (liWork<iWorkMin && liWork>0) {
179 ASSERT(0);
180 return -2;
181 }
182
183 if (lWork==0) {
184 work.resize(workMin);
185 }
186 if (liWork==0) {
187 iwork.resize(iWorkMin);
188 }
189
190 work(1) = workMin;
191 iwork(1) = iWorkMin;
192
193 IndexType n1 = m;
194 IndexType n2 = n - m;
195 IndexType nn = n1*n2;
196
197 bool pair;
198 IndexType ks;
199 ElementType scale;
200
201 //
202 // Quick return if possible.
203 //
204 if (m==n || m==0) {
205 if (wantS) {
206 s = One;
207 }
208 if (wantSP) {
209 sep = lan(OneNorm, T);
210 }
211 goto DONE;
212 }
213 //
214 // Collect the selected blocks at the top-left corner of T.
215 //
216 ks = 0;
217 pair = false;
218 for (IndexType k=1; k<=n; ++k) {
219 if (pair) {
220 pair = false;
221 } else {
222 bool swap = select(k);
223 if (k<n) {
224 if (T(k+1,k)!=Zero) {
225 pair = true;
226 swap = swap || select(k+1);
227 }
228 }
229 if (swap) {
230 ++ks;
231 //
232 // Swap the K-th block to position KS.
233 //
234 IndexType iErr = 0;
235 IndexType kk = k;
236 if (k!=ks) {
237 trexc(computeQ, T, Q, kk, ks, work(_(1,n)));
238 }
239 if (iErr==1 || iErr==2) {
240 //
241 // Blocks too close to swap: exit.
242 //
243 info = 1;
244 if (wantS) {
245 s = Zero;
246 }
247 if (wantSP) {
248 sep = Zero;
249 }
250 goto DONE;
251 }
252 if (pair) {
253 ++ks;
254 }
255 }
256 }
257 }
258
259 if (wantS) {
260 GeMatrixView<ElementType> Work(n1, n2, work, n1);
261 //
262 // Solve Sylvester equation for R:
263 //
264 // T11*R - R*T22 = scale*T12
265 //
266 Work = T(_(1,n1),_(n1+1,n1+n2));
267 auto T11 = T(_(1,n1),_(1,n1));
268 auto T22 = T(_(n1+1,n1+n2),_(n1+1,n1+n2));
269
270 trsyl(NoTrans, NoTrans, -1, T11, T22, Work, scale);
271 //
272 // Estimate the reciprocal of the condition number of the cluster
273 // of eigenvalues.
274 //
275 ElementType rNorm = lan(FrobeniusNorm, Work);
276 if (rNorm==Zero) {
277 s = One;
278 } else {
279 s = scale / (sqrt(scale*scale/rNorm + rNorm)*sqrt(rNorm));
280 }
281 }
282
283 if (wantSP) {
284 //
285 // Estimate sep(T11,T22).
286 //
287 ElementType est = Zero;
288 IndexType kase = 0;
289
290 do {
291 auto _v = work(_(nn+1,nn+nn));
292 auto _x = work(_(1,nn));
293 auto _isgn = iwork(_(1,nn));
294
295 lacn2(_v, _x, _isgn, est, kase, isave);
296 if (kase==0) {
297 break;
298 }
299
300 GeMatrixView<ElementType> Work(n1, n2, work, n1);
301 auto T11 = T(_(1,n1),_(1,n1));
302 auto T22 = T(_(n1+1,n1+n2),_(n1+1,n1+n2));
303
304 if (kase==1) {
305 //
306 // Solve T11*R - R*T22 = scale*X.
307 //
308 trsyl(NoTrans, NoTrans, -1, T11, T22, Work, scale);
309 } else {
310 //
311 // Solve T11**T*R - R*T22**T = scale*X.
312 //
313 trsyl(Trans, Trans, -1, T11, T22, Work, scale);
314 }
315 } while (true);
316
317 sep = scale / est;
318 }
319
320 DONE:
321
322 //
323 // Store the output eigenvalues in WR and WI.
324 //
325 for (IndexType k=1; k<=n; ++k) {
326 wr(k) = T(k,k);
327 wi(k) = Zero;
328 }
329 for (IndexType k=1; k<=n-1; ++k) {
330 if (T(k+1,k)!=Zero) {
331 wi(k) = sqrt(abs(T(k,k+1)))*sqrt(abs(T(k+1,k)));
332 wi(k+1) = -wi(k);
333 }
334 }
335
336 work(1) = workMin;
337 iwork(1) = iWorkMin;
338 return info;
339 }
340
341 //== interface for native lapack ===============================================
342
343 #ifdef CHECK_CXXLAPACK
344
345 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
346 typename IndexType, typename S, typename SEP,
347 typename WORK, typename IWORK>
348 IndexType
349 trsen_native(TRSEN::Job job,
350 bool computeQ,
351 const DenseVector<SELECT> &select,
352 GeMatrix<MT> &T,
353 GeMatrix<MQ> &Q,
354 DenseVector<WR> &wr,
355 DenseVector<WI> &wi,
356 IndexType &m,
357 S &s,
358 SEP &sep,
359 DenseVector<WORK> &work,
360 DenseVector<IWORK> &iwork)
361 {
362 typedef typename GeMatrix<MT>::ElementType ElementType;
363
364 const char JOB = job;
365 const char COMPQ = (computeQ) ? 'V' : 'N';
366 const INTEGER N = T.numRows();
367 const INTEGER LDT = T.leadingDimension();
368 const INTEGER LDQ = Q.leadingDimension();
369 INTEGER _M = m;
370 DOUBLE _S = s;
371 DOUBLE _SEP = sep;
372 const INTEGER LWORK = work.length();
373 INTEGER LIWORK = iwork.length();
374 INTEGER INFO;
375
376 ASSERT(JOB=='N' || JOB=='E' || JOB=='V' || JOB=='B');
377
378 DenseVector<Array<LOGICAL> > _select(N);
379 for (IndexType i=1; i<=N; ++i) {
380 _select(i) = select(i);
381 }
382
383 DenseVector<Array<INTEGER> > _iwork(iwork.length());
384 for (IndexType i=1; i<=iwork.length(); ++i) {
385 _iwork(i) = iwork(i);
386 }
387
388 if (IsSame<ElementType,DOUBLE>::value) {
389 LAPACK_IMPL(dtrsen)(&JOB,
390 &COMPQ,
391 _select.data(),
392 &N,
393 T.data(),
394 &LDT,
395 Q.data(),
396 &LDQ,
397 wr.data(),
398 wi.data(),
399 &_M,
400 &_S,
401 &_SEP,
402 work.data(),
403 &LWORK,
404 _iwork.data(),
405 &LIWORK,
406 &INFO);
407 } else {
408 ASSERT(0);
409 }
410 if (INFO<0) {
411 std::cerr << "dtrsen: INFO = " << INFO << std::endl;
412 }
413 ASSERT(INFO>=0);
414
415 m = _M;
416 s = _S;
417 sep = _SEP;
418 for (IndexType i=1; i<=iwork.length(); ++i) {
419 iwork(i) = _iwork(i);
420 }
421
422 return INFO;
423 }
424
425 #endif // CHECK_CXXLAPACK
426
427 //== public interface ==========================================================
428
429 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
430 typename IndexType, typename S, typename SEP,
431 typename WORK, typename IWORK>
432 IndexType
433 trsen(TRSEN::Job job,
434 bool computeQ,
435 const DenseVector<SELECT> &select,
436 GeMatrix<MT> &T,
437 GeMatrix<MQ> &Q,
438 DenseVector<WR> &wr,
439 DenseVector<WI> &wi,
440 IndexType &m,
441 S &s,
442 SEP &sep,
443 DenseVector<WORK> &work,
444 DenseVector<IWORK> &iwork)
445 {
446 LAPACK_DEBUG_OUT("trsen");
447 //
448 // Test the input parameters
449 //
450 # ifndef NDEBUG
451 ASSERT(T.firstRow()==1);
452 ASSERT(T.firstCol()==1);
453 ASSERT(T.numRows()==T.numCols());
454
455 const IndexType n = T.numRows();
456 if (computeQ) {
457 ASSERT(Q.firstRow()==1);
458 ASSERT(Q.firstCol()==1);
459 ASSERT(Q.numRows()==n);
460 ASSERT(Q.numCols()==n);
461 }
462
463 ASSERT(wr.firstIndex()==1);
464 ASSERT(wr.length()==n);
465
466 ASSERT(wi.firstIndex()==1);
467 ASSERT(wi.length()==n);
468 # endif
469
470 # ifdef CHECK_CXXLAPACK
471 //
472 // Make copies of output arguments
473 //
474 const typename GeMatrix<MT>::NoView T_org = T;
475 const typename GeMatrix<MQ>::NoView Q_org = Q;
476 const typename DenseVector<WR>::NoView wr_org = wr;
477 const typename DenseVector<WI>::NoView wi_org = wi;
478 const IndexType m_org = m;
479 const S s_org = s;
480 const SEP sep_org = sep;
481 const typename DenseVector<WORK>::NoView work_org = work;
482 const typename DenseVector<IWORK>::NoView iwork_org = iwork;
483 # endif
484
485 //
486 // Call implementation
487 //
488 IndexType info = trsen_generic(job, computeQ, select, T, Q, wr, wi,
489 m, s, sep, work, iwork);
490
491 # ifdef CHECK_CXXLAPACK
492 //
493 // Make copies of results computed by generic implementation
494 //
495 const typename GeMatrix<MT>::NoView T_generic = T;
496 const typename GeMatrix<MQ>::NoView Q_generic = Q;
497 const typename DenseVector<WR>::NoView wr_generic = wr;
498 const typename DenseVector<WI>::NoView wi_generic = wi;
499 const IndexType m_generic = m;
500 const S s_generic = s;
501 const SEP sep_generic = sep;
502 const typename DenseVector<WORK>::NoView work_generic = work;
503 const typename DenseVector<IWORK>::NoView iwork_generic = iwork;
504
505 //
506 // restore output arguments
507 //
508 T = T_org;
509 Q = Q_org;
510 wr = wr_org;
511 wi = wi_org;
512 m = m_org;
513 s = s_org;
514 sep = sep_org;
515 work = work_org;
516 iwork = iwork_org;
517
518 //
519 // Compare generic results with results from the native implementation
520 //
521 IndexType _info = trsen_native(job, computeQ, select, T, Q, wr, wi,
522 m, s, sep, work, iwork);
523 bool failed = false;
524 if (! isIdentical(T_generic, T, "T_generic", "T")) {
525 std::cerr << "CXXLAPACK: T_generic = " << T_generic << std::endl;
526 std::cerr << "F77LAPACK: T = " << T << std::endl;
527 failed = true;
528 }
529 if (! isIdentical(Q_generic, Q, "Q_generic", "Q")) {
530 std::cerr << "CXXLAPACK: Q_generic = " << Q_generic << std::endl;
531 std::cerr << "F77LAPACK: Q = " << Q << std::endl;
532 failed = true;
533 }
534 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
535 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
536 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
537 failed = true;
538 }
539 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
540 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
541 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
542 failed = true;
543 }
544 if (! isIdentical(m_generic, m, "m_generic", "m")) {
545 std::cerr << "CXXLAPACK: m_generic = " << m_generic << std::endl;
546 std::cerr << "F77LAPACK: m = " << m << std::endl;
547 failed = true;
548 }
549 if (! isIdentical(s_generic, s, "s_generic", "s")) {
550 std::cerr << "CXXLAPACK: s_generic = " << s_generic << std::endl;
551 std::cerr << "F77LAPACK: s = " << s << std::endl;
552 failed = true;
553 }
554 if (! isIdentical(sep_generic, sep, "sep_generic", "sep")) {
555 std::cerr << "CXXLAPACK: sep_generic = " << sep_generic << std::endl;
556 std::cerr << "F77LAPACK: sep = " << sep << std::endl;
557 failed = true;
558 }
559 if (! isIdentical(work_generic, work, "work_generic", "work")) {
560 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
561 std::cerr << "F77LAPACK: work = " << work << std::endl;
562 failed = true;
563 }
564 if (! isIdentical(iwork_generic, iwork, "iwork_generic", "iwork")) {
565 std::cerr << "CXXLAPACK: iwork_generic = "
566 << iwork_generic << std::endl;
567 std::cerr << "F77LAPACK: iwork = " << iwork << std::endl;
568 failed = true;
569 }
570 if (! isIdentical(info, _info, "info", "_info")) {
571 std::cerr << "CXXLAPACK: info = " << info << std::endl;
572 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
573 failed = true;
574 }
575
576 if (failed) {
577 std::cerr << "error in: trsen.tcc" << std::endl;
578 ASSERT(0);
579 } else {
580 // std::cerr << "passed: trsen.tcc" << std::endl;
581 }
582 # endif
583
584 return info;
585 }
586
587 //-- forwarding ----------------------------------------------------------------
588
589 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
590 typename IndexType, typename S, typename SEP,
591 typename WORK, typename IWORK>
592 IndexType
593 trsen(TRSEN::Job job,
594 bool computeQ,
595 const DenseVector<SELECT> &select,
596 MT &&T,
597 MQ &&Q,
598 WR &&wr,
599 WI &&wi,
600 IndexType &&m,
601 S &&s,
602 SEP &&sep,
603 WORK &&work,
604 IWORK &&iwork)
605 {
606 trsen(job, computeQ, select, T, Q, wr, wi, m, s, sep, work, iwork);
607 }
608
609 } } // namespace lapack, flens
610
611 #endif // FLENS_LAPACK_EIG_TRSEN_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 DTRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI,
36 $ M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO )
37 *
38 * -- LAPACK 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
45 #ifndef FLENS_LAPACK_EIG_TRSEN_TCC
46 #define FLENS_LAPACK_EIG_TRSEN_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 template <typename SELECT, typename MT, typename IndexType>
56 Pair<IndexType>
57 trsen_generic_wsq(TRSEN::Job job,
58 const DenseVector<SELECT> &select,
59 GeMatrix<MT> &T,
60 IndexType &m)
61 {
62 using std::max;
63
64 typedef typename GeMatrix<MT>::ElementType ElementType;
65 const ElementType Zero(0);
66
67 const IndexType n = T.numRows();
68
69 //
70 // Decode and test the input parameters
71 //
72 const bool wantBH = (job==TRSEN::Both);
73 const bool wantSP = (job==TRSEN::InvariantSubspaceOnly) || wantBH;
74 //
75 // Set M to the dimension of the specified invariant subspace,
76 // and test LWORK and LIWORK.
77 //
78 m = 0;
79 bool pair = false;
80 for (IndexType k=1; k<=n; ++k) {
81 if (pair) {
82 pair = false;
83 } else {
84 if (k<n) {
85 if (T(k+1,k)==Zero) {
86 if (select(k)) {
87 ++m;
88 }
89 } else {
90 pair = true;
91 if (select(k) || select(k+1)) {
92 m += 2;
93 }
94 }
95 } else {
96 if (select(n)) {
97 ++m;
98 }
99 }
100 }
101 }
102
103 IndexType n1 = m;
104 IndexType n2 = n - m;
105 IndexType nn = n1*n2;
106
107 IndexType workMin = 0,
108 iWorkMin = 0;
109
110 if (wantSP) {
111 workMin = max(IndexType(1), 2*nn);
112 iWorkMin = max(IndexType(1), nn);
113 } else if (job==TRSEN::None) {
114 workMin = max(IndexType(1), n);
115 iWorkMin = 1;
116 } else if (job==TRSEN::EigenvaluesOnly) {
117 workMin = max(IndexType(1), nn);
118 iWorkMin = 1;
119 }
120
121 return Pair<IndexType>(workMin, iWorkMin);
122 }
123
124 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
125 typename IndexType, typename S, typename SEP,
126 typename WORK, typename IWORK>
127 IndexType
128 trsen_generic(TRSEN::Job job,
129 bool computeQ,
130 const DenseVector<SELECT> &select,
131 GeMatrix<MT> &T,
132 GeMatrix<MQ> &Q,
133 DenseVector<WR> &wr,
134 DenseVector<WI> &wi,
135 IndexType &m,
136 S &s,
137 SEP &sep,
138 DenseVector<WORK> &work,
139 DenseVector<IWORK> &iwork)
140 {
141 using std::abs;
142
143 typedef typename GeMatrix<MT>::ElementType ElementType;
144 const ElementType Zero(0), One(1);
145
146 const IndexType n = T.numRows();
147 const Underscore<IndexType> _;
148
149 IndexType info = 0;
150
151 //
152 // .. Local Arrays ..
153 // this array is used to save variables between calls to lacn2
154 //
155 IndexType _isaveData[3] = {0, 0, 0};
156 DenseVectorView<IndexType>
157 isave = typename DenseVectorView<IndexType>::Engine(3, _isaveData);
158 //
159 // Decode and test the input parameters
160 //
161 const bool wantBH = (job==TRSEN::Both);
162 const bool wantS = (job==TRSEN::EigenvaluesOnly) || wantBH;
163 const bool wantSP = (job==TRSEN::InvariantSubspaceOnly) || wantBH;
164 //
165 // Set M to the dimension of the specified invariant subspace,
166 // and compute minimal size of WORK and IWORK
167 //
168 auto wsq = trsen_generic_wsq(job, select, T, m);
169 const IndexType workMin = wsq.first;
170 const IndexType iWorkMin = wsq.second;
171
172 IndexType lWork = work.length();
173 IndexType liWork = iwork.length();
174
175 if (lWork<workMin && lWork>0) {
176 ASSERT(0);
177 return -1;
178 } else if (liWork<iWorkMin && liWork>0) {
179 ASSERT(0);
180 return -2;
181 }
182
183 if (lWork==0) {
184 work.resize(workMin);
185 }
186 if (liWork==0) {
187 iwork.resize(iWorkMin);
188 }
189
190 work(1) = workMin;
191 iwork(1) = iWorkMin;
192
193 IndexType n1 = m;
194 IndexType n2 = n - m;
195 IndexType nn = n1*n2;
196
197 bool pair;
198 IndexType ks;
199 ElementType scale;
200
201 //
202 // Quick return if possible.
203 //
204 if (m==n || m==0) {
205 if (wantS) {
206 s = One;
207 }
208 if (wantSP) {
209 sep = lan(OneNorm, T);
210 }
211 goto DONE;
212 }
213 //
214 // Collect the selected blocks at the top-left corner of T.
215 //
216 ks = 0;
217 pair = false;
218 for (IndexType k=1; k<=n; ++k) {
219 if (pair) {
220 pair = false;
221 } else {
222 bool swap = select(k);
223 if (k<n) {
224 if (T(k+1,k)!=Zero) {
225 pair = true;
226 swap = swap || select(k+1);
227 }
228 }
229 if (swap) {
230 ++ks;
231 //
232 // Swap the K-th block to position KS.
233 //
234 IndexType iErr = 0;
235 IndexType kk = k;
236 if (k!=ks) {
237 trexc(computeQ, T, Q, kk, ks, work(_(1,n)));
238 }
239 if (iErr==1 || iErr==2) {
240 //
241 // Blocks too close to swap: exit.
242 //
243 info = 1;
244 if (wantS) {
245 s = Zero;
246 }
247 if (wantSP) {
248 sep = Zero;
249 }
250 goto DONE;
251 }
252 if (pair) {
253 ++ks;
254 }
255 }
256 }
257 }
258
259 if (wantS) {
260 GeMatrixView<ElementType> Work(n1, n2, work, n1);
261 //
262 // Solve Sylvester equation for R:
263 //
264 // T11*R - R*T22 = scale*T12
265 //
266 Work = T(_(1,n1),_(n1+1,n1+n2));
267 auto T11 = T(_(1,n1),_(1,n1));
268 auto T22 = T(_(n1+1,n1+n2),_(n1+1,n1+n2));
269
270 trsyl(NoTrans, NoTrans, -1, T11, T22, Work, scale);
271 //
272 // Estimate the reciprocal of the condition number of the cluster
273 // of eigenvalues.
274 //
275 ElementType rNorm = lan(FrobeniusNorm, Work);
276 if (rNorm==Zero) {
277 s = One;
278 } else {
279 s = scale / (sqrt(scale*scale/rNorm + rNorm)*sqrt(rNorm));
280 }
281 }
282
283 if (wantSP) {
284 //
285 // Estimate sep(T11,T22).
286 //
287 ElementType est = Zero;
288 IndexType kase = 0;
289
290 do {
291 auto _v = work(_(nn+1,nn+nn));
292 auto _x = work(_(1,nn));
293 auto _isgn = iwork(_(1,nn));
294
295 lacn2(_v, _x, _isgn, est, kase, isave);
296 if (kase==0) {
297 break;
298 }
299
300 GeMatrixView<ElementType> Work(n1, n2, work, n1);
301 auto T11 = T(_(1,n1),_(1,n1));
302 auto T22 = T(_(n1+1,n1+n2),_(n1+1,n1+n2));
303
304 if (kase==1) {
305 //
306 // Solve T11*R - R*T22 = scale*X.
307 //
308 trsyl(NoTrans, NoTrans, -1, T11, T22, Work, scale);
309 } else {
310 //
311 // Solve T11**T*R - R*T22**T = scale*X.
312 //
313 trsyl(Trans, Trans, -1, T11, T22, Work, scale);
314 }
315 } while (true);
316
317 sep = scale / est;
318 }
319
320 DONE:
321
322 //
323 // Store the output eigenvalues in WR and WI.
324 //
325 for (IndexType k=1; k<=n; ++k) {
326 wr(k) = T(k,k);
327 wi(k) = Zero;
328 }
329 for (IndexType k=1; k<=n-1; ++k) {
330 if (T(k+1,k)!=Zero) {
331 wi(k) = sqrt(abs(T(k,k+1)))*sqrt(abs(T(k+1,k)));
332 wi(k+1) = -wi(k);
333 }
334 }
335
336 work(1) = workMin;
337 iwork(1) = iWorkMin;
338 return info;
339 }
340
341 //== interface for native lapack ===============================================
342
343 #ifdef CHECK_CXXLAPACK
344
345 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
346 typename IndexType, typename S, typename SEP,
347 typename WORK, typename IWORK>
348 IndexType
349 trsen_native(TRSEN::Job job,
350 bool computeQ,
351 const DenseVector<SELECT> &select,
352 GeMatrix<MT> &T,
353 GeMatrix<MQ> &Q,
354 DenseVector<WR> &wr,
355 DenseVector<WI> &wi,
356 IndexType &m,
357 S &s,
358 SEP &sep,
359 DenseVector<WORK> &work,
360 DenseVector<IWORK> &iwork)
361 {
362 typedef typename GeMatrix<MT>::ElementType ElementType;
363
364 const char JOB = job;
365 const char COMPQ = (computeQ) ? 'V' : 'N';
366 const INTEGER N = T.numRows();
367 const INTEGER LDT = T.leadingDimension();
368 const INTEGER LDQ = Q.leadingDimension();
369 INTEGER _M = m;
370 DOUBLE _S = s;
371 DOUBLE _SEP = sep;
372 const INTEGER LWORK = work.length();
373 INTEGER LIWORK = iwork.length();
374 INTEGER INFO;
375
376 ASSERT(JOB=='N' || JOB=='E' || JOB=='V' || JOB=='B');
377
378 DenseVector<Array<LOGICAL> > _select(N);
379 for (IndexType i=1; i<=N; ++i) {
380 _select(i) = select(i);
381 }
382
383 DenseVector<Array<INTEGER> > _iwork(iwork.length());
384 for (IndexType i=1; i<=iwork.length(); ++i) {
385 _iwork(i) = iwork(i);
386 }
387
388 if (IsSame<ElementType,DOUBLE>::value) {
389 LAPACK_IMPL(dtrsen)(&JOB,
390 &COMPQ,
391 _select.data(),
392 &N,
393 T.data(),
394 &LDT,
395 Q.data(),
396 &LDQ,
397 wr.data(),
398 wi.data(),
399 &_M,
400 &_S,
401 &_SEP,
402 work.data(),
403 &LWORK,
404 _iwork.data(),
405 &LIWORK,
406 &INFO);
407 } else {
408 ASSERT(0);
409 }
410 if (INFO<0) {
411 std::cerr << "dtrsen: INFO = " << INFO << std::endl;
412 }
413 ASSERT(INFO>=0);
414
415 m = _M;
416 s = _S;
417 sep = _SEP;
418 for (IndexType i=1; i<=iwork.length(); ++i) {
419 iwork(i) = _iwork(i);
420 }
421
422 return INFO;
423 }
424
425 #endif // CHECK_CXXLAPACK
426
427 //== public interface ==========================================================
428
429 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
430 typename IndexType, typename S, typename SEP,
431 typename WORK, typename IWORK>
432 IndexType
433 trsen(TRSEN::Job job,
434 bool computeQ,
435 const DenseVector<SELECT> &select,
436 GeMatrix<MT> &T,
437 GeMatrix<MQ> &Q,
438 DenseVector<WR> &wr,
439 DenseVector<WI> &wi,
440 IndexType &m,
441 S &s,
442 SEP &sep,
443 DenseVector<WORK> &work,
444 DenseVector<IWORK> &iwork)
445 {
446 LAPACK_DEBUG_OUT("trsen");
447 //
448 // Test the input parameters
449 //
450 # ifndef NDEBUG
451 ASSERT(T.firstRow()==1);
452 ASSERT(T.firstCol()==1);
453 ASSERT(T.numRows()==T.numCols());
454
455 const IndexType n = T.numRows();
456 if (computeQ) {
457 ASSERT(Q.firstRow()==1);
458 ASSERT(Q.firstCol()==1);
459 ASSERT(Q.numRows()==n);
460 ASSERT(Q.numCols()==n);
461 }
462
463 ASSERT(wr.firstIndex()==1);
464 ASSERT(wr.length()==n);
465
466 ASSERT(wi.firstIndex()==1);
467 ASSERT(wi.length()==n);
468 # endif
469
470 # ifdef CHECK_CXXLAPACK
471 //
472 // Make copies of output arguments
473 //
474 const typename GeMatrix<MT>::NoView T_org = T;
475 const typename GeMatrix<MQ>::NoView Q_org = Q;
476 const typename DenseVector<WR>::NoView wr_org = wr;
477 const typename DenseVector<WI>::NoView wi_org = wi;
478 const IndexType m_org = m;
479 const S s_org = s;
480 const SEP sep_org = sep;
481 const typename DenseVector<WORK>::NoView work_org = work;
482 const typename DenseVector<IWORK>::NoView iwork_org = iwork;
483 # endif
484
485 //
486 // Call implementation
487 //
488 IndexType info = trsen_generic(job, computeQ, select, T, Q, wr, wi,
489 m, s, sep, work, iwork);
490
491 # ifdef CHECK_CXXLAPACK
492 //
493 // Make copies of results computed by generic implementation
494 //
495 const typename GeMatrix<MT>::NoView T_generic = T;
496 const typename GeMatrix<MQ>::NoView Q_generic = Q;
497 const typename DenseVector<WR>::NoView wr_generic = wr;
498 const typename DenseVector<WI>::NoView wi_generic = wi;
499 const IndexType m_generic = m;
500 const S s_generic = s;
501 const SEP sep_generic = sep;
502 const typename DenseVector<WORK>::NoView work_generic = work;
503 const typename DenseVector<IWORK>::NoView iwork_generic = iwork;
504
505 //
506 // restore output arguments
507 //
508 T = T_org;
509 Q = Q_org;
510 wr = wr_org;
511 wi = wi_org;
512 m = m_org;
513 s = s_org;
514 sep = sep_org;
515 work = work_org;
516 iwork = iwork_org;
517
518 //
519 // Compare generic results with results from the native implementation
520 //
521 IndexType _info = trsen_native(job, computeQ, select, T, Q, wr, wi,
522 m, s, sep, work, iwork);
523 bool failed = false;
524 if (! isIdentical(T_generic, T, "T_generic", "T")) {
525 std::cerr << "CXXLAPACK: T_generic = " << T_generic << std::endl;
526 std::cerr << "F77LAPACK: T = " << T << std::endl;
527 failed = true;
528 }
529 if (! isIdentical(Q_generic, Q, "Q_generic", "Q")) {
530 std::cerr << "CXXLAPACK: Q_generic = " << Q_generic << std::endl;
531 std::cerr << "F77LAPACK: Q = " << Q << std::endl;
532 failed = true;
533 }
534 if (! isIdentical(wr_generic, wr, "wr_generic", "wr")) {
535 std::cerr << "CXXLAPACK: wr_generic = " << wr_generic << std::endl;
536 std::cerr << "F77LAPACK: wr = " << wr << std::endl;
537 failed = true;
538 }
539 if (! isIdentical(wi_generic, wi, "wi_generic", "wi")) {
540 std::cerr << "CXXLAPACK: wi_generic = " << wi_generic << std::endl;
541 std::cerr << "F77LAPACK: wi = " << wi << std::endl;
542 failed = true;
543 }
544 if (! isIdentical(m_generic, m, "m_generic", "m")) {
545 std::cerr << "CXXLAPACK: m_generic = " << m_generic << std::endl;
546 std::cerr << "F77LAPACK: m = " << m << std::endl;
547 failed = true;
548 }
549 if (! isIdentical(s_generic, s, "s_generic", "s")) {
550 std::cerr << "CXXLAPACK: s_generic = " << s_generic << std::endl;
551 std::cerr << "F77LAPACK: s = " << s << std::endl;
552 failed = true;
553 }
554 if (! isIdentical(sep_generic, sep, "sep_generic", "sep")) {
555 std::cerr << "CXXLAPACK: sep_generic = " << sep_generic << std::endl;
556 std::cerr << "F77LAPACK: sep = " << sep << std::endl;
557 failed = true;
558 }
559 if (! isIdentical(work_generic, work, "work_generic", "work")) {
560 std::cerr << "CXXLAPACK: work_generic = " << work_generic << std::endl;
561 std::cerr << "F77LAPACK: work = " << work << std::endl;
562 failed = true;
563 }
564 if (! isIdentical(iwork_generic, iwork, "iwork_generic", "iwork")) {
565 std::cerr << "CXXLAPACK: iwork_generic = "
566 << iwork_generic << std::endl;
567 std::cerr << "F77LAPACK: iwork = " << iwork << std::endl;
568 failed = true;
569 }
570 if (! isIdentical(info, _info, "info", "_info")) {
571 std::cerr << "CXXLAPACK: info = " << info << std::endl;
572 std::cerr << "F77LAPACK: _info = " << _info << std::endl;
573 failed = true;
574 }
575
576 if (failed) {
577 std::cerr << "error in: trsen.tcc" << std::endl;
578 ASSERT(0);
579 } else {
580 // std::cerr << "passed: trsen.tcc" << std::endl;
581 }
582 # endif
583
584 return info;
585 }
586
587 //-- forwarding ----------------------------------------------------------------
588
589 template <typename SELECT, typename MT, typename MQ, typename WR, typename WI,
590 typename IndexType, typename S, typename SEP,
591 typename WORK, typename IWORK>
592 IndexType
593 trsen(TRSEN::Job job,
594 bool computeQ,
595 const DenseVector<SELECT> &select,
596 MT &&T,
597 MQ &&Q,
598 WR &&wr,
599 WI &&wi,
600 IndexType &&m,
601 S &&s,
602 SEP &&sep,
603 WORK &&work,
604 IWORK &&iwork)
605 {
606 trsen(job, computeQ, select, T, Q, wr, wi, m, s, sep, work, iwork);
607 }
608
609 } } // namespace lapack, flens
610
611 #endif // FLENS_LAPACK_EIG_TRSEN_TCC