
//#define WANT_STREAM

#include "include.h"

#include "newmat.h"

#ifdef use_namespace
using namespace NEWMAT;
#endif


/**************************** test program ******************************/

void Print(const Matrix& X);
void Print(const UpperTriangularMatrix& X);
void Print(const DiagonalMatrix& X);
void Print(const SymmetricMatrix& X);
void Print(const LowerTriangularMatrix& X);
void Clean(Matrix&, Real);

void trymat3()
{
//   cout << "\nThird test of Matrix package\n";
   Tracer et("Third test of Matrix package");
   Tracer::PrintTrace(true);

   int i,j;

   {
      Tracer et1("Stage 1");
      SymmetricMatrix S(7);
      for (i=1;i<=7;i++) for (j=1;j<=i;j++) S(i,j)=i*i+j;
		S=-S+2.0;

      DiagonalMatrix D(7);
      for (i=1;i<=7;i++) D(i,i)=S(i,i);

      Matrix M4(7,7); { M4=D+(D+4.0); M4=M4-D*2.0;  M4=M4-4.0; Print(M4); }
      SymmetricMatrix S2=D; Matrix M2=S2;  { M2=-D+M2; Print(M2); }
      UpperTriangularMatrix U2=D; { M2=U2; M2=D-M2; Print(M2); }
      LowerTriangularMatrix L2=D; { M2=L2; M2=D-M2; Print(M2); }
      M2=D; M2=M2-D; Print(M2);
      for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i-j;
      U2=L2.t(); D=D.t(); S=S.t();
      M4=(L2-1.0)+(U2+1.0)-D-S; Print(M4);
      M4=(-L2+1.0)+(-U2-1.0)+D+S; Print(M4);
   }

   {
      Tracer et1("Stage 2");
      DiagonalMatrix D(6);
      for (i=1;i<=6;i++) D(i,i)=i*3.0+i*i+2.0;
      UpperTriangularMatrix U2(7); LowerTriangularMatrix L2(7);
      for (i=1;i<=7;i++) for (j=1;j<=i;j++) L2(i,j)=2.0-i*i+j;
		{ U2=L2.t(); }
      DiagonalMatrix D1(7); for (i=1;i<=7;i++) D1(i,i)=(i-2)*(i-4);
      Matrix M2(6,7);
      for (i=1;i<=6;i++) for (j=1;j<=7;j++) M2(i,j)=2.0+i*j+i*i+2.0*j*j;
      Matrix MD=D; SymmetricMatrix MD1(1); MD1=D1;
      Matrix MX=MD*M2*MD1 - D*(M2*D1); Print(MX);
      MX=MD*M2*MD1 - (D*M2)*D1; Print(MX);
      {
         D.ReDimension(7); for (i=1;i<=7;i++) D(i,i)=i*3.0+i*i+2.0;
         LowerTriangularMatrix LD(1); LD=D;
         UpperTriangularMatrix UD(1); UD=D;
         M2=U2; M2=LD*M2*MD1 - D*(U2*D1); Print(M2);
         M2=U2; M2=UD*M2*MD1 - (D*U2)*D1; Print(M2);
         M2=L2; M2=LD*M2*MD1 - D*(L2*D1); Print(M2);
         M2=L2; M2=UD*M2*MD1 - (D*L2)*D1; Print(M2);
      }
   }

   {
      Tracer et1("Stage 3");
      // test inverse * scalar
      DiagonalMatrix D(6);
      for (i=1;i<=6;i++) D(i)=i*i;
      DiagonalMatrix E = D.i() * 4.0;
      DiagonalMatrix I(6); I = 1.0;
      E=D*E-I*4.0; Print(E);
      E = D.i() / 0.25; E=D*E-I*4.0; Print(E);
   }
   {
      Tracer et1("Stage 4");
      Matrix sigma(3,3); Matrix sigmaI(3,3);
      sigma = 0; sigma(1,1) = 1.0; sigma(2,2) = 1.0; sigma(3,3) = 1.0;
      sigmaI = sigma.i();
      sigmaI -= sigma;  Clean(sigmaI, 0.000000001); Print(sigmaI);
   }
   {
      Tracer et1("Stage 5");
      Matrix X(5,5); DiagonalMatrix DM(5);
      for (int i=1; i<=5; i++) for (int j=1; j<=5; j++)
         X(i,j) = (23*i+59*j) % 43;
      DM << 1 << 8 << -7 << 2 << 3;
      Matrix Y = X.i() * DM; Y = X * Y - DM;
      Clean(Y, 0.000000001); Print(Y);
   }

//   cout << "\nEnd of third test\n";
}

