00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__householder_h
00013 #define __deal2__householder_h
00014
00015
00016 #include <cmath>
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/lac/full_matrix.h>
00019 #include <deal.II/lac/vector_memory.h>
00020
00021 #include <vector>
00022
00023 DEAL_II_NAMESPACE_OPEN
00024
00025
00026
00027 template<typename number> class Vector;
00028
00029
00050 template<typename number>
00051 class Householder : private FullMatrix<number>
00052 {
00053 public:
00057 Householder ();
00058
00063 template<typename number2>
00064 Householder (const FullMatrix<number2>&);
00065
00070 template<typename number2>
00071 void
00072 initialize (const FullMatrix<number2>&);
00073
00091 template<typename number2>
00092 double least_squares (Vector<number2> &dst,
00093 const Vector<number2> &src) const;
00094
00099 template<typename number2>
00100 double least_squares (BlockVector<number2> &dst,
00101 const BlockVector<number2> &src) const;
00102
00108 template<class VECTOR>
00109 void vmult (VECTOR &dst, const VECTOR &src) const;
00110
00111 template<class VECTOR>
00112 void Tvmult (VECTOR &dst, const VECTOR &src) const;
00113
00114
00115 private:
00121 std::vector<number> diagonal;
00122 };
00123
00126 #ifndef DOXYGEN
00127
00128
00129
00130
00131 template <typename number>
00132 Householder<number>::Householder()
00133 {}
00134
00135
00136
00137 template <typename number>
00138 template <typename number2>
00139 void
00140 Householder<number>::initialize(const FullMatrix<number2>& M)
00141 {
00142 const unsigned int m = M.n_rows(), n = M.n_cols();
00143 this->reinit(m, n);
00144 this->fill(M);
00145 Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00146 diagonal.resize(m);
00147
00148
00149 Assert (this->n_cols() <= this->n_rows(),
00150 ExcDimensionMismatch(this->n_cols(), this->n_rows()));
00151
00152 for (unsigned int j=0 ; j<n ; ++j)
00153 {
00154 number2 sigma = 0;
00155 unsigned int i;
00156
00157 for (i=j ; i<m ; ++i)
00158 sigma += this->el(i,j)*this->el(i,j);
00159
00160
00161 if (std::fabs(sigma) < 1.e-15) return;
00162
00163 number2 s = (this->el(j,j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
00164
00165 number2 beta = std::sqrt(1./(sigma-s*this->el(j,j)));
00166
00167
00168
00169
00170 diagonal[j] = beta*(this->el(j,j) - s);
00171 this->el(j,j) = s;
00172
00173 for (i=j+1 ; i<m ; ++i)
00174 this->el(i,j) *= beta;
00175
00176
00177
00178
00179 for (unsigned int k=j+1 ; k<n ; ++k)
00180 {
00181 number2 sum = diagonal[j]*this->el(j,k);
00182 for (i=j+1 ; i<m ; ++i)
00183 sum += this->el(i,j)*this->el(i,k);
00184
00185 this->el(j,k) -= sum*this->diagonal[j];
00186 for (i=j+1 ; i<m ; ++i)
00187 this->el(i,k) -= sum*this->el(i,j);
00188 }
00189 }
00190 }
00191
00192
00193 template <typename number>
00194 template <typename number2>
00195 Householder<number>::Householder(const FullMatrix<number2>& M)
00196 {
00197 initialize(M);
00198 }
00199
00200
00201 template <typename number>
00202 template <typename number2>
00203 double
00204 Householder<number>::least_squares (Vector<number2>& dst,
00205 const Vector<number2>& src) const
00206 {
00207 Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00208 AssertDimension(dst.size(), this->n());
00209 AssertDimension(src.size(), this->m());
00210
00211 const unsigned int m = this->m(), n = this->n();
00212
00213 GrowingVectorMemory<Vector<number2> > mem;
00214 Vector<number2>* aux = mem.alloc();
00215 aux->reinit(src, true);
00216 *aux = src;
00217
00218
00219
00220
00221 for (unsigned int j=0;j<n;++j)
00222 {
00223
00224 number2 sum = diagonal[j]* (*aux)(j);
00225 for (unsigned int i=j+1 ; i<m ; ++i)
00226 sum += this->el(i,j)*(*aux)(i);
00227
00228 (*aux)(j) -= sum*diagonal[j];
00229 for (unsigned int i=j+1 ; i<m ; ++i)
00230 (*aux)(i) -= sum*this->el(i,j);
00231 }
00232
00233 number2 sum = 0.;
00234 for (unsigned int i=n ; i<m ; ++i)
00235 sum += (*aux)(i) * (*aux)(i);
00236
00237 this->backward(dst, *aux);
00238
00239 mem.free(aux);
00240
00241 return std::sqrt(sum);
00242 }
00243
00244 template <typename number>
00245 template <typename number2>
00246 double
00247 Householder<number>::least_squares (BlockVector<number2>& dst,
00248 const BlockVector<number2>& src) const
00249 {
00250 Assert (!this->empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
00251 AssertDimension(dst.size(), this->n());
00252 AssertDimension(src.size(), this->m());
00253
00254 const unsigned int m = this->m(), n = this->n();
00255
00256 GrowingVectorMemory<BlockVector<number2> > mem;
00257 BlockVector<number2>* aux = mem.alloc();
00258 aux->reinit(src, true);
00259 *aux = src;
00260
00261
00262
00263
00264 for (unsigned int j=0;j<n;++j)
00265 {
00266
00267 number2 sum = diagonal[j]* (*aux)(j);
00268 for (unsigned int i=j+1 ; i<m ; ++i)
00269 sum += this->el(i,j)*(*aux)(i);
00270
00271 (*aux)(j) -= sum*diagonal[j];
00272 for (unsigned int i=j+1 ; i<m ; ++i)
00273 (*aux)(i) -= sum*this->el(i,j);
00274 }
00275
00276 number2 sum = 0.;
00277 for (unsigned int i=n ; i<m ; ++i)
00278 sum += (*aux)(i) * (*aux)(i);
00279
00280
00281
00282 Vector<number2> v_dst, v_aux;
00283 v_dst = dst;
00284 v_aux = *aux;
00285
00286 this->backward(v_dst, v_aux);
00287
00288
00289 dst = v_dst;
00290
00291 mem.free(aux);
00292
00293 return std::sqrt(sum);
00294 }
00295
00296
00297 template <typename number>
00298 template <class VECTOR>
00299 void
00300 Householder<number>::vmult (VECTOR &dst, const VECTOR &src) const
00301 {
00302 least_squares (dst, src);
00303 }
00304
00305
00306 template <typename number>
00307 template <class VECTOR>
00308 void
00309 Householder<number>::Tvmult (VECTOR &, const VECTOR &) const
00310 {
00311 Assert(false, ExcNotImplemented());
00312 }
00313
00314
00315
00316 #endif // DOXYGEN
00317
00318 DEAL_II_NAMESPACE_CLOSE
00319
00320 #endif
00321
00322