Home
Browse Files
  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 DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
 36      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
 37      $                   INFO )
 38  *
 39  *  -- LAPACK routine (version 3.3.1) --
 40  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 41  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 42  *  -- April 2011                                                      --
 43  */
 44 
 45 #ifndef FLENS_LAPACK_EIG_TRSNA_TCC
 46 #define FLENS_LAPACK_EIG_TRSNA_TCC 1
 47 
 48 #include <flens/aux/aux.h>
 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 VSELECT, typename MT, typename MVL, typename MVR,
 56           typename VS, typename VSEP, typename M, typename MM,
 57           typename MWORK, typename VIWORK>
 58 void
 59 trsna_generic(TRSNA::Job                    job,
 60               TRSNA::HowMany                howMany,
 61               const DenseVector<VSELECT>    &select,
 62               const GeMatrix<MT>            &T,
 63               const GeMatrix<MVL>           &_VL,
 64               const GeMatrix<MVR>           &_VR,
 65               DenseVector<VS>               &s,
 66               DenseVector<VSEP>             &sep,
 67               const MM                      &mm,
 68               M                             &m,
 69               GeMatrix<MWORK>               &Work,
 70               DenseVector<VIWORK>           &iWork)
 71 {
 72     using std::abs;
 73 
 74     typedef typename GeMatrix<MT>::ElementType  ElementType;
 75     typedef typename GeMatrix<MT>::IndexType    IndexType;
 76 
 77     const Underscore<IndexType> _;
 78 
 79     const IndexType n = T.numRows();
 80 
 81     const ElementType Zero(0), One(1), Two(2);
 82 //
 83 //  Local Arrays
 84 //
 85     IndexType iSaveData[3] = {0, 0, 0};
 86     DenseVectorView<IndexType>
 87         iSave = typename DenseVectorView<IndexType>::Engine(3, iSaveData, 1);
 88 //
 89 //  Decode and test the input parameters
 90 //
 91     const bool wantBH = (job==TRSNA::Both);
 92     const bool wantS = (job==TRSNA::EigenvaluesOnly) || wantBH;
 93     const bool wantSP = (job==TRSNA::EigenvectorsOnly) || wantBH;
 94 
 95     const bool someCon = (howMany==TRSNA::Selected);
 96 //
 97 //  Set M to the number of eigenpairs for which condition numbers
 98 //  are required, and test MM.
 99 //
100     if (someCon) {
101         m = 0;
102         bool pair = false;
103         for (IndexType k=1; k<=n; ++k) {
104             if (pair) {
105                 pair= false;
106             } else {
107                 if (k<n) {
108                     if (T(k+1,k)==Zero) {
109                         if (select(k)) {
110                             ++m;
111                         }
112                     } else {
113                         pair = true;
114                         if (select(k) || select(k+1)) {
115                             m += 2;
116                         }
117                     }
118                 } else {
119                     if (select(n)) {
120                         ++m;
121                     }
122                 }
123             }
124         }
125     } else {
126         m = n;
127     }
128 
129     if (mm<m) {
130         ASSERT(0);
131     }
132 
133     if (wantSP) {
134         ASSERT(_VL.numCols()>=m);
135         ASSERT(_VR.numCols()>=m);
136     }
137     // TODO: if one forgets to make this auto views const you get
138     //       some error that is hard to understand for newbies ...
139     //       Idea: disallow the creation of non-const views from
140     //       const matrices/vectors.
141     const auto VL = _VL(_,_(1,m));
142     const auto VR = _VR(_,_(1,m));
143 
144 //
145 //  Quick return if possible
146 //
147     if (n==0) {
148         return;
149     }
150 //
151     if (n==1) {
152         if (someCon) {
153             if (!select(1)) {
154                 return;
155             }
156         }
157         if (wantS) {
158             s(1) = One;
159         }
160         if (wantSP) {
161             sep(1) = abs(T(1,1));
162         }
163         return;
164     }
165 //
166 //  Get machine constants
167 //
168     const ElementType eps = lamch<ElementType>(Precision);
169     ElementType smallNum = lamch<ElementType>(SafeMin) / eps;
170     ElementType bigNum = One/smallNum;
171     labad(smallNum, bigNum);
172 
173     IndexType ks = 0;
174     bool  pair = false;
175     for (IndexType k=1; k<=n; ++k) {
176 //
177 //      Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
178 //
179         if (pair) {
180             pair = false;
181             continue;
182         } else {
183             if (k<n) {
184                 pair = (T(k+1,k)!=Zero);
185             }
186         }
187 //
188 //      Determine whether condition numbers are required for the k-th
189 //      eigenpair.
190 //
191         if (someCon) {
192             if (pair) {
193                 if (!select(k) && !select(k+1)) {
194                     continue;
195                 }
196             } else {
197                 if (!select(k)) {
198                     continue;
199                 }
200             }
201         }
202 
203         ++ks;
204 
205         if (wantS) {
206 //
207 //          Compute the reciprocal condition number of the k-th
208 //          eigenvalue.
209 //
210             if (!pair) {
211 //
212 //              Real eigenvalue.
213 //
214                 const ElementType prod = VR(_,ks) * VL(_,ks);
215                 const ElementType rNrm = blas::nrm2(VR(_,ks));
216                 const ElementType lNrm = blas::nrm2(VL(_,ks));
217                 s(ks) = abs(prod) / (rNrm*lNrm);
218             } else {
219 //
220 //              Complex eigenvalue.
221 //
222                 const ElementType prod1 = VR(_,ks) * VL(_,ks)
223                                         + VR(_,ks+1) * VL(_,ks+1);
224                 const ElementType prod2 = VL(_,ks) * VR(_,ks+1)
225                                         - VL(_,ks+1) * VR(_,ks);
226                 const ElementType rNrm = lapy2(blas::nrm2(VR(_,ks)),
227                                                blas::nrm2(VR(_,ks+1)));
228                 const ElementType lNrm = lapy2(blas::nrm2(VL(_,ks)),
229                                                blas::nrm2(VL(_,ks+1)));
230                 const ElementType cond = lapy2(prod1, prod2) / (rNrm*lNrm);
231                 s(ks)   = cond;
232                 s(ks+1) = cond;
233             }
234         }
235 
236         if (wantSP) {
237 //
238 //          Estimate the reciprocal condition number of the k-th
239 //          eigenvector.
240 //
241 //          Copy the matrix T to the array WORK and swap the diagonal
242 //          block beginning at T(k,k) to the (1,1) position.
243 //
244             auto _T = Work(_,_(1,n));
245             _T = T;
246             IndexType iFirst =k;
247             IndexType iLast = 1;
248             IndexType iErr = trexc(false, _T, _T, iFirst, iLast, Work(_,n+1));
249 
250             ElementType est, mu, scale;
251             IndexType   n2, nn;
252 
253             if (iErr==1 || iErr==2) {
254 //
255 //              Could not swap because blocks not well separated
256 //
257                 scale = One;
258                 est   = bigNum;
259             } else {
260 //
261 //              Reordering successful
262 //
263                 if (Work(2,1)==Zero) {
264 //
265 //                  Form C = T22 - lambda*I in WORK(2:N,2:N).
266 //
267                     for (IndexType i=2; i<=n; ++i) {
268                         Work(i,i) -= Work(1,1);
269                     }
270                     n2 = 1;
271                     nn = n - 1;
272                 } else {
273 //
274 //                  Triangularize the 2 by 2 block by unitary
275 //                  transformation U = [  cs   i*ss ]
276 //                                     [ i*ss   cs  ].
277 //                  such that the (1,1) position of WORK is complex
278 //                  eigenvalue lambda with positive imaginary part. (2,2)
279 //                  position of WORK is the complex eigenvalue lambda
280 //                  with negative imaginary  part.
281 //
282                     mu = sqrt(abs(Work(1,2))) * sqrt(abs(Work(2,1)));
283                     const ElementType delta = lapy2(mu, Work(2,1));
284                     const ElementType cs = mu / delta;
285                     const ElementType sn = -Work(2,1) / delta;
286 //
287 //                  Form
288 //
289 //                  C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
290 //                                           [   mu                     ]
291 //                                           [         ..               ]
292 //                                           [             ..           ]
293 //                                           [                  mu      ]
294 //                  where C**T is transpose of matrix C,
295 //                  and RWORK is stored starting in the N+1-st column of
296 //                  WORK.
297 //
298                     Work(2,_(3,n)) *= cs;
299                     for (IndexType j=3; j<=n; ++j) {
300                         Work(j,j) -= Work(1,1);
301                     }
302                     Work(2,2) = Zero;
303 
304                     Work(1,n+1) =Two*mu;
305                     for (IndexType i=2; i<=n-1; ++i) {
306                         Work(i,n+1) = sn*Work(1,i+1);
307                     }
308                     n2 = 2;
309                     nn = 2*(n-1);
310                 }
311 //
312 //              Estimate norm(inv(C**T))
313 //
314                 est = Zero;
315                 IndexType kase = 0;
316                 do {
317                     auto _v = Work(_,_(n+2,n+3)).vectorView(1,nn);
318                     auto _x = Work(_,_(n+4,n+5)).vectorView(1,nn);
319                     auto _iSgn = iWork(_(1,nn));
320                     lacn2(_v, _x, _iSgn, est, kase, iSave);
321                     if (kase==0) {
322                         break;
323                     } else {
324                         auto T = Work(_(2,n),_(2,n));
325                         auto b = Work(_(1,n-1),n+1);
326                         auto x = Work(_,_(n+4,n+5)).vectorView(1,2*(n-1));
327                         auto w = Work(_(1,n-1),n+6);
328                         ElementType dummyMu;
329                         if (kase==1) {
330                             if (n2==1) {
331 //
332 //                              Real eigenvalue: solve C**T*x = scale*c.
333 //
334                                 laqtr(true, true, T, b, dummyMu, scale, x, w);
335                             } else {
336 //
337 //                              Complex eigenvalue: solve
338 //                              C**T*(p+iq) = scale*(c+id) in real arithmetic.
339 //
340                                 laqtr(true, false, T, b, mu, scale, x, w);
341                             }
342                         } else {
343                             if (n2==1) {
344 //
345 //                              Real eigenvalue: solve C*x = scale*c.
346 //
347                                 laqtr(false, true, T, w, dummyMu, scale, x, w);
348                             } else {
349 //
350 //                              Complex eigenvalue: solve
351 //                              C*(p+iq) = scale*(c+id) in real arithmetic.
352 //
353                                 laqtr(false, false, T, b, mu, scale, x, w);
354 
355                             }
356                         }
357                     }
358                 } while (true);
359             }
360 
361             sep(ks) = scale / max(est, smallNum);
362             if (pair) {
363                 sep(ks+1) = sep(ks);
364             }
365         }
366 
367         if (pair) {
368             ++ks;
369         }
370     }
371 }
372 
373 //== interface for native lapack ===============================================
374 
375 #ifdef CHECK_CXXLAPACK
376 
377 template <typename VSELECT, typename MT, typename MVL, typename MVR,
378           typename VS, typename VSEP, typename M, typename MM,
379           typename MWORK, typename VIWORK>
380 void
381 trsna_native(TRSNA::Job                    job,
382              TRSNA::HowMany                howMany,
383              const DenseVector<VSELECT>    &select,
384              const GeMatrix<MT>            &T,
385              const GeMatrix<MVL>           &VL,
386              const GeMatrix<MVR>           &VR,
387              DenseVector<VS>               &s,
388              DenseVector<VSEP>             &sep,
389              const MM                      &mm,
390              M                             &m,
391              GeMatrix<MWORK>               &Work,
392              DenseVector<VIWORK>           &iWork)
393 {
394     typedef typename GeMatrix<MT>::ElementType   ElementType;
395     typedef typename GeMatrix<MT>::IndexType     IndexType;
396 
397     const char       JOB    = char(job);
398     const char       HOWMNY = char(howMany);
399     LOGICAL          *SELECT = 0;
400     const INTEGER    N = T.numRows();
401     const INTEGER    LDT = T.leadingDimension();
402     const INTEGER    LDVL = VL.leadingDimension();
403     const INTEGER    LDVR = VR.leadingDimension();
404     const INTEGER    _MM = mm;
405     INTEGER          _M = m;
406     INTEGER          LDWORK = Work.leadingDimension();
407     INTEGER          INFO;
408 
409     if (howMany==TRSNA::Selected) {
410         SELECT = new LOGICAL[N];
411         for (INTEGER i=1; i<=N; ++i) {
412             SELECT[i] = select(i);
413         }
414     }
415 
416     DenseVector<Array<INTEGER> >    _iWork(2*(N-1));
417     ASSERT(_iWork.length()>=iWork.length());
418     for (IndexType i=1; i<=iWork.length(); ++i) {
419         _iWork(i) = iWork(i);
420     }
421 
422     if (IsSame<ElementType,DOUBLE>::value) {
423         LAPACK_IMPL(dtrsna)(&JOB,
424                             &HOWMNY,
425                             SELECT,
426                             &N,
427                             T.data(),
428                             &LDT,
429                             VL.data(),
430                             &LDVL,
431                             VR.data(),
432                             &LDVR,
433                             s.data(),
434                             sep.data(),
435                             &_MM,
436                             &_M,
437                             Work.data(),
438                             &LDWORK,
439                             _iWork.data(),
440                             &INFO);
441     } else {
442         ASSERT(0);
443     }
444     ASSERT(INFO==0);
445 
446     if (howMany==TRSNA::Selected) {
447         ASSERT(SELECT);
448         delete [] SELECT;
449     }
450 
451     m = _M;
452 
453     for (IndexType i=1; i<=iWork.length(); ++i) {
454         iWork(i) = _iWork(i);
455     }
456 
457 }
458 
459 #endif // CHECK_CXXLAPACK
460 
461 //== public interface ==========================================================
462 
463 template <typename VSELECT, typename MT, typename MVL, typename MVR,
464           typename VS, typename VSEP, typename MM, typename M,
465           typename MWORK, typename VIWORK>
466 void
467 trsna(TRSNA::Job                    job,
468       TRSNA::HowMany                howMany,
469       const DenseVector<VSELECT>    &select,
470       const GeMatrix<MT>            &T,
471       const GeMatrix<MVL>           &VL,
472       const GeMatrix<MVR>           &VR,
473       DenseVector<VS>               &s,
474       DenseVector<VSEP>             &sep,
475       const MM                      &mm,
476       M                             &m,
477       GeMatrix<MWORK>               &Work,
478       DenseVector<VIWORK>           &iWork)
479 {
480     typedef typename GeMatrix<MT>::IndexType IndexType;
481 
482     const IndexType n = T.numRows();
483 
484 #   ifndef NDEBUG
485     ASSERT(T.firstRow()==1);
486     ASSERT(T.firstCol()==1);
487     ASSERT(T.numRows()==T.numCols());
488 
489     if (howMany!=TRSNA::All) {
490         ASSERT(select.firstIndex()==1);
491         ASSERT(select.length()==n);
492     }
493 
494     if (job!=TRSNA::EigenvectorsOnly) {
495         ASSERT(VL.firstRow()==1);
496         ASSERT(VL.firstCol()==1);
497         ASSERT(VL.numRows()==n);
498     }
499 
500     if (job!=TRSNA::EigenvectorsOnly) {
501         ASSERT(VR.firstRow()==1);
502         ASSERT(VR.firstCol()==1);
503         ASSERT(VR.numRows()==n);
504     }
505 
506     ASSERT(s.firstIndex()==1);
507     ASSERT(s.length()==mm);
508 
509     ASSERT(sep.firstIndex()==1);
510     ASSERT(sep.length()==mm);
511 
512     if (job!=TRSNA::EigenvaluesOnly) {
513         ASSERT(Work.firstRow()==1);
514         ASSERT(Work.firstCol()==1);
515         ASSERT(Work.numRows()==n);
516         ASSERT(Work.numCols()==n+6);
517 
518         ASSERT(iWork.firstIndex()==1);
519         ASSERT(iWork.length()==2*(n-1));
520     }
521 #   endif
522 
523 //
524 //  Make copies of output arguments
525 //
526 #   ifdef CHECK_CXXLAPACK
527     typename DenseVector<VS>::NoView      s_org     = s;
528     typename DenseVector<VSEP>::NoView    sep_org   = sep;
529     M                                     m_org     = m;
530     typename GeMatrix<MWORK>::NoView      Work_org  = Work;
531     typename DenseVector<VIWORK>::NoView  iWork_org = iWork;
532 #   endif
533 
534     trsna_generic(job, howMany, select, T, VL, VR, s, sep, mm, m, Work, iWork);
535 
536 #   ifdef CHECK_CXXLAPACK
537 //
538 //  Make copies of results computed by the generic implementation
539 //
540     typename DenseVector<VS>::NoView      s_generic     = s;
541     typename DenseVector<VSEP>::NoView    sep_generic   = sep;
542     M                                     m_generic     = m;
543     typename GeMatrix<MWORK>::NoView      Work_generic  = Work;
544     typename DenseVector<VIWORK>::NoView  iWork_generic = iWork;
545 
546 //
547 //  restore output arguments
548 //
549     s       = s_org;
550     sep     = sep_org;
551     m       = m_org;
552     Work    = Work_org;
553     iWork   = iWork_org;
554 
555     trsna_native(job, howMany, select, T, VL, VR, s, sep, mm, m, Work, iWork);
556 
557     bool failed = false;
558     if (! isIdentical(s_generic, s, "s_generic", "s")) {
559         std::cerr << "CXXLAPACK: s_generic = " << s_generic << std::endl;
560         std::cerr << "F77LAPACK: s = " << s << std::endl;
561         failed = true;
562     }
563 
564     if (! isIdentical(sep_generic, sep, "sep_generic", "sep")) {
565         std::cerr << "CXXLAPACK: sep_generic = " << sep_generic << std::endl;
566         std::cerr << "F77LAPACK: sep = " << sep << std::endl;
567         failed = true;
568     }
569 
570     if (! isIdentical(m_generic, m, "m_generic", "m")) {
571         std::cerr << "CXXLAPACK: m_generic = " << m_generic << std::endl;
572         std::cerr << "F77LAPACK: m = " << m << std::endl;
573         failed = true;
574     }
575 
576     if (! isIdentical(Work_generic, Work, "Work_generic", "Work")) {
577         std::cerr << "CXXLAPACK: Work_generic = " << Work_generic << std::endl;
578         std::cerr << "F77LAPACK: Work = " << Work << std::endl;
579         failed = true;
580     }
581 
582     if (! isIdentical(iWork_generic, iWork, "iWork_generic", "iWork")) {
583         std::cerr << "CXXLAPACK: iWork_generic = "
584                   << iWork_generic << std::endl;
585         std::cerr << "F77LAPACK: iWork = " << iWork << std::endl;
586         failed = true;
587     }
588 
589     if (failed) {
590         std::cerr << "n = " << n << std::endl;
591         ASSERT(0);
592     } else {
593 //        std::cerr << "passed: trsna.tcc" << std::endl;
594     }
595 #   endif
596 }
597 
598 //-- forwarding ----------------------------------------------------------------
599 template <typename VSELECT, typename MT, typename MVL, typename MVR,
600           typename VS, typename VSEP, typename MM, typename M,
601           typename MWORK, typename VIWORK>
602 void
603 trsna(TRSNA::Job                    job,
604       TRSNA::HowMany                howMany,
605       const VSELECT                 &select,
606       const MT                      &T,
607       const MVL                     &VL,
608       const MVR                     &VR,
609       VS                            &&s,
610       VSEP                          &&sep,
611       const MM                      &mm,
612       M                             &&m,
613       MWORK                         &&Work,
614       VIWORK                        &&iWork)
615 {
616     CHECKPOINT_ENTER;
617     trsna(job, howMany, select, T, VL, VR, s, sep, mm, m, Work, iWork);
618     CHECKPOINT_LEAVE;
619 }
620 
621 } } // namespace lapack, flens
622 
623 #endif // FLENS_LAPACK_EIG_TRSNA_TCC
Home