Home Browse Files
       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
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
#include <cxxstd/iostream.h>
///
///  With header __flens.cxx__ all of FLENS gets included.
///
///  :links:  __flens.cxx__ -> file:flens/flens.cxx
#include <flens/flens.cxx>

using namespace std;
using namespace flens;

typedef double   T;

//
//  Traits for creating one based matrix/vector views
//
namespace flens {

template <typename A>
struct OneBased
{
};

template <typename MA>
struct OneBased<GeMatrix<MA> >
{
    typedef typename MA::ElementType      T;
    typedef typename MA::IndexType        IndexType;
    typedef IndexOptions<IndexType, 1>    IndexBase;

    static const StorageOrder order = MA::order;

    typedef GeMatrix<FullStorageView<T, order, IndexBase> >         View;
    typedef GeMatrix<ConstFullStorageView<T, order, IndexBase> >    ConstView;
};

template <typename VX>
struct OneBased<DenseVector<VX> >
{
    typedef typename VX::ElementType      T;
    typedef typename VX::IndexType        IndexType;
    typedef IndexOptions<IndexType, 1>    IndexBase;

    typedef DenseVector<ArrayView<T, IndexBase> >         View;
    typedef DenseVector<ConstArrayView<T, IndexBase> >    ConstView;
};

} // namespace flens

//
//  LAPACK wrapper for non-one based indices
//
namespace flens { namespace mylapack {

template <typename MA, typename VPIV>
typename RestrictTo<IsGeMatrix<MA>::value
                 && IsIntegerDenseVector<VPIV>::value,
         typename RemoveRef<MA>::Type::IndexType>::Type
trf(MA &&A, VPIV &&piv)
{
///
/// Remove references from rvalue types
///
    typedef typename RemoveRef<MA>::Type    MatrixA;
    typedef typename MatrixA::IndexType     IndexType;
    typedef typename RemoveRef<VPIV>::Type  VectorPiv;

///
/// Create views of the arguments
///
    typename OneBased<MatrixA>::View    A_    = A;
    typename OneBased<VectorPiv>::View  piv_  = piv;

    A_.changeIndexBase(1,1);
    piv_.changeIndexBase(1);

///
/// Make the views one-based
///
    IndexType info = lapack::trf(A_, piv_);

    const IndexType diff = piv.firstIndex() - piv_.firstIndex();

    for (IndexType i=1; i<=piv_.length(); ++i) {
        piv_(i) += diff;
    }

    return info;
}

} } // namespace mylapack, flens

int
main()
{
    ///
    ///  Define an index-option type for zero based indices
    ///
    typedef IndexOptions<int, 0>                           ZeroBased;

    ///
    ///  Define some convenient typedefs for the matrix type
    ///  of our system of linear equations.
    ///
    typedef GeMatrix<FullStorage<T, ColMajor, ZeroBased> > Matrix;

    ///
    ///  We also need an extra vector type for the pivots.  The type of the
    ///  pivots is taken for the system matrix.
    ///
    typedef Matrix::IndexType                              IndexType;
    typedef DenseVector<Array<IndexType, ZeroBased> >      IndexVector;

    ///
    ///  Set up the baby problem ...
    ///
    const IndexType m = 4,
                    n = 4;

    ///
    /// Zero based matrix and vector
    ///
    Matrix            Ab(m, n);
    IndexVector       piv(m);

    Ab = 2,  4,  4,  4,
         4,  2,  4,  4,
         4,  4,  2,  4,
         4,  4,  4,  2;

    cout << "Ab = " << Ab << endl;

    ///
    /// Compute the $LU$ factorization with __lapack::trf__
    ///
    mylapack::trf(Ab, piv);

    cout << "Ab.firstRow() = " << Ab.firstRow() << endl;
    cout << "Ab.firstCol() = " << Ab.firstCol() << endl;

    cout << "Ab = " << Ab << endl;
    cout << "piv = " << piv << endl;
}





































































































Home