Skip to content

GitLab

Commit 72622075 authored by Andreas Pfeiffer's avatar Andreas Pfeiffer
Browse files

merged changes from 1.8.2.1 (double[] -> std::vector), fixed bug in sub

parent e08afde6
Expand all Hide whitespace changes
Inline Side-by-side
Matrix/Matrix/DiagMatrix.h
View file @ 72622075
// -*- C++ -*-
// CLASSDOC OFF
// $Id: DiagMatrix.h,v 1.3 2003/10/23 21:29:50 garren Exp $
// $Id: DiagMatrix.h,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
// CLASSDOC ON
//
......@@ -51,6 +51,8 @@
#pragma interface
#endif
#include <vector>
#include "CLHEP/Matrix/defs.h"
#include "CLHEP/Matrix/GenMatrix.h"
......@@ -212,7 +214,7 @@ private:
friend HepMatrix operator*(const HepMatrix &m1, const HepDiagMatrix &m2);
friend HepVector operator*(const HepDiagMatrix &m1, const HepVector &m2);
double *m;
std::vector<double,Alloc<double,25> > m;
int nrow;
#if defined(__sun) || !defined(__GNUG__)
//
......
Matrix/Matrix/DiagMatrix.icc
View file @ 72622075
// -*- C++ -*-
// $Id: DiagMatrix.icc,v 1.2 2003/07/18 05:31:48 garren Exp $
// $Id: DiagMatrix.icc,v 1.2.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -55,7 +55,7 @@ inline double & HepDiagMatrix::fast(int row,int col)
if (row != col)
error("Index error in HepDiagMatrix::fast(i,j): i != j");
return *(m+(col-1));
return *(m.begin()+(col-1));
}
inline const double & HepDiagMatrix::fast(int row,int col) const
......@@ -65,7 +65,7 @@ inline const double & HepDiagMatrix::fast(int row,int col) const
error("Range error in HepDiagMatrix::fast()");
#endif
if (row == col) {
return *(m+(col-1));
return *(m.begin()+(col-1));
} else {
#if defined(__sun) || !defined(__GNUG__)
//
......
Matrix/Matrix/GenMatrix.h
View file @ 72622075
// -*- C++ -*-
// CLASSDOC OFF
// $Id: GenMatrix.h,v 1.3 2003/10/23 21:29:50 garren Exp $
// $Id: GenMatrix.h,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
// CLASSDOC ON
//
......@@ -49,6 +49,8 @@
#pragma interface
#endif
#include <vector>
#include <iostream>
#include "CLHEP/Matrix/defs.h"
......@@ -63,10 +65,44 @@ class HepGenMatrix;
* @ingroup matrix
*/
class HepGenMatrix {
public:
virtual ~HepGenMatrix() {}
template <class T, size_t size> class Alloc
{
public:
typedef T value_type;
typedef size_t size_type;
typedef ptrdiff_t difference_type;
typedef T* pointer;
typedef const T* const_pointer;
typedef T& reference;
typedef const T& const_reference;
pointer address(reference r) const { return &r; }
const_pointer address(const_reference r) const { return &r; }
Alloc() throw() {}
Alloc(const Alloc<T,size>&) throw() {}
~Alloc() throw() {}
pointer allocate(size_type n ) { if( n <= size ) return pool; else return new T[n]; }
void deallocate(pointer p, size_type n) { if (p == pool ) return; delete [] p; }
void construct(pointer p, const T& val ) { new(p) T(val); }
void destroy(pointer p) { p->~T(); }
size_type max_size() const throw() { size_type c = (size_type)(-1) /sizeof(T); return (0 < c ? c : 1); }
template<class O> struct rebind { typedef Alloc<O,size> other; };
private:
T pool[size];
};
typedef std::vector<double,Alloc<double,25> >::iterator mIter;
typedef std::vector<double,Alloc<double,25> >::const_iterator mcIter;
virtual int num_row() const = 0;
virtual int num_col() const = 0;
......@@ -103,7 +139,7 @@ public:
// ** Note that the indexing starts from [0][0]. **
inline static void swap(int&,int&);
inline static void swap(double *&, double *&);
inline static void swap(std::vector<double,Alloc<double,25> >&, std::vector<double,Alloc<double,25> >&);
virtual bool operator== ( const HepGenMatrix& ) const;
// equality operator for matrices (BaBar)
......@@ -131,8 +167,8 @@ private:
friend class HepGenMatrix_row;
friend class HepGenMatrix_row_const;
double data_array[size_max];
//-ap: removed this as it is taken over by the std::vector<double>
//-ap double data_array[size_max];
};
double norm(const HepGenMatrix &m);
......
Matrix/Matrix/GenMatrix.icc
View file @ 72622075
// -*- C++ -*-
// $Id: GenMatrix.icc,v 1.2 2003/07/18 05:31:48 garren Exp $
// $Id: GenMatrix.icc,v 1.2.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -41,7 +41,9 @@ namespace CLHEP {
// swap
//
inline void HepGenMatrix::swap(int &i,int &j) {int t=i;i=j;j=t;}
inline void HepGenMatrix::swap(double * &i,double * &j) {double * t=i;i=j;j=t;}
inline void HepGenMatrix::swap(std::vector<double,Alloc<double,25> >& i, std::vector<double,Alloc<double,25> >& j) {
std::vector<double,Alloc<double,25> > t=i;i=j;j=t;
}
//
// operator [] (I cannot make it virtual because return types are different.)
......
Matrix/Matrix/Matrix.h
View file @ 72622075
// -*- C++ -*-
// CLASSDOC OFF
// $Id: Matrix.h,v 1.3 2003/10/23 21:29:50 garren Exp $
// $Id: Matrix.h,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
// CLASSDOC ON
//
......@@ -228,6 +228,8 @@
#pragma interface
#endif
#include <vector>
#include "CLHEP/Matrix/defs.h"
#include "CLHEP/Matrix/GenMatrix.h"
......@@ -426,7 +428,7 @@ private:
int dfinv_matrix(int *ir);
// invert the matrix. See CERNLIB DFINV.
double *m;
std::vector<double,Alloc<double,25> > m;
int nrow, ncol;
int size;
};
......
Matrix/Matrix/Matrix.icc
View file @ 72622075
// -*- C++ -*-
// $Id: Matrix.icc,v 1.2 2003/07/18 05:31:48 garren Exp $
// $Id: Matrix.icc,v 1.2.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -56,7 +56,7 @@ inline double & HepMatrix::operator()(int row, int col)
if(row<1 || row>num_row() || col<1 || col>num_col())
error("Range error in HepMatrix::operator()");
#endif
return *(m+(row-1)*ncol+col-1);
return *(m.begin()+(row-1)*ncol+col-1);
}
inline const double & HepMatrix::operator()(int row, int col) const
......@@ -65,7 +65,7 @@ inline const double & HepMatrix::operator()(int row, int col) const
if(row<1 || row>num_row() || col<1 || col>num_col())
error("Range error in HepMatrix::operator()");
#endif
return *(m+(row-1)*ncol+col-1);
return *(m.begin()+(row-1)*ncol+col-1);
}
inline HepMatrix::HepMatrix_row HepMatrix::operator[] (int r)
......@@ -95,7 +95,7 @@ inline double &HepMatrix::HepMatrix_row::operator[](int c) {
if (_r<0 || _r>=_a.num_row() || c<0 || c>=_a.num_col())
HepGenMatrix::error("Range error in HepMatrix::operator[][]");
#endif
return *(_a.m+_r*_a.ncol+c);
return *(_a.m.begin()+_r*_a.ncol+c);
}
inline const double &HepMatrix::HepMatrix_row_const::operator[](int c) const
......@@ -104,7 +104,7 @@ inline const double &HepMatrix::HepMatrix_row_const::operator[](int c) const
if (_r<0 || _r>=_a.num_row() || c<0 || c>=_a.num_col())
HepGenMatrix::error("Range error in HepMatrix::operator[][]");
#endif
return *(_a.m+_r*_a.ncol+c);
return *(_a.m.begin()+_r*_a.ncol+c);
}
inline HepMatrix::HepMatrix_row::HepMatrix_row(HepMatrix&a,int r)
......@@ -122,12 +122,14 @@ inline HepMatrix::HepMatrix_row_const::HepMatrix_row_const
// This function swaps two Matrices without doing a full copy.
inline void swap(HepMatrix &m1,HepMatrix &m2) {
HepGenMatrix::swap(m1.m,m2.m);
/*** commented
HepGenMatrix::swap(m1.nrow,m2.nrow);
HepGenMatrix::swap(m1.ncol,m2.ncol);
HepGenMatrix::swap(m1.size,m2.size);
*/
}
inline HepMatrix HepMatrix::inverse(int &ierr) const
/*-ap inline */ HepMatrix HepMatrix::inverse(int &ierr) const
#ifdef HEP_GNU_OPTIMIZED_RETURN
return mTmp(*this);
{
......
Matrix/Matrix/SymMatrix.h
View file @ 72622075
// -*- C++ -*-
// CLASSDOC OFF
// $Id: SymMatrix.h,v 1.3 2003/10/23 21:29:50 garren Exp $
// $Id: SymMatrix.h,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
// CLASSDOC ON
//
......@@ -109,6 +109,8 @@
#pragma interface
#endif
#include <vector>
#include "CLHEP/Matrix/defs.h"
#include "CLHEP/Matrix/GenMatrix.h"
......@@ -298,7 +300,7 @@ private:
friend HepSymMatrix vT_times_v(const HepVector &v);
// Returns v * v.T();
double *m;
std::vector<double,Alloc<double,25> > m;
int nrow;
int size; // total number of elements
......
Matrix/Matrix/SymMatrix.icc
View file @ 72622075
// -*- C++ -*-
// $Id: SymMatrix.icc,v 1.2 2003/07/18 05:31:48 garren Exp $
// $Id: SymMatrix.icc,v 1.2.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -55,7 +55,7 @@ inline double & HepSymMatrix::fast(int row,int col)
if(row<1||row>num_row() || col<1||col>num_col())
error("Range error in HepSymMatrix::fast()");
#endif
return *(m+(row*(row-1))/2+(col-1));
return *(m.begin()+(row*(row-1))/2+(col-1));
}
inline const double & HepSymMatrix::fast(int row,int col) const
{
......@@ -63,7 +63,7 @@ inline const double & HepSymMatrix::fast(int row,int col) const
if(row<1||row>num_row() || col<1||col>num_col())
error("Range error in HepSymMatrix::fast()");
#endif
return *(m+(row*(row-1))/2+(col-1));
return *(m.begin()+(row*(row-1))/2+(col-1));
}
inline double & HepSymMatrix::operator()(int row, int col)
......@@ -105,9 +105,9 @@ inline double &HepSymMatrix::HepSymMatrix_row::operator[](int c)
error("Range error in HepSymMatrix::operator[][]");
#endif
if (_r >= c ) {
return *(_a.m + (_r+1)*_r/2 + c);
return *(_a.m.begin() + (_r+1)*_r/2 + c);
} else {
return *(_a.m + (c+1)*c/2 + _r);
return *(_a.m.begin() + (c+1)*c/2 + _r);
}
}
......@@ -119,9 +119,9 @@ HepSymMatrix::HepSymMatrix_row_const::operator[](int c) const
error("Range error in HepSymMatrix::operator[][]");
#endif
if (_r >= c ) {
return *(_a.m + (_r+1)*_r/2 + c);
return *(_a.m.begin() + (_r+1)*_r/2 + c);
} else {
return *(_a.m + (c+1)*c/2 + _r);
return *(_a.m.begin() + (c+1)*c/2 + _r);
}
}
......
Matrix/Matrix/Vector.h
View file @ 72622075
// -*- C++ -*-
// CLASSDOC OFF
// $Id: Vector.h,v 1.3 2003/10/23 21:29:50 garren Exp $
// $Id: Vector.h,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
// CLASSDOC ON
//
......@@ -196,7 +196,7 @@ private:
friend HepSymMatrix vT_times_v(const HepVector &v);
friend HepVector qr_solve(HepMatrix *, const HepVector &);
double *m;
std::vector<double,Alloc<double,25> > m;
int nrow;
};
......
Matrix/Matrix/Vector.icc
View file @ 72622075
// -*- C++ -*-
// $Id: Vector.icc,v 1.3 2003/07/18 05:31:48 garren Exp $
// $Id: Vector.icc,v 1.3.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -64,7 +64,7 @@ inline double & HepVector::operator()(int row)
error("Range error in HepVector::operator()");
#endif
return *(m+row-1);
return *(m.begin()+row-1);
}
inline const double & HepVector::operator()(int row) const
{
......@@ -73,7 +73,7 @@ inline const double & HepVector::operator()(int row) const
error("Range error in HepVector::operator()");
#endif
return *(m+row-1);
return *(m.begin()+row-1);
}
inline double & HepVector::operator[](int row)
{
......@@ -82,7 +82,7 @@ inline double & HepVector::operator[](int row)
error("Range error in HepVector::operator[]");
#endif
return *(m+row);
return *(m.begin()+row);
}
inline const double & HepVector::operator[](int row) const
{
......@@ -91,7 +91,7 @@ inline const double & HepVector::operator[](int row) const
error("Range error in HepVector::operator[]");
#endif
return *(m+row);
return *(m.begin()+row);
}
#ifdef MATRIX_BOUND_CHECK
......@@ -104,7 +104,7 @@ inline double & HepVector::operator()(int row, int)
{
#endif
return *(m+(row-1));
return *(m.begin()+(row-1));
}
#ifdef MATRIX_BOUND_CHECK
......@@ -117,7 +117,7 @@ inline const double & HepVector::operator()(int row, int) const
{
#endif
return *(m+(row-1));
return *(m.begin()+(row-1));
}
} // namespace CLHEP
......
Matrix/src/DiagMatrix.cc
View file @ 72622075
// -*- C++ -*-
// $Id: DiagMatrix.cc,v 1.4 2003/08/13 20:00:12 garren Exp $
// $Id: DiagMatrix.cc,v 1.4.2.1 2004/08/25 18:37:41 pfeiffer Exp $
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
......@@ -59,21 +59,21 @@ namespace CLHEP {
// Simple operation for all elements
#define SIMPLE_UOP(OPER) \
register double *a=m; \
register double *e=m+num_size(); \
register HepMatrix::mIter a=m.begin(); \
register HepMatrix::mIter e=m.begin()+num_size(); \
for(;a<e; a++) (*a) OPER t;
#define SIMPLE_BOP(OPER) \
register double *a=m; \
register double *b=m2.m; \
register double *e=m+num_size(); \
register HepMatrix::mIter a=m.begin(); \
register HepMatrix::mcIter b=m2.m.begin(); \
register HepMatrix::mIter e=m.begin()+num_size(); \
for(;a<e; a++, b++) (*a) OPER (*b);
#define SIMPLE_TOP(OPER) \
register double *a=m1.m; \
register double *b=m2.m; \
register double *t=mret.m; \
register double *e=m1.m+m1.nrow; \
register HepMatrix::mcIter a=m1.m.begin(); \
register HepMatrix::mcIter b=m2.m.begin(); \
register HepMatrix::mIter t=mret.m.begin(); \
register HepMatrix::mcIter e=m1.m.begin()+m1.nrow; \
for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
#define CHK_DIM_2(r1,r2,c1,c2,fun) \
......@@ -100,26 +100,23 @@ const double HepDiagMatrix::zero = 0;
// Constructors. (Default constructors are inlined and in .icc file)
HepDiagMatrix::HepDiagMatrix(int p)
: nrow(p)
: m(std::vector<double,Alloc<double,25> >(p)), nrow(p)
{
m = new_m(nrow);
}
HepDiagMatrix::HepDiagMatrix(int p, int init)
: nrow(p)
{
m = new_m(nrow);
: m(std::vector<double,Alloc<double,25> >(p)), nrow(p)
{
switch(init)
{
case 0:
memset(m, 0, nrow * sizeof(double));
m.assign(nrow,0);
break;
case 1:
{
double *a=m;
double *b=m + p;
HepMatrix::mIter a=m.begin();
HepMatrix::mIter b=m.begin() + p;
for( ; a<b; a++) *a = 1.0;
break;
}
......@@ -129,25 +126,22 @@ HepDiagMatrix::HepDiagMatrix(int p, int init)
}
HepDiagMatrix::HepDiagMatrix(int p, HepRandom &r)
: nrow(p)
: m(std::vector<double,Alloc<double,25> >(p)), nrow(p)
{
m = new_m(nrow);
double *a = m;
double *b = m + num_size();
HepMatrix::mIter a = m.begin();
HepMatrix::mIter b = m.begin() + num_size();
for(;a<b;a++) *a = r();
}
//
// Destructor
//
HepDiagMatrix::~HepDiagMatrix() {
delete_m(nrow, m);
}
HepDiagMatrix::HepDiagMatrix(const HepDiagMatrix &m1)
: nrow(m1.nrow)
: m(std::vector<double,Alloc<double,25> >(m1.nrow)), nrow(m1.nrow)
{
m = new_m(nrow);
memcpy(m, m1.m, nrow*sizeof(double));
m = m1.m;
}
//
......@@ -166,9 +160,9 @@ return mret(max_row-min_row+1);
#endif
if(max_row > num_row())
error("HepDiagMatrix::sub: Index out of range");
double *a = mret.m;
double *b = m + min_row - 1;
double *e = mret.m + mret.num_row();
HepMatrix::mIter a = mret.m.begin();
HepMatrix::mcIter b = m.begin() + min_row - 1;
HepMatrix::mIter e = mret.m.begin() + mret.num_row();
for(;a<e;) *(a++) = *(b++);
return mret;
}
......@@ -179,9 +173,9 @@ HepDiagMatrix HepDiagMatrix::sub(int min_row, int max_row)
HepDiagMatrix mret(max_row-min_row+1);
if(max_row > num_row())
error("HepDiagMatrix::sub: Index out of range");
double *a = mret.m;
double *b = m + min_row - 1;
double *e = mret.m + mret.num_row();
HepMatrix::mIter a = mret.m.begin();
HepMatrix::mIter b = m.begin() + min_row - 1;
HepMatrix::mIter e = mret.m.begin() + mret.num_row();
for(;a<e;) *(a++) = *(b++);
return mret;
}
......@@ -191,9 +185,9 @@ void HepDiagMatrix::sub(int row,const HepDiagMatrix &m1)
{
if(row <1 || row+m1.num_row()-1 > num_row() )
error("HepDiagMatrix::sub: Index out of range");
double *a = m1.m;
double *b = m + row - 1;
double *e = m1.m + m1.num_row();
HepMatrix::mcIter a = m1.m.begin();
HepMatrix::mIter b = m.begin() + row - 1;
HepMatrix::mcIter e = m1.m.begin() + m1.num_row();
for(;a<e;) *(b++) = *(a++);
}
......@@ -224,9 +218,9 @@ HepDiagMatrix HepDiagMatrix::operator- () const
{
HepDiagMatrix m2(nrow);
#endif
register double *a=m;
register double *b=m2.m;
register double *e=m+num_size();
register HepMatrix::mcIter a=m.begin();
register HepMatrix::mIter b=m2.m.begin();
register HepMatrix::mcIter e=m.begin()+num_size();
for(;a<e; a++, b++) (*b) = -(*a);
return m2;
}
......@@ -421,10 +415,10 @@ HepMatrix operator*(const HepMatrix &m1,const HepDiagMatrix &m2)
HepMatrix mret(m1.num_row(),m2.num_col());
#endif
CHK_DIM_1(m1.num_col(),m2.num_row(),*);
double *mit1=m1.m;
double *mir=mret.m;
HepMatrix::mcIter mit1=m1.m.begin();
HepMatrix::mIter mir=mret.m.begin();
for(int irow=1;irow<=m1.num_row();irow++) {
double *mcc = m2.m;
HepMatrix::mcIter mcc = m2.m.begin();
for(int icol=1;icol<=m1.num_col();icol++) {
*(mir++) = *(mit1++) * (*(mcc++));
}
......@@ -441,9 +435,9 @@ HepMatrix operator*(const HepDiagMatrix &m1,const HepMatrix &m2)
HepMatrix mret(m1.num_row(),m2.num_col());
#endif
CHK_DIM_1(m1.num_col(),m2.num_row(),*);
double *mit1=m2.m;
double *mir=mret.m;
double *mrr = m1.m;
HepMatrix::mcIter mit1=m2.m.begin();
HepMatrix::mIter mir=mret.m.begin();
HepMatrix::mcIter mrr = m1.m.begin();
for(int irow=1;irow<=m2.num_row();irow++) {
for(int icol=1;icol<=m2.num_col();icol++) {
*(mir++) = *(mit1++) * (*mrr);
......@@ -462,10 +456,10 @@ HepDiagMatrix operator*(const HepDiagMatrix &m1,const HepDiagMatrix &m2)
HepDiagMatrix mret(m1.num_row());
#endif
CHK_DIM_1(m1.num_col(),m2.num_row(),*);
double *a = mret.m;
double *b = m1.m;
double *c = m2.m;
double *e = mret.m + m1.num_col();
HepMatrix::mIter a = mret.m.begin();
HepMatrix::mcIter b = m1.m.begin();
HepMatrix::mcIter c = m2.m.begin();
HepMatrix::mIter e = mret.m.begin() + m1.num_col();
for(;a<e;) *(a++) = *(b++) * (*(c++));
return mret;
}
......@@ -479,7 +473,8 @@ HepVector operator*(const HepDiagMatrix &m1,const HepVector &m2)
HepVector mret(m1.num_row());
#endif
CHK_DIM_1(m1.num_col(),m2.num_row(),*);
double *mir=mret.m, *mi1 = m1.m, *mi2 = m2.m;
HepGenMatrix::mIter mir=mret.m.begin();
HepGenMatrix::mcIter mi1 = m1.m.begin(), mi2 = m2.m.begin();
for(int icol=1;icol<=m1.num_col();icol++) {
*(mir++) = *(mi1++) * *(mi2++);
}
......@@ -494,8 +489,8 @@ HepMatrix & HepMatrix::operator+=(const HepDiagMatrix &m2)
{
CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
int n = num_row();
double *mrr = m;
double *mr = m2.m;
mIter mrr = m.begin();
HepMatrix::mcIter mr = m2.m.begin();
for(int r=1;r<=n;r++) {
*mrr += *(mr++);
mrr += (n+1);
......@@ -506,8 +501,8 @@ HepMatrix & HepMatrix::operator+=(const HepDiagMatrix &m2)
HepSymMatrix & HepSymMatrix::operator+=(const HepDiagMatrix &m2)
{
CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=);
register double *a=m;
register double *b=m2.m;
register HepMatrix::mIter a=m.begin();
register HepMatrix::mcIter b=m2.m.begin();
for(int i=1;i<=num_row();i++) {
*a += *(b++);
a += (i+1);
......@@ -526,8 +521,8 @@ HepMatrix & HepMatrix::operator-=(const HepDiagMatrix &m2)
{
CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=);
int n = num_row();
double *mrr = m;
double *mr = m2.m;
mIter mrr = m.begin();
HepMatrix::mcIter mr = m2.m.begin();
for(int r=1;r<=n;r++) {
*mrr -= *(mr++);
mrr += (n+1);
......@@ -538,8 +533,8 @@ HepMatrix & HepMatrix::operator-=(const HepDiagMatrix &m2)
HepSymMatrix & HepSymMatrix::operator-=(const HepDiagMatrix &m2