load "lapack"
int N=3,M=3;
real[int,int] A(N,M),B(M,N),AxB(N,N),BxA(M,M),AT(M,N); 
A=1; 
A(0,0)=0;
A(0,2)=0;
A(1,0)=0;
A(1,1)=0;
B=3; 
B(0,1)=0;
B(0,2)=0;
B(1,0)=0;
B(1,0)=0;
//
// Matrix operations  product, transpose
//
AxB=A*B;// works if lapack is loaded
AT=A';
cout << " double array A = " << A << endl;
cout << " double array B = " << B << endl;
cout << " double array product A*B works if lapack is loaded = " << AxB << endl;
cout << " double array transpose A'  " <<AT<< endl;


real[int,int] P(3,3),C(3,4),D(3,4),E(4,4),Ct(4,3),T(4,3),Tt(3,4);
real[int] a(3),b(4),d(4);
a=[1,2,3];
b=[4,5,6,7];
C=1;
C(2:3,1)=0;
cout << " double array C = " << C << endl;
cout << " first row of C = " << C(0,:) << endl;
d=C(0,:);
cout << " first row of C = " << d <<"size "<<d.n<< endl;
cout << " single array dim 3x1 displayed as a row 1x3 "<<endl;
cout<<" a = " << a <<"size "<<a.n<< endl;
cout << " single array dim 4x1 displayed as a row 1x4 "<<endl;
cout<<" b = " << b <<"size "<<b.n<< endl;
D=a*b';
cout << " D=a*b' = " << D << endl;
a=C*b;
cout << " (dim 3x4)*dim(4x1)=(dim 3x1) C*b  size = " << a << endl;
b=D*d;
cout << "  (dim 3x4)*dim(4x1)=(dim 3x1) D*d  size = " << b << endl;
P=1;
cout << "   P  size = " << P << endl;
cout << "   C  size = " << C << endl;
D=P*C;
cout << "  (dim 3x3)*dim(3x4)=(dim 3x4) P*C  size = " << D << endl;
T=1;
Tt=T';
cout << "  T  size = " << T << endl;
cout << "  T'  size = " << Tt << endl;

//E=Ct*D;