diff -p --exclude=Makefile.am --exclude=Makefile.in istl/basearray.hh ../dune-istl-1.1/istl/basearray.hh *** istl/basearray.hh Wed Aug 15 16:06:05 2007 --- ../dune-istl-1.1/istl/basearray.hh Thu Apr 3 03:59:22 2008 *************** namespace Dune { *** 18,24 **** /** \brief A simple array container for objects of type B ! Implements: - iterator access - const_iterator access --- 18,24 ---- /** \brief A simple array container for objects of type B ! Implement. - iterator access - const_iterator access *************** namespace Dune { *** 78,84 **** { public: //! \brief The unqualified value type ! typedef typename Dune::RemoveConst::Type ValueType; friend class BidirectionalIteratorFacade, const ValueType>; friend class BidirectionalIteratorFacade, ValueType>; --- 78,84 ---- { public: //! \brief The unqualified value type ! typedef typename remove_const::type ValueType; friend class BidirectionalIteratorFacade, const ValueType>; friend class BidirectionalIteratorFacade, ValueType>; *************** namespace Dune { *** 553,559 **** { public: //! \brief The unqualified value type ! typedef typename Dune::RemoveConst::Type ValueType; friend class BidirectionalIteratorFacade, const ValueType>; friend class BidirectionalIteratorFacade, ValueType>; --- 553,559 ---- { public: //! \brief The unqualified value type ! typedef typename remove_const::type ValueType; friend class BidirectionalIteratorFacade, const ValueType>; friend class BidirectionalIteratorFacade, ValueType>; diff -p --exclude=Makefile.am --exclude=Makefile.in istl/bcrsmatrix.hh ../dune-istl-1.1/istl/bcrsmatrix.hh *** istl/bcrsmatrix.hh Wed Aug 15 16:06:05 2007 --- ../dune-istl-1.1/istl/bcrsmatrix.hh Thu Apr 3 03:59:22 2008 *************** *** 1,7 **** #ifndef DUNE_BCRSMATRIX_HH #define DUNE_BCRSMATRIX_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_BCRSMATRIX_HH #define DUNE_BCRSMATRIX_HH ! #include #include #include #include *************** namespace Dune { *** 141,148 **** --- 141,153 ---- B[3][3] = 8; \endcode */ + #ifdef DUNE_EXPRESSIONTEMPLATES + template + class BCRSMatrix : public ExprTmpl::Matrix< BCRSMatrix > + #else template class BCRSMatrix + #endif { private: enum BuildStage{ *************** namespace Dune { *** 245,251 **** public: //! \brief The unqualified value type ! typedef typename Dune::RemoveConst::Type ValueType; friend class RandomAccessIteratorFacade, const ValueType>; friend class RandomAccessIteratorFacade, ValueType>; --- 250,256 ---- public: //! \brief The unqualified value type ! typedef typename remove_const::type ValueType; friend class RandomAccessIteratorFacade, const ValueType>; friend class RandomAccessIteratorFacade, ValueType>; *************** namespace Dune { *** 686,700 **** r[i].setsize(s); } ! //! increment size of row i by 1 ! void incrementrowsize (size_type i) { if (build_mode!=random) DUNE_THROW(ISTLError,"requires random build mode"); if (ready) DUNE_THROW(ISTLError,"matrix row sizes already built up"); ! r[i].setsize(r[i].getsize()+1); } //! indicate that size of all rows is defined --- 691,715 ---- r[i].setsize(s); } ! //! get current number of indices in row i ! size_type getrowsize (size_type i) const ! { ! #ifdef DUNE_ISTL_WITH_CHECKING ! if (r==0) DUNE_THROW(ISTLError,"row not initialized yet"); ! if (i>=n) DUNE_THROW(ISTLError,"index out of range"); ! #endif ! return r[i].getsize(); ! } ! ! //! increment size of row i by s (1 by default) ! void incrementrowsize (size_type i, size_type s = 1) { if (build_mode!=random) DUNE_THROW(ISTLError,"requires random build mode"); if (ready) DUNE_THROW(ISTLError,"matrix row sizes already built up"); ! r[i].setsize(r[i].getsize()+s); } //! indicate that size of all rows is defined *************** namespace Dune { *** 1138,1143 **** --- 1153,1164 ---- return m; } + //! number of blocks that are stored (the number of blocks that possibly are nonzero) + size_type nonzeroes () const + { + return nnz; + } + /** \brief row dimension of block r \bug Does not count empty rows (FlySpray #7) */ *************** namespace Dune { *** 1361,1366 **** --- 1382,1406 ---- }; + #ifdef DUNE_EXPRESSIONTEMPLATES + template + struct FieldType< BCRSMatrix > + { + typedef typename FieldType::type type; + }; + + template + struct BlockType< BCRSMatrix > + { + typedef B type; + }; + template + struct RowType< BCRSMatrix > + { + typedef CompressedBlockVectorWindow type; + }; + #endif + /** @} end documentation */ } // end namespace diff -p --exclude=Makefile.am --exclude=Makefile.in istl/bvector.hh ../dune-istl-1.1/istl/bvector.hh *** istl/bvector.hh Tue Jul 10 13:35:45 2007 --- ../dune-istl-1.1/istl/bvector.hh Thu Apr 3 03:59:22 2008 *************** *** 1,13 **** #ifndef DUNE_BVECTOR_HH #define DUNE_BVECTOR_HH ! #include #include #include "istlexception.hh" #include "allocator.hh" #include "basearray.hh" /*! \file \brief This file implements a vector space as a tensor product of --- 1,17 ---- #ifndef DUNE_BVECTOR_HH #define DUNE_BVECTOR_HH ! #include #include #include "istlexception.hh" #include "allocator.hh" #include "basearray.hh" + #ifdef DUNE_EXPRESSIONTEMPLATES + #include + #endif + /*! \file \brief This file implements a vector space as a tensor product of *************** namespace Dune { *** 68,73 **** --- 72,78 ---- } //===== vector space arithmetic + #ifndef DUNE_EXPRESSIONTEMPLATES //! vector space addition block_vector_unmanaged& operator+= (const block_vector_unmanaged& y) { *************** namespace Dune { *** 101,107 **** for (int i=0; in; ++i) (*this)[i] /= k; return *this; } ! //! vector space axpy operation block_vector_unmanaged& axpy (const field_type& a, const block_vector_unmanaged& y) { --- 106,112 ---- for (int i=0; in; ++i) (*this)[i] /= k; return *this; } ! #endif //! vector space axpy operation block_vector_unmanaged& axpy (const field_type& a, const block_vector_unmanaged& y) { *************** namespace Dune { *** 128,133 **** --- 133,139 ---- //===== norms + #ifndef DUNE_EXPRESSIONTEMPLATES //! one norm (sum over absolute values of entries) double one_norm () const { *************** namespace Dune { *** 175,180 **** --- 181,187 ---- for (size_type i=0; in; ++i) max = std::max(max,(*this)[i].infinity_norm_real()); return max; } + #endif //===== sizes *************** namespace Dune { *** 214,221 **** --- 221,234 ---- Setting the compile time switch DUNE_ISTL_WITH_CHECKING enables error checking. */ + #ifdef DUNE_EXPRESSIONTEMPLATES + template + class BlockVector : public block_vector_unmanaged , + public Dune::ExprTmpl::Vector< Dune::BlockVector > + #else template class BlockVector : public block_vector_unmanaged + #endif { public: *************** namespace Dune { *** 371,376 **** --- 384,540 ---- this->n=size; } + + + + #ifdef DUNE_EXPRESSIONTEMPLATES + //! random access to blocks + B& operator[] (size_type i) + { + return block_vector_unmanaged::operator[](i); + } + + //! same for read only access + const B& operator[] (size_type i) const + { + return block_vector_unmanaged::operator[](i); + } + + //! dimension of the vector space + size_type N () const + { + return block_vector_unmanaged::N(); + } + + BlockVector (const BlockVector& a) { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Copy Constructor BlockVector\n"; + #endif + capacity_ = a.capacity_; + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + this->n = a.N(); + // copy data + assignFrom(a); + } + template + BlockVector (Dune::ExprTmpl::Expression op) { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Copy Constructor Expression\n"; + #endif + capacity_ = op.N(); + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + this->n = op.N(); + assignFrom(op); + } + template + BlockVector (const Dune::ExprTmpl::Vector & op) { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Copy Constructor Vector\n"; + #endif + capacity_ = op.N(); + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + this->n = op.N(); + assignFrom(op); + } + BlockVector& operator = (const BlockVector& a) { + if (&a!=this) // check if this and a are different objects + { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT + << "BlockVector Assignment Operator BlockVector\n"; + #endif + // adjust size of vector + if (capacity_!=a.capacity_) // check if size is different + { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT + << "BlockVector Resize to size(BlockVector)\n"; + #endif + if (capacity_>0) A::template free(this->p); // delete old memory + capacity_ = a.capacity_; + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + } + this->n = a.N(); + // copy data + assignFrom(a); + } + return assignFrom(a); + } + template + BlockVector& operator = (Dune::ExprTmpl::Expression op) { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Assignment Operator Expression\n"; + #endif + // adjust size of vector + if (capacity_ < op.N()) // check if size is different + { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Resize to size(Expression)\n"; + #endif + if (capacity_>0) A::template free(this->p); // delete old memory + capacity_ = op.N(); + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + } + this->n = op.N(); + return assignFrom(op); + } + template + BlockVector& operator = (const Dune::ExprTmpl::Vector & op) { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Assignment Operator Vector\n"; + #endif + // adjust size of vector + if (capacity_ < op.N()) // check if size is different + { + #ifdef DUNE_VVERBOSE + Dune::dvverb << INDENT << "BlockVector Resize to size(Vector)\n"; + #endif + if (capacity_>0) A::template free(this->p); // delete old memory + capacity_ = op.N(); + if (capacity_>0) + this->p = A::template malloc(capacity_); + else + { + this->p = 0; + capacity_ = 0; + } + } + this->n = op.N(); + return assignFrom(op); + } + #endif + + #ifndef DUNE_EXPRESSIONTEMPLATES //! copy constructor BlockVector (const BlockVector& a) : block_vector_unmanaged(a) *************** namespace Dune { *** 413,418 **** --- 577,583 ---- // and copy elements for (size_type i=0; in; i++) this->p[i]=a.p[i]; } + #endif //! free dynamic memory ~BlockVector () *************** namespace Dune { *** 421,426 **** --- 586,592 ---- } //! assignment + #ifndef DUNE_EXPRESSIONTEMPLATES BlockVector& operator= (const BlockVector& a) { if (&a!=this) // check if this and a are different objects *************** namespace Dune { *** 452,457 **** --- 618,624 ---- // forward to regular assignement operator return this->operator=(static_cast(a)); } + #endif //! assign from scalar BlockVector& operator= (const field_type& k) *************** namespace Dune { *** 989,994 **** --- 1156,1174 ---- } }; + #ifdef DUNE_EXPRESSIONTEMPLATES + template + struct FieldType< BlockVector > + { + typedef typename FieldType::type type; + }; + template + struct BlockType< BlockVector > + { + typedef B type; + }; + #endif + } // end namespace #endif diff -p --exclude=Makefile.am --exclude=Makefile.in istl/communicator.hh ../dune-istl-1.1/istl/communicator.hh *** istl/communicator.hh Thu May 10 13:12:23 2007 --- ../dune-istl-1.1/istl/communicator.hh Thu Apr 3 03:59:23 2008 *************** *** 1,4 **** ! // $Id: communicator.hh 776 2007-05-10 11:12:19Z mblatt $ #ifndef DUNE_COMMUNICATOR #define DUNE_COMMUNICATOR --- 1,4 ---- ! // $Id: communicator.hh 772 2007-05-08 08:56:30Z christi $ #ifndef DUNE_COMMUNICATOR #define DUNE_COMMUNICATOR *************** namespace Dune *** 153,168 **** template class FieldVector; ! template class BlockVector; template ! struct CommPolicy, A> > { ! typedef BlockVector, A> Type; ! typedef typename Type::block_type IndexedType; ! typedef SizeOne IndexedTypeFlag; static const void* getAddress(const Type& v, int i); --- 153,168 ---- template class FieldVector; ! template class VariableBlockVector; template ! struct CommPolicy, A> > { ! typedef VariableBlockVector, A> Type; ! typedef typename Type::B IndexedType; ! typedef VariableSize IndexedTypeFlag; static const void* getAddress(const Type& v, int i); *************** namespace Dune *** 468,474 **** * @param interface The interface that defines what indices are to be communicated. */ template ! typename EnableIf::IndexedTypeFlag>::value, void>::Type build(const Interface& interface); /** --- 468,474 ---- * @param interface The interface that defines what indices are to be communicated. */ template ! typename EnableIf::IndexedTypeFlag>::value, void>::Type build(const Interface& interface); /** *************** namespace Dune *** 898,912 **** } template ! inline const void* CommPolicy, A> >::getAddress(const Type& v, int index) { ! return &(v[index]); } template ! inline int CommPolicy, A> >::getSize(const Type& v, int index) { ! return 1; } --- 898,912 ---- } template ! inline const void* CommPolicy, A> >::getAddress(const Type& v, int index) { ! return &(v[index][0]); } template ! inline int CommPolicy, A> >::getSize(const Type& v, int index) { ! return v[index].getsize(); } *************** namespace Dune *** 1132,1138 **** template template ! typename EnableIf::IndexedTypeFlag>::value, void>::Type BufferedCommunicator::build(const Interface& interface) { typedef typename Interface::InformationMap::const_iterator const_iterator; --- 1132,1138 ---- template template ! typename EnableIf::IndexedTypeFlag>::value, void>::Type BufferedCommunicator::build(const Interface& interface) { typedef typename Interface::InformationMap::const_iterator const_iterator; diff -p --exclude=Makefile.am --exclude=Makefile.in istl/gsetc.hh ../dune-istl-1.1/istl/gsetc.hh *** istl/gsetc.hh Thu May 10 11:35:26 2007 --- ../dune-istl-1.1/istl/gsetc.hh Thu Apr 3 03:59:23 2008 *************** *** 1,7 **** #ifndef DUNE_GSETC_HH #define DUNE_GSETC_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_GSETC_HH #define DUNE_GSETC_HH ! #include #include #include #include diff -p --exclude=Makefile.am --exclude=Makefile.in istl/ilu.hh ../dune-istl-1.1/istl/ilu.hh *** istl/ilu.hh Thu May 10 11:35:11 2007 --- ../dune-istl-1.1/istl/ilu.hh Thu Apr 3 03:59:22 2008 *************** *** 1,7 **** #ifndef DUNE_ILU_HH #define DUNE_ILU_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_ILU_HH #define DUNE_ILU_HH ! #include #include #include #include diff -p --exclude=Makefile.am --exclude=Makefile.in istl/indexset.hh ../dune-istl-1.1/istl/indexset.hh *** istl/indexset.hh Thu May 10 13:12:23 2007 --- ../dune-istl-1.1/istl/indexset.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: indexset.hh 776 2007-05-10 11:12:19Z mblatt $ #ifndef DUNE_INDEXSET_HH #define DUNE_INDEXSET_HH --- 1,4 ---- ! // $Id: indexset.hh 771 2007-05-07 13:35:57Z mblatt $ #ifndef DUNE_INDEXSET_HH #define DUNE_INDEXSET_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/indicessyncer.hh ../dune-istl-1.1/istl/indicessyncer.hh *** istl/indicessyncer.hh Fri Mar 16 15:04:04 2007 --- ../dune-istl-1.1/istl/indicessyncer.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: indicessyncer.hh 727 2007-03-16 14:05:10Z mblatt $ #ifndef DUNE_INDICESSYNCER_HH #define DUNE_INDICESSYNCER_HH --- 1,4 ---- ! // $Id: indicessyncer.hh 726 2007-03-16 13:53:17Z mblatt $ #ifndef DUNE_INDICESSYNCER_HH #define DUNE_INDICESSYNCER_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/interface.hh ../dune-istl-1.1/istl/interface.hh *** istl/interface.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/interface.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: interface.hh 764 2007-05-03 12:39:19Z mblatt $ #ifndef DUNE_INTERFACE_HH #define DUNE_INTERFACE_HH --- 1,4 ---- ! // $Id: interface.hh 762 2007-05-01 17:01:20Z mblatt $ #ifndef DUNE_INTERFACE_HH #define DUNE_INTERFACE_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/io.hh ../dune-istl-1.1/istl/io.hh *** istl/io.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/io.hh Thu Apr 3 03:59:23 2008 *************** *** 1,7 **** #ifndef DUNE_ISTLIO_HH #define DUNE_ISTLIO_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_ISTLIO_HH #define DUNE_ISTLIO_HH ! #include #include #include #include *************** namespace Dune { *** 302,307 **** --- 302,308 ---- // reset the output format s.flags(oldflags); + s.precision(oldprec); } /** \brief Helper method for the writeMatrixToMatlab routine. diff -p --exclude=Makefile.am --exclude=Makefile.in istl/mpitraits.hh ../dune-istl-1.1/istl/mpitraits.hh *** istl/mpitraits.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/mpitraits.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: mpitraits.hh 764 2007-05-03 12:39:19Z mblatt $ #ifndef DUNE_MPITRAITS_HH #define DUNE_MPITRAITS_HH --- 1,4 ---- ! // $Id: mpitraits.hh 761 2007-05-01 10:48:34Z mblatt $ #ifndef DUNE_MPITRAITS_HH #define DUNE_MPITRAITS_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/operators.hh ../dune-istl-1.1/istl/operators.hh *** istl/operators.hh Thu May 10 11:35:30 2007 --- ../dune-istl-1.1/istl/operators.hh Thu Apr 3 03:59:23 2008 *************** *** 1,7 **** #ifndef DUNE_ISTLOPERATORS_HH #define DUNE_ISTLOPERATORS_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_ISTLOPERATORS_HH #define DUNE_ISTLOPERATORS_HH ! #include #include #include #include diff -p --exclude=Makefile.am --exclude=Makefile.in istl/owneroverlapcopy.hh ../dune-istl-1.1/istl/owneroverlapcopy.hh *** istl/owneroverlapcopy.hh Wed Aug 15 16:06:46 2007 --- ../dune-istl-1.1/istl/owneroverlapcopy.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: owneroverlapcopy.hh 803 2007-08-15 14:06:45Z mblatt $ #ifndef DUNE_OWNEROVERLAPCOPY_HH #define DUNE_OWNEROVERLAPCOPY_HH --- 1,4 ---- ! // $Id: owneroverlapcopy.hh 783 2007-08-03 08:17:56Z sander $ #ifndef DUNE_OWNEROVERLAPCOPY_HH #define DUNE_OWNEROVERLAPCOPY_HH *************** *** 9,15 **** #include #include ! #include"math.h" // MPI header #if HAVE_MPI --- 9,15 ---- #include #include ! #include"cmath" // MPI header #if HAVE_MPI *************** *** 17,23 **** #endif ! #include #include #if HAVE_MPI --- 17,23 ---- #endif ! #include #include #if HAVE_MPI *************** namespace Dune { *** 81,94 **** * The triple consists of the global index and the local * index and an attribute */ ! typedef tripel IndexTripel; /** * @brief A triple describing a remote index. * * The triple consists of a process number and the global index and * the attribute of the index at the remote process. */ ! typedef tripel RemoteIndexTripel; /** * @brief Add a new index triple to the set of local indices. --- 81,94 ---- * The triple consists of the global index and the local * index and an attribute */ ! typedef Tuple IndexTripel; /** * @brief A triple describing a remote index. * * The triple consists of a process number and the global index and * the attribute of the index at the remote process. */ ! typedef Tuple RemoteIndexTripel; /** * @brief Add a new index triple to the set of local indices. *************** namespace Dune { *** 97,105 **** */ void addLocalIndex (const IndexTripel& x) { ! if (x.third!=OwnerOverlapCopyAttributeSet::owner && ! x.third!=OwnerOverlapCopyAttributeSet::overlap && ! x.third!=OwnerOverlapCopyAttributeSet::copy) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); localindices.insert(x); } --- 97,105 ---- */ void addLocalIndex (const IndexTripel& x) { ! if (Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::owner && ! Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::overlap && ! Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::copy) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); localindices.insert(x); } *************** namespace Dune { *** 111,119 **** */ void addRemoteIndex (const RemoteIndexTripel& x) { ! if (x.third!=OwnerOverlapCopyAttributeSet::owner && ! x.third!=OwnerOverlapCopyAttributeSet::overlap && ! x.third!=OwnerOverlapCopyAttributeSet::copy) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); remoteindices.insert(x); } --- 111,119 ---- */ void addRemoteIndex (const RemoteIndexTripel& x) { ! if (Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::owner && ! Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::overlap && ! Element<2>::get(x)!=OwnerOverlapCopyAttributeSet::copy) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); remoteindices.insert(x); } *************** namespace Dune { *** 438,450 **** pis.beginResize(); for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i) { ! if (i->third==OwnerOverlapCopyAttributeSet::owner) ! pis.add(i->first,LI(i->second,OwnerOverlapCopyAttributeSet::owner,true)); ! if (i->third==OwnerOverlapCopyAttributeSet::overlap) ! pis.add(i->first,LI(i->second,OwnerOverlapCopyAttributeSet::overlap,true)); ! if (i->third==OwnerOverlapCopyAttributeSet::copy) ! pis.add(i->first,LI(i->second,OwnerOverlapCopyAttributeSet::copy,true)); ! // std::cout << cc.rank() << ": adding index " << i->first << " " << i->second << " " << i->third << std::endl; } pis.endResize(); --- 438,450 ---- pis.beginResize(); for (localindex_iterator i=indexinfo.localIndices().begin(); i!=indexinfo.localIndices().end(); ++i) { ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::owner) ! pis.add(Element<0>::get(*i),LI(Element<1>::get(*i),OwnerOverlapCopyAttributeSet::owner,true)); ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::overlap) ! pis.add(Element<0>::get(*i),LI(Element<1>::get(*i),OwnerOverlapCopyAttributeSet::overlap,true)); ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::copy) ! pis.add(Element<0>::get(*i),LI(Element<1>::get(*i),OwnerOverlapCopyAttributeSet::copy,true)); ! // std::cout << cc.rank() << ": adding index " << Element<0>::get(*i) << " " << Element<1>::get(*i) << " " << Element<2>::get(*i) << std::endl; } pis.endResize(); *************** namespace Dune { *** 454,485 **** if (indexinfo.remoteIndices().size()>0) { remoteindex_iterator i=indexinfo.remoteIndices().begin(); ! int p = i->first; RILM modifier = ri.template getModifier(p); typename PIS::const_iterator pi=pis.begin(); for ( ; i!=indexinfo.remoteIndices().end(); ++i) { // handle processor change ! if (p!=i->first) { ! p = i->first; modifier = ri.template getModifier(p); pi=pis.begin(); } // position to correct entry in parallel index set ! while (pi->global()!=i->second && pi!=pis.end()) ++pi; if (pi==pis.end()) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); // insert entry ! // std::cout << cc.rank() << ": adding remote index " << i->first << " " << i->second << " " << i->third << std::endl; ! if (i->third==OwnerOverlapCopyAttributeSet::owner) modifier.insert(RX(OwnerOverlapCopyAttributeSet::owner,&(*pi))); ! if (i->third==OwnerOverlapCopyAttributeSet::overlap) modifier.insert(RX(OwnerOverlapCopyAttributeSet::overlap,&(*pi))); ! if (i->third==OwnerOverlapCopyAttributeSet::copy) modifier.insert(RX(OwnerOverlapCopyAttributeSet::copy,&(*pi))); } } --- 454,485 ---- if (indexinfo.remoteIndices().size()>0) { remoteindex_iterator i=indexinfo.remoteIndices().begin(); ! int p = Element<0>::get(*i); RILM modifier = ri.template getModifier(p); typename PIS::const_iterator pi=pis.begin(); for ( ; i!=indexinfo.remoteIndices().end(); ++i) { // handle processor change ! if (p!=Element<0>::get(*i)) { ! p = Element<0>::get(*i); modifier = ri.template getModifier(p); pi=pis.begin(); } // position to correct entry in parallel index set ! while (pi->global()!=Element<1>::get(*i) && pi!=pis.end()) ++pi; if (pi==pis.end()) DUNE_THROW(ISTLError,"OwnerOverlapCopyCommunication: global index not in index set"); // insert entry ! // std::cout << cc.rank() << ": adding remote index " << Element<0>::get(*i) << " " << Element<1>::get(*i) << " " << Element<2>::get(*i) << std::endl; ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::owner) modifier.insert(RX(OwnerOverlapCopyAttributeSet::owner,&(*pi))); ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::overlap) modifier.insert(RX(OwnerOverlapCopyAttributeSet::overlap,&(*pi))); ! if (Element<2>::get(*i)==OwnerOverlapCopyAttributeSet::copy) modifier.insert(RX(OwnerOverlapCopyAttributeSet::copy,&(*pi))); } } Common subdirectories: istl/paamg and ../dune-istl-1.1/istl/paamg Only in ../dune-istl-1.1/istl/: pardiso.hh diff -p --exclude=Makefile.am --exclude=Makefile.in istl/plocalindex.hh ../dune-istl-1.1/istl/plocalindex.hh *** istl/plocalindex.hh Thu May 10 13:12:23 2007 --- ../dune-istl-1.1/istl/plocalindex.hh Thu Apr 3 03:59:23 2008 *************** *** 1,4 **** ! // $Id: plocalindex.hh 776 2007-05-10 11:12:19Z mblatt $ #ifndef DUNE_PLOCALINDEX_HH #define DUNE_PLOCALINDEX_HH --- 1,4 ---- ! // $Id: plocalindex.hh 775 2007-05-08 16:39:05Z sander $ #ifndef DUNE_PLOCALINDEX_HH #define DUNE_PLOCALINDEX_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/preconditioners.hh ../dune-istl-1.1/istl/preconditioners.hh *** istl/preconditioners.hh Thu May 10 11:35:27 2007 --- ../dune-istl-1.1/istl/preconditioners.hh Thu Apr 3 03:59:23 2008 *************** *** 1,7 **** #ifndef DUNE_PRECONDITIONERS_HH #define DUNE_PRECONDITIONERS_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_PRECONDITIONERS_HH #define DUNE_PRECONDITIONERS_HH ! #include #include #include #include *************** namespace Dune { *** 78,86 **** A solver solves a linear operator equation A(x)=b by applying one or several steps of the preconditioner. The method pre() ! is called before before the first apply operation. ! x and b are right hand side and solution vector of the linear ! system. It may. e.g., scale the system, allocate memory or compute a (I)LU decomposition. Note: The ILU decomposition could also be computed in the constructor or with a separate method of the derived method if several --- 78,86 ---- A solver solves a linear operator equation A(x)=b by applying one or several steps of the preconditioner. The method pre() ! is called before the first apply operation. ! b and x are right hand side and solution vector of the linear ! system respectively. It may. e.g., scale the system, allocate memory or compute a (I)LU decomposition. Note: The ILU decomposition could also be computed in the constructor or with a separate method of the derived method if several *************** namespace Dune { *** 248,256 **** */ virtual void apply (X& v, const Y& d) { ! for (int i=0; i<_n; i++){ ! bsorf(_A_,v,d,_w,BL()); ! } } /*! --- 248,275 ---- */ virtual void apply (X& v, const Y& d) { ! this->template apply(v,d); ! } ! ! /*! ! \brief Apply the preconditioner in a special direction. ! ! The template parameter forward indications the direction ! the smoother is applied. If true The application is ! started at the lowest index in the vector v, if false at ! the highest index of vector v. ! */ ! template ! void apply(X& v, const Y& d) ! { ! if(forward) ! for (int i=0; i<_n; i++){ ! bsorf(_A_,v,d,_w,BL()); ! } ! else ! for (int i=0; i<_n; i++){ ! bsorb(_A_,v,d,_w,BL()); ! } } /*! diff -p --exclude=Makefile.am --exclude=Makefile.in istl/remoteindices.hh ../dune-istl-1.1/istl/remoteindices.hh *** istl/remoteindices.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/remoteindices.hh Thu Apr 3 03:59:23 2008 *************** *** 1,4 **** ! // $Id: remoteindices.hh 764 2007-05-03 12:39:19Z mblatt $ #ifndef DUNE_REMOTEINDICES_HH #define DUNE_REMOTEINDICES_HH --- 1,4 ---- ! // $Id: remoteindices.hh 762 2007-05-01 17:01:20Z mblatt $ #ifndef DUNE_REMOTEINDICES_HH #define DUNE_REMOTEINDICES_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/scalarproducts.hh ../dune-istl-1.1/istl/scalarproducts.hh *** istl/scalarproducts.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/scalarproducts.hh Thu Apr 3 03:59:22 2008 *************** *** 1,7 **** #ifndef DUNE_SCALARPRODUCTS_HH #define DUNE_SCALARPRODUCTS_HH ! #include #include #include #include --- 1,7 ---- #ifndef DUNE_SCALARPRODUCTS_HH #define DUNE_SCALARPRODUCTS_HH ! #include #include #include #include *************** sequential case are provided. *** 34,40 **** Krylov space methods need to compute scalar products and norms (for convergence test only). These methods have to know about the ! underlying data decomposition. For the sequantial case adefault implementation is provided. */ template --- 34,40 ---- Krylov space methods need to compute scalar products and norms (for convergence test only). These methods have to know about the ! underlying data decomposition. For the sequential case a default implementation is provided. */ template diff -p --exclude=Makefile.am --exclude=Makefile.in istl/schwarz.hh ../dune-istl-1.1/istl/schwarz.hh *** istl/schwarz.hh Thu May 10 11:35:06 2007 --- ../dune-istl-1.1/istl/schwarz.hh Thu Apr 3 03:59:22 2008 *************** *** 6,13 **** #include // STL vector class #include ! #include // Yes, we do some math here ! #include // There is nothing better than sprintf #include // for timing measurements #include --- 6,12 ---- #include // STL vector class #include ! #include // Yes, we do some math here #include // for timing measurements #include *************** namespace Dune { *** 273,279 **** namespace Amg { template class ConstructionTraits; ! }; /** * @brief Block parallel preconditioner. --- 272,278 ---- namespace Amg { template class ConstructionTraits; ! } /** * @brief Block parallel preconditioner. diff -p --exclude=Makefile.am --exclude=Makefile.in istl/selection.hh ../dune-istl-1.1/istl/selection.hh *** istl/selection.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/selection.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: selection.hh 764 2007-05-03 12:39:19Z mblatt $ #ifndef DUNE_SELECTION_HH #define DUNE_SELECTION_HH --- 1,4 ---- ! // $Id: selection.hh 761 2007-05-01 10:48:34Z mblatt $ #ifndef DUNE_SELECTION_HH #define DUNE_SELECTION_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/solvercategory.hh ../dune-istl-1.1/istl/solvercategory.hh *** istl/solvercategory.hh Thu May 10 11:31:21 2007 --- ../dune-istl-1.1/istl/solvercategory.hh Thu Apr 3 03:59:22 2008 *************** *** 1,4 **** ! // $Id: solvercategory.hh 764 2007-05-03 12:39:19Z mblatt $ #ifndef DUNE_SOLVERCATEGORY_HH #define DUNE_SOLVERCATEGORY_HH --- 1,4 ---- ! // $Id: solvercategory.hh 751 2007-04-20 13:59:25Z mblatt $ #ifndef DUNE_SOLVERCATEGORY_HH #define DUNE_SOLVERCATEGORY_HH diff -p --exclude=Makefile.am --exclude=Makefile.in istl/solvers.hh ../dune-istl-1.1/istl/solvers.hh *** istl/solvers.hh Thu May 10 11:35:19 2007 --- ../dune-istl-1.1/istl/solvers.hh Thu Apr 3 03:59:23 2008 *************** *** 1,12 **** #ifndef DUNE_SOLVERS_HH #define DUNE_SOLVERS_HH ! #include #include #include #include #include - #include // there is nothing better than printf #include "istlexception.hh" #include "operators.hh" --- 1,11 ---- #ifndef DUNE_SOLVERS_HH #define DUNE_SOLVERS_HH ! #include #include #include #include #include #include "istlexception.hh" #include "operators.hh" *************** namespace Dune { *** 20,26 **** @ingroup ISTL */ /** @addtogroup ISTL_Solvers ! @{ */ --- 19,25 ---- @ingroup ISTL */ /** @addtogroup ISTL_Solvers ! @{ */ *************** namespace Dune { *** 43,77 **** struct InverseOperatorResult { /** \brief Default constructor */ ! InverseOperatorResult () ! { ! clear(); ! } /** \brief Resets all data */ ! void clear () ! { ! iterations = 0; ! reduction = 0; ! converged = false; ! conv_rate = 1; ! elapsed = 0; ! } /** \brief Number of iterations */ ! int iterations; /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */ ! double reduction; /** \brief True if convergence criterion has been met */ ! bool converged; /** \brief Convergence rate (average reduction per step) */ ! double conv_rate; /** \brief Elapsed time in seconds */ ! double elapsed; }; --- 42,76 ---- struct InverseOperatorResult { /** \brief Default constructor */ ! InverseOperatorResult () ! { ! clear(); ! } /** \brief Resets all data */ ! void clear () ! { ! iterations = 0; ! reduction = 0; ! converged = false; ! conv_rate = 1; ! elapsed = 0; ! } /** \brief Number of iterations */ ! int iterations; /** \brief Reduction achieved: \f$ \|b-A(x^n)\|/\|b-A(x^0)\|\f$ */ ! double reduction; /** \brief True if convergence criterion has been met */ ! bool converged; /** \brief Convergence rate (average reduction per step) */ ! double conv_rate; /** \brief Elapsed time in seconds */ ! double elapsed; }; *************** namespace Dune { *** 90,112 **** template class InverseOperator { public: ! //! \brief Type of the domain of the operator we are the inverse of. typedef X domain_type; ! //! \brief Type of the range of the operator we are the inverse of. typedef Y range_type; /** \brief The field type of the operator. */ typedef typename X::field_type field_type; /** ! \brief Apply inverse operator, ! ! \warning Note: right hand side b may be overwritten! ! ! \param x The left hand side to store the result in. ! \param b The right hand side ! \param res Object to store the statistics about applying the operator. */ virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0; --- 89,111 ---- template class InverseOperator { public: ! //! \brief Type of the domain of the operator to be inverted. typedef X domain_type; ! //! \brief Type of the range of the operator to be inverted. typedef Y range_type; /** \brief The field type of the operator. */ typedef typename X::field_type field_type; /** ! \brief Apply inverse operator, ! ! \warning Note: right hand side b may be overwritten! ! ! \param x The left hand side to store the result in. ! \param b The right hand side ! \param res Object to store the statistics about applying the operator. */ virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0; *************** namespace Dune { *** 124,129 **** --- 123,163 ---- //! \brief Destructor virtual ~InverseOperator () {} + + protected: + // spacing values + enum { iterationSpacing = 5 , normSpacing = 16 }; + + //! helper function for printing header of solver output + void printHeader(std::ostream& s) const + { + s << std::setw(iterationSpacing) << " Iter"; + s << std::setw(normSpacing) << "Defect"; + s << std::setw(normSpacing) << "Rate" << std::endl; + } + + //! helper function for printing solver output + template + void printOutput(std::ostream& s, + const int iter, + const DataType& norm, + const DataType& norm_old) const + { + const DataType rate = norm/norm_old; + s << std::setw(iterationSpacing) << iter << " "; + s << std::setw(normSpacing) << norm << " "; + s << std::setw(normSpacing) << rate << std::endl; + } + + //! helper function for printing solver output + template + void printOutput(std::ostream& s, + const int iter, + const DataType& norm) const + { + s << std::setw(iterationSpacing) << iter << " "; + s << std::setw(normSpacing) << norm << std::endl; + } }; *************** namespace Dune { *** 142,152 **** template class LoopSolver : public InverseOperator { public: ! //! \brief The domain type of the operator we are the inverse for. typedef X domain_type; ! //! \brief The range type of the operator we are the inverse for. typedef X range_type; ! //! \brief The field type of the operator we are the inverse for. typedef typename X::field_type field_type; /*! --- 176,186 ---- template class LoopSolver : public InverseOperator { public: ! //! \brief The domain type of the operator that we do the inverse for. typedef X domain_type; ! //! \brief The range type of the operator that we do the inverse for. typedef X range_type; ! //! \brief The field type of the operator that we do the inverse for. typedef typename X::field_type field_type; /*! *************** namespace Dune { *** 170,176 **** */ template LoopSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 204,210 ---- */ template LoopSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 178,205 **** } /** ! \brief Set up loop solver ! ! \param op The operator we solve. ! \param sp The scalar product to use, e. g. SeqScalarproduct. ! \param prec The preconditioner to apply in each iteration of the loop. ! Has to inherit from Preconditioner. ! \param reduction The relative defect reduction to achieve when applying ! the operator. ! \param maxit The maximum number of iteration steps allowed when applying ! the operator. ! \param verbose The verbosity level. ! ! Verbose levels are: !
    !
  • 0 : print nothing
    • !
  • 1 : print initial and final defect and statistics
    • !
  • 2 : print line for each iteration
    • !
*/ template LoopSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 212,239 ---- } /** ! \brief Set up loop solver ! ! \param op The operator we solve. ! \param sp The scalar product to use, e. g. SeqScalarproduct. ! \param prec The preconditioner to apply in each iteration of the loop. ! Has to inherit from Preconditioner. ! \param reduction The relative defect reduction to achieve when applying ! the operator. ! \param maxit The maximum number of iteration steps allowed when applying ! the operator. ! \param verbose The verbosity level. ! ! Verbose levels are: !
    !
  • 0 : print nothing
    • !
  • 1 : print initial and final defect and statistics
    • !
  • 2 : print line for each iteration
    • !
*/ template LoopSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 207,293 **** } ! //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&) ! virtual void apply (X& x, X& b, InverseOperatorResult& res) ! { ! // clear solver statistics ! res.clear(); ! ! // start a timer ! Timer watch; ! ! // prepare preconditioner ! _prec.pre(x,b); ! ! // overwrite b with defect ! _op.applyscaleadd(-1,x,b); ! ! // compute norm, \todo parallelization ! double def0 = _sp.norm(b); ! ! // printing ! if (_verbose>0) ! { ! printf("=== LoopSolver\n"); ! if (_verbose>1) printf(" Iter Defect Rate\n"); ! if (_verbose>1) printf("%5d %12.4E\n",0,def0); ! } ! ! // allocate correction vector ! X v(x); ! ! // iteration loop ! int i=1; double def=def0; ! for ( ; i<=_maxit; i++ ) ! { ! v = 0; // clear correction ! _prec.apply(v,b); // apply preconditioner ! x += v; // update solution ! _op.applyscaleadd(-1,v,b); // update defect ! double defnew=_sp.norm(b);// comp defect norm ! if (_verbose>1) // print ! printf("%5d %12.4E %12.4g\n",i,defnew,defnew/def); ! def = defnew; // update norm ! if (def0) ! printf("=== rate=%g, T=%g, TIT=%g, IT=%d\n",res.conv_rate,res.elapsed,res.elapsed/i,i); ! } ! ! //! \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) ! virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) ! { ! _reduction = reduction; ! (*this).apply(x,b,res); ! } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; --- 241,336 ---- } ! //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&) ! virtual void apply (X& x, X& b, InverseOperatorResult& res) ! { ! // clear solver statistics ! res.clear(); ! ! // start a timer ! Timer watch; ! ! // prepare preconditioner ! _prec.pre(x,b); ! ! // overwrite b with defect ! _op.applyscaleadd(-1,x,b); ! ! // compute norm, \todo parallelization ! double def0 = _sp.norm(b); ! ! // printing ! if (_verbose>0) ! { ! std::cout << "=== LoopSolver" << std::endl; ! if (_verbose>1) ! { ! this->printHeader(std::cout); ! this->printOutput(std::cout,0,def0); ! } ! } ! ! // allocate correction vector ! X v(x); ! ! // iteration loop ! int i=1; double def=def0; ! for ( ; i<=_maxit; i++ ) ! { ! v = 0; // clear correction ! _prec.apply(v,b); // apply preconditioner ! x += v; // update solution ! _op.applyscaleadd(-1,v,b); // update defect ! double defnew=_sp.norm(b);// comp defect norm ! if (_verbose>1) // print ! this->printOutput(std::cout,i,defnew,def); ! //std::cout << i << " " << defnew << " " << defnew/def << std::endl; ! def = defnew; // update norm ! if (defprintOutput(std::cout,i,def); ! ! // postprocess preconditioner ! _prec.post(x); ! ! // fill statistics ! res.iterations = i; ! res.reduction = def/def0; ! res.conv_rate = pow(res.reduction,1.0/i); ! res.elapsed = watch.elapsed(); ! ! // final print ! if (_verbose>0) ! { ! std::cout << "=== rate=" << res.conv_rate ! << ", T=" << res.elapsed ! << ", TIT=" << res.elapsed/i ! << ", IT=" << i << std::endl; ! } ! } ! ! //! \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) ! virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) ! { ! _reduction = reduction; ! (*this).apply(x,b,res); ! } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; *************** namespace Dune { *** 296,306 **** template class GradientSolver : public InverseOperator { public: ! //! \brief The domain type of the operator we are the inverse for. typedef X domain_type; ! //! \brief The range type of the operator we are the inverse for. typedef X range_type; ! //! \brief The field type of the operator we are the inverse for. typedef typename X::field_type field_type; /*! --- 339,349 ---- template class GradientSolver : public InverseOperator { public: ! //! \brief The domain type of the operator that we do the inverse for. typedef X domain_type; ! //! \brief The range type of the operator that we do the inverse for. typedef X range_type; ! //! \brief The field type of the operator that we do the inverse for. typedef typename X::field_type field_type; /*! *************** namespace Dune { *** 310,316 **** */ template GradientSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 353,359 ---- */ template GradientSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 323,329 **** */ template GradientSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 366,372 ---- */ template GradientSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 344,390 **** X p(x); // create local vectors X q(b); ! double def0 = _sp.norm(b);// compute norm if (_verbose>0) // printing ! { ! printf("=== GradientSolver\n"); ! if (_verbose>1) printf(" Iter Defect Rate\n"); ! if (_verbose>1) printf("%5d %12.4E\n",0,def0); ! } int i=1; double def=def0; // loop variables field_type lambda; ! for ( ; i<=_maxit; i++ ) ! { ! p = 0; // clear correction ! _prec.apply(p,b); // apply preconditioner ! _op.apply(p,q); // q=Ap ! lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization ! x.axpy(lambda,p); // update solution ! b.axpy(-lambda,q); // update defect ! ! double defnew=_sp.norm(b);// comp defect norm ! if (_verbose>1) // print ! printf("%5d %12.4E %12.4g\n",i,defnew,defnew/def); ! def = defnew; // update norm ! if (def0) // final print ! printf("=== rate=%g, T=%g, TIT=%g, IT=%d\n",res.conv_rate,res.elapsed,res.elapsed/i,i); } /*! --- 387,441 ---- X p(x); // create local vectors X q(b); ! double def0 = _sp.norm(b);// compute norm if (_verbose>0) // printing ! { ! std::cout << "=== GradientSolver" << std::endl; ! if (_verbose>1) ! { ! this->printHeader(std::cout); ! this->printOutput(std::cout,0,def0); ! } ! } int i=1; double def=def0; // loop variables field_type lambda; ! for ( ; i<=_maxit; i++ ) ! { ! p = 0; // clear correction ! _prec.apply(p,b); // apply preconditioner ! _op.apply(p,q); // q=Ap ! lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization ! x.axpy(lambda,p); // update solution ! b.axpy(-lambda,q); // update defect ! ! double defnew=_sp.norm(b);// comp defect norm ! if (_verbose>1) // print ! this->printOutput(std::cout,i,defnew,def); ! ! def = defnew; // update norm ! if (defprintOutput(std::cout,i,def); ! _prec.post(x); // postprocess preconditioner res.iterations = i; // fill statistics res.reduction = def/def0; res.conv_rate = pow(res.reduction,1.0/i); res.elapsed = watch.elapsed(); if (_verbose>0) // final print ! std::cout << "=== rate=" << res.conv_rate ! << ", T=" << res.elapsed ! << ", TIT=" << res.elapsed/i ! << ", IT=" << i << std::endl; } /*! *************** namespace Dune { *** 399,411 **** } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; --- 450,462 ---- } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; *************** namespace Dune { *** 414,424 **** template class CGSolver : public InverseOperator { public: ! //! \brief The domain type of the operator we are the inverse for. typedef X domain_type; ! //! \brief The range type of the operator we are the inverse for. typedef X range_type; ! //! \brief The field type of the operator we are the inverse for. typedef typename X::field_type field_type; /*! --- 465,475 ---- template class CGSolver : public InverseOperator { public: ! //! \brief The domain type of the operator to be inverted. typedef X domain_type; ! //! \brief The range type of the operator to be inverted. typedef X range_type; ! //! \brief The field type of the operator to be inverted. typedef typename X::field_type field_type; /*! *************** namespace Dune { *** 453,557 **** */ virtual void apply (X& x, X& b, InverseOperatorResult& res) { ! res.clear(); // clear solver statistics ! Timer watch; // start a timer ! _prec.pre(x,b); // prepare preconditioner ! _op.applyscaleadd(-1,x,b); // overwrite b with defect ! ! X p(x); // the search direction ! X q(x); // a temporary vector ! ! double def0 = _sp.norm(b);// compute norm ! if (def0<1E-30) // convergence check ! { ! res.converged = true; ! res.iterations = 0; // fill statistics ! res.reduction = 0; ! res.conv_rate = 0; ! res.elapsed=0; ! if (_verbose>0) // final print ! printf("=== rate=%g, T=%g, TIT=%g, IT=%d\n",res.conv_rate,res.elapsed,res.elapsed,0); ! return; ! } ! ! if (_verbose>0) // printing ! { ! printf("=== CGSolver\n"); ! if (_verbose>1) printf(" Iter Defect Rate\n"); ! if (_verbose>1) printf("%5d %12.4E\n",0,def0); ! } ! ! // some local variables ! double def=def0; // loop variables ! field_type rho,rholast,lambda,alpha,beta; ! ! // determine initial search direction ! p = 0; // clear correction ! _prec.apply(p,b); // apply preconditioner ! rholast = _sp.dot(p,b); // orthogonalization ! ! // the loop ! int i=1; ! for ( ; i<=_maxit; i++ ) ! { ! // minimize in given search direction p ! _op.apply(p,q); // q=Ap ! alpha = _sp.dot(p,q); // scalar product ! lambda = rholast/alpha; // minimization ! x.axpy(lambda,p); // update solution ! b.axpy(-lambda,q); // update defect ! ! // convergence test ! double defnew=_sp.norm(b);// comp defect norm ! if (_verbose>1) // print ! printf("%5d %12.4E %12.4E\n",i,defnew,defnew/def); ! def = defnew; // update norm ! if (def0) // final print ! printf("=== rate=%g, T=%g, TIT=%g, IT=%d\n",res.conv_rate,res.elapsed,res.elapsed/i,i); ! } /*! \brief Apply inverse operator with given reduction factor. \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) */ ! virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) ! { ! _reduction = reduction; ! (*this).apply(x,b,res); ! } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; --- 504,621 ---- */ virtual void apply (X& x, X& b, InverseOperatorResult& res) { ! res.clear(); // clear solver statistics ! Timer watch; // start a timer ! _prec.pre(x,b); // prepare preconditioner ! _op.applyscaleadd(-1,x,b); // overwrite b with defect ! ! X p(x); // the search direction ! X q(x); // a temporary vector ! ! double def0 = _sp.norm(b);// compute norm ! if (def0<1E-30) // convergence check ! { ! res.converged = true; ! res.iterations = 0; // fill statistics ! res.reduction = 0; ! res.conv_rate = 0; ! res.elapsed=0; ! if (_verbose>0) // final print ! std::cout << "=== rate=" << res.conv_rate ! << ", T=" << res.elapsed << ", TIT=" << res.elapsed ! << ", IT=0" << std::endl; ! return; ! } ! ! if (_verbose>0) // printing ! { ! std::cout << "=== CGSolver" << std::endl; ! if (_verbose>1) { ! this->printHeader(std::cout); ! this->printOutput(std::cout,0,def0); ! } ! } ! ! // some local variables ! double def=def0; // loop variables ! field_type rho,rholast,lambda,alpha,beta; ! ! // determine initial search direction ! p = 0; // clear correction ! _prec.apply(p,b); // apply preconditioner ! rholast = _sp.dot(p,b); // orthogonalization ! ! // the loop ! int i=1; ! for ( ; i<=_maxit; i++ ) ! { ! // minimize in given search direction p ! _op.apply(p,q); // q=Ap ! alpha = _sp.dot(p,q); // scalar product ! lambda = rholast/alpha; // minimization ! x.axpy(lambda,p); // update solution ! b.axpy(-lambda,q); // update defect ! ! // convergence test ! double defnew=_sp.norm(b);// comp defect norm ! ! if (_verbose>1) // print ! this->printOutput(std::cout,i,defnew,def); ! ! def = defnew; // update norm ! if (defprintOutput(std::cout,i,def); ! ! _prec.post(x); // postprocess preconditioner ! res.iterations = i; // fill statistics ! res.reduction = def/def0; ! res.conv_rate = pow(res.reduction,1.0/i); ! res.elapsed = watch.elapsed(); ! ! if (_verbose>0) // final print ! { ! std::cout << "=== rate=" << res.conv_rate ! << ", T=" << res.elapsed ! << ", TIT=" << res.elapsed/i ! << ", IT=" << i << std::endl; ! } ! } /*! \brief Apply inverse operator with given reduction factor. \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&) */ ! virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) ! { ! _reduction = reduction; ! (*this).apply(x,b,res); ! } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; *************** namespace Dune { *** 560,570 **** template class BiCGSTABSolver : public InverseOperator { public: ! //! \brief The domain type of the operator we are the inverse for. typedef X domain_type; ! //! \brief The range type of the operator we are the inverse for. typedef X range_type; ! //! \brief The field type of the operator we are the inverse for. typedef typename X::field_type field_type; /*! --- 624,634 ---- template class BiCGSTABSolver : public InverseOperator { public: ! //! \brief The domain type of the operator to be inverted. typedef X domain_type; ! //! \brief The range type of the operator to be inverted. typedef X range_type; ! //! \brief The field type of the operator to be inverted typedef typename X::field_type field_type; /*! *************** namespace Dune { *** 574,580 **** */ template BiCGSTABSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 638,644 ---- */ template BiCGSTABSolver (L& op, P& prec, ! double reduction, int maxit, int verbose) : ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 587,593 **** */ template BiCGSTABSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); --- 651,657 ---- */ template BiCGSTABSolver (L& op, S& sp, P& prec, ! double reduction, int maxit, int verbose) : _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) { IsTrue< static_cast(L::category) == static_cast(P::category) >::yes(); *************** namespace Dune { *** 601,678 **** */ virtual void apply (X& x, X& b, InverseOperatorResult& res) { ! const double EPSILON=1e-80; ! int it; ! field_type rho, rho_new, alpha, beta, h, omega; ! field_type norm, norm_old, norm_0; ! // ! // get vectors and matrix ! // ! X& r=b; ! X p(x); ! X v(x); ! X t(x); ! X y(x); ! X rt(x); ! ! // ! // begin iteration ! // ! // r = r - Ax; rt = r ! res.clear(); // clear solver statistics ! Timer watch; // start a timer ! _op.applyscaleadd(-1,x,r); // overwrite b with defect ! _prec.pre(x,r); // prepare preconditioner ! rt=r; ! norm = norm_old = norm_0 = _sp.norm(r); ! p=0; ! v=0; ! rho = 1; ! alpha = 1; ! omega = 1; ! ! if (_verbose>0) // printing ! { ! printf("=== BiCGSTABSolver\n"); ! if (_verbose>1) printf(" Iter Defect Rate\n"); ! if (_verbose>1) printf("%5d %12.4E\n",0,norm_0); ! } ! ! if ( norm < (_reduction * norm_0) || norm<1E-30) ! { ! res.converged = 1; ! _prec.post(x); // postprocess preconditioner ! res.iterations = 0; // fill statistics ! res.reduction = 0; ! res.conv_rate = 0; ! res.elapsed = watch.elapsed(); ! return; ! } ! ! // ! // iteration ! // ! it = 0; ! while ( true ) ! { ! // ! // preprocess, set vecsizes etc. ! // ! // rho_new = < rt , r > ! rho_new = _sp.dot(rt,r); ! // look if breakdown occured ! if (std::abs(rho) <= EPSILON) DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho " << rho << " <= EPSILON " << EPSILON << " after " << it << " iterations"); --- 665,747 ---- */ virtual void apply (X& x, X& b, InverseOperatorResult& res) { ! const double EPSILON=1e-80; ! int it; ! field_type rho, rho_new, alpha, beta, h, omega; ! field_type norm, norm_old, norm_0; ! // ! // get vectors and matrix ! // ! X& r=b; ! X p(x); ! X v(x); ! X t(x); ! X y(x); ! X rt(x); ! ! // ! // begin iteration ! // ! // r = r - Ax; rt = r ! res.clear(); // clear solver statistics ! Timer watch; // start a timer ! _op.applyscaleadd(-1,x,r); // overwrite b with defect ! _prec.pre(x,r); // prepare preconditioner ! rt=r; ! norm = norm_old = norm_0 = _sp.norm(r); ! p=0; ! v=0; ! ! rho = 1; ! alpha = 1; ! omega = 1; ! ! if (_verbose>0) // printing ! { ! std::cout << "=== BiCGSTABSolver" << std::endl; ! if (_verbose>1) ! { ! this->printHeader(std::cout); ! this->printOutput(std::cout,0,norm_0); ! //std::cout << " Iter Defect Rate" << std::endl; ! //std::cout << " 0" << std::setw(14) << norm_0 << std::endl; ! } ! } ! ! if ( norm < (_reduction * norm_0) || norm<1E-30) ! { ! res.converged = 1; ! _prec.post(x); // postprocess preconditioner ! res.iterations = 0; // fill statistics ! res.reduction = 0; ! res.conv_rate = 0; ! res.elapsed = watch.elapsed(); ! return; ! } ! // ! // iteration ! // ! it = 0; ! while ( true ) ! { ! // ! // preprocess, set vecsizes etc. ! // ! // rho_new = < rt , r > ! rho_new = _sp.dot(rt,r); ! // look if breakdown occured ! if (std::abs(rho) <= EPSILON) DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho " << rho << " <= EPSILON " << EPSILON << " after " << it << " iterations"); *************** namespace Dune { *** 682,791 **** << " after " << it << " iterations"); ! if (it==0) ! p = r; ! else ! { ! beta = ( rho_new / rho ) * ( alpha / omega ); ! p.axpy(-omega,v); // p = r + beta (p - omega*v) ! p *= beta; ! p += r; ! } ! ! // y = W^-1 * p ! y = 0; ! _prec.apply(y,p); // apply preconditioner ! ! // v = A * y ! _op.apply(y,v); ! ! // alpha = rho_new / < rt, v > ! h = _sp.dot(rt,v); ! ! if ( std::abs(h) < EPSILON ) ! DUNE_THROW(ISTLError,"h=0 in BiCGSTAB"); ! ! alpha = rho_new / h; ! ! // apply first correction to x ! // x <- x + alpha y ! x.axpy(alpha,y); ! ! // r = r - alpha*v ! r.axpy(-alpha,v); ! ! // ! // test stop criteria ! // ! ! it++; ! norm = _sp.norm(r); ! ! if (_verbose>1) // print ! printf("%5d %12.4E %12.4g\n",it,norm,norm/norm_old); ! ! if ( norm < (_reduction * norm_0) ) ! { ! res.converged = 1; ! break; ! } ! if (it >= _maxit) break; ! norm_old = norm; ! // y = W^-1 * r ! y = 0; ! _prec.apply(y,r); ! // t = A * y ! _op.apply(y,t); ! ! // omega = < t, r > / < t, t > ! omega = _sp.dot(t,r)/_sp.dot(t,t); ! // apply second correction to x ! // x <- x + omega y ! x.axpy(omega,y); ! // r = s - omega*t (remember : r = s) ! r.axpy(-omega,t); ! rho = rho_new; ! // ! // test stop criteria ! // ! it++; ! norm = _sp.norm(r); ! if (_verbose>1) // print ! printf("%5d %12.4E %12.4g\n",it,norm,norm/norm_old); ! if ( norm < (_reduction * norm_0) || norm<1E-30) ! { ! res.converged = 1; ! break; ! } ! if (it >= _maxit) break; ! norm_old = norm; ! }// while ! if (_verbose==1) // printing for non verbose ! printf("%5d %12.4E\n",it,norm); ! _prec.post(x); // postprocess preconditioner ! res.iterations = it; // fill statistics ! res.reduction = norm/norm_0; ! res.conv_rate = pow(res.reduction,1.0/it); ! res.elapsed = watch.elapsed(); ! if (_verbose>0) // final print ! printf("=== rate=%g, T=%g, TIT=%g, IT=%d\n",res.conv_rate,res.elapsed,res.elapsed/it,it); } /*! --- 751,868 ---- << " after " << it << " iterations"); ! if (it==0) ! p = r; ! else ! { ! beta = ( rho_new / rho ) * ( alpha / omega ); ! p.axpy(-omega,v); // p = r + beta (p - omega*v) ! p *= beta; ! p += r; ! } ! ! // y = W^-1 * p ! y = 0; ! _prec.apply(y,p); // apply preconditioner ! // v = A * y ! _op.apply(y,v); ! ! // alpha = rho_new / < rt, v > ! h = _sp.dot(rt,v); ! ! if ( std::abs(h) < EPSILON ) ! DUNE_THROW(ISTLError,"h=0 in BiCGSTAB"); ! ! alpha = rho_new / h; ! ! // apply first correction to x ! // x <- x + alpha y ! x.axpy(alpha,y); ! ! // r = r - alpha*v ! r.axpy(-alpha,v); ! ! // ! // test stop criteria ! // ! ! it++; ! norm = _sp.norm(r); ! ! if (_verbose>1) // print ! { ! this->printOutput(std::cout,it,norm,norm_old); ! } ! ! if ( norm < (_reduction * norm_0) ) ! { ! res.converged = 1; ! break; ! } ! ! if (it >= _maxit) break; ! norm_old = norm; ! // y = W^-1 * r ! y = 0; ! _prec.apply(y,r); ! // t = A * y ! _op.apply(y,t); ! ! // omega = < t, r > / < t, t > ! omega = _sp.dot(t,r)/_sp.dot(t,t); ! // apply second correction to x ! // x <- x + omega y ! x.axpy(omega,y); ! // r = s - omega*t (remember : r = s) ! r.axpy(-omega,t); ! rho = rho_new; ! // ! // test stop criteria ! // ! it++; ! norm = _sp.norm(r); ! if (_verbose > 1) // print ! { ! this->printOutput(std::cout,it,norm,norm_old); ! } ! if ( norm < (_reduction * norm_0) || norm<1E-30) ! { ! res.converged = 1; ! break; ! } ! if (it >= _maxit) break; ! norm_old = norm; ! }// while ! ! if (_verbose==1) // printing for non verbose ! this->printOutput(std::cout,it,norm); ! _prec.post(x); // postprocess preconditioner ! res.iterations = it; // fill statistics ! res.reduction = norm/norm_0; ! res.conv_rate = pow(res.reduction,1.0/it); ! res.elapsed = watch.elapsed(); ! if (_verbose>0) // final print ! std::cout << "=== rate=" << res.conv_rate ! << ", T=" << res.elapsed ! << ", TIT=" << res.elapsed/it ! << ", IT=" << it << std::endl; } /*! *************** namespace Dune { *** 800,812 **** } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; --- 877,889 ---- } private: ! SeqScalarProduct ssp; ! LinearOperator& _op; ! Preconditioner& _prec; ! ScalarProduct& _sp; ! double _reduction; ! int _maxit; ! int _verbose; }; Only in ../dune-istl-1.1/istl/: superlu.hh Only in ../dune-istl-1.1/istl/: supermatrix.hh Common subdirectories: istl/test and ../dune-istl-1.1/istl/test Common subdirectories: istl/tutorial and ../dune-istl-1.1/istl/tutorial Only in ../dune-istl-1.1/istl/: vbcrsmatrix.hh diff -p --exclude=Makefile.am --exclude=Makefile.in istl/vbvector.hh ../dune-istl-1.1/istl/vbvector.hh *** istl/vbvector.hh Wed Aug 15 16:06:05 2007 --- ../dune-istl-1.1/istl/vbvector.hh Thu Apr 3 03:59:22 2008 *************** *** 1,7 **** #ifndef DUNE_VBVECTOR_HH #define DUNE_VBVECTOR_HH ! #include #include #include --- 1,7 ---- #ifndef DUNE_VBVECTOR_HH #define DUNE_VBVECTOR_HH ! #include #include #include