/* ----------------------------------------------------------------------------

realmat: realmat.cc

realmat is a C++ matrix class. The header file realmat.h describes its use.

Copyright (C) 1990.

A. Ronald Gallant
P.O. Box 5513 
Raleigh NC 27650-5513 
USA   

Permission to use, copy, modify, and distribute this software and its 
documentation for any purpose and without fee is hereby granted, provided 
that the above copyright notice appear in all copies and that both that 
copyright notice and this permission notice appear in supporting 
documentation.  

This software is provided "as is" without any expressed or implied warranty.

---------------------------------------------------------------------------- */
#include "realmat.h"

realmat::realmat(INTEGER r, INTEGER c, REAL* a)
{
  if (r<=0) error("Error, realmat::realmat, Number of rows not positive");
  if (c<=0) error("Error, realmat::realmat, Number of columns not positive");
  rows=r; 
  cols=c; 
  len=r*c; 
  x = a;
}


void realmat::realmat_constructor(INTEGER r, INTEGER c, REAL fill_value)
{
  if (r<=0) error("Error, realmat::realmat, Number of rows not positive");
  if (c<=0) error("Error, realmat::realmat, Number of columns not positive");
  rows=r; 
  cols=c; 
  len=r*c; 
  x = new REAL[len];
  if (x == 0) error("Error, realmat::realmat, Operator new failed");
  REAL* top = &(x[len-1]);
  REAL* t = x;
  while (t <= top) *t++ = fill_value;
}

realmat::realmat(realmat& a)
{
    rows=a.rows;
    cols=a.cols;
    if(rows<=0) 
      error("Error, realmat::realmat, Number of rows not positive");
    if(cols<=0) 
      error("Error, realmat::realmat, Number of columns not positive");
    len=a.rows*a.cols; 
    x = new REAL[len];
    if (x == 0) error("Error, realmat::realmat, Operator new failed");
    REAL* top = &(x[len-1]);
    REAL* t = x;
    REAL* u = a.x;
    while (t <= top) *t++ = *u++;
}

void realmat::resize_constructor(INTEGER r, INTEGER c, REAL fill_value)
{
  if (r<=0) error("Error, realmat::resize, Number of rows not positive");
  if (c<=0) error("Error, realmat::resize, Number of columns not positive");
  delete x;
  rows=r; 
  cols=c; 
  len=r*c; 
  x = new REAL[len];
  if (x == 0) error("Error, realmat::resize, Operator new failed");
  REAL* top = &(x[len-1]);
  REAL* t = x;
  while (t <= top) *t++ = fill_value;
}

realmat& realmat::operator=(realmat& a)
{
  if (this != &a) {
    delete x;
    rows=a.rows;
    cols=a.cols;
    len=a.len;
    x = new REAL[len];
    if (x == 0) error("Error, realmat::operator=, Operator new failed");
    REAL* top = &(x[len-1]);
    REAL* t = x;
    REAL* u = a.x;
    while (t <= top) *t++ = *u++;
  }
  return *this;
}

REAL& realmat::check1(INTEGER i)
{
  if ((1<=i) && (i<=len))
    return x[i-1];
  else
    error("Error, realmat::check1, Index out of range");
    exit(1); //This keeps the compiler from complaining;
}

REAL& realmat::check2(INTEGER i, INTEGER j)
{
  if ((1<=i) && (1<=j) && (i<=rows) && (j<=cols))
    return x[i + rows*j - rows - 1];  //  return x[rows*(j-1)+i-1]
  else
    error("Error, realmat::check2, Index out of range");
    exit(1); //This keeps the compiler from complaining;
}

realmat T(realmat& a)
{
  INTEGER newrows = a.cols;
  INTEGER newcols = a.rows;
  INTEGER newlen = newrows*newcols;
  REAL*   newx = new REAL[newlen];
  if (newx == 0) a.error("Error, realmat, T, Operator new failed");
  INTEGER i,j;
  for (j = 0; j < a.cols; j++)
  for (i = 0; i < a.rows; i++)
    newx[j + newrows*i] = a.x[i + a.rows*j];
  realmat d(newrows,newcols,newx);
  return d;
}

realmat operator+(realmat& a, realmat& b)
{
  if ((a.cols == b.cols) && (a.rows == b.rows)) {
    REAL* newx = new REAL[a.len];
    if (newx == 0) a.error("Error, realmat, operator+, Operator new failed");
    REAL* top = &(newx[a.len-1]);
    REAL* t = newx;
    REAL* u = a.x;
    REAL* v = b.x;
    while (t <= top) *t++ = *u++ + *v++;
    realmat d(a.rows,a.cols,newx);
    return d;
  }
  else
    a.error("Error, realmat, operator+, Matrices not conformable."); 
    exit(1); //This keeps the compiler from complaining;
}

realmat operator+(realmat& a)
{
  REAL* newx = new REAL[a.len];
  if (newx == 0) a.error("Error, realmat, operator+, Operator new failed");
  REAL* top = &(newx[a.len-1]);
  REAL* t = newx;
  REAL* u = a.x;
  while (t <= top) *t++ = *u++;
  realmat d(a.rows,a.cols,newx);
  return d;
}

realmat operator-(realmat& a, realmat& b)
{
  if ((a.cols == b.cols) && (a.rows == b.rows)) {
    REAL* newx = new REAL[a.len];
    if (newx == 0) a.error("Error, realmat, operator-, Operator new failed");
    REAL* top = &(newx[a.len-1]);
    REAL* t = newx;
    REAL* u = a.x;
    REAL* v = b.x;
    while (t <= top) *t++ = *u++ - *v++;
    realmat d(a.rows,a.cols,newx);
    return d;
  }
  else
    a.error("Error, realmat, operator-, Matrices not conformable."); 
    exit(1); //This keeps the compiler from complaining;
}

realmat operator-(realmat& a)
{
  REAL* newx = new REAL[a.len];
  if (newx == 0) a.error("Error, realmat, operator-, Operator new failed");
  REAL* top = &(newx[a.len-1]);
  REAL* t = newx;
  REAL* u = a.x;
  while (t <= top) *t++ = - *u++;
  realmat d(a.rows,a.cols,newx);
  return d;
}

realmat operator*(realmat& a, realmat& b)
{
  if (a.cols == b.rows) {
    INTEGER newlen = a.rows*b.cols;
    REAL*   newx = new REAL[newlen];
    if (newx == 0) a.error("Error, realmat, operator*, Operator new failed");
    REAL    zero = 0.0;
    REAL*   top = &(newx[newlen-1]);
    REAL*   t = newx;
    while (t <= top) *t++ = zero;
    INTEGER i,j,k;
    for (j = 0; j < b.cols; j++)
    for (k = 0; k < a.cols; k++)
    for (i = 0; i < a.rows; i++)
      newx[i+a.rows*j] += a.x[i+a.rows*k] * b.x[k+b.rows*j];
    realmat d(a.rows,b.cols,newx);
    return d;
  }
  else {
    a.error("Error, realmat, operator*, Matrices not conformable."); 
    exit(1); //this keeps compiler from complaining
  }
}

realmat operator*(REAL& a, realmat& b)
{
  REAL*   newx = new REAL[b.len];
  if (newx == 0) b.error("Error, realmat, operator*, Operator new failed");
  REAL* top = &(newx[b.len-1]);
  REAL* t = newx;
  REAL* u = b.x;
  while (t <= top) *t++ = a * *u++;
  realmat d(b.rows,b.cols,newx);
  return d;
}

realmat operator*(INTEGER& a, realmat& b)
{
  REAL* newx = new REAL[b.len];
  if (newx == 0) b.error("Error, realmat, operator*, Operator new failed");
  REAL  f = a;
  REAL* top = &(newx[b.len-1]);
  REAL* t = newx;
  REAL* u = b.x;
  while (t <= top) *t++ = f * *u++;
  realmat d(b.rows,b.cols,newx);
  return d;
}

void default_realmat_error_handler(const char* msg)
{
  cerr << msg << "\n";
  exit(1);
}

ONE_ARG_ERROR_HANDLER_T realmat_error_handler = default_realmat_error_handler;

ONE_ARG_ERROR_HANDLER_T set_realmat_error_handler(ONE_ARG_ERROR_HANDLER_T f)
{
  ONE_ARG_ERROR_HANDLER_T old = realmat_error_handler;
  realmat_error_handler = f;
  return old;
}

void realmat::error(const char* msg)
{
  (*realmat_error_handler)(msg);
}

ostream& operator<<(ostream& stream, realmat& a)
{
  char    line[256], eol[2] = {'\n','\0'}, bl[3] = {' ','\n','\0'};
  char    *next,*save,*fcode;
  REAL    f,af;
  INTEGER i,j,linesize,maxcol,pad,start,stop,r,c;

  r=rows(a);
  c=cols(a);

  ((LINESIZE<72)||(LINESIZE>133)) ? (linesize=133) : (linesize=LINESIZE);
  maxcol = (linesize - 8)/12;
  if (c < maxcol) maxcol = c;
  pad = (linesize - 8 - 12*maxcol)/2 + 1;
  memset(line,' ',pad);
  save = line + pad;

  start=1;
  do {
    stop = start - 1 + maxcol;
    if (stop > c) stop = c;

    stream << bl;
    stream << bl;

    next = save;
    memset(next,' ',6);
    next = next + 6;

    for (j = start; j <= stop; j++) {
      if (j < 1000)      {sprintf(next,"      Col%3i",j); next = next+12;}
      else if (j <32760) {sprintf(next,"      C%5i"  ,j); next = next+12;}
      else               {sprintf(next,"      TooBig"  ); next = next+12;}
    }                                   

    memcpy(next,eol,2);
                          
    stream << line;
    stream << bl;

    for (i = 1; i <= r; i++) {
      next = save;
      if (i < 1000)      {sprintf(next,"Row%3i",i); next = next+6;}
      else if (i <32760) {sprintf(next,"R%5i"  ,i); next = next+6;}
      else               {sprintf(next,"TooBig"  ); next = next+6;}
      for (j = start; j <= stop; j++) {
        fcode = "%12.3e";
        f = a.elem(i,j);
        af = f; if (f < 0.E0)  af = -f ;
        if (af < 1.E+8 ) fcode = "%12.0f";
        if (af < 1.E+5 ) fcode = "%12.1f";
        if (af < 1.E+4 ) fcode = "%12.2f";
        if (af < 1.E+3 ) fcode = "%12.3f";
        if (af < 1.E+2 ) fcode = "%12.4f";
        if (af < 1.E+1 ) fcode = "%12.5f";
        if (af < 1.E+0 ) fcode = "%12.6f";
        if (af < 1.E-1 ) fcode = "%12.7f";
        if (af < 1.E-2 ) fcode = "%12.8f";
        if (af < 1.E-4 ) fcode = "%12.3e";
        if (af < 1.E-30) fcode = "%12.1f";
        sprintf(next,fcode,f); next = next+12;
      }
      memcpy(next,eol,2);
      stream << line;
    }
    start = stop + 1;
  } while (stop < c);

  return stream;
}

realmat invpsd(realmat& a, REAL eps)
{
  if (a.cols != a.rows) 
    a.error("Error, realmat, invpsd, Matrix not square.");
  INTEGER i;
  for (i=0; i<a.rows; i++)
    if (a.x[i+a.rows*i] < 0) 
      a.error("Error, realmat, invpsd, Matrix not positive semi-definite.");

  INTEGER newrows = a.cols;
  INTEGER newcols = a.rows;
  INTEGER newlen = newrows*newcols;
  REAL*   newx = new REAL[newlen];
  if (newx == 0) a.error("Error, realmat, invpsd, Operator new failed");

  REAL* top = &(newx[a.len-1]);
  REAL* t = newx;
  REAL* u = a.x;
    while (t <= top) *t++ = *u++;

  REAL* s = new REAL[newrows];
  if (s == 0) a.error("Error, realmat, invpsd, Operator new failed");

  dcond(newx,newrows,s,0);
  INTEGER ier = dsweep(newx,newrows,eps);
  dcond(newx,newrows,s,1);

  ier++;  //This stops a compiler warning.

  delete s;

  realmat d(newrows,newcols,newx);
  return d;
}