From 2c558291112a26d48c99e5483b05f0f3fc8ad729 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 28 May 2009 09:55:05 +0000 Subject: [PATCH] Make some small optimization in vmult of FullMatrix. git-svn-id: https://svn.dealii.org/trunk@18884 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/full_matrix.templates.h | 66 +++++++------------ deal.II/lac/include/lac/householder.h | 29 ++++---- .../lac/include/lac/sparse_matrix.templates.h | 4 +- 3 files changed, 44 insertions(+), 55 deletions(-) diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 936884e238..4b7744d5e0 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -198,27 +198,17 @@ FullMatrix::vmult (Vector& dst, Assert (&src != &dst, ExcSourceEqualsDestination()); const number* e = this->data(); - const unsigned int size_m = m(), - size_n = n(); - if (!adding) - { - for (unsigned int i=0; i*>(&src))(0); + const unsigned int size_m = m(), size_n = n(); + for (unsigned int i=0; i::Tvmult (Vector &dst, const bool adding) const { Assert (!this->empty(), ExcEmptyMatrix()); - + Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m())); Assert (&src != &dst, ExcSourceEqualsDestination()); - const unsigned int size_m = m(), - size_n = n(); + const number* e = this->data(); + number2* dst_ptr = &dst(0); + const unsigned int size_m = m(), size_n = n(); + // zero out data if we are not adding if (!adding) + for (unsigned int j=0; j::Householder() {} + template template void Householder::initialize(const FullMatrix& M) { - this->reinit(M.n_rows(), M.n_cols()); + const unsigned int m = M.n_rows(), n = M.n_cols(); + this->reinit(m, n); this->fill(M); - diagonal.resize(M.n_rows()); Assert (!this->empty(), typename FullMatrix::ExcEmptyMatrix()); + diagonal.resize(m); + // m > n, src.n() = m Assert (this->n_cols() <= this->n_rows(), ExcDimensionMismatch(this->n_cols(), this->n_rows())); - for (unsigned int j=0 ; jn() ; ++j) + for (unsigned int j=0 ; jm() ; ++i) + for (i=j ; iel(i,j)*this->el(i,j); // We are ready if the column is // empty. Are we? @@ -148,20 +151,20 @@ Householder::initialize(const FullMatrix& M) diagonal[j] = beta*(this->el(j,j) - s); this->el(j,j) = s; - for (i=j+1 ; im() ; ++i) + for (i=j+1 ; iel(i,j) *= beta; // For all subsequent columns do // the Householder reflexion - for (unsigned int k=j+1 ; kn() ; ++k) + for (unsigned int k=j+1 ; kel(j,k); - for (i=j+1 ; im() ; ++i) + for (i=j+1 ; iel(i,j)*this->el(i,k); this->el(j,k) -= sum*this->diagonal[j]; - for (i=j+1 ; im() ; ++i) + for (i=j+1 ; iel(i,k) -= sum*this->el(i,j); } } @@ -185,6 +188,8 @@ Householder::least_squares (Vector& dst, Assert (!this->empty(), typename FullMatrix::ExcEmptyMatrix()); AssertDimension(dst.size(), this->n()); AssertDimension(src.size(), this->m()); + + const unsigned int m = this->m(), n = this->n(); GrowingVectorMemory > mem; Vector* aux = mem.alloc(); @@ -194,20 +199,20 @@ Householder::least_squares (Vector& dst, // Multiply Q_n ... Q_2 Q_1 src // Where Q_i = I-v_i v_i^T - for (unsigned int j=0;jn();++j) + for (unsigned int j=0;jm() ; ++i) + for (unsigned int i=j+1 ; iel(i,j)*(*aux)(i); // dst -= v * sum (*aux)(j) -= sum*diagonal[j]; - for (unsigned int i=j+1 ; im() ; ++i) + for (unsigned int i=j+1 ; iel(i,j); } // Compute norm of residual number2 sum = 0.; - for (unsigned int i=this->n() ; im() ; ++i) + for (unsigned int i=n ; ibackward(dst, *aux); diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index df34ec9546..5b7ee61e36 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -374,11 +374,11 @@ namespace internal else for (unsigned int row=begin_row; row