← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/la/SparsityPattern.cpp

  • Committer: Anders Logg
  • Date: 2008-05-15 16:02:44 UTC
  • Revision ID: logg@simula.no-20080515160244-e2a4wkfxfzmts16s
Add new class SingularSolver:

This class provides a linear solver for singular linear systems
Ax = b where A has a one-dimensional null-space (kernel). This
may happen for example when solving Poisson's equation with
pure Neumann boundary conditions.

The solver attempts to create an extended non-singular system
by adding the constraint [1, 1, 1, ...]^T x = 0.

If an optional mass matrix M is supplied, the solver attempts
to create an extended non-singular system by adding the
constraint m^T x = 0 where m is the lumped mass matrix. This
corresponds to setting the average (integral) of the finite
element function with coefficients x to zero.

The solver makes not attempt to check that the null-space is
indeed one-dimensional. It is also assumed that the system
Ax = b retains its sparsity pattern between calls to solve().

Show diffs side-by-side

added added

removed removed

Lines of Context:
dolfin/la/SparsityPattern.cpp
1
1
// Copyright (C) 2007 Garth N. Wells.
2
2
// Licensed under the GNU LGPL Version 2.1.
3
3
//
4
 
// Modified by Magnus Vikstrom 2008.
 
4
// Modified by Magnus Vikstrom, 2008.
 
5
// Modified by Anders Logg, 2008.
5
6
//
6
7
// First added:  2007-03-13
7
 
// Last changed: 2008-01-24
 
8
// Last changed: 2008-05-15
8
9
 
9
10
#include <dolfin/log/dolfin_log.h>
10
11
#include "SparsityPattern.h"
15
16
using namespace dolfin;
16
17
 
17
18
//-----------------------------------------------------------------------------
18
 
SparsityPattern::SparsityPattern() 
 
19
SparsityPattern::SparsityPattern(uint M, uint N) : range(0)
 
20
{
 
21
  uint dims[2];
 
22
  dims[0] = M;
 
23
  dims[1] = N;
 
24
  init(2, dims);
 
25
}
 
26
//-----------------------------------------------------------------------------
 
27
SparsityPattern::SparsityPattern(uint M) : range(0)
 
28
{
 
29
  uint dims[2];
 
30
  dims[0] = M;
 
31
  dims[1] = 0;
 
32
  init(1, dims);
 
33
}
 
34
//-----------------------------------------------------------------------------
 
35
  SparsityPattern::SparsityPattern() : range(0)
19
36
{
20
37
  dim[0] = 0;
21
38
  dim[1] = 0;
25
42
//-----------------------------------------------------------------------------
26
43
SparsityPattern::~SparsityPattern()
27
44
{
28
 
  //Do nothing
 
45
  // Do nothing
29
46
}
30
47
//-----------------------------------------------------------------------------
31
48
void SparsityPattern::init(uint rank, const uint* dims)
51
68
  initRange();
52
69
}
53
70
//-----------------------------------------------------------------------------
54
 
void SparsityPattern::insert(const uint* num_rows, const uint * const * rows)
 
71
void SparsityPattern::insert(uint m, const uint* rows, uint n, const uint* cols)
55
72
{ 
56
 
  for (unsigned int i = 0; i<num_rows[0];++i)
57
 
    for (unsigned int j = 0; j<num_rows[1];++j)
58
 
      sparsity_pattern[rows[0][i]].insert(rows[1][j]);
 
73
  for (uint i = 0; i < m; ++i)
 
74
    for (uint j = 0; j < n; ++j)
 
75
      sparsity_pattern[rows[i]].insert(cols[j]);
59
76
}
60
77
//-----------------------------------------------------------------------------
61
78
void SparsityPattern::pinsert(const uint* num_rows, const uint * const * rows)
171
188
  for(uint p=0; p<num_procs; ++p)
172
189
    range[p+1] = range[p] + dim[0]/num_procs + ((dim[0]%num_procs) > p ? 1 : 0);
173
190
}
 
191
//-----------------------------------------------------------------------------