0

Blitz++: Power iteration

Tensor Power iteration using Blitz++. The following code is too dirty.

#include <iostream>

#include <blitz/array.h>
#include <random/uniform.h>
#include <time.h>
#include <ctime>
#include <algorithm>
#ifdef BZ_HAVE_STD
#include <fstream>
#else
#include <fstream.h>
#endif

blitz::firstIndex i1;
blitz::secondIndex i2;
blitz::thirdIndex i3;
blitz::fourthIndex i4;
blitz::fifthIndex i5;
blitz::sixthIndex i6;
blitz::seventhIndex i7;
blitz::eighthIndex i8;
blitz::ninthIndex i9;
blitz::tenthIndex i10;
blitz::eleventhIndex i11;

template <int index>
blitz::Array<float,1> TVMIteration(
    blitz::Array<float,index>& T,
    blitz::Array<float,1>& v
)
{
    blitz::Array<float,1> Tv(v.size());
    if(index==11)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11)*v(i11), i11)*v(i10), i10)*v(i9), i9)*v(i8), i8)*v(i7), i7)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==10)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10)*v(i10), i10)*v(i9), i9)*v(i8), i8)*v(i7), i7)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==9)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6,i7,i8,i9)*v(i9), i9)*v(i8), i8)*v(i7), i7)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==8)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6,i7,i8)*v(i8), i8)*v(i7), i7)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==7)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6,i7)*v(i7), i7)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==6)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5,i6)*v(i6), i6)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==5)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4,i5)*v(i5), i5)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==4)
    {
        Tv = blitz::sum(blitz::sum(blitz::sum(T(i1,i2,i3,i4)*v(i4), i4)*v(i3), i3)*v(i2), i2);
    }
    else if(index==3)
    {
        Tv = blitz::sum(blitz::sum(T(i1,i2,i3)*v(i3), i3)*v(i2), i2);
    }
    else if(index==2)
    {
        Tv = blitz::sum(T(i1,i2)*v(i2), i2);
    }
    return Tv;
}

template <int rankT>
void PM(
    blitz::Array<float,1>& x
)
{
    std::cout << "O(T) = " << rankT << std::endl;
    int size = x.size();
    blitz::Array<float, rankT> T;
    if(rankT == 2)
    {
        T.resize(size,size);
        T = x(i1) * x(i2);
    }
    else if(rankT == 3)
    {
        T.resize(size,size,size);
        T = x(i1) * x(i2) * x(i3);
    }
    else if(rankT == 4)
    {
        T.resize(size,size,size,size);
        T = x(i1) * x(i2) * x(i3) * x(i4);
    }
    else if(rankT == 5)
    {
        T.resize(size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5);
    }
    else if(rankT == 6)
    {
        T.resize(size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6);
    }
    else if(rankT == 7)
    {
        T.resize(size,size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6)*x(i7);
    }
    else if(rankT == 8)
    {
        T.resize(size,size,size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6)*x(i7)*x(i8);
    }
    else if(rankT == 9)
    {
        T.resize(size,size,size,size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6)*x(i7)*x(i8)*x(i9);
    }
    else if(rankT == 10)
    {
        T.resize(size,size,size,size,size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6)*x(i7)*x(i8)*x(i9)*x(i10);
    }
    else if(rankT == 11)
    {
        T.resize(size,size,size,size,size,size,size,size,size,size,size);
        T = x(i1)*x(i2)*x(i3)*x(i4)*x(i5)*x(i6)*x(i7)*x(i8)*x(i9)*x(i10)*x(i11);
    }

    blitz::Array<float,1> v(size), vTmp(size), vInit(size);
    ranlib::Uniform<float> r;
    r.seed((unsigned int)time(0));
    for(int n = 0; n < size; ++n)
    {
        vInit(n) = r.random();
    }

    clock_t t0, t1;
    float vDiff;
    int n;

    /// TVMIteration
    t0 = clock();
    vDiff = 1000.0;
    n = 0;
    v = vInit;
    blitz::Array<float,1> T_(size);
    while((n<1000) && (vDiff > std::numeric_limits<float>::epsilon()))
    {
        T_ = TVMIteration<rankT>(T,v);
        vTmp = T_;
        vTmp /= sum(vTmp);
        vDiff = sum(v-vTmp);
        v = vTmp;
        std::cout << "  ite(" << ++n-1 << "), diff = " << vDiff << std::endl;
    }
    t1 = clock();
    std::cout << "  TVMIteration: " << (double)(t1-t0) / CLOCKS_PER_SEC << " sec." << std::endl;

//    std::cout << "x = " << x << std::endl;
//    std::cout << "v = " << v << std::endl;
}

int main()
{
    int sizeVec = 60;
    blitz::Array<float,1> x(sizeVec);
    ranlib::Uniform<float> r;
    r.seed((unsigned int)time(0));
    for(int n = 0; n < sizeVec; ++n)
    {
        x(n) = r.random();
    }
    x /= sum(x);
    std::cout << "x = " << x << std::endl;

    PM<2>(x);
    PM<3>(x);
    PM<4>(x);
    PM<5>(x);
//    PM<6>(x);
    return 0;
}

Advertisements