QuOc

 

preconditioner.h

Go to the documentation of this file.
00001 /* Copyright (c) 2001-2010 AG Rumpf, INS, Universitaet Bonn                      */
00002 /*                                                                               */
00003 /* The contents of this file are subject to the terms of the Common Development  */
00004 /* and Distribution License Version 1.0 (the "License"); you may not use         */
00005 /* this file except in compliance with the License. You may obtain a copy of     */
00006 /* the License at http://www.sun.com/cddl/                                       */
00007 /*                                                                               */
00008 /* Software distributed under the License is distributed on an "AS IS" basis,    */
00009 /* WITHOUT WARRANTY OF ANY KIND, either express or implied.                      */
00010 /*                                                                               */
00011 #ifndef __PRECONDITIONER_H
00012 #define __PRECONDITIONER_H
00013 
00014 #include <solver.h>
00015 
00016 namespace aol {
00017 
00018 /** Basis class for template specialization */
00019 template< typename VectorType >
00020 class DiagonalPreconditioner;
00021 
00022 
00023 /**
00024  * \brief A class for diagonally preconditioning matrices (and matrices only)
00025  * \author Schwen
00026  * \ingroup solver
00027  */
00028 template< typename RealType >
00029 class DiagonalPreconditioner < aol::Vector<RealType> > : public Op< aol::Vector<RealType> > {
00030 protected:
00031   aol::DiagonalMatrix< RealType > _diagInv;
00032 
00033 public:
00034   template< typename MatrixType >
00035   explicit DiagonalPreconditioner ( const MatrixType &Matrix ) : _diagInv ( Matrix.getNumRows() ) {
00036     _diagInv.setToInverseDiagonalOf ( Matrix );
00037   }
00038 
00039   virtual void applyAdd ( const Vector<RealType> &Arg, Vector<RealType> &Dest ) const {
00040     _diagInv.applyAdd ( Arg, Dest );
00041   }
00042 
00043   virtual void apply ( const Vector<RealType> &Arg, Vector<RealType> &Dest ) const {
00044     _diagInv.apply ( Arg, Dest );
00045   }
00046   
00047   void getPreconditionerMatrix ( aol::DiagonalMatrix<RealType> &mat ) const {
00048     mat = _diagInv;
00049   }
00050 
00051   const aol::DiagonalMatrix<RealType> &getPreconditionerMatrixReference ( ) const {
00052     return _diagInv;
00053   }
00054 };
00055 
00056 
00057 /**
00058  * \brief A class for diagonally preconditioning BlockMatrices
00059  * \author Teusner, Berkels
00060  * \ingroup solver
00061  *
00062  * The diagonal entries of the BlockOp may contain NULL pointers.
00063  * In this case the corresponding block of the preconditioner is
00064  * set to the IdentityOp and getPreconditionerBlockMatrix can't be
00065  * called.
00066  */
00067 template< typename RealType >
00068 class DiagonalPreconditioner < aol::MultiVector<RealType> > : public BlockOp< RealType, aol::Op< aol::Vector<RealType> > > {
00069   aol::IdentityOp<aol::Vector<RealType> > _identity;
00070 
00071 public:
00072   template< typename BlockMatrixType >
00073   explicit DiagonalPreconditioner ( const BlockMatrixType &BlockMat ) : aol::BlockOp< RealType, aol::Op< aol::Vector<RealType> > > ( BlockMat.getNumRows(), BlockMat.getNumCols() ) {
00074     for ( int i = 0; i < BlockMat.getNumRows(); i++ ) {
00075       if ( BlockMat.getPointer ( i, i ) != NULL )
00076         this->setReference ( i, i, * ( new DiagonalPreconditioner< aol::Vector<RealType> > ( BlockMat.getReference ( i, i ) ) ) );
00077       else
00078         this->setReference ( i, i, _identity );
00079     }
00080   }
00081 
00082 public:
00083   ~DiagonalPreconditioner () {
00084     for ( int i = 0; i < this->getNumRows(); i++ ) {
00085       if ( this->getPointer ( i, i ) != &_identity )
00086         delete ( this->getPointer ( i, i ) );
00087     }
00088   }
00089 
00090   void getPreconditionerBlockMatrix ( aol::BlockMatrix < aol::DiagonalMatrix< RealType > > &blockMat ) const {
00091     for ( int i = 0; i < this->getNumRows(); i++ ) {
00092       const DiagonalPreconditioner< aol::Vector<RealType> >* pTemp = dynamic_cast< const DiagonalPreconditioner< aol::Vector<RealType> >* >( this->getPointer(i, i) );
00093       if ( pTemp )
00094         pTemp->getPreconditionerMatrix( blockMat.getReference(i, i) );
00095       else
00096         throw Exception ( "Entry that is not of type DiagonalPreconditioner<RealType> found (most likely an IdenityOp), can't get the preconditioner matrix for this.\n", __FILE__, __LINE__ );
00097     }
00098   }
00099 
00100 };
00101 
00102 
00103 /** Basis class for template specialization */
00104 template< typename VectorType >
00105 class GeometricScalingPreconditioner;
00106 
00107 
00108 /**
00109  * \brief A class for preconditioning matrices using geometric scaling (division by l2 norm of rows, cf. Gordon & Gordon, arXiv:0812.2769v2 [cs.MS])
00110  * \author Schwen
00111  * \ingroup solver
00112  */
00113 template< typename RealType >
00114 class GeometricScalingPreconditioner < aol::Vector<RealType> > : public Op< aol::Vector<RealType> > {
00115 protected:
00116   aol::DiagonalMatrix< RealType > _diagInv;
00117 
00118 public:
00119   template< typename MatrixType >
00120   explicit GeometricScalingPreconditioner ( const MatrixType &Matrix ) : _diagInv ( Matrix.getNumRows() ) {
00121     for ( int i = 0; i < Matrix.getNumRows(); ++i ) {
00122       RealType l2sum = aol::NumberTrait<RealType>::zero;
00123       vector<typename Row<RealType>::RowEntry > vec;
00124       Matrix.makeRowEntries ( vec, i );
00125       for ( typename vector<typename Row<RealType>::RowEntry >::const_iterator it = vec.begin(); it != vec.end(); ++it ) {
00126         l2sum += aol::Sqr ( it->value );
00127       }
00128 
00129       // make sure we do not divide by zero
00130       if ( l2sum == aol::NumberTrait<RealType>::zero )
00131         l2sum = aol::NumberTrait<RealType>::one;
00132 
00133       _diagInv.set ( i, i, aol::NumberTrait<RealType>::one / sqrt ( l2sum ) );
00134     }
00135   }
00136 
00137   virtual void applyAdd ( const Vector<RealType> &Arg, Vector<RealType> &Dest ) const {
00138     _diagInv.applyAdd ( Arg, Dest );
00139   }
00140 
00141   virtual void apply ( const Vector<RealType> &Arg, Vector<RealType> &Dest ) const {
00142     _diagInv.apply ( Arg, Dest );
00143   }
00144   
00145   void getPreconditionerMatrix ( aol::DiagonalMatrix<RealType> &mat ) const {
00146     mat = _diagInv;
00147   }
00148 };
00149 
00150 
00151 /**
00152  * \brief A class for preconditioning block matrices using geometric scaling (division by l2 norm of unblocked rows)
00153  * \author Schwen
00154  * \ingroup solver
00155  */
00156 template< typename RealType >
00157 class GeometricScalingPreconditioner < aol::MultiVector<RealType> > : public BlockOp< RealType, aol::Op< aol::Vector<RealType> > > {
00158 public:
00159   template< typename BlockMatrixType >
00160   explicit GeometricScalingPreconditioner ( const BlockMatrixType &BlockMat ) : aol::BlockOp< RealType, aol::Op< aol::Vector<RealType> > > ( BlockMat.getNumRows(), BlockMat.getNumCols() ) {
00161     for ( int i = 0; i < BlockMat.getNumRows(); i++ ) {
00162       DiagonalMatrix<RealType> *pCurrentDiag = new DiagonalMatrix< RealType > ( BlockMat.getReference ( i, i ).getNumRows() );
00163       this->setPointer ( i, i, pCurrentDiag );
00164 
00165       for ( int ii = 0; ii < pCurrentDiag->getNumRows(); ++ii ) {
00166         RealType l2sum = aol::NumberTrait<RealType>::zero;
00167         vector<typename Row<RealType>::RowEntry > vec;
00168         BlockMat.makeUnblockedRowEntries ( vec, i, ii );
00169         
00170         for ( typename vector<typename Row<RealType>::RowEntry >::const_iterator it = vec.begin(); it != vec.end(); ++it ) {
00171           l2sum += aol::Sqr ( it->value );
00172         }
00173         
00174         // make sure we do not divide by zero
00175         if ( l2sum == aol::NumberTrait<RealType>::zero )
00176           l2sum = aol::NumberTrait<RealType>::one;
00177         
00178         pCurrentDiag->set ( ii, ii, aol::NumberTrait<RealType>::one / sqrt ( l2sum ) );
00179       }
00180     }
00181   }
00182 
00183 public:
00184   ~GeometricScalingPreconditioner () {
00185     for ( int i = 0; i < this->getNumRows(); i++ ) {
00186       delete ( this->getPointer ( i, i ) );
00187     }
00188   }
00189 
00190   void getPreconditionerBlockMatrix ( aol::BlockMatrix < aol::DiagonalMatrix< RealType > > &blockMat ) const {
00191     for ( int i = 0; i < this->getNumRows(); i++ ) {
00192       blockMat.getReference ( i, i ) = dynamic_cast< const aol::DiagonalMatrix<RealType>& > ( this->getReference ( i, i ) );
00193     }
00194   }
00195 
00196 };
00197 
00198 
00199 /** Basis class for template specialization */
00200 template< typename VectorType, typename MatrixType >
00201 class SSORPreconditioner;
00202 
00203 
00204 /** \brief A class for preconditioning sparse matrices by SSOR.
00205  *  Represents for a given matrix \f$A\f$ and its decomposition
00206  *  \f$A=L+D+U\f$ into the lower triangular, diagonal, and upper triangular part
00207  *  the operator \f$\left(\frac{D}\omega+L\right)\frac\omega{2-\omega}D^{-1}\left(\frac{D}\omega+U\right)\f$,
00208  *  which is equivalent to computing one single SSOR iterate starting with the initial vector 0.
00209  *  Note, the SOR-iteration to solve the system \f$Ax=b\f$ is given by
00210  *  \f$\left(\frac{D}\omega+L\right)x^{m+1}=\frac{2-\omega}\omega Dx^m-\left(\frac{D}\omega+U\right)x^m+b\f$,
00211  *  the backward SOR-iteration by
00212  *  \f$\left(\frac{D}\omega+U\right)x^{m+1}=\frac{2-\omega}\omega Dx^m-\left(\frac{D}\omega+L\right)x^m+b\f$,
00213  *  and the SSOR-iteration by alternatingly applying the SOR-iteration first and then the backward SOR-iteration.
00214  *  \author Droske
00215  *  \ingroup solver
00216  */
00217 template <typename RealType, typename MatrixType>
00218 class SSORPreconditioner< aol::Vector<RealType>, MatrixType > : public Op<aol::Vector<RealType> > {
00219 protected:
00220   const MatrixType &_sparse;
00221   RealType _omega;
00222 public:
00223   SSORPreconditioner ( const MatrixType &SparseMatrix, RealType Omega = 1.2 )
00224       : _sparse ( SparseMatrix ), _omega ( Omega ) {}
00225 
00226   virtual void applyAdd ( const aol::Vector<RealType> &Arg, aol::Vector<RealType> &Dest ) const {
00227     aol::Vector<RealType> tmp ( Dest.size() );
00228     apply ( Arg, tmp );
00229     Dest += tmp;
00230   }
00231 
00232   virtual void apply ( const aol::Vector<RealType> &Arg, aol::Vector<RealType> &Dest ) const {
00233     aol::Vector<RealType> &x = Dest;
00234     const aol::Vector<RealType> &b = Arg;
00235 
00236     vector<typename Row<RealType>::RowEntry > vec;
00237     for ( int i = 0; i < x.size(); i++ ) {
00238       _sparse.makeRowEntries ( vec, i );
00239 
00240 #ifdef CHECK_ZERO_DIAG
00241       RealType v = 0., diag = 1.;
00242 #else
00243       RealType v = 0., diag = 0.;
00244 #endif
00245       for ( typename vector<typename Row<RealType>::RowEntry >::iterator it = vec.begin(); it != vec.end(); ++it ) {
00246 #ifdef BOUNDS_CHECK
00247         if ( it->col >= x.size() || it->col < 0 ) {
00248           cerr << "error: it-> col = " <<  it-> col << endl;
00249         }
00250 #endif
00251         if ( it->col == i ) {
00252           diag = it->value;
00253 #ifdef CHECK_ZERO_DIAG
00254           if ( Abs ( diag ) < 1e-20 ) {
00255             diag = 1.;
00256           }
00257 #endif
00258         } else if ( it->col < i ) {
00259           v += it->value * x[it->col];
00260         }
00261       }
00262       x[i] = ( b[i] - _omega * v ) / diag;
00263     }
00264 
00265     for ( int i = x.size() - 1; i >= 0; --i ) {
00266       _sparse.makeRowEntries ( vec, i );
00267 
00268 #ifdef CHECK_ZERO_DIAG
00269       RealType v = 0., diag = 1.;
00270 #else
00271       RealType v = 0., diag = 0.;
00272 #endif
00273       for ( typename vector<typename Row<RealType>::RowEntry >::iterator it = vec.begin(); it != vec.end(); ++it ) {
00274         if ( it->col == i ) {
00275           diag = it->value;
00276 #ifdef CHECK_ZERO_DIAG
00277           if ( Abs ( diag ) < 1e-20 ) {
00278             diag = 1.;
00279           }
00280 #endif
00281           x[i] *= diag;
00282         } else if ( it->col > i ) {
00283           v += it->value * x[it->col];
00284         }
00285       }
00286       x[i] = ( x[i] - _omega * v ) / diag;
00287     }
00288   }
00289 };
00290 
00291 /** \brief SSOR preconditioner for block matrices.
00292  *  Represents for a given matrix \f$A\f$ and its decomposition
00293  *  \f$A=L+D+U\f$ into the lower triangular, diagonal, and upper triangular part
00294  *  the operator \f$\left(\frac{D}\omega+L\right)\frac\omega{2-\omega}D^{-1}\left(\frac{D}\omega+U\right)\f$,
00295  *  which is equivalent to computing one single SSOR iterate starting with the initial vector 0.
00296  *  Note, the SOR-iteration to solve the system \f$Ax=b\f$ is given by
00297  *  \f$\left(\frac{D}\omega+L\right)x^{m+1}=\frac{2-\omega}\omega Dx^m-\left(\frac{D}\omega+U\right)x^m+b\f$,
00298  *  the backward SOR-iteration by
00299  *  \f$\left(\frac{D}\omega+U\right)x^{m+1}=\frac{2-\omega}\omega Dx^m-\left(\frac{D}\omega+L\right)x^m+b\f$,
00300  *  and the SSOR-iteration by alternatingly applying the SOR-iteration first and then the backward SOR-iteration.
00301  *  \author Wirth
00302  *  \ingroup solver
00303  */
00304 template <typename RealType, typename MatrixType>
00305 class SSORPreconditioner< aol::MultiVector<RealType>, MatrixType > :
00306 public aol::Op<MultiVector<RealType> > {
00307 protected:
00308   const aol::BlockOpBase<RealType, MatrixType> &_blockMat;
00309   RealType _omega;
00310 public:
00311   SSORPreconditioner ( const aol::BlockOpBase<RealType, MatrixType> &BlockMat, RealType Omega = 1.2 ) :
00312     _blockMat( BlockMat ),
00313     _omega( Omega ) {}
00314 
00315   virtual void applyAdd ( const MultiVector<RealType> &Arg, MultiVector<RealType> &Dest ) const {
00316     MultiVector<RealType> tmp ( Dest, aol::STRUCT_COPY );
00317     apply ( Arg, tmp );
00318     Dest += tmp;
00319   }
00320 
00321   virtual void apply ( const MultiVector<RealType> &Arg, MultiVector<RealType> &Dest ) const {
00322     // the following notation is used: L,D,U - lower triangular, diagonal, upper triangular part of _blockMat;
00323     // b - Arg; A = L+D+U
00324 
00325     vector<typename Row<RealType>::RowEntry > matRow;
00326 
00327     // solve (\omega L+D)y=b for y by forward substitution and write y into Dest
00328     for ( int k = 0; k < _blockMat.getNumRows(); k++ ) {
00329       // compute \sum_{j=1}^{k-1}A_{kj}y_j
00330       Vector<RealType> sum( Dest[k], aol::STRUCT_COPY );
00331       for ( int j = 0; j < k; j++ )
00332         _blockMat.getReference( k, j ).applyAdd( Dest[j], sum );
00333       // solve (\omega L_{kk}+D_{kk})y_k=b_k-\omega\sum_{j=1}^{k-1}A_{kj}y_j by forward substitution
00334       for ( int i = 0; i < Dest[k].size(); i++ ) {
00335         _blockMat.getReference( k, k ).makeRowEntries( matRow, i );
00336         RealType v = 0., diag = 1.;
00337         for ( typename vector<typename Row<RealType>::RowEntry >::iterator it = matRow.begin(); it != matRow.end(); ++it )
00338           if ( it->col < i )
00339             v += it->value * Dest[k][it->col];
00340           else if ( it->col == i )
00341             diag = it->value;
00342         Dest[k][i] = ( Arg[k][i] - _omega * ( sum[i] + v ) ) / diag;
00343       }
00344     }
00345 
00346     // solve (D+\omega U)x=Dy for x by backward substitution and write x into Dest
00347     for ( int k = _blockMat.getNumRows() - 1; k >= 0; k-- ) {
00348       // compute \sum_{j=k+1}^{n}A_{kj}x_j
00349       Vector<RealType> sum( Dest[k], aol::STRUCT_COPY );
00350       for ( int j = k+1; j < _blockMat.getNumCols(); j++ )
00351         _blockMat.getReference( k, j ).applyAdd( Dest[j], sum );
00352       // solve (D_{kk}+\omega U_{kk})x_k=D_{kk}y_k-\omega\sum_{j=k+1}^{n}A_{kj}x_j by backward substitution
00353       for ( int i = Dest[k].size() - 1; i >= 0; --i ) {
00354         _blockMat.getReference( k, k ).makeRowEntries( matRow, i );
00355         RealType v = 0., diag = 1.;
00356         for ( typename vector<typename Row<RealType>::RowEntry >::iterator it = matRow.begin(); it != matRow.end(); ++it )
00357           if ( it->col > i )
00358             v += it->value * Dest[k][it->col];
00359           else if ( it->col == i )
00360             diag = it->value;
00361         Dest[k][i] = ( diag * Dest[k][i] - _omega * ( sum[i] + v ) ) / diag;
00362       }
00363     }
00364   }
00365 };
00366 
00367 
00368 /** SSOR preconditioning which is naively implemented by performing
00369  *  one forward and one backward Block-Gauss-Seidel sweep. This does
00370  *  not have optimal complexity.
00371  *  \todo think about initialization of Dest (OS,BW)
00372  *  \author Schwen
00373  *  \ingroup solver
00374  */
00375 template< typename VectorType, typename MatrixType >
00376 class GaussSeidelPreconditioner : public Op< VectorType > {
00377 protected:
00378   typedef typename VectorType::DataType DataType;
00379   aol::GaussSeidelSweeper< DataType, VectorType, MatrixType > _sweeper;
00380   aol::GaussSeidelSweepingMode _firstSweepMode, _secondSweepMode;
00381   
00382 public:
00383   explicit GaussSeidelPreconditioner ( const MatrixType &Matrix, const DataType Relax = 1.2 ) : _sweeper ( Matrix, Relax ), _firstSweepMode ( aol::GAUSS_SEIDEL_FORWARD ), _secondSweepMode ( aol::GAUSS_SEIDEL_BACKWARD ) {
00384   }
00385 
00386   void setSweepingModes ( const aol::GaussSeidelSweepingMode firstMode, const aol::GaussSeidelSweepingMode secondMode ) {
00387     _firstSweepMode = firstMode;
00388     _secondSweepMode = secondMode;
00389   }
00390 
00391   virtual void applyAdd ( const VectorType &Arg, VectorType &Dest ) const {
00392     VectorType tmp ( Dest, aol::DEEP_COPY );
00393     // ? VectorType tmp ( Dest, aol::STRUCT_COPY );
00394     apply ( Arg, tmp );
00395     Dest += tmp;
00396   }
00397 
00398   virtual void apply ( const VectorType &Arg, VectorType &Dest ) const {
00399     // ? Dest.setZero();
00400     _sweeper.apply ( Arg, Dest, _firstSweepMode );
00401     _sweeper.apply ( Arg, Dest, _secondSweepMode );
00402   }
00403 };
00404 
00405 
00406 /** A class for preconditioning sparse matrices by ILU(0).
00407  *  ATTN: The preconditioner must compute a decomposition (in the constructor), so it uses a copy of the matrix
00408  *  The matrix must be filled at time of the precoditioners construction
00409  *  In the contructor get (,) is used for the needed column values, values are written via set/add (,)
00410  *  \author unknown
00411  *  \ingroup solver
00412  */
00413 template <typename RealType, typename MatrixType>
00414 class ILU0Preconditioner : public Op<aol::Vector<RealType> > {
00415 
00416 protected:
00417   MatrixType _decomp;
00418 
00419 public:
00420   ILU0Preconditioner ( const MatrixType& matrix )
00421       : _decomp ( matrix ) {
00422     vector<typename Row<RealType>::RowEntry > rowi;
00423     vector<typename Row<RealType>::RowEntry > rowk;
00424 
00425     for ( int i = 1; i < _decomp.getNumRows (); ++i ) { // Nothing is done in first row
00426       _decomp.makeSortedRowEntries ( rowi, i ); // Rowwise elimination
00427 
00428       for ( typename vector<typename Row<RealType>::RowEntry >::iterator ikentry = rowi.begin ();
00429             ikentry != rowi.end () && ikentry->col < i; ++ikentry ) { // Only left from diagonal, RowEntries are sorted
00430         int k = ikentry->col;
00431 
00432         RealType ikval = ikentry->value;
00433         _decomp.makeSortedRowEntries ( rowk, k );
00434 
00435         typename vector<typename Row<RealType>::RowEntry >::iterator kjentry = rowk.begin ();
00436 
00437         // Find index kk
00438         while ( kjentry != rowk.end () && kjentry->col < k ) ++kjentry;
00439         if ( kjentry == rowk.end () || kjentry->col != k ) throw Exception ( "Pivot zero in ILU", __FILE__, __LINE__ );
00440 
00441         // Eliminate index ik, store factor
00442         ikval /= kjentry->value;
00443         _decomp.set ( i, k, ikval );
00444 
00445         typename vector<typename Row<RealType>::RowEntry >::iterator ikentrynext = ikentry; ++ikentrynext; // Only right from current column
00446         for ( typename vector<typename Row<RealType>::RowEntry >::iterator ijentry = ikentrynext; ijentry != rowi.end (); ++ijentry ) {
00447           int j = ijentry->col;
00448 
00449           // Find index kj
00450           while ( kjentry != rowk.end () && kjentry->col < j ) ++kjentry;
00451           if ( kjentry == rowk.end () || kjentry->col != j ) continue; // Zero entry
00452           RealType kjval = kjentry->value;
00453           RealType newijval = - ikval * kjval;
00454           _decomp.add ( i, j, newijval );
00455           ijentry->value += newijval;
00456         }
00457       }
00458     }
00459   }
00460 
00461   virtual ~ILU0Preconditioner () {}
00462 
00463   virtual void applyAdd ( const Vector<RealType>& arg, Vector<RealType>& dest ) const {
00464     Vector<RealType> tmp ( dest.size () );
00465     apply ( arg, tmp );
00466     dest += tmp;
00467   }
00468 
00469   virtual void apply ( const Vector<RealType> &arg, Vector<RealType> &dest ) const {
00470     vector<typename Row<RealType>::RowEntry > row;
00471 
00472     // Forwards-substitution
00473     for ( int i = 0; i < _decomp.getNumRows (); ++i ) {
00474       _decomp.makeSortedRowEntries ( row, i ); // Sorted to allow early break
00475       typename vector<typename Row<RealType>::RowEntry >::iterator end = row.end ();
00476       RealType sum = 0;
00477 
00478       for ( typename vector<typename Row<RealType>::RowEntry >::iterator kentry = row.begin (); kentry != end && kentry->col < i; ++kentry ) {
00479         int k = kentry->col;
00480 
00481         sum += dest [k] * kentry->value;
00482       }
00483       dest [i] = arg [i] - sum;
00484     }
00485 
00486     // Backwards-substitution
00487     for ( int i = _decomp.getNumRows () - 1; i >= 0; --i ) {
00488       _decomp.makeRowEntries ( row, i ); // Need not be sorted
00489       RealType sum = 0;
00490       RealType diag = 0;
00491 
00492       // Strangely this crashes with the intel-CC if using a normal iterator loop (like above)
00493       for ( int kentry = 0, end = row.size (); kentry < end; kentry++ ) {
00494         const int k = row[kentry].col;
00495 
00496         if ( k < i ) continue; // U part
00497         if ( k == i ) {
00498           diag = row [kentry].value;
00499           continue;
00500         }
00501         sum += dest [k] * row[kentry].value;
00502       }
00503       dest [i] -= sum; dest [i] /= diag;
00504     }
00505   }
00506 };
00507 
00508 template <typename RealType, typename BlockMatrixType>
00509 class ILU0BlockPreconditioner : public Op<aol::MultiVector<RealType> > {
00510 
00511 protected:
00512   BlockMatrixType _decomp;
00513 
00514 public:
00515   ILU0BlockPreconditioner ( const BlockMatrixType& matrix )
00516       : _decomp ( matrix ) {
00517     vector<typename Row<RealType>::RowEntry > rowi;
00518     vector<typename Row<RealType>::RowEntry > rowk;
00519 
00520     for ( int i = 0; i < _decomp.getNumRows (); ++i ) { // Nothing is done in first row
00521       for ( int locRow = 0; locRow < _decomp.getReference(0, 0).getNumRows(); locRow++ ) {
00522         if ( !(i==0 && locRow==0) ) {
00523           _decomp.makeUnblockedSortedRowEntries( rowi, i, locRow );
00524 
00525           for ( typename vector<typename Row<RealType>::RowEntry >::iterator ikentry = rowi.begin ();
00526                 ikentry != rowi.end () && ikentry->col < i*_decomp.getReference(0, 0).getNumRows()+locRow; ++ikentry ) { // Only left from diagonal, RowEntries are sorted
00527             int k = ikentry->col;
00528 
00529             RealType ikval = ikentry->value;
00530             _decomp.makeUnblockedSortedRowEntries( rowk, floor(k/_decomp.getReference(0, 0).getNumRows()), k%_decomp.getReference(0, 0).getNumRows() );
00531 
00532             typename vector<typename Row<RealType>::RowEntry >::iterator kjentry = rowk.begin ();
00533 
00534             // Find index kk
00535             while ( kjentry != rowk.end () && kjentry->col < k ) ++kjentry;
00536             if ( kjentry == rowk.end () || kjentry->col != k ) throw Exception ( "Pivot zero in ILU", __FILE__, __LINE__ );
00537 
00538             // Eliminate index ik, store factor
00539             ikval /= kjentry->value;
00540 
00541             _decomp.getReference(i, floor(k/_decomp.getReference(0, 0).getNumCols())).set ( locRow, k%_decomp.getReference(0, 0).getNumCols(), ikval );
00542 
00543             typename vector<typename Row<RealType>::RowEntry >::iterator ikentrynext = ikentry; ++ikentrynext; // Only right from current column
00544             for ( typename vector<typename Row<RealType>::RowEntry >::iterator ijentry = ikentrynext; ijentry != rowi.end (); ++ijentry ) {
00545               int j = ijentry->col;
00546 
00547               // Find index kj
00548               while ( kjentry != rowk.end () && kjentry->col < j ) ++kjentry;
00549               if ( kjentry == rowk.end () || kjentry->col != j ) continue; // Zero entry
00550               RealType kjval = kjentry->value;
00551               RealType newijval = - ikval * kjval;
00552               _decomp.getReference(i, floor(j/_decomp.getReference(0, 0).getNumCols())).add ( locRow, j%_decomp.getReference(0, 0).getNumCols(), newijval );
00553               ijentry->value += newijval;
00554             }
00555           }
00556         }
00557       }
00558     }
00559 //   _decomp.getPointer(0,0)->print(cerr);
00560 // _decomp.getPointer(0,1)->print(cerr);
00561 // _decomp.getPointer(1,0)->print(cerr);
00562 // _decomp.getPointer(1,1)->print(cerr);
00563   }
00564 
00565   virtual ~ILU0BlockPreconditioner () {}
00566 
00567   virtual void applyAdd ( const MultiVector<RealType>& arg, MultiVector<RealType>& dest ) const {
00568     MultiVector<RealType> tmp ( dest, STRUCT_COPY );
00569     apply ( arg, tmp );
00570     dest += tmp;
00571   }
00572 
00573   virtual void apply ( const MultiVector<RealType> &arg, MultiVector<RealType> &dest ) const {
00574     vector<typename Row<RealType>::RowEntry > row;
00575 
00576     // Forwards-substitution
00577     for ( int i = 0; i < _decomp.getNumRows (); ++i ) { 
00578       for ( int locRow = 0; locRow < _decomp.getReference(0, 0).getNumRows(); locRow++ ) {
00579         _decomp.makeUnblockedSortedRowEntries( row, i, locRow );
00580 
00581         typename vector<typename Row<RealType>::RowEntry >::iterator end = row.end ();
00582         RealType sum = 0;
00583 
00584         for ( typename vector<typename Row<RealType>::RowEntry >::iterator kentry = row.begin (); kentry != end && kentry->col < i*_decomp.getReference(0, 0).getNumRows()+locRow; ++kentry ) {
00585           int k = kentry->col;
00586 
00587           sum += dest [floor(k/_decomp.getReference(0, 0).getNumRows())][k%_decomp.getReference(0, 0).getNumRows()] * kentry->value;
00588         }
00589         dest [i][locRow] = arg [i][locRow] - sum;
00590       }
00591     }
00592 
00593     // Backwards-substitution
00594     for ( int i = _decomp.getNumRows () - 1; i >= 0; --i ) {
00595       for ( int locRow = _decomp.getReference(0, 0).getNumRows() - 1; locRow >= 0; --locRow ) {
00596         _decomp.makeUnblockedRowEntries ( row, i, locRow ); // Need not be sorted
00597         RealType sum = 0;
00598         RealType diag = 0;
00599 
00600         // Strangely this crashes with the intel-CC if using a normal iterator loop (like above)
00601         for ( int kentry = 0, end = row.size (); kentry < end; kentry++ ) {
00602           const int k = row[kentry].col;
00603 
00604           if ( k < i*_decomp.getReference(0, 0).getNumRows()+locRow ) continue; // U part
00605           if ( k == i*_decomp.getReference(0, 0).getNumRows()+locRow ) {
00606             diag = row [kentry].value;
00607             continue;
00608           }
00609           sum += dest [floor(k/_decomp.getReference(0, 0).getNumRows())][k%_decomp.getReference(0, 0).getNumRows()] * row[kentry].value;
00610         }
00611         dest [i][locRow] -= sum; 
00612         dest [i][locRow] /= diag;
00613       }
00614     }
00615   }
00616 };
00617 
00618 
00619 /** \brief A class for block-diagonally preconditioning BlockMatrices with diagonal blocks all the same size and all set.
00620  *  \author Schwen
00621  *  \ingroup solver
00622  *  This preconditioner emulates rearranging the unknowns such that x, y and z component of function at one grid point are 
00623  *  subsequent unknowns in the system and then inverts d x d (d = 2, 3) blocks along the diagonal, then the rearranging is undone.
00624  */
00625 
00626 template< typename RealType, qc::Dimension Dim >
00627 class BlockDiagonalPreconditioner : public Op< aol::MultiVector< RealType > > {
00628 
00629 protected:
00630   aol::BlockMatrix< aol::DiagonalMatrix<RealType> > _blockDiagBlockInv;
00631 
00632 public:
00633   template< typename BlockMatrixType >
00634   explicit BlockDiagonalPreconditioner ( const BlockMatrixType &BlockMat ) :
00635     _blockDiagBlockInv ( Dim, Dim, BlockMat.getReference ( 0, 0 ).getNumRows(), BlockMat.getReference ( 0, 0 ).getNumCols() ) {
00636 
00637     setBlockDiagBlockInv ( BlockMat, _blockDiagBlockInv );
00638   }
00639 
00640 
00641   virtual void applyAdd ( const MultiVector<RealType> &Arg, MultiVector<RealType> &Dest ) const {
00642     _blockDiagBlockInv.applyAdd ( Arg, Dest );
00643   }
00644 
00645   virtual void apply ( const MultiVector<RealType> &Arg, MultiVector<RealType> &Dest ) const {
00646     _blockDiagBlockInv.apply ( Arg, Dest );
00647   }
00648 
00649   void getPreconditionerBlockMatrix ( aol::BlockMatrix < aol::DiagonalMatrix< RealType > > &BlockMat ) const {
00650     if ( BlockMat.getNumCols() != Dim || BlockMat.getNumRows() != Dim )
00651       throw aol::Exception ( "aol::BlockDiagonalPreconditioner: illegal number of blocks", __FILE__, __LINE__ );
00652 
00653     for ( int i = 0; i < Dim; ++i ) {
00654       for ( int j = 0; j < Dim; ++j ) {
00655         BlockMat.getReference(i, j) = _blockDiagBlockInv.getReference(i, j);
00656       }
00657     }
00658   }
00659 
00660   const aol::BlockMatrix< aol::DiagonalMatrix<RealType> >& getPreconditionerBlockMatrixRef ( ) const {
00661     return ( _blockDiagBlockInv );
00662   }
00663 
00664   template< typename BlockMatrixType >
00665   static void setBlockDiagBlockInv ( const BlockMatrixType &BlockMat, aol::BlockMatrix< aol::DiagonalMatrix<RealType> > &bdbInv ) {
00666     if ( BlockMat.getNumCols() != Dim || BlockMat.getNumRows() != Dim )
00667       throw aol::Exception ( "aol::BlockDiagonalPreconditioner::setBlockDiagBlockInv: illegal number of blocks", __FILE__, __LINE__ );
00668 
00669     const int n = BlockMat.getReference ( 0, 0 ).getNumRows(); // we assume squared blocks of all the same size and rely on bounds checking
00670 
00671     typename aol::MatDimTrait<RealType, Dim>::MatType diagBlock, diagBlockInv;
00672     
00673     for ( int i = 0; i < n; ++i ) {
00674       for ( short j = 0; j < Dim; ++j ) {      
00675         for ( short k = 0; k < Dim; ++k ) {
00676           diagBlock.set ( j, k, BlockMat.getReference ( j, k ).getDiag ( i ) );
00677         }
00678       }
00679       
00680       if ( diagBlock.det() != aol::NumberTrait<RealType>::zero ) {
00681         diagBlockInv = diagBlock.inverse();
00682       } else {
00683         diagBlockInv.setIdentity();
00684       }
00685       
00686       for ( short j = 0; j < Dim; ++j ) {      
00687         for ( short k = 0; k < Dim; ++k ) {
00688           bdbInv.getReference ( j, k ).set( i, i, diagBlockInv.get ( j, k ) );
00689         }
00690       }
00691     }
00692 
00693   }
00694 
00695 };
00696 
00697 /** Block-SSOR preconditioning which is naively implemented by
00698  *  performing one forward and one backward relaxed Block-Gauss-Seidel
00699  *  sweep. This does not have optimal complexity.
00700  *  Currently, it only works for 3D (due to the BlockGaussSeidelSweeper).
00701  *  Can be parallelized by using parallelized Gauss-Seidel iterations
00702  *  (which is no longer SSOR!)
00703  *  \author Schwen
00704  *  \ingroup solver
00705  */
00706 template< typename VectorType, typename BlockOpType >
00707 class BlockGaussSeidelPreconditioner : public Op< VectorType > {
00708 protected:
00709   typedef typename VectorType::DataType DataType;
00710   aol::BlockGaussSeidelSweeper< DataType, BlockOpType > _sweeper;
00711   aol::GaussSeidelSweepingMode _firstSweepMode, _secondSweepMode;
00712   
00713 public:
00714   explicit BlockGaussSeidelPreconditioner ( const BlockOpType &BlockMatrix, const DataType Relax = 1.2 ) : _sweeper ( BlockMatrix, Relax ), _firstSweepMode ( aol::GAUSS_SEIDEL_FORWARD ), _secondSweepMode ( aol::GAUSS_SEIDEL_BACKWARD ) {
00715   }
00716 
00717   void setSweepingModes ( const aol::GaussSeidelSweepingMode firstMode, const aol::GaussSeidelSweepingMode secondMode ) {
00718     _firstSweepMode = firstMode;
00719     _secondSweepMode = secondMode;
00720   }
00721 
00722   virtual void applyAdd ( const VectorType &Arg, VectorType &Dest ) const {
00723     VectorType tmp ( Dest, aol::DEEP_COPY );
00724     // ? VectorType tmp ( Dest, aol::STRUCT_COPY );
00725     apply ( Arg, tmp );
00726     Dest += tmp;
00727   }
00728 
00729   virtual void apply ( const VectorType &Arg, VectorType &Dest ) const {
00730     // ? Dest.setZero()
00731     _sweeper.apply ( Arg, Dest, _firstSweepMode );
00732     _sweeper.apply ( Arg, Dest, _secondSweepMode );
00733   }
00734 };
00735 
00736 
00737 /** Block variant of the SSOR preconditioner.
00738  *  \author Schwen (based on code by Droske and Wirth)
00739  *  \ingroup solver
00740  */
00741 template< typename RealType, typename BlockMatrixType, qc::Dimension Dim >
00742 class BlockSSORPreconditioner : public Op< aol::MultiVector< RealType > > {
00743 protected:
00744   std::vector< typename aol::MatDimTrait<RealType,Dim>::MatType > _invDiags; // for cacheing inverses of the diagonal blocks
00745   const BlockMatrixType &_bMat;
00746   RealType _omega;
00747 
00748 public:
00749   BlockSSORPreconditioner ( const BlockMatrixType &BlockMatrix, const RealType Omega = 1.2f )
00750     : _invDiags ( ), _bMat ( BlockMatrix ), _omega ( Omega ) {
00751     
00752     cacheDiagInverses ( _bMat, _invDiags );
00753  
00754   }
00755 
00756   virtual void applyAdd ( const aol::MultiVector<RealType> &Arg, aol::MultiVector<RealType> &Dest ) const {
00757     aol::MultiVector<RealType> tmp ( Dest, aol::STRUCT_COPY );
00758     apply ( Arg, tmp );
00759     Dest += tmp;
00760   }
00761 
00762   virtual void apply ( const aol::MultiVector<RealType> &Arg, aol::MultiVector<RealType> &Dest ) const {
00763 
00764     for ( int i = 0; i < Dest[0].size(); ++i ) { // assuming all have the same size
00765       aol::Vec<Dim,RealType> v;
00766       for ( short a = 0; a < Dim; ++a ) {
00767         for ( short b = 0; b < Dim; ++b ) {
00768           std::vector< typename aol::Row<RealType>::RowEntry > vec;
00769           this->_bMat.getReference ( a, b ).makeSortedRowEntries ( vec, i );
00770           
00771           for ( typename std::vector<typename aol::Row<RealType>::RowEntry >::iterator it = vec.begin(); ( it != vec.end() ) && ( it->col < i ); ++it ) {
00772             v[ a ] += it->value * Dest[ b ][ it->col ];
00773           }
00774         }
00775       }
00776       aol::Vec<Dim,RealType> tmp, tmpprod;
00777       for ( short a = 0; a < Dim; ++a ) 
00778         tmp[a] = Arg[a][i];
00779       tmp -= _omega * v;
00780       _invDiags[i].mult ( tmp, tmpprod );
00781       for ( short a = 0; a < Dim; ++a ) 
00782         Dest[a][i] = tmpprod[a];
00783     }
00784 
00785     for ( int i = Dest[0].size() - 1; i >= 0; --i ) {
00786       aol::Vec<Dim,RealType> v;
00787       for ( short a = 0; a < Dim; ++a ) {
00788         for ( short b = 0; b < Dim; ++b ) {
00789           std::vector< typename aol::Row<RealType>::RowEntry > vec;
00790           this->_bMat.getReference ( a, b ).makeRowEntries ( vec, i ); // sortedness not necessary here
00791           
00792           for ( typename std::vector<typename aol::Row<RealType>::RowEntry >::iterator it = vec.begin(); it != vec.end(); ++it ) {
00793             if ( it->col > i ) {
00794               v[ a ] += it->value * Dest[ b ][ it->col ];
00795             }
00796           }
00797         }
00798       }
00799       v *= _omega;
00800       aol::Vec<Dim,RealType> tmpmlt;
00801       _invDiags[i].mult ( v, tmpmlt );
00802       for ( short a = 0; a < Dim; ++a ) 
00803         Dest[a][i] -= tmpmlt[a];
00804     }
00805   }
00806 
00807   static void cacheDiagInverses ( const BlockMatrixType &BlockMat, std::vector< typename aol::MatDimTrait<RealType,Dim>::MatType > &invDiags ) {
00808     const int n = BlockMat.getReference(0,0).getNumRows(); // assuming this is correct ...
00809     invDiags.resize ( n );
00810     
00811     for ( int row = 0; row < n; ++row ) {
00812       typename aol::MatDimTrait<RealType,Dim>::MatType diag;
00813       
00814       for ( short a = 0; a < Dim; ++a ) {
00815         for ( short b = 0; b < Dim; ++b ) {
00816           std::vector< typename aol::Row<RealType>::RowEntry > vec;
00817           BlockMat.getReference ( a, b ).makeRowEntries ( vec, row );
00818           for ( typename std::vector<typename aol::Row<RealType>::RowEntry >::iterator it = vec.begin(); it != vec.end(); ++it ) {
00819             if ( it->col == row ) {
00820               diag[a][b] = it->value;                 
00821               break;
00822             }
00823           }
00824         }
00825       }
00826       
00827       if ( diag.det() != 0.0 )
00828         invDiags[row] = diag.inverse(); 
00829       else
00830         invDiags[row].setIdentity();
00831     }
00832   }
00833 };
00834 
00835 
00836 } // end namespace
00837 
00838 #endif

Generated on Fri Sep 9 2011 21:07:54 for QuocMesh by doxygen 1.7.1