cout << " ---------------------------------------------------------------  " << endl;
cout << " ----------------------ARRAY+MATRIX TESTS ----------------------  " << endl;
cout << " ---------------------------------------------------------------  " << endl;
real[int,int] R(4,4);
R=0;// All entries of R set to zero
// Set nonzero some entries
R(0,0)=2;R(0,1)=-3;R(1,0)=7;R(1,1)=2;R(1,2)=11;R(2,1)= 5;R(2,2)=2;R(3,2)=2;R(1,3)=2;

cout << " double array R. Size  " << R << endl;
// Set R sparse matrix format
matrix A=R;
cout << "R as sparse matrix and called A " <<  endl;
cout << "............................... " <<  endl;
cout << A << endl;
cout << "............................... " <<  endl;
int[int] IA,JA;
real[int] ZA;
[IA,JA,ZA]=A ;
cout <<" A=[IA,JA,ZA] and then "<<endl;
cout <<" IA non-zero row indices from sparse matrix  A " <<  endl;
cout << "............................... " <<  endl;
cout << IA << endl;
cout << "............................... " <<  endl;
cout << "JA non-zero column  indices from sparse matrix  A " <<  endl;
cout << "............................... " <<  endl;
cout << JA << endl;
cout << "............................... " <<  endl;
cout << "ZA non-zero entries from sparse matrix  A " <<  endl;
cout << "............................... " <<  endl;
cout << ZA << endl;
cout << "...................................................................... " <<  endl;
cout << "..Observe the last entry in IA, JA and ZA to save the correct size !!! " <<  endl;
cout << "...................................................................... " <<  endl;
cout << " Recover A from IA,JA,ZA via  B=[IA,JA,ZA]; " <<  endl;
matrix B=[IA,JA,ZA];
cout << "sparse matrix B " <<  endl;
cout << "............................... " <<  endl;
cout << B << endl;
cout << "............................... " <<  endl;
cout << "size of A " <<A.n<<  endl;
int[int] IB(A.n);
//
cout<<" Construct an index vector IB to remove one column or row of A "<<endl;
cout<<" Vector IB has the size (number of rows or columns) of A "<<endl;
cout<<" and -1 at the row or column we want to remove "<<endl;
cout<<" The rest are renumbered from 0 to remaining rows (or columns) "<<endl;
cout<<" For instance, size of A is "<<A.n<<endl;
//
IB(0)=0;
IB(1)=-1;
IB(2)=1;
IB(3)=2;

cout << "... IB index vector to remove second column and row from A ... " <<  endl;
cout<<IB<<endl;

matrix D;
D=A(IB^-1,IB^-1);// To use only entries not in row or column 2
cout << "sparse matrix A without row and column 2 " <<  endl;
cout << "............................... " <<  endl;
cout << D << endl;
cout << "............................... " <<  endl;

int[int] ID,JD;
real[int] ZD;
[ID,JD,ZD]=D ;
real[int,int] DD(D.n,D.m);
DD=0;
for (int i=0;i<ID.n;i++)
DD(ID(i),JD(i))=D(ID(i),JD(i));
cout << " double array D is R without row and column 2  " << endl;
cout << " double array D  size " << DD << endl;
cout << " Recall that double array R was: size " << R << endl;