Blitz ++ matrix multiplication (tensor operation) Example

// Sort by mongokin

// Blitz ++ tensor calculation example

/***************************************************************************** * matmult.cpp     Blitz++ tensor notation example ***************************************************************************** * This example illustrates the tensor-like notation provided by Blitz++. */#include <blitz/array.h>#include <iostream>using namespace blitz;int main(){    // Create two 4x4 arrays.  We want them to look like matrices, so    // we'll make the valid index range 1..4 (rather than 0..3 which is    // the default).    Range r(1,4);    Array<float,2> A(r,r), B(r,r);    // The first will be a Hilbert matrix:    //    // a   =   1    //  ij   -----    //       i+j-1    //    // Blitz++ provides a set of types { firstIndex, secondIndex, ... }    // which act as placeholders for indices.  These can be used directly    // in expressions.  For example, we can fill out the A matrix like this:    firstIndex i;    // Placeholder for the first index    secondIndex j;   // Placeholder for the second index    A = 1.0 / (i+j-1);    cout << "A = " << A << endl;    // A = 4 x 4    //         1       0.5  0.333333      0.25    //       0.5  0.333333      0.25       0.2    //  0.333333      0.25       0.2  0.166667    //      0.25       0.2  0.166667  0.142857    // Now the A matrix has each element equal to a_ij = 1/(i+j-1).    // The matrix B will be the permutation matrix    //    // [ 0 0 0 1 ]    // [ 0 0 1 0 ]    // [ 0 1 0 0 ]    // [ 1 0 0 0 ]    //    // Here are two ways of filling out B:    B = (i == (5-j));         // Using an equation -- a bit cryptic    cout << "B = " << B << endl;    // B = 4 x 4    //         0         0         0         1    //         0         0         1         0    //         0         1         0         0    //         1         0         0         0    B = 0, 0, 0, 1,           // Using an initializer list        0, 0, 1, 0,                   0, 1, 0, 0,        1, 0, 0, 0;    cout << "B = " << B << endl;    // Now some examples of tensor-like notation.    Array<float,3> C(r,r,r);  // A three-dimensional array: 1..4, 1..4, 1..4    thirdIndex k;             // Placeholder for the third index    // This expression will set    //    // c    = a   * b    //  ijk    ik    kj //   C = A(i,k) * B(k,j);

Cout <"c =" <C <Endl; // in real tensor notation, the repeated K index wocould imply a // contraction (or summation) along K. in blitz ++, you must explicitly // indicate contractions using the sum (expr, index) function: array <float, 2> D (R, R ); D = sum (a (I, k) * B (K, J), k); // indicator shrinkage, calculate the Matrix Product // The above expression computes the Matrix Product of A and B. cout <"d =" <D <Endl; // d = 4x4 // 0.25 0.333333 0.5 1 // 0.2 0.25 0.333333 0.5 // 0.166667 0.2 0.25 0.333333 0.142857 // indices like I, j, K can be used in any order in an expression. // For example, the following computes a Kronecker Product of A and B, // But permutes the indices along the way: array <float, 4> E (R, r); // a four-dimen1_array fourthindex L; // placeholder for the fourth index E = a (L, j) * B (K, I); // metric Rotation

//cout << "E = " << E << endl;    // Now let's fill out a two-dimensional array with a radially symmetric    // decaying sinusoid.    int N = 64;                   // Size of array: N x N    Array<float,2> F(N,N);    float midpoint = (N-1)/2.;    int cycles = 3;    float omega = 2.0 * M_PI * cycles / double(N);    float tau = - 10.0 / N;    F = cos(omega * sqrt(pow2(i-midpoint) + pow2(j-midpoint)))        * exp(tau * sqrt(pow2(i-midpoint) + pow2(j-midpoint)));

//cout << "F = " << F << endl;    return 0;}

// Output

A = 4 x 4
[         1       0.5  0.333333      0.25 
        0.5  0.333333      0.25       0.2 
   0.333333      0.25       0.2  0.166667 
       0.25       0.2  0.166667  0.142857 ]

B = 4x4
[0 0 0 1
0 0 1 0
0 1 0 0
1 0 0 0]

B = 4x4
[0 0 0 1
0 0 1 0
0 1 0 0
1 0 0 0]

D = 4x4
[0.25 0.333333 0.5 1
0.2 0.25 0.333333 0.5
0.166667 0.2 0.25 0.333333
0.142857 0.166667 0.2 0.25]