gejsv.cc (flens/lapack/interface/src/gejsv.cc)

       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122

#define STR(x)      #x
#define STRING(x)   STR(x)

#include <flens/lapack/interface/include/config.h>

namespace flens { namespace lapack {

extern "C" {

//-- dgejsv --------------------------------------------------------------------
void
LAPACK_DECL(dgejsv)(const char       *JOBA,
                    const char       *JOBU,
                    const char       *JOBV,
                    const char       *JOBR,
                    const char       *JOBT,
                    const char       *JOBP,
                    const INTEGER    *M,
                    const INTEGER    *N,
                    DOUBLE           *A,
                    const INTEGER    *LDA,
                    DOUBLE           *SVA,
                    DOUBLE           *U,
                    const INTEGER    *LDU,
                    DOUBLE           *V,
                    const INTEGER    *LDV,
                    DOUBLE           *WORK,
                    const INTEGER    *LWORK,
                    INTEGER          *IWORK,
                    INTEGER          *INFO)
{
//
//  Test the input parameters so that we pass LAPACK error checks
//
    const bool lsvec  = (*JOBU=='U') || (*JOBU=='F');
    const bool jracc  = (*JOBV=='J');
    const bool rsvec  = (*JOBV=='V') || jracc;
    const bool rowpiv = (*JOBA=='F') || (*JOBA=='G');
    const bool l2rank = (*JOBA=='R');
    const bool l2aber = (*JOBA=='A');
    const bool errest = (*JOBA=='E') || (*JOBA=='G');
    const bool l2tran = (*JOBT=='T');
    const bool l2kill = (*JOBR=='R');
    const bool defr   = (*JOBR=='N');
    const bool l2pert = (*JOBP=='P');

    const int m = *M;
    const int n = *N;

    *INFO = 0;
    if (!(rowpiv || l2rank || l2aber || errest || *JOBA=='C')) {
        *INFO = -1;
    } else if (!(lsvec || *JOBU=='N' || *JOBU=='W')) {
        *INFO = -2;
    } else if (!(rsvec || *JOBV=='N' || *JOBV=='W') || (jracc && !lsvec)) {
        *INFO = -3;
    } else if (!(l2kill || defr)) {
        *INFO = -4;
    } else if (!(l2tran || *JOBT=='N')) {
        *INFO = -5;
    } else if (!(l2pert || *JOBP=='N')) {
        *INFO = -6;
    } else if (m<0) {
        *INFO = -7;
    } else if ((n<0) || (n>m)) {
        *INFO = -8;
    } else if (*LDA<m) {
        *INFO = -10;
    } else if (lsvec && *LDU<m) {
        *INFO = -13;
    } else if (rsvec && *LDV<n) {
        *INFO = -14;
    } else if ((!(lsvec || rsvec || errest) || (*LWORK < max(7,4*n+1,2*m+n)))
        || (!(lsvec || rsvec) || errest || (*LWORK < max(7,4*n+n*n,2*m+n)))
        || (lsvec || (!rsvec) || (*LWORK < max(7,2*m+n,4*n+1)))
        || (rsvec || (!lsvec) || (*LWORK < max(7,2*m+n,4*n+1)))
        || (lsvec || rsvec || (!jracc) || (*LWORK<max(2*m+n,6*n+2*n*n)))
        || (lsvec || rsvec || jracc || *LWORK<max(2*m+n,4*n+n*n,2*n+n*n+6)))
    {
        *INFO = -17;
    }

    if (*INFO!=0) {
        *INFO = -(*INFO);
        LAPACK_ERROR("DGEJSV", INFO);
        *INFO = -(*INFO);
        return;
    }

//
//  Call FLENS implementation
//
    JSV::Accuracy  accuracy    = JSV::Accuracy(*JOBA);
    JSV::JobU  jobU            = JSV::JobU(*JOBU);
    JSV::JobV  jobV            = JSV::JobV(*JOBV);
    const bool restrictedRange = (*JOBR=='R');
    const bool considerTransA  = (*JOBT=='T');
    const bool perturb         = (*JOBR=='P');

    DGeMatrixView     _A       = DFSView(m, n, A, *LDA);
    DDenseVectorView  _sva     = DArrayView(n, SVA, INTEGER(1));
    DGeMatrixView     _U       = DFSView(m, n, U, *LDU);
    DGeMatrixView     _V       = DFSView(n, n, V, *LDV);
    DDenseVectorView  _work    = DArrayView(*LWORK, WORK, INTEGER(1));
    IDenseVectorView  _iwork   = IArrayView(m+3*n, IWORK, INTEGER(1));

    *INFO = jsv(accuracy, jobU, jobV,
                restrictedRange, considerTransA, perturb,
                _A, _sva, _U, _V, _work, _iwork);

    if (*INFO<0) {
        *INFO = -(*INFO);
        LAPACK_ERROR("DGEJSV", INFO);
        *INFO = -(*INFO);
        return;
    }
}

} // extern "C"

} } // namespace lapack, flens