Newsgroups: comp.sources.misc
From: robert@am.dsir.govt.nz (Robert Davies)
Subject:  v26i090:  newmat04 - A matrix package in C++, Part04/05
Message-ID: <1991Nov30.205353.4021@sparky.imd.sterling.com>
X-Md4-Signature: 9805326b081925173c11d229565edab0
Date: Sat, 30 Nov 1991 20:53:53 GMT
Approved: kent@sparky.imd.sterling.com

Submitted-by: robert@am.dsir.govt.nz (Robert Davies)
Posting-number: Volume 26, Issue 90
Archive-name: newmat04/part04
Environment: C++
Supersedes: newmat02: Volume 21, Issue 49-52

#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  newmat5.cxx newmat6.cxx newmat7.cxx newmat8.cxx
#   newmat9.cxx newmatrm.cxx
# Wrapped by robert@kea on Sat Nov 30 18:29:43 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'newmat5.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmat5.cxx'\"
else
echo shar: Extracting \"'newmat5.cxx'\" \(7775 characters\)
sed "s/^X//" >'newmat5.cxx' <<'END_OF_FILE'
X//$$ newmat5.cxx         Transpose, evaluate etc
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X#include "include.hxx"
X
X#include "newmat.hxx"
X#include "newmatrc.hxx"
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,5); ExeCount++; }
X
X#define REPORT {}
X
X
X/************************ carry out operations ******************************/
X
X
XGeneralMatrix* GeneralMatrix::Transpose(MatrixType mt)
X{
X   if ((int)mt==0) mt = Type().t();           // type of transposed matrix
X   GeneralMatrix* gm1 = mt.New(ncols,nrows); gm1->ReleaseAndDelete();
X
X   if (mt == Type().t())
X   {
X      REPORT
X      for (int i=0; i<ncols; i++)
X      {
X	 MatrixRow mr(gm1, StoreOnExit+DirectPart, i);
X         MatrixCol mc(this, mr.Store(), LoadOnEntry, i);
X      }
X   }
X   else
X   {
X      REPORT
X      MatrixRow mr(this, LoadOnEntry);
X      MatrixCol mc(gm1, StoreOnExit+DirectPart);
X      int i = nrows;
X      while (i--) { mc.Copy(mr); mr.Next(); mc.Next(); }
X   }
X   tDelete(); return gm1;
X}
X
XGeneralMatrix* SymmetricMatrix::Transpose(MatrixType mt)
X{ REPORT  return Evaluate(mt); }
X
X
XGeneralMatrix* DiagonalMatrix::Transpose(MatrixType mt)
X{ REPORT return Evaluate(mt); }
X
XBOOL GeneralMatrix::IsZero() const
X{
X   REPORT
X   real* s=store; int i=storage;
X   while (i--) { if (*s++) return FALSE; }
X   return TRUE;
X}
X
XGeneralMatrix* ColumnVector::Transpose(MatrixType mt)
X{
X   REPORT
X   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
X   gmx->nrows = 1; gmx->ncols = gmx->storage = storage;
X   return BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* RowVector::Transpose(MatrixType mt)
X{
X   REPORT
X   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
X   gmx->ncols = 1; gmx->nrows = gmx->storage = storage;
X   return BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* GeneralMatrix::Evaluate(MatrixType mt)
X{
X   if ((int)mt==0) { REPORT return this; }
X   if (mt==this->Type()) { REPORT return this; }
X   REPORT
X   GeneralMatrix* gmx = mt.New(nrows,ncols);
X   MatrixRow mr(this, LoadOnEntry);
X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X   int i=nrows;
X   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
X   tDelete(); gmx->ReleaseAndDelete(); return gmx;
X}
X
XGeneralMatrix* ConstMatrix::Evaluate(MatrixType mt)
X{
X   if ((int)mt==0 || mt==cgm->Type()) { REPORT return (GeneralMatrix*)cgm; }
X
X   REPORT
X   GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols());
X   MatrixRow mr((GeneralMatrix*)cgm, LoadOnEntry);//assume won't change this
X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X   int i=cgm->Nrows();
X   while (i--) { mrx.Copy(mr); mrx.Next(); mr.Next(); }
X   gmx->ReleaseAndDelete(); return gmx; // no tDelete
X}
X
XGeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm=bm->Evaluate();
X   if ((int)mt==0)
X      { MatrixType mteqel(MatrixType::EqEl); mt = gm->Type()+mteqel; }
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   if (!(mt==gm->Type()))
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc);
X      MatrixRow mr(gm, LoadOnEntry); 
X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X      while (nr--) { mrx.Add(mr,f); mrx.Next(); mr.Next(); }
X      gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
X   }
X   else if (gm->reuse()) { REPORT gm->Add(f); return gm; }
X   else
X   {
X      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc);
X      gmy->ReleaseAndDelete(); gmy->Add(gm,f); return gmy;
X   }
X}
X
XGeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm=bm->Evaluate();
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   if ((int)mt==0) mt = gm->Type();
X   if (mt==gm->Type())
X   {
X      if (gm->reuse()) { REPORT gm->Multiply(f); return gm; }
X      else
X      {
X         REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc);
X         gmx->ReleaseAndDelete(); gmx->Multiply(gm,f); return gmx;
X      }
X   }
X   else
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc);
X      MatrixRow mr(gm, LoadOnEntry); 
X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X      while (nr--) { mrx.Multiply(mr,f); mrx.Next(); mr.Next(); }
X      gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
X   }
X}
X
XGeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm=bm->Evaluate();
X   if ((int)mt==0) mt = gm->Type();
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   if (mt==gm->Type())
X   {
X      if (gm->reuse()) { REPORT gm->Negate(); return gm; }
X      else
X      {
X         REPORT
X         GeneralMatrix* gmx = gm->Type().New(nr,nc); gmx->ReleaseAndDelete();
X         gmx->Negate(gm); return gmx;
X      }
X   }
X   else
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc);
X      MatrixRow mr(gm, LoadOnEntry); 
X      MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X      while (nr--) { mrx.Negate(mr); mrx.Next(); mr.Next(); }
X      gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
X   }
X}   
X
XGeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   GeneralMatrix* gm=bm->Evaluate();
X   if ((int)mt==0) mt = gm->Type().t();
X   GeneralMatrix* gmx=gm->Transpose(mt);
X   return gmx;
X}
X   
XGeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm = bm->Evaluate();
X   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
X   gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
X   return gm->BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm = bm->Evaluate();
X   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
X   gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
X   return gm->BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm = bm->Evaluate();
X   GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
X   gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
X   return gm->BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
X{
X   GeneralMatrix* gm = bm->Evaluate();
X   GeneralMatrix* gmx = new Matrix; MatrixErrorNoSpace(gmx);
X   gmx->nrows = nr; gmx->ncols = nc; gmx->storage = gm->storage;
X   if (nr*nc != gmx->storage)
X      MatrixError("Incompatible dimensions in CopyToMatrix");
X   return gm->BorrowStore(gmx,mt);
X}
X
XGeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   GeneralMatrix* gm = bm->Evaluate();
X   if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
X      MatrixError("SubMatrix dimension error");
X   if (int(mt)==0) mt = MatrixType::Rect;
X   GeneralMatrix* gmx = mt.New(row_number, col_number);
X   int i = row_number;
X   MatrixRow mr(gm, LoadOnEntry, row_skip); 
X   MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X   MatrixRowCol sub;
X   while (i--)
X   {
X      mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
X      mrx.Copy(sub); mrx.Next(); mr.Next();
X   }
X   gmx->ReleaseAndDelete(); gm->tDelete(); return gmx;
X}   
X
XGeneralMatrix* ReturnMatrix::Evaluate(MatrixType mt)
X{ return gm->Evaluate(mt); }
X
X
Xvoid GeneralMatrix::Add(GeneralMatrix* gm1, real f)
X{
X   REPORT
X   register real* s1=gm1->store;
X   register real* s=store; register int i=storage;
X   while (i--) *s++ = *s1++ + f;
X}
X   
Xvoid GeneralMatrix::Add(real f)
X{
X   REPORT
X   register real* s=store; register int i=storage; while (i--) *s++ += f;
X}
X   
Xvoid GeneralMatrix::Negate(GeneralMatrix* gm1)
X{
X   // change sign of elements
X   REPORT
X   real* s1=gm1->store; real* s=store; register int i=storage;
X   while(i--) *s++ = -(*s1++);
X}
X   
Xvoid GeneralMatrix::Negate()
X{
X   REPORT
X   real* s=store; register int i=storage; while(i--) { *s = -(*s); s++; } 
X}
X   
Xvoid GeneralMatrix::Multiply(GeneralMatrix* gm1, real f)
X{
X   REPORT
X   register real* s1=gm1->store;
X   register real* s=store; register int i=storage;
X   while (i--) *s++ = *s1++ * f;
X}
X   
Xvoid GeneralMatrix::Multiply(real f)
X{
X   REPORT
X   register real* s=store; register int i=storage; while (i--) *s++ *= f;
X}
X   
END_OF_FILE
if test 7775 -ne `wc -c <'newmat5.cxx'`; then
    echo shar: \"'newmat5.cxx'\" unpacked with wrong size!
fi
# end of 'newmat5.cxx'
fi
if test -f 'newmat6.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmat6.cxx'\"
else
echo shar: Extracting \"'newmat6.cxx'\" \(7084 characters\)
sed "s/^X//" >'newmat6.cxx' <<'END_OF_FILE'
X//$$ newmat6.cxx            Operators, element access, submatrices
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X#include "include.hxx"
X
X#include "newmat.hxx"
X#include "newmatrc.hxx"
X
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,6); ExeCount++; }
X
X#define REPORT {}
X
X/*************************** general utilities *************************/
X
Xstatic int tristore(int n)                      // els in triangular matrix
X{ return (n*(n+1))/2; }
X
X
X/****************************** operators *******************************/
X
Xreal& Matrix::operator()(int m, int n)
X{
X   if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
X   return store[(m-1)*ncols+n-1];
X}
X
Xreal& SymmetricMatrix::operator()(int m, int n)
X{
X   if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
X   if (m>=n) return store[tristore(m-1)+n-1];
X   else return store[tristore(n-1)+m-1];
X}
X
Xreal& UpperTriangularMatrix::operator()(int m, int n)
X{
X   if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
X   return store[(m-1)*ncols+n-1-tristore(m-1)];
X}
X
Xreal& LowerTriangularMatrix::operator()(int m, int n)
X{
X   if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
X   return store[tristore(m-1)+n-1];
X}
X
Xreal& DiagonalMatrix::operator()(int m, int n)
X{
X   if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
X   return store[n-1];
X}
X
Xreal& DiagonalMatrix::operator()(int m)
X{
X   if (m<=0 || m>nrows) MatrixError("Index out of range");
X   return store[m-1];
X}
X
Xreal& ColumnVector::operator()(int m)
X{
X   if (m<=0 || m> nrows) MatrixError("Index out of range");
X   return store[m-1];
X}
X
Xreal& RowVector::operator()(int n)
X{
X   if (n<=0 || n> ncols) MatrixError("Index out of range");
X   return store[n-1];
X}
X
X#ifndef __ZTC__
X
Xreal Matrix::operator()(int m, int n) const
X{
X   if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
X   return store[(m-1)*ncols+n-1];
X}
X
Xreal SymmetricMatrix::operator()(int m, int n) const
X{
X   if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
X   if (m>=n) return store[tristore(m-1)+n-1];
X   else return store[tristore(n-1)+m-1];
X}
X
Xreal UpperTriangularMatrix::operator()(int m, int n) const
X{
X   if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
X   return store[(m-1)*ncols+n-1-tristore(m-1)];
X}
X
Xreal LowerTriangularMatrix::operator()(int m, int n) const
X{
X   if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
X   return store[tristore(m-1)+n-1];
X}
X
Xreal DiagonalMatrix::operator()(int m, int n) const
X{
X   if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
X   return store[n-1];
X}
X
Xreal DiagonalMatrix::operator()(int m) const
X{
X   if (m<=0 || m>nrows) MatrixError("Index out of range");
X   return store[m-1];
X}
X
Xreal ColumnVector::operator()(int m) const
X{
X   if (m<=0 || m> nrows) MatrixError("Index out of range");
X   return store[m-1];
X}
X
Xreal RowVector::operator()(int n) const
X{
X   if (n<=0 || n> ncols) MatrixError("Index out of range");
X   return store[n-1];
X}
X
X#endif
X
XBaseMatrix::operator real()
X{
X   REPORT
X   GeneralMatrix* gm = Evaluate();
X   if (gm->nrows!=1 || gm->ncols!=1)
X      MatrixError("Attempt to convert non 1x1 matrix to scalar");
X   real x = *(gm->store); gm->tDelete(); return x;
X}
X
XAddedMatrix BaseMatrix::operator+(BaseMatrix& bm)
X{ REPORT return AddedMatrix(this, &bm); }
X
XMultipliedMatrix BaseMatrix::operator*(BaseMatrix& bm)
X{ REPORT return MultipliedMatrix(this, &bm); }
X
XSolvedMatrix InvertedMatrix::operator*(BaseMatrix& bmx)
X{ REPORT return SolvedMatrix(bm, &bmx); }
X
XSubtractedMatrix BaseMatrix::operator-(BaseMatrix& bm)
X{ REPORT return SubtractedMatrix(this, &bm); }
X
XShiftedMatrix BaseMatrix::operator+(real f)
X{ REPORT return ShiftedMatrix(this, f); }
X
XScaledMatrix BaseMatrix::operator*(real f)
X{ REPORT return ScaledMatrix(this, f); }
X
XScaledMatrix BaseMatrix::operator/(real f)
X{ REPORT return ScaledMatrix(this, 1.0/f); }
X
XShiftedMatrix BaseMatrix::operator-(real f)
X{ REPORT return ShiftedMatrix(this, -f); }
X
XTransposedMatrix BaseMatrix::t() { REPORT return TransposedMatrix(this); }
X
XNegatedMatrix BaseMatrix::operator-() { REPORT return NegatedMatrix(this); }
X
XInvertedMatrix BaseMatrix::i() { REPORT return InvertedMatrix(this); }
X
XConstMatrix GeneralMatrix::c() const
X{
X   if (tag != -1) MatrixError(".c() applied to temporary matrix");
X   REPORT return ConstMatrix(this);
X}
X
XRowedMatrix BaseMatrix::CopyToRow() { REPORT return RowedMatrix(this); }
X
XColedMatrix BaseMatrix::CopyToColumn() { REPORT return ColedMatrix(this); }
X
XDiagedMatrix BaseMatrix::CopyToDiagonal() { REPORT return DiagedMatrix(this); }
X
XMatedMatrix BaseMatrix::CopyToMatrix(int nrx, int ncx)
X{ REPORT return MatedMatrix(this,nrx,ncx); }
X
Xvoid GeneralMatrix::operator=(real f)
X{ REPORT int i=storage; real* s=store; while (i--) { *s++ = f; } }
X
Xvoid Matrix::operator=(BaseMatrix& X)
X{ REPORT CheckConversion(X); Eq(X,MatrixType::Rect); } 
X
Xvoid RowVector::operator=(BaseMatrix& X)
X{
X   REPORT CheckConversion(X); Eq(X,MatrixType::RowV);
X   if (nrows!=1) MatrixError("Illegal conversion to row vector");
X}
X
Xvoid ColumnVector::operator=(BaseMatrix& X)
X{
X   REPORT CheckConversion(X); Eq(X,MatrixType::ColV);
X   if (ncols!=1) MatrixError("Illegal conversion to column vector");
X}
X
Xvoid SymmetricMatrix::operator=(BaseMatrix& X)
X{ REPORT CheckConversion(X); Eq(X,MatrixType::Sym); }
X 
Xvoid UpperTriangularMatrix::operator=(BaseMatrix& X)
X{ REPORT CheckConversion(X); Eq(X,MatrixType::UT); }
X
Xvoid LowerTriangularMatrix::operator=(BaseMatrix& X)
X{ REPORT CheckConversion(X); Eq(X,MatrixType::LT); }
X
Xvoid DiagonalMatrix::operator=(BaseMatrix& X)
X{ REPORT CheckConversion(X); Eq(X,MatrixType::Diag); }
X
Xvoid GeneralMatrix::operator<<(const real* r)
X{
X   REPORT
X   int i = storage; real* s=store;
X   while(i--) *s++ = *r++;
X}
X
X
X/************************* element access *********************************/
X
Xreal& Matrix::element(int m, int n)
X{
X   if (m<0 || m>= nrows || n<0 || n>= ncols) MatrixError("Index out of range");
X   return store[m*ncols+n];
X}
X
Xreal& SymmetricMatrix::element(int m, int n)
X{
X   if (m<0 || n<0 || m >= nrows || n>=ncols) MatrixError("Index out of range");
X   if (m>=n) return store[tristore(m)+n];
X   else return store[tristore(n)+m];
X}
X
Xreal& UpperTriangularMatrix::element(int m, int n)
X{
X   if (m<0 || n<m || n>=ncols) MatrixError("Index out of range");
X   return store[m*ncols+n-tristore(m)];
X}
X
Xreal& LowerTriangularMatrix::element(int m, int n)
X{
X   if (n<0 || m<n || m>=nrows) MatrixError("Index out of range");
X   return store[tristore(m)+n];
X}
X
Xreal& DiagonalMatrix::element(int m, int n)
X{
X   if (n<0 || m!=n || m>=nrows || n>=ncols) MatrixError("Index out of range");
X   return store[n];
X}
X
Xreal& DiagonalMatrix::element(int m)
X{
X   if (m<0 || m>=nrows) MatrixError("Index out of range");
X   return store[m];
X}
X
Xreal& ColumnVector::element(int m)
X{
X   if (m<0 || m>= nrows) MatrixError("Index out of range");
X   return store[m];
X}
X
Xreal& RowVector::element(int n)
X{
X   if (n<0 || n>= ncols) MatrixError("Index out of range");
X   return store[n];
X}
X
END_OF_FILE
if test 7084 -ne `wc -c <'newmat6.cxx'`; then
    echo shar: \"'newmat6.cxx'\" unpacked with wrong size!
fi
# end of 'newmat6.cxx'
fi
if test -f 'newmat7.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmat7.cxx'\"
else
echo shar: Extracting \"'newmat7.cxx'\" \(11535 characters\)
sed "s/^X//" >'newmat7.cxx' <<'END_OF_FILE'
X//$$ newmat7.cxx     Invert, solve, binary operations
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X#include "include.hxx"
X
X#include "newmat.hxx"
X#include "newmatrc.hxx"
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,7); ExeCount++; }
X
X#define REPORT {}
X
X
X/***************************** solve routines ******************************/
X
XGeneralMatrix* GeneralMatrix::MakeSolver()
X{
X   REPORT
X   GeneralMatrix* gm = new CroutMatrix(*this);
X   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
X}
X
XGeneralMatrix* Matrix::MakeSolver()
X{
X   REPORT
X   GeneralMatrix* gm = new CroutMatrix(*this);
X   MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
X}
X
Xvoid CroutMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
X{
X   REPORT
X   real* el = mcin.store; int i = mcin.skip;
X   while (i--) *el++ = 0.0;
X   el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
X   while (i--) *el++ = 0.0;
X   lubksb(mcin.store, mcout.skip);
X}
X
X
X// Do we need check for entirely zero output?
X
Xvoid UpperTriangularMatrix::Solver(MatrixRowCol& mcout,
X   const MatrixRowCol& mcin)
X{
X   REPORT
X   real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
X   while (i-- > 0) *elx++ = 0.0;
X   int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; real* el = elx;
X   int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
X   while (j-- > 0) *elx++ = 0.0;
X   real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
X   while (i-- > 0)
X   {
X      elx = el; real sum = 0.0; int jx = j++; Ael -= nc;
X      while (jx--) sum += *(--Ael) * *(--elx);
X      elx--; *elx = (*elx - sum) / *(--Ael);
X   }
X}
X
Xvoid LowerTriangularMatrix::Solver(MatrixRowCol& mcout,
X   const MatrixRowCol& mcin)
X{
X   REPORT
X   real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
X   while (i-- > 0) *elx++ = 0.0;
X   int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
X   int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
X   while (j-- > 0) *elx++ = 0.0;
X   real* el = mcin.store+nc; real* Ael = store + (nc*(nc+1))/2; j = 0;
X   while (i-- > 0)
X   {
X      elx = el; real sum = 0.0; int jx = j++; Ael += nc;
X      while (jx--) sum += *Ael++ * *elx++;
X      *elx = (*elx - sum) / *Ael++;
X   }
X}
X
X/******************* carry out binary operations *************************/
X
Xstatic void Add(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
Xstatic void Add(GeneralMatrix*,GeneralMatrix*);
X
Xstatic void Subtract(GeneralMatrix*,GeneralMatrix*, GeneralMatrix*);
Xstatic void Subtract(GeneralMatrix*,GeneralMatrix*);
Xstatic void ReverseSubtract(GeneralMatrix*,GeneralMatrix*);
X
Xstatic GeneralMatrix* GeneralAdd(GeneralMatrix*,GeneralMatrix*,MatrixType);
Xstatic GeneralMatrix* GeneralSub(GeneralMatrix*,GeneralMatrix*,MatrixType);
Xstatic GeneralMatrix* GeneralMult(GeneralMatrix*,GeneralMatrix*,MatrixType);
Xstatic GeneralMatrix* GeneralSolv(GeneralMatrix*,GeneralMatrix*,MatrixType);
X
XGeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
X{ REPORT return GeneralAdd(bm1->Evaluate(), bm2->Evaluate(), mt); }
X
XGeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
X{ REPORT return GeneralSub(bm1->Evaluate(), bm2->Evaluate(), mt); }
X
XGeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   GeneralMatrix* gm2 = bm2->Evaluate();
X   MatrixType mtsym(MatrixType::Sym);
X   if (gm2->Type() == mtsym) gm2 = gm2->Evaluate(MatrixType::Rect);
X   return GeneralMult(bm1->Evaluate(), gm2, mt);
X}
X
XGeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
X{ REPORT return GeneralSolv(bm1->Evaluate(), bm2->Evaluate(), mt); }
X
X// routines for adding or subtracting matrices of identical storage structure
X
Xstatic void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   REPORT
X   register real* s1=gm1->Store(); register real* s2=gm2->Store();
X   register real* s=gm->Store(); register int i=gm->Storage();
X   while (i--) *s++ = *s1++ + *s2++;
X}
X   
Xstatic void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   register real* s2=gm2->Store();
X   register real* s=gm->Store(); register int i=gm->Storage();
X   while (i--) *s++ += *s2++;
X}
X
Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   REPORT
X   register real* s1=gm1->Store(); register real* s2=gm2->Store();
X   register real* s=gm->Store(); register int i=gm->Storage();
X   while (i--) *s++ = *s1++ - *s2++;
X}
X
Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   register real* s2=gm2->Store();
X   register real* s=gm->Store(); register int i=gm->Storage();
X   while (i--) *s++ -= *s2++;
X}
X
Xstatic void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   register real* s2=gm2->Store();
X   register real* s=gm->Store(); register int i=gm->Storage();
X   while (i--) { *s = *s2++ - *s; s++; }
X}
X
Xstatic GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MatrixType mtx)
X{
X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
X      MatrixError("Incompatible dimensions");
X   if (!mtx) mtx = gm1->Type() + gm2->Type();
X   if (mtx == gm1->Type() && mtx == gm2->Type())
X					     // modify when band or sparse
X					     // matrices are included.
X   {
X      if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
X      else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
X      else
X      {
X         REPORT
X         GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
X         gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
X      }
X   }
X   else
X   {
X      if (gm1->Type()==mtx && gm1->reuse() )    // must have type test first
X      {
X	 REPORT
X	 MatrixRow mr1(gm1, StoreOnExit+LoadOnEntry+DirectPart);
X	 MatrixRow mr2(gm2, LoadOnEntry);
X	 while (nr--) { mr1.Add(mr2); mr1.Next(); mr2.Next(); }
X         gm2->tDelete(); return gm1;
X      }
X      else if (gm2->Type()==mtx && gm2->reuse() )
X      {
X         REPORT
X	 MatrixRow mr1(gm1, LoadOnEntry);
X	 MatrixRow mr2(gm2, StoreOnExit+LoadOnEntry+DirectPart);
X	 while (nr--) { mr2.Add(mr1); mr1.Next(); mr2.Next(); }
X	 if (gm1->Type()!=mtx) gm1->tDelete();
X         return gm2;
X      }
X      else
X      {
X         REPORT
X	 GeneralMatrix* gmx = mtx.New(nr,nc);
X         MatrixRow mr1(gm1, LoadOnEntry);
X	 MatrixRow mr2(gm2, LoadOnEntry);
X	 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X	 while (nr--)
X	 { mrx.Add(mr1,mr2); mr1.Next(); mr2.Next(); mrx.Next(); }
X	 if (gm1->Type()!=mtx) gm1->tDelete();
X	 if (gm2->Type()!=mtx) gm2->tDelete();
X         gmx->ReleaseAndDelete(); return gmx;
X      }
X   }
X}
X
X
Xstatic GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MatrixType mtx)
X{
X   if (!mtx) mtx = gm1->Type() - gm2->Type();
X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
X      MatrixError("Incompatible dimensions");
X   if (mtx == gm1->Type() && mtx == gm2->Type())
X					     // modify when band or sparse
X                                             // matrices are included.
X   {
X      if (gm1->reuse())
X      { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
X      else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
X      else
X      {
X         REPORT
X	 GeneralMatrix* gmx = gm1->Type().New(nr,nc); gmx->ReleaseAndDelete();
X         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
X      }
X   }
X   else
X   {
X      if (gm1->Type() == mtx && gm1->reuse() )
X      {
X         REPORT
X	 MatrixRow mr1(gm1, LoadOnEntry+StoreOnExit+DirectPart);
X	 MatrixRow mr2(gm2, LoadOnEntry);
X	 while (nr--) { mr1.Sub(mr2); mr1.Next(); mr2.Next(); }
X         gm2->tDelete(); return gm1;
X      }
X      else if (gm2->Type() == mtx && gm2->reuse() )
X      {
X         REPORT
X         MatrixRow mr1(gm1, LoadOnEntry);
X	 MatrixRow mr2(gm2, LoadOnEntry+StoreOnExit+DirectPart);
X	 while (nr--) { mr2.RevSub(mr1); mr1.Next(); mr2.Next(); }
X	 if (gm1->Type()!=mtx) gm1->tDelete();
X	 return gm2;
X      }
X      else
X      {
X         REPORT
X	 GeneralMatrix* gmx = mtx.New(nr,nc);
X         MatrixRow mr1(gm1, LoadOnEntry);
X         MatrixRow mr2(gm2, LoadOnEntry);
X	 MatrixRow mrx(gmx, StoreOnExit+DirectPart);
X	 while (nr--)
X	    { mrx.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mrx.Next(); }
X	 if (gm1->Type()!=mtx) gm1->tDelete();
X	 if (gm2->Type()!=mtx) gm2->tDelete();
X	 gmx->ReleaseAndDelete(); return gmx;
X      }
X   }
X}
X
X/*
Xstatic GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MatrixType mtx)
X{
X   REPORT
X   if (!mtx) mtx = gm1->Type() * gm2->Type();
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
X   GeneralMatrix* gmx = mtx.New(nr,nc);
X
X   MatrixCol mcx(gmx, StoreOnExit+DirectPart);
X   MatrixCol mc2(gm2, LoadOnEntry);
X   while (nc--)
X   {
X      MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
X      real* el = mcx();                              // pointer to an element
X      int n = mcx.Storage();
X      while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
X      mc2.Next(); mcx.Next();
X   }
X   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
X}
X*/
X
Xstatic GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MatrixType mtx)
X{
X   REPORT
X   if (!mtx) mtx = gm1->Type() * gm2->Type();
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
X   GeneralMatrix* gmx = mtx.New(nr,nc);
X
X   real* el = gmx->Store(); int n = gmx->Storage();
X   while (n--) *el++ = 0.0;
X   MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
X   MatrixRow mr1(gm1, LoadOnEntry);
X   while (nr--)
X   {
X      MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
X      el = mr1();                              // pointer to an element
X      n = mr1.Storage();
X      while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
X      mr1.Next(); mrx.Next();
X   }
X   gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
X}
X
Xstatic GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MatrixType mtx)
X{
X   REPORT
X   if (!mtx) mtx = gm1->Type().i()*gm2->Type();
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows()) MatrixError("Incompatible dimensions");
X   GeneralMatrix* gmx = mtx.New(nr,nc);
X
X   real* r = new real [nr]; MatrixErrorNoSpace(r);
X   MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
X   MatrixCol mc2(gm2, r, LoadOnEntry);
X   GeneralMatrix* gms = gm1->MakeSolver();
X   while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
X   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
X   delete [nr] r;
X   return gmx;
X}
X
XGeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
X{
X   // Matrix Inversion - use solve routines
X   REPORT
X   GeneralMatrix* gm=bm->Evaluate();
X   int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
X   return GeneralSolv(gm,&I,mtx);
X}
X
X
X/*************************** norm functions ****************************/
X
Xreal BaseMatrix::Norm1()
X{
X   // maximum of sum of absolute values of a column
X   REPORT
X   GeneralMatrix* gm = Evaluate(); int nc = gm->Ncols(); real value = 0.0;
X   MatrixCol mc(gm, LoadOnEntry);
X   while (nc--)
X      { real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
X   gm->tDelete(); return value;
X}
X
Xreal BaseMatrix::NormInfinity()
X{
X   // maximum of sum of absolute values of a row
X   REPORT
X   GeneralMatrix* gm = Evaluate(); int nr = gm->Nrows(); real value = 0.0;
X   MatrixRow mr(gm, LoadOnEntry);
X   while (nr--)
X      { real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
X   gm->tDelete(); return value;
X}
END_OF_FILE
if test 11535 -ne `wc -c <'newmat7.cxx'`; then
    echo shar: \"'newmat7.cxx'\" unpacked with wrong size!
fi
# end of 'newmat7.cxx'
fi
if test -f 'newmat8.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmat8.cxx'\"
else
echo shar: Extracting \"'newmat8.cxx'\" \(6243 characters\)
sed "s/^X//" >'newmat8.cxx' <<'END_OF_FILE'
X//$$ newmat8.cxx         Advanced LU transform, scalar functions
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X#define WANT_MATH
X
X#include "include.hxx"
X
X#include "newmatap.hxx"
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,8); ExeCount++; }
X
X#define REPORT {}
X
X
X/************************** LU transformation ****************************/
X
Xvoid CroutMatrix::ludcmp()
X// LU decomposition - from numerical recipes in C
X{
X   REPORT
X   int i,j;
X
X   real* vv=new real [nrows]; MatrixErrorNoSpace(vv);
X   real* a;
X
X   a=store;
X   for (i=0;i<nrows;i++)
X   {
X      real big=0.0;
X      j=nrows; while (j--) { real temp=fabs(*a++); if (temp > big) big=temp; }
X      if (big == 0.0) MatrixError("Singular matrix");
X      vv[i]=1.0/big;
X   }
X
X   real* aj=store;
X   for (j=0;j<nrows;j++)
X   {
X      real* ai=store;
X      for (i=0;i<j;i++)
X      {
X         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
X         int k=i; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
X         *ajx = sum; ai += nrows;
X      }
X
X      real big=0.0; int imax;
X      for (i=j;i<nrows;i++)
X      {
X         real sum=*(ai+j); real* aix=ai; real* ajx=aj;
X         int k=j; while (k--) { sum -= *aix++ * *ajx; ajx += nrows; }
X         *aix = sum; ai += nrows;
X         real dum=vv[i]*fabs(sum); if (dum >= big) { big=dum; imax=i; }
X      }
X
X      if (j != imax)
X      {
X         real* amax=store+imax*nrows; real* ajx=store+j*nrows;
X         int k=nrows; while(k--) { real dum=*amax; *amax++=*ajx; *ajx++=dum; }
X         d=!d; vv[imax]=vv[j];
X      }
X
X      indx[j]=imax; ai=aj+j*nrows;
X      if (*ai == 0.0) MatrixError("Singular matrix");
X      real dum=1.0/(*ai);
X      i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
X
X      aj++;
X   }
X   delete [nrows] vv;
X}
X
Xvoid CroutMatrix::lubksb(real* B, int mini)
X{
X   REPORT
X   int i,j; int ii=-1; real* ai=store;
X
X   for (i=0;i<nrows;i++)
X   {
X      int ip=indx[i]; real sum=B[ip]; B[ip]=B[i];
X      if (ii>=0)
X      {
X         real* aij=ai+ii; real* bj=B+ii; j=i-ii;
X         while (j--) { sum -= *aij++ * *bj; bj++; }
X      }
X      else if (sum) ii=i;
X      B[i]=sum; ai += nrows;
X   }
X
X   for (i=nrows-1;i>=mini;i--)
X   {
X      real* bj=B+i; ai -= nrows; real* ajx=ai+i; real sum=*bj; bj++;
X      j=nrows-i; while(--j) { sum -= *(++ajx) * *bj; bj++; }
X      B[i]=sum/(*(ai+i));
X   }
X}
X
X
X/****************************** scalar functions ****************************/
X
Xinline real square(real x) { return x*x; }
X
Xreal GeneralMatrix::SumSquare()
X{
X   REPORT
X   real sum = 0.0; int i = storage; real* s = store;
X   while (i--) sum += square(*s++);
X   tDelete(); return sum;
X}
X
Xreal GeneralMatrix::SumAbsoluteValue()
X{
X   REPORT
X   real sum = 0.0; int i = storage; real* s = store;
X   while (i--) sum += fabs(*s++);
X   tDelete(); return sum;
X}
X
Xreal GeneralMatrix::MaximumAbsoluteValue()
X{
X   REPORT
X   real maxval = 0.0; int i = storage; real* s = store;
X   while (i--) { real a = fabs(*s++); if (maxval < a) maxval = a; }
X   tDelete(); return maxval;
X}
X
Xreal SymmetricMatrix::SumSquare()
X{
X   REPORT
X   real sum1 = 0.0; real sum2 = 0.0; real* s = store; int nr = nrows;
X   for (int i = 0; i<nr; i++)
X   {
X      int j = i;
X      while (j--) sum2 += square(*s++);
X      sum1 += square(*s++);
X   }
X   tDelete(); return sum1 + 2.0 * sum2;
X}
X
Xreal SymmetricMatrix::SumAbsoluteValue()
X{
X   REPORT
X   real sum1 = 0.0; real sum2 = 0.0; real* s = store; int nr = nrows;
X   for (int i = 0; i<nr; i++)
X   {
X      int j = i;
X      while (j--) sum2 += fabs(*s++);
X      sum1 += fabs(*s++);
X   }
X   tDelete(); return sum1 + 2.0 * sum2;
X}
X
Xreal BaseMatrix::SumSquare()
X{
X   REPORT GeneralMatrix* gm = Evaluate();
X   real s = gm->SumSquare(); return s;
X}
X
Xreal BaseMatrix::SumAbsoluteValue()
X{
X   REPORT GeneralMatrix* gm = Evaluate();
X   real s = gm->SumAbsoluteValue(); return s;
X}
X
Xreal BaseMatrix::MaximumAbsoluteValue()
X{
X   REPORT GeneralMatrix* gm = Evaluate();
X   real s = gm->MaximumAbsoluteValue(); return s;
X}
X
Xreal Matrix::Trace()
X{
X   REPORT
X   int i = nrows; int d = i+1;
X   if (i != ncols) MatrixError("trace of non-square matrix");
X   real sum = 0.0; real* s = store;
X   while (i--) { sum += *s; s += d; }
X   tDelete(); return sum;
X}
X
Xreal DiagonalMatrix::Trace()
X{
X   REPORT
X   int i = nrows; real sum = 0.0; real* s = store;
X   while (i--) sum += *s++;
X   tDelete(); return sum;
X}
X
Xreal SymmetricMatrix::Trace()
X{
X   REPORT
X   int i = nrows; real sum = 0.0; real* s = store; int j = 2;
X   while (i--) { sum += *s; s += j++; }
X   tDelete(); return sum;
X}
X
Xreal LowerTriangularMatrix::Trace()
X{
X   REPORT
X   int i = nrows; real sum = 0.0; real* s = store; int j = 2;
X   while (i--) { sum += *s; s += j++; }
X   tDelete(); return sum;
X}
X
Xreal UpperTriangularMatrix::Trace()
X{
X   REPORT
X   int i = nrows; real sum = 0.0; real* s = store;
X   while (i) { sum += *s; s += i--; }
X   tDelete(); return sum;
X}
X
Xreal BaseMatrix::Trace()
X{
X   REPORT
X   GeneralMatrix* gm = Evaluate(MatrixType::Diag);
X   real sum = gm->Trace(); return sum;
X}
X
Xvoid LogAndSign::operator*=(real x)
X{
X   if (x > 0.0) { log_value += log(x); }
X   else if (x < 0.0) { log_value += log(-x); sign = -sign; }
X   else sign = 0;
X}
X
Xreal LogAndSign::Value() { return sign * exp(log_value); }
X
XLogAndSign DiagonalMatrix::LogDeterminant()
X{
X   REPORT
X   int i = nrows; LogAndSign sum; real* s = store;
X   while (i--) sum *= *s++;
X   tDelete(); return sum;
X}
X
XLogAndSign LowerTriangularMatrix::LogDeterminant()
X{
X   REPORT
X   int i = nrows; LogAndSign sum; real* s = store; int j = 2;
X   while (i--) { sum *= *s; s += j++; }
X   tDelete(); return sum;
X}
X
XLogAndSign UpperTriangularMatrix::LogDeterminant()
X{
X   REPORT
X   int i = nrows; LogAndSign sum; real* s = store;
X   while (i) { sum *= *s; s += i--; }
X   tDelete(); return sum;
X}
X
XLogAndSign BaseMatrix::LogDeterminant()
X{
X   REPORT GeneralMatrix* gm = Evaluate();
X   LogAndSign sum = gm->LogDeterminant(); return sum;
X}
X
XLogAndSign GeneralMatrix::LogDeterminant()
X{
X   REPORT
X   if (nrows != ncols) MatrixError("determinant of non-square matrix"); 
X   CroutMatrix C(*this); return C.LogDeterminant();
X}
X
XLogAndSign CroutMatrix::LogDeterminant()
X{
X   REPORT
X   int i = nrows; int dd = i+1; LogAndSign sum; real* s = store;
X   while (i--) { sum *= *s; s += dd; }
X   if (!d) sum.ChangeSign(); return sum;
X
X}
END_OF_FILE
if test 6243 -ne `wc -c <'newmat8.cxx'`; then
    echo shar: \"'newmat8.cxx'\" unpacked with wrong size!
fi
# end of 'newmat8.cxx'
fi
if test -f 'newmat9.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmat9.cxx'\"
else
echo shar: Extracting \"'newmat9.cxx'\" \(946 characters\)
sed "s/^X//" >'newmat9.cxx' <<'END_OF_FILE'
X//$$ newmat9.cxx         Input and output
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X
X#define WANT_STREAM
X
X#include "include.hxx"
X
X#include "newmat.hxx"
X#include "newmatrc.hxx"
X#include "newmatio.hxx"
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,9); ExeCount++; }
X
X#define REPORT {}
X
Xostream& operator<<(ostream& s, BaseMatrix& X)
X{
X   GeneralMatrix* gm = X.Evaluate(); operator<<(s, *gm);
X   gm->tDelete(); return s;
X}
X
X
Xostream& operator<<(ostream& s, const GeneralMatrix& X)
X{
X   int w = s.width();  int nr = X.Nrows();  long f = s.flags();
X   s.setf(ios::fixed, ios::floatfield);
X   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
X   for (int i=1; i<=nr; i++)
X   {
X      int skip = mr.skip;  int storage = mr.storage;
X      real* store = mr.store+skip;  skip *= w+1;
X      while (skip--) s << " ";
X      while (storage--) s << setw(w) << *store++ << " ";
X      mr.Next();  s << "\n";
X   }
X   s << flush;  s.flags(f); return s;
X}
X
END_OF_FILE
if test 946 -ne `wc -c <'newmat9.cxx'`; then
    echo shar: \"'newmat9.cxx'\" unpacked with wrong size!
fi
# end of 'newmat9.cxx'
fi
if test -f 'newmatrm.cxx' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'newmatrm.cxx'\"
else
echo shar: Extracting \"'newmatrm.cxx'\" \(3141 characters\)
sed "s/^X//" >'newmatrm.cxx' <<'END_OF_FILE'
X//$$newmatrm.cxx                         rectangular matrix operations
X
X// Copyright (C) 1991: R B Davies and DSIR
X
X#include "include.hxx"
X
X#include "newmat.hxx"
X#include "newmatrm.hxx"
X
X
X// operations on rectangular matrices
X
X
Xvoid RectMatrixRow::Reset (const Matrix& M, int row, int skip, int length)
X{
X   RectMatrixRowCol::Reset
X      ( M.Store()+row*M.Ncols()+skip, length, 1, M.Ncols() );
X}
X
Xvoid RectMatrixRow::Reset (const Matrix& M, int row)
X{
X   RectMatrixRowCol::Reset( M.Store()+row*M.Ncols(), M.Ncols(), 1, M.Ncols() );
X}
X
Xvoid RectMatrixCol::Reset (const Matrix& M, int skip, int col, int length)
X{
X   RectMatrixRowCol::Reset
X      ( M.Store()+col+skip*M.Ncols(), length, M.Ncols(), 1 );
X}
X
Xvoid RectMatrixCol::Reset (const Matrix& M, int col)
X   { RectMatrixRowCol::Reset( M.Store()+col, M.Nrows(), M.Ncols(), 1 ); }
X
X
Xreal RectMatrixRowCol::SumSquare() const
X{
X   long_real sum = 0.0; int i = n; real* s = store; int d = spacing;
X   while (i--) { sum += (long_real)*s * *s; s += d; }
X   return (real)sum;
X}
X
Xreal RectMatrixRowCol::operator*(const RectMatrixRowCol& rmrc) const
X{
X   long_real sum = 0.0; int i = n;
X   real* s = store; int d = spacing;
X   real* s1 = rmrc.store; int d1 = rmrc.spacing;
X   if (i!=rmrc.n) MatrixError("Row/column length error in *");
X   while (i--) { sum += (long_real)*s * *s1; s += d; s1 += d1; }
X   return (real)sum;
X}
X
Xvoid RectMatrixRowCol::AddScaled(const RectMatrixRowCol& rmrc, real r)
X{
X   int i = n; real* s = store; int d = spacing;
X   real* s1 = rmrc.store; int d1 = rmrc.spacing;
X   if (i!=rmrc.n) MatrixError("Row/column length error in AddScaled");
X   while (i--) { *s += *s1 * r; s += d; s1 += d1; }
X}
X
Xvoid RectMatrixRowCol::Divide(const RectMatrixRowCol& rmrc, real r)
X{
X   int i = n; real* s = store; int d = spacing;
X   real* s1 = rmrc.store; int d1 = rmrc.spacing;
X   if (i!=rmrc.n) MatrixError("Row/column length error in Divide");
X   while (i--) { *s = *s1 / r; s += d; s1 += d1; }
X}
X
Xvoid RectMatrixRowCol::Divide(real r)
X{
X   int i = n; real* s = store; int d = spacing;
X   while (i--) { *s /= r; s += d; }
X}
X
Xvoid RectMatrixRowCol::Negate()
X{
X   int i = n; real* s = store; int d = spacing;
X   while (i--) { *s = - *s; s += d; }
X}
X
Xvoid RectMatrixRowCol::Zero()
X{
X   int i = n; real* s = store; int d = spacing;
X   while (i--) { *s = 0.0; s += d; }
X}
X
Xvoid ComplexScale(RectMatrixCol& U, RectMatrixCol& V, real x, real y)
X{
X   int n = U.n;
X   if (n != V.n) MatrixError("Columns different length in ComplexScale");
X   real* u = U.store; real* v = V.store; 
X   int su = U.spacing; int sv = V.spacing;
X   while (n--)
X   {
X      real z = *u * x - *v * y;  *v =  *u * y + *v * x;  *u = z;
X      u += su;  v += sv;
X   }
X}
X
Xvoid Rotate(RectMatrixCol& U, RectMatrixCol& V, real tau, real s)
X{
X   //  (U, V) = (U, V) * (c, s)  where  tau = s/(1+c), c^2 + s^2 = 1
X   int n = U.n;
X   if (n != V.n) MatrixError("Columns different length in Rotate");
X   real* u = U.store; real* v = V.store; 
X   int su = U.spacing; int sv = V.spacing;
X   while (n--)
X   {
X      real zu = *u; real zv = *v;
X      *u -= s * (zv + zu * tau); *v += s * (zu - zv * tau);
X      u += su;  v += sv;
X   }
X}
X
X
END_OF_FILE
if test 3141 -ne `wc -c <'newmatrm.cxx'`; then
    echo shar: \"'newmatrm.cxx'\" unpacked with wrong size!
fi
# end of 'newmatrm.cxx'
fi
echo shar: End of shell archive.
exit 0

exit 0 # Just in case...
-- 
Kent Landfield                   INTERNET: kent@sparky.IMD.Sterling.COM
Sterling Software, IMD           UUCP:     uunet!sparky!kent
Phone:    (402) 291-8300         FAX:      (402) 291-4362
Please send comp.sources.misc-related mail to kent@uunet.uu.net.