~ubuntu-branches/ubuntu/saucy/deal.ii/saucy

« back to all changes in this revision

Viewing changes to lac/include/lac/petsc_block_vector.h

  • Committer: Package Import Robot
  • Author(s): "Adam C. Powell, IV", Adam C. Powell, IV, Christophe Trophime
  • Date: 2012-02-21 06:57:30 UTC
  • mfrom: (3.1.7 sid)
  • Revision ID: package-import@ubuntu.com-20120221065730-r2iz70lg557wcd2e
Tags: 7.1.0-1
[ Adam C. Powell, IV ]
* New upstream (closes: #652057).
* Updated to use PETSc and SLEPc 3.2, and forward-ported all patches.
* Removed NetCDF Build-Depends because it uses serial HDF5.
* Made Sacado cmath patch work with new configure.
* Moved -dev package symlink lines in rules to arch all section.

[ Christophe Trophime ]
* debian/rules:
   - add dh_strip --dbg-package to generate dbg package (closes: #652058)
   - add .install files to simplify rules
* Add support for mumps, arpack (closes: #637655)
* Add patch for slepc 3.2 (closes: #659245)

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
//---------------------------------------------------------------------------
2
 
//    $Id: petsc_block_vector.h 20934 2010-04-02 13:09:20Z bangerth $
3
 
//    Version: $Name$
4
 
//
5
 
//    Copyright (C) 2004, 2005, 2006, 2007, 2010 by the deal.II authors
6
 
//
7
 
//    This file is subject to QPL and may not be  distributed
8
 
//    without copyright and license information. Please refer
9
 
//    to the file deal.II/doc/license.html for the  text  and
10
 
//    further information on this license.
11
 
//
12
 
//---------------------------------------------------------------------------
13
 
#ifndef __deal2__petsc_block_vector_h
14
 
#define __deal2__petsc_block_vector_h
15
 
 
16
 
 
17
 
#include <base/config.h>
18
 
 
19
 
#ifdef DEAL_II_USE_PETSC
20
 
 
21
 
#  include <lac/petsc_vector.h>
22
 
#  include <lac/petsc_parallel_block_vector.h>
23
 
#  include <lac/block_indices.h>
24
 
#  include <lac/block_vector_base.h>
25
 
#  include <lac/exceptions.h>
26
 
 
27
 
DEAL_II_NAMESPACE_OPEN
28
 
 
29
 
 
30
 
 
31
 
namespace PETScWrappers
32
 
{
33
 
/*! @addtogroup PETScWrappers
34
 
 *@{
35
 
 */
36
 
 
37
 
/**
38
 
 * An implementation of block vectors based on the vector class implemented in
39
 
 * PETScWrappers. While the base class provides for most of the interface,
40
 
 * this class handles the actual allocation of vectors and provides functions
41
 
 * that are specific to the underlying vector type.
42
 
 *
43
 
 * @ingroup Vectors
44
 
 * @see @ref GlossBlockLA "Block (linear algebra)"
45
 
 * @author Wolfgang Bangerth, 2004
46
 
 */
47
 
  class BlockVector : public BlockVectorBase<Vector>
48
 
  {
49
 
    public:
50
 
                                       /**
51
 
                                        * Typedef the base class for simpler
52
 
                                        * access to its own typedefs.
53
 
                                        */
54
 
      typedef BlockVectorBase<Vector> BaseClass;
55
 
    
56
 
                                       /**
57
 
                                        * Typedef the type of the underlying
58
 
                                        * vector.
59
 
                                        */
60
 
      typedef BaseClass::BlockType  BlockType;
61
 
 
62
 
                                       /**
63
 
                                        * Import the typedefs from the base
64
 
                                        * class.
65
 
                                        */
66
 
      typedef BaseClass::value_type      value_type;
67
 
      typedef BaseClass::pointer         pointer;
68
 
      typedef BaseClass::const_pointer   const_pointer;
69
 
      typedef BaseClass::reference       reference;
70
 
      typedef BaseClass::const_reference const_reference;
71
 
      typedef BaseClass::size_type       size_type;
72
 
      typedef BaseClass::iterator        iterator;
73
 
      typedef BaseClass::const_iterator  const_iterator;
74
 
 
75
 
                                       /**
76
 
                                        *  Constructor. There are three
77
 
                                        *  ways to use this
78
 
                                        *  constructor. First, without
79
 
                                        *  any arguments, it generates
80
 
                                        *  an object with no
81
 
                                        *  blocks. Given one argument,
82
 
                                        *  it initializes <tt>num_blocks</tt>
83
 
                                        *  blocks, but these blocks have
84
 
                                        *  size zero. The third variant
85
 
                                        *  finally initializes all
86
 
                                        *  blocks to the same size
87
 
                                        *  <tt>block_size</tt>.
88
 
                                        *
89
 
                                        *  Confer the other constructor
90
 
                                        *  further down if you intend to
91
 
                                        *  use blocks of different
92
 
                                        *  sizes.
93
 
                                        */
94
 
      explicit BlockVector (const unsigned int num_blocks = 0,
95
 
                            const unsigned int block_size = 0);
96
 
    
97
 
                                       /**
98
 
                                        * Copy-Constructor. Dimension set to
99
 
                                        * that of V, all components are copied
100
 
                                        * from V
101
 
                                        */
102
 
      BlockVector (const BlockVector  &V);
103
 
 
104
 
                                       /**
105
 
                                        * Copy-constructor: copy the values
106
 
                                        * from a PETSc wrapper parallel block
107
 
                                        * vector class.
108
 
                                        * 
109
 
                                        *
110
 
                                        * Note that due to the communication
111
 
                                        * model of MPI, @em all processes have
112
 
                                        * to actually perform this operation,
113
 
                                        * even if they do not use the
114
 
                                        * result. It is not sufficient if only
115
 
                                        * one processor tries to copy the
116
 
                                        * elements from the other processors
117
 
                                        * over to its own process space.
118
 
                                        */
119
 
      explicit BlockVector (const MPI::BlockVector &v);
120
 
 
121
 
                                       /**
122
 
                                        * Constructor. Set the number of
123
 
                                        * blocks to <tt>n.size()</tt> and
124
 
                                        * initialize each block with
125
 
                                        * <tt>n[i]</tt> zero elements.
126
 
                                        */
127
 
      BlockVector (const std::vector<unsigned int> &n);
128
 
 
129
 
                                       /**
130
 
                                        * Constructor. Set the number of
131
 
                                        * blocks to
132
 
                                        * <tt>n.size()</tt>. Initialize the
133
 
                                        * vector with the elements
134
 
                                        * pointed to by the range of
135
 
                                        * iterators given as second and
136
 
                                        * third argument. Apart from the
137
 
                                        * first argument, this
138
 
                                        * constructor is in complete
139
 
                                        * analogy to the respective
140
 
                                        * constructor of the
141
 
                                        * <tt>std::vector</tt> class, but the
142
 
                                        * first argument is needed in
143
 
                                        * order to know how to subdivide
144
 
                                        * the block vector into
145
 
                                        * different blocks.
146
 
                                        */
147
 
      template <typename InputIterator>
148
 
      BlockVector (const std::vector<unsigned int> &n,
149
 
                   const InputIterator              first,
150
 
                   const InputIterator              end);
151
 
    
152
 
                                       /**
153
 
                                        * Destructor. Clears memory
154
 
                                        */
155
 
      ~BlockVector ();
156
 
 
157
 
                                       /**
158
 
                                        * Copy operator: fill all components of
159
 
                                        * the vector with the given scalar
160
 
                                        * value.
161
 
                                        */
162
 
      BlockVector & operator = (const value_type s);
163
 
 
164
 
                                       /**
165
 
                                        * Copy operator for arguments of the
166
 
                                        * same type.
167
 
                                        */
168
 
      BlockVector &
169
 
      operator= (const BlockVector &V);
170
 
 
171
 
                                       /**
172
 
                                        * Copy all the elements of the
173
 
                                        * parallel block vector @p v into this
174
 
                                        * local vector. Note that due to the
175
 
                                        * communication model of MPI, @em all
176
 
                                        * processes have to actually perform
177
 
                                        * this operation, even if they do not
178
 
                                        * use the result. It is not sufficient
179
 
                                        * if only one processor tries to copy
180
 
                                        * the elements from the other
181
 
                                        * processors over to its own process
182
 
                                        * space.
183
 
                                        */
184
 
      BlockVector &
185
 
      operator = (const MPI::BlockVector &v);
186
 
 
187
 
                                       /**
188
 
                                        * Reinitialize the BlockVector to
189
 
                                        * contain <tt>num_blocks</tt> blocks of
190
 
                                        * size <tt>block_size</tt> each.
191
 
                                        *
192
 
                                        * If <tt>fast==false</tt>, the vector
193
 
                                        * is filled with zeros.
194
 
                                        */
195
 
      void reinit (const unsigned int num_blocks,
196
 
                   const unsigned int block_size,
197
 
                   const bool fast = false);
198
 
  
199
 
                                       /**
200
 
                                        * Reinitialize the BlockVector such
201
 
                                        * that it contains
202
 
                                        * <tt>block_sizes.size()</tt>
203
 
                                        * blocks. Each block is reinitialized
204
 
                                        * to dimension
205
 
                                        * <tt>block_sizes[i]</tt>.
206
 
                                        *
207
 
                                        * If the number of blocks is the
208
 
                                        * same as before this function
209
 
                                        * was called, all vectors remain
210
 
                                        * the same and reinit() is
211
 
                                        * called for each vector.
212
 
                                        *
213
 
                                        * If <tt>fast==false</tt>, the vector
214
 
                                        * is filled with zeros.
215
 
                                        *
216
 
                                        * Note that you must call this
217
 
                                        * (or the other reinit()
218
 
                                        * functions) function, rather
219
 
                                        * than calling the reinit()
220
 
                                        * functions of an individual
221
 
                                        * block, to allow the block
222
 
                                        * vector to update its caches of
223
 
                                        * vector sizes. If you call
224
 
                                        * reinit() on one of the
225
 
                                        * blocks, then subsequent
226
 
                                        * actions on this object may
227
 
                                        * yield unpredictable results
228
 
                                        * since they may be routed to
229
 
                                        * the wrong block.
230
 
                                        */ 
231
 
      void reinit (const std::vector<unsigned int> &N,
232
 
                   const bool                       fast=false);
233
 
    
234
 
                                       /**
235
 
                                        * Change the dimension to that
236
 
                                        * of the vector <tt>V</tt>. The same
237
 
                                        * applies as for the other
238
 
                                        * reinit() function.
239
 
                                        *
240
 
                                        * The elements of <tt>V</tt> are not
241
 
                                        * copied, i.e.  this function is
242
 
                                        * the same as calling <tt>reinit
243
 
                                        * (V.size(), fast)</tt>.
244
 
                                        *
245
 
                                        * Note that you must call this
246
 
                                        * (or the other reinit()
247
 
                                        * functions) function, rather
248
 
                                        * than calling the reinit()
249
 
                                        * functions of an individual
250
 
                                        * block, to allow the block
251
 
                                        * vector to update its caches of
252
 
                                        * vector sizes. If you call
253
 
                                        * reinit() of one of the
254
 
                                        * blocks, then subsequent
255
 
                                        * actions of this object may
256
 
                                        * yield unpredictable results
257
 
                                        * since they may be routed to
258
 
                                        * the wrong block.
259
 
                                        */
260
 
      void reinit (const BlockVector &V,
261
 
                   const bool         fast=false);
262
 
    
263
 
                                       /**
264
 
                                        * Swap the contents of this
265
 
                                        * vector and the other vector
266
 
                                        * <tt>v</tt>. One could do this
267
 
                                        * operation with a temporary
268
 
                                        * variable and copying over the
269
 
                                        * data elements, but this
270
 
                                        * function is significantly more
271
 
                                        * efficient since it only swaps
272
 
                                        * the pointers to the data of
273
 
                                        * the two vectors and therefore
274
 
                                        * does not need to allocate
275
 
                                        * temporary storage and move
276
 
                                        * data around.
277
 
                                        *
278
 
                                        * Limitation: right now this
279
 
                                        * function only works if both
280
 
                                        * vectors have the same number
281
 
                                        * of blocks. If needed, the
282
 
                                        * numbers of blocks should be
283
 
                                        * exchanged, too.
284
 
                                        *
285
 
                                        * This function is analog to the
286
 
                                        * the swap() function of all C++
287
 
                                        * standard containers. Also,
288
 
                                        * there is a global function
289
 
                                        * swap(u,v) that simply calls
290
 
                                        * <tt>u.swap(v)</tt>, again in analogy
291
 
                                        * to standard functions.
292
 
                                        */
293
 
      void swap (BlockVector &v);    
294
 
 
295
 
                                     /**
296
 
                                      * Print to a stream.
297
 
                                      */
298
 
      void print (std::ostream       &out,
299
 
                  const unsigned int  precision = 3,
300
 
                  const bool          scientific = true,
301
 
                  const bool          across = true) const;
302
 
 
303
 
                                       /** @addtogroup Exceptions
304
 
                                        * @{ */
305
 
 
306
 
                                       /**
307
 
                                        * Exception
308
 
                                        */
309
 
      DeclException0 (ExcIteratorRangeDoesNotMatchVectorSize);
310
 
                                       ///@}
311
 
  };
312
 
 
313
 
/*@}*/
314
 
 
315
 
/*----------------------- Inline functions ----------------------------------*/
316
 
 
317
 
 
318
 
 
319
 
  inline
320
 
  BlockVector::BlockVector (const unsigned int n_blocks,
321
 
                            const unsigned int block_size)
322
 
  {
323
 
    reinit (n_blocks, block_size);
324
 
  }
325
 
 
326
 
 
327
 
 
328
 
  inline
329
 
  BlockVector::BlockVector (const std::vector<unsigned int> &n)
330
 
  {
331
 
    reinit (n, false);
332
 
  }
333
 
 
334
 
 
335
 
  inline
336
 
  BlockVector::BlockVector (const BlockVector& v)
337
 
                  :
338
 
                  BlockVectorBase<Vector > ()
339
 
  {
340
 
    this->components.resize (v.n_blocks());
341
 
    block_indices = v.block_indices;
342
 
  
343
 
    for (unsigned int i=0; i<this->n_blocks(); ++i)
344
 
      this->components[i] = v.components[i];
345
 
  }
346
 
 
347
 
 
348
 
 
349
 
  inline
350
 
  BlockVector::BlockVector (const MPI::BlockVector& v)
351
 
                  :
352
 
                  BlockVectorBase<Vector > ()
353
 
  {
354
 
    this->components.resize (v.get_block_indices().size());
355
 
    block_indices = v.get_block_indices();
356
 
  
357
 
    for (unsigned int i=0; i<this->n_blocks(); ++i)
358
 
      this->components[i] = v.block(i);
359
 
  }
360
 
 
361
 
  
362
 
 
363
 
  template <typename InputIterator>
364
 
  BlockVector::BlockVector (const std::vector<unsigned int> &n,
365
 
                            const InputIterator              first,
366
 
                            const InputIterator              end)
367
 
  {
368
 
                                     // first set sizes of blocks, but
369
 
                                     // don't initialize them as we will
370
 
                                     // copy elements soon
371
 
    reinit (n, true);
372
 
    InputIterator start = first;
373
 
    for (unsigned int b=0; b<n.size(); ++b)
374
 
      {
375
 
        InputIterator end = start;
376
 
        std::advance (end, static_cast<signed int>(n[b]));
377
 
 
378
 
        for (unsigned int i=0; i<n[b]; ++i, ++start)
379
 
          this->block(b)(i) = *start;
380
 
      }
381
 
    Assert (start == end, ExcIteratorRangeDoesNotMatchVectorSize());
382
 
  }
383
 
 
384
 
 
385
 
 
386
 
  inline
387
 
  BlockVector &
388
 
  BlockVector::operator = (const value_type s)
389
 
  {
390
 
    BaseClass::operator = (s);
391
 
    return *this;
392
 
  }
393
 
 
394
 
 
395
 
 
396
 
  inline
397
 
  BlockVector &
398
 
  BlockVector::operator = (const BlockVector &v)
399
 
  {
400
 
    BaseClass::operator = (v);
401
 
    return *this;
402
 
  }
403
 
 
404
 
 
405
 
 
406
 
  inline
407
 
  BlockVector &
408
 
  BlockVector::operator = (const MPI::BlockVector &v)
409
 
  {
410
 
    BaseClass::operator = (v);
411
 
    return *this;
412
 
  }
413
 
 
414
 
 
415
 
 
416
 
  inline
417
 
  BlockVector::~BlockVector ()
418
 
  {}
419
 
 
420
 
  
421
 
  inline
422
 
  void
423
 
  BlockVector::reinit (const unsigned int n_bl,
424
 
                       const unsigned int bl_sz,
425
 
                       const bool         fast)
426
 
  {
427
 
    std::vector<unsigned int> n(n_bl, bl_sz);
428
 
    reinit(n, fast);
429
 
  }
430
 
 
431
 
  
432
 
 
433
 
  inline
434
 
  void
435
 
  BlockVector::reinit (const std::vector<unsigned int> &n,
436
 
                       const bool                       fast)
437
 
  {
438
 
    block_indices.reinit (n);
439
 
    if (this->components.size() != this->n_blocks())
440
 
      this->components.resize(this->n_blocks());
441
 
  
442
 
    for (unsigned int i=0; i<this->n_blocks(); ++i)
443
 
      this->components[i].reinit(n[i], fast);
444
 
  }
445
 
 
446
 
 
447
 
  inline
448
 
  void
449
 
  BlockVector::reinit (const BlockVector& v,
450
 
                       const bool fast)
451
 
  {
452
 
    block_indices = v.get_block_indices();
453
 
    if (this->components.size() != this->n_blocks())
454
 
      this->components.resize(this->n_blocks());
455
 
  
456
 
    for (unsigned int i=0;i<this->n_blocks();++i)
457
 
      block(i).reinit(v.block(i), fast);
458
 
  }
459
 
 
460
 
  
461
 
 
462
 
  inline
463
 
  void
464
 
  BlockVector::swap (BlockVector &v)
465
 
  {
466
 
    Assert (this->n_blocks() == v.n_blocks(),
467
 
            ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
468
 
  
469
 
    for (unsigned int i=0; i<this->n_blocks(); ++i)
470
 
      this->components[i].swap (v.components[i]);
471
 
    ::dealii::swap (this->block_indices, v.block_indices);
472
 
  }
473
 
 
474
 
 
475
 
 
476
 
  inline
477
 
  void 
478
 
  BlockVector::print (std::ostream       &out,
479
 
                      const unsigned int  precision,
480
 
                      const bool          scientific,
481
 
                      const bool          across) const
482
 
  {
483
 
    for (unsigned int i=0;i<this->n_blocks();++i)
484
 
      {
485
 
        if (across)
486
 
          out << 'C' << i << ':';
487
 
        else
488
 
          out << "Component " << i << std::endl;
489
 
        this->components[i].print(out, precision, scientific, across);
490
 
      }
491
 
  }
492
 
 
493
 
 
494
 
  
495
 
 
496
 
/**
497
 
 * Global function which overloads the default implementation
498
 
 * of the C++ standard library which uses a temporary object. The
499
 
 * function simply exchanges the data of the two vectors.
500
 
 *
501
 
 * @relates PETScWrappers::BlockVector
502
 
 * @author Wolfgang Bangerth, 2000
503
 
 */
504
 
  inline
505
 
  void swap (BlockVector &u,
506
 
             BlockVector &v)
507
 
  {
508
 
    u.swap (v);
509
 
  }
510
 
  
511
 
}
512
 
 
513
 
 
514
 
DEAL_II_NAMESPACE_CLOSE
515
 
 
516
 
#endif  // DEAL_II_USE_PETSC
517
 
 
518
 
#endif