Blitz++: Power iteration

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

[code lang="cpp"]

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*1; 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; }

[/code]

*1:unsigned int)time(0