newmat3.cpp

matrix library for linux and windos

         // while (i--) { *RowCopy++ = *Mstore; Mstore += ++row; }
         if (i) for (;;)
            { *RowCopy++ = *Mstore; if (!(--i)) break; Mstore += ++row; }
      }
   }
}

void SymmetricMatrix::GetCol(MatrixRowCol& mrc)
{
   // do not allow StoreHere
   if (+(mrc.cw*StoreHere))
      Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));

   int col=mrc.rowcol; mrc.length=nrows;
   REPORT
   mrc.skip=0;
   if (+(mrc.cw*DirectPart))    // actually get row ??
      { REPORT mrc.storage=col+1; mrc.data=store+(col*(col+1))/2; }
   else
   {
      // do not allow StoreOnExit and !DirectPart
      if (+(mrc.cw*StoreOnExit))
         Throw(InternalException("SymmetricMatrix::GetCol(MatrixRowCol&)"));

      mrc.storage=ncols; Real* ColCopy;
      if ( +(mrc.cw*HaveStore)) { REPORT ColCopy = mrc.data; }
      else
      {
         REPORT                                      // not accessed
         ColCopy = new Real [ncols]; MatrixErrorNoSpace(ColCopy);
         MONITOR_REAL_NEW("Make (SymGetCol)",ncols,ColCopy)
         mrc.cw += HaveStore; mrc.data = ColCopy;
      }
      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* Mstore = store+(col*(col+1))/2; int i = col;
         while (i--) *ColCopy++ = *Mstore++;
         i = ncols-col;
         // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
         if (i) for (;;)
            { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
      }
   }
}

void SymmetricMatrix::GetCol(MatrixColX& mrc)
{
   int col=mrc.rowcol; mrc.length=nrows;
   if (+(mrc.cw*DirectPart))
   {
      REPORT
      mrc.skip=col; int i=nrows-col; mrc.storage=i;
      mrc.data = mrc.store+col;
      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT                           // not accessed
         Real* ColCopy = mrc.data;
         Real* Mstore = store+(col*(col+3))/2;
         // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
         if (i) for (;;)
            { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
      }
   }
   else
   {
      REPORT
      // do not allow StoreOnExit and !DirectPart
      if (+(mrc.cw*StoreOnExit))
         Throw(InternalException("SymmetricMatrix::GetCol(MatrixColX&)"));

      mrc.skip=0; mrc.storage=ncols;
      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* ColCopy = mrc.data;
         Real* Mstore = store+(col*(col+1))/2; int i = col;
         while (i--) *ColCopy++ = *Mstore++;
         i = ncols-col;
         // while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
         if (i) for (;;)
            { *ColCopy++ = *Mstore; if (!(--i)) break; Mstore += ++col; }
      }
   }
}

// Do not need RestoreRow because we do not allow !DirectPart && StoreOnExit

void SymmetricMatrix::RestoreCol(MatrixColX& mrc)
{
   REPORT
   // Really do restore column
   int col=mrc.rowcol; Real* Cstore = mrc.data;
   Real* Mstore = store+(col*(col+3))/2; int i = nrows-col;
   // while (i--) { *Mstore = *Cstore++; Mstore+= ++col; }
   if (i) for (;;)
      { *Mstore = *Cstore++; if (!(--i)) break; Mstore+= ++col; }

}

// routines for row vector

void RowVector::GetCol(MatrixRowCol& mrc)
{
   REPORT
   // do not allow StoreHere
   if (+(mrc.cw*StoreHere))
      Throw(InternalException("RowVector::GetCol(MatrixRowCol&)"));

   mrc.skip=0; mrc.storage=1; mrc.length=nrows; mrc.data = store+mrc.rowcol;

}

void RowVector::GetCol(MatrixColX& mrc)
{
   REPORT
   mrc.skip=0; mrc.storage=1; mrc.length=nrows;
   if (mrc.cw >= LoadOnEntry)
      { REPORT *(mrc.data) = *(store+mrc.rowcol); }

}

void RowVector::NextCol(MatrixRowCol& mrc)
{ REPORT mrc.rowcol++; mrc.data++; }

void RowVector::NextCol(MatrixColX& mrc)
{
   if (+(mrc.cw*StoreOnExit)) { REPORT *(store+mrc.rowcol)=*(mrc.data); }

   mrc.rowcol++;
   if (mrc.rowcol < ncols)
   {
      if (+(mrc.cw*LoadOnEntry)) { REPORT *(mrc.data)=*(store+mrc.rowcol); }
   }
   else { REPORT mrc.cw -= StoreOnExit; }
}

void RowVector::RestoreCol(MatrixColX& mrc)
   { REPORT *(store+mrc.rowcol)=*(mrc.data); }           // not accessed


// routines for band matrices

void BandMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   int r = mrc.rowcol; int w = lower+1+upper; mrc.length=ncols;
   int s = r-lower;
   if (s<0) { mrc.data = store+(r*w-s); w += s; s = 0; }
   else mrc.data = store+r*w;
   mrc.skip = s; s += w-ncols; if (s>0) w -= s; mrc.storage = w;
}

// should make special versions of this for upper and lower band matrices

void BandMatrix::NextRow(MatrixRowCol& mrc)
{
   REPORT
   int r = ++mrc.rowcol;
   if (r<=lower) { mrc.storage++; mrc.data += lower+upper; }
   else  { mrc.skip++; mrc.data += lower+upper+1; }
   if (r>=ncols-upper) mrc.storage--;
}

void BandMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   int c = mrc.rowcol; int n = lower+upper; int w = n+1;
   mrc.length=nrows; Real* ColCopy;
   int b; int s = c-upper;
   if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
   mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
   if ( +(mrc.cw*(StoreHere+HaveStore)) )
      { REPORT ColCopy = mrc.data; }
   else
   {
      REPORT
      ColCopy = new Real [n+1]; MatrixErrorNoSpace(ColCopy);
      MONITOR_REAL_NEW("Make (BMGetCol)",n+1,ColCopy)
      mrc.cw += HaveStore; mrc.data = ColCopy;
   }

   if (+(mrc.cw*LoadOnEntry))
   {
      REPORT
      Real* Mstore = store+b;
      // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
      if (w) for (;;)
         { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
   }
}

void BandMatrix::GetCol(MatrixColX& mrc)
{
   REPORT
   int c = mrc.rowcol; int n = lower+upper; int w = n+1;
   mrc.length=nrows; int b; int s = c-upper;
   if (s<=0) { w += s; s = 0; b = c+lower; } else b = s*w+n;
   mrc.skip = s; s += w-nrows; if (s>0) w -= s; mrc.storage = w;
   mrc.data = mrc.store+mrc.skip;

   if (+(mrc.cw*LoadOnEntry))
   {
      REPORT
      Real* ColCopy = mrc.data; Real* Mstore = store+b;
      // while (w--) { *ColCopy++ = *Mstore; Mstore+=n; }
      if (w) for (;;)
         { *ColCopy++ = *Mstore; if (!(--w)) break; Mstore+=n; }
   }
}

void BandMatrix::RestoreCol(MatrixRowCol& mrc)
{
   REPORT
   int c = mrc.rowcol; int n = lower+upper; int s = c-upper;
   Real* Mstore = store + ((s<=0) ? c+lower : s*n+s+n);
   Real* Cstore = mrc.data;
   int w = mrc.storage;
   // while (w--) { *Mstore = *Cstore++; Mstore += n; }
   if (w) for (;;)
      { *Mstore = *Cstore++; if (!(--w)) break; Mstore += n; }
}

// routines for symmetric band matrix

void SymmetricBandMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   int r=mrc.rowcol; int s = r-lower; int w1 = lower+1; int o = r*w1;
   mrc.length = ncols;
   if (s<0) { w1 += s; o -= s; s = 0; }
   mrc.skip = s;

   if (+(mrc.cw*DirectPart))
      { REPORT  mrc.data = store+o; mrc.storage = w1; }
   else
   {
      // do not allow StoreOnExit and !DirectPart
      if (+(mrc.cw*StoreOnExit))
         Throw(InternalException("SymmetricBandMatrix::GetRow(MatrixRowCol&)"));
      int w = w1+lower; s += w-ncols; Real* RowCopy;
      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;
      if (!(mrc.cw*HaveStore))
      {
         REPORT
         RowCopy = new Real [2*lower+1]; MatrixErrorNoSpace(RowCopy);
         MONITOR_REAL_NEW("Make (SmBGetRow)",2*lower+1,RowCopy)
         mrc.cw += HaveStore; mrc.data = RowCopy;
      }
      else { REPORT  RowCopy = mrc.data; }

      if (+(mrc.cw*LoadOnEntry))
      {
	      REPORT
         Real* Mstore = store+o;
         while (w1--) *RowCopy++ = *Mstore++;
         Mstore--;
         while (w2--) { Mstore += lower; *RowCopy++ = *Mstore; }
      }
   }
}

void SymmetricBandMatrix::GetCol(MatrixRowCol& mrc)
{
   // do not allow StoreHere
   if (+(mrc.cw*StoreHere))
      Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));

   int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
   REPORT
   int s = c-lower; int o = c*w1;
   if (s<0) { w1 += s; o -= s; s = 0; }
   mrc.skip = s;

   if (+(mrc.cw*DirectPart))
   { REPORT  mrc.data = store+o; mrc.storage = w1; }
   else
   {
      // do not allow StoreOnExit and !DirectPart
      if (+(mrc.cw*StoreOnExit))
         Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixRowCol&)"));
      int w = w1+lower; s += w-ncols; Real* ColCopy;
      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;

      if ( +(mrc.cw*HaveStore) ) { REPORT ColCopy = mrc.data; }
      else
      {
         REPORT ColCopy = new Real [2*lower+1]; MatrixErrorNoSpace(ColCopy);
         MONITOR_REAL_NEW("Make (SmBGetCol)",2*lower+1,ColCopy)
         mrc.cw += HaveStore; mrc.data = ColCopy;
      }

      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* Mstore = store+o;
         while (w1--) *ColCopy++ = *Mstore++;
         Mstore--;
         while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
      }
   }
}

void SymmetricBandMatrix::GetCol(MatrixColX& mrc)
{
   int c=mrc.rowcol; int w1 = lower+1; mrc.length=nrows;
   if (+(mrc.cw*DirectPart))
   {
      REPORT
      int b = c*w1+lower;
      mrc.skip = c; c += w1-nrows; w1 -= c; mrc.storage = w1;
      Real* ColCopy = mrc.data = mrc.store+mrc.skip;
      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* Mstore = store+b;
         // while (w1--) { *ColCopy++ = *Mstore; Mstore += lower; }
         if (w1) for (;;)
            { *ColCopy++ = *Mstore; if (!(--w1)) break; Mstore += lower; }
      }
   }
   else
   {
      REPORT
      // do not allow StoreOnExit and !DirectPart
      if (+(mrc.cw*StoreOnExit))
         Throw(InternalException("SymmetricBandMatrix::GetCol(MatrixColX&)"));
      int s = c-lower; int o = c*w1;
      if (s<0) { w1 += s; o -= s; s = 0; }
      mrc.skip = s;

      int w = w1+lower; s += w-ncols;
      if (s>0) w -= s; mrc.storage = w; int w2 = w-w1;

      Real* ColCopy = mrc.data = mrc.store+mrc.skip;

      if (+(mrc.cw*LoadOnEntry))
      {
         REPORT
         Real* Mstore = store+o;
         while (w1--) *ColCopy++ = *Mstore++;
         Mstore--;
         while (w2--) { Mstore += lower; *ColCopy++ = *Mstore; }
      }

   }
}

void SymmetricBandMatrix::RestoreCol(MatrixColX& mrc)
{
   REPORT
   int c = mrc.rowcol;
   Real* Mstore = store + c*lower+c+lower;
   Real* Cstore = mrc.data; int w = mrc.storage;
   // while (w--) { *Mstore = *Cstore++; Mstore += lower; }
   if (w) for (;;)
      { *Mstore = *Cstore++; if (!(--w)) break; Mstore += lower; }
}

// routines for identity matrix

void IdentityMatrix::GetRow(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=mrc.rowcol; mrc.storage=1; mrc.data=store; mrc.length=ncols;
}

void IdentityMatrix::GetCol(MatrixRowCol& mrc)
{
   REPORT
   mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
   if (+(mrc.cw*StoreHere))              // should not happen
      Throw(InternalException("IdentityMatrix::GetCol(MatrixRowCol&)"));
   else  { REPORT mrc.data=store; }
}

void IdentityMatrix::GetCol(MatrixColX& mrc)
{
   REPORT
   mrc.skip=mrc.rowcol; mrc.storage=1; mrc.length=nrows;
   mrc.data = mrc.store+mrc.skip; *(mrc.data)=*store;
}

void IdentityMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }

void IdentityMatrix::NextCol(MatrixRowCol& mrc) { REPORT mrc.IncrId(); }

void IdentityMatrix::NextCol(MatrixColX& mrc)
{
   REPORT
   if (+(mrc.cw*StoreOnExit)) { REPORT *store=*(mrc.data); }
   mrc.IncrDiag();            // must increase mrc.data so need IncrDiag
   int t1 = +(mrc.cw*LoadOnEntry);
   if (t1 && mrc.rowcol < ncols) { REPORT *(mrc.data)=*store; }
}




// *************************** destructors *******************************

MatrixRowCol::~MatrixRowCol()
{
   if (+(cw*HaveStore))
   {
      MONITOR_REAL_DELETE("Free    (RowCol)",-1,data)  // do not know length
      delete [] data;
   }
}

MatrixRow::~MatrixRow() { if (+(cw*StoreOnExit)) gm->RestoreRow(*this); }

MatrixCol::~MatrixCol() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }

MatrixColX::~MatrixColX() { if (+(cw*StoreOnExit)) gm->RestoreCol(*this); }

#ifdef use_namespace
}
#endif