#include <boost/numeric/ublas/matrix.hpp>

#include <boost/numeric/ublas/io.hpp>

using namespace boost::numeric::ublas;

int main(int argc, char *argv[])

{

matrix<matrix<double> > AA(3,3), BB(3,3); // or matrix<matrix<double> > AA(3,3,0), BB(3,3,0);

matrix<double> CC(3,3);

for (unsigned int i=0; i!=3; i++)

for (unsigned int j=0; j!=3; j++)

CC(i,j) = i + j + 1 + i*j;

std::cout << CC << std::endl;

AA(0,0)=CC;

AA(1,1)=CC;

AA(2,2)=CC;

std::cout << AA << std::endl;

for (unsigned int i=0; i!=3; i++)

for (unsigned int j=0; j!=3; j++)

CC(i,j) = 2*i + 0.1*j + 0.75 + i*j*i*0.22;

BB(0,0)=CC;

BB(1,1)=CC;

BB(2,2)=CC;

std::cout << BB << std::endl;

std::cout << AA+BB << std::endl;

return 0;

}