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