MTL 带状矩阵乘法

//MTL 带状矩阵乘法
//整理　by RobinKin
//
#include <iostream>
#include "mtl/matrix.h"
#include "mtl/mtl.h"
#include "mtl/utils.h"
#include "mtl/linalg_vec.h"
/*
Sample Output
Array A in packed form
[
[1,1,1,],
[2,2,2,],
[3,3,3,],
[4,4,],
[5,],
]
x
[1,2,3,4,5,]
y
[1,1,1,1,1,]
Ax + y
[15,39,75,75,53,]
*/
using namespace mtl;
typedef matrix< double, banded<>, banded<>, row_major>::type Matrix;
typedef dense1D<double> Vector;
int
main()
{
const int M = 5;
const int N = 5;

Matrix A(M, N, 0, 2);
Vector x(N);
Vector y(M);
// 1 1 1 1 1
//
// 2 2 2 2 1
//
// A = 3 3 3 x = 3 y = 1
//
// 4 4 4 1
//
// 5 5 1
int row = 1;
//赋值
Matrix::iterator ri = A.begin();
while (ri != A.end()) {
Matrix::Row::iterator i = (*ri).begin();
while (i != (*ri).end()) {
*i = row;
++i;
}
++row;
++ri;
}
for (int i = 0; i < N; ++i)
x[i] = i + 1;
mtl::set_value(y, 1);

std::cout << "Array A in packed form" << std::endl;
print_row(A);
std::cout << "x" << std::endl;
print_vector(x);
std::cout << "y" << std::endl;
print_vector(y);
//y=2*x*A+3*y
mult(A, scaled(x, 2), scaled(y, 3), y);
std::cout << "Ax + y" << std::endl;
print_vector(y);
return 0;
}