//整理 by RobinKin
//Blitz++ 张量计算的示例
/*****************************************************************************
* 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 would 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);//指标收缩, 计算矩阵积
// The above expression computes the matrix product of A and B.
cout << "D = " << D << endl;
// D = 4 x 4
// 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
// 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,r,r); // A four-dimensional array
fourthIndex l; // Placeholder for the fourth index
E = A(l,j) * B(k,i);//指标轮换
//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;
}
//输出
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 = 4 x 4
[ 0 0 0 1
0 0 1 0
0 1 0 0
1 0 0 0 ]
B = 4 x 4
[ 0 0 0 1
0 0 1 0
0 1 0 0
1 0 0 0 ]
D = 4 x 4
[ 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 ]
