Skip to content

GitLab

Commit f27b20ea authored by Lynn Garren's avatar Lynn Garren
Browse files

fix post incremented iterators

parent e0a87672
Expand all Hide whitespace changes
Inline Side-by-side
Matrix/ChangeLog
View file @ f27b20ea
2008-07-16 Lynn Garren <garren@fnal.gov>
* src/SymMatrix.cc, Matrix.cc, DiagMatrix.cc, MatrixLinear.cc:
Iterators were set to a nonexistent memory location in many places.
Even though the iterators were not used, this operation is not allowed
and causes problems with the newer VC++ compilers.
In some cases, a more efficient rewrite was possible.
==============================
01.05.08 Release CLHEP-2.0.3.3
==============================
......
Matrix/src/DiagMatrix.cc
View file @ f27b20ea
......@@ -490,7 +490,7 @@ HepMatrix & HepMatrix::operator+=(const HepDiagMatrix &m2)
HepMatrix::mcIter mr = m2.m.begin();
for(int r=1;r<=n;r++) {
*mrr += *(mr++);
mrr += (n+1);
if(r<n) mrr += (n+1);
}
return (*this);
}
......@@ -502,7 +502,7 @@ HepSymMatrix & HepSymMatrix::operator+=(const HepDiagMatrix &m2)
register HepMatrix::mcIter b=m2.m.begin();
for(int i=1;i<=num_row();i++) {
*a += *(b++);
a += (i+1);
if(i<num_row()) a += (i+1);
}
return (*this);
}
......@@ -522,7 +522,7 @@ HepMatrix & HepMatrix::operator-=(const HepDiagMatrix &m2)
HepMatrix::mcIter mr = m2.m.begin();
for(int r=1;r<=n;r++) {
*mrr -= *(mr++);
mrr += (n+1);
if(r<n) mrr += (n+1);
}
return (*this);
}
......@@ -534,7 +534,7 @@ HepSymMatrix & HepSymMatrix::operator-=(const HepDiagMatrix &m2)
register HepMatrix::mcIter b=m2.m.begin();
for(int i=1;i<=num_row();i++) {
*a -= *(b++);
a += (i+1);
if(i<num_row()) a += (i+1);
}
return (*this);
}
......@@ -573,7 +573,7 @@ HepMatrix & HepMatrix::operator=(const HepDiagMatrix &m1)
HepMatrix::mcIter mr = m1.m.begin();
for(int r=1;r<=n;r++) {
*mrr = *(mr++);
mrr += (n+1);
if(r<n) mrr += (n+1);
}
return (*this);
}
......@@ -640,7 +640,7 @@ void HepDiagMatrix::assign (const HepMatrix &m1)
HepMatrix::mIter b = m.begin();
for(int r=1;r<=nrow;r++) {
*(b++) = *a;
a += (nrow+1);
if(r<nrow) a += (nrow+1);
}
}
......@@ -655,7 +655,7 @@ void HepDiagMatrix::assign(const HepSymMatrix &m1)
HepMatrix::mIter b = m.begin();
for(int r=1;r<=nrow;r++) {
*(b++) = *a;
a += (r+1);
if(r<nrow) a += (r+1);
}
}
......
Matrix/src/Matrix.cc
View file @ f27b20ea
......@@ -112,8 +112,7 @@ HepMatrix::HepMatrix(int p,int q,int init)
{
if ( ncol == nrow ) {
mIter a = m.begin();
mIter b = m.end();
for( ; a<b; a+=(ncol+1)) *a = 1.0;
for( int step=0; step < size; step+=(ncol+1) ) *(a+step) = 1.0;
} else {
error("Invalid dimension in HepMatrix(int,int,1).");
}
......@@ -152,23 +151,18 @@ HepMatrix::HepMatrix(const HepSymMatrix &m1)
{
size = nrow * ncol;
int n = ncol;
mcIter sjk = m1.m.begin();
mIter m1j = m.begin();
mIter mj = m.begin();
// j >= k
for(int j=1;j<=nrow;j++) {
mIter mjk = mj;
mIter mkj = m1j;
for(int k=1;k<=j;k++) {
*(mjk++) = *sjk;
if(j!=k) *mkj = *sjk;
sjk++;
mkj += n;
}
mj += n;
m1j++;
}
for(int j=0; j!=nrow; ++j) {
for(int k=0; k<=j; ++k) {
m[j*ncol+k] = *sjk;
// we could copy the diagonal element twice or check
// doing the check may be a tiny bit faster,
// so we choose that option for now
if(k!=j) m[k*nrow+j] = *sjk;
++sjk;
}
}
}
HepMatrix::HepMatrix(const HepDiagMatrix &m1)
......@@ -177,11 +171,11 @@ HepMatrix::HepMatrix(const HepDiagMatrix &m1)
size = nrow * ncol;
int n = num_row();
mIter mrr = m.begin();
mIter mrr;
mcIter mr = m1.m.begin();
for(int r=1;r<=n;r++) {
for(int r=0;r<n;r++) {
mrr = m.begin()+(n+1)*r;
*mrr = *(mr++);
mrr += (n+1);
}
}
......@@ -214,13 +208,13 @@ return mret(max_row-min_row+1,max_col-min_col+1);
mIter a = mret.m.begin();
int nc = num_col();
mcIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
for(int irow=1; irow<=mret.num_row(); irow++) {
int rowsize = mret.num_row();
for(int irow=1; irow<=rowsize; ++irow) {
mcIter brc = b1;
for(int icol=1; icol<=mret.num_col(); icol++) {
for(int icol=0; icol<mret.num_col(); ++icol) {
*(a++) = *(brc++);
}
b1 += nc;
if(irow<rowsize) b1 += nc;
}
return mret;
}
......@@ -233,13 +227,13 @@ void HepMatrix::sub(int row,int col,const HepMatrix &m1)
mcIter a = m1.m.begin();
int nc = num_col();
mIter b1 = m.begin() + (row - 1) * nc + col - 1;
for(int irow=1; irow<=m1.num_row(); irow++) {
int rowsize = m1.num_row();
for(int irow=1; irow<=rowsize; ++irow) {
mIter brc = b1;
for(int icol=1; icol<=m1.num_col(); icol++) {
for(int icol=0; icol<m1.num_col(); ++icol) {
*(brc++) = *(a++);
}
b1 += nc;
if(irow<rowsize) b1 += nc;
}
}
......@@ -469,14 +463,14 @@ return mret(ncol,nrow);
{
HepMatrix mret(ncol,nrow);
#endif
register mcIter pl = m.end();
register mcIter pme = m.begin();
register mIter pt = mret.m.begin();
register mIter ptl = mret.m.end();
for (; pme < pl; pme++, pt+=nrow)
{
if (pt >= ptl) pt -= (size-1);
(*pt) = (*pme);
for( int nr=0; nr<nrow; ++nr) {
for( int nc=0; nc<ncol; ++nc) {
pt = mret.m.begin() + nr + nrow*nc;
(*pt) = (*pme);
++pme;
}
}
return mret;
}
......@@ -515,12 +509,13 @@ int HepMatrix::dfinv_matrix(int *ir) {
*m21 = -(*m22) * (*m11) * (*m21);
*m12 = -(*m12);
if (n>2) {
mIter mi = m.begin() + 2 * n;
mIter mii= m.begin() + 2 * n + 2;
mIter mimim = m.begin() + n + 1;
mIter mimim = m11 + n + 1;
for (int i=3;i<=n;i++) {
// calculate these to avoid pointing off the end of the storage array
mIter mi = m11 + (i-1) * n;
mIter mii= m11 + (i-1) * n + i - 1;
int im2 = i - 2;
mIter mj = m.begin();
mIter mj = m11;
mIter mji = mj + i - 1;
mIter mij = mi;
for (int j=1;j<=im2;j++) {
......@@ -535,50 +530,50 @@ int HepMatrix::dfinv_matrix(int *ir) {
s32 += (*(mjkp++)) * (*mkpi);
mkj += n;
mkpi += n;
}
} // for k
*mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
*mji = -s32;
mj += n;
mji += n;
mij++;
}
} // for j
*(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
*(mimim+1) = -(*(mimim+1));
mi += n;
mimim += (n+1);
mii += (n+1);
}
}
mIter mi = m.begin();
mIter mii = m.begin();
} // for i
} // n>2
mIter mi = m11;
mIter mii = m11;
for (int i=1;i<n;i++) {
int ni = n - i;
mIter mij = mi;
int j;
for (j=1; j<=i;j++) {
s33 = *mij;
register mIter mikj = mi + n + j - 1;
// change initial definition of mikj to avoid pointing off the end of the storage array
register mIter mikj = mi + j - 1;
register mIter miik = mii + 1;
mIter min_end = mi + n;
for (;miik<min_end;) {
s33 += (*mikj) * (*(miik++));
// iterate by n as we enter the loop to avoid pointing off the end of the storage array
mikj += n;
s33 += (*mikj) * (*(miik++));
}
*(mij++) = s33;
}
for (j=1;j<=ni;j++) {
s34 = 0.0;
mIter miik = mii + j;
mIter mikij = mii + j * n + j;
for (int k=j;k<=ni;k++) {
// calculate mikij here to avoid pointing off the end of the storage array
mIter mikij = mii + k * n + j;
s34 += *mikij * (*(miik++));
mikij += n;
}
*(mii+j) = s34;
}
mi += n;
mii += (n+1);
}
} // for i
int nxch = ir[n];
if (nxch==0) return 0;
for (int mm=1;mm<=nxch;mm++) {
......@@ -586,18 +581,17 @@ int HepMatrix::dfinv_matrix(int *ir) {
int ij = ir[k];
int i = ij >> 12;
int j = ij%4096;
mIter mki = m.begin() + i - 1;
mIter mkj = m.begin() + j - 1;
for (k=1; k<=n;k++) {
// avoid setting the iterator beyond the end of the storage vector
mIter mki = m11 + (k-1)*n + i - 1;
mIter mkj = m11 + (k-1)*n + j - 1;
// 2/24/05 David Sachs fix of improper swap bug that was present
// for many years:
double ti = *mki; // 2/24/05
*mki = *mkj;
*mkj = ti; // 2/24/05
mki += n;
mkj += n;
}
}
} // for mm
return 0;
}
......@@ -631,15 +625,14 @@ int HepMatrix::dfact_matrix(double &det, int *ir) {
int k = j;
p = (fabs(*mjj));
if (j!=n) {
mIter mij = mj + n + j - 1;
for (int i=j+1;i<=n;i++) {
q = (fabs(*(mij)));
// replace mij with calculation of position
for (int i=j+1;i<n;i++) {
q = (fabs(*(mj + n*(i-j) + j - 1)));
if (q > p) {
k = i;
p = q;
}
mij += n;
}
} // for i
if (k==j) {
if (p <= epsilon) {
det = 0;
......@@ -649,7 +642,7 @@ int HepMatrix::dfact_matrix(double &det, int *ir) {
}
det = -det; // in this case the sign of the determinant
// must not change. So I change it twice.
}
} // k==j
mIter mjl = mj;
mIter mkl = m.begin() + (k-1)*n;
for (int l=1;l<=n;l++) {
......@@ -659,14 +652,14 @@ int HepMatrix::dfact_matrix(double &det, int *ir) {
}
nxch = nxch + 1; // this makes the determinant change its sign
ir[nxch] = (((j)<<12)+(k));
} else {
} else { // j!=n
if (p <= epsilon) {
det = 0.0;
ifail = imposs;
jfail = jrange;
return ifail;
}
}
} // j!=n
det *= *mjj;
*mjj = 1.0 / *mjj;
t = (fabs(det));
......@@ -677,11 +670,12 @@ int HepMatrix::dfact_matrix(double &det, int *ir) {
det = 1.0;
if (jfail==jrange) jfail = jover;
}
// calculate mk and mkjp so we don't point off the end of the vector
if (j!=n) {
mIter mk = mj + n;
mIter mkjp = mk + j;
mIter mjk = mj + j;
for (k=j+1;k<=n;k++) {
mIter mk = mj + n*(k-j);
mIter mkjp = mk + j;
s11 = - (*mjk);
s12 = - (*mkjp);
if (j!=1) {
......@@ -694,17 +688,18 @@ int HepMatrix::dfact_matrix(double &det, int *ir) {
s12 += (*mijp) * (*(mki++));
mik += n;
mijp += n;
}
}
} // for i
} // j!=1
*(mjk++) = -s11 * (*mjj);
*(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
mk += n;
mkjp += n;
}
} // for k
} // j!=n
// avoid setting the iterator beyond the end of the vector
if(j!=n) {
mj += n;
mjj += (n+1);
}
mj += n;
mjj += (n+1);
}
} // for j
if (nxch%2==1) det = -det;
if (jfail !=jrange) det = 0.0;
ir[n] = nxch;
......
Matrix/src/MatrixLinear.cc
View file @ f27b20ea
......@@ -102,8 +102,10 @@ void back_solve(const HepMatrix &R, HepVector *b)
(*br)-=(*(Rrc++))*(*(bc++));
}
(*br)/=(*Rrr);
br--;
Rrr -= n+1;
if(r>1) {
br--;
Rrr -= n+1;
}
}
}
......@@ -128,11 +130,13 @@ void back_solve(const HepMatrix &R, HepMatrix *b)
HepMatrix::mcIter Rrc = Rrr+1;
for (int c=r+1;c<=b->num_row();c++) {
(*bri)-= (*(Rrc++)) * (*bci);
bci += nc;
if(c<b->num_row()) bci += nc;
}
(*bri)/= (*Rrr);
Rrr -= (n+1);
bri -= nc;
if(r>1) {
Rrr -= (n+1);
bri -= nc;
}
}
bbi++;
}
......@@ -154,8 +158,10 @@ void col_givens(HepMatrix *A, double c,double s,
for (int j=row_min;j<=row_max;j++) {
double tau1=(*Ajk1); double tau2=(*Ajk2);
(*Ajk1)=c*tau1-s*tau2;(*Ajk2)=s*tau1+c*tau2;
Ajk1 += n;
Ajk2 += n;
if(j<row_max) {
Ajk1 += n;
Ajk2 += n;
}
}
}
......@@ -194,7 +200,7 @@ void col_house(HepMatrix *a,const HepMatrix &v,double vnormsq,
vp += nv;
}
wptr++;
acrb += n;
if(c<a->num_col()) acrb += n;
}
w*=beta;
......@@ -210,7 +216,7 @@ void col_house(HepMatrix *a,const HepMatrix &v,double vnormsq,
vp += nv;
}
wptr++;
arcb += n;
if(r<a->num_row()) arcb += n;
}
}
......@@ -232,7 +238,7 @@ double condition(const HepSymMatrix &m)
for (int i=2; i<=n; i++) {
if (max<fabs(*mii)) max=fabs(*mii);
if (min>fabs(*mii)) min=fabs(*mii);
mii += i+1;
if(i<n) mii += i+1;
}
return max/min;
}
......@@ -282,10 +288,10 @@ void diag_step(HepSymMatrix *t,int begin,int end)
(*(tkp2k+1))=bq*c;
x=(*tkp1k);
z=(*tkp2k);
tkk += k+1;
tkp1k += k+2;
}
tkk += k+1;
tkp1k += k+2;
tkp2k += k+3;
if(k<end-2) tkp2k += k+3;
}
}
......@@ -324,10 +330,10 @@ void diag_step(HepSymMatrix *t,HepMatrix *u,int begin,int end)
(*(tkp2k+1))=bq*c;
x=(*tkp1k);
z=(*(tkp2k));
tkk += k+1;
tkp1k += k+2;
}
tkk += k+1;
tkp1k += k+2;
tkp2k += k+3;
if(k<end-2) tkp2k += k+3;
}
}
......@@ -352,8 +358,10 @@ HepMatrix diagonalize(HepSymMatrix *s)
tolerance*(fabs(*sii)+fabs(*(sip1i+1)))) {
(*sip1i)=0;
}
sii += i+1;
sip1i += i+2;
if(i<end-1) {
sii += i+1;
sip1i += i+2;
}
}
while(begin<end && s->fast(begin+1,begin) ==0) begin++;
while(end>begin && s->fast(end,end-1) ==0) end--;
......@@ -423,7 +431,7 @@ void house_with_update(HepMatrix *a,int row,int col)
int r;
for (r=row;r<=a->num_row();r++) {
(*(vp++))=(*arc);
arc += n;
if(r<a->num_row()) arc += n;
}
double normsq=v.normsq();
double norm=sqrt(normsq);
......@@ -431,12 +439,12 @@ void house_with_update(HepMatrix *a,int row,int col)
v(1)+=sign((*a)(row,col))*norm;
normsq+=v(1)*v(1);
(*a)(row,col)=-sign((*a)(row,col))*norm;
arc = a->m.begin() + row * n + (col-1);
for (r=row+1;r<=a->num_row();r++) {
(*arc)=0;
arc += n;
}
if (row<a->num_row()) {
arc = a->m.begin() + row * n + (col-1);
for (r=row+1;r<=a->num_row();r++) {
(*arc)=0;
if(r<a->num_row()) arc += n;
}
row_house(a,v,normsq,row,col+1);
}
}
......@@ -452,8 +460,10 @@ void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col)
for (r=row;r<=a->num_row();r++) {
(*vrc)=(*arc);
normsq+=(*vrc)*(*vrc);
vrc += nv;
arc += na;
if(r<a->num_row()) {
vrc += nv;
arc += na;
}
}
double norm=sqrt(normsq);
vrc = v->m.begin() + (row-1) * nv + (col-1);
......@@ -461,13 +471,14 @@ void house_with_update(HepMatrix *a,HepMatrix *v,int row,int col)
(*vrc)+=sign((*a)(row,col))*norm;
normsq+=(*vrc)*(*vrc);
(*a)(row,col)=-sign((*a)(row,col))*norm;
arc = a->m.begin() + row * na + (col-1);
for (r=row+1;r<=a->num_row();r++) {
(*arc)=0;
arc += na;
}
if (row<a->num_row())
if (row<a->num_row()) {
arc = a->m.begin() + row * na + (col-1);
for (r=row+1;r<=a->num_row();r++) {
(*arc)=0;
if(r<a->num_row()) arc += na;
}
row_house(a,*v,normsq,row,col+1,row,col);
}
}
/* -----------------------------------------------------------------------
house_with_update2 Version: 1.00 Date Last Changed: 1/28/93
......@@ -487,8 +498,10 @@ void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col)
{
(*vrc)=(*arc);
normsq+=(*vrc)*(*vrc);
arc += r;
vrc += nv;
if(r<a->num_row()) {
arc += r;
vrc += nv;
}
}
double norm=sqrt(normsq);
vrc = v->m.begin() + (row-1) * nv + (col-1);
......@@ -498,7 +511,7 @@ void house_with_update2(HepSymMatrix *a,HepMatrix *v,int row,int col)
arc += row;
for (r=row+1;r<=a->num_row();r++) {
(*arc)=0;
arc += r;
if(r<a->num_row()) arc += r;
}
}
......@@ -642,7 +655,7 @@ void row_house(HepMatrix *a,const HepVector &v,double vnormsq,
HepMatrix::mIter arc = arcb;
for (int r=row;r<=a->num_row();r++) {
(*wptr)+=(*arc)*(*(vp++));
arc += na;
if(r<a->num_row()) arc += na;
}
wptr++;
arcb++;
......@@ -660,7 +673,7 @@ void row_house(HepMatrix *a,const HepVector &v,double vnormsq,
(*(arc++))+=(*vp)*(*(wptr++));
}
vp++;
arcb += na;
if(r<a->num_row()) arcb += na;
}
}
......@@ -681,8 +694,10 @@ void row_house(HepMatrix *a,const HepMatrix &v,double vnormsq,
HepMatrix::mcIter vpc = vpcb;
for (int r=row;r<=a->num_row();r++) {
(*wptr)+=(*arc)*(*vpc);
arc += na;
vpc += nv;
if(r<a->num_row()) {
arc += na;
vpc += nv;
}
}
wptr++;
arcb++;
......@@ -697,8 +712,10 @@ void row_house(HepMatrix *a,const HepMatrix &v,double vnormsq,
for (c=col;c<=a->num_col();c++) {
(*(arc++))+=(*vpc)*(*(wptr++));
}
arcb += na;
vpc += nv;
if(r<a->num_row()) {
arcb += na;
vpc += nv;
}