← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to dolfin/la/SingularSolver.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/SingularSolver.cpp
 
1
// Copyright (C) 2008 Anders Logg.
 
2
// Licensed under the GNU LGPL Version 2.1.
 
3
//
 
4
// First added:  2008-05-15
 
5
// Last changed: 2008-05-15
 
6
 
 
7
#include <dolfin/common/Array.h>
 
8
#include "LinearAlgebraFactory.h"
 
9
#include "GenericMatrix.h"
 
10
#include "GenericVector.h"
 
11
#include "SparsityPattern.h"
 
12
#include "SingularSolver.h"
 
13
 
 
14
using namespace dolfin;
 
15
 
 
16
//-----------------------------------------------------------------------------
 
17
SingularSolver::SingularSolver(SolverType solver_type,
 
18
                               PreconditionerType pc_type)
 
19
  : Parametrized(), linear_solver(solver_type, pc_type), B(0), y(0), c(0)
 
20
{
 
21
  // Set parameters for linear solver
 
22
  linear_solver.set("parent", *this);
 
23
}
 
24
//-----------------------------------------------------------------------------
 
25
SingularSolver::~SingularSolver()
 
26
{
 
27
  delete B;
 
28
  delete y;
 
29
  delete c;
 
30
}
 
31
//-----------------------------------------------------------------------------
 
32
dolfin::uint SingularSolver::solve(const GenericMatrix& A,
 
33
                                   GenericVector& x, const GenericVector& b)
 
34
{
 
35
  message("Solving singular system...");
 
36
 
 
37
  // Initialize data structures for extended system
 
38
  init(A);
 
39
 
 
40
  // Create extended system
 
41
  create(A, b, 0);
 
42
  
 
43
  // Solve extended system
 
44
  const uint num_iterations = linear_solver.solve(*B, *y, *c);
 
45
 
 
46
  // Extract solution
 
47
  x.init(y->size() - 1);
 
48
  real* vals = new real[y->size()];
 
49
  y->get(vals);
 
50
  x.set(vals);
 
51
  delete [] vals;
 
52
 
 
53
  return num_iterations;
 
54
}
 
55
//-----------------------------------------------------------------------------
 
56
dolfin::uint SingularSolver::solve(const GenericMatrix& A,
 
57
                                   GenericVector& x, const GenericVector& b,
 
58
                                   const GenericMatrix& M)
 
59
{
 
60
  message("Solving singular system...");
 
61
 
 
62
  // Initialize data structures for extended system
 
63
  init(A);
 
64
 
 
65
  // Create extended system
 
66
  create(A, b, &M);
 
67
 
 
68
  // Solve extended system
 
69
  const uint num_iterations = linear_solver.solve(*B, *y, *c);
 
70
 
 
71
  // Extract solution
 
72
  x.init(y->size() - 1);
 
73
  real* vals = new real[y->size()];
 
74
  y->get(vals);
 
75
  x.set(vals);
 
76
  delete [] vals;
 
77
 
 
78
  return num_iterations;
 
79
}
 
80
//-----------------------------------------------------------------------------
 
81
void SingularSolver::init(const GenericMatrix& A)
 
82
{
 
83
  // Check size of system
 
84
  if (A.size(0) != A.size(1))
 
85
    error("Matrix must be square.");
 
86
  if (A.size(0) == 0)
 
87
    error("Matrix size must be non-zero.");
 
88
 
 
89
  // Get dimension
 
90
  const uint N = A.size(0);
 
91
 
 
92
  // Check if we have already initialized system
 
93
  if (B && B->size(0) == N + 1 && B->size(1) == N + 1)
 
94
    return;
 
95
 
 
96
  cout << "Initializing" << endl;
 
97
  
 
98
  // Delete any old data
 
99
  delete B;
 
100
  delete y;
 
101
  delete c;
 
102
 
 
103
  // Create sparsity pattern for B
 
104
  SparsityPattern s(N + 1, N + 1);
 
105
 
 
106
  // Copy sparsity pattern for A and last column
 
107
  Array<uint> columns;
 
108
  Array<real> dummy;
 
109
  for (uint i = 0; i < N; i++)
 
110
  {
 
111
    // Get row
 
112
    A.getrow(i, columns, dummy);
 
113
 
 
114
    // Copy columns to array
 
115
    const uint num_cols = columns.size() + 1;
 
116
    uint* cols = new uint[num_cols];
 
117
    for (uint j = 0; j < columns.size(); j++)
 
118
      cols[j] = columns[j];
 
119
 
 
120
    // Add last entry
 
121
    cols[num_cols - 1] = N;
 
122
 
 
123
    // Insert into sparsity pattern
 
124
    s.insert(1, &i, num_cols, cols);
 
125
 
 
126
    // Delete temporary array
 
127
    delete [] cols;
 
128
  }
 
129
 
 
130
  // Add last row
 
131
  const uint num_cols = N;
 
132
  uint* cols = new uint[num_cols];
 
133
  for (uint j = 0; j < num_cols; j++)
 
134
    cols[j] = j;
 
135
  const uint row = N;
 
136
  s.insert(1, &row, num_cols, cols);
 
137
  delete [] cols;
 
138
 
 
139
  // Create matrix and vector
 
140
  B = A.factory().createMatrix();
 
141
  y = A.factory().createVector();
 
142
  c = A.factory().createVector();
 
143
  B->init(s);
 
144
  y->init(N + 1);
 
145
  c->init(N + 1);
 
146
}
 
147
//-----------------------------------------------------------------------------
 
148
void SingularSolver::create(const GenericMatrix& A, const GenericVector& b,
 
149
                            const GenericMatrix* M)
 
150
{
 
151
  dolfin_assert(B);
 
152
  dolfin_assert(c);
 
153
 
 
154
  // Reset matrix
 
155
  B->zero();
 
156
 
 
157
  // Copy rows from A into B
 
158
  const uint N = A.size(0);
 
159
  Array<uint> columns;
 
160
  Array<real> values;
 
161
  for (uint i = 0; i < N; i++)
 
162
  {
 
163
    A.getrow(i, columns, values);
 
164
    B->setrow(i, columns, values);
 
165
  }
 
166
 
 
167
  // Compute lumped mass matrix
 
168
  columns.resize(N);
 
169
  values.resize(N);
 
170
  if (M)
 
171
  {
 
172
    GenericVector* ones = A.factory().createVector();
 
173
    GenericVector* z = A.factory().createVector();
 
174
    ones->init(N);
 
175
    z->init(N);
 
176
    *ones = 1.0;
 
177
    A.mult(*ones, *z);
 
178
    for (uint i = 0; i < N; i++)
 
179
    {
 
180
      columns[i] = i;
 
181
      values[i] = (*z)[i];
 
182
    }
 
183
    delete ones;
 
184
    delete z;
 
185
  }
 
186
  else
 
187
  {
 
188
    for (uint i = 0; i < N; i++)
 
189
    {
 
190
      columns[i] = i;
 
191
      values[i] = 1.0;
 
192
    }
 
193
  }
 
194
 
 
195
  // Add last row
 
196
  B->setrow(N, columns, values);
 
197
 
 
198
  // Add last column
 
199
  for (uint i = 0; i < N; i++)
 
200
    B->set(&values[i], 1, &i, 1, &N);
 
201
 
 
202
  // Copy values from b into c
 
203
  real* vals = new real[N + 1];
 
204
  b.get(vals);
 
205
  vals[N] = 0.0;
 
206
  c->set(vals);
 
207
  delete [] vals;
 
208
 
 
209
  // Apply changes
 
210
  B->apply();
 
211
  c->apply();
 
212
}
 
213
//-----------------------------------------------------------------------------