load "lapack"
real[int,int] tab(1,7);
tab=1;
tab(3:5,0)=0;
cout << " the single array in double array format tab = " << tab << endl;
matrix tab1;
tab1=tab;
cout << " the sparse vector tab = " << tab1 << endl;

int N=3,M=4;
real[int,int] A(N,M),B(M,N),C(N,M),AxB(N,N),AmC(N,M),BxA(M,M),AxBt(N,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;
C=-1;
C(0:1,2)=0;

//
// Matrix operation product
//
AxB=A*B;// works if lapack is loaded
BxA=B*A;
cout << " double array A (dim 3x4) = " << A << endl;
cout << " double array B (dim 4x3) = " << B << endl;
cout << " double array product A*B (dim 3x3)= " << AxB << endl;
cout << " double array product B*A (dim 4x4)= " << BxA << endl;
cout << " double array  A  " <<A<< endl;
//B=A';
//cout << " double array transpose A' only works the square case " <<B<< endl;
//cout <<"non square case the size of B is changed"<<endl;
matrix SA,SB,SAxSB,SAt,SC,SAmSC;
SA=A;
SB=B;
SC=C;
cout << " SA = " << SA << endl;
cout << " SB = " << SB << endl;
cout << " C  = " << C << endl;
cout << " SC = " << SC << endl;

SAmSC=SA+(-SC);
AmC=A-C;
SAxSB=SA*SB;
SAt=SA';

cout << " sparse matrix diff SA-SC = " << SAmSC << endl;
cout << " the same as matrix difference A-C = " << AmC << endl;
cout << " sparse matrix prod SA*SB = " << SAxSB << endl;
cout << " the same as matrix product A*B = " << AxB << endl;
cout << " sparse matrix transpose SA' = " << SAt << endl;

real[int,int] T(4,7),Tt(7,4);
T=10;
T(0,:)=1;T(:,3)=0;
cout << " the matrix T = " << T << endl;

matrix ST,STpSA,STt;
ST=T;
STt=ST';
cout << " the sparse matrix ST = " << ST << endl;
cout << " the sparse matrix STt = " << STt << endl;
cout << " +++++++++++++++++++++++++++++++++++++++++  " << endl;
cout << " adding sparse matrices of different sizes  " << endl;
cout << " takes largest of both sizes and add entries with common indices "<<endl;
cout << " +++++++++++++++++++++++++++++++++++++++++  " << endl;
cout << " the sparse matrix SA = " << SA << endl;
cout << " the sparse matrix ST = " << ST << endl;
STpSA=ST+SA;
cout << " the sparse matrix ST+SA = " << STpSA << endl;
STt=ST';
cout << " the sparse matrix ST = " << ST << endl;
cout << " the sparse matrix ST transpose IS WORKING in sparse matrix format= " << STt << endl;
real[int] v(7);
v=1;
real[int] test=ST*v;
cout << "  ST  = " << ST<< endl;
cout << "  v  = " << v<< endl;
cout << "4x7 sparse matrix times a 7x1 vector ST*v  = " << test << endl;