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
#include <cassert>
#include <iostream>
#include <utility>
#include <flens/flens.cxx>

using namespace flens;
using namespace std;

template <typename MA>
std::pair<long, long>
eigenvalues_wsq(const GeMatrix<MA> &A)
{
    // Compute and return minimal and optimal workspace sizes
    return std::pair<long,long>(A.numRows()*4, A.numRows()*256);
}

template <typename MA, typename VV, typename VWS>
void
eigenvalues(const GeMatrix<MA> &A, DenseVector<VV> &v, DenseVector<VWS> &ws)
{
    auto wsq = eigenvalues_wsq(A);
    if (ws.length()==0) {
        ws.resize(wsq.second);
    }
    assert(ws.length()>=wsq.first);

    //
    // ... compute eigenvalues v of A ... here the actual work happens!!!
    //
}

template <typename MA, typename VV>
void
eigenvalues(const GeMatrix<MA> &A, DenseVector<VV> &v)
{
    DenseVector<VV> ws;

    eigenvalues(A, v, ws);
}

int
main()
{
    typedef GeMatrix<FullStorage<double> >  DGeMatrix;
    typedef DenseVector<Array<double> >     DDenseVector;

    // Use case 1: Use optimal workspace (performance for memory).  We have
    //             a typical case where the same problem size is involved.  So
    //             the workspace can be reused.
    {
        DGeMatrix     A(100, 100);
        DDenseVector  v(100);

        auto workspace = eigenvalues_wsq(A);
        DDenseVector  ws(workspace.second);

        for (int i=0; i<666; ++i) {
            fillRandom(A);
            eigenvalues(A, v, ws);
        }
    }

    // Use case 2: Use minimal workspace (memory for performance)
    {
        DGeMatrix     A(100, 100);
        DDenseVector  v(100);

        auto workspace = eigenvalues_wsq(A);
        DDenseVector  ws(workspace.first);

        for (int i=0; i<666; ++i) {
            fillRandom(A);
            eigenvalues(A, v, ws);
        }
    }

    // Use case 3: Allocate (optimal) workspace on first call.
    {
        DGeMatrix     A(100, 100);
        DDenseVector  v(100);
        DDenseVector  ws;

        for (int i=0; i<666; ++i) {
            fillRandom(A);
            eigenvalues(A, v, ws);
        }
    }

    // NO-Use case 4: Workspace gets allocated for each call!!!
    {
        DGeMatrix     A(100, 100);
        DDenseVector  v(100);

        for (int i=0; i<666; ++i) {
            fillRandom(A);
            eigenvalues(A, v);
        }
    }

    // Use case 5: Ok, workspace is needed only once.
    {
        DGeMatrix     A(100, 100);
        DDenseVector  v(100);

        fillRandom(A);
        eigenvalues(A, v);
    }
}





































































































Home