//$$ newmat1.cxx   Utilities and matrix type functions

// Copyright (C) 1991: R B Davies and DSIR

#define WANT_STREAM                     // to include stream package

#include "include.hxx"

#include "newmat.hxx"


/**************************** error handlers *******************************/

void MatrixError(char* message)
{
   cerr << "\nError detected by matrix package: " << message << "\n";
   exit(0);
}

void MatrixErrorNoSpace(void* v)
// give error message if v is null
{
   if (v) return;
   cerr << "\nError detected by matrix package: no space\n";
   exit(0);
}

/************************* ExeCounter functions *****************************/


// the next statment may need to be deleted form some compilers

int ExeCounter::nreports = 0;

ExeCounter::ExeCounter(int xl, int xf) : line(xl), fileid(xf), nexe(0) {}

ExeCounter::~ExeCounter()
{
   nreports++;
   cout << nreports << "  " << fileid << "  " << line << "  " << nexe << "\n";
}
 

/************************* MatrixType functions *****************************/


BOOL MatrixType::operator==(const MatrixType& t) const
{ return (type == t.type); }

BOOL MatrixType::operator!=(const MatrixType& t) const
{ return (type != t.type); }

MatrixType MatrixType::operator+(const MatrixType& mt) const
{
   Type type1=mt.type;
   if (type==UnSp || type1==UnSp) return UnSp;
   if (type==type1) return type;   // won't work for orthogonal
   if (type==ColV || type1==ColV) return ColV;
   if (type==RowV || type1==RowV) return RowV;
   if (type==EqEl && (type1==Sym || type1==Diag)) return Sym;
   if (type1==EqEl && (type==Sym || type==Diag)) return Sym;
   if (type==EqEl || type1==EqEl) return Rect;
   if (type==Diag) return type1;   // won't work for assym
   if (type1==Diag) return type;
   return Rect;
}

MatrixType MatrixType::operator-(const MatrixType& mt) const
{
   return *this+mt;        // won't work for orthogonal or pos def
}

MatrixType MatrixType::operator*(const MatrixType& mt) const
{
   Type type1=mt.type;
   if (type==UnSp || type1==UnSp) return UnSp;
   if (type1==ColV) { if (type==RowV) return Diag; else return ColV; }
   if (type==RowV) return RowV;
   if (type==Sym || type1==Sym) return Rect;
   if (type==type1) return type;    // won't work for pos def or assym
   if (type==EqEl || type1==EqEl) return Rect;
   if (type==Diag) return type1;
   if (type1==Diag) return type;
   return Rect;
}

BOOL MatrixType::operator>=(const MatrixType& mt) const
{
   Type type1=mt.type;
   if (type==type1) return TRUE;    // by definition
   if (type==Rect || type==RowV || type==ColV) return TRUE;
   if (type==EqEl) return FALSE;
   if (type1==Diag) return TRUE;    // won't work for pos def or assym
   return FALSE;
}

MatrixType MatrixType::i() const
{
   switch (type)
   {
   case UnSp:  return UnSp;
   case UT:    return UT;
   case LT:    return LT;
   case Rect:  return Rect;
   case Sym:   return Sym;
   case Diag:  return Diag;
   case RowV:  return Diag;
   case ColV:  return Diag;
   case EqEl:  return Diag;
   case Crout: return UnSp;
   }
}

MatrixType MatrixType::t() const
{
   switch (type)
   {
   case UnSp:  return UnSp;
   case UT:    return LT;
   case LT:    return UT;
   case Rect:  return Rect;
   case Sym:   return Sym;
   case Diag:  return Diag;
   case RowV:  return ColV;
   case ColV:  return RowV;
   case EqEl:  return EqEl;
   case Crout: return UnSp;
   }
}

MatrixType MatrixType::sub() const
{
   switch (type)
   {
   case UnSp:  return UnSp;
   case UT:    return Rect;
   case LT:    return Rect;
   case Rect:  return Rect;
   case Sym:   return Rect;
   case Diag:  return Rect;
   case RowV:  return RowV;
   case ColV:  return ColV;
   case EqEl:  return EqEl;
   case Crout: return UnSp;
   }
}

MatrixType MatrixType::ssub() const
{
   return *this;                  // won't work for selection matrix
}

MatrixType::operator char*() const
{
   switch (type)
   {
   case UnSp:  return "UnSp ";
   case UT:    return "UT   ";
   case LT:    return "LT   ";
   case Rect:  return "Rect ";
   case Sym:   return "Sym  ";
   case Diag:  return "Diag ";
   case RowV:  return "RowV ";
   case ColV:  return "ColV ";
   case EqEl:  return "EqEl ";
   case Crout: return "Crout";
   }
}

GeneralMatrix* MatrixType::New() const
{
   GeneralMatrix* gm;
   switch (type)
   {
   case UnSp:  MatrixError("Can't build matrix type UnSp");
   case UT:    gm = new UpperTriangularMatrix; break;
   case LT:    gm = new LowerTriangularMatrix; break;
   case Rect:  gm = new Matrix; break;
   case Sym:   gm = new SymmetricMatrix; break;
   case Diag:  gm = new DiagonalMatrix; break;
   case RowV:  gm = new RowVector; break;
   case ColV:  gm = new ColumnVector; break;
   case EqEl:  MatrixError("Can't build matrix type EqEl");
   case Crout: MatrixError("Can't build matrix type Crout");
   }
   MatrixErrorNoSpace(gm);
   gm->Protect(); return gm;
}

GeneralMatrix* MatrixType::New(int nr, int nc) const
{
   GeneralMatrix* gm;
   switch (type)
   {
   case UnSp:  MatrixError("Can't build matrix type UnSp");
   case UT:    gm = new UpperTriangularMatrix(nr); break;
   case LT:    gm = new LowerTriangularMatrix(nr); break;
   case Rect:  gm = new Matrix(nr, nc); break;
   case Sym:   gm = new SymmetricMatrix(nr); break;
   case Diag:  gm = new DiagonalMatrix(nr); break;
   case RowV:  if (nr!=1) MatrixError("Illegal Row Vector");
	       gm = new RowVector(nc); break;
   case ColV:  if (nc!=1) MatrixError("Illegal Column Vector");
	       gm = new ColumnVector(nr); break;
   case EqEl:  MatrixError("Can't build matrix type EqEl");
   case Crout: MatrixError("Can't build matrix type Crout");
   }
   MatrixErrorNoSpace(gm);
   gm->Protect(); return gm;
}

void TestTypeAdd()
{
   MatrixType list[] = { MatrixType::UT,
                         MatrixType::LT,
                         MatrixType::Rect,
                         MatrixType::Sym,
                         MatrixType::Diag,
                         MatrixType::RowV,
                         MatrixType::ColV,
			 MatrixType::EqEl };

   cout << "+     ";
   for (int i=0; i<MatrixType::nTypes(); i++) cout << (char*)list[i] << " ";
   cout << "\n";
   for (i=0; i<MatrixType::nTypes(); i++)
   {
      cout << (char*)(list[i]) << " ";
      for (int j=0; j<MatrixType::nTypes(); j++)
	 cout << (char*)(list[j]+list[i]) << " ";
      cout << "\n";
   }
   cout << "\n";
}

void TestTypeMult()
{
   MatrixType list[] = { MatrixType::UT,
                         MatrixType::LT,
                         MatrixType::Rect,
                         MatrixType::Sym,
                         MatrixType::Diag,
                         MatrixType::RowV,
                         MatrixType::ColV,
			 MatrixType::EqEl };
   cout << "*     ";
   for (int i=0; i<MatrixType::nTypes(); i++)
      cout << (char*)list[i] << " ";
   cout << "\n";
   for (i=0; i<MatrixType::nTypes(); i++)
   {
      cout << (char*)list[i] << " ";
      for (int j=0; j<MatrixType::nTypes(); j++)
	 cout << (char*)(list[j]*list[i]) << " ";
      cout << "\n";
   }
   cout << "\n";
}

void TestTypeOrder()
{
   MatrixType list[] = { MatrixType::UT,
                         MatrixType::LT,
                         MatrixType::Rect,
                         MatrixType::Sym,
                         MatrixType::Diag,
                         MatrixType::RowV,
                         MatrixType::ColV,
			 MatrixType::EqEl };
   cout << ">=    ";
   for (int i = 0; i<MatrixType::nTypes(); i++)
      cout << (char*)list[i] << " ";
   cout << "\n";
   for (i=0; i<MatrixType::nTypes(); i++)
   {
      cout << (char*)list[i] << " ";
      for (int j=0; j<MatrixType::nTypes(); j++)
	 cout << ((list[j]>=list[i]) ? "Yes   " : "No    ");
      cout << "\n";
   }
   cout << "\n";
}
