Source Code 1994 March/Source_Code_CD-ROM_Walnut_Creek_March_1994.iso/compsrcs/misc/volume34/newmat07/part05

Text File | 1993-01-10 | 58.0 KB | 2,173 lines

Newsgroups: comp.sources.misc
From: [e-mail] (Robert Davies)
Subject: v34i111:  newmat07 - A matrix package in C++, Part05/08
Message-ID: <[e-mail]>
X-Md4-Signature: 95c3bfe3e634dc03ab1b546ccee8137a
Date: Mon, 11 Jan 1993 15:32:15 GMT
Approved: [e-mail]

Submitted-by: [e-mail] (Robert Davies)
Posting-number: Volume 34, Issue 111
Archive-name: newmat07/part05
Environment: C++
Supersedes: newmat06: Volume 34, Issue 7-13

#! /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 archive 5 (of 8)."
# Contents:  newmat5.cxx newmat6.cxx newmat7.cxx newmat8.cxx
#   newmat9.cxx newmatrm.cxx
# Wrapped by robert@kea on Sun Jan 10 23:58:07 1993
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'\" \(11567 characters\)
sed "s/^X//" >'newmat5.cxx' <<'END_OF_FILE'
X//$$ newmat5.cxx         Transpose, evaluate etc
X
X// Copyright (C) 1991,2,3: R B Davies
X
X#include "include.h"
X
X#include "newmat.h"
X#include "newmatrc.h"
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(TransposedMatrix* tm, MatrixType mt)
X{
X   GeneralMatrix* gm1;
X
X   if (Compare(Type().t(),mt))
X   {
X      REPORT
X      gm1 = mt.New(ncols,nrows,tm);
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      gm1 = mt.New(ncols,nrows,tm);
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(); gm1->ReleaseAndDelete(); return gm1;
X}
X
XGeneralMatrix* SymmetricMatrix::Transpose(TransposedMatrix*, MatrixType mt)
X{ REPORT  return Evaluate(mt); }
X
X
XGeneralMatrix* DiagonalMatrix::Transpose(TransposedMatrix*, MatrixType mt)
X{ REPORT return Evaluate(mt); }
X
XBoolean 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(TransposedMatrix*, 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(TransposedMatrix*, 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 (Compare(this->Type(),mt)) { REPORT return this; }
X   REPORT
X   GeneralMatrix* gmx = mt.New(nrows,ncols,this);
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 (Compare(cgm->Type(),mt))
X   {
X      REPORT
X#ifdef TEMPS_DESTROYED_QUICKLY
X      GeneralMatrix* gmx = (GeneralMatrix*)cgm; delete this; return gmx;
X#else
X      return (GeneralMatrix*)cgm;
X#endif
X   }
X   REPORT
X   GeneralMatrix* gmx = cgm->Type().New(cgm->Nrows(),cgm->Ncols(),this);
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();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   delete this;
X#endif
X   return gmx; // no tDelete
X}
X
XGeneralMatrix* ShiftedMatrix::Evaluate(MatrixType mt)
X{
X   gm=((BaseMatrix*&)bm)->Evaluate();
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   Compare(gm->Type().AddEqualEl(),mt);
X   if (!(mt==gm->Type()))
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc,this);
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();
X#ifdef TEMPS_DESTROYED_QUICKLY
X      delete this;
X#endif
X      return gmx;
X   }
X   else if (gm->reuse())
X   {
X      REPORT gm->Add(f);
X#ifdef TEMPS_DESTROYED_QUICKLY
X      GeneralMatrix* gmx = gm; delete this; return gmx;
X#else
X      return gm;
X#endif
X   }
X   else
X   {
X      REPORT GeneralMatrix* gmy = gm->Type().New(nr,nc,this);
X      gmy->ReleaseAndDelete(); gmy->Add(gm,f);
X#ifdef TEMPS_DESTROYED_QUICKLY
X      delete this;
X#endif
X      return gmy;
X   }
X}
X
XGeneralMatrix* ScaledMatrix::Evaluate(MatrixType mt)
X{
X   gm=((BaseMatrix*&)bm)->Evaluate();
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   if (Compare(gm->Type(),mt))
X   {
X      if (gm->reuse())
X      {
X         REPORT gm->Multiply(f);
X#ifdef TEMPS_DESTROYED_QUICKLY
X         GeneralMatrix* gmx = gm; delete this; return gmx;
X#else
X         return gm;
X#endif
X      }
X      else
X      {
X         REPORT GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
X         gmx->ReleaseAndDelete(); gmx->Multiply(gm,f);
X#ifdef TEMPS_DESTROYED_QUICKLY
X         delete this;
X#endif
X         return gmx;
X      }
X   }
X   else
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc,this);
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();
X#ifdef TEMPS_DESTROYED_QUICKLY
X      delete this;
X#endif
X      return gmx;
X   }
X}
X
XGeneralMatrix* NegatedMatrix::Evaluate(MatrixType mt)
X{
X   gm=((BaseMatrix*&)bm)->Evaluate();
X   int nr=gm->Nrows(); int nc=gm->Ncols();
X   if (Compare(gm->Type(),mt))
X   {
X      if (gm->reuse())
X      {
X         REPORT gm->Negate();
X#ifdef TEMPS_DESTROYED_QUICKLY
X         GeneralMatrix* gmx = gm; delete this; return gmx;
X#else
X         return gm;
X#endif
X      }
X      else
X      {
X         REPORT
X         GeneralMatrix* gmx = gm->Type().New(nr,nc,this);
X         gmx->ReleaseAndDelete(); gmx->Negate(gm);
X#ifdef TEMPS_DESTROYED_QUICKLY
X         delete this;
X#endif
X         return gmx;
X      }
X   }
X   else
X   {
X      REPORT
X      GeneralMatrix* gmx = mt.New(nr,nc,this);
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();
X#ifdef TEMPS_DESTROYED_QUICKLY
X      delete this;
X#endif
X      return gmx;
X   }
X}
X
XGeneralMatrix* TransposedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   gm=((BaseMatrix*&)bm)->Evaluate();
X   Compare(gm->Type().t(),mt);
X   GeneralMatrix* gmx=gm->Transpose(this, mt);
X#ifdef TEMPS_DESTROYED_QUICKLY
X   delete this;
X#endif
X   return gmx;
X}
X
XGeneralMatrix* RowedMatrix::Evaluate(MatrixType mt)
X{
X   gm = ((BaseMatrix*&)bm)->Evaluate();
X   GeneralMatrix* gmx = new RowVector; MatrixErrorNoSpace(gmx);
X   gmx->nrows = 1; gmx->ncols = gmx->storage = gm->storage;
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
X#else
X   return gm->BorrowStore(gmx,mt);
X#endif
X}
X
XGeneralMatrix* ColedMatrix::Evaluate(MatrixType mt)
X{
X   gm = ((BaseMatrix*&)bm)->Evaluate();
X   GeneralMatrix* gmx = new ColumnVector; MatrixErrorNoSpace(gmx);
X   gmx->ncols = 1; gmx->nrows = gmx->storage = gm->storage;
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
X#else
X   return gm->BorrowStore(gmx,mt);
X#endif
X}
X
XGeneralMatrix* DiagedMatrix::Evaluate(MatrixType mt)
X{
X   gm = ((BaseMatrix*&)bm)->Evaluate();
X   GeneralMatrix* gmx = new DiagonalMatrix; MatrixErrorNoSpace(gmx);
X   gmx->nrows = gmx->ncols = gmx->storage = gm->storage;
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
X#else
X   return gm->BorrowStore(gmx,mt);
X#endif
X}
X
XGeneralMatrix* MatedMatrix::Evaluate(MatrixType mt)
X{
X   Tracer tr("MatedMatrix::Evaluate");
X   gm = ((BaseMatrix*&)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      Throw(IncompatibleDimensionsException());
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmy = gm; delete this; return gmy->BorrowStore(gmx,mt);
X#else
X   return gm->BorrowStore(gmx,mt);
X#endif
X}
X
XGeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   Tracer tr("SubMatrix(evaluate)");
X   gm = ((BaseMatrix*&)bm)->Evaluate();
X   if (row_number < 0) row_number = gm->Nrows();
X   if (col_number < 0) col_number = gm->Ncols();
X   if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
X      Throw(SubMatrixDimensionException());
X   if (IsSym) Compare(gm->Type().ssub(), mt);
X   else Compare(gm->Type().sub(), mt);
X   GeneralMatrix* gmx = mt.New(row_number, col_number, this);
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();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   delete this;
X#endif
X   return gmx;
X}
X
X
XGeneralMatrix* ReturnMatrixX::Evaluate(MatrixType mt)
X{
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx = gm; delete this; return gmx->Evaluate(mt);
X#else
X   return gm->Evaluate(mt);
X#endif
X}
X
X
Xvoid GeneralMatrix::Add(GeneralMatrix* gm1, Real f)
X{
X   REPORT
X   Real* s1=gm1->store; Real* s=store; int i=(storage >> 2);
X   while (i--)
X   { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; }
X   i = storage & 3; while (i--) *s++ = *s1++ + f;
X}
X
Xvoid GeneralMatrix::Add(Real f)
X{
X   REPORT
X   Real* s=store; int i=(storage >> 2);
X   while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; }
X   i = storage & 3; 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; int i=(storage >> 2);
X   while (i--)
X   { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); }
X   i = storage & 3; while(i--) *s++ = -(*s1++);
X}
X
Xvoid GeneralMatrix::Negate()
X{
X   REPORT
X   Real* s=store; int i=(storage >> 2);
X   while (i--)
X   { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; }
X   i = storage & 3; while(i--) { *s = -(*s); s++; }
X}
X
Xvoid GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f)
X{
X   REPORT
X   Real* s1=gm1->store; Real* s=store;  int i=(storage >> 2);
X   while (i--)
X   { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; }
X   i = storage & 3; while (i--) *s++ = *s1++ * f;
X}
X
Xvoid GeneralMatrix::Multiply(Real f)
X{
X   REPORT
X   Real* s=store; int i=(storage >> 2);
X   while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; }
X   i = storage & 3; while (i--) *s++ *= f;
X}
X
X
X/************************ MatrixInput routines ****************************/
X
Xint MatrixInput::n = 0;       // number values still to be read
XReal* MatrixInput::r;         // pointer to last thing read
Xint MatrixInput::depth = 0;   // number of objects of this class in existence
X
XMatrixInput MatrixInput::operator<<(Real f)
X{
X   if (!(n--))
X   { depth=0;  n=0; Throw(ProgramException("List of values too long")); }
X   *(++r) = f;
X   return MatrixInput();
X}
X
X
XMatrixInput BandMatrix::operator<<(Real)
X{
X   Throw(ProgramException("Cannot use list read with a BandMatrix"));
X   return MatrixInput();
X}
X
Xvoid BandMatrix::operator<<(const Real*)
X{ Throw(ProgramException("Cannot use array read with a BandMatrix")); }
X
XMatrixInput GeneralMatrix::operator<<(Real f)
X{
X   if (MatrixInput::n)
X   {
X      MatrixInput::depth=0;            // so we can recover
X      MatrixInput::n=0;                // so we can recover
X      Throw(ProgramException("A list of values was too short"));
X   }
X   MatrixInput::n = Storage();
X   if (MatrixInput::n<=0)
X      Throw(ProgramException("Loading data to zero dimension matrix"));
X   MatrixInput::r = Store(); *(MatrixInput::r) = f; MatrixInput::n--;
X   return MatrixInput();
X}
X
XMatrixInput::~MatrixInput()
X{
X   if (depth && --depth == 0 && n != 0)
X   {
X      depth = 0; n = 0;                // so we can recover
X      Throw(ProgramException("A list of values was too short"));
X   }
X}
X
END_OF_FILE
if test 11567 -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'\" \(13261 characters\)
sed "s/^X//" >'newmat6.cxx' <<'END_OF_FILE'
X//$$ newmat6.cxx            Operators, element access, submatrices
X
X// Copyright (C) 1991,2,3: R B Davies
X
X#include "include.h"
X
X#include "newmat.h"
X#include "newmatrc.h"
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
X   return store[n-1];
X}
X
XReal& DiagonalMatrix::operator()(int m)
X{
X   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
X   return store[m-1];
X}
X
XReal& ColumnVector::operator()(int m)
X{
X   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
X   return store[m-1];
X}
X
XReal& RowVector::operator()(int n)
X{
X   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
X   return store[n-1];
X}
X
XReal& BandMatrix::operator()(int m, int n)
X{
X   int w = upper+lower+1; int i = lower+n-m;
X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this));
X   return store[w*(m-1)+i];
X}
X
XReal& UpperBandMatrix::operator()(int m, int n)
X{
X   int w = upper+1; int i = n-m;
X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this));
X   return store[w*(m-1)+i];
X}
X
XReal& LowerBandMatrix::operator()(int m, int n)
X{
X   int w = lower+1; int i = lower+n-m;
X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this));
X   return store[w*(m-1)+i];
X}
X
XReal& SymmetricBandMatrix::operator()(int m, int n)
X{
X   int w = lower+1;
X   if (m>=n)
X   {
X      int i = lower+n-m;
X      if ( m>nrows || n<=0 || i<0 )
X         Throw(IndexException(m,n,*this));
X      return store[w*(m-1)+i];
X   }
X   else
X   {
X      int i = lower+m-n;
X      if ( n>nrows || m<=0 || i<0 )
X         Throw(IndexException(m,n,*this));
X      return store[w*(n-1)+i];
X   }
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
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)
X      Throw(IndexException(m,n,*this));
X   return store[n-1];
X}
X
XReal DiagonalMatrix::operator()(int m) const
X{
X   if (m<=0 || m>nrows) Throw(IndexException(m,*this));
X   return store[m-1];
X}
X
XReal ColumnVector::operator()(int m) const
X{
X   if (m<=0 || m> nrows) Throw(IndexException(m,*this));
X   return store[m-1];
X}
X
XReal RowVector::operator()(int n) const
X{
X   if (n<=0 || n> ncols) Throw(IndexException(n,*this));
X   return store[n-1];
X}
X
XReal BandMatrix::operator()(int m, int n) const
X{
X   int w = upper+lower+1; int i = lower+n-m;
X   if (m<=0 || m>nrows || n<=0 || n>ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this));
X   return store[w*(m-1)+i];
X}
X
XReal SymmetricBandMatrix::operator()(int m, int n) const
X{
X   int w = lower+1;
X   if (m>=n)
X   {
X      int i = lower+n-m;
X      if ( m>nrows || n<=0 || i<0 )
X         Throw(IndexException(m,n,*this));
X      return store[w*(m-1)+i];
X   }
X   else
X   {
X      int i = lower+m-n;
X      if ( n>nrows || m<=0 || i<0 )
X         Throw(IndexException(m,n,*this));
X      return store[w*(n-1)+i];
X   }
X}
X
X#endif
X
XReal BaseMatrix::AsScalar() const
X{
X   REPORT
X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   if (gm->nrows!=1 || gm->ncols!=1)
X   {
X      Try
X      {
X         Tracer tr("AsScalar");
X     Throw(ProgramException("Cannot convert to scalar", *gm));
X      }
X      CatchAll { gm->tDelete(); Throw(); }
X   }
X   Real x = *(gm->store); gm->tDelete(); return x;
X}
X
X#ifdef TEMPS_DESTROYED_QUICKLY
X
XAddedMatrix& BaseMatrix::operator+(const BaseMatrix& bm) const
X{
X   REPORT
X   AddedMatrix* x = new AddedMatrix(this, &bm);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XMultipliedMatrix& BaseMatrix::operator*(const BaseMatrix& bm) const
X{
X   REPORT
X   MultipliedMatrix* x = new MultipliedMatrix(this, &bm);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
X//SolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx) const
XSolvedMatrix& InvertedMatrix::operator*(const BaseMatrix& bmx)
X{
X   REPORT
X   SolvedMatrix* x;
X   Try { x = new SolvedMatrix(bm, &bmx); MatrixErrorNoSpace(x); }
X   CatchAll { delete this; Throw(); }
X   delete this;                // since we are using bm rather than this
X   return *x;
X}
X
XSubtractedMatrix& BaseMatrix::operator-(const BaseMatrix& bm) const
X{
X   REPORT
X   SubtractedMatrix* x = new SubtractedMatrix(this, &bm);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XShiftedMatrix& BaseMatrix::operator+(Real f) const
X{
X   REPORT
X   ShiftedMatrix* x = new ShiftedMatrix(this, f);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XScaledMatrix& BaseMatrix::operator*(Real f) const
X{
X   REPORT
X   ScaledMatrix* x = new ScaledMatrix(this, f);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XScaledMatrix& BaseMatrix::operator/(Real f) const
X{
X   REPORT
X   ScaledMatrix* x = new ScaledMatrix(this, 1.0/f);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XShiftedMatrix& BaseMatrix::operator-(Real f) const
X{
X   REPORT
X   ShiftedMatrix* x = new ShiftedMatrix(this, -f);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XTransposedMatrix& BaseMatrix::t() const
X{
X   REPORT
X   TransposedMatrix* x = new TransposedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XNegatedMatrix& BaseMatrix::operator-() const
X{
X   REPORT
X   NegatedMatrix* x = new NegatedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XInvertedMatrix& BaseMatrix::i() const
X{
X   REPORT
X   InvertedMatrix* x = new InvertedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XConstMatrix& GeneralMatrix::c() const
X{
X   if (tag != -1)
X      Throw(ProgramException(".c() applied to temporary matrix"));
X   REPORT
X   ConstMatrix* x = new ConstMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XRowedMatrix& BaseMatrix::AsRow() const
X{
X   REPORT
X   RowedMatrix* x = new RowedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XColedMatrix& BaseMatrix::AsColumn() const
X{
X   REPORT
X   ColedMatrix* x = new ColedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XDiagedMatrix& BaseMatrix::AsDiagonal() const
X{
X   REPORT
X   DiagedMatrix* x = new DiagedMatrix(this);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
XMatedMatrix& BaseMatrix::AsMatrix(int nrx, int ncx) const
X{
X   REPORT
X   MatedMatrix* x = new MatedMatrix(this,nrx,ncx);
X   MatrixErrorNoSpace(x);
X   return *x;
X}
X
X#else
X
XAddedMatrix BaseMatrix::operator+(const BaseMatrix& bm) const
X{ REPORT return AddedMatrix(this, &bm); }
X
XMultipliedMatrix BaseMatrix::operator*(const BaseMatrix& bm) const
X{ REPORT return MultipliedMatrix(this, &bm); }
X
XSolvedMatrix InvertedMatrix::operator*(const BaseMatrix& bmx) const
X{ REPORT return SolvedMatrix(bm, &bmx); }
X
XSubtractedMatrix BaseMatrix::operator-(const BaseMatrix& bm) const
X{ REPORT return SubtractedMatrix(this, &bm); }
X
XShiftedMatrix BaseMatrix::operator+(Real f) const
X{ REPORT return ShiftedMatrix(this, f); }
X
XScaledMatrix BaseMatrix::operator*(Real f) const
X{ REPORT return ScaledMatrix(this, f); }
X
XScaledMatrix BaseMatrix::operator/(Real f) const
X{ REPORT return ScaledMatrix(this, 1.0/f); }
X
XShiftedMatrix BaseMatrix::operator-(Real f) const
X{ REPORT return ShiftedMatrix(this, -f); }
X
XTransposedMatrix BaseMatrix::t() const
X{ REPORT return TransposedMatrix(this); }
X
XNegatedMatrix BaseMatrix::operator-() const
X{ REPORT return NegatedMatrix(this); }
X
XInvertedMatrix BaseMatrix::i() const
X{ REPORT return InvertedMatrix(this); }
X
XConstMatrix GeneralMatrix::c() const
X{
X   if (tag != -1)
X      Throw(ProgramException(".c() applied to temporary matrix"));
X   REPORT return ConstMatrix(this);
X}
X
XRowedMatrix BaseMatrix::AsRow() const
X{ REPORT return RowedMatrix(this); }
X
XColedMatrix BaseMatrix::AsColumn() const
X{ REPORT return ColedMatrix(this); }
X
XDiagedMatrix BaseMatrix::AsDiagonal() const
X{ REPORT return DiagedMatrix(this); }
X
XMatedMatrix BaseMatrix::AsMatrix(int nrx, int ncx) const
X{ REPORT return MatedMatrix(this,nrx,ncx); }
X
X#endif
X
Xvoid GeneralMatrix::operator=(Real f)
X{ REPORT int i=storage; Real* s=store; while (i--) { *s++ = f; } }
X
Xvoid Matrix::operator=(const BaseMatrix& X)
X{
X   REPORT //CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::Rt);
X}
X
Xvoid RowVector::operator=(const BaseMatrix& X)
X{
X   REPORT  // CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::RV);
X   if (nrows!=1)
X   {
X      Tracer tr("RowVector(=)");
X      Throw(VectorException(*this));
X   }
X}
X
Xvoid ColumnVector::operator=(const BaseMatrix& X)
X{
X   REPORT //CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::CV);
X   if (ncols!=1)
X   {
X      Tracer tr("ColumnVector(=)");
X      Throw(VectorException(*this));
X   }
X}
X
Xvoid SymmetricMatrix::operator=(const BaseMatrix& X)
X{
X   REPORT // CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::Sm);
X}
X
Xvoid UpperTriangularMatrix::operator=(const BaseMatrix& X)
X{
X   REPORT //CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::UT);
X}
X
Xvoid LowerTriangularMatrix::operator=(const BaseMatrix& X)
X{
X   REPORT //CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::LT);
X}
X
Xvoid DiagonalMatrix::operator=(const BaseMatrix& X)
X{
X   REPORT // CheckConversion(X);
X   MatrixConversionCheck mcc;
X   Eq(X,MatrixType::Dg);
X}
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)
X      Throw(IndexException(m,n,*this,TRUE));
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)
X      Throw(IndexException(m,n,*this,TRUE));
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)
X      Throw(IndexException(m,n,*this,TRUE));
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)
X      Throw(IndexException(m,n,*this,TRUE));
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)
X      Throw(IndexException(m,n,*this,TRUE));
X   return store[n];
X}
X
XReal& DiagonalMatrix::element(int m)
X{
X   if (m<0 || m>=nrows) Throw(IndexException(m,*this,TRUE));
X   return store[m];
X}
X
XReal& ColumnVector::element(int m)
X{
X   if (m<0 || m>= nrows) Throw(IndexException(m,*this,TRUE));
X   return store[m];
X}
X
XReal& RowVector::element(int n)
X{
X   if (n<0 || n>= ncols)  Throw(IndexException(n,*this,TRUE));
X   return store[n];
X}
X
XReal& BandMatrix::element(int m, int n)
X{
X   int w = upper+lower+1; int i = lower+n-m;
X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this,TRUE));
X   return store[w*m+i];
X}
X
XReal& UpperBandMatrix::element(int m, int n)
X{
X   int w = upper+1; int i = n-m;
X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this,TRUE));
X   return store[w*m+i];
X}
X
XReal& LowerBandMatrix::element(int m, int n)
X{
X   int w = lower+1; int i = lower+n-m;
X   if (m<0 || m>= nrows || n<0 || n>= ncols || i<0 || i>=w)
X      Throw(IndexException(m,n,*this,TRUE));
X   return store[w*m+i];
X}
X
XReal& SymmetricBandMatrix::element(int m, int n)
X{
X   int w = lower+1;
X   if (m>=n)
X   {
X      int i = lower+n-m;
X      if ( m>=nrows || n<0 || i<0 )
X         Throw(IndexException(m,n,*this,TRUE));
X      return store[w*m+i];
X   }
X   else
X   {
X      int i = lower+m-n;
X      if ( n>=nrows || m<0 || i<0 )
X         Throw(IndexException(m,n,*this,TRUE));
X      return store[w*n+i];
X   }
X}
X
END_OF_FILE
if test 13261 -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'\" \(16089 characters\)
sed "s/^X//" >'newmat7.cxx' <<'END_OF_FILE'
X//$$ newmat7.cxx     Invert, solve, binary operations
X
X// Copyright (C) 1991,2,3: R B Davies
X
X#include "include.h"
X
X#include "newmat.h"
X#include "newmatrc.h"
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 GeneralMatrix*
X   GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
Xstatic GeneralMatrix*
X   GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
Xstatic GeneralMatrix*
X   GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
Xstatic GeneralMatrix*
X   GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
X
XGeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   gm1=((BaseMatrix*&)bm1)->Evaluate();
X   gm2=((BaseMatrix*&)bm2)->Evaluate();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx;
X   Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
X   CatchAll { delete this; Throw(); }
X   delete this; return gmx;
X#else
X   return GeneralAdd(gm1,gm2,this,mt);
X#endif
X}
X
XGeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   gm1=((BaseMatrix*&)bm1)->Evaluate();
X   gm2=((BaseMatrix*&)bm2)->Evaluate();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx;
X   Try { gmx = GeneralSub(gm1,gm2,this,mt); }
X   CatchAll { delete this; Throw(); }
X   delete this; return gmx;
X#else
X   return GeneralSub(gm1,gm2,this,mt);
X#endif
X}
X
XGeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   gm2 = ((BaseMatrix*&)bm2)->Evaluate();
X   gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
X   gm1=((BaseMatrix*&)bm1)->Evaluate();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx;
X   Try { gmx = GeneralMult(gm1, gm2, this, mt); }
X   CatchAll { delete this; Throw(); }
X   delete this; return gmx;
X#else
X   return GeneralMult(gm1, gm2, this, mt);
X#endif
X}
X
XGeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
X{
X   REPORT
X   gm1=((BaseMatrix*&)bm1)->Evaluate();
X   gm2=((BaseMatrix*&)bm2)->Evaluate();
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx;
X   Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
X   CatchAll { delete this; Throw(); }
X   delete this; return gmx;
X#else
X   return GeneralSolv(gm1,gm2,this,mt);
X#endif
X}
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   Real* s1=gm1->Store(); Real* s2=gm2->Store();
X   Real* s=gm->Store(); int i=gm->Storage() >> 2;
X   while (i--)
X   {
X       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
X       *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
X   }
X   i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
X}
X
Xstatic void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
X   while (i--)
X   { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
X   i=gm->Storage() & 3; while (i--) *s++ += *s2++;
X}
X
Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   REPORT
X   Real* s1=gm1->Store(); Real* s2=gm2->Store();
X   Real* s=gm->Store(); int i=gm->Storage() >> 2;
X   while (i--)
X   {
X       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
X       *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
X   }
X   i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
X}
X
Xstatic void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
X   while (i--)
X   { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
X   i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
X}
X
Xstatic void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   REPORT
X   Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
X   while (i--)
X   {
X      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
X      *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
X   }
X   i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
X}
X
X// routines for adding or subtracting matrices of different storage structure
X
Xstatic void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   int nr = gm->Nrows();
X   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
X   MatrixRow mr(gm, StoreOnExit+DirectPart);
X   while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
X}
X
Xstatic void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
X// Add into first argument
X{
X   int nr = gm->Nrows();
X   MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
X   MatrixRow mr2(gm2, LoadOnEntry);
X   while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
X}
X
Xstatic void SubtractDS
X   (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   int nr = gm->Nrows();
X   MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
X   MatrixRow mr(gm, StoreOnExit+DirectPart);
X   while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
X}
X
Xstatic void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   int nr = gm->Nrows();
X   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
X   MatrixRow mr2(gm2, LoadOnEntry);
X   while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
X}
X
Xstatic void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
X{
X   int nr = gm->Nrows();
X   MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
X   MatrixRow mr2(gm2, LoadOnEntry);
X   while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
X}
X
X#ifdef __GNUG__
Xvoid AddedMatrix::SelectVersion
X   (MatrixType mtx, int& c1, int& c2) const
X#else
Xvoid AddedMatrix::SelectVersion
X   (MatrixType mtx, Boolean& c1, Boolean& c2) const
X#endif
X// for determining version of add and subtract routines
X// will need to modify if further matrix structures are introduced
X{
X   MatrixBandWidth bm1 = gm1->BandWidth();
X   MatrixBandWidth bm2 = gm2->BandWidth();
X   MatrixBandWidth bmx = bm1 + bm2;
X   c1 = (mtx == gm1->Type()) && (bmx == bm1);
X   c2 = (mtx == gm2->Type()) && (bmx == bm2);
X}
X
Xstatic GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   AddedMatrix* am, MatrixType mtx)
X{
X   Tracer tr("GeneralAdd");
X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
X      Throw(IncompatibleDimensionsException());
X   Compare(gm1->Type() + gm2->Type(),mtx);
X#ifdef __GNUG__
X   int c1,c2; am->SelectVersion(mtx,c1,c2);
X#else
X   Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
X#endif
X   if (c1 && c2)
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 GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
X         gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
X      }
X   }
X   else
X   {
X      if (c1 && gm1->reuse() )               // must have type test first
X      { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
X      else if (c2 && gm2->reuse() )
X      { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
X      else
X      {
X         REPORT
X     GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
X     if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
X         gmx->ReleaseAndDelete(); return gmx;
X      }
X   }
X}
X
X
Xstatic GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   SubtractedMatrix* sm, MatrixType mtx)
X{
X   Tracer tr("GeneralSub");
X   Compare(gm1->Type() + gm2->Type(),mtx);
X   int nr=gm1->Nrows(); int nc=gm1->Ncols();
X   if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
X      Throw(IncompatibleDimensionsException());
X#ifdef __GNUG__
X   int c1,c2; sm->SelectVersion(mtx,c1,c2);
X#else
X   Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
X#endif
X   if (c1 && c2)
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,sm);
X         gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
X      }
X   }
X   else
X   {
X      if ( c1 && gm1->reuse() )
X      { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
X      else if ( c2 && gm2->reuse() )
X      {
X         REPORT
X         ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
X      }
X      else
X      {
X         REPORT
X     GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
X     if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
X     gmx->ReleaseAndDelete(); return gmx;
X      }
X   }
X}
X
Xstatic GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MultipliedMatrix* mm, MatrixType mtx)
X{
X   REPORT
X   Tracer tr("GeneralMult1");
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows())
X      Throw(IncompatibleDimensionsException());
X   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
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
Xstatic GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MultipliedMatrix* mm, MatrixType mtx)
X{
X   // version that accesses by row only - not good for thin matrices
X   // or column vectors in right hand term. Needs fixing
X   REPORT
X   Tracer tr("GeneralMult2");
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows())
X      Throw(IncompatibleDimensionsException());
X   GeneralMatrix* gmx = mtx.New(nr,nc,mm);
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* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
X{
X   // matrix multiplication for type Matrix only
X   REPORT
X   Tracer tr("MatrixMult");
X
X   int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
X   if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
X
X   Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
X
X   Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
X
X   if (ncr)
X   {
X      while (nr--)
X      {
X         Real* s2x = s2; int j = ncr;
X         Real* sx = s; Real f = *s1++; int k = nc;
X         while (k--) *sx++ = f * *s2x++;
X         while (--j)
X            { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
X         s = sx;
X      }
X   }
X   else *gm = 0.0;
X
X   gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
X}
X
Xstatic GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   MultipliedMatrix* mm, MatrixType mtx)
X{
X   if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
X   else
X   {
X      Compare(gm1->Type() * gm2->Type(),mtx);
X      int nr = gm2->Nrows(); int nc = gm2->Ncols();
X      if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
X      else return GeneralMult2(gm1, gm2, mm, mtx);
X   }
X}
X
Xstatic GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
X   BaseMatrix* sm, MatrixType mtx)
X{
X   REPORT
X   Tracer tr("GeneralSolv");
X   Compare(gm1->Type().i() * gm2->Type(),mtx);
X   int nr=gm1->Nrows(); int nc=gm2->Ncols();
X   if (gm1->Ncols() !=gm2->Nrows())
X      Throw(IncompatibleDimensionsException());
X   GeneralMatrix* gmx = mtx.New(nr,nc,sm);
X
X   Real* r = new Real [nr]; MatrixErrorNoSpace(r);
X#ifndef ATandT
X   MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
X                           // deleted for ATandT, to balance deletion below
X#endif
X   MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
X   MatrixCol mc2(gm2, r, LoadOnEntry);
X   GeneralMatrix* gms = gm1->MakeSolver();
X   Try
X   {
X      while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
X   }
X   CatchAll
X   {
X      gms->tDelete(); delete gmx; gm2->tDelete();
X#ifndef ATandT
X      MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
X                          // ATandT version 2.1 gives an internal error
X#endif
X#ifdef Version21
X      delete [] r;
X#else
X      delete [nr] r;
X#endif
X      Throw();
X   }
X   gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
X#ifndef ATandT
X   MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
X                          // ATandT version 2.1 gives an internal error
X#endif
X#ifdef Version21
X   delete [] r;
X#else
X   delete [nr] r;
X#endif
X   return gmx;
X}
X
XGeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
X{
X   // Matrix Inversion - use solve routines
X   REPORT
X   gm=((BaseMatrix*&)bm)->Evaluate();
X   int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
X#ifdef TEMPS_DESTROYED_QUICKLY
X   GeneralMatrix* gmx;
X   Try { gmx = GeneralSolv(gm,&I,this,mtx); }
X   CatchAll { delete this; Throw(); }
X   delete this; return gmx;
X#else
X   return GeneralSolv(gm,&I,this,mtx);
X#endif
X}
X
X
X/*************************** norm functions ****************************/
X
XReal BaseMatrix::Norm1() const
X{
X   // maximum of sum of absolute values of a column
X   REPORT
X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   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() const
X{
X   // maximum of sum of absolute values of a row
X   REPORT
X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   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 16089 -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'\" \(8288 characters\)
sed "s/^X//" >'newmat8.cxx' <<'END_OF_FILE'
X//$$ newmat8.cxx         Advanced LU transform, scalar functions
X
X// Copyright (C) 1991,2,3: R B Davies
X
X#define WANT_MATH
X
X#include "include.h"
X
X#include "newmatap.h"
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   Tracer trace("Crout(ludcmp)");
X   int i,j;
X
X   Real* vv=new Real [nrows]; MatrixErrorNoSpace(vv);
X   MONITOR_REAL_NEW("Make  (CroutMat)",nrows,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)
X      {
X         MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
X#ifdef Version21
X         sing = TRUE; delete [] vv; return;
X#else
X         sing = TRUE; delete [nrows] vv; return;
X#endif
X      }
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)
X      {
X         MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
X#ifdef Version21
X         sing = TRUE; delete [] vv; return;
X#else
X         sing = TRUE; delete [nrows] vv; return;
X#endif
X      }
X      Real dum=1.0/(*ai);
X      i=nrows-j; while (--i) { ai += nrows; *ai *= dum; }
X
X      aj++;
X   }
X   MONITOR_REAL_DELETE("Delete (CroutMt)",nrows,vv)
X#ifdef Version21
X   delete [] vv;
X#else
X   delete [nrows] vv;
X#endif
X}
X
Xvoid CroutMatrix::lubksb(Real* B, int mini)
X{
X   REPORT
X   Tracer trace("Crout(lubksb)");
X   if (sing) Throw(SingularException(*this));
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() const
X{
X   REPORT
X   Real sum = 0.0; int i = storage; Real* s = store;
X   while (i--) sum += square(*s++);
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal GeneralMatrix::SumAbsoluteValue() const
X{
X   REPORT
X   Real sum = 0.0; int i = storage; Real* s = store;
X   while (i--) sum += fabs(*s++);
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal GeneralMatrix::MaximumAbsoluteValue() const
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   ((GeneralMatrix&)*this).tDelete(); return maxval;
X}
X
XReal SymmetricMatrix::SumSquare() const
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   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
X}
X
XReal SymmetricMatrix::SumAbsoluteValue() const
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   ((GeneralMatrix&)*this).tDelete(); return sum1 + 2.0 * sum2;
X}
X
XReal BaseMatrix::SumSquare() const
X{
X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   Real s = gm->SumSquare(); return s;
X}
X
XReal BaseMatrix::SumAbsoluteValue() const
X{
X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   Real s = gm->SumAbsoluteValue(); return s;
X}
X
XReal BaseMatrix::MaximumAbsoluteValue() const
X{
X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   Real s = gm->MaximumAbsoluteValue(); return s;
X}
X
XReal Matrix::Trace() const
X{
X   REPORT
X   Tracer trace("Trace");
X   int i = nrows; int d = i+1;
X   if (i != ncols) Throw(NotSquareException(*this));
X   Real sum = 0.0; Real* s = store;
X   while (i--) { sum += *s; s += d; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal DiagonalMatrix::Trace() const
X{
X   REPORT
X   int i = nrows; Real sum = 0.0; Real* s = store;
X   while (i--) sum += *s++;
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal SymmetricMatrix::Trace() const
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   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal LowerTriangularMatrix::Trace() const
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   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal UpperTriangularMatrix::Trace() const
X{
X   REPORT
X   int i = nrows; Real sum = 0.0; Real* s = store;
X   while (i) { sum += *s; s += i--; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal BandMatrix::Trace() const
X{
X   REPORT
X   int i = nrows; int w = lower+upper+1;
X   Real sum = 0.0; Real* s = store+lower;
X   while (i--) { sum += *s; s += w; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal SymmetricBandMatrix::Trace() const
X{
X   REPORT
X   int i = nrows; int w = lower+1;
X   Real sum = 0.0; Real* s = store+lower;
X   while (i--) { sum += *s; s += w; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XReal BaseMatrix::Trace() const
X{
X   REPORT
X   GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate(MatrixType::Dg);
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() const { return sign * exp(log_value); }
X
XLogAndSign::LogAndSign(Real f)
X{
X   if (f == 0.0) { log_value = 0.0; sign = 0; return; }
X   else if (f < 0.0) { sign = -1; f = -f; }
X   else sign = 1;
X   log_value = log(f);
X}
X
XLogAndSign DiagonalMatrix::LogDeterminant() const
X{
X   REPORT
X   int i = nrows; LogAndSign sum; Real* s = store;
X   while (i--) sum *= *s++;
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XLogAndSign LowerTriangularMatrix::LogDeterminant() const
X{
X   REPORT
X   int i = nrows; LogAndSign sum; Real* s = store; int j = 2;
X   while (i--) { sum *= *s; s += j++; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XLogAndSign UpperTriangularMatrix::LogDeterminant() const
X{
X   REPORT
X   int i = nrows; LogAndSign sum; Real* s = store;
X   while (i) { sum *= *s; s += i--; }
X   ((GeneralMatrix&)*this).tDelete(); return sum;
X}
X
XLogAndSign BaseMatrix::LogDeterminant() const
X{
X   REPORT GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
X   LogAndSign sum = gm->LogDeterminant(); return sum;
X}
X
XLogAndSign GeneralMatrix::LogDeterminant() const
X{
X   REPORT
X   Tracer tr("Determinant");
X   if (nrows != ncols) Throw(NotSquareException(*this));
X   CroutMatrix C(*this); return C.LogDeterminant();
X}
X
XLogAndSign CroutMatrix::LogDeterminant() const
X{
X   REPORT
X   if (sing) return 0.0;
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}
X
XLinearEquationSolver::LinearEquationSolver(const BaseMatrix& bm)
X: gm( ( ((BaseMatrix&)bm).Evaluate() )->MakeSolver() )
X{
X   if (gm==&bm) { REPORT  gm = gm->Image(); }
X   // want a copy if  *gm is actually bm
X   else { REPORT  gm->Protect(); }
X}
X
END_OF_FILE
if test 8288 -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'\" \(954 characters\)
sed "s/^X//" >'newmat9.cxx' <<'END_OF_FILE'
X//$$ newmat9.cxx         Input and output
X
X// Copyright (C) 1991,2,3: R B Davies
X
X
X#define WANT_STREAM
X
X#include "include.h"
X
X#include "newmat.h"
X#include "newmatrc.h"
X#include "newmatio.h"
X
X//#define REPORT { static ExeCounter ExeCount(__LINE__,9); ExeCount++; }
X
X#define REPORT {}
X
Xostream& operator<<(ostream& s, const BaseMatrix& X)
X{
X   GeneralMatrix* gm = ((BaseMatrix&)X).Evaluate(); operator<<(s, *gm);
X   gm->tDelete(); return s;
X}
X
X
Xostream& operator<<(ostream& s, const GeneralMatrix& X)
X{
X   MatrixRow mr((GeneralMatrix*)&X, LoadOnEntry);
X   int w = s.width();  int nr = X.Nrows();  long f = s.flags();
X   s.setf(ios::fixed, ios::floatfield);
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 954 -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'\" \(3388 characters\)
sed "s/^X//" >'newmatrm.cxx' <<'END_OF_FILE'
X//$$newmatrm.cxx                         rectangular matrix operations
X
X// Copyright (C) 1991,2,3: R B Davies
X
X#include "include.h"
X
X#include "newmat.h"
X#include "newmatrm.h"
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)
X   {
X      Tracer tr("newmatrm");
X      Throw(InternalException("Dimensions differ in *"));
X   }
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)
X   {
X      Tracer tr("newmatrm");
X      Throw(InternalException("Dimensions differ in AddScaled"));
X   }
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)
X   {
X      Tracer tr("newmatrm");
X      Throw(InternalException("Dimensions differ in Divide"));
X   }
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)
X   {
X      Tracer tr("newmatrm");
X      Throw(InternalException("Dimensions differ in ComplexScale"));
X   }
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)
X   {
X      Tracer tr("newmatrm");
X      Throw(InternalException("Dimensions differ in Rotate"));
X   }
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 3388 -ne `wc -c <'newmatrm.cxx'`; then
    echo shar: \"'newmatrm.cxx'\" unpacked with wrong size!
fi
# end of 'newmatrm.cxx'
fi
echo shar: End of archive 5 \(of 8\).
cp /dev/null ark5isdone
MISSING=""
for I in 1 2 3 4 5 6 7 8 ; do
    if test ! -f ark${I}isdone ; then
    MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 8 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0

exit 0 # Just in case...