Home
Browse Files
  1 #include <flens/lapack/interface/include/config.h>
  2 
  3 
  4 namespace flens { namespace lapack {
  5 
  6 template <typename T>
  7 struct SelectFunction
  8 {
  9     typedef LOGICAL (* Function)(const T *, const T *);
 10 
 11     SelectFunction(Function _select)
 12         : select(_select)
 13     {
 14     }
 15 
 16     bool
 17     operator()(const T &a, const T &b)
 18     {
 19         return select(&a, &b);
 20     }
 21 
 22     Function  select;
 23 };
 24 
 25 
 26 extern "C" {
 27 
 28 //-- dgeesx --------------------------------------------------------------------
 29 void
 30 LAPACK_DECL(dgeesx)(const char       *JOBVS,
 31                     const char       *SORT,
 32                     LOGICAL          (*SELECT)(const DOUBLE *, const DOUBLE *),
 33                     const char       *SENSE,
 34                     const INTEGER    *N,
 35                     DOUBLE           *A,
 36                     const INTEGER    *LDA,
 37                     INTEGER          *SDIM,
 38                     DOUBLE           *WR,
 39                     DOUBLE           *WI,
 40                     DOUBLE           *VS,
 41                     const INTEGER    *LDVS,
 42                     DOUBLE           *RCONDE,
 43                     DOUBLE           *RCONDV,
 44                     DOUBLE           *WORK,
 45                     const INTEGER    *LWORK,
 46                     INTEGER          *IWORK,
 47                     const INTEGER    *LIWORK,
 48                     LOGICAL          *BWORK,
 49                     INTEGER          *INFO)
 50 {
 51     DEBUG_FLENS_LAPACK("dgeesx");
 52 
 53     using std::max;
 54     using std::min;
 55 //
 56 //  Test the input parameters so that we pass LAPACK error checks
 57 //
 58     *INFO = 0;
 59     const bool lQuery = (*LWORK==-1) || (*LIWORK==-1);
 60     const bool wantVS = (*JOBVS=='V');
 61     const bool wantST = (*SORT=='S');
 62 
 63     const bool wantSN = (*SENSE=='N');
 64     const bool wantSE = (*SENSE=='E');
 65     const bool wantSV = (*SENSE=='V');
 66     const bool wantSB = (*SENSE=='B');
 67 
 68     if ((!wantVS) && (*JOBVS!='N')) {
 69         *INFO = 1;
 70     } else if ((!wantST) && (*SORT!='N')) {
 71         *INFO = 2;
 72     } else if (!(wantSN || wantSE || wantSV || wantSB)
 73             || ( !wantST && !wantSN))
 74     {
 75         *INFO = 4;
 76     } else if (*N<0) {
 77         *INFO = 5;
 78     } else if (*LDA<max(INTEGER(1),*N)) {
 79         *INFO = 7;
 80     } else if (*LDVS<1 || (wantVS && *LDVS<*N)) {
 81         *INFO = 12;
 82     }
 83 
 84     if (*INFO!=0) {
 85         LAPACK_ERROR("DGEESX", INFO);
 86         *INFO = -(*INFO);
 87         return;
 88     }
 89 
 90 //
 91 //  Setup FLENS matrix/vector types
 92 //
 93     typedef DGeMatrixView::IndexType IndexType;
 94 
 95     SENSE::Sense sense = getFlensLapackEnum<SENSE::Sense>(*SENSE);
 96 
 97     DGeMatrixView       _A      = DFSView(*N, *N, A, *LDA);
 98     IndexType           _SDIM   = *SDIM;
 99     DDenseVectorView    _WR     = DArrayView(*N, WR, 1);
100     DDenseVectorView    _WI     = DArrayView(*N, WI, 1);
101     DGeMatrixView       _VS     = DFSView(*N, *N, VS, *LDVS);
102     DDenseVectorView    _WORK   = DArrayView(*LWORK, WORK, 1);
103 
104     IDenseVector        _IWORK(*LIWORK);
105     for (INTEGER i=1; i<*LIWORK; ++i) {
106         _IWORK(i) = IWORK[i-1];
107     }
108 
109     BDenseVector        _BWORK(*N);
110     for (INTEGER i=1; i<*N; ++i) {
111         _BWORK(i) = BWORK[i-1];
112     }
113 
114 //
115 //  Test if work has at least minimal worksize
116 //
117     auto wsQuery = esx_wsq(wantVS, sense, _A);
118     INTEGER minWork = wsQuery.first;
119 
120     if (*LWORK<minWork && !lQuery) {
121         *INFO = 16;
122     } else if (*LIWORK<1 && !lQuery) {
123         *INFO = 18;
124     }
125 
126     if (*INFO!=0) {
127         LAPACK_ERROR("DGEESX", INFO);
128         *INFO = -(*INFO);
129         return;
130     }
131 //
132 //  Call FLENS implementation
133 //
134     SelectFunction<DOUBLE> select(SELECT);
135 
136     esx(wantVS, wantST, select, sense, _A, _SDIM, _WR, _WI, _VS,
137         *RCONDE, *RCONDV, _WORK, _IWORK, _BWORK);
138 
139     *SDIM   = _SDIM;
140     for (INTEGER i=1; i<*LIWORK; ++i) {
141         IWORK[i-1] = _IWORK(i);
142     }
143     for (INTEGER i=1; i<*N; ++i) {
144         BWORK[i-1] = _BWORK(i);
145     }
146 }
147 
148 } // extern "C"
149 
150 } } // namespace lapack, flens
Home