diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..d3b286e8 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,27 @@ +default_language_version: + python: python3 +ci: + autofix_commit_msg: | + [pre-commit.ci] auto fixes from pre-commit hooks + autofix_prs: true + autoupdate_branch: 'pre-commit-autoupdate' + autoupdate_commit_msg: '[pre-commit.ci] pre-commit autoupdate' + autoupdate_schedule: monthly + skip: [no-commit-to-branch] + submodules: false +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.6.0 + hooks: + - id: check-yaml + exclude: 'meta.yaml' + - id: end-of-file-fixer + - id: trailing-whitespace + exclude: '\.(rst|txt)$' + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.4.0 + hooks: + - id: no-commit-to-branch + name: Prevent Commit to Main Branch + args: ["--branch", "main"] + stages: [pre-commit] diff --git a/CHANGELOG.md b/CHANGELOG.md index b6c6ccb1..c0a09067 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,13 +3,13 @@ ## Version 2024.1.1, - 2024-02-06 ### Added -- Make sure Molecule::BuildConnectivityTable() always list all atoms, +- Make sure Molecule::BuildConnectivityTable() always list all atoms, even if their connectivity list is empty (pathological case) ## Version 2022.1.4, - 2022-12-03 ### Added -- Add UnitCell::ChangeSpaceGroup(), which updates lattice parameter symmetry constraints. +- Add UnitCell::ChangeSpaceGroup(), which updates lattice parameter symmetry constraints. ### Changed - Take into account spacegroup clock in UnitCell::GetLatticePar and InitMatrices @@ -20,8 +20,8 @@ ### Changed -- the list of HKL reflections will now be automatically be re-generated - for a PowderPatternDiffraction when the Crystal's spacegroup changes, +- the list of HKL reflections will now be automatically be re-generated + for a PowderPatternDiffraction when the Crystal's spacegroup changes, or the lattice parameters are modified by more than 0.5% - Fixed the powder pattern indexing test @@ -44,7 +44,7 @@ ### Changed -- Add relative_length_tolerance and absolute_angle_tolerance_degree to +- Add relative_length_tolerance and absolute_angle_tolerance_degree to SpaceGroupExplorer::Run() and RunAll() - Crystal::XMLInput(): add a hook to re-use atomic scattering power when mDeleteSubObjInDestructor is False @@ -72,13 +72,13 @@ ### Added - Move SpaceGroupExplorer in a separate class for non-GUI access. Allow keeping or not the tested spacegroup, or the best solution. - Sort solutions by the GoF multiplied by the ratio of the number + Sort solutions by the GoF multiplied by the ratio of the number of non-extinct reflections in the spacegoup relative to P1. - MonteCarloObj: add public access to AutoLSQ option - OptimizationObj: add access to Options by number or name - RefinableObj: provide access to options by name - Add STL-type methods (begin, end, size) for ObjRegistry and Molecule objects -- Add Crystal.GetFormula(). Use formula to automatically name Crystal and +- Add Crystal.GetFormula(). Use formula to automatically name Crystal and DiffractionDataSingleCrystal when imported from CIF and no name is given ### Changed diff --git a/src/ObjCryst/ObjCryst/CIF.cpp b/src/ObjCryst/ObjCryst/CIF.cpp index f17a2bec..8fd79cb9 100644 --- a/src/ObjCryst/ObjCryst/CIF.cpp +++ b/src/ObjCryst/ObjCryst/CIF.cpp @@ -1008,7 +1008,7 @@ Crystal* CreateCrystalFromCIF(CIF &cif,const bool verbose,const bool checkSymAsX (*fpObjCrystInformUser)("CIF: Opening CIF"); Chronometer chrono; chrono.start(); - + // If oneScatteringPowerPerElement==true, we hold this to compute the average Biso per element std::map > vElementBiso; @@ -1184,7 +1184,7 @@ Crystal* CreateCrystalFromCIF(CIF &cif,const bool verbose,const bool checkSymAsX // Try to set name from CIF. If that fails, the computed formula will be used at the end if(pos->second.mName!="") pCryst->SetName(pos->second.mName); else if(pos->second.mFormula!="") pCryst->SetName(pos->second.mFormula); - + const float t1=chrono.seconds(); (*fpObjCrystInformUser)((boost::format("CIF: Create Crystal:%s(%s)(dt=%6.3fs)")%pCryst->GetName() % pCryst->GetSpaceGroup().GetName() % t1).str()); diff --git a/src/ObjCryst/ObjCryst/Crystal.h b/src/ObjCryst/ObjCryst/Crystal.h index c60bd9af..a32e07c6 100644 --- a/src/ObjCryst/ObjCryst/Crystal.h +++ b/src/ObjCryst/ObjCryst/Crystal.h @@ -189,7 +189,7 @@ class Crystal:public UnitCell /// \todo one function to print on one line and a PrintLong() function /// \param os the stream to which the information is outputed (default=cout) void Print(ostream &os=cout) const; - + /// Formula with atoms in alphabetic order std::string GetFormula() const; /// Weight for the crystal formula, in atomic units (g/mol). This should be diff --git a/src/ObjCryst/ObjCryst/Molecule.cpp b/src/ObjCryst/ObjCryst/Molecule.cpp index b6ffbf2f..09909ed0 100644 --- a/src/ObjCryst/ObjCryst/Molecule.cpp +++ b/src/ObjCryst/ObjCryst/Molecule.cpp @@ -7962,7 +7962,7 @@ Molecule *ZScatterer2Molecule(ZScatterer *scatt) z0+=z; mol->AddAtom(x,y,z,scatt->GetZAtomRegistry().GetObj(i).GetScatteringPower(), scatt->GetZAtomRegistry().GetObj(i).GetName()); - + #if 0 if(i>0) { @@ -8014,7 +8014,7 @@ Molecule *ZScatterer2Molecule(ZScatterer *scatt) #endif mol->GetAtom(i).SetOccupancy(scatt->GetZAtomRegistry().GetObj(i).GetOccupancy()); } - + CrystVector_REAL x(nb),y(nb),z(nb),radius(nb); vector > scattpow(nb); for(unsigned int i=0;iRemovePar(&this->GetPar(mScaleFactor.data()+mPowderPatternComponentRegistry.GetNb()-1)); this->Print(); } - + this->RemoveSubRefObj(comp); comp.DeRegisterClient(*this); mClockPowderPatternCalc.Reset(); mClockIntegratedFactorsPrep.Reset(); mPowderPatternComponentRegistry.DeRegister(comp); - + // Shift scale factors unsigned int i=0; for(unsigned int i=0;iGetNbPowderPatternComponent();i++) @@ -6968,7 +6968,7 @@ SPGScore SpaceGroupExplorer::Run(const cctbx::sgtbx::space_group &spg, const boo lsq.PrepareRefParList(true); const unsigned int saved_par = lsq.GetCompiledRefinedObj().CreateParamSet("SpaceGroupExplorer saved parameters"); lsq.GetCompiledRefinedObj().SaveParamSet(saved_par); - + Chronometer chrono; for(unsigned int j=0;jGetNbReflBelowMaxSinThetaOvLambda(); if(verbose>0) cout << boost::format(" Rwp= %5.2f%% GoF=%9.2f nGoF =%9.2f (%3u reflections, %3u extinct)") % rw % gof % ngof % nbrefl % nbextinct446<GetCrystal().Init(a,b,c,d,e,f,spghm,name); @@ -7059,7 +7059,7 @@ void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, c const REAL relative_length_tolerance, const REAL absolute_angle_tolerance_degree) { Crystal *pCrystal=&(mpDiff->GetCrystal()); - + // Initial lattice parameters & spg const REAL a=pCrystal->GetLatticePar(0), b=pCrystal->GetLatticePar(1), @@ -7086,7 +7086,7 @@ void SpaceGroupExplorer::RunAll(const bool fitprofile_all, const bool verbose, c mvSPG.clear(); mvSPGExtinctionFingerprint.clear(); - + // Nb refl below max sin(theta/lambda) for p1, to compute nGoF unsigned int nb_refl_p1=1; @@ -7170,7 +7170,7 @@ REAL SpaceGroupExplorer::GetP1IntegratedGoF() // <(mP1IntegratedProfileMin, mP1IntegratedProfileMax,mP1IntegratedProfileMax)<(mpDiff->GetH(), mpDiff->GetK(), mpDiff->GetL(), mpDiff->GetFhklCalcSq()); REAL integratedChi2=0.; REAL integratedChi2LikeNorm=0.; @@ -7199,5 +7199,5 @@ REAL SpaceGroupExplorer::GetP1IntegratedGoF() } return chi2 / nbpoint; } - + }//namespace ObjCryst diff --git a/src/ObjCryst/ObjCryst/PowderPattern.h b/src/ObjCryst/ObjCryst/PowderPattern.h index 0a907c0a..86bf3a4b 100644 --- a/src/ObjCryst/ObjCryst/PowderPattern.h +++ b/src/ObjCryst/ObjCryst/PowderPattern.h @@ -1097,7 +1097,7 @@ class PowderPattern : public RefinableObj /// This is mutable because generally we use the 'best' scale factor, but /// it should not be... mutable CrystVector_REAL mScaleFactor; - + /// Mu*R parameter for cylinder absorption correction REAL mMuR; @@ -1189,7 +1189,7 @@ struct SPGScore string hm; /// Rw factor, from PowderPattern::GetRw() REAL rw; - /// Goodness-of-fit, from PowderPattern::GetChi2() normalised by + /// Goodness-of-fit, from PowderPattern::GetChi2() normalised by REAL gof; /// Normalised & integrated goodness-of-fit = iGoF * (nb_refl) / (nb_refl in P1), /// where iGoF is the integrated goodness-of-fit, computed using the integration @@ -1289,7 +1289,7 @@ class SpaceGroupExplorer /// Map extinction fingerprint std::map,SPGScore> mvSPGExtinctionFingerprint; }; - + }//namespace ObjCryst #endif // _OBJCRYST_POWDERPATTERN_H_ diff --git a/src/ObjCryst/RefinableObj/GlobalOptimObj.cpp b/src/ObjCryst/RefinableObj/GlobalOptimObj.cpp index bf5dd2be..84079fac 100644 --- a/src/ObjCryst/RefinableObj/GlobalOptimObj.cpp +++ b/src/ObjCryst/RefinableObj/GlobalOptimObj.cpp @@ -2299,16 +2299,16 @@ void MonteCarloObj::InitLSQ(const bool useFullPowderPatternProfile) } // Only refine structural parameters (excepting parameters already fixed) and scale factor mLSQ.PrepareRefParList(true); - + // Intensity corrections can be refined std::list vIntCorrPar; for(int i=0; iIsDescendantFromOrSameAs(gpRefParTypeScattDataCorrInt) && mLSQ.GetCompiledRefinedObj().GetPar(i).IsFixed()==false) vIntCorrPar.push_back(&mLSQ.GetCompiledRefinedObj().GetPar(i)); - + mLSQ.SetParIsFixed(gpRefParTypeScattData,true); mLSQ.SetParIsFixed(gpRefParTypeScattDataScale,false); - + for(std::list::iterator pos=vIntCorrPar.begin();pos!=vIntCorrPar.end();pos++) (*pos)->SetIsFixed(false); mLSQ.SetParIsFixed(gpRefParTypeUnitCell,true); diff --git a/src/ObjCryst/RefinableObj/RefinableObj.cpp b/src/ObjCryst/RefinableObj/RefinableObj.cpp index 52ce27c3..e0520789 100644 --- a/src/ObjCryst/RefinableObj/RefinableObj.cpp +++ b/src/ObjCryst/RefinableObj/RefinableObj.cpp @@ -931,10 +931,10 @@ template void ObjRegistry::DeRegister(T &obj) if(0!=mpWXRegistry) mpWXRegistry->Remove(obj.WXGet()); #endif mvpRegistry.erase(pos); - + typename list::iterator pos2=find(mvpRegistryList.begin(),mvpRegistryList.end(),&obj); mvpRegistryList.erase(pos2); - + mListClock.Click(); VFN_DEBUG_EXIT("ObjRegistry("<Protect(); return gm; } diff --git a/src/newmat/newmat4.cpp b/src/newmat/newmat4.cpp index 5c5cc41d..34f6acae 100644 --- a/src/newmat/newmat4.cpp +++ b/src/newmat/newmat4.cpp @@ -200,12 +200,12 @@ CroutMatrix::CroutMatrix(const BaseMatrix& m) void CroutMatrix::get_aux(CroutMatrix& X) { X.d = d; X.sing = sing; - if (tag_val == 0 || tag_val == 1) // reuse the array + if (tag_val == 0 || tag_val == 1) // reuse the array { REPORT X.indx = indx; indx = 0; d = true; sing = true; return; } else if (nrows_val == 0) { REPORT indx = 0; d = true; sing = true; return; } else // copy the array - { + { REPORT Tracer tr("CroutMatrix::get_aux"); int *ix = new int [nrows_val]; MatrixErrorNoSpace(ix); @@ -259,7 +259,7 @@ ReturnMatrix GeneralMatrix::for_return() const Matrix::Matrix(Real a, int m, int n) : GeneralMatrix(m * n) { REPORT nrows_val=m; ncols_val=n; operator=(a); } - + Matrix::Matrix(const Real* a, int m, int n) : GeneralMatrix(m * n) { REPORT nrows_val=m; ncols_val=n; *this << a; } @@ -433,7 +433,7 @@ void Matrix::resize_keep(int nr, int nc) { Tracer tr("Matrix::resize_keep"); if (nr == nrows_val && nc == ncols_val) { REPORT return; } - + if (nr <= nrows_val && nc <= ncols_val) { REPORT @@ -456,7 +456,7 @@ void Matrix::resize_keep(int nr, int nc) X.submatrix(1,nr,1,nc) = submatrix(1,nr,1,nc); swap(X); } -} +} void SquareMatrix::resize_keep(int nr) { @@ -483,7 +483,7 @@ void SquareMatrix::resize_keep(int nr, int nc) if (nr != nc) Throw(NotSquareException(*this)); resize_keep(nr); } - + void SymmetricMatrix::resize_keep(int nr) { @@ -501,7 +501,7 @@ void SymmetricMatrix::resize_keep(int nr) X.sym_submatrix(1,nrows_val) = *this; swap(X); } -} +} void UpperTriangularMatrix::resize_keep(int nr) { @@ -519,7 +519,7 @@ void UpperTriangularMatrix::resize_keep(int nr) X.sym_submatrix(1,nrows_val) = *this; swap(X); } -} +} void LowerTriangularMatrix::resize_keep(int nr) { @@ -537,7 +537,7 @@ void LowerTriangularMatrix::resize_keep(int nr) X.sym_submatrix(1,nrows_val) = *this; swap(X); } -} +} void DiagonalMatrix::resize_keep(int nr) { @@ -555,7 +555,7 @@ void DiagonalMatrix::resize_keep(int nr) X.sym_submatrix(1,nrows_val) = *this; swap(X); } -} +} void RowVector::resize_keep(int nc) { @@ -599,7 +599,7 @@ void ColumnVector::resize_keep(int nr) X.rows(1,nrows_val) = *this; swap(X); } -} +} void ColumnVector::resize_keep(int nr, int nc) { @@ -683,7 +683,7 @@ MatrixBandWidth BandMatrix::bandwidth() const MatrixBandWidth BandLUMatrix::bandwidth() const { REPORT return MatrixBandWidth(m1,m2); } - + MatrixBandWidth GenericMatrix::bandwidth()const { REPORT return gm->bandwidth(); } @@ -1271,7 +1271,7 @@ void GeneralMatrix::swap(GeneralMatrix& gm) t = storage; storage = gm.storage; gm.storage = t; Real* s = store; store = gm.store; gm.store = s; } - + void nricMatrix::swap(nricMatrix& gm) { REPORT @@ -1339,7 +1339,7 @@ RealStarStar::RealStarStar(Matrix& A) MatrixErrorNoSpace(a); Real* d = A.data(); for (int i = 0; i < m; ++i) a[i] = d + i * n; -} +} ConstRealStarStar::ConstRealStarStar(const Matrix& A) { @@ -1351,7 +1351,7 @@ ConstRealStarStar::ConstRealStarStar(const Matrix& A) MatrixErrorNoSpace(a); const Real* d = A.data(); for (int i = 0; i < m; ++i) a[i] = d + i * n; -} +} diff --git a/src/newmat/newmat5.cpp b/src/newmat/newmat5.cpp index 703fb00e..83b495b9 100644 --- a/src/newmat/newmat5.cpp +++ b/src/newmat/newmat5.cpp @@ -303,7 +303,7 @@ GeneralMatrix* GetSubMatrix::Evaluate(MatrixType mt) else Compare(gm->Type().sub(), mt); GeneralMatrix* gmx = mt.New(row_number, col_number, this); int i = row_number; - MatrixRow mr(gm, LoadOnEntry, row_skip); + MatrixRow mr(gm, LoadOnEntry, row_skip); MatrixRow mrx(gmx, StoreOnExit+DirectPart); MatrixRowCol sub; while (i--) @@ -330,7 +330,7 @@ void GeneralMatrix::Add(GeneralMatrix* gm1, Real f) { *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; *s++ = *s1++ + f; } i = storage & 3; while (i--) *s++ = *s1++ + f; } - + void GeneralMatrix::Add(Real f) { REPORT @@ -338,7 +338,7 @@ void GeneralMatrix::Add(Real f) while (i--) { *s++ += f; *s++ += f; *s++ += f; *s++ += f; } i = storage & 3; while (i--) *s++ += f; } - + void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f) { REPORT @@ -347,7 +347,7 @@ void GeneralMatrix::NegAdd(GeneralMatrix* gm1, Real f) { *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; *s++ = f - *s1++; } i = storage & 3; while (i--) *s++ = f - *s1++; } - + void GeneralMatrix::NegAdd(Real f) { REPORT @@ -359,7 +359,7 @@ void GeneralMatrix::NegAdd(Real f) } i = storage & 3; while (i--) { *s = f - *s; s++; } } - + void GeneralMatrix::Negate(GeneralMatrix* gm1) { // change sign of elements @@ -369,7 +369,7 @@ void GeneralMatrix::Negate(GeneralMatrix* gm1) { *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); *s++ = -(*s1++); } i = storage & 3; while(i--) *s++ = -(*s1++); } - + void GeneralMatrix::Negate() { REPORT @@ -378,7 +378,7 @@ void GeneralMatrix::Negate() { *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; *s = -(*s); s++; } i = storage & 3; while(i--) { *s = -(*s); s++; } } - + void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f) { REPORT @@ -387,7 +387,7 @@ void GeneralMatrix::Multiply(GeneralMatrix* gm1, Real f) { *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; *s++ = *s1++ * f; } i = storage & 3; while (i--) *s++ = *s1++ * f; } - + void GeneralMatrix::Multiply(Real f) { REPORT @@ -395,7 +395,7 @@ void GeneralMatrix::Multiply(Real f) while (i--) { *s++ *= f; *s++ *= f; *s++ *= f; *s++ *= f; } i = storage & 3; while (i--) *s++ *= f; } - + /************************ MatrixInput routines ****************************/ diff --git a/src/newmat/newmat6.cpp b/src/newmat/newmat6.cpp index eda39091..575ff08f 100644 --- a/src/newmat/newmat6.cpp +++ b/src/newmat/newmat6.cpp @@ -351,7 +351,7 @@ void Matrix::operator=(const BaseMatrix& X) REPORT //CheckConversion(X); // MatrixConversionCheck mcc; Eq(X,MatrixType::Rt); -} +} void SquareMatrix::operator=(const BaseMatrix& X) { @@ -432,7 +432,7 @@ void CroutMatrix::operator=(const CroutMatrix& gm) ((CroutMatrix&)gm).get_aux(*this); Eq(gm); } - + diff --git a/src/newmat/newmat7.cpp b/src/newmat/newmat7.cpp index 5befa5d9..17bd7783 100644 --- a/src/newmat/newmat7.cpp +++ b/src/newmat/newmat7.cpp @@ -1033,7 +1033,7 @@ ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B) Real* bn = b+n; Real* bn2 = bn+n; Real* cn = c+n; Real* cn2 = cn+n; - int i = n; + int i = n; while (i--) { *c++ = *an * *bn2 - *an2 * *bn; @@ -1050,4 +1050,3 @@ ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B) #endif ///@} - diff --git a/src/newmat/newmat8.cpp b/src/newmat/newmat8.cpp index f48edc6b..d3796a03 100644 --- a/src/newmat/newmat8.cpp +++ b/src/newmat/newmat8.cpp @@ -750,7 +750,7 @@ ReturnMatrix BaseMatrix::sum_square_rows() const int s = mr.Storage(); Real* in = mr.Data(); while (s--) sum += square(*in++); - ssq(i) = sum; + ssq(i) = sum; mr.Next(); } } @@ -795,7 +795,7 @@ ReturnMatrix BaseMatrix::sum_rows() const int s = mr.Storage(); Real* in = mr.Data(); while (s--) sum += *in++; - sum_vec(i) = sum; + sum_vec(i) = sum; mr.Next(); } } diff --git a/src/newmat/newmatap.h b/src/newmat/newmatap.h index dc88ac14..e6884d8d 100644 --- a/src/newmat/newmatap.h +++ b/src/newmat/newmatap.h @@ -25,13 +25,13 @@ void QRZT(const Matrix&, Matrix&, Matrix&); void QRZ(Matrix&, UpperTriangularMatrix&); -void QRZ(const Matrix&, Matrix&, Matrix&); +void QRZ(const Matrix&, Matrix&, Matrix&); inline void QRZT(Matrix& X, Matrix& Y, LowerTriangularMatrix& L, Matrix& M) - { QRZT(X, L); QRZT(X, Y, M); } + { QRZT(X, L); QRZT(X, Y, M); } inline void QRZ(Matrix& X, Matrix& Y, UpperTriangularMatrix& U, Matrix& M) - { QRZ(X, U); QRZ(X, Y, M); } + { QRZ(X, U); QRZ(X, Y, M); } inline void HHDecompose(Matrix& X, LowerTriangularMatrix& L) { QRZT(X,L); } @@ -60,7 +60,7 @@ inline void UpdateQRZ(UpperTriangularMatrix& X, UpperTriangularMatrix& U) inline void UpdateQRZ(const UpperTriangularMatrix& X, Matrix& MX, Matrix& MU) { updateQRZ(X, MX, MU); } - + inline void UpdateQRZ(const Matrix& X, Matrix& MX, Matrix& MU) { updateQRZ(X, MX, MU); } @@ -68,7 +68,7 @@ inline void UpdateQRZ(const Matrix& X, Matrix& MX, Matrix& MU) // Matrix A's first n columns are orthonormal // so A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix. // Fill out the remaining columns of A to make them orthonormal -// so A.t() * A is the identity matrix +// so A.t() * A is the identity matrix void extend_orthonormal(Matrix& A, int n); @@ -99,9 +99,9 @@ inline void RightCircularUpdateCholesky(UpperTriangularMatrix &chol, // a LEFT circular shift of the rows and columns from // 1,...,k-1,k,k+1,...l,l+1,...,p to // 1,...,k-1,k+1,...l,k,l+1,...,p to -void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l); +void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l); inline void LeftCircularUpdateCholesky(UpperTriangularMatrix &chol, - int k, int l) { left_circular_update_Cholesky(chol, k, l); } + int k, int l) { left_circular_update_Cholesky(chol, k, l); } void SVD(const Matrix&, DiagonalMatrix&, Matrix&, Matrix&, @@ -268,5 +268,3 @@ ReturnMatrix Helmert_transpose(const Matrix& Y, bool full=false); ///@} - - diff --git a/src/newmat/newmatex.cpp b/src/newmat/newmatex.cpp index d00496ec..0df2657d 100644 --- a/src/newmat/newmatex.cpp +++ b/src/newmat/newmatex.cpp @@ -322,4 +322,3 @@ void BaseMatrix::IEQND() const #ifdef use_namespace } #endif - diff --git a/src/newmat/newmatnl.cpp b/src/newmat/newmatnl.cpp index 56c52ea2..35c4b46d 100644 --- a/src/newmat/newmatnl.cpp +++ b/src/newmat/newmatnl.cpp @@ -237,7 +237,7 @@ void MLE_D_FI::Fit(ColumnVector& Parameters) FindMaximum2::Fit(Parameters,Lim); cout << "\nConverged" << endl; } - + void MLE_D_FI::MakeCovariance() { if (Covariance.Nrows()==0) @@ -252,7 +252,7 @@ void MLE_D_FI::MakeCovariance() void MLE_D_FI::GetStandardErrors(ColumnVector& SEX) { MakeCovariance(); SEX = SE.AsColumn(); } - + void MLE_D_FI::GetCorrelations(SymmetricMatrix& Corr) { MakeCovariance(); Corr << SE.i() * Covariance * SE.i(); } diff --git a/src/newmat/newmatnl.h b/src/newmat/newmatnl.h index 0a415845..7477d010 100644 --- a/src/newmat/newmatnl.h +++ b/src/newmat/newmatnl.h @@ -249,7 +249,7 @@ class NonLinearLeastSquares : public FindMaximum2 // The next class is the prototype class for calculating the // log-likelihood. -// I assume first derivatives are available and something like the +// I assume first derivatives are available and something like the // Fisher Information or variance/covariance matrix of the first // derivatives or minus the matrix of second derivatives is // available. This matrix must be positive definite. @@ -324,4 +324,3 @@ class MLE_D_FI : public FindMaximum2 ///@} - diff --git a/src/newmat/newmatrm.cpp b/src/newmat/newmatrm.cpp index 726d10f1..66e5e27f 100644 --- a/src/newmat/newmatrm.cpp +++ b/src/newmat/newmatrm.cpp @@ -141,7 +141,7 @@ void ComplexScale(RectMatrixCol& U, RectMatrixCol& V, Real x, Real y) Tracer tr("newmatrm"); Throw(InternalException("Dimensions differ in ComplexScale")); } - Real* u = U.store; Real* v = V.store; + Real* u = U.store; Real* v = V.store; int su = U.spacing; int sv = V.spacing; //while (n--) //{ diff --git a/src/newmat/newmatrm.h b/src/newmat/newmatrm.h index d8494d73..ad64e18f 100644 --- a/src/newmat/newmatrm.h +++ b/src/newmat/newmatrm.h @@ -130,7 +130,7 @@ inline void GivensRotation(Real cGivens, Real sGivens, Real& x, Real& y) Real tmp1 = -sGivens * x + cGivens * y; x = tmp0; y = tmp1; } - + inline void GivensRotationR(Real cGivens, Real sGivens, Real& x, Real& y) { // also change sign of y @@ -138,7 +138,7 @@ inline void GivensRotationR(Real cGivens, Real sGivens, Real& x, Real& y) Real tmp0 = cGivens * x + sGivens * y; Real tmp1 = sGivens * x - cGivens * y; x = tmp0; y = tmp1; -} +} diff --git a/src/newmat/nl_ex.cpp b/src/newmat/nl_ex.cpp index 25c9a012..66cf73b7 100644 --- a/src/newmat/nl_ex.cpp +++ b/src/newmat/nl_ex.cpp @@ -95,7 +95,7 @@ int my_main() #ifdef DO_FREE_CHECK FreeCheck::Status(); #endif - + return 0; } @@ -113,13 +113,10 @@ int main() } CatchAll { - cout << "\nProgram fails - exception generated\n\n"; + cout << "\nProgram fails - exception generated\n\n"; } return 0; } ///@} - - - diff --git a/src/newmat/nl_ex.txt b/src/newmat/nl_ex.txt index 46dcbf4b..4982562c 100644 --- a/src/newmat/nl_ex.txt +++ b/src/newmat/nl_ex.txt @@ -39,4 +39,3 @@ Var/cov 0.5295 -0.0324 -0.0250 -0.0324 0.6668 -0.0096 -0.0250 -0.0096 0.0014 - diff --git a/src/newmat/nm11.htm b/src/newmat/nm11.htm index f34829ad..d920b4d6 100644 --- a/src/newmat/nm11.htm +++ b/src/newmat/nm11.htm @@ -1,7 +1,7 @@ - + - + @@ -27,7 +27,7 @@

Documentation for newmat11, a matrix library in C++

3. Reference manual
4. Error messages - + 5. Design of the library
6. Function summary
7. Change History
@@ -37,7 +37,7 @@

Documentation for newmat11, a matrix library in C++

 

This is the how to use documentation for newmat plus some background information on its design.

-

There is additional support material on my web site. +

There is additional support material on my web site.

Navigation:  This page is arranged in sections, sub-sections and sub-sub-sections; four cross-references are given at the top @@ -72,7 +72,7 @@

1.1 Conditions of use

I place no restrictions on the use of newmat except that I take no liability for any problems that may arise from its use, distribution or other dealings with it.

-You can use it in your commercial projects (as well as your non-commercial +You can use it in your commercial projects (as well as your non-commercial projects).

You can make and distribute modified or merged versions. You can @@ -91,7 +91,7 @@

1.1 Conditions of use

your own risk. I (Robert Davies) take no responsibility for any errors or omissions in this package or for any misfortune that may befall you or others as a result of your use, distribution or other dealings with it.
-

Please report bugs to me at robert at statsresearch.co.nz   +

Please report bugs to me at robert at statsresearch.co.nz   [replace at by you-know-what-character in the email address].

When reporting a bug please tell me which C++ compiler you are using, and what version. Also give me details of your computer. And tell me which version @@ -197,7 +197,7 @@

1.3 Is this the library for operators so you can write things like X = A * (B + C);
  • need a variety of types of matrices (but not sparse)
  • need only one element type (float or double)
    • -
  • work with matrices in the range 10 x 10 up to what can be stored in memory +
  • work with matrices in the range 10 x 10 up to what can be stored in memory
  • tolerate a moderately large but not huge package
  • need high quality but not necessarily the latest numerical methods.
    • @@ -295,14 +295,14 @@

    2.1 Overview

    used for testing, example.cpp is an example and sl_ex.cpp, nl_ex.cpp and garch.cpp are examples of the non-linear solve and optimisation routines. A -demonstration and test of the exception mechanism is in test_exc.cpp. The files +demonstration and test of the exception mechanism is in test_exc.cpp. The files nm_ex1.cpp, nm_ex2.cpp and nm_ex3.cpp contain more simple examples.

    I include a number of make files for compiling the example and the test package. See the section on make files for details. Alternatively, with the PC compilers, its pretty quick just to load all the files in the interactive environments by pointing and clicking.

    Use the large or win32 console model when you are using a PC. Do not -outline inline functions. You may need to increase the stack size on +outline inline functions. You may need to increase the stack size on older operating systems or compilers.

    Your source files that access the newmat will need to #include one or more of the following files.

    @@ -338,18 +338,18 @@

    2.1 Overview

    2.2 Make files

    next - skip - up - start

    -

    I have included make files for compiling my test and example programs for -various versions of CC, Borland, Microsoft, Open Watcom, Intel and Gnu compilers. You +

    I have included make files for compiling my test and example programs for +various versions of CC, Borland, Microsoft, Open Watcom, Intel and Gnu compilers. You can generate make files for a number of other compilers with my genmake utility. Make files provide a way of compiling your programs without using the IDE that comes with PC compilers. See the files section for details.

    PC

    -

    I include make files for various versions of Borland, Microsoft and Intel -compilers. With the -Borland compiler you will need to edit it to show where +

    I include make files for various versions of Borland, Microsoft and Intel +compilers. With the +Borland compiler you will need to edit it to show where you have stored your Borland compiler. For make files for other compilers use my -genmake utility. To compile my test and example +genmake utility. To compile my test and example programs use Borland 5.5 (Builder 5) use

       make -f nm_b55.mak

    or with Borland 5.6 (Builder 6) use

    @@ -400,7 +400,7 @@

    2.3 Customising

    compiler supported exceptions if and only if you have set the option on your compiler to recognise exceptions. Disabling exceptions sometimes helps with compilers that are incompatible with my exception simulation scheme.

    -

    Activate the appropriate statement to make the element type float or +

    Activate the appropriate statement to make the element type float or double. I suggest you leave it at double.

    The option DO_FREE_CHECK is used for tracking memory leaks and normally should not be activated.

    @@ -408,7 +408,7 @@

    2.3 Customising

    element access. Note that this does not change the starting point for indices when you are using round brackets for accessing elements or selecting submatrices. It does enable you to use C style square -brackets. This also activates additional constructors for Matrix, ColumnVector +brackets. This also activates additional constructors for Matrix, ColumnVector and RowVector to simplify use with Numerical Recipes in C++.

    Activate #define use_namespace if you want to use namespaces. Do this only if you are sure your compiler @@ -417,7 +417,7 @@

    2.3 Customising

    missing links.

    Activate #define _STANDARD_ to use the standard names for the included files and to find the floating point precision data using the floating -point standard. This will work with most recent compilers and is done +point standard. This will work with most recent compilers and is done automatically for Borland, Gnu, Intel and Microsoft compilers.

    If you haven't defined _STANDARD_ and are using a compiler that include.h does not recognise and you want to pick up the floating point @@ -427,7 +427,7 @@

    2.3 Customising

    particular.

    There is a line

       //#define set_unix_options
    -

    You can activate this if you are using a Linux or Unix system. It is not used +

    You can activate this if you are using a Linux or Unix system. It is not used by Newmat but is used by some of my other programs.

    2.4 Compilers

    next - skip - @@ -449,7 +449,7 @@

    2.4 Compilers

    I have tested this library on a number of compilers. Here are the levels of success and any special considerations. In most cases I have chosen code that works under all the compilers I have access to, but I have had to include some -specific work-arounds for some compilers. For the newest PC versions, I use a Pentium +specific work-arounds for some compilers. For the newest PC versions, I use a Pentium 4 computer running windows XP or Linux (Red Hat workstation version). The Unix versions are on a Sun Sparc station. Thanks to Victoria University for access to the Sparc.

    @@ -462,7 +462,7 @@

    2.4.1 AT&T

    next - skip - up - start

    The AT&T compiler used to be available on a wide variety of Unix -workstations. I don't know if anyone still uses it. However the AT&T options are +workstations. I don't know if anyone still uses it. However the AT&T options are the default if your compiler is not recognised.

    AT&T C++ 2.1; 3.0.1 on a Sun: Previous versions worked on these compilers, which I no longer have access to.

    @@ -473,25 +473,25 @@

    2.4.2 Borland

    next - skip - up - start

    Newer compilers

    -

    Borland Builder version 8. My tests here have been compiling from a -make file using nm_b58.mak and making a console program. This works fine -except that my values of LnMinimum and LnMaximum in precisio.h +

    Borland Builder version 8. My tests here have been compiling from a +make file using nm_b58.mak and making a console program. This works fine +except that my values of LnMinimum and LnMaximum in precisio.h are not correct.

    -

    Borland Builder version 6: My tests have been on the personal -version. See the notes for version 5. If you are -compiling with a make file you can use nm_b56.mak as a model. You can set -the newmat options to use namespace and the standard library. If you are +

    Borland Builder version 6: My tests have been on the personal +version. See the notes for version 5. If you are +compiling with a make file you can use nm_b56.mak as a model. You can set +the newmat options to use namespace and the standard library. If you are compiling a GUI program you may need to comment out the line defining -TypeDefException in include.h. I don't believe exceptions work -completely correctly in either version 5 or version 6. However, this does not +TypeDefException in include.h. I don't believe exceptions work +completely correctly in either version 5 or version 6. However, this does not seem to be a problem with my use of them in newmat.

    Borland Builder version 5: This works fine in console mode and no special editing of the source codes is required. I haven't tested it in GUI mode. You can set the newmat options to use namespace and the standard library. You -should turn off the Borland option to use pre-compiled headers. There -are notes on compiling with the IDE on my website. +should turn off the Borland option to use pre-compiled headers. There +are notes on compiling with the IDE on my website. Alternatively you can use the nm_b55.mak make file.

    -

    Borland Builder version 4: I have successfully used this on older +

    Borland Builder version 4: I have successfully used this on older versions of newmat using the console wizard (menu item file/new - select new tab). Use compiler exceptions. Suppose you are compiling my test program tmt. Rename my @@ -508,10 +508,10 @@

    Newer compilers

    Older compilers

    Borland C++ 5.02:

    I am not longer checking compatibility with this compiler.

    -

    Use the large or 32 bit flat model. If you are not debugging, turn off the +

    Use the large or 32 bit flat model. If you are not debugging, turn off the options that collect debugging information. Use my simulated exceptions.

    When running my test program under ms-dos you may run out of memory. Either -compile the test routine to run under easywin or use simulated exceptions +compile the test routine to run under easywin or use simulated exceptions rather than the built in exceptions.

    If you can, upgrade to windows 95 or window NT and use the 32 bit console model.

    @@ -519,23 +519,23 @@

    Older compilers

    less than 64K bytes in length (90x90 for a rectangular matrix if you are using double as your element type). Otherwise your program will crash without warning or explanation. You will need to break the tmt set of test files into several parts to get the program +HREF="#testing">tmt set of test files into several parts to get the program to fit into your computer and run without stack overflow.

    You can generate make files for versions 5 with my genmake utility.

    Borland C++ 3 and 4.

    -

    The program used to compile in version 3.1 if you enable the simulated booleans -- comment out the line #define bool_LIB 0 in include.h -and use the simulated exceptions. I haven't checked the latest versions -of Newmat. The main test program is too large to run +

    The program used to compile in version 3.1 if you enable the simulated booleans +- comment out the line #define bool_LIB 0 in include.h +and use the simulated exceptions. I haven't checked the latest versions +of Newmat. The main test program is too large to run unless you break it up into several parts. I haven't tried it under version 4.

    2.4.3 Gnu G++

    next - skip - up - start

    Gnu G++ 3, 4 (Linux), 3 (Sun): These work OK. If you are using a much earlier version -see if you can upgrade. It  used to work with 2.95 and 2.96 but I don't -have access to these now. You can't use standard with the 2.9X -versions. The namespace option worked with 2.96 on Linux but not with 2.95 on +see if you can upgrade. It  used to work with 2.95 and 2.96 but I don't +have access to these now. You can't use standard with the 2.9X +versions. The namespace option worked with 2.96 on Linux but not with 2.95 on the Sun. Standard is automatically turned on with the 3.X.

    This version of Newmat is not compatible with versions 2.6 or earlier.

    2.4.4 HP-UX

    @@ -552,11 +552,11 @@

    2.4.5 Intel

    next - skip - up - start

    -

    Newmat works correctly with the Intel 10 C++ compiler for Windows and Linux. (Not tested for the other versions). Standard is -automatically turned on for both the Linux versions and Windows versions. If +

    Newmat works correctly with the Intel 10 C++ compiler for Windows and Linux. (Not tested for the other versions). Standard is +automatically turned on for both the Linux versions and Windows versions. If this causes a problem for the version you are using you can find the lines in -include.h that control this and comment them out. Note that the Intel compiler for -Linux is free for non-commercial use. (One of the versions of 8.1 gave a +include.h that control this and comment them out. Note that the Intel compiler for +Linux is free for non-commercial use. (One of the versions of 8.1 gave a warning message every time I had something like Real x; ... if (x==0.0) ..., which was often. This is now seems to be fixed.)

    @@ -566,15 +566,15 @@

    2.4.6 Microsoft

    Newer versions

    See my web site for instructions how to work Microsoft's IDE.

    -

    Microsoft Visual C++ 7, 7.1, 8: These work OK. All my tests have -been in console mode. You can turn on my namespace option. Standard is turned on +

    Microsoft Visual C++ 7, 7.1, 8: These work OK. All my tests have +been in console mode. You can turn on my namespace option. Standard is turned on by default for these versions.

    -

    Microsoft Visual C++ 6: Get the latest service pack. I have tried this +

    Microsoft Visual C++ 6: Get the latest service pack. I have tried this in console mode and it seems to work satisfactorily. Use the compiler supported exceptions. You may be able to use the namespace and standard options. If you want to work under MFC -you may need to #include "stdafx.h" at the beginning of each .cpp file +you may need to #include "stdafx.h" at the beginning of each .cpp file (or turn off precompiled headers).

    -

    Microsoft Visual C++ 5: I have tried this in console mode on previous +

    Microsoft Visual C++ 5: I have tried this in console mode on previous versions of Newmat. It seems to work satisfactorily. There may be a problem with namespace (fixed by Service Pack 3?). Turn optimisation @@ -587,7 +587,7 @@

    Older versions

    2.4.7 Sun

    next - skip - up - start

    -

    Sun C++ (version = ?): This seems to work fine with +

    Sun C++ (version = ?): This seems to work fine with compiler supported exceptions. Sun C++ (version 5): There was a problem with exceptions. If you use my simulated exceptions the non-linear optimisation programs hang. If you use the compiler @@ -596,48 +596,48 @@

    2.4.7 Sun

    2.4.8 Watcom

    next - skip - up - start

    -

    Open Watcom (version 1.7a): this works. You can set the standard -option in include.h. The scientific and fixed +

    Open Watcom (version 1.7a): this works. You can set the standard +option in include.h. The scientific and fixed manipulators don't work.

    -

    Watcom C++ (version 10a): this used to work, I don't know if it works +

    Watcom C++ (version 10a): this used to work, I don't know if it works now.

    2.5 Updating from previous versions

    next - skip - up - start

    -

    Newmat11 - if you are upgrading from earlier versions note the +

    Newmat11 - if you are upgrading from earlier versions note the following:

    • The class Exception in myexcept.h has been replaced by - BaseException and a typedef statement included so programs that use - the Exception class will still work. If the Exception class was + BaseException and a typedef statement included so programs that use + the Exception class will still work. If the Exception class was causing a problem comment out the line defining TypeDefException in include.h.
      • -
    • Newmat11 does not support the TEMP_DESTROYED_QUICKLY options so won't work +
    • Newmat11 does not support the TEMP_DESTROYED_QUICKLY options so won't work with very old versions of Gnu G++
    • Old AT&T work-arounds are removed
      • -
    • The simulated Booleans class is now stored at the end of include.h. In +
    • The simulated Booleans class is now stored at the end of include.h. In general, I am not testing with compilers that don't support bool.
      • -
    • QRZ and QRZT no longer throw an exception if they generate a singular +
    • QRZ and QRZT no longer throw an exception if they generate a singular triangular matrix
    • Converting functions to lower case. Both the old and new names will work. See the list of functions for new and old names
    • nm_i5.mak, nm_il5.mak replaced by nm_i8.mak, nm_il8.mak
    • extra file nm_misc.cpp to be included in compilation
      • -
    • scientific and fixed manipulators now recognised in output +
    • scientific and fixed manipulators now recognised in output statements. If you want fixed output and have been using scientific in a previous output statement, you may need to include the - fixed manipulator. (Doesn't work in Visual C++, version 6 and open + fixed manipulator. (Doesn't work in Visual C++, version 6 and open Watcom).

    Newmat10 includes new maxima, minima, determinant, dot product and Frobenius norm functions, a faster FFT, revised make files for GCC -and CC compilers, several corrections, new ReSize function, IdentityMatrix +and CC compilers, several corrections, new ReSize function, IdentityMatrix and Kronecker Product. Singular values from SVD are sorted. The program files include a new file, newfft.cpp, so you -will need to include this in the list of files in your IDE and make files. There +will need to include this in the list of files in your IDE and make files. There is also a new test file tmtm.cpp. Pointer arithmetic now mostly meets requirements of standard. You can use << to load data into rows @@ -650,7 +650,7 @@

    2.5 Updating from previous
  • Boolean, TRUE, FALSE are now bool, true, false. See customising if your compiler supports the bool class.
    • -
  • ReDimension is now ReSize. +
  • ReDimension is now ReSize.
  • The simulated exception package has been updated.
    • @@ -692,7 +692,7 @@

    2.5 Updating from previous Edit the file include.h to select one of these three options. Don't simulate exceptions if you have set your compiler's option to implement exceptions. -
  • New QR decomposition functions. +
  • New QR decomposition functions.
  • A non-linear least squares class.
    • @@ -700,7 +700,7 @@

    2.5 Updating from previous
  • Concatenation and elementwise multiplication.
    • -
  • A new GenericMatrix class. +
  • A new GenericMatrix class.
  • Sum function.
  • Some of the make files reorganised. @@ -754,8 +754,8 @@

    2.6 Catching exceptions

    return 0; }

    Or see my file nm_ex1.cpp for an easy way of organising this.

    -

    If you are using a GUI version rather a console version of the program you -will need to catch the exception and display the error message in a pop-up +

    If you are using a GUI version rather a console version of the program you +will need to catch the exception and display the error message in a pop-up window.

    If you are using my simulated exceptions or have set the disable exceptions option in include.h then uncaught exceptions automatically print the @@ -767,13 +767,13 @@

    2.6 Catching exceptions

    2.7 Examples

    next - skip - up - start

    -

    I include a number of example files. See the sections on make -files and on compilers for information about compiling +

    I include a number of example files. See the sections on make +files and on compilers for information about compiling them.

    -

    Invert matrix: nm_ex1.cpp. Load values into a 4x4 matrix; invert it +

    Invert matrix: nm_ex1.cpp. Load values into a 4x4 matrix; invert it and check the result. The output is in nm_ex1.txt.

    -

    Eigenvalues and eigenvectors of Hilbert matrix: nm_ex2.cpp. Calculate the -eigenvalues and eigenvectors of a 7x7 Hilbert matrix. The output is in +

    Eigenvalues and eigenvectors of Hilbert matrix: nm_ex2.cpp. Calculate the +eigenvalues and eigenvectors of a 7x7 Hilbert matrix. The output is in nm_ex2.txt.

    Values in precisio.h: nm_ex3.cpp.

    Linear regression example: example.cpp. This gives a linear @@ -784,7 +784,7 @@

    2.7 Examples

    and the use of the DO_FREE_CHECK option.

    Other example files are nl_ex.cpp and garch.cpp for demonstrating the non-linear fitting routines, sl_ex for demonstrating -the solve function and test_exc for demonstrating and testing exception +the solve function and test_exc for demonstrating and testing exception handling.

    2.8 Testing

    next - skip - @@ -815,8 +815,8 @@

    2.8 Testing

    different sized items in different places, or might not allocate items consecutively. Or they might mix the items with memory blocks from other programs. Nevertheless, I seem to get consistent answers from some of the -compilers I work with, so I think this is a worthwhile test. The compilers that -the test seems to work for include the Borland compilers, Microsoft VC++ 6 , +compilers I work with, so I think this is a worthwhile test. The compilers that +the test seems to work for include the Borland compilers, Microsoft VC++ 6 , Watcom, and Gnu 2.96 for Linux.

    If the DO_FREE_CHECK option in include.h is activated, the program checks that each new is balanced with exactly one @@ -844,41 +844,41 @@

    2.9 Bugs

    2.10 Problem areas

    next - skip - up - start

    -

    This section lists parts of Newmat which users (including me) have found -difficult or unnatural. Also see the newmat FAQ list on +

    This section lists parts of Newmat which users (including me) have found +difficult or unnatural. Also see the newmat FAQ list on my website.

    Invert, element access, matrix multiply etc causes the program to crash

    -

    Newmat throws an exception when it detects an error. This can cause a program +

    Newmat throws an exception when it detects an error. This can cause a program crash unless the exception is caught with a catch statement. See catching exceptions.

    1x1 matrix not automatically converted to a Real

    Use the as_scalar() member function or the -dotproduct() function to take the dot +dotproduct() function to take the dot product of two vectors.

    Constructors do not initialise elements

    -

    For example, Matrix A(4,5); does not initialise the elements of A. Use the +

    For example, Matrix A(4,5); does not initialise the elements of A. Use the statement A=0.0 to set the values to zero.

    resize does not initialise elements

    -

    For example, A.resize(5,6); does not set the elements of A. +

    For example, A.resize(5,6); does not set the elements of A. If you want to keep values use resize_keep. See resize.

    Setting Matrix to a scalar sets all the values

    -

    A(1,3) = 0.0; sets one element of a Matrix to zero. A = 0.0; sets all the -elements to zero. This is very convenient but also a source of error that is - +

    A(1,3) = 0.0; sets one element of a Matrix to zero. A = 0.0; sets all the +elements to zero. This is very convenient but also a source of error that is + hard to see if you wanted A(1,3) = 0.0; but left out the element details.

    Symmetry not detected automatically

    -

    For example, SymmetricMatrix SM = A.t() * A; will fail. Use SymmetricMatrix +

    For example, SymmetricMatrix SM = A.t() * A; will fail. Use SymmetricMatrix SM; SM << A.t() * A;

    << does not work with constructors

    @@ -887,7 +887,7 @@

    2.10 Problem areas

    Multiple multiplication may be inefficient

    -

    For example, A * B * X where A and B are matrices and X is a column vector is +

    For example, A * B * X where A and B are matrices and X is a column vector is likely to be much slower than A * (B * X).

     

    @@ -1381,7 +1381,7 @@

    3. Reference manual

    3.17 Output
    3.18 Accessing unspecified type - + 3.19 Cholesky decomposition
    3.20 QR decomposition
    3.21 Singular value decomposition
    @@ -1422,7 +1422,7 @@

    3.1 Constructors

    DiagonalMatrix D(n); -

    Band matrices need to include bandwidth information in their constructors. +

    Band matrices need to include bandwidth information in their constructors. 

        BandMatrix BM(n, lower, upper);
     UpperBandMatrix UB(n, upper);
@@ -1432,7 +1432,7 @@ 
    3.1 Constructors
    
 
 
    The integers upper and lower are the number of non-zero
 diagonals above and below the diagonal (excluding the diagonal)
-respectively.  The UpperBandMatrix and LowerBandMatrix are upper and lower 
+respectively.  The UpperBandMatrix and LowerBandMatrix are upper and lower
 triangular band matrices. So an UpperBandMatrix is essentially a BandMatrix with
 lower = 0 and a LowerBandMatrix is a BandMatrix with upper = 0.
    
 
    The RowVector and ColumnVector types take just one argument in their
@@ -1447,10 +1447,10 @@ 
    3.1 Constructors
    
 
        Matrix A(m, n); A = 0.0;
 
    
 
-
    The IdentityMatrix takes one argument in its constructor specifying its 
+
    The IdentityMatrix takes one argument in its constructor specifying its
 dimension.
    
 
        IdentityMatrix I(n);
    
-
    The value of the diagonal elements is set to 1 by default, but you can 
+
    The value of the diagonal elements is set to 1 by default, but you can
 change this value as with other matrix types. 
    
 
 
    You can also construct vectors and matrices without specifying the
@@ -1468,7 +1468,7 @@ 
    3.1 Constructors
    
 
 
    Only conversions that don't lose information are supported - eg you cannot
 convert an upper triangular matrix into a diagonal matrix using =. 
    
-
    3.2 Accessing elements 
+
    3.2 Accessing elements
 
    
 
    next - skip -
 up - start
    
@@ -1496,7 +1496,7 @@ 
    3.2 Accessing elements
 think of the storage as being in the upper triangle (but it does matter
 in some other situations such as when entering data).
    
 
    The IdentityMatrix type does not support element access.
    
-
    For interfacing with traditional C functions that involve one and two 
+
    For interfacing with traditional C functions that involve one and two
 dimensional arrays see accessing C functions.
    
 
    3.3 Assignment and copying
    
 
    next - skip -
@@ -1559,24 +1559,24 @@ 
    3.3 Assignment and copying
    
     Matrix A(m,n); A = r;
    -To swap the values in two matrices A and B use one of the +To swap the values in two matrices A and B use one of the following expressions
       A.swap(B);
       swap(A,B);
    -The matrices A and B must be of the same type. This can be any +The matrices A and B must be of the same type. This can be any of the matrix types listed in the section on constructors, CroutMatrix, BandLUMatrix or -GenericMatrix. Swap works by switching pointers and does not do any actual copying +GenericMatrix. Swap works by switching pointers and does not do any actual copying of the main data arrays.

    Notes:

      -
    • When you do a matrix assignment to another matrix or matrix expression with -either = or << the original data array associated with the matrix - being assigned to is destroyed -even if there is no change in length. See the section on storage. +
    • When you do a matrix assignment to another matrix or matrix expression with +either = or << the original data array associated with the matrix + being assigned to is destroyed +even if there is no change in length. See the section on storage. This means, in particular, that pointers to matrix elements - e.g. -Real* a; a = &(A(1,1)); become invalid. If you want avoid this you - can use -Inject rather than =. But remember that you may need to zero +Real* a; a = &(A(1,1)); become invalid. If you want avoid this you + can use +Inject rather than =. But remember that you may need to zero the matrix first.
    @@ -1649,7 +1649,7 @@

    3.5 Unary operators

    RowVector X = A.sum_square_columns(); // sum of squares of // elements of each column -

    See the function summary list for the older +

    See the function summary list for the older depreciated function name.

    @@ -1706,34 +1706,34 @@

    3.6 Binary operators

    submatrices.
  • See the section on submatrices on using a submatrix on the RHS of an expression.
    • -
  • crossproduct - assumes A and B are both RowVectors of length 3 or both -ColumnVectors of length 3. Result is a Matrix of same dimension as A or B (will -automatically convert to RowVector if A and B are RowVectors and ColumnVector if +
  • crossproduct - assumes A and B are both RowVectors of length 3 or both +ColumnVectors of length 3. Result is a Matrix of same dimension as A or B (will +automatically convert to RowVector if A and B are RowVectors and ColumnVector if they are ColumnVectors).
    • -
  • crossproduct_rows - assumes A and B are of type Matrix with the same number +
  • crossproduct_rows - assumes A and B are of type Matrix with the same number of rows and 3 columns. Does a cross product on corresponding pairs of rows.
    • -
  • crossproduct_columns - assumes A and B are of type Matrix with  the same -number of columns and 3 rows. Does a cross product on corresponding pairs of +
  • crossproduct_columns - assumes A and B are of type Matrix with  the same +number of columns and 3 rows. Does a cross product on corresponding pairs of columns.
  • Two matrices are equal if their difference is zero. They may be of different types. For the CroutMatrix or BandLUMatrix they must be of the same type and have all their elements equal. This is not a very useful operator and is included for compatibility with some container templates.
    • -
  • The inequality operators are included for compatibility with some versions +
  • The inequality operators are included for compatibility with some versions of the standard template library. If actually called, they will throw an exception. So don't try to sort a list of matrices.
  • A row vector multiplied by a column vector yields a 1x1 matrix, not a Real. To get a Real result use either as_scalar or dotproduct.
    • -
  • The result from Kronecker product, KP(A, B), possesses an attribute such as -upper triangular, lower triangular, band, symmetric, diagonal if both of the -matrices A and B have the attribute. If A is band and B is a square type -(eg SquareMatrix, Diagonal, LowerTriangularMatrix etc) then +
  • The result from Kronecker product, KP(A, B), possesses an attribute such as +upper triangular, lower triangular, band, symmetric, diagonal if both of the +matrices A and B have the attribute. If A is band and B is a square type +(eg SquareMatrix, Diagonal, LowerTriangularMatrix etc) then the result is band.
    • -
  • Elementwise product is also known as the Schur product or the Hadamard +
  • Elementwise product is also known as the Schur product or the Hadamard product.
    • -
  • See the function summary list for the older +
  • See the function summary list for the older depreciated function names.

    • Remember that the product of symmetric matrices is not @@ -1755,14 +1755,14 @@

    3.7 Matrix and scalar

    up - start

    The following expressions multiply the elements of a matrix A by a scalar f: A * f or f * A . Likewise one can divide -the elements of a matrix A by a scalar f: A / f . +the elements of a matrix A by a scalar f: A / f .

    The expressions A + f and A - f add or subtract a rectangular matrix of the same dimension as A with elements equal to f to or from the matrix A .

    The expression f + A is an alternative to A + f. The expression f - A subtracts matrix A from a rectangular matrix -of the same dimension as A and with elements equal to f . +of the same dimension as A and with elements equal to f .

    The expression A += f replaces A by A + f. Operators -=, *=, /= are similarly defined.

    @@ -1788,10 +1788,10 @@

    3.8 Scalar functions of a mt.value() to get a string (UT, LT, Rect, Sym, Diag, Band, UB, LB, Crout, BndLU) showing the type (Vector types are returned as Rect).

    MatrixBandWidth has member functions upper() and -lower() for finding the upper and lower bandwidths (number of diagonals -above and below the diagonal, both zero for a diagonal matrix). For non-band +lower() for finding the upper and lower bandwidths (number of diagonals +above and below the diagonal, both zero for a diagonal matrix). For non-band matrices -1 is returned for both these values.

    -

    See the function summary list for the older +

    See the function summary list for the older depreciated function names.

    3.9 Scalar functions of a matrix - maximum & minimum

    @@ -1848,7 +1848,7 @@

    3.10 Scalar functions of a // interpreted as vectors -

    See the function summary list for the older +

    See the function summary list for the older depreciated function names.

    A.log_determinant() returns a value of type LogAndSign. If ld is of @@ -1875,7 +1875,7 @@

    3.10 Scalar functions of a

    The versions sum(A), sum_square(A), sum_absolute_value(A), trace(A), log_determinant(A), determinant(A), norm1(A), norm_infinity(A), norm_Frobenius(A) can be used in place of A.sum(), A.sum_square(), A.sum_absolute_value(), -A.trace(), A.log_determinant(), A.determinant(A), A.norm1(), A.norm_infinity(), A.norm_Frobenius(). +A.trace(), A.log_determinant(), A.determinant(A), A.norm1(), A.norm_infinity(), A.norm_Frobenius().

    3.11 Submatrices

    next - skip - @@ -1885,7 +1885,7 @@

    3.11 Submatrices

    This selects a submatrix from A. The arguments fr,lr,fc,lc are the first row, last row, -first column, last column of the submatrix with the numbering beginning at 1. +first column, last column of the submatrix with the numbering beginning at 1.

    I allow lr = fr-1 or lc = fc-1 or to indicate that a matrix of zero rows or columns is to be returned.

    @@ -1904,7 +1904,7 @@

    3.11 Submatrices

    A.column(f) // select single column -

    See the function summary list for the older +

    See the function summary list for the older depreciated function names.

    In each case f and l mean the first and last row or column @@ -1941,8 +1941,8 @@

    3.11 Submatrices

    memory.

    Do not use the += and -= operations with a submatrix of a SymmetricMatrix or BandSymmetricMatrix on the LHS.

    -

    You can't pass a submatrix (or any of its variants) as a reference -non-constant matrix in a function argument. For example, the following will not +

    You can't pass a submatrix (or any of its variants) as a reference +non-constant matrix in a function argument. For example, the following will not work:

       void YourFunction(Matrix& A);
    ...
@@ -1954,9 +1954,9 @@ 
    3.11 Submatrices
    
 
    3.12 Change dimensions
    
 
    next - skip -
 up - start
    
-
    The following operations change the dimensions of a matrix. The values of the 
-elements are lost, for resize. resize_keep keeps element 
-values and zeros the elements in the matrix with the new size that are not in 
+
    The following operations change the dimensions of a matrix. The values of the
+elements are lost, for resize. resize_keep keeps element
+values and zeros the elements in the matrix with the new size that are not in
 the matrix with the old size.
    
 
        A.resize(nrows,ncols);        // for type Matrix or nricMatrix
     A.resize(n);                  // for other types, except Band
@@ -1976,9 +1976,9 @@ 
    3.12 Change dimensions
    
 B. This includes the band-width in the case of a band matrix. It is an
 error for A to be a band matrix and B not a band matrix (or
 diagonal matrix). 
    
-
    Remember that resize destroys values. If you want to resize, but 
+
    Remember that resize destroys values. If you want to resize, but
 keep the values in the bit that is left use resize_keep.
    
-
    See the function summary list for the older 
+
    See the function summary list for the older
 depreciated function names.
    
 
 
    3.13 Change type
    
@@ -1996,7 +1996,7 @@ 
    3.13 Change type
    
 
 
    The expression A.as_scalar() is used to convert a 1 x 1 matrix to a
 scalar.
    
-
    See the function summary list for the older 
+
    See the function summary list for the older
 depreciated function names.
    
 
    3.14 Multiple matrix solve
    
 
    next - skip -
@@ -2016,7 +2016,7 @@ 
    3.14 Multiple matrix solve
    
 A. In these cases you probably want to avoid repeating the LU
 decomposition of A for each solve or determinant calculation. 
    
 
    If A is a square or symmetric matrix use 
    
-
        ColumnVector p, q;    
+
        ColumnVector p, q;
     ...
     CroutMatrix X = A;                // carries out LU decomposition
     ColumnVector Ap = X.i()*p; ColumnVector Aq = X.i()*q;
@@ -2024,7 +2024,7 @@ 
    3.14 Multiple matrix solve
    
 
    
 
 
    rather than 
    
-
        ColumnVector p, q;    
+
        ColumnVector p, q;
     ...
     ColumnVector Ap = A.i()*p; ColumnVector Aq = A.i()*q;
     LogAndSign ld = A.log_determinant();
@@ -2035,12 +2035,12 @@ 
    3.14 Multiple matrix solve
    
 
        BandLUMatrix X = A;               // carries out LU decomposition
 
    
 
-
    A CroutMatrix or BandLUMatrix can be copied and you can have a constructor 
+
    A CroutMatrix or BandLUMatrix can be copied and you can have a constructor
 with no parameters (use = to give it values). They work with
 release(), release_and_delete and
-GenericMatrix and ReturnMatrix. You 
+GenericMatrix and ReturnMatrix. You
 can't do any other manipulation apart from taking the inverse or solving with
-i(), or finding the determinant or log determinant. See the function summary list 
+i(), or finding the determinant or log determinant. See the function summary list
 for accessing the internals of a CroutMatrix or BandLUMatrix.
    
 
    You can alternatively use 
    
 
        LinearEquationSolver X = A;
@@ -2048,8 +2048,8 @@ 
    3.14 Multiple matrix solve
    
 
 
    This will choose the most appropriate decomposition of A. That is,
 the band form if A is banded; the Crout decomposition if A is
-square or symmetric and no decomposition if A is triangular or 
-diagonal. It doesn't know about positive definite matrices so won't use 
+square or symmetric and no decomposition if A is triangular or
+diagonal. It doesn't know about positive definite matrices so won't use
 Cholesky.
    
 
    3.15 Memory management
    
 
    next - skip -
@@ -2091,7 +2091,7 @@ 
    3.15 Memory management
    
     }
 
    
 
-
    If your compiler objects to this code, replace the return statements with 
+
    If your compiler objects to this code, replace the return statements with
 
    
 
        return A.for_return();
 
    
@@ -2109,7 +2109,7 @@ 
    3.15 Memory management
    
        ......
        A.release(); return A;   // will compile but could give wrong answers.
     }
    
-
    since with some compilers A might retain its released status 
+
    since with some compilers A might retain its released status
 after being returned.
    
 
    You can also use .release() or ->release_and_delete() to
 allow a matrix expression to recycle space. Suppose you call 
    
@@ -2122,7 +2122,7 @@ 
    3.15 Memory management
    
 can be used to improve efficiency and reduce the use of memory. 
    
 
    Use ->release_and_delete for matrices created by new if you want
 to completely delete the matrix after it is accessed.
    
-
    See the function summary list for the older 
+
    See the function summary list for the older
 depreciated function names.
    
 
    3.16 Efficiency
    
 
    next - skip -
@@ -2148,7 +2148,7 @@ 
    3.17 Output
    
 
    
 
 
    This will work only with systems that support the standard input/output
-routines including manipulators. The scientific and fixed manipulators won't 
+routines including manipulators. The scientific and fixed manipulators won't
 work with Visual C++, version 6.
    
 
    You need to #include the files iostream.h,
 iomanip.h, newmatio.h in your C++ source files that use this facility. The
@@ -2186,11 +2186,11 @@ 
    3.18 Unspecified type
    
 
    
 
 
    A GenericMatrix matrix can be used anywhere where a matrix expression can be
-used and also on the left hand side of an =. You can pass any type of 
+used and also on the left hand side of an =. You can pass any type of
 matrix to a const
 GenericMatrix& argument in a function. However most scalar functions
 including nrows(), ncols(), type() and element access do not work with it. Nor
-does the ReturnMatrix construct. swap does work with objects of type 
+does the ReturnMatrix construct. swap does work with objects of type
 GenericMatrix. See also the paragraph on LinearEquationSolver. 
    
 
    An alternative and less flexible approach is to use BaseMatrix or
@@ -2235,13 +2235,13 @@ 
    3.19 Cholesky
 
    
 
 
    If S is a symmetric band matrix then L is a band matrix
-and an alternative procedure is provided for carrying out the decomposition: 
+and an alternative procedure is provided for carrying out the decomposition:
 
    
 
        SymmetricBandMatrix S;
     ......
     LowerBandMatrix L = Cholesky(S);
    
 
-
    See section 3.32 on updating a Cholesky 
+
    See section 3.32 on updating a Cholesky
 decomposition.
    
 
    
 
@@ -2264,7 +2264,7 @@ 
    3.20 QR decomposition
    
          s    t
    -

    where U is upper triangular (the R of the QR transform). That is +

    where U is upper triangular (the R of the QR transform). That is 

          Q  / 0   0 \  =  / U   M \
          \ X   Y /     \ 0   Z /
    @@ -2299,12 +2299,12 @@

    3.20 QR decomposition

    The are also routines QRZT(X, L), QRZT(X, Y, M) and QRZT(X, Y, L, M) which do the same decomposition on the transposes of all these matrices. QRZT -replaces the routines HHDecompose in earlier versions of newmat. +replaces the routines HHDecompose in earlier versions of newmat. HHDecompose is still defined but just calls QRZT.

    For an example of the use of this decomposition see the file example.cpp.

    -

    See the section on updating a Cholesky decomposition +

    See the section on updating a Cholesky decomposition for updating U.

    Alternatively to update U or L with a new block of data, X, one can use

        updateQRZ(X, U);
@@ -2313,10 +2313,10 @@ 
    3.20 QR decomposition
    
 
        updateQRZT(X, L);
    
 
    You can update the M matrix with
    
 
        updateQRZ(X, Y, M);
    
-
    At present, there is no corresponding function for updateQRZT. Also 
+
    At present, there is no corresponding function for updateQRZT. Also
 note, at present, updateQRZ(X, Y, M) is not optimised for accessing
 contiguous locations.
    
-
    If you have carried out QRZ transforms on two blocks of data, X1, X2, 
+
    If you have carried out QRZ transforms on two blocks of data, X1, X2,
 (same number of columns in each) with the data to be fitted in Y1, 
 Y2, you can combine them as follows
    
 
        // QRZ transforms of two blocks of data
@@ -2327,28 +2327,28 @@ 
    3.20 QR decomposition
    
     updateQRZ(U1, U2); updateQRZ(U1, M1, M2);
    
 
    Results are in U2 and M2. The values in U1 and 
 M1 are modified so use copies if you still need the original values.
    
-
    The functions updateQRZ(U1, U2) and updateQRZ(U1, M1, M2) 
-are not optimised for accessing contiguous locations. 
-This is probably not an issue, since usually all the data will fit into cache 
+
    The functions updateQRZ(U1, U2) and updateQRZ(U1, M1, M2)
+are not optimised for accessing contiguous locations.
+This is probably not an issue, since usually all the data will fit into cache
 memory.  
    
 
    Notes on updateQRZ and updateQRZT:
    
 
      
-  
    • Use these routines if your data is arriving in blocks or there is too much to 
+  
    • Use these routines if your data is arriving in blocks or there is too much to
 fit into memory.
      • 
   
    • Use QRZ or QRZT on the first block of data or you can use 
-  updateQRZ or 
+  updateQRZ or
 updateQRZT if U or L is correctly dimensioned and set to zero.
      • 
-  
    • The block of data, X, will be modified by these routines but the 
+  
    • The block of data, X, will be modified by these routines but the
 modified data is probably not useful for further analysis.
      • 
-  
    • The signs of some of the rows of U or columns of L maybe the reverse of what you would 
-expect (i.e. you get the negative of what you get from QRZ or QRZT applied to 
+  
    • The signs of some of the rows of U or columns of L maybe the reverse of what you would
+expect (i.e. you get the negative of what you get from QRZ or QRZT applied to
 	the whole dataset). Usually this will not be a problem.
      • 
-	
    • The second version of updateQRZ is sometimes a more useful way of 
-	combining blocks, especially if you want to use the reduce function 
-	in Intel's Threaded Building Blocks to parallelise QRZ. (Contact me 
-	for what you need to do to make Newmat compatible with Threaded Building 
+	
    • The second version of updateQRZ is sometimes a more useful way of
+	combining blocks, especially if you want to use the reduce function
+	in Intel's Threaded Building Blocks to parallelise QRZ. (Contact me
+	for what you need to do to make Newmat compatible with Threaded Building
 	Blocks).
      • 
-  
    • See the function summary list for the older 
+  
    • See the function summary list for the older
 depreciated function names.
      • 
 
    
 
     
    
@@ -2390,8 +2390,8 @@ 
    3.21 Singular value
 number of the singular values are identical one can apply an orthogonal
 transformation to the corresponding columns of U and the corresponding
 columns of V.
    
-
    If m < n apply the SVD transform to the transpose 
-of A and swap U and V. If necessary, extend U to a square matrix as described 
+
    If m < n apply the SVD transform to the transpose
+of A and swap U and V. If necessary, extend U to a square matrix as described
 above.
    
 
    3.22 Eigenvalue
 decomposition
    
@@ -2435,7 +2435,7 @@ 
    3.22 Eigenvalue
 
    Remember that an eigenvalue decomposition is not completely
 unique - see the comments about the SVD
 decomposition.
    
-
    See the function summary list for the older 
+
    See the function summary list for the older
 depreciated function names. 
    
 
    3.23 Sorting
    
 
    next - skip -
@@ -2452,9 +2452,9 @@ 
    3.23 Sorting
    
 can do) an exception is thrown. 
    
 
    You will get incorrect results if you try to sort a band matrix - but why
 would you want to sort a band matrix? 
    
-
    See the function summary list for the older 
+
    See the function summary list for the older
 depreciated function names.
    
-
    3.24 Fast Fourier transform 
+
    3.24 Fast Fourier transform
 
    
 
    next - skip -
 up - start
    
@@ -2476,7 +2476,7 @@ 
    3.24 Fast Fourier transform
 Y(j+1). Likewise h[k] is complex and stored in
 F(k+1) and G(k+1). The fast Fourier algorithm takes order 
 n log(n) operations (for good values of n) rather than
-n**2 that straight evaluation (see the file tmtf.cpp) takes. 
+n**2 that straight evaluation (see the file tmtf.cpp) takes.
 
    
 
    I use one of two methods: 
    
 
      
@@ -2507,7 +2507,7 @@ 
      3.24 Fast Fourier transform
    FFT2I(F, G, X, Y);                      // inverse, X=F and Y=G are OK

    where X, Y, F, G are of type Matrix. X and Y are the real and imaginary input matrices; F -and G are the real and imaginary output matrices. The dimensions of Y +and G are the real and imaginary output matrices. The dimensions of Y must be the same as those of X and should be the product of numbers less than about 10 for fast execution.

    The formula is

    @@ -2518,7 +2518,7 @@

    3.24 Fast Fourier transform

    where z[j,k] is stored complex and stored in X(j+1,k+1) and Y(j+1,k+1) and X and Y have dimension m x n. Likewise h[p,q] is complex and stored in -F(p+1,q+1) and G(p+1,q+1). +F(p+1,q+1) and G(p+1,q+1).

     

    3.25 Fast trigonometric @@ -2589,8 +2589,8 @@

    3.26 Interface to Numerical Recipes in C

    next - skip - up - start

    -

    This section refers to Numerical Recipes in C. This section is not -relevant to Numerical Recipes in C++. I'll put a note on the website soon +

    This section refers to Numerical Recipes in C. This section is not +relevant to Numerical Recipes in C++. I'll put a note on the website soon on how to interface with Numerical Recipes in C++.

    This package can be used with the vectors and matrices defined in Numerical Recipes in C. You need to edit the routines in Numerical @@ -2658,10 +2658,10 @@

    3.27 Exceptions

    Tracer class. Also see the section on error messages for additional notes on the messages generated by the exceptions.

    -

    You can also print this information in a catch clause by printing Exception::what(). +

    You can also print this information in a catch clause by printing Exception::what().

    If you are using compiler supported exceptions then see the section -on catching exceptions.  +on catching exceptions. 

    See the file test_exc.cpp as an example of catching an exception and printing the error message.

    @@ -2672,10 +2672,10 @@

    3.27 Exceptions

    Catch or CatchAll block to determine the action.

    The library includes the alternatives of using the inbuilt exceptions provided by a compiler, simulating exceptions, or disabling exceptions. See -customising for selecting the correct exception option. +customising for selecting the correct exception option.

    The rest of this section describes my partial simulation of exceptions for -compilers which do not support C++ exceptions. Skip the rest of this section and +compilers which do not support C++ exceptions. Skip the rest of this section and the next section if you are using compiler supported exceptions. I use Carlos Vidal's article in the September 1992 C Users Journal as a starting point.

    Newmat does a partial clean up of memory following throwing an exception - @@ -2729,7 +2729,7 @@

    3.28 Cleanup after an

    To overcome this problem newmat keeps a list of all the currently declared matrices and its exception mechanism will return heap memory when you do a Throw and Catch.

    -

    However it will not return heap memory from objects from other packages. +

    However it will not return heap memory from objects from other packages.

    If you want the mechanism to work with another class you will have to do four things:

    @@ -2750,7 +2750,7 @@

    3.28 Cleanup after an -

    Use CleanUp() rather than cleanup() since this is what is +

    Use CleanUp() rather than cleanup() since this is what is defined in class Janitor.

    Note that the function CleanUp() does somewhat the same duties as the destructor. However CleanUp() has to do the cleaning for @@ -2769,13 +2769,13 @@

    3.29 Non-linear in the non-linear least squares.

    Files newmatnl.h, newmatnl.cpp contain a series of classes for non-linear least squares and maximum likelihood. These classes work on very well-behaved -functions but need upgrading for less well-behaved functions. I haven't followed -the usual practice of inflating the values of the diagonal elements of the -matrix of second derivatives. I originally thought I could avoid this if my -program had a good line search. But this was wrong and when I use this program -on all but the most well-behaved problems I run the fit first with the diagonal -elements inflated by a factor of 2 to 5 and the critical value for the stopping -criterion set to something like 50. Then rerun with with no inflation factor and +functions but need upgrading for less well-behaved functions. I haven't followed +the usual practice of inflating the values of the diagonal elements of the +matrix of second derivatives. I originally thought I could avoid this if my +program had a good line search. But this was wrong and when I use this program +on all but the most well-behaved problems I run the fit first with the diagonal +elements inflated by a factor of 2 to 5 and the critical value for the stopping +criterion set to something like 50. Then rerun with with no inflation factor and critical value 0.0001.

    Documentation for both of these is in the definition files. Simple examples are in sl_ex.cpp, nl_ex.cpp and garch.cpp.

    @@ -2794,7 +2794,7 @@

    3.30 Standard template
  • Make sure the option DO_FREE_CHECK is not turned on.
  • You can store only one type of matrix in a container. If you want to use a -variety of types use the GenericMatrix type or store pointers to the matrices. +variety of types use the GenericMatrix type or store pointers to the matrices.
  • The vector and deque container templates like to copy their elements. For the vector container this happens when you insert an element anywhere except at @@ -2818,7 +2818,7 @@

    3.31 Namespace

    at the beginning of any file that needs to access all my libraries.

    -

    This works correctly with Borland C++ version 5 and +

    This works correctly with Borland C++ version 5 and Builder 5 and 6.

    Microsoft Visual C++ version 5 works in my example and test files, but fails with apparently insignificant changes (it may be more @@ -2826,7 +2826,7 @@

    3.31 Namespace

    "newmatap.h", but no other newmat include file, then also #include "newmatio.h". It seems to work with Microsoft Visual C++ version 6 if you have applied at least service pack 2.

    -

    My use of namespace does not work with Gnu g++ version +

    My use of namespace does not work with Gnu g++ version 2.8.1 but does work with versions 3.x.

    I have defined the following namespaces:

      @@ -2850,8 +2850,8 @@

      3.32 Updating the Cholesky matrix

         UpperTriangularMatrix U;
    QRZ(X, U);

      See sections 3.19 and 3.20.

      -

      Now suppose we want to append an extra row to X or delete a row from X -or rearrange the columns of X. The functions described here allow you to +

      Now suppose we want to append an extra row to X or delete a row from X +or rearrange the columns of X. The functions described here allow you to update U without  recalculating it.

      @@ -2859,33 +2859,33 @@

      3.32 Updating the Cholesky matrix

      downdate_Cholesky(U, x); // x is a RowVector right_circular_update_Cholesky(U, j, k); // j and k are integers left_circular_update_Cholesky(U, j, k); // j and k are integers -

      update_Cholesky carries out the modification of U corresponding to appending +

      update_Cholesky carries out the modification of U corresponding to appending an extra row x to X.

      -

      downdate_Cholesky carries out the modification corresponding -to removing a row x from X. A ProgramException exception is +

      downdate_Cholesky carries out the modification corresponding +to removing a row x from X. A ProgramException exception is thrown if the modification fails.

      -

      right_circular_update_Cholesky supposes that columns j,j+1,...k of X are +

      right_circular_update_Cholesky supposes that columns j,j+1,...k of X are replaced by columns k,j,j+1,...k-1.

      -

      left_circular_update_Cholesky supposes that columns j,j+1,...k of X are replaced +

      left_circular_update_Cholesky supposes that columns j,j+1,...k of X are replaced by columns j+1,...k,j.

      -

      These functions are based on a contribution from Nick Bennett of -Schlumberger-Doll Research. See also the LINPACK User's Guide, Chapter 10, +

      These functions are based on a contribution from Nick Bennett of +Schlumberger-Doll Research. See also the LINPACK User's Guide, Chapter 10, Dongarra et. al., SIAM, Philadelphia, 1979.

      -

      Where you want to append a number of new rows consider using the update +

      Where you want to append a number of new rows consider using the update routine in section 3.20.

      -

      See the function summary list for the older +

      See the function summary list for the older depreciated function names.

      @@ -2896,21 +2896,21 @@

      3.33 Accessing C functions

      up - start

      -

      You have a C function that uses one and two dimensional arrays as +

      You have a C function that uses one and two dimensional arrays as vectors and matrices. You want to access it from Newmat.

      -

      One dimensional arrays are easy. Set up a ColumnVector, RowVector +

      One dimensional arrays are easy. Set up a ColumnVector, RowVector or DiagonalMatrix with the correct dimension and where the function has a double* argument enter X.data() where X denotes the -ColumnVector, RowVector or DiagonalMatrix. (I am assuming you +ColumnVector, RowVector or DiagonalMatrix. (I am assuming you have left Real being typedefed as a double). If you have a const double* argument use X.const_data().

      -

      You can't do this with two dimensional arrays where you have a double** -argument. Newmat includes classes RealStarStar and ConstRealStarStar for this situation. -To access a Matrix A from a function c_function(double** a) +

      You can't do this with two dimensional arrays where you have a double** +argument. Newmat includes classes RealStarStar and ConstRealStarStar for this situation. +To access a Matrix A from a function c_function(double** a) use either

      @@ -2921,15 +2921,15 @@

      3.33 Accessing C functions

      If the argument is const double** use ConstRealStarStar.

      -

      The Matrix A must be correctly dimensioned and must not be -resized or set equal to another matrix between setting up the RealStarStar -object and its use in the function. Also the following construction will not +

      The Matrix A must be correctly dimensioned and must not be +resized or set equal to another matrix between setting up the RealStarStar +object and its use in the function. Also the following construction will not work

         double** a = RealStarStar(A);      // wrong
    c_function(a);
      -

      since the RealStarStar structure will have been destroyed before you +

      since the RealStarStar structure will have been destroyed before you get to the second line.

      @@ -2940,8 +2940,8 @@

      3.34 Simple integer array class

      up - start

      -

      This is primarily for use within Newmat. You can set up a simple array -of integers with the SimpleIntArray class. Here are the descriptions of +

      This is primarily for use within Newmat. You can set up a simple array +of integers with the SimpleIntArray class. Here are the descriptions of the constructors and functions.

      @@ -2952,18 +2952,18 @@

      3.34 Simple integer array class

      SimpleIntArray A(n); - Constructs int array of length n - + Constructs int array of length n - individual elements are not initialised A = i; - sets values to i where i is an int + sets values to i where i is an int variable A = B; - sets values to those of B where -B is + sets values to those of B where +B is a SimpleIntArray, change size if necessary. @@ -2972,17 +2972,17 @@

      3.34 Simple integer array class

      A.resize(n); - change the length of A, don't keep + change the length of A, don't keep values A.resize_keep(n); - change the length of A, do keep + change the length of A, do keep values; if length is being increased set new elements to zero. int x = A[i]; - element access; i runs from 0 to + element access; i runs from 0 to n-1 @@ -2991,12 +2991,12 @@

      3.34 Simple integer array class

      const int* d = A.data(); - get beginning of data array as const + get beginning of data array as const int* if A is const const int* d = A.const_data(); - get beginning of data array as const + get beginning of data array as const int* @@ -3010,12 +3010,12 @@

      3.34 Simple integer array class

      Extend orthonormal set of columns

    • next - skip - up - start

    -

    Suppose a Matrix A's first n columns are orthonormal so that +

    Suppose a Matrix A's first n columns are orthonormal so that A.Columns(1,n).t() * A.Columns(1,n) is the identity matrix. Suppose -we want to fill out the remaining columns of A to make them orthonormal so +we want to fill out the remaining columns of A to make them orthonormal so that A.t() * A is the identity matrix. Then use the function

       extend_orthonormal(A, n);
    -

    Matrix A is then replaced by the matrix with the additional +

    Matrix A is then replaced by the matrix with the additional columns.

    Use this function to extend U from the QRZ or SVD decompositions to form a square (orthogonal) matrix.

    @@ -3023,19 +3023,19 @@

    3.34 Simple integer array class

    • the first n columns of A must be orthonormal
    • n must be less or equal to the number of columns of A
      • -
    • the number of columns of A must be less than or equal to the number of +
    • the number of columns of A must be less than or equal to the number of rows of A

    3.36 Miscellaneous functions

    next - skip - up - start

    -

    The section includes some miscellaneous functions I have needed for my work. So +

    The section includes some miscellaneous functions I have needed for my work. So far there are only the Helmert transforms.

    Helmert transforms

    -

    This section refers to the Helmert transform used in some statistical -packages for extracting contrasts. It is different from the Helmert transform -used in geodesy. The version of the transform I am going to use is based on the +

    This section refers to the Helmert transform used in some statistical +packages for extracting contrasts. It is different from the Helmert transform +used in geodesy. The version of the transform I am going to use is based on the n x n matrix

    | 1/sqrt(1*2)                       | * | -1  1              |
 |     1/sqrt(2*3)                   |   | -1 -1  2           |
@@ -3043,14 +3043,14 @@ 
    Helmert transforms
    
 |               ...                 |   |  ...               |
 |                  1/(sqrt(n-1)*n)  |   | -1 -1  ...  -1 n-1 |
 |                         1/sqrt(n) |   |  1  1  ...       1 |
    -

    You can interpret multiplying a column vector by this matrix as follows. Form -the k-th element in the resulting column vector by taking the k+1-th -element and subtracting the average of the previous k elements. Then -multiply by sqrt((k+1)/k) so we get an orthonormal transform. The -last element is just the sum divided by sqrt(n). Usually one will omit +

    You can interpret multiplying a column vector by this matrix as follows. Form +the k-th element in the resulting column vector by taking the k+1-th +element and subtracting the average of the previous k elements. Then +multiply by sqrt((k+1)/k) so we get an orthonormal transform. The +last element is just the sum divided by sqrt(n). Usually one will omit the last element since we want just the contrasts.

    -

    Here are the functions. X and Y are ColumnVectors, A -and B are matrices, b is a boolean, j, n are +

    Here are the functions. X and Y are ColumnVectors, A +and B are matrices, b is a boolean, j, n are integers, x is a Real.

    @@ -3058,7 +3058,7 @@

    Helmert transforms

    A = Helmert(n, b);
    A = Helmert(n); - Return the n x n + Return the n x n Helmert matrix or the (n-1) x n version with the last row omitted. @@ -3069,36 +3069,36 @@

    Helmert transforms

    Y = Helmert(n,j,b);
    Y = Helmert(n,j); - Return j-th column of Helmert + Return j-th column of Helmert matrix as a ColumnVector. X = Helmert_transpose(Y,b);
    X = Helmert_transpose(Y); - Multiply by the transpose of the + Multiply by the transpose of the Helmert matrix. x = Helmert_transpose(Y,j,b);
    x = Helmert_transpose(Y,j); - Multiply by the transpose of the + Multiply by the transpose of the Helmert matrix, return just the j-th element. B = Helmert(A,b);
    B = Helmert(A); - Multiply a Matrix by the Helmert + Multiply a Matrix by the Helmert matrix. A = Helmert_transpose(B,b);
    A = Helmert_transpose(B); - Multiply a Matrix by the transpose of + Multiply a Matrix by the transpose of the Helmert matrix. -

    If b is true the full Helmert matrix is used, if it is false or +

    If b is true the full Helmert matrix is used, if it is false or omitted the n-th row (or n-column of the transpose) is omitted.

     

    @@ -3181,7 +3181,7 @@

    Some general comments

    By safety, I mean getting the right answer, and not causing crashes or damage to the computer system.

    By usability, I mean being easy to learn and use, including not being -too complicated, being intuitive, saving the users' time, being nice to use. +too complicated, being intuitive, saving the users' time, being nice to use.

    Efficiency means minimising the use of computer memory and time.

    In the early days of computers the emphasis was on efficiency. But computer @@ -3238,53 +3238,53 @@

    Levels of access

    safety is reduced but you can get better efficiency.

    Some performance data

    This section looks at the performance on newmat for simple sums, -comparing it with C code and with a simple array program. +comparing it with C code and with a simple array program. (This is now very old data).

    The following table lists the time (in seconds) for carrying out the operations X=A+B;, X=A+B+C;, X=A+B+C+D;, X=A+B+C+D+E; where X,A,B,C,D,E are of type ColumnVector, with a variety of programs. I am -using Microsoft VC++, version 6 in console mode under windows 2000 on a PC with +using Microsoft VC++, version 6 in console mode under windows 2000 on a PC with a 1 ghz Pentium III and 512 mbytes of memory. 

        length    iters. newmat      C C-res.  subs.  array
 X = A + B
-         2   5000000   27.8    0.3    8.8    1.9    9.5 
-        20    500000    3.0    0.3    1.1    1.9    1.2 
-       200     50000    0.5    0.3    0.4    1.9    0.3 
-      2000      5000    0.4    0.3    0.4    2.0    1.0 
-     20000       500    4.5    4.5    4.5    6.7    4.4 
-    200000        50    5.2    4.7    5.5    5.8    5.2 
+         2   5000000   27.8    0.3    8.8    1.9    9.5
+        20    500000    3.0    0.3    1.1    1.9    1.2
+       200     50000    0.5    0.3    0.4    1.9    0.3
+      2000      5000    0.4    0.3    0.4    2.0    1.0
+     20000       500    4.5    4.5    4.5    6.7    4.4
+    200000        50    5.2    4.7    5.5    5.8    5.2
 
 X = A + B + C
-         2   5000000   36.6    0.4    8.9    2.5   12.2 
-        20    500000    4.0    0.4    1.2    2.5    1.6 
-       200     50000    0.8    0.3    0.5    2.5    0.5 
-      2000      5000    3.6    4.4    4.6    9.0    4.4 
-     20000       500    6.8    5.4    5.4    9.6    6.8 
-    200000        50    8.6    6.0    6.7    7.1    8.6 
+         2   5000000   36.6    0.4    8.9    2.5   12.2
+        20    500000    4.0    0.4    1.2    2.5    1.6
+       200     50000    0.8    0.3    0.5    2.5    0.5
+      2000      5000    3.6    4.4    4.6    9.0    4.4
+     20000       500    6.8    5.4    5.4    9.6    6.8
+    200000        50    8.6    6.0    6.7    7.1    8.6
 
 X = A + B + C + D
-         2   5000000   44.0    0.7    9.3    3.1   14.6 
-        20    500000    4.9    0.6    1.5    3.1    1.9 
-       200     50000    1.0    0.6    0.8    3.2    0.8 
-      2000      5000    5.6    6.6    6.8   11.5    5.9 
-     20000       500    9.0    6.7    6.8   11.0    8.5 
-    200000        50   11.9    7.1    7.9    9.5   12.0 
+         2   5000000   44.0    0.7    9.3    3.1   14.6
+        20    500000    4.9    0.6    1.5    3.1    1.9
+       200     50000    1.0    0.6    0.8    3.2    0.8
+      2000      5000    5.6    6.6    6.8   11.5    5.9
+     20000       500    9.0    6.7    6.8   11.0    8.5
+    200000        50   11.9    7.1    7.9    9.5   12.0
 
 X = A + B + C + D + E
-         2   5000000   50.6    1.0    9.5    3.8   17.1 
-        20    500000    5.7    0.8    1.7    3.9    2.4 
-       200     50000    1.3    0.9    1.0    3.9    1.0 
-      2000      5000    7.0    8.3    8.2   13.8    7.1 
-     20000       500   11.5    8.1    8.4   13.2   11.0 
-    200000        50   15.2    8.7    9.5   12.4   15.4 
+         2   5000000   50.6    1.0    9.5    3.8   17.1
+        20    500000    5.7    0.8    1.7    3.9    2.4
+       200     50000    1.3    0.9    1.0    3.9    1.0
+      2000      5000    7.0    8.3    8.2   13.8    7.1
+     20000       500   11.5    8.1    8.4   13.2   11.0
+    200000        50   15.2    8.7    9.5   12.4   15.4
    -

    I have graphed the results +

    I have graphed the results and included rather more array lengths.

    The first column gives the lengths of the arrays, the second the number of iterations and the remaining columns the total time required in seconds. If the only thing that consumed time was the double precision addition then the -numbers within each block of the table would be the same. The summation is +numbers within each block of the table would be the same. The summation is repeated 5 times within each loop, for example:

       for (i=1; i<=m; ++i)
    {
@@ -3294,30 +3294,30 @@ 
    Some performance data
    
 
    The column labelled newmat is using the standard newmat add. The column labelled C uses the usual C method: while
 (j1--) *x1++ = *a1++ + *b1++; . The following column also includes an
 X.resize() in the outer loop to correspond to the reassignment of
-memory that newmat would do. In the next column the calculation is using 
+memory that newmat would do. In the next column the calculation is using
 the usual C style for loop
 and accessing the elements using newmat subscripts such as
 A(i). The final column is the time taken by a
-simple array package. This uses an alternative method for avoiding temporaries 
+simple array package. This uses an alternative method for avoiding temporaries
 and unnecessary copies that does not involve runtime tests. It does its sums in blocks of 4 and copies in blocks of
 8 in the same way that newmat does. 
    
 
    Here are my conclusions. 
    
 
      
 
    • Newmat does very badly for length 2 and doesn't do  well for
 length 20. There is  a lot of code in newmat for
-determining which sum algorithm to use and it is not surprising that this 
+determining which sum algorithm to use and it is not surprising that this
 impacts on performance for small lengths.
 However the array program is also having difficulty with length 2 so it
 is unlikely that the problem could be completely eliminated. 
      • 
 
    • For arrays of length 2000 or longer newmat is doing about as well as
 C and slightly better than C with resize in the X=A+B table. For the
 other two tables it tends to be slower, but not dramatically so. 
      • 
-
    • It is really important for fast processing with the Pentium III to stay 
+
    • It is really important for fast processing with the Pentium III to stay
 within the Pentium cache.
      • 
 
    • Addition using the newmat subscripts, while considerably slower than
 the others, is still surprisingly good for the longer arrays. 
      • 
 
    • The array program and newmat are similar for
-lengths 2000 or higher (the longer times for the array program for the longest 
+lengths 2000 or higher (the longer times for the array program for the longest
 arrays shown on the graph are probably a quirk of the timing program). 
      • 
 
    
 
    In summary: for the situation considered here, newmat is doing very
@@ -3454,7 +3454,7 @@ 
    What element types?
    
 
    It might also be worth including symmetric and triangular matrices with
 extra precision elements (double or long double) to be used for storage only
 and with a minimum of operations defined. These would be used for accumulating
-the results of sums of squares and product matrices or multi-stage QR 
+the results of sums of squares and product matrices or multi-stage QR
 decompositions. 
    
 
    Allow matrix expressions
    
 
    I want to be able to write matrix expressions the way I would on paper. So
@@ -3479,9 +3479,9 @@ 
    Naming convention
    
 
    Exceptions to the general rule are the functions for transpose and inverse.
 To make matrix expressions more like the corresponding mathematical formulae, I
 have used the single letter abbreviations, t() and i(). 
    
-
    I am now switching to using lowercase for  functions with individual 
-words separated by "_". This is following the convention in 
-the standard library and I think it looks neater. Class names will remain with 
+
    I am now switching to using lowercase for  functions with individual
+words separated by "_". This is following the convention in
+the standard library and I think it looks neater. Class names will remain with
 individual words being capitalised.
    
 
    Row and column index ranges
    
 
    In mathematical work matrix subscripts usually start at one. In C, array
@@ -3521,70 +3521,70 @@ 
    5.4 Data storage
    
 
    next - skip -
 up - start 
    
 
    The stack and heap
    
-
    To understand how newmat stores matrices you need to know a little bit 
+
    To understand how newmat stores matrices you need to know a little bit
 about the heap and stack.
    
-
    The data values of variables or objects in a C++ program are stored in either 
-of two sections of memory called the stack and the heap. Sometimes 
+
    The data values of variables or objects in a C++ program are stored in either
+of two sections of memory called the stack and the heap. Sometimes
 there is more than one heap to cater for different sized variables.
    
 
    If you declare an automatic variable
    
 
       int x;
    
-
    then the value of x is stored on the stack. As you declare more 
-variables the stack gets bigger. When you exit a block (i.e a section of code 
-delimited by curly brackets {...}) the memory used by the automatic 
+
    then the value of x is stored on the stack. As you declare more
+variables the stack gets bigger. When you exit a block (i.e a section of code
+delimited by curly brackets {...}) the memory used by the automatic
 variables declared in the block is released and the stack shrinks.
    
 
    When you declare a variable with new, for example,
    
 
       int* y = new int;
    
-
    the pointer y is stored on the stack but the value it is 
-pointing to is stored on the heap. Memory on the heap is not 
+
    the pointer y is stored on the stack but the value it is
+pointing to is stored on the heap. Memory on the heap is not
 released until the program explicitly does this with a delete statement
    
 
       delete y;
    
 
    or the program exits.
    
-
    On the stack, variables and objects are is always added to the end of 
-the stack and are removed in the reverse order to that in which they are 
+
    On the stack, variables and objects are is always added to the end of
+the stack and are removed in the reverse order to that in which they are
 added - that is the last on will be the first off. This is not the case with the 
-heap, where the variables and objects can be removed in any order. So one 
-can get alternating pieces of used and unused memory. When a new variable or 
-object is declared on the heap the system needs to search for piece of 
-unused memory large enough to hold it. This means that storing on the heap 
-will usually be a slower process than storing on the stack. There is also 
-likely to be waste space on the heap because of gaps between the used 
+heap, where the variables and objects can be removed in any order. So one
+can get alternating pieces of used and unused memory. When a new variable or
+object is declared on the heap the system needs to search for piece of
+unused memory large enough to hold it. This means that storing on the heap
+will usually be a slower process than storing on the stack. There is also
+likely to be waste space on the heap because of gaps between the used
 blocks of memory that are too small for the next object you want to store on the
-heap. There is also the possibility of wasting space if you forget to 
-remove a variable or object on the heap even though you have finished 
-using it. However, the stack is usually limited to holding small objects 
-with size known at compile time. Large objects, objects whose size you don't 
-know at compile time, and objects that you want to persist after the end of the 
+heap. There is also the possibility of wasting space if you forget to
+remove a variable or object on the heap even though you have finished
+using it. However, the stack is usually limited to holding small objects
+with size known at compile time. Large objects, objects whose size you don't
+know at compile time, and objects that you want to persist after the end of the
 block need to be stored on the heap.
    
-
    In C++, the constructor/destructor system enables one to build 
-complicated objects such as matrices that behave as automatic variables stored 
-on the stack, so the programmer doesn't have to worry about deleting them 
-at the end of the block, but which really utilise the heap for storing 
+
    In C++, the constructor/destructor system enables one to build
+complicated objects such as matrices that behave as automatic variables stored
+on the stack, so the programmer doesn't have to worry about deleting them
+at the end of the block, but which really utilise the heap for storing
 their data.
    
 
    Structure of matrix objects
    
 
    Each matrix object contains the basic information such as the number of rows
 and columns, the amount of memory used, a status variable and a pointer to the data array which is on
 the heap. So if you declare a matrix
    
 
       Matrix A(1000,1000);
    
-
    there is an small amount of memory used on the stack for storing the numbers 
-of rows and columns, the amount of  memory used, the status variable and 
-the pointer together with 1,000,000 Real locations stored on the heap. 
+
    there is an small amount of memory used on the stack for storing the numbers
+of rows and columns, the amount of  memory used, the status variable and
+the pointer together with 1,000,000 Real locations stored on the heap.
 When you exit the block in which A is declared, the heap memory used by
-A is automatically returned to the system, as well as the memory used on 
+A is automatically returned to the system, as well as the memory used on
 the stack.
    
 
    Of course, if you use new to declare a matrix
    
 
       Matrix* B = new Matrix(1000,1000);
    
-
    both the information about the size and the actual data are stored on heap 
+
    both the information about the size and the actual data are stored on heap
 and not deleted until the program exits or you do an explicit delete:
    
 
       delete B;
    
-
    If you carry out an assignment with = or << or do a 
-resize() the data array currently associated with a matrix is destroyed and 
+
    If you carry out an assignment with = or << or do a
+resize() the data array currently associated with a matrix is destroyed and
 a new array generated. For example
    
 
       Matrix A(1000,1000);
    Matrix B(50, 50);
    ... put some values in B
    A = B;
    
-
    At the last step the heap memory associated with A is returned to the 
-system and a new block of heap memory is assigned to contain the new values. 
+
    At the last step the heap memory associated with A is returned to the
+system and a new block of heap memory is assigned to contain the new values.
 This happens even if there is no change in the amount of memory required. 
    
 
    One block or several
    
 
    The elements of the matrix are stored as a single array. Alternatives would
@@ -3766,9 +3766,9 @@ 
    5.6 Memory management -
 Could be improved?
 
 
-
    This is now all rather out of date. With Pentiums, at least, the important 
-requirement for speed seems to be to minimise transfers between the RAM memory 
-and the on-chip memory. There isn't much you can do about add and subtract, but 
+
    This is now all rather out of date. With Pentiums, at least, the important
+requirement for speed seems to be to minimise transfers between the RAM memory
+and the on-chip memory. There isn't much you can do about add and subtract, but
 there lots of possibilities for some of the other operations.
    
 
    5.7 Evaluation of expressions -
 lazy evaluation
    
@@ -3880,7 +3880,7 @@ 
    5.8 How to overcome an
 matrices of the same type then there is no need to access the individual rows
 or columns and a faster general algorithm is appropriate. 
    
 
    Generally the method works well. However symmetric matrices are not always
-handled very efficiently (yet) since complete rows are not stored explicitly. 
+handled very efficiently (yet) since complete rows are not stored explicitly.
 
    
 
    The original version of the package did not use this access by row or column
 method and provided the multitude of algorithms for the combination of
@@ -3888,7 +3888,7 @@ 
    5.8 How to overcome an
 longer than the present one when providing the same facilities with 5 distinct
 types of matrices. It would have been very difficult to increase the number of
 matrix types in the original version. Apparently 4 to 5 types is about the
-break even point for switching to the approach adopted in the present package. 
+break even point for switching to the approach adopted in the present package.
 
    
 
    However it must also be admitted that there is a substantial overhead in the
 approach adopted in the present package for small matrices. The test program
@@ -3905,10 +3905,10 @@ 
    5.9 Destruction of
 
    Versions before version 5 of newmat did not work correctly with Gnu C++
 (version 5 or earlier). This was because the tree structure used to represent a
 matrix expression was set up on the stack.  Early versions of Gnu C++ destroyed temporary structures as soon as the
-function that accesses them finished.  To overcome this problem, there was an 
-option to store the temporaries forming the tree structure on the heap (created 
-with new) and to delete them explicitly. Now that the C++ standards committee 
-has said that temporary structures should not be destroyed before a statement 
+function that accesses them finished.  To overcome this problem, there was an
+option to store the temporaries forming the tree structure on the heap (created
+with new) and to delete them explicitly. Now that the C++ standards committee
+has said that temporary structures should not be destroyed before a statement
 finishes, I have deleted this option.
    
 
    5.10 A calculus of matrix
 types
    
@@ -4016,9 +4016,9 @@ 
    5.11 Pointer arithmetic
    
    }
    -

    which is more complicated and consequently introduces more chance of error. +

    which is more complicated and consequently introduces more chance of error.

    -

    5.12 Error handling +

    5.12 Error handling

    next - skip - up - start

    @@ -4027,7 +4027,7 @@

    5.12 Error handling newmat08 was released (in 1995), compiler exception handling in the compilers I had access to was unreliable. I recommended you used my simulated exceptions. In 1997 compiler supported exceptions seemed to work on a variety of -compilers - but not all compilers. This was still true in 2001. One compiler +compilers - but not all compilers. This was still true in 2001. One compiler company was still having problems in 2003 (not sure about 2004). Try using the compiler supported exceptions if you have a recent compiler, but if you are getting strange crashes or errors try going back to my simulated exceptions.

    @@ -4076,7 +4076,7 @@

    5.14 Complex matrices

    You can simulate most complex operations by representing Z = X + iY by 

        /  X   Y \
-    \ -Y   X / 
+    \ -Y   X /

    Most matrix operations will simulate the corresponding complex operation, @@ -4094,7 +4094,7 @@

    6. Function summary

    up - start

    -

    6.1 Member functions for matrices +

    6.1 Member functions for matrices and matrix expressions
    6.2 Member functions for matrices
    6.3 Operators
    @@ -4104,15 +4104,15 @@

    6. Function summary

    This section lists member and global functions for matrices defined in -newmat.h. Where there are alternative names the lower-case non-capitalised +newmat.h. Where there are alternative names the lower-case non-capitalised versions are the preferred ones.

    6.1 Member functions for matrices and matrix expressions

    next - skip - up - start

    -

    Member functions for matrices and matrix expressions. These do not apply to -CroutMatrix and BandLUMatrix +

    Member functions for matrices and matrix expressions. These do not apply to +CroutMatrix and BandLUMatrix unless explicitly noted.

    @@ -4211,8 +4211,8 @@

    6.1 Member functions for matrices and matri select a range of columns of a matrix - Scalar functions - maxima & - minima

    (also global versions of maximum(), minimum(),  + Scalar functions - maxima & + minima

    (also global versions of maximum(), minimum(),  maximum_absolute_value(), mimimum_absolute_value()) .maximum_absolute_value()
    .    MaximumAbsoluteValue() @@ -4331,7 +4331,7 @@

    6.1 Member functions for matrices and matri

     

    6.2 Member functions for matrices

    next - skip - up - start

    -

    Member functions for matrices but not matrix expressions. These do not apply +

    Member functions for matrices but not matrix expressions. These do not apply to CroutMatrix and BandLUMatrix unless explicitly noted.

    @@ -4340,7 +4340,7 @@

    6.2 Member functions for matrices

    - @@ -4357,11 +4357,11 @@

    6.2 Member functions for matrices

    - - @@ -4380,7 +4380,7 @@

    6.2 Member functions for matrices

    - @@ -4402,7 +4402,7 @@

    6.2 Member functions for matrices

    - @@ -4414,12 +4414,12 @@

    6.2 Member functions for matrices

    - - @@ -4452,7 +4452,7 @@

    6.2 Member functions for matrices

    - @@ -4472,7 +4472,7 @@

    6.2 Member functions for matrices

    - @@ -4484,7 +4484,7 @@

    6.2 Member functions for matrices

    6.3 Operators

    next - skip - up - start

    -

    Operators for matrices and matrix expressions. These do not apply to CroutMatrix and BandLUMatrix +

    Operators for matrices and matrix expressions. These do not apply to CroutMatrix and BandLUMatrix unless explicitly noted.

    @@ -4567,7 +4567,7 @@

    6.3 Operators

    - @@ -4599,7 +4599,7 @@

    6.4 Global functions - newmat.h

    - @@ -4675,7 +4675,7 @@

    6.5 Global functions - newmatap.h

    - @@ -4686,19 +4686,19 @@

    6.5 Global functions - newmatap.h

    - - - @@ -4707,7 +4707,7 @@

    6.5 Global functions - newmatap.h

    - @@ -4717,39 +4717,39 @@

    6.5 Global functions - newmatap.h

    - - - - - - - @@ -4758,57 +4758,57 @@

    6.5 Global functions - newmatap.h

    - - - - - - - - - - - - - - - @@ -4823,45 +4823,45 @@

    6.5 Global functions - newmatap.h

    - - - - - - - - - @@ -4870,7 +4870,7 @@

    6.5 Global functions - newmatap.h

    - @@ -4879,7 +4879,7 @@

    6.5 Global functions - newmatap.h

    - @@ -4888,7 +4888,7 @@

    6.5 Global functions - newmatap.h

    - @@ -4910,12 +4910,12 @@

    6.5 Global functions - newmatap.h

    - -
    description
    Element access

    (see +

    		Element access

    (see also operators)

    		 .element(int,int) access element - subscripts start at 0
    .swap(Matrix&)swap bodies of two matrices of same type (also global + swap bodies of two matrices of same type (also global version, also works with CroutMatrix and BandLUMatrix)
    Scalar functions - size and + Scalar functions - size and shape

    (also work with CroutMatrix and BandLUMatrix)

    		 .type()
    .Type()
    .size()
    .Storage()    		number of stored elements (including unused elements in band + number of stored elements (including unused elements in band matrices)
    .const_data_indx()constant pointer to row swap array (CroutMatrix and BandLUMatrix + constant pointer to row swap array (CroutMatrix and BandLUMatrix only)
    .is_singular()
    .IsSingular()    		test for exact singularity (CroutMatrix and BandLUMatrix + test for exact singularity (CroutMatrix and BandLUMatrix only)
    .even_exchanges()true if there have been an even number of row exchanges (CroutMatrix and BandLUMatrix + true if there have been an even number of row exchanges (CroutMatrix and BandLUMatrix only)
    .resize(int,int)
    .ReSize(int,int)    		change the dimensions (non-square matrices, triangular band + change the dimensions (non-square matrices, triangular band matrices and symmetric band matrices)
    .resize_keep(int)change the dimensions, keep values (vectors and square matrices, + change the dimensions, keep values (vectors and square matrices, not band)
    =copy Real + copy Real to all elements
    description
    Binary operators

    (see +

    		Binary operators

    (see also operators)

    		 SP(const BaseMatrix&, const BaseMatrix&) element-wise product of two matrices
    updateQRZT(Matrix&, LowerTriangularMatrix&)
    UpdateQRZT(Matrix&, LowerTriangularMatrix&)    		add extra columns to transposed version of QRZ + add extra columns to transposed version of QRZ transform
    updateQRZ(const Matrix&, Matrix&, Matrix&)
    UpdateQRZ(const Matrix&, Matrix&, Matrix&)    		combine the results of the solve parts of QRZ + combine the results of the solve parts of QRZ transforms on two blocks of data
    updateQRZ(UpperTriangularMatrix&, UpperTriangularMatrix&)
    UpdateQRZ(UpperTriangularMatrix&, UpperTriangularMatrix&)    		combine the results of QRZ transforms on two + combine the results of QRZ transforms on two blocks of data
    updateQRZ(const UpperTriangularMatrix&, Matrix&, Matrix&)
    UpdateQRZ(const UpperTriangularMatrix&, Matrix&, Matrix&)    		combine the results of the solve parts of QRZ + combine the results of the solve parts of QRZ transforms on two blocks of data
    extend a set of orthonormal columns
    Cholesky + Cholesky decomposition Cholesky(const SymmetricMatrix&) Cholesky decomposition of symmetric matrix Cholesky decomposition of symmetric band matrix
    Update Cholesky + Update Cholesky decompositionupdate_Cholesky (UpperTriangularMatrix&, + update_Cholesky (UpperTriangularMatrix&, RowVector)
    - UpdateCholesky (UpperTriangularMatrix&, + UpdateCholesky (UpperTriangularMatrix&, RowVector)    		 add extra row to Cholesky/QR decomposition
    downdate_Cholesky (UpperTriangularMatrix&, + downdate_Cholesky (UpperTriangularMatrix&, RowVector)
    - DowndateCholesky (UpperTriangularMatrix&, + DowndateCholesky (UpperTriangularMatrix&, RowVector)    		 remove row from Cholesky/QR decomposition
    right_circular_update_Cholesky (UpperTriangularMatrix&, + right_circular_update_Cholesky (UpperTriangularMatrix&, int, int)
    - RightCircularUpdateCholesky (UpperTriangularMatrix&, + RightCircularUpdateCholesky (UpperTriangularMatrix&, int, int)    		 rearrange columns in Cholesky/QR decomposition
    left_circular_update_Cholesky (UpperTriangularMatrix&, + left_circular_update_Cholesky (UpperTriangularMatrix&, int, int)
    - LeftCircularUpdateCholesky (UpperTriangularMatrix&, + LeftCircularUpdateCholesky (UpperTriangularMatrix&, int, int)    		 rearrange columns in Cholesky/QR decomposition
    Singular value + Singular value decompositionSVD(const Matrix&, DiagonalMatrix&, + SVD(const Matrix&, DiagonalMatrix&, Matrix&, Matrix&, bool, bool) singular value decomposition - get U and V
    SVD decomposition - get just singular values
    SVD(const Matrix& A, DiagonalMatrix& D, + SVD(const Matrix& A, DiagonalMatrix& D, Matrix&, bool) SVD decomposition - get U
    Eigenvalue + Eigenvalue decomposition of a symmetric matrixJacobi(const SymmetricMatrix&, + Jacobi(const SymmetricMatrix&, DiagonalMatrix&) Jacobi eigenvalue decomposition - get only eigenvalues
    Jacobi(const SymmetricMatrix&, + Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&)Jacobi eigenvalue decomposition - get only + Jacobi eigenvalue decomposition - get only eigenvalues
    Jacobi(const SymmetricMatrix&, + Jacobi(const SymmetricMatrix&, DiagonalMatrix&, Matrix&)Jacobi eigenvalue decomposition - get + Jacobi eigenvalue decomposition - get eigenvalues and eigenvectors
    Jacobi(const SymmetricMatrix&, + Jacobi(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&, Matrix&, bool)Jacobi eigenvalue decomposition - get + Jacobi eigenvalue decomposition - get eigenvalues and eigenvectors
    eigenvalues(const SymmetricMatrix&, + eigenvalues(const SymmetricMatrix&, DiagonalMatrix&)
    - EigenValues(const SymmetricMatrix&, + EigenValues(const SymmetricMatrix&, DiagonalMatrix&)    		Householder eigenvalue decomposition - get only + Householder eigenvalue decomposition - get only eigenvalues
    eigenvalues(const SymmetricMatrix&, + eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&)
    - EigenValues(const SymmetricMatrix&, + EigenValues(const SymmetricMatrix&, DiagonalMatrix&, SymmetricMatrix&)    		Householder eigenvalue decomposition with back + Householder eigenvalue decomposition with back transform - get only eigenvalues
    eigenvalues(const SymmetricMatrix&, + eigenvalues(const SymmetricMatrix&, DiagonalMatrix&, Matrix&)
    - EigenValues(const SymmetricMatrix&, + EigenValues(const SymmetricMatrix&, DiagonalMatrix&, Matrix&)    		Householder eigenvalue decomposition - get + Householder eigenvalue decomposition - get eigenvalues and eigenvectors
    descending sort
    Fast Fourier + Fast Fourier transformFFT(const ColumnVector&, const + FFT(const ColumnVector&, const ColumnVector&, ColumnVector&, ColumnVector&) fast Fourier transform
    FFTI(const ColumnVector&, const + FFTI(const ColumnVector&, const ColumnVector&, ColumnVector&, ColumnVector&) fast Fourier transform - inverse
    RealFFT(const ColumnVector&, ColumnVector&, + RealFFT(const ColumnVector&, ColumnVector&, ColumnVector&) FFT of real vector
    RealFFTI(const ColumnVector&, const + RealFFTI(const ColumnVector&, const ColumnVector&, ColumnVector&) FFT of real vector - inverse
    FFT2(const Matrix& U, const Matrix& V, + FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y) two dimensional FFT
    FFT2I(const Matrix& U, const Matrix& V, + FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y) two dimensional FFT - inverse
    Fast + Fast trigonometric transform DCT_II(const ColumnVector&, ColumnVector&) type II discrete cosine transform
    DCT_II_inverse(const ColumnVector&, + DCT_II_inverse(const ColumnVector&, ColumnVector&) type II discrete cosine transform - inverse
    type II discrete sine transform
    DST_II_inverse(const ColumnVector&, + DST_II_inverse(const ColumnVector&, ColumnVector&) type II discrete sine transform - inverse
    discrete cosine transform
    DCT_inverse(const ColumnVector&, + DCT_inverse(const ColumnVector&, ColumnVector&) discrete cosine transform - inverse
    discrete sine transform
    DST_inverse(const ColumnVector&, + DST_inverse(const ColumnVector&, ColumnVector&) discrete sine transform - inverse
    Helmert_transpose(const ColumnVector&, bool=false)
    Helmert_transpose(const Matrix&, bool=false)    		multiply by transpose of Helmert transform + multiply by transpose of Helmert transform matrix
    Helmert_transpose(const ColumnVector&, int, bool=false)multiply by transpose of Helmert transform + multiply by transpose of Helmert transform matrix, return one element of result
    @@ -5030,24 +5030,24 @@

    7. Change history

    next - skip - up - start

    Newmat11 - November, 2008:

    -

    Remove work-arounds for older compilers, Borland Builder 6 and -Open Watcom compatibility, SquareMatrix, -load from array of ints, crossproducts, Cholesky and QRZ update functions, swap -functions, FFT2, access to arrays in traditional C functions, SimpleIntArray -class, compatibility with Numerical Recipes in C++, sum_rows(), sum_columns(), sum_squares_rows() and -sum_squares_columns() functions, extend_orthogonal function, resize_keep -function, speed-ups and -bug-fixes, change to lower case for functions, can copy CroutMatrix, -BandLUMatrix, Helmert transform, compatibility fix for G++ 4.1, start inserting +

    Remove work-arounds for older compilers, Borland Builder 6 and +Open Watcom compatibility, SquareMatrix, +load from array of ints, crossproducts, Cholesky and QRZ update functions, swap +functions, FFT2, access to arrays in traditional C functions, SimpleIntArray +class, compatibility with Numerical Recipes in C++, sum_rows(), sum_columns(), sum_squares_rows() and +sum_squares_columns() functions, extend_orthogonal function, resize_keep +function, speed-ups and +bug-fixes, change to lower case for functions, can copy CroutMatrix, +BandLUMatrix, Helmert transform, compatibility fix for G++ 4.1, start inserting comments for Doxygen, SP_eq, scientific format for output, fix for 64 bits.

    Newmat10A - October, 2002, Newmat10B - January 2005:

    -

    Fix error in Kronecker product; fixes for Intel and GCC3 +

    Fix error in Kronecker product; fixes for Intel and GCC3 compilers.

    Newmat10 - January, 2002:

    Improve compatibility with GCC, fix errors in FFT and GenericMatrix, update simulated exceptions, maxima, minima, determinant, dot product and Frobenius norm functions, update make files for CC and GCC, faster -FFT, A.ReSize(B), fix pointer arithmetic, << for loading +FFT, A.ReSize(B), fix pointer arithmetic, << for loading data into rows, IdentityMatrix, Kronecker product, sort singular values.

    Newmat09 - September, 1997:

    Operator ==, !=, +=, -=, @@ -5087,7 +5087,7 @@

    7. Change history

    Newmat05 - June 1992:

    For private release only

    Newmat04 - December 1991:

    -

    Fix problem with G++1.40, some extra documentation +

    Fix problem with G++1.40, some extra documentation

    Newmat03 - November 1991:

    Col and Cols become Column and Columns. Added Sort, SVD, @@ -5142,4 +5142,4 @@

    8. Problem report form

     

     

    - \ No newline at end of file + diff --git a/src/newmat/nm_b55.mak b/src/newmat/nm_b55.mak index c7c45030..98911c2d 100644 --- a/src/newmat/nm_b55.mak +++ b/src/newmat/nm_b55.mak @@ -14,7 +14,7 @@ CC = bcc32 -W- -v- -H- -3 -N -Og -Oi -Ov -f -I$(INCLUDEPATH) .cpp.obj: $(CC) -c {$< } -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -231,4 +231,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_b56.mak b/src/newmat/nm_b56.mak index 8414f19b..0a372d3d 100644 --- a/src/newmat/nm_b56.mak +++ b/src/newmat/nm_b56.mak @@ -14,7 +14,7 @@ CC = bcc32 -W- -v- -H- -6 -Og -Oi -Ov -I$(INCLUDEPATH) .cpp.obj: $(CC) -c {$< } -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -231,4 +231,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_cc.mak b/src/newmat/nm_cc.mak index cb3a7354..84e78004 100644 --- a/src/newmat/nm_cc.mak +++ b/src/newmat/nm_cc.mak @@ -10,11 +10,11 @@ PRE = ./ .cpp.o: rm -f $*.cxx - ln $*.cpp $*.cxx + ln $*.cpp $*.cxx $(CXX) $(CXXFLAGS) -c $*.cxx - rm $*.cxx + rm $*.cxx -everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch +everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch newmat_lobj = newmat1.o newmat2.o newmat3.o newmat4.o newmat5.o newmat6.o newmat7.o newmat8.o newmatex.o bandmat.o submat.o myexcept.o cholesky.o evalue.o fft.o hholder.o jacobi.o newfft.o sort.o svd.o nm_misc.o newmatrm.o newmat9.o @@ -214,4 +214,3 @@ sl_ex.txx: sl_ex garch.txx: garch $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_ex1.cpp b/src/newmat/nm_ex1.cpp index 05bb45c3..28b7d666 100644 --- a/src/newmat/nm_ex1.cpp +++ b/src/newmat/nm_ex1.cpp @@ -11,7 +11,7 @@ // should not be required for this example // included because it seems to help MS VC6 // when you have namespace turned on - + #include "newmatio.h" // newmat headers including output functions #ifdef use_namespace @@ -25,7 +25,7 @@ int my_main() // called by main() // declare a matrix Matrix X(4,4); - + // load values row by row X.row(1) << 3.7 << -2.1 << 7.4 << -1.0; X.row(2) << 4.1 << 0.0 << 3.9 << 4.0; @@ -35,16 +35,16 @@ int my_main() // called by main() // print the matrix cout << "Matrix X" << endl; cout << setw(15) << setprecision(8) << X << endl; - + // calculate its inverse and print it Matrix Y = X.i(); cout << "Inverse of X" << endl; cout << setw(15) << setprecision(8) << Y << endl; - + // multiply X by its inverse and print the result (should be near identity) cout << "X * inverse of X" << endl; cout << setw(15) << setprecision(8) << (X * Y) << endl; - + return 0; } @@ -60,4 +60,3 @@ int main() } ///@} - diff --git a/src/newmat/nm_ex1.txt b/src/newmat/nm_ex1.txt index 53eab140..9da0390b 100644 --- a/src/newmat/nm_ex1.txt +++ b/src/newmat/nm_ex1.txt @@ -15,4 +15,3 @@ X * inverse of X 0.00000000 1.00000000 0.00000000 0.00000000 -0.00000000 -0.00000000 1.00000000 -0.00000000 -0.00000000 0.00000000 0.00000000 1.00000000 - diff --git a/src/newmat/nm_ex2.cpp b/src/newmat/nm_ex2.cpp index d0887a8a..339e8225 100644 --- a/src/newmat/nm_ex2.cpp +++ b/src/newmat/nm_ex2.cpp @@ -12,7 +12,7 @@ /// /// The dimensions of this matrix are not large enough for there to be numerical /// problems but we will be able to see that wide range of values of the -/// eigenvalues. +/// eigenvalues. #define WANT_STREAM // include iostream and iomanipulators @@ -27,13 +27,13 @@ using namespace RBD_LIBRARIES; int my_main() // called by main() { Tracer tr("my_main "); // for tracking exceptions - + int n = 7; // this is the order we will work with int i, j; // declare a matrix SymmetricMatrix H(n); - + // load values for Hilbert matrix for (i = 1; i <= n; ++i) for (j = 1; j <= i; ++j) H(i, j) = 1.0 / (i + j - 1); @@ -41,7 +41,7 @@ int my_main() // called by main() // print the matrix cout << "SymmetricMatrix H" << endl; cout << setw(10) << setprecision(7) << H << endl; - + // calculate its eigenvalues and eigenvectors and print them Matrix U; DiagonalMatrix D; eigenvalues(H, D, U); @@ -51,13 +51,13 @@ int my_main() // called by main() cout << setw(10) << setprecision(7) << U << endl; // check orthogonality - cout << "U * U.t() (should be near identity)" << endl; + cout << "U * U.t() (should be near identity)" << endl; cout << setw(10) << setprecision(7) << (U * U.t()) << endl; - + // check decomposition - cout << "U * D * U.t() (should be near H)" << endl; + cout << "U * D * U.t() (should be near H)" << endl; cout << setw(10) << setprecision(7) << (U * D * U.t()) << endl; - + return 0; } diff --git a/src/newmat/nm_ex2.txt b/src/newmat/nm_ex2.txt index fb988220..f552b8a9 100644 --- a/src/newmat/nm_ex2.txt +++ b/src/newmat/nm_ex2.txt @@ -42,4 +42,3 @@ U * D * U.t() (should be near H) 0.2000000 0.1666667 0.1428571 0.1250000 0.1111111 0.1000000 0.0909091 0.1666667 0.1428571 0.1250000 0.1111111 0.1000000 0.0909091 0.0833333 0.1428571 0.1250000 0.1111111 0.1000000 0.0909091 0.0833333 0.0769231 - diff --git a/src/newmat/nm_ex3.cpp b/src/newmat/nm_ex3.cpp index 8c3309fd..a5fb9ed0 100644 --- a/src/newmat/nm_ex3.cpp +++ b/src/newmat/nm_ex3.cpp @@ -28,24 +28,24 @@ int my_main() // called by main() cout << "Minimum " << FloatingPointPrecision::Minimum() << endl; cout << "LnMinimum " << FloatingPointPrecision::LnMinimum() << endl; cout << endl; - + ColumnVector values(5), sel_values(4); - + values(1) = sel_values(1) = FloatingPointPrecision::Epsilon(); values(2) = FloatingPointPrecision::Maximum(); values(3) = sel_values(2) = FloatingPointPrecision::LnMaximum(); values(4) = sel_values(3) = FloatingPointPrecision::Minimum(); values(5) = sel_values(4) = FloatingPointPrecision::LnMinimum(); - + cout << "Load values into a ColumnVector and print in scientific format" << endl << endl; cout << "Values in scientific format" << endl; cout << setw(15) << setprecision(5) << scientific << values << endl; - + cout << "Print in fixed format" << endl << endl; cout << "Selected values in fixed format (omit Maximum)" << endl; cout << setw(15) << setprecision(5) << fixed << sel_values << endl; - + cout << "Test Epsilon" << endl << endl; Real x = 1.0 + FloatingPointPrecision::Epsilon(); cout << "This value should be non-zero" << endl; @@ -56,7 +56,7 @@ int my_main() // called by main() cout << "(1.0 + Epsilon / 2) - 1.0" << setw(15) << setprecision(5) << scientific << (Real)(x - 1.0) << endl; cout << endl; - + cout << "Test Minimum" << endl << endl; cout << "This value should be non-zero" << endl; cout << "Minimum " << setw(15) << setprecision(5) << @@ -70,11 +70,11 @@ int my_main() // called by main() cout << "Minimum * Epsilon / 2.0 " << setw(15) << setprecision(5) << scientific << x << endl; cout << endl; - + cout << "Test LnMaximum and LnMinimum" << endl << endl; x = exp(FloatingPointPrecision::LnMaximum()); cout << "This should print a value close to Maximum" << endl; - cout << "Fails with some compilers - needs fixing" << endl; + cout << "Fails with some compilers - needs fixing" << endl; cout << "exp(LnMaximum) " << setw(15) << setprecision(5) << scientific << x << endl; x = exp(FloatingPointPrecision::LnMinimum()); @@ -82,9 +82,9 @@ int my_main() // called by main() cout << "exp(LnMinimum) " << setw(15) << setprecision(5) << scientific << x << endl; cout << endl; - - - + + + return 0; } @@ -100,4 +100,3 @@ int main() } ///@} - diff --git a/src/newmat/nm_ex3.txt b/src/newmat/nm_ex3.txt index 6bb79ff4..de926135 100644 --- a/src/newmat/nm_ex3.txt +++ b/src/newmat/nm_ex3.txt @@ -47,4 +47,3 @@ Fails with some compilers - needs fixing exp(LnMaximum) 1.79769e+308 This should print a value close to Minimum, not zero exp(LnMinimum) 2.22507e-308 - diff --git a/src/newmat/nm_gnu.mak b/src/newmat/nm_gnu.mak index aca0335a..7441f0b2 100644 --- a/src/newmat/nm_gnu.mak +++ b/src/newmat/nm_gnu.mak @@ -9,7 +9,7 @@ MINOR = 0 %.o: %.cpp $(CXX) $(CXXFLAGS) -c $*.cpp -everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch +everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch newmat_lobj = newmat1.o newmat2.o newmat3.o newmat4.o newmat5.o newmat6.o newmat7.o newmat8.o newmatex.o bandmat.o submat.o myexcept.o cholesky.o evalue.o fft.o hholder.o jacobi.o newfft.o sort.o svd.o nm_misc.o newmatrm.o newmat9.o @@ -209,4 +209,3 @@ sl_ex.txx: sl_ex garch.txx: garch $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_i10.mak b/src/newmat/nm_i10.mak index 6fcf660f..08b476ee 100644 --- a/src/newmat/nm_i10.mak +++ b/src/newmat/nm_i10.mak @@ -1,5 +1,5 @@ -conlibs = +conlibs = DIFF = sdiff PRE = @@ -9,7 +9,7 @@ PRE = .cpp.obj: icl -c -GX -GR -Ge -GS -Qprec -Qprec_div -nologo -Qlong_double $*.cpp -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -208,4 +208,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_i8.mak b/src/newmat/nm_i8.mak index 30b36ba8..f8a06024 100644 --- a/src/newmat/nm_i8.mak +++ b/src/newmat/nm_i8.mak @@ -9,7 +9,7 @@ PRE = .cpp.obj: icl -c -GX -GR -Ge -GS -Qprec -Qprec_div -nologo -Qlong_double $*.cpp -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -208,4 +208,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_il8.mak b/src/newmat/nm_il8.mak index 4d174121..c5f49d5f 100644 --- a/src/newmat/nm_il8.mak +++ b/src/newmat/nm_il8.mak @@ -9,7 +9,7 @@ MINOR = 0 %.o: %.cpp $(CXX) $(CXXFLAGS) -c $*.cpp -everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch +everything: tmt example nm_ex1 nm_ex2 nm_ex3 test_exc nl_ex sl_ex garch newmat_lobj = newmat1.o newmat2.o newmat3.o newmat4.o newmat5.o newmat6.o newmat7.o newmat8.o newmatex.o bandmat.o submat.o myexcept.o cholesky.o evalue.o fft.o hholder.o jacobi.o newfft.o sort.o svd.o nm_misc.o newmatrm.o newmat9.o @@ -209,4 +209,3 @@ sl_ex.txx: sl_ex garch.txx: garch $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_m6.mak b/src/newmat/nm_m6.mak index 97931f36..036d4322 100644 --- a/src/newmat/nm_m6.mak +++ b/src/newmat/nm_m6.mak @@ -9,7 +9,7 @@ PRE = .cpp.obj: cl -c -W3 -Ox -GX $*.cpp -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -208,4 +208,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_m8.mak b/src/newmat/nm_m8.mak index 5b29da94..a56efa3d 100644 --- a/src/newmat/nm_m8.mak +++ b/src/newmat/nm_m8.mak @@ -9,7 +9,7 @@ PRE = .cpp.obj: cl -c -W3 -Ox -EHsc $*.cpp -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -208,4 +208,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/nm_misc.cpp b/src/newmat/nm_misc.cpp index 611aaad2..10116d2c 100644 --- a/src/newmat/nm_misc.cpp +++ b/src/newmat/nm_misc.cpp @@ -27,7 +27,7 @@ ReturnMatrix Helmert(int n, bool full) Tracer et("Helmert "); if (n <= 0) Throw(ProgramException("Dimension <= 0 ")); Matrix H; - + if (full) H.resize(n,n); else H.resize(n-1,n); H = 0.0; for (int i = 1; i < n; ++i) @@ -37,7 +37,7 @@ ReturnMatrix Helmert(int n, bool full) } if (full) { H.row(n) = 1.0 / sqrt((Real)n); } H.release(); return H.for_return(); -} +} @@ -55,7 +55,7 @@ ReturnMatrix Helmert(const ColumnVector& X, bool full) { sum += X(i); Y(i) = (i * X(i+1) - sum) / sqrt((Real)i * (i+1)); } if (full) { sum += X(n); Y(n) = sum / sqrt((Real)n); } Y.release(); return Y.for_return(); -} +} // same as above for X a ColumnVector, length n, element j = 1; otherwise 0 ReturnMatrix Helmert(int n, int j, bool full) @@ -71,7 +71,7 @@ ReturnMatrix Helmert(int n, int j, bool full) for (int i = j; i < n; ++i) Y(i) = - 1.0 / sqrt((Real)i * (i+1)); if (full) Y(n) = 1.0 / sqrt((Real)n); Y.release(); return Y.for_return(); -} +} ReturnMatrix Helmert_transpose(const ColumnVector& Y, bool full) { diff --git a/src/newmat/nm_ow.mak b/src/newmat/nm_ow.mak index 0e69d234..6c3af9af 100644 --- a/src/newmat/nm_ow.mak +++ b/src/newmat/nm_ow.mak @@ -8,7 +8,7 @@ PRE = -everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe +everything: tmt.exe example.exe nm_ex1.exe nm_ex2.exe nm_ex3.exe test_exc.exe nl_ex.exe sl_ex.exe garch.exe newmat_lobj = newmat1.obj newmat2.obj newmat3.obj newmat4.obj newmat5.obj newmat6.obj newmat7.obj newmat8.obj newmatex.obj bandmat.obj submat.obj myexcept.obj cholesky.obj evalue.obj fft.obj hholder.obj jacobi.obj newfft.obj sort.obj svd.obj nm_misc.obj newmatrm.obj newmat9.obj @@ -209,4 +209,3 @@ sl_ex.txx: sl_ex.exe garch.txx: garch.exe $(PRE)garch > garch.txx $(DIFF) garch.txt garch.txx - diff --git a/src/newmat/precisio.h b/src/newmat/precisio.h index 141e1555..e22a1e3a 100644 --- a/src/newmat/precisio.h +++ b/src/newmat/precisio.h @@ -25,7 +25,7 @@ namespace NEWMAT { #endif using namespace std; - + /// Floating point precision. class FloatingPointPrecision { @@ -52,7 +52,7 @@ class FloatingPointPrecision { return (Real)log(Maximum()); } static Real Minimum() // minimum positive value - { return numeric_limits::min(); } + { return numeric_limits::min(); } static int MinimumDecimalExponent() // minimum decimal exponent { return numeric_limits::min_exponent10; } diff --git a/src/newmat/rbd.css b/src/newmat/rbd.css index d512592b..7db67ff6 100644 --- a/src/newmat/rbd.css +++ b/src/newmat/rbd.css @@ -17,4 +17,4 @@ a:link { color: blue } a:visited { color: green } a:active { color: red } table { font-size: 10pt } -.question { color: navy; } \ No newline at end of file +.question { color: navy; } diff --git a/src/newmat/readme.txt b/src/newmat/readme.txt index 0293c78f..f6aad1b2 100644 --- a/src/newmat/readme.txt +++ b/src/newmat/readme.txt @@ -8,5 +8,3 @@ Documentation is in nm11.htm. This library is freeware. See the documentation for details. To contact the author please email robert@statsresearch.co.nz - - diff --git a/src/newmat/sl_ex.cpp b/src/newmat/sl_ex.cpp index 485b6e29..599f3a32 100644 --- a/src/newmat/sl_ex.cpp +++ b/src/newmat/sl_ex.cpp @@ -3,7 +3,7 @@ /// \file sl_ex.cpp /// Test one-dimensional solve function. -/// This is an example of the use of solution to find the cube root of +/// This is an example of the use of solution to find the cube root of /// the integers -10 to 10 /// /// you will need to compile and link solution.cpp and my_except.cpp @@ -50,7 +50,7 @@ int main() } CatchAll { - cout << "\nProgram fails - exception generated\n\n"; + cout << "\nProgram fails - exception generated\n\n"; } return 0; } diff --git a/src/newmat/solution.h b/src/newmat/solution.h index b752e18a..cfcbb61c 100644 --- a/src/newmat/solution.h +++ b/src/newmat/solution.h @@ -96,4 +96,3 @@ class OneDimSolve ///@} - diff --git a/src/newmat/submat.cpp b/src/newmat/submat.cpp index f3cc95d8..2074f5c4 100644 --- a/src/newmat/submat.cpp +++ b/src/newmat/submat.cpp @@ -256,7 +256,7 @@ void GetSubMatrix::operator+=(const BaseMatrix& bmx) SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate(); if (row_number != gmx->Nrows() || col_number != gmx->Ncols()) Throw(IncompatibleDimensionsException()); - if (gm->type().is_symmetric() && + if (gm->type().is_symmetric() && ( ! gmx->type().is_symmetric() || row_skip != col_skip) ) Throw(ProgramException("Illegal operation on symmetric")); MatrixRow mrx(gmx, LoadOnEntry); @@ -289,7 +289,7 @@ void GetSubMatrix::SP_eq(const BaseMatrix& bmx) SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate(); if (row_number != gmx->Nrows() || col_number != gmx->Ncols()) Throw(IncompatibleDimensionsException()); - if (gm->type().is_symmetric() && + if (gm->type().is_symmetric() && ( ! gmx->type().is_symmetric() || row_skip != col_skip) ) Throw(ProgramException("Illegal operation on symmetric")); MatrixRow mrx(gmx, LoadOnEntry); @@ -321,7 +321,7 @@ void GetSubMatrix::operator-=(const BaseMatrix& bmx) SetUpLHS(); gmx = ((BaseMatrix&)bmx).Evaluate(); if (row_number != gmx->Nrows() || col_number != gmx->Ncols()) Throw(IncompatibleDimensionsException()); - if (gm->type().is_symmetric() && + if (gm->type().is_symmetric() && ( ! gmx->type().is_symmetric() || row_skip != col_skip) ) Throw(ProgramException("Illegal operation on symmetric")); MatrixRow mrx(gmx, LoadOnEntry); @@ -398,4 +398,3 @@ void GetSubMatrix::operator*=(Real r) #endif ///@} - diff --git a/src/newmat/test_exc.txt b/src/newmat/test_exc.txt index a0be2118..69a9c5df 100644 --- a/src/newmat/test_exc.txt +++ b/src/newmat/test_exc.txt @@ -207,4 +207,3 @@ compilers - see documentation) Checking for lost memory (large block): 9191360 9191360 - ok Checking for lost memory (small block): 9191360 9191360 - ok - diff --git a/src/newmat/tmt.cpp b/src/newmat/tmt.cpp index 3b4af857..db25ce92 100644 --- a/src/newmat/tmt.cpp +++ b/src/newmat/tmt.cpp @@ -194,21 +194,21 @@ void FillWithValues(MultWithCarry&MWC, Matrix& M) for (int i = 1; i <= nr; ++i) for (int j = 1; j <= nc; ++ j) M(i, j) = MWC.Next(); } - + void FillWithValues(MultWithCarry&MWC, UpperTriangularMatrix& M) { int nr = M.nrows(); for (int i = 1; i <= nr; ++i) for (int j = i; j <= nr; ++ j) M(i, j) = MWC.Next(); } - + void FillWithValues(MultWithCarry&MWC, LowerTriangularMatrix& M) { int nr = M.nrows(); for (int i = 1; i <= nr; ++i) for (int j = 1; j <= i; ++ j) M(i, j) = MWC.Next(); } - + void FillWithValues(MultWithCarry&MWC, DiagonalMatrix& M) { int nr = M.nrows(); @@ -222,8 +222,8 @@ void FillWithValues(MultWithCarry&MWC, SymmetricMatrix& M) for (int i = 1; i <= nr; ++i) for (int j = 1; j <= i; ++ j) M(i, j) = MWC.Next(); } - - + + #ifdef use_namespace } @@ -271,8 +271,8 @@ int main() TestTypeAdd(); TestTypeMult(); TestTypeConcat(); TestTypeSP(); TestTypeKP(); TestTypeOrder(); - - Try { + + Try { trymat1(); trymat2(); trymat3(); @@ -554,6 +554,3 @@ time_lapse::~time_lapse() ///@} - - - diff --git a/src/newmat/tmt.h b/src/newmat/tmt.h index 4233d6d6..d9305452 100644 --- a/src/newmat/tmt.h +++ b/src/newmat/tmt.h @@ -38,11 +38,11 @@ class MultWithCarry }; // fill a matrix with values from the MultWithCarry random number generator -void FillWithValues(MultWithCarry& MWC, Matrix& M); -void FillWithValues(MultWithCarry& MWC, UpperTriangularMatrix& M); -void FillWithValues(MultWithCarry& MWC, LowerTriangularMatrix& M); -void FillWithValues(MultWithCarry& MWC, DiagonalMatrix& M); -void FillWithValues(MultWithCarry& MWC, SymmetricMatrix& M); +void FillWithValues(MultWithCarry& MWC, Matrix& M); +void FillWithValues(MultWithCarry& MWC, UpperTriangularMatrix& M); +void FillWithValues(MultWithCarry& MWC, LowerTriangularMatrix& M); +void FillWithValues(MultWithCarry& MWC, DiagonalMatrix& M); +void FillWithValues(MultWithCarry& MWC, SymmetricMatrix& M); void Print(const Matrix& X); diff --git a/src/newmat/tmt.txt b/src/newmat/tmt.txt index 76a79f79..55ad1126 100644 --- a/src/newmat/tmt.txt +++ b/src/newmat/tmt.txt @@ -6591,4 +6591,3 @@ Number of matrices tested = 1604 Number of non-zero matrices (should be 1) = 1 Elapsed (processor) time = 0.36 seconds - diff --git a/src/newmat/tmt1.cpp b/src/newmat/tmt1.cpp index 80d76ded..054b5db8 100644 --- a/src/newmat/tmt1.cpp +++ b/src/newmat/tmt1.cpp @@ -137,15 +137,15 @@ void trymat1() Matrix M1X = BM1; Matrix M2X = BM2; Matrix MX = BM; MX -= M1X + M2X; Print(MX); MX = BM1; MX += BM2; MX -= M1X; MX -= M2X; Print(MX); - SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t(); + SymmetricBandMatrix SM1; SM1 << BM1 * BM1.t(); SymmetricBandMatrix SM2; SM2 << BM2 * BM2.t(); SM1 *= 5.5; M1X *= M1X.t(); M1X *= 5.5; M2X *= M2X.t(); SM1 -= SM2; M1 = SM1 - M1X + M2X; Print(M1); - M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1); - M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1); - M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1); - + M1 = BM1; BM1 *= SM1; M1 = M1 * SM1 - BM1; Print(M1); + M1 = BM1; BM1 -= SM1; M1 = M1 - SM1 - BM1; Print(M1); + M1 = BM1; BM1 += SM1; M1 = M1 + SM1 - BM1; Print(M1); + } { Tracer et1("Stage 6"); @@ -271,16 +271,16 @@ void trymat1() ColumnVector cv(10); for (int i = 1; i <= 10; ++i) cv(i) = mwc.Next(); cv *= 100.0; - cv(1) -= 6.27874; - cv(2) -= 42.1718; - cv(3) -= 80.2854; - cv(4) -= 12.961; - cv(5) -= 17.7499; - cv(6) -= 13.2657; - cv(7) -= 50.4923; - cv(8) -= 26.095; - cv(9) -= 57.9147; - cv(10) -= 30.1778; + cv(1) -= 6.27874; + cv(2) -= 42.1718; + cv(3) -= 80.2854; + cv(4) -= 12.961; + cv(5) -= 17.7499; + cv(6) -= 13.2657; + cv(7) -= 50.4923; + cv(8) -= 26.095; + cv(9) -= 57.9147; + cv(10) -= 30.1778; Clean(cv, 0.0001); Print(cv); } diff --git a/src/newmat/tmt2.cpp b/src/newmat/tmt2.cpp index 443aea74..c435e675 100644 --- a/src/newmat/tmt2.cpp +++ b/src/newmat/tmt2.cpp @@ -285,7 +285,7 @@ void trymat2() I.ReSize(3); J.ReSize(3); K = I * J; N = K - I; Print(N); K << D; N = K - D; Print(N); } - + // test add integer { Matrix X(2,3); @@ -296,15 +296,15 @@ void trymat2() X += (-10); X -= Y; Print(X); - + // also test f suffix X << 5.25f << 7.75f << 1.25f << 9.00f << 1.00f << 2.50f; X -= Y; Print(X); - + } - - + + // cout << "\nEnd of second test\n"; diff --git a/src/newmat/tmt4.cpp b/src/newmat/tmt4.cpp index d2ed96bf..fed955b5 100644 --- a/src/newmat/tmt4.cpp +++ b/src/newmat/tmt4.cpp @@ -157,7 +157,7 @@ void trymat4() D1(1) = 18; D1(2) = 23; D1(3) = 31; D1(4) = 17; D1 -= D; Print(D1); } - + { Tracer et1("Stage 8"); // test swap functions @@ -173,13 +173,13 @@ void trymat4() D = A - B1; Print(D); Print(A); // now should be zero D.ReSize(1,2); D(1,1) = a; D(1,2) = b; Print(D); - + A.ReSize(20,20); FillWithValues(MWC, A); - + UpperTriangularMatrix UA; UA << A; UpperTriangularMatrix UA1 = UA; UpperTriangularMatrix UB; swap(UA, UB); Print(UA); UB -= UA1; Print(UB); - + LowerTriangularMatrix LA; LA << A; LowerTriangularMatrix LA1 = LA; LowerTriangularMatrix LB; swap(LB, LA); Print(LA); LB -= LA1; Print(LB); @@ -187,63 +187,63 @@ void trymat4() SymmetricMatrix SA; SA << A; SymmetricMatrix SA1 = SA; SymmetricMatrix SB; swap(SA, SB); Print(SA); SB -= SA1; Print(SB); - + DiagonalMatrix DA; DA << A; DiagonalMatrix DA1 = DA; DiagonalMatrix DB; swap(DB, DA); Print(DA); DB -= DA1; Print(DB); - + RowVector RVA = A.Row(1); RowVector RVA1 = RVA; RowVector RVB; swap(RVB, RVA); Print(RVA); RVB -= RVA1; Print(RVB); - + ColumnVector CVA = A.Column(1); ColumnVector CVA1 = CVA; ColumnVector CVB; swap(CVA, CVB); Print(CVA); CVB -= CVA1; Print(CVB); - + BandMatrix BA(20, 7, 4); BA.Inject(A); BandMatrix BA1 = BA; BandMatrix BB; swap(BA, BB); D = BA; Print(D); BB -= BA1; D = BB; Print(D); - + LowerBandMatrix LBA(20, 6); LBA.Inject(A); LowerBandMatrix LBA1 = LBA; LowerBandMatrix LBB; swap(LBB, LBA); D = LBA; Print(D); LBB -= LBA1; D = LBB; Print(D); - + UpperBandMatrix UBA(20, 9); UBA.Inject(A); UpperBandMatrix UBA1 = UBA; UpperBandMatrix UBB; swap(UBA, UBB); D = UBA; Print(D); UBB -= UBA1; D = UBB; Print(D); - + SymmetricBandMatrix SBA(20, 4); SBA.Inject(A); SymmetricBandMatrix SBA1 = SBA; SymmetricBandMatrix SBB; - + swap(SBB, SBA); D = SBA; Print(D); SBB -= SBA1; D = SBB; Print(D); - + B.ReSize(10,10); FillWithValues(MWC, B); - + CroutMatrix CA = A; IdentityMatrix IA(20); CroutMatrix CB = B; IdentityMatrix IB(10); swap(CA, CB); swap(IA, IB); D = CA.i() * B - IA; Clean(D,0.00000001); Print(D); D = CB.i() * A - IB; Clean(D,0.00000001); Print(D); - + BA.ReSize(20, 5, 7); BA.Inject(A); BandLUMatrix BLUA = BA; BB.ReSize(10, 3, 4); BB.Inject(B); BandLUMatrix BLUB = BB; swap(BLUA, BLUB); D = BLUA.i() * BB - IA; Clean(D,0.00000001); Print(D); D = BLUB.i() * BA - IB; Clean(D,0.00000001); Print(D); - + SBA.ReSize(20, 5); SBA.Inject(A); BandLUMatrix SBLUA = SBA; SBB.ReSize(10, 3); SBB.Inject(B); BandLUMatrix SBLUB = SBB; swap(SBLUA, SBLUB); D = SBLUA.i() * SBB - IA; Clean(D,0.00000001); Print(D); D = SBLUB.i() * SBA - IB; Clean(D,0.00000001); Print(D); - + UA << A; GenericMatrix GUA = UA; GenericMatrix GB = B; swap(GUA, GB); D = GB - UA; Print(D); D = B - GUA; Print(D); - + } // cout << "\nEnd of fourth test\n"; diff --git a/src/newmat/tmt6.cpp b/src/newmat/tmt6.cpp index 8defa2ad..f8423b72 100644 --- a/src/newmat/tmt6.cpp +++ b/src/newmat/tmt6.cpp @@ -79,7 +79,7 @@ void trymat6() for (i=1;i<=6;i++) for (j=i;j<=6;j++) S(i,j)=i*i+2.0+j; Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S; Matrix M(6,6); - for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0; + for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0; { Tracer et1("Stage 1"); Print(Matrix(MS+(-MS))); diff --git a/src/newmat/tmt7.cpp b/src/newmat/tmt7.cpp index c4e43fa3..4ece58cc 100644 --- a/src/newmat/tmt7.cpp +++ b/src/newmat/tmt7.cpp @@ -29,7 +29,7 @@ void c_matrix_multiply(int p, int q, int r, c[i][k] = sum; } } - + void trymat7() @@ -50,7 +50,7 @@ void trymat7() Matrix MD=D; Matrix ML=L; Matrix MU=U; Matrix MS=S; Matrix M(6,6); - for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0; + for (i=1;i<=6;i++) for (j=1;j<=6;j++) M(i,j)=i*j+i*i-10.0; { Tracer et1("Stage 1"); Print(Matrix((S-M)-(MS-M))); @@ -99,7 +99,7 @@ void trymat7() CV << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75; D << 2 << 6 <<3 << 8 << -4 << 17.5 << 2 << 1 << -2 << 5 << 3.75; X = CV.AsDiagonal(); X = X-D; Print(X); - SymmetricBandMatrix SB1(11,7); SB1 = 5; + SymmetricBandMatrix SB1(11,7); SB1 = 5; SymmetricBandMatrix SB2 = SB1 + D; X.ReSize(11,11); X=0; for (i=1;i<=11;i++) for (j=1;j<=11;j++) @@ -145,12 +145,12 @@ void trymat7() for (i = 0; i < 5; i++) if (F[i] != 0) Test(9)=1; Print(Test); } - + { // testing RealStarStar Tracer et("Stage 6"); MultWithCarry MWC; - + Matrix A(10, 12), B(12, 15), C(10, 15); FillWithValues(MWC, A); FillWithValues(MWC, B); ConstRealStarStar a(A); @@ -165,13 +165,13 @@ void trymat7() ConstRealStarStar(A),ConstRealStarStar(B),RealStarStar(C)); X = C - A * B; Clean(X,0.00000001); Print(X); } - + { // testing resize_keep Tracer et("Stage 7"); Matrix X, Y; MultWithCarry MWC; - + X.resize(20,35); FillWithValues(MWC, X); Matrix M(20,35); M = X; X = M.submatrix(1,15,1,25); @@ -181,14 +181,14 @@ void trymat7() M.resize_keep(29,27); Y -= M; Print(Y); M.resize_keep(0,5); M.resize_keep(10,10); Print(M); M.resize_keep(15,0); M.resize_keep(10,10); Print(M); - + X.resize(20,35); FillWithValues(MWC, X); M = X; M.resize_keep(38,17); Y.resize(38,17); Y = 0; Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17); Y -= M; Print(Y); - + X.resize(40,12); FillWithValues(MWC, X); M = X; M.resize_keep(38,17); @@ -207,23 +207,23 @@ void trymat7() nM.resize_keep(29,27); Y -= nM; Print(Y); nM.resize_keep(0,5); nM.resize_keep(10,10); Print(nM); nM.resize_keep(15,0); nM.resize_keep(10,10); Print(nM); - + X.resize(20,35); FillWithValues(MWC, X); nM = X; nM.resize_keep(38,17); Y.resize(38,17); Y = 0; Y.submatrix(1,20,1,17) = X.submatrix(1,20,1,17); Y -= nM; Print(Y); - + X.resize(40,12); FillWithValues(MWC, X); nM = X; nM.resize_keep(38,17); Y.resize(38,17); Y = 0; Y.submatrix(1,38,1,12) = X.submatrix(1,38,1,12); - Y -= nM; Print(Y); + Y -= nM; Print(Y); + +#endif -#endif - X.resize(20,20); FillWithValues(MWC, X); SquareMatrix SQM(20); SQM << X; X = SQM.sym_submatrix(1,13); @@ -232,7 +232,7 @@ void trymat7() Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X; SQM.resize_keep(23,23); Y -= SQM; Print(Y); SQM.resize_keep(0); SQM.resize_keep(50); Print(SQM); - + X.resize(20,20); FillWithValues(MWC, X); SymmetricMatrix SM(20); SM << X; X = SM.sym_submatrix(1,13); @@ -241,7 +241,7 @@ void trymat7() Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X; SM.resize_keep(23); Y -= SM; Print(Y); SM.resize_keep(0); SM.resize_keep(50); Print(SM); - + X.resize(20,20); FillWithValues(MWC, X); LowerTriangularMatrix LT(20); LT << X; X = LT.sym_submatrix(1,13); @@ -250,7 +250,7 @@ void trymat7() Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X; LT.resize_keep(23); Y -= LT; Print(Y); LT.resize_keep(0); LT.resize_keep(50); Print(LT); - + X.resize(20,20); FillWithValues(MWC, X); UpperTriangularMatrix UT(20); UT << X; X = UT.sym_submatrix(1,13); @@ -259,7 +259,7 @@ void trymat7() Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X; UT.resize_keep(23); Y -= UT; Print(Y); UT.resize_keep(0); UT.resize_keep(50); Print(UT); - + X.resize(20,20); FillWithValues(MWC, X); DiagonalMatrix DM(20); DM << X; X = DM.sym_submatrix(1,13); @@ -268,7 +268,7 @@ void trymat7() Y.resize(23,23); Y = 0; Y.sym_submatrix(1,13) = X; DM.resize_keep(23); Y -= DM; Print(Y); DM.resize_keep(0); DM.resize_keep(50); Print(DM); - + X.resize(1,20); FillWithValues(MWC, X); RowVector RV(20); RV << X; X = RV.columns(1,13); @@ -277,7 +277,7 @@ void trymat7() Y.resize(1,23); Y = 0; Y.columns(1,13) = X; RV.resize_keep(1,23); Y -= RV; Print(Y); RV.resize_keep(0); RV.resize_keep(50); Print(RV); - + X.resize(20,1); FillWithValues(MWC, X); ColumnVector CV(20); CV << X; X = CV.rows(1,13); @@ -286,9 +286,9 @@ void trymat7() Y.resize(23,1); Y = 0; Y.rows(1,13) = X; CV.resize_keep(23,1); Y -= CV; Print(Y); CV.resize_keep(0); CV.resize_keep(50); Print(CV); - - - } + + + } // cout << "\nEnd of seventh test\n"; diff --git a/src/newmat/tmt8.cpp b/src/newmat/tmt8.cpp index 60b5a16b..f369e053 100644 --- a/src/newmat/tmt8.cpp +++ b/src/newmat/tmt8.cpp @@ -145,7 +145,7 @@ void trymat8() M.ReSize(MX); DCR( D.nric(), C.nric(), 5, R.nric(), 9, M.nric() ); M -= MX; Print(M); - + // test swap nricMatrix A(3,4); nricMatrix B(4,5); A.Row(1) << 2 << 7 << 3 << 6; @@ -164,7 +164,7 @@ void trymat8() for (int k = 1; k <= 4; ++k) X.nric()[i][j] += A1.nric()[i][k] * B1.nric()[k][j]; } - X1 -= X; Print(X1); + X1 -= X; Print(X1); } #endif { diff --git a/src/newmat/tmt9.cpp b/src/newmat/tmt9.cpp index bef46361..fc4fd4cc 100644 --- a/src/newmat/tmt9.cpp +++ b/src/newmat/tmt9.cpp @@ -41,7 +41,7 @@ void trymat9() { Tracer et1("Stage 2"); Matrix DXE = D.i() * X * E; - DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE); + DXE = E.i() * DXE.t() * D - X.t(); Clean(DXE, 0.00000001); Print(DXE); E=D; for (i=1;i<=7;i++) E(i,i)=i*3+1; } DiagonalMatrix F=D; @@ -216,7 +216,7 @@ void trymat9() { Tracer et1("Stage 8"); - // test copying CroutMatrix and BandLUMatrix - see also tmtd.cpp + // test copying CroutMatrix and BandLUMatrix - see also tmtd.cpp MultWithCarry MWC; SymmetricBandMatrix SBM(50, 10); for (int i = 1; i <= 50; ++i) for (int j = 1; j <= i; ++j) @@ -241,7 +241,7 @@ void trymat9() { Tracer et1("Stage 9"); - // do it again odd matrix size + // do it again odd matrix size MultWithCarry MWC; SymmetricBandMatrix SBM(51, 10); for (int i = 1; i <= 51; ++i) for (int j = 1; j <= i; ++j) @@ -263,7 +263,7 @@ void trymat9() D(4) = z.sign() - x.sign(); Clean(D,0.00000001); Print(D); } - + { Tracer et1("Stage 10"); // testing equal again @@ -280,12 +280,12 @@ void trymat9() B *= 3; A = B; X = A.i() * B - IdentityMatrix(10); Clean(X,0.00000001); Print(X); - } - - - - - + } + + + + + // cout << "\nEnd of ninth test\n"; } diff --git a/src/newmat/tmtc.cpp b/src/newmat/tmtc.cpp index a88f4f82..2d485957 100644 --- a/src/newmat/tmtc.cpp +++ b/src/newmat/tmtc.cpp @@ -184,7 +184,7 @@ void trymatc() { RowVector RV = 5 * M1.Row(i); M.Row(i) -= RV; } M += M1; Print(M); - + } { @@ -205,13 +205,13 @@ void trymatc() nricMatrix N2 = N; nricMatrix Y = N; Print(Y); Y = N1 - N2; Print(Y); - + N = M1 / 2; N1 = N * 2; N.Release(); N2 = N * 2; Y = N; Print(N); Y = (N1 - M1) | (N2 - M1); Print(Y); #endif } - + { Tracer et("Stage 13"); // test sum of squares of rows or columns @@ -272,9 +272,9 @@ void trymatc() A1.ReSize(10,0); X.ReSize(10,1); X = 0; X -= A1.sum_square_rows(); Clean(X, 0.00000000000001); Print(X); X.ReSize(1,0); X -= A1.sum_square_columns(); Clean(X, 0.00000000000001); Print(X); - + } - + { Tracer et("Stage 14"); // test extend orthonormal @@ -289,7 +289,7 @@ void trymatc() // Check orthogonality X = A.t() * A - IdentityMatrix(5); Clean(X, 0.000000001); Print(X); - // Try orthonality extend + // Try orthonality extend SquareMatrix A1(20); A1.Columns(1,5) = A; extend_orthonormal(A1,5); @@ -297,10 +297,10 @@ void trymatc() X = A - A1.Columns(1,5); Print(X); // Check orthogonality X = A1.t() * A1 - IdentityMatrix(20); - Clean(X, 0.000000001); Print(X); + Clean(X, 0.000000001); Print(X); X = A1 * A1.t() - IdentityMatrix(20); Clean(X, 0.000000001); Print(X); - // Test with smaller number of columns + // Test with smaller number of columns Matrix A2(20,15); A2.Columns(1,5) = A; extend_orthonormal(A2,5); @@ -333,7 +333,7 @@ void trymatc() B -= A; Clean(B, 0.000000001); Print(B); } - + { Tracer et("Stage 16"); // test sum of rows or columns @@ -396,9 +396,9 @@ void trymatc() A1.ReSize(10,0); X.ReSize(10,1); X = 0; X -= A1.sum_rows(); Print(X); X.ReSize(1,0); X -= A1.sum_columns(); Print(X); - + } - + { Tracer et("Stage 17"); // SP_eq on submatrices @@ -408,20 +408,20 @@ void trymatc() Matrix X = A; X.submatrix(5,21,3,13).SP_eq(B); Matrix Y = A; Y.submatrix(5,21,3,13) = SP(Y.submatrix(5,21,3,13),B); Y -= X; Print(Y); - + UpperTriangularMatrix UT(33); FillWithValues(mwc, UT); UpperTriangularMatrix UTX = UT; UTX.submatrix(5,21,3,13).SP_eq(B); Y = UT; Y.submatrix(5,21,3,13) = SP(Y.submatrix(5,21,3,13),B); Y -= UTX; Print(Y); - + UT.resize(10); FillWithValues(mwc, UT); X = A; X.submatrix(5,14,3,12).SP_eq(UT); Y = A; Y.submatrix(5,14,3,12) = SP(Y.submatrix(5,14,3,12),UT); Y -= X; Print(Y); - + SymmetricMatrix SM(30), SM1(10); FillWithValues(mwc, SM); FillWithValues(mwc, SM1); SymmetricMatrix SMX = SM; @@ -429,12 +429,12 @@ void trymatc() //SMX.sym_submatrix(5,14) += SM1; //Y = SM; Y.submatrix(5,14,5,14) = Y.submatrix(5,14,5,14) + SM1; //Y -= SMX; Print(Y); - + SMX = SM; SMX.submatrix(5,14,5,14).SP_eq(SM1); Y = SM; Y.submatrix(5,14,5,14) = SP(Y.submatrix(5,14,5,14),SM1); Y -= SMX; Print(Y); - + BandMatrix BM(50, 10, 12); X.resize(50,50); FillWithValues(mwc, X); BM.inject(X); X = BM; @@ -449,17 +449,17 @@ void trymatc() BM.submatrix(2,8,4,13).SP_eq(B); X.submatrix(2,8,4,13) = SP(X.submatrix(2,8,4,13), B); X -= BM; Print(X); - + X.resize(50,50); FillWithValues(mwc, X); BM.inject(X); X = BM; B.resize(21, 23); FillWithValues(mwc, B); BM.submatrix(2,22,4,26).SP_eq(B); X.submatrix(2,22,4,26) = SP(X.submatrix(2,22,4,26), B); X -= BM; Print(X); - + } - - { + + { Tracer et("Stage 18"); // += real on submatrices of symmetric matrices // these don't work either @@ -470,27 +470,27 @@ void trymatc() //SMX.sym_submatrix(7,11) += 5; //Matrix X = SM; X.sym_submatrix(7,11) += 5; //X -= SMX; Print(X); - + } - + { Tracer et("Stage 19"); // SQ_eq on matrices MultWithCarry mwc; - + Matrix X(10,13), Y(10,13); FillWithValues(mwc, X); FillWithValues(mwc, Y); Matrix Z = X; Z.SP_eq(Y); Z -= SP(X, Y); Print(Z); - + Z = X; Z.SP_eq(2 * Y); Z -= 2 * SP(X, Y); Print(Z); - + Z = X; Z.SP_eq(Z); Z -= SP(X, X); Print(Z); - + GenericMatrix GX = X, GY = Y; test_SP_eq(GX, GY); test_SP_eq(GX, GX); - + X.resize(17,17); Y.resize(17,17); FillWithValues(mwc, X); FillWithValues(mwc, Y); UpperTriangularMatrix UT; UT << Y; @@ -504,28 +504,28 @@ void trymatc() UpperTriangularMatrix UX(19), UY(19); FillWithValues(mwc, UX); FillWithValues(mwc, UY); - + GX = UX; GY = UY; test_SP_eq(GX, GY); test_SP_eq(GX, GY.t()); UpperTriangularMatrix UZ = UX; UZ.SP_eq(UY); UZ -= SP(UX, UY); Print(UZ); - + UZ = UX; UZ.SP_eq(2 * UY); UZ -= 2 * SP(UX, UY); Print(UZ); - + UT << UY; LT << UY; D << UY; UpperTriangularMatrix UZUT = UX; UZUT.SP_eq(UT); UpperTriangularMatrix UZLT = UX; UZLT.SP_eq(LT); UpperTriangularMatrix UZD = UX; UZD.SP_eq(D); Z = UZUT + UZLT - UZD - SP(UX, UY); Clean(Z, 0.000000001); Print(Z); - + SymmetricMatrix SM(13), SM1(13), SMX, SM2; FillWithValues(mwc, SM); FillWithValues(mwc, SM1); SMX = SP(SM, SM1); SM2 = SM; SM2.SP_eq(SM1); SM2 -= SMX; Print(SM2); GX = SM; GY = SM1; test_SP_eq(GX, GY); - + X.resize(13,13); FillWithValues(mwc, X); Z = SP(X,SM); SM2 = SM; X.SP_eq(SM); X -= Z; Print(X); @@ -533,7 +533,7 @@ void trymatc() FillWithValues(mwc, X); GX = X; GY = SM; test_SP_eq(GX, GY); GX = SM; GY = X; test_SP_eq(GX, GY); - + BandMatrix BMX(13, 5,9); BMX.inject(X); SymmetricBandMatrix SBM(13,6); SBM.inject(SM); GX = BMX; GY = SBM; test_SP_eq(GX, GY); @@ -542,11 +542,11 @@ void trymatc() BandMatrix BMX2 = SBM; GX = BMX; GY = BMX2; test_SP_eq(GX, GY); BMX1 = BMX; BMX1.SP_eq(BMX2); Z = BMX1 - SP(BMX,BMX2); Print(Z); - - + + } - - + + // cout << "\nEnd of twelfth test\n"; } diff --git a/src/newmat/tmtd.cpp b/src/newmat/tmtd.cpp index edf574e9..6ff3cd74 100644 --- a/src/newmat/tmtd.cpp +++ b/src/newmat/tmtd.cpp @@ -137,9 +137,9 @@ void QRZ_checkT(Matrix& X, Matrix& Y, UpperTriangularMatrix& U, Matrix& M) Matrix XT = X.t(), YT = Y.t(), MT, BX; LowerTriangularMatrix L; QRZ(X, Y, U, M); QRZT(XT, YT, L, MT); - BX = U - L.t(); Clean(BX, 0.000000001); Print(BX); - BX = M - MT.t(); Clean(BX, 0.000000001); Print(BX); -} + BX = U - L.t(); Clean(BX, 0.000000001); Print(BX); + BX = M - MT.t(); Clean(BX, 0.000000001); Print(BX); +} @@ -289,7 +289,7 @@ void trymatd() M = A * Cholesky(B); M = M * M.t() - A * B * A; Clean(M, 0.000000001); Print(M); } - + { Tracer et1("Stage 6"); int N=49; @@ -334,8 +334,8 @@ void trymatd() D(6) = B1.nrows(); D(7) = B1.ncols(); D(8) = B1.size(); Clean(C, 0.000000001); Print(C); // see if we get an implicit invert - B1 = -A; - D(9) = determinant(B1) / x + 1; // odd number of rows - sign will change + B1 = -A; + D(9) = determinant(B1) / x + 1; // odd number of rows - sign will change C = -A * Inverter2(B1) - IdentityMatrix(7); Clean(C, 0.000000001); Print(C); // check for_return @@ -368,7 +368,7 @@ void trymatd() const int* i2 = B2.const_data_indx(); D(25) = (d != d1 ? 0 : 1) + (d1 == d2 ? 0 : 1) + (i != i1 ? 0 : 1) + (i1 == i2 ? 0 : 1); - + Clean(D, 0.000000001); Print(D); } @@ -403,8 +403,8 @@ void trymatd() D(6) = B1.nrows(); D(7) = B1.ncols(); D(8) = B1.size(); Clean(C, 0.000000001); Print(C); // see if we get an implicit invert - B1 = -A; - D(9) = determinant(B1) / x + 1; // odd number of rows - sign will change + B1 = -A; + D(9) = determinant(B1) / x + 1; // odd number of rows - sign will change C = -A * Inverter2(B1) - IdentityMatrix(7); Clean(C, 0.000000001); Print(C); // check for_return @@ -471,7 +471,7 @@ void trymatd() QRZ(X, U2); Matrix Diff = U1 - U2; Clean(Diff, 0.000000001); Print(Diff); - // Try adding new row to X and updating triangular matrix + // Try adding new row to X and updating triangular matrix RowVector NewRow(10); for (j = 1; j <= 10; ++j) NewRow(j) = 2.0 * (mwc.Next() - 0.5); UpdateCholesky(U2, NewRow); @@ -492,7 +492,7 @@ void trymatd() CircularShift(X, 1,6); CircularShift(X, 6,10); } - + { Tracer et1("Stage 10"); // Try updating QRZ, QRZT decomposition @@ -507,9 +507,9 @@ void trymatd() tuqrz2.Reset(); tuqrz2.ClearRow(5); tuqrz2.ClearRow(6); tuqrz2.DoTest(); tuqrz2.Reset(); tuqrz2.ClearRow(15); tuqrz2.DoTest(); TestUpdateQRZ tuqrz3(5, 0, 10, 0); tuqrz3.DoTest(); - + } - + { Tracer et1("Stage 11"); // combining three decompositions @@ -537,14 +537,14 @@ void trymatd() fix_signs(U3); BX = U - U3; Clean(BX, 0.000000001); Print(BX); - + // same again with some columns set to zero FillWithValues(MWC, X1); FillWithValues(MWC, Y1); FillWithValues(MWC, X2); FillWithValues(MWC, Y2); FillWithValues(MWC, X3); FillWithValues(MWC, Y3); X1.columns(1,8) = 0; X2.columns(5,13) = 0; - X3.columns(10,20) = 0; + X3.columns(10,20) = 0; X = X1 & X2 & X3; Y = Y1 & Y2 & Y3; QRZ(X1, Y1, U1, M1); QRZ(X2, Y2, U2, M2); @@ -560,14 +560,14 @@ void trymatd() fix_signs(U3); BX = U - U3; Clean(BX, 0.000000001); Print(BX); - + // same again with some duplicate columns FillWithValues(MWC, X1); FillWithValues(MWC, Y1); FillWithValues(MWC, X2); FillWithValues(MWC, Y2); FillWithValues(MWC, X3); FillWithValues(MWC, Y3); X1.columns(4,5) = X1.columns(1,2); X2.columns(7,12) = X2.columns(1,6); - X3.columns(11,20) = X3.columns(1,10); + X3.columns(11,20) = X3.columns(1,10); X = X1 & X2 & X3; Y = Y1 & Y2 & Y3; QRZ(X1, Y1, U1, M1); QRZ(X2, Y2, U2, M2); @@ -583,9 +583,9 @@ void trymatd() fix_signs(U3); BX = U - U3; Clean(BX, 0.000000001); Print(BX); - - - + + + // Try with incremental update FillWithValues(MWC, X1); FillWithValues(MWC, Y1); FillWithValues(MWC, X2); FillWithValues(MWC, Y2); @@ -608,7 +608,7 @@ void trymatd() FillWithValues(MWC, X3); FillWithValues(MWC, Y3); X1.columns(1,8) = 0; X2.columns(5,13) = 0; - X3.columns(10,20) = 0; + X3.columns(10,20) = 0; X = X1 & X2 & X3, Y = Y1 & Y2 & Y3; QRZ(X1, Y1, U1, M1); UpdateQRZ(X2, U1); UpdateQRZ(X2, Y2, M1); @@ -620,15 +620,15 @@ void trymatd() fix_signs(U1); BX = U - U1; Clean(BX, 0.000000001); Print(BX); - - + + // same again with some duplicate columns FillWithValues(MWC, X1); FillWithValues(MWC, Y1); FillWithValues(MWC, X2); FillWithValues(MWC, Y2); FillWithValues(MWC, X3); FillWithValues(MWC, Y3); X1.columns(4,5) = X1.columns(1,2); X2.columns(7,12) = X2.columns(1,6); - X3.columns(11,20) = X3.columns(1,10); + X3.columns(11,20) = X3.columns(1,10); X = X1 & X2 & X3, Y = Y1 & Y2 & Y3; QRZ(X1, Y1, U1, M1); UpdateQRZ(X2, U1); UpdateQRZ(X2, Y2, M1); @@ -640,7 +640,7 @@ void trymatd() fix_signs(U1); BX = U - U1; Clean(BX, 0.000000001); Print(BX); - + // Try with incremental update using QRZT LowerTriangularMatrix L, L1; X1.resize(21, 90); X2.resize(21, 250); X3.resize(21, 113); @@ -662,7 +662,7 @@ void trymatd() FillWithValues(MWC, X3); X1.rows(1,8) = 0; X2.rows(5,13) = 0; - X3.rows(10,21) = 0; + X3.rows(10,21) = 0; X = X1 | X2 | X3; QRZT(X1, L1); UpdateQRZT(X2, L1); @@ -671,15 +671,15 @@ void trymatd() fix_signs(L1); BX = L - L1; Clean(BX, 0.000000001); Print(BX); - - + + // same again with some duplicate columns FillWithValues(MWC, X1); FillWithValues(MWC, X2); FillWithValues(MWC, X3); X1.columns(4,5) = X1.columns(1,2); X2.columns(7,12) = X2.columns(1,6); - X3.columns(12,21) = X3.columns(1,10); + X3.columns(12,21) = X3.columns(1,10); X = X1 | X2 | X3; QRZT(X1, L1); UpdateQRZT(X2, L1); @@ -688,9 +688,9 @@ void trymatd() fix_signs(L1); BX = L - L1; Clean(BX, 0.000000001); Print(BX); - } - - + } + + // cout << "\nEnd of Thirteenth test\n"; } @@ -764,6 +764,6 @@ void fix_signs(LowerTriangularMatrix& L) L.column(i) = -L.column(i); } } - + ///@} diff --git a/src/newmat/tmte.cpp b/src/newmat/tmte.cpp index c3e46ad1..215e50a4 100644 --- a/src/newmat/tmte.cpp +++ b/src/newmat/tmte.cpp @@ -502,7 +502,7 @@ void trymate() Clean(X, 0.000000001); Print(X); Y = V * V.t() - IdentityMatrix(2); Clean(Y, 0.000000001); Print(Y); - + } diff --git a/src/newmat/tmtf.cpp b/src/newmat/tmtf.cpp index 57d21c20..8ddd640b 100644 --- a/src/newmat/tmtf.cpp +++ b/src/newmat/tmtf.cpp @@ -303,7 +303,7 @@ static void test9(int m, int n) Clean(X,0.000000001); Clean(Y,0.000000001); Print(X); Print(Y); FFT2I(X1, Y1, X1, Y1); X1 -= A; Y1 -= B; - Clean(X1,0.000000001); Clean(Y1,0.000000001); Print(X1); Print(Y1); + Clean(X1,0.000000001); Clean(Y1,0.000000001); Print(X1); Print(Y1); } @@ -367,7 +367,7 @@ void trymatf() test9(1,16); test9(16,1); test9(4,3); test9(4,4); test9(4,5); test9(5,3); - + // now do the same thing all over again forcing use of old FFT FFT_Controller::OnlyOldFFT = true; @@ -417,7 +417,7 @@ void trymatf() // compare DCT, DST and slow versions test8(2); test8(26); test8(32); test8(18); - + // compare FFT2 with slow version // don't redo these diff --git a/src/newmat/tmth.cpp b/src/newmat/tmth.cpp index 0d3100e4..ef77cbdd 100644 --- a/src/newmat/tmth.cpp +++ b/src/newmat/tmth.cpp @@ -766,7 +766,7 @@ void trymath() RV(2) = BM.BandWidth().Upper(); Print(RV); } - + { // Helmert multiplies Tracer et1("Stage 9"); @@ -779,22 +779,22 @@ void trymath() Matrix H = Helmert(9); H -= Y.t(); Clean(H,0.000000001); Print(H); Matrix Z = Helmert(Y) - I; Clean(Z,0.000000001); Print(Z); - + Matrix A(9, 8); for (i = 1; i <= 9; ++i) for (j = 1; j <= 8; ++j) A(i, j) = Helmert_transpose(X.column(j), i); - A -= Y; Clean(A,0.000000001); Print(A); - + A -= Y; Clean(A,0.000000001); Print(A); + X = I; Y = Helmert_transpose(X, true); H = Helmert(8, true); H -= Y.t(); Clean(H,0.000000001); Print(H); Z = Helmert(Y, true) - I; Clean(Z,0.000000001); Print(Z); - + A.resize(8, 8); for (i = 1; i <= 8; ++i) for (j = 1; j <= 8; ++j) A(i, j) = Helmert_transpose(X.column(j), i, true); - A -= Y; Clean(A,0.000000001); Print(A); + A -= Y; Clean(A,0.000000001); Print(A); @@ -804,7 +804,7 @@ void trymath() H = Helmert(9, true); H -= Y; Clean(H,0.000000001); Print(H); Z = Helmert_transpose(Y, true) - I; Clean(Z,0.000000001); Print(Z); - + A.ReSize(9, 9); for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i, true); A -= Y; Clean(A,0.000000001); Print(A); @@ -813,10 +813,10 @@ void trymath() A.ReSize(8, 9); for (i = 1; i <= 9; ++i) A.Column(i) = Helmert(9, i); A -= Y; Clean(A,0.000000001); Print(A); - + ColumnVector Twos(100); Twos = 2; ColumnVector CV = Helmert(Twos); Clean(CV,0.000000001); Print(CV); - + X.resize(25,30); FillWithValues(MCW, X); Y = Helmert(X); @@ -824,7 +824,7 @@ void trymath() Clean(Z,0.000000001); Print(Z); Z = Helmert(X,true).row(25) - X.sum_columns() / 5.0; Clean(Z,0.000000001); Print(Z); - + I.resize(15); X = I; Z = Helmert_transpose(X, true) - Helmert(X, true).t(); @@ -832,15 +832,15 @@ void trymath() I.resize(14); Y = I; Z = Helmert(X) - Helmert_transpose(Y).t(); Clean(Z,0.000000001); Print(Z); - - - + + + } - - - - - + + + + + // cout << "\nEnd of Seventeenth test\n"; } diff --git a/src/newmat/tmti.cpp b/src/newmat/tmti.cpp index c74057e6..577e3e15 100644 --- a/src/newmat/tmti.cpp +++ b/src/newmat/tmti.cpp @@ -237,4 +237,3 @@ void trymati() ///@} - diff --git a/src/newmat/tmtj.cpp b/src/newmat/tmtj.cpp index 9463f0e9..8e5a9977 100644 --- a/src/newmat/tmtj.cpp +++ b/src/newmat/tmtj.cpp @@ -142,7 +142,7 @@ void trymatj() SymmetricBandMatrix SB(7,3); SymmetricMatrix S; S << (B+B.t()); SB.Inject(S); A = BA; S = SB; - Matrix X; + Matrix X; X = SP(BA,SB); X=X-SP(A,S); Print(X); X = SP(BA*2,SB); X=X-SP(A,S*2); Print(X); X = SP(BA,SB/4); X=X-SP(A/4,S); Print(X); diff --git a/src/newmat/tmtk.cpp b/src/newmat/tmtk.cpp index 4a57bc22..f1db617e 100644 --- a/src/newmat/tmtk.cpp +++ b/src/newmat/tmtk.cpp @@ -61,7 +61,7 @@ void trymatk() A1.CleanUp(); B1.CleanUp(); cout << "LowerTriangularMatrix\n"; - LowerTriangularMatrix A2(35), B2(35); + LowerTriangularMatrix A2(35), B2(35); for (i=0; i<35; i++) for (j=0; j<=i; j++) { A2[i][j] = i+100*j; B2(i+1,j+1) = i+100*j; } X = A2 - B2; Print(X); Y = X; @@ -176,7 +176,7 @@ void trymatk() } Print(X); Print(Y); Aa.CleanUp(); Ba.CleanUp(); - + // test special constructors used in Numerical Recipes for C++ Real a[] = {1.2, 5.6, 7.9, 3.8, 4.5, 1.3, 5.2, 9.9, 2.1, 4.7, 0.0, 1.6 }; diff --git a/src/newmat/tmtm.cpp b/src/newmat/tmtm.cpp index 057a28c3..9e103d63 100644 --- a/src/newmat/tmtm.cpp +++ b/src/newmat/tmtm.cpp @@ -264,7 +264,7 @@ void trymatm() M = AB; C -= M; Print(C); C << M; C += -M; Print(C); - + } @@ -272,8 +272,3 @@ void trymatm() ///@} - - - - -