This file is indexed.

/usr/include/dune/istl/novlpschwarz.hh is in libdune-istl-dev 2.4.1-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_NOVLPSCHWARZ_HH
#define DUNE_ISTL_NOVLPSCHWARZ_HH

#include <iostream>              // for input/output to shell
#include <fstream>               // for input/output to files
#include <vector>                // STL vector class
#include <sstream>

#include <cmath>                // Yes, we do some math here

#include <dune/common/timer.hh>

#include "io.hh"
#include "bvector.hh"
#include "vbvector.hh"
#include "bcrsmatrix.hh"
#include "io.hh"
#include "gsetc.hh"
#include "ilu.hh"
#include "operators.hh"
#include "solvers.hh"
#include "preconditioners.hh"
#include "scalarproducts.hh"
#include "owneroverlapcopy.hh"

namespace Dune {

  /**
   * @defgroup ISTL_Parallel Parallel Solvers
   * @ingroup ISTL_Solvers
   * Instead of using parallel data structures (matrices and vectors) that
   * (implicitly) know the data distribution and communication patterns,
   * there is a clear separation of the parallel data composition together
   *  with the communication APIs from the data structures. This allows for
   * implementing overlapping and nonoverlapping domain decompositions as
   * well as data parallel parallelisation approaches.
   *
   * The \ref ISTL_Solvers "solvers" can easily be turned into parallel solvers
   * initializing them with matching parallel subclasses of the the base classes
   * ScalarProduct, Preconditioner and LinearOperator.
   *
   * The information of the data distribution is provided by OwnerOverlapCopyCommunication
   * of \ref ISTL_Comm "communication API".
   *
   * Currently only data parallel versions are shipped with dune-istl. Domain
   * decomposition can be found in module dune-dd.
   */
  /**
     @addtogroup ISTL_Operators
     @{
   */

  /**
   * @brief A nonoverlapping operator with communication object.
   */
  template<class M, class X, class Y, class C>
  class NonoverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
  {
  public:
    //! \brief The type of the matrix we operate on.
    typedef M matrix_type;
    //! \brief The type of the domain.
    typedef X domain_type;
    //! \brief The type of the range.
    typedef Y range_type;
    //! \brief The field type of the range
    typedef typename X::field_type field_type;
    //! \brief The type of the communication object
    typedef C communication_type;

    typedef typename C::PIS PIS;
    typedef typename C::RI RI;
    typedef typename RI::RemoteIndexList RIL;
    typedef typename RI::const_iterator RIIterator;
    typedef typename RIL::const_iterator RILIterator;
    typedef typename M::ConstColIterator ColIterator;
    typedef typename M::ConstRowIterator RowIterator;
    typedef std::multimap<int,int> MM;
    typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
    typedef typename RIMap::iterator RIMapit;

    enum {
      //! \brief The solver category.
      category=SolverCategory::nonoverlapping
    };

    /**
     * @brief constructor: just store a reference to a matrix.
     *
     * @param A The assembled matrix.
     * @param com The communication object for syncing owner and copy
     * data points. (E.~g. OwnerOverlapCommunication )
     */
    NonoverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
      : _A_(A), communication(com), buildcomm(true)
    {}

    //! apply operator to x:  \f$ y = A(x) \f$
    virtual void apply (const X& x, Y& y) const
    {
      y = 0;
      novlp_op_apply(x,y,1);
      communication.addOwnerCopyToOwnerCopy(y,y);
    }

    //! apply operator to x, scale and add:  \f$ y = y + \alpha A(x) \f$
    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
    {
      // only apply communication to alpha*A*x to make it consistent,
      // y already has to be consistent.
      Y y1(y);
      y = 0;
      novlp_op_apply(x,y,alpha);
      communication.addOwnerCopyToOwnerCopy(y,y);
      y += y1;
    }

    //! get matrix via *
    virtual const matrix_type& getmat () const
    {
      return _A_;
    }

    void novlp_op_apply (const X& x, Y& y, field_type alpha) const
    {
      //get index sets
      const PIS& pis=communication.indexSet();
      const RI& ri = communication.remoteIndices();

      // at the beginning make a multimap "bordercontribution".
      // process has i and j as border dofs but is not the owner
      // => only contribute to Ax if i,j is in bordercontribution
      if (buildcomm == true) {

        // set up mask vector
        if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())) {
          mask.resize(x.size());
          for (typename std::vector<double>::size_type i=0; i<mask.size(); i++)
            mask[i] = 1;
          for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
            if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
              mask[i->local().local()] = 0;
            else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
              mask[i->local().local()] = 2;
        }

        for (MM::iterator iter = bordercontribution.begin();
             iter != bordercontribution.end(); ++iter)
          bordercontribution.erase(iter);
        std::map<int,int> owner; //key: local index i, value: process, that owns i
        RIMap rimap;

        // for each local index make multimap rimap:
        // key: local index i, data: pair of process that knows i and pointer to RI entry
        for (RowIterator i = _A_.begin(); i != _A_.end(); ++i)
          if (mask[i.index()] == 0)
            for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote) {
              RIL& ril = *(remote->second.first);
              for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
                if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
                  if (rindex->localIndexPair().local().local() == i.index()) {
                    rimap.insert
                      (std::make_pair(i.index(),
                                      std::pair<int,RILIterator>(remote->first, rindex)));
                    if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
                      owner.insert(std::make_pair(i.index(),remote->first));
                  }
            }

        int iowner = 0;
        for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
          if (mask[i.index()] == 0) {
            std::map<int,int>::iterator it = owner.find(i.index());
            iowner = it->second;
            std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
            for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
              if (mask[j.index()] == 0) {
                bool flag = true;
                for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi) {
                  std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
                  for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
                    if (foundj->second.first == foundi->second.first)
                      if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
                          || foundj->second.first == iowner
                          || foundj->second.first  < communication.communicator().rank()) {
                        flag = false;
                        continue;
                      }
                  if (flag == false)
                    continue;
                }
                // donĀ“t contribute to Ax if
                // 1. the owner of j has i as interior/border dof
                // 2. iowner has j as interior/border dof
                // 3. there is another process with smaller rank that has i and j
                // as interor/border dofs
                // if the owner of j does not have i as interior/border dof,
                // it will not be taken into account
                if (flag==true)
                  bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
              }
            }
          }
        }
        buildcomm = false;
      }

      //compute alpha*A*x nonoverlapping case
      for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) {
        if (mask[i.index()] == 0) {
          //dof doesn't belong to process but is border (not ghost)
          for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) {
            if (mask[j.index()] == 1) //j is owner => then sum entries
              (*j).usmv(alpha,x[j.index()],y[i.index()]);
            else if (mask[j.index()] == 0) {
              std::pair<MM::iterator, MM::iterator> itp =
                bordercontribution.equal_range(i.index());
              for (MM::iterator it = itp.first; it != itp.second; ++it)
                if ((*it).second == (int)j.index())
                  (*j).usmv(alpha,x[j.index()],y[i.index()]);
            }
          }
        }
        else if (mask[i.index()] == 1) {
          for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j)
            if (mask[j.index()] != 2)
              (*j).usmv(alpha,x[j.index()],y[i.index()]);
        }
      }
    }

  private:
    const matrix_type& _A_;
    const communication_type& communication;
    mutable bool buildcomm;
    mutable std::vector<double> mask;
    mutable std::multimap<int,int>  bordercontribution;
  };

  /** @} */

  /**
   * @addtogroup ISTL_SP
   * @{
   */
  /**
   * \brief Nonoverlapping Scalar Product with communication object.
   *
   * Consistent vectors in interior and border are assumed.
   */
  template<class X, class C>
  class NonoverlappingSchwarzScalarProduct : public ScalarProduct<X>
  {
  public:
    //! \brief The type of the domain.
    typedef X domain_type;
    //!  \brief The type of the range
    typedef typename X::field_type field_type;
    //! \brief The type of the communication object
    typedef C communication_type;

    //! define the category
    enum {category=SolverCategory::nonoverlapping};

    /*! \brief Constructor
     * \param com The communication object for syncing owner and copy
     * data points. (E.~g. OwnerOverlapCommunication )
     */
    NonoverlappingSchwarzScalarProduct (const communication_type& com)
      : communication(com)
    {}

    /*! \brief Dot product of two vectors.
       It is assumed that the vectors are consistent on the interior+border
       partition.
     */
    virtual field_type dot (const X& x, const X& y)
    {
      field_type result;
      communication.dot(x,y,result);
      return result;
    }

    /*! \brief Norm of a right-hand side vector.
       The vector must be consistent on the interior+border partition
     */
    virtual double norm (const X& x)
    {
      return communication.norm(x);
    }

    /*! \brief make additive vector consistent
     */
    void make_consistent (X& x) const
    {
      communication.copyOwnerToAll(x,x);
    }

  private:
    const communication_type& communication;
  };

  template<class X, class C>
  struct ScalarProductChooser<X,C,SolverCategory::nonoverlapping>
  {
    /** @brief The type of the scalar product for the nonoverlapping case. */
    typedef NonoverlappingSchwarzScalarProduct<X,C> ScalarProduct;
    /** @brief The type of the communication object to use. */
    typedef C communication_type;

    enum {
      /** @brief The solver category. */
      solverCategory=SolverCategory::nonoverlapping
    };

    static ScalarProduct* construct(const communication_type& comm)
    {
      return new ScalarProduct(comm);
    }
  };

  namespace Amg
  {
    template<class T> class ConstructionTraits;
  }

  /**
   * @brief Nonoverlapping parallel preconditioner.
   *
   * This is essentially a wrapper that take a sequential
   * preconditoner. In each step the sequential preconditioner
   * is applied and then all owner data points are updated on
   * all other processes.
   */

  template<class C, class P>
  class NonoverlappingBlockPreconditioner
    : public Dune::Preconditioner<typename P::domain_type,typename P::range_type> {
    friend class Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
  public:
    //! \brief The domain type of the preconditioner.
    typedef typename P::domain_type domain_type;
    //! \brief The range type of the preconditioner.
    typedef typename P::range_type range_type;
    //! \brief The type of the communication object.
    typedef C communication_type;

    // define the category
    enum {
      //! \brief The category the preconditioner is part of.
      category=SolverCategory::nonoverlapping
    };

    /*! \brief Constructor.

       constructor gets all parameters to operate the prec.
       \param prec The sequential preconditioner.
       \param c The communication object for syncing owner and copy
       data points. (E.~g. OwnerOverlapCommunication )
     */
    NonoverlappingBlockPreconditioner (P& prec, const communication_type& c)
      : preconditioner(prec), communication(c)
    {}

    /*!
       \brief Prepare the preconditioner.

       \copydoc Preconditioner::pre(domain_type&,range_type&)
     */
    virtual void pre (domain_type& x, range_type& b)
    {
      preconditioner.pre(x,b);
    }

    /*!
       \brief Apply the preconditioner

       \copydoc Preconditioner::apply(domain_type&,const range_type&)
     */
    virtual void apply (domain_type& v, const range_type& d)
    {
      // block preconditioner equivalent to WrappedPreconditioner from
      // pdelab/backend/ovlpistsolverbackend.hh,
      // but not to BlockPreconditioner from schwarz.hh
      preconditioner.apply(v,d);
      communication.addOwnerCopyToOwnerCopy(v,v);
    }

    /*!
       \brief Clean up.

       \copydoc Preconditioner::post(domain_type&)
     */
    virtual void post (domain_type& x)
    {
      preconditioner.post(x);
    }

  private:
    //! \brief a sequential preconditioner
    P& preconditioner;

    //! \brief the communication object
    const communication_type& communication;
  };

  /** @} end documentation */

} // end namespace

#endif