//---------------------------------------------------------------------------
// N.V.Shokhirev
// created: 20040726
// modified: 20041106
//---------------------------------------------------------------------------

#pragma hdrstop

#include "uLinAlgTests.h"
#include "dynarrays.h"
#include "DynLinAlg.h"
#include "MatUtils.h"
#include "dynarrutils.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

String LinAlgTest()
{
  AnsiString s = "", ss = "LinAlg Test";
  bool fail = false;
  int n = 5;

  FArr1D a(n);
  FArr1D b(n);
  FArr1D c(n);
  FArr1D d(n);
  FArr1D e(n);
  FArr2D q(n,n);

// #1 strid test
  RandArr(a,0.0,1.0);
  RandArr(b,0.0,1.0);
  RandArr(d,3.0,4.0);
  RandArr(c,-1.0,1.0);
  TridTo2D(b, d, a, q);
  e = c;
  c = MxV(q, e);
  strid(b, d, a, c);
  if ( ! Equal(c, e) ) { fail = true;   s += " strid |"; };

// #2 GaussD test
  FArr2D m(n,n+1);
  RandArr(q,-1.0,1.0);
  ShiftDiag(q, 4.0);
  c = MxV(q, a);  // a is from the test #1
  for(int i =1; i<=n; i++)
  {
    for(int j =1; j<=n; j++) m(i,j) = q(i,j);
    m(i,n+1) = c(i);
  };
  fail = !GaussD(m);
  if (fail) s += " GaussD - singular |"; ;
  for(int i =1; i<=n;i++) c(i) = m(i,n+1);
  if ( ! Equal(c, a) ) { fail = true;   s += " GaussD |"; };

// #3 GaussP test
  c = MxV(q, a);  // a is from the test #1, q - from #2
  for(int i =1; i<=n; i++)
  {
    for(int j =1; j<=n; j++) m(i,j) = q(i,j);
    m(i,n+1) = c(i);
  };
  fail = !GaussP(m);
  if (fail) s += " GaussP - singular |";
  for(int i =1; i<=n;i++) c(i) = m(i,n+1);
  if ( ! Equal(c, a) ) { fail = true;   s += " GaussP |"; };

// #4 JacobiIter test
  b = MxV(q, a);  // a is from the test #1, q - from #2
  RandArr(d,-0.1,0.1);
  for(int i =1; i<=n;i++) c(i) = a(i) + d(i);
  fail = !JacobiIter(q, c, b,1.0e-6);
  if (fail) s += " JacobiIter diverged |";
  if ( ! Equal(c, a, 1.0e-6) ) { fail = true;   s += " JacobiIter |"; };

  if (fail) ss += " FAIL:" + s;
  else ss += " - OK";
  return ss;
};

String CGTest()
{
  AnsiString s = "", ss = "Conjugate Gradients Test";
  bool fail = false;
  int n = 5;

  FArr2D a(n,n);
  FArr1D b(n);
  FArr1D x(n);
  FArr1D e(n);

  RandArr(a,1.0,2.0, 1);
  ShiftDiag(a, 4.0);
  RandArr(e,-1.0,1.0);
  b = MxV(a, e);
  fail = !ConjGrad(a, b, x, 1.0e-8);
  if (fail) {s += " ConjGrad: Iter > IterMax |"; };
  if ( ! Equal(x, e,1.1e-4) ) { fail = true;   s += " ConjGrad |"; };

  int m = 6;
  FArr2D c(m,n);
  FArr1D bb(m);
//   C(m,n), b(m), x(n)
  RandArr(c,1.0,2.0);
  RandArr(e,-1.0,1.0);
  bb = MxV(c, e);
  x.fill(0.0);
  fail = !ConjGradLS(c, bb, x, 1.0e-8);
  if (fail) {s += " ConjGradLS: Iter > IterMax |"; };
  if ( ! Equal(x, e,1.1e-4) ) { fail = true;   s += " ConjGradLS |"; };

  if (fail) ss += " FAIL:" + s;
  else ss += " - OK";
  return ss;
};