From 0ae59302dcdd90c48ea302b0b21128562753440c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 6 May 2002 07:32:20 +0000 Subject: [PATCH] Qualify access of member variables of template dependent base classes by this->, as requested by the aCC compiler. git-svn-id: https://svn.dealii.org/trunk@5804 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/point.h | 22 +- deal.II/lac/include/lac/full_matrix.h | 4 +- .../lac/include/lac/full_matrix.templates.h | 327 ++++++++++-------- .../lac/include/lac/sparse_ilu.templates.h | 55 +-- 4 files changed, 218 insertions(+), 190 deletions(-) diff --git a/deal.II/base/include/base/point.h b/deal.II/base/include/base/point.h index 6d27eba533..20c3c2492d 100644 --- a/deal.II/base/include/base/point.h +++ b/deal.II/base/include/base/point.h @@ -200,7 +200,7 @@ inline Point::Point (const double x) { Assert (dim==1, ExcInvalidConstructorCalled()); - values[0] = x; + this->values[0] = x; }; @@ -210,8 +210,8 @@ inline Point::Point (const double x, const double y) { Assert (dim==2, ExcInvalidConstructorCalled()); - values[0] = x; - values[1] = y; + this->values[0] = x; + this->values[1] = y; }; @@ -221,9 +221,9 @@ inline Point::Point (const double x, const double y, const double z) { Assert (dim==3, ExcInvalidConstructorCalled()); - values[0] = x; - values[1] = y; - values[2] = z; + this->values[0] = x; + this->values[1] = y; + this->values[2] = z; }; @@ -233,7 +233,7 @@ inline double Point::operator () (const unsigned int index) const { Assert (indexvalues[index]; }; @@ -243,7 +243,7 @@ inline double & Point::operator () (const unsigned int index) { Assert (indexvalues[index]; }; @@ -272,7 +272,7 @@ Point Point::operator - () const { Point result; for (unsigned int i=0; ivalues[i] = -this->values[i]; return result; }; @@ -302,7 +302,7 @@ double Point::square () const { double q=0; for (unsigned int i=0; ivalues[i] * this->values[i]; return q; }; @@ -314,7 +314,7 @@ double Point::distance (const Point &p) const double sum=0; for (unsigned int i=0; ivalues[i]-p(i); sum += diff*diff; } diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index 5c841f3b70..01b33d8eb0 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -735,7 +735,7 @@ template inline unsigned int FullMatrix::m() const { - return n_rows(); + return this->n_rows(); }; @@ -744,7 +744,7 @@ template inline unsigned int FullMatrix::n() const { - return n_cols(); + return this->n_cols(); }; diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index bba5d9850c..8734282943 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -75,10 +75,10 @@ template bool FullMatrix::all_zero () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); - const number *p = data(), - *e = data() + n()*m(); + const number *p = this->data(), + *e = this->data() + n()*m(); while (p!=e) if (*p++ != 0.0) return false; @@ -109,7 +109,7 @@ FullMatrix::vmult (Vector& dst, const Vector& src, const bool adding) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); @@ -119,11 +119,19 @@ FullMatrix::vmult (Vector& dst, number2 s; number2 s0,s1,s2; s = src(0); - s0 = s*data()[0]; s1 = s*data()[3]; s2 = s*data()[6]; + s0 = s*this->data()[0]; + s1 = s*this->data()[3]; + s2 = s*this->data()[6]; + s = src(1); - s0 += s*data()[1]; s1 += s*data()[4]; s2 += s*data()[7]; + s0 += s*this->data()[1]; + s1 += s*this->data()[4]; + s2 += s*this->data()[7]; + s = src(2); - s0 += s*data()[2]; s1 += s*data()[5]; s2 += s*data()[8]; + s0 += s*this->data()[2]; + s1 += s*this->data()[5]; + s2 += s*this->data()[8]; if (!adding) { @@ -143,13 +151,28 @@ FullMatrix::vmult (Vector& dst, number2 s; number2 s0,s1,s2,s3; s = src(0); - s0 = s*data()[0]; s1 = s*data()[4]; s2 = s*data()[8]; s3 = s*data()[12]; + s0 = s*this->data()[0]; + s1 = s*this->data()[4]; + s2 = s*this->data()[8]; + s3 = s*this->data()[12]; + s = src(1); - s0 += s*data()[1]; s1 += s*data()[5]; s2 += s*data()[9]; s3 += s*data()[13]; + s0 += s*this->data()[1]; + s1 += s*this->data()[5]; + s2 += s*this->data()[9]; + s3 += s*this->data()[13]; + s = src(2); - s0 += s*data()[2]; s1 += s*data()[6]; s2 += s*data()[10]; s3 += s*data()[14]; + s0 += s*this->data()[2]; + s1 += s*this->data()[6]; + s2 += s*this->data()[10]; + s3 += s*this->data()[14]; + s = src(3); - s0 += s*data()[3]; s1 += s*data()[7]; s2 += s*data()[11]; s3 += s*data()[15]; + s0 += s*this->data()[3]; + s1 += s*this->data()[7]; + s2 += s*this->data()[11]; + s3 += s*this->data()[15]; if (!adding) { @@ -171,29 +194,29 @@ FullMatrix::vmult (Vector& dst, number2 s; number2 s0,s1,s2,s3,s4,s5,s6,s7; s = src(0); - s0 = s*data()[0]; s1 = s*data()[8]; s2 = s*data()[16]; s3 = s*data()[24]; - s4 = s*data()[32]; s5 = s*data()[40]; s6 = s*data()[48]; s7 = s*data()[56]; + s0 = s*this->data()[0]; s1 = s*this->data()[8]; s2 = s*this->data()[16]; s3 = s*this->data()[24]; + s4 = s*this->data()[32]; s5 = s*this->data()[40]; s6 = s*this->data()[48]; s7 = s*this->data()[56]; s = src(1); - s0 += s*data()[1]; s1 += s*data()[9]; s2 += s*data()[17]; s3 += s*data()[25]; - s4 += s*data()[33]; s5 += s*data()[41]; s6 += s*data()[49]; s7 += s*data()[57]; + s0 += s*this->data()[1]; s1 += s*this->data()[9]; s2 += s*this->data()[17]; s3 += s*this->data()[25]; + s4 += s*this->data()[33]; s5 += s*this->data()[41]; s6 += s*this->data()[49]; s7 += s*this->data()[57]; s = src(2); - s0 += s*data()[2]; s1 += s*data()[10]; s2 += s*data()[18]; s3 += s*data()[26]; - s4 += s*data()[34]; s5 += s*data()[42]; s6 += s*data()[50]; s7 += s*data()[58]; + s0 += s*this->data()[2]; s1 += s*this->data()[10]; s2 += s*this->data()[18]; s3 += s*this->data()[26]; + s4 += s*this->data()[34]; s5 += s*this->data()[42]; s6 += s*this->data()[50]; s7 += s*this->data()[58]; s = src(3); - s0 += s*data()[3]; s1 += s*data()[11]; s2 += s*data()[19]; s3 += s*data()[27]; - s4 += s*data()[35]; s5 += s*data()[43]; s6 += s*data()[51]; s7 += s*data()[59]; + s0 += s*this->data()[3]; s1 += s*this->data()[11]; s2 += s*this->data()[19]; s3 += s*this->data()[27]; + s4 += s*this->data()[35]; s5 += s*this->data()[43]; s6 += s*this->data()[51]; s7 += s*this->data()[59]; s = src(4); - s0 += s*data()[4]; s1 += s*data()[12]; s2 += s*data()[20]; s3 += s*data()[28]; - s4 += s*data()[36]; s5 += s*data()[44]; s6 += s*data()[52]; s7 += s*data()[60]; + s0 += s*this->data()[4]; s1 += s*this->data()[12]; s2 += s*this->data()[20]; s3 += s*this->data()[28]; + s4 += s*this->data()[36]; s5 += s*this->data()[44]; s6 += s*this->data()[52]; s7 += s*this->data()[60]; s = src(5); - s0 += s*data()[5]; s1 += s*data()[13]; s2 += s*data()[21]; s3 += s*data()[29]; - s4 += s*data()[37]; s5 += s*data()[45]; s6 += s*data()[53]; s7 += s*data()[61]; + s0 += s*this->data()[5]; s1 += s*this->data()[13]; s2 += s*this->data()[21]; s3 += s*this->data()[29]; + s4 += s*this->data()[37]; s5 += s*this->data()[45]; s6 += s*this->data()[53]; s7 += s*this->data()[61]; s = src(6); - s0 += s*data()[6]; s1 += s*data()[14]; s2 += s*data()[22]; s3 += s*data()[30]; - s4 += s*data()[38]; s5 += s*data()[46]; s6 += s*data()[54]; s7 += s*data()[62]; + s0 += s*this->data()[6]; s1 += s*this->data()[14]; s2 += s*this->data()[22]; s3 += s*this->data()[30]; + s4 += s*this->data()[38]; s5 += s*this->data()[46]; s6 += s*this->data()[54]; s7 += s*this->data()[62]; s = src(7); - s0 += s*data()[7]; s1 += s*data()[15]; s2 += s*data()[23]; s3 += s*data()[31]; - s4 += s*data()[39]; s5 += s*data()[47]; s6 += s*data()[55]; s7 += s*data()[63]; + s0 += s*this->data()[7]; s1 += s*this->data()[15]; s2 += s*this->data()[23]; s3 += s*this->data()[31]; + s4 += s*this->data()[39]; s5 += s*this->data()[47]; s6 += s*this->data()[55]; s7 += s*this->data()[63]; if (!adding) { @@ -220,7 +243,7 @@ FullMatrix::vmult (Vector& dst, } else { - const number* e = data(); + const number* e = this->data(); const unsigned int size_m = m(), size_n = n(); if (!adding) @@ -254,7 +277,7 @@ void FullMatrix::Tvmult (Vector &dst, const Vector &src, const bool adding) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m())); @@ -291,7 +314,7 @@ double FullMatrix::residual (Vector& dst, const Vector& src, const Vector& right) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); @@ -318,7 +341,7 @@ template void FullMatrix::forward (Vector &dst, const Vector &src) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n())); @@ -341,7 +364,7 @@ template void FullMatrix::backward (Vector &dst, const Vector &src) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); unsigned int j; unsigned int nu = (m()::fill_permutation (const FullMatrix &src, /* void FullMatrix::fill (const number2* entries) */ /* { */ /* if (n_cols()*n_rows() != 0) */ -/* std::copy (entries, entries+n_rows()*n_cols(), data()); */ +/* std::copy (entries, entries+n_rows()*n_cols(), this->data()); */ /* } */ @@ -404,7 +427,7 @@ void FullMatrix::add_row (const unsigned int i, const number s, const unsigned int j) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); for (unsigned int k=0; k::add_row (const unsigned int i, const number t, const unsigned int k) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); const unsigned int size_m = m(); for (unsigned l=0; l void FullMatrix::add_col (const unsigned int i, const number s, const unsigned int j) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); for (unsigned int k=0; k::add_col (const unsigned int i, const number s, const unsigned int j, const number t, const unsigned int k) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); for (unsigned int l=0; l::add_col (const unsigned int i, const number s, template void FullMatrix::swap_row (const unsigned int i, const unsigned int j) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number s; for (unsigned int k=0; k::swap_row (const unsigned int i, const unsigned int j) template void FullMatrix::swap_col (const unsigned int i, const unsigned int j) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number s; for (unsigned int k=0; k::swap_col (const unsigned int i, const unsigned int j) template void FullMatrix::diagadd (const number src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (m() == n(), ExcDimensionMismatch(m(),n())); for (unsigned int i=0; i::mmult (FullMatrix &dst, const FullMatrix &src, const bool adding) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m())); Assert (dst.n() == src.n(), ExcDimensionMismatch(dst.n(), src.n())); Assert (dst.m() == m(), ExcDimensionMismatch(m(), dst.m())); @@ -525,7 +548,7 @@ void FullMatrix::Tmmult (FullMatrix &dst, const FullMatrix &src, const bool adding) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); Assert (n() == dst.m(), ExcDimensionMismatch(n(), dst.m())); Assert (src.n() == dst.n(), ExcDimensionMismatch(src.n(), dst.n())); @@ -556,14 +579,14 @@ template template number2 FullMatrix::matrix_norm_square (const Vector &v) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size())); Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); number2 sum = 0.; const unsigned int n_rows = m(); - const number *val_ptr = data(); + const number *val_ptr = this->data(); const number2 *v_ptr; for (unsigned int row=0; row number2 FullMatrix::matrix_scalar_product (const Vector &u, const Vector &v) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert(m() == u.size(), ExcDimensionMismatch(m(),v.size())); Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); @@ -594,7 +617,7 @@ number2 FullMatrix::matrix_scalar_product (const Vector &u, number2 sum = 0.; const unsigned int n_rows = m(); const unsigned int n_cols = n(); - const number *val_ptr = data(); + const number *val_ptr = this->data(); const number2 *v_ptr; for (unsigned int row=0; row::symmetrize () template number FullMatrix::l1_norm () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number sum=0, max=0; const unsigned int n_rows = m(), n_cols = n(); @@ -652,7 +675,7 @@ number FullMatrix::l1_norm () const template number FullMatrix::linfty_norm () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number sum=0, max=0; const unsigned int n_rows = m(), n_cols = n(); @@ -675,7 +698,7 @@ FullMatrix::print (std::ostream &s, const unsigned int w, const unsigned int p) const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); for (unsigned int i=0; i void FullMatrix::add (const number s,const FullMatrix& src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - number* val = const_cast (data()); + number* val = const_cast (this->data()); const number2* srcval = src.data(); if ((n()==3) && (m()==3)) @@ -816,12 +839,12 @@ template void FullMatrix::add_diag (const number s, const FullMatrix& src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - number* val = const_cast (data()); + number* val = const_cast (this->data()); const number2* srcval = src.data(); if ((n()==3) && (m()==3)) @@ -941,13 +964,13 @@ template void FullMatrix::Tadd (const number s, const FullMatrix& src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (m() == n(), ExcNotQuadratic()); Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - number* val = const_cast (data()); + number* val = const_cast (this->data()); const number2* srcval = src.data(); if ((n()==3) && (m()==3)) @@ -1073,9 +1096,9 @@ FullMatrix::operator == (const FullMatrix &M) const // empty, or of same size and with // same values, if they shall be // equal - bool result = (data()==0) && (M.data()==0); + bool result = (this->data()==0) && (M.data()==0); result = result || ((m()==M.m()) && (n()==M.n()) && - std::equal (data(), data()+m()*n(), + std::equal (this->data(), this->data()+m()*n(), M.data())); return result; @@ -1086,25 +1109,25 @@ template double FullMatrix::determinant () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (n_cols() == n_rows(), ExcDimensionMismatch(n_cols(), n_rows())); Assert ((n_cols()>=1) && (n_cols()<=3), ExcNotImplemented(n_cols())); - switch (n_cols()) + switch (this->n_cols()) { case 1: - return el(0,0); + return this->el(0,0); case 2: - return el(0,0)*el(1,1) - el(1,0)*el(0,1); + return this->el(0,0)*this->el(1,1) - this->el(1,0)*this->el(0,1); case 3: - return (el(0,0)*el(1,1)*el(2,2) - -el(0,0)*el(1,2)*el(2,1) - -el(1,0)*el(0,1)*el(2,2) - +el(1,0)*el(0,2)*el(2,1) - +el(2,0)*el(0,1)*el(1,2) - -el(2,0)*el(0,2)*el(1,1)); + return (this->el(0,0)*this->el(1,1)*this->el(2,2) + -this->el(0,0)*this->el(1,2)*this->el(2,1) + -this->el(1,0)*this->el(0,1)*this->el(2,2) + +this->el(1,0)*this->el(0,2)*this->el(2,1) + +this->el(2,0)*this->el(0,1)*this->el(1,2) + -this->el(2,0)*this->el(0,2)*this->el(1,1)); default: return 0; }; @@ -1115,11 +1138,11 @@ template number FullMatrix::norm2 () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number s = 0.; - for (unsigned int i=0;in_rows()*this->n_cols(); ++i) + s += this->data()[i]*this->data()[i]; return std::sqrt(s); } @@ -1128,12 +1151,12 @@ template number FullMatrix::relative_symmetry_norm2 () const { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); number s = 0.; number a = 0.; - for (unsigned int i=0;in_rows(); ++i) + for (unsigned int j=0; jn_cols(); ++j) { a += ((*this)(i,j)-(*this)(j,i))*((*this)(i,j)-(*this)(j,i)); s += (*this)(i,j)*(*this)(i,j); @@ -1149,7 +1172,7 @@ template void FullMatrix::invert (const FullMatrix &M) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (n_cols() == n_rows(), ExcNotQuadratic()); Assert (n_cols() == M.n_cols(), @@ -1157,20 +1180,20 @@ FullMatrix::invert (const FullMatrix &M) Assert (n_rows() == M.n_rows(), ExcDimensionMismatch(n_rows(),M.n_rows())); - switch (n_cols()) + switch (this->n_cols()) { case 1: - el(0,0) = 1.0/M.el(0,0); + this->el(0,0) = 1.0/M.el(0,0); return; case 2: // this is Maple output, // thus a bit unstructured { const number t4 = 1.0/(M.el(0,0)*M.el(1,1)-M.el(0,1)*M.el(1,0)); - el(0,0) = M.el(1,1)*t4; - el(0,1) = -M.el(0,1)*t4; - el(1,0) = -M.el(1,0)*t4; - el(1,1) = M.el(0,0)*t4; + this->el(0,0) = M.el(1,1)*t4; + this->el(0,1) = -M.el(0,1)*t4; + this->el(1,0) = -M.el(1,0)*t4; + this->el(1,1) = M.el(0,0)*t4; return; }; @@ -1184,15 +1207,15 @@ FullMatrix::invert (const FullMatrix &M) t04 = M.el(0,2)*M.el(2,0), t07 = 1.0/(t4*M.el(2,2)-t6*M.el(2,1)-t8*M.el(2,2)+ t00*M.el(2,1)+t01*M.el(1,2)-t04*M.el(1,1)); - el(0,0) = (M.el(1,1)*M.el(2,2)-M.el(1,2)*M.el(2,1))*t07; - el(0,1) = -(M.el(0,1)*M.el(2,2)-M.el(0,2)*M.el(2,1))*t07; - el(0,2) = -(-M.el(0,1)*M.el(1,2)+M.el(0,2)*M.el(1,1))*t07; - el(1,0) = -(M.el(1,0)*M.el(2,2)-M.el(1,2)*M.el(2,0))*t07; - el(1,1) = (M.el(0,0)*M.el(2,2)-t04)*t07; - el(1,2) = -(t6-t00)*t07; - el(2,0) = -(-M.el(1,0)*M.el(2,1)+M.el(1,1)*M.el(2,0))*t07; - el(2,1) = -(M.el(0,0)*M.el(2,1)-t01)*t07; - el(2,2) = (t4-t8)*t07; + this->el(0,0) = (M.el(1,1)*M.el(2,2)-M.el(1,2)*M.el(2,1))*t07; + this->el(0,1) = -(M.el(0,1)*M.el(2,2)-M.el(0,2)*M.el(2,1))*t07; + this->el(0,2) = -(-M.el(0,1)*M.el(1,2)+M.el(0,2)*M.el(1,1))*t07; + this->el(1,0) = -(M.el(1,0)*M.el(2,2)-M.el(1,2)*M.el(2,0))*t07; + this->el(1,1) = (M.el(0,0)*M.el(2,2)-t04)*t07; + this->el(1,2) = -(t6-t00)*t07; + this->el(2,0) = -(-M.el(1,0)*M.el(2,1)+M.el(1,1)*M.el(2,0))*t07; + this->el(2,1) = -(M.el(0,0)*M.el(2,1)-t01)*t07; + this->el(2,2) = (t4-t8)*t07; return; }; @@ -1260,40 +1283,40 @@ FullMatrix::invert (const FullMatrix &M) const number t131 = M.el(0,0)*M.el(1,3); const number t133 = M.el(1,0)*M.el(0,2); const number t135 = M.el(1,0)*M.el(0,3); - el(0,0) = (M.el(1,1)*M.el(2,2)*M.el(3,3)-M.el(1,1)*M.el(2,3)*M.el(3,2)- - M.el(2,1)*M.el(1,2)*M.el(3,3)+M.el(2,1)*M.el(1,3)*M.el(3,2)+ - M.el(3,1)*M.el(1,2)*M.el(2,3)-M.el(3,1)*M.el(1,3)*M.el(2,2))*t65; - el(0,1) = -(M.el(0,1)*M.el(2,2)*M.el(3,3)-M.el(0,1)*M.el(2,3)*M.el(3,2)- - t71*M.el(3,3)+t73*M.el(3,2)+t75*M.el(2,3)-t77*M.el(2,2))*t65; - el(0,2) = (t81*M.el(3,3)-t83*M.el(3,2)-t85*M.el(3,3)+t87*M.el(3,2)+ - t75*M.el(1,3)-t77*M.el(1,2))*t65; - el(0,3) = -(t81*M.el(2,3)-t83*M.el(2,2)-t85*M.el(2,3)+t87*M.el(2,2)+ - t71*M.el(1,3)-t73*M.el(1,2))*t65; - el(1,0) = -(t101*M.el(3,3)-t103*M.el(3,2)-t105*M.el(3,3)+t107*M.el(3,2)+ - t109*M.el(2,3)-t111*M.el(2,2))*t65; - el(1,1) = (t115*M.el(3,3)-t117*M.el(3,2)-t119*M.el(3,3)+t121*M.el(3,2)+ - t123*M.el(2,3)-t125*M.el(2,2))*t65; - el(1,2) = -(t129*M.el(3,3)-t131*M.el(3,2)-t133*M.el(3,3)+t135*M.el(3,2)+ - t123*M.el(1,3)-t125*M.el(1,2))*t65; - el(1,3) = (t129*M.el(2,3)-t131*M.el(2,2)-t133*M.el(2,3)+t135*M.el(2,2)+ - t119*M.el(1,3)-t121*M.el(1,2))*t65; - el(2,0) = (t32*M.el(3,3)-t103*M.el(3,1)-t46*M.el(3,3)+t107*M.el(3,1)+ - t57*M.el(2,3)-t111*M.el(2,1))*t65; - el(2,1) = -(t19*M.el(3,3)-t117*M.el(3,1)-t43*M.el(3,3)+t121*M.el(3,1)+ - t54*M.el(2,3)-t125*M.el(2,1))*t65; - el(2,2) = (t14*M.el(3,3)-t131*M.el(3,1)-t29*M.el(3,3)+t135*M.el(3,1)+ - t54*M.el(1,3)-t125*M.el(1,1))*t65; - el(2,3) = -(t14*M.el(2,3)-t131*M.el(2,1)-t29*M.el(2,3)+t135*M.el(2,1)+ - t43*M.el(1,3)-t121*M.el(1,1))*t65; - el(3,0) = -(t32*M.el(3,2)-t101*M.el(3,1)-t46*M.el(3,2)+t105*M.el(3,1)+ - t57*M.el(2,2)-t109*M.el(2,1))*t65; - el(3,1) = (t19*M.el(3,2)-t115*M.el(3,1)-t43*M.el(3,2)+t119*M.el(3,1)+ - t54*M.el(2,2)-t123*M.el(2,1))*t65; - el(3,2) = -(t14*M.el(3,2)-t129*M.el(3,1)-t29*M.el(3,2)+t133*M.el(3,1)+ - t54*M.el(1,2)-t123*M.el(1,1))*t65; - el(3,3) = (t14*M.el(2,2)-t129*M.el(2,1)-t29*M.el(2,2)+t133*M.el(2,1)+ - t43*M.el(1,2)-t119*M.el(1,1))*t65; - + this->el(0,0) = (M.el(1,1)*M.el(2,2)*M.el(3,3)-M.el(1,1)*M.el(2,3)*M.el(3,2)- + M.el(2,1)*M.el(1,2)*M.el(3,3)+M.el(2,1)*M.el(1,3)*M.el(3,2)+ + M.el(3,1)*M.el(1,2)*M.el(2,3)-M.el(3,1)*M.el(1,3)*M.el(2,2))*t65; + this->el(0,1) = -(M.el(0,1)*M.el(2,2)*M.el(3,3)-M.el(0,1)*M.el(2,3)*M.el(3,2)- + t71*M.el(3,3)+t73*M.el(3,2)+t75*M.el(2,3)-t77*M.el(2,2))*t65; + this->el(0,2) = (t81*M.el(3,3)-t83*M.el(3,2)-t85*M.el(3,3)+t87*M.el(3,2)+ + t75*M.el(1,3)-t77*M.el(1,2))*t65; + this->el(0,3) = -(t81*M.el(2,3)-t83*M.el(2,2)-t85*M.el(2,3)+t87*M.el(2,2)+ + t71*M.el(1,3)-t73*M.el(1,2))*t65; + this->el(1,0) = -(t101*M.el(3,3)-t103*M.el(3,2)-t105*M.el(3,3)+t107*M.el(3,2)+ + t109*M.el(2,3)-t111*M.el(2,2))*t65; + this->el(1,1) = (t115*M.el(3,3)-t117*M.el(3,2)-t119*M.el(3,3)+t121*M.el(3,2)+ + t123*M.el(2,3)-t125*M.el(2,2))*t65; + this->el(1,2) = -(t129*M.el(3,3)-t131*M.el(3,2)-t133*M.el(3,3)+t135*M.el(3,2)+ + t123*M.el(1,3)-t125*M.el(1,2))*t65; + this->el(1,3) = (t129*M.el(2,3)-t131*M.el(2,2)-t133*M.el(2,3)+t135*M.el(2,2)+ + t119*M.el(1,3)-t121*M.el(1,2))*t65; + this->el(2,0) = (t32*M.el(3,3)-t103*M.el(3,1)-t46*M.el(3,3)+t107*M.el(3,1)+ + t57*M.el(2,3)-t111*M.el(2,1))*t65; + this->el(2,1) = -(t19*M.el(3,3)-t117*M.el(3,1)-t43*M.el(3,3)+t121*M.el(3,1)+ + t54*M.el(2,3)-t125*M.el(2,1))*t65; + this->el(2,2) = (t14*M.el(3,3)-t131*M.el(3,1)-t29*M.el(3,3)+t135*M.el(3,1)+ + t54*M.el(1,3)-t125*M.el(1,1))*t65; + this->el(2,3) = -(t14*M.el(2,3)-t131*M.el(2,1)-t29*M.el(2,3)+t135*M.el(2,1)+ + t43*M.el(1,3)-t121*M.el(1,1))*t65; + this->el(3,0) = -(t32*M.el(3,2)-t101*M.el(3,1)-t46*M.el(3,2)+t105*M.el(3,1)+ + t57*M.el(2,2)-t109*M.el(2,1))*t65; + this->el(3,1) = (t19*M.el(3,2)-t115*M.el(3,1)-t43*M.el(3,2)+t119*M.el(3,1)+ + t54*M.el(2,2)-t123*M.el(2,1))*t65; + this->el(3,2) = -(t14*M.el(3,2)-t129*M.el(3,1)-t29*M.el(3,2)+t133*M.el(3,1)+ + t54*M.el(1,2)-t123*M.el(1,1))*t65; + this->el(3,3) = (t14*M.el(2,2)-t129*M.el(2,1)-t29*M.el(2,2)+t133*M.el(2,1)+ + t43*M.el(1,2)-t119*M.el(1,1))*t65; + break; } @@ -1324,7 +1347,7 @@ FullMatrix::precondition_Jacobi (Vector &dst, const somenumber *src_ptr = src.begin(); for (unsigned int i=0; iel(i,i); }; @@ -1340,7 +1363,7 @@ FullMatrix::print_formatted (std::ostream &out, { unsigned int width = width_; - Assert ((data() != 0) || (n_cols()+n_rows()==0), + Assert ((this->data() != 0) || (n_cols()+n_rows()==0), ExcInternalError()); // set output format, but store old @@ -1362,9 +1385,9 @@ FullMatrix::print_formatted (std::ostream &out, for (unsigned int i=0; iel(i,j) != 0) out << std::setw(width) - << el(i,j) * denominator << ' '; + << this->el(i,j) * denominator << ' '; else out << std::setw(width) << zero_string << ' '; out << std::endl; @@ -1382,7 +1405,7 @@ template void FullMatrix::gauss_jordan() { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); Assert (n_cols() == n_rows(), ExcNotQuadratic()); // Gauss-Jordan-Algorithmus @@ -1398,7 +1421,7 @@ FullMatrix::gauss_jordan() // regular double diagonal_sum = 0; for (unsigned int i=0; iel(i,i)); const double typical_diagonal_element = diagonal_sum/N; std::vector p(N); @@ -1411,13 +1434,13 @@ FullMatrix::gauss_jordan() // part of the line on and // right of the diagonal for // the largest element - number max = std::fabs(el(j,j)); + number max = std::fabs(this->el(j,j)); unsigned int r = j; for (unsigned int i=j+1; i max) + if (std::fabs(this->el(i,j)) > max) { - max = std::fabs(el(i,j)); + max = std::fabs(this->el(i,j)); r = i; } } @@ -1430,38 +1453,38 @@ FullMatrix::gauss_jordan() if (r>j) { for (unsigned int k=0; kel(j,k), this->el(r,k)); std::swap (p[j], p[r]); } // transformation - const number hr = 1./el(j,j); - el(j,j) = hr; + const number hr = 1./this->el(j,j); + this->el(j,j) = hr; for (unsigned int k=0; kel(i,k) -= this->el(i,j)*this->el(j,k)*hr; } } for (unsigned int i=0; iel(i,j) *= hr; + this->el(j,i) *= -hr; } - el(j,j) = hr; + this->el(j,j) = hr; } // column interchange std::vector hv(N); for (unsigned int i=0; iel(i,k); for (unsigned int k=0; kel(i,k) = hv[k]; } } @@ -1472,7 +1495,7 @@ template void FullMatrix::householder(Vector& src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); // m > n, src.n() = m Assert (n_cols() <= n_rows(), ExcDimensionMismatch(n_cols(), n_rows())); @@ -1482,30 +1505,30 @@ FullMatrix::householder(Vector& src) { number2 sigma = 0; unsigned int i; - for (i=j ; iel(i,j)*this->el(i,j); if (std::fabs(sigma) < 1.e-15) return; - number2 s = el(j,j); + number2 s = this->el(j,j); s = (s<0) ? std::sqrt(sigma) : -std::sqrt(sigma); number2 dj = s; - number2 beta = 1./(s*el(j,j)-sigma); - el(j,j) -= s; + number2 beta = 1./(s*this->el(j,j)-sigma); + this->el(j,j) -= s; for (unsigned int k=j+1 ; kel(i,j)*this->el(i,k); sum *= beta; - for (i=j ; iel(i,k) += sum*this->el(i,j); } number2 sum = 0.; - for (i=j ; iel(i,j)*src(i); sum *= beta; - for (i=j ; iel(i,j); + this->el(j,j) = dj; } } @@ -1515,7 +1538,7 @@ template double FullMatrix::least_squares(Vector& dst, Vector& src) { - Assert (data() != 0, ExcEmptyMatrix()); + Assert (this->data() != 0, ExcEmptyMatrix()); // m > n, m = src.n, n = dst.n diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 2728d5aec7..70f1762f9f 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -74,8 +74,8 @@ void SparseILU::decompose (const SparseMatrix &matrix, if (true) { // preset the elements - std::fill_n (&global_entry(0), - n_nonzero_elements(), + std::fill_n (&this->global_entry(0), + this->n_nonzero_elements(), 0); // note: pointers to the sparsity @@ -85,31 +85,34 @@ void SparseILU::decompose (const SparseMatrix &matrix, const unsigned int * const column_numbers = matrix.get_sparsity_pattern().get_column_numbers(); - for (unsigned int row=0; rowm(); ++row) for (const unsigned int * col = &column_numbers[rowstart_indices[row]]; col != &column_numbers[rowstart_indices[row+1]]; ++col) set (row, *col, matrix.global_entry(col-column_numbers)); }; if (strengthen_diagonal > 0) - for (unsigned int row=0; rowm(); ++row) { // get the length of the row // (without the diagonal element) - const unsigned int rowlength = get_sparsity_pattern().get_rowstart_indices()[row+1] - -get_sparsity_pattern().get_rowstart_indices()[row] - -1; + const unsigned int + rowlength = (this->get_sparsity_pattern().get_rowstart_indices()[row+1] + - + this->get_sparsity_pattern().get_rowstart_indices()[row] + - + 1); // get the global index of the first // non-diagonal element in this row const unsigned int rowstart - = get_sparsity_pattern().get_rowstart_indices()[row] + 1; - number * const diagonal_element = &global_entry(rowstart-1); + = this->get_sparsity_pattern().get_rowstart_indices()[row] + 1; + number * const diagonal_element = &this->global_entry(rowstart-1); number rowsum = 0; for (unsigned int global_index=rowstart; global_indexglobal_entry(global_index)); *diagonal_element += strengthen_diagonal * rowsum; }; @@ -117,7 +120,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, // now work only on this // matrix - const SparsityPattern &sparsity = get_sparsity_pattern(); + const SparsityPattern &sparsity = this->get_sparsity_pattern(); const unsigned int * const rowstart_indices = sparsity.get_rowstart_indices(); const unsigned int * const column_numbers = sparsity.get_column_numbers(); @@ -138,7 +141,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, // i := row - for (unsigned int row=1; rowm(); ++row) { // invert diagonal element of the // previous row. this is a hack, @@ -148,10 +151,10 @@ void SparseILU::decompose (const SparseMatrix &matrix, // and since it makes the backward // step when applying the decomposition // significantly faster - AssertThrow((global_entry(rowstart_indices[row-1]) !=0), + AssertThrow((this->global_entry(this->rowstart_indices[row-1]) !=0), ExcDivideByZero()); - global_entry (rowstart_indices[row-1]) + this->global_entry (rowstart_indices[row-1]) = 1./global_entry (rowstart_indices[row-1]); // let k run over all lower-left @@ -168,7 +171,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, for (const unsigned int * col_ptr = first_of_row; col_ptr!=first_after_diagonal; ++col_ptr) { const unsigned int global_index_ik = col_ptr-column_numbers; - global_entry(global_index_ik) *= diag_element(*col_ptr); + this->global_entry(global_index_ik) *= this->diag_element(*col_ptr); // now do the inner loop over // j. note that we need to do @@ -180,8 +183,8 @@ void SparseILU::decompose (const SparseMatrix &matrix, const int global_index_ki = sparsity(*col_ptr,row); if (global_index_ki != -1) - diag_element(row) -= global_entry(global_index_ik) * - global_entry(global_index_ki); + this->diag_element(row) -= this->global_entry(global_index_ik) * + this->global_entry(global_index_ki); for (const unsigned int * j = col_ptr+1; j<&column_numbers[rowstart_indices[row+1]]; @@ -202,8 +205,8 @@ void SparseILU::decompose (const SparseMatrix &matrix, global_index_kj = sparsity(*col_ptr,*j); if ((global_index_ij != -1) && (global_index_kj != -1)) - global_entry(global_index_ij) -= global_entry(global_index_ik) * - global_entry(global_index_kj); + this->global_entry(global_index_ij) -= this->global_entry(global_index_ik) * + this->global_entry(global_index_kj); }; }; }; @@ -212,7 +215,7 @@ void SparseILU::decompose (const SparseMatrix &matrix, // element still has to be inverted // because the for-loop doesn't do // it... - diag_element(m()-1)=1./diag_element(m()-1); + this->diag_element(this->m()-1) = 1./this->diag_element(m()-1); /* OLD CODE, rather crude first implementation with an algorithm taken @@ -271,8 +274,10 @@ void SparseILU::apply_decomposition (Vector &dst, Assert (dst.size() == m(), ExcSizeMismatch(dst.size(), m())); const unsigned int N=dst.size(); - const unsigned int * const rowstart_indices = get_sparsity_pattern().get_rowstart_indices(); - const unsigned int * const column_numbers = get_sparsity_pattern().get_column_numbers(); + const unsigned int * const rowstart_indices + = this->get_sparsity_pattern().get_rowstart_indices(); + const unsigned int * const column_numbers + = this->get_sparsity_pattern().get_column_numbers(); // solve LUx=b in two steps: // first Ly = b, then // Ux = y @@ -299,7 +304,7 @@ void SparseILU::apply_decomposition (Vector &dst, row); for (const unsigned int * col=rowstart; col!=first_after_diagonal; ++col) - dst(row) -= global_entry (col-column_numbers) * dst(*col); + dst(row) -= this->global_entry (col-column_numbers) * dst(*col); }; // now the backward solve. same @@ -322,12 +327,12 @@ void SparseILU::apply_decomposition (Vector &dst, static_cast(row)); for (const unsigned int * col=first_after_diagonal; col!=rowend; ++col) - dst(row) -= global_entry (col-column_numbers) * dst(*col); + dst(row) -= this->global_entry (col-column_numbers) * dst(*col); // scale by the diagonal element. // note that the diagonal element // was stored inverted - dst(row) *= diag_element(row); + dst(row) *= this->diag_element(row); }; }; -- 2.39.5