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 '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; inrows = 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 || nncols) 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 || mnrows) 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 || nncols) 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 || mnrows) 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=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=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 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= 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=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; iSumSquare(); 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.