← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to sandbox/la/hypre/main.cpp

  • Committer: Anders Logg
  • Date: 2008-05-22 06:46:27 UTC
  • mto: (2668.9.6 trunk)
  • mto: This revision was merged to the branch mainline in revision 2670.
  • Revision ID: logg@simula.no-20080522064627-rahrtqj3jqxnlc2j
More work on BCs

Show diffs side-by-side

added added

removed removed

Lines of Context:
sandbox/la/hypre/main.cpp
4
4
// Simple test program for solving a singular system
5
5
// with GMRES + AMG (PETSc + Hypre BoomerAMG)
6
6
//
7
 
// This system has a zero eigenvalue with corresponding
8
 
// eigenvector v = (1, 1, 1, ...). This gives a compatibility
9
 
// condition (corresponding to the integral of the right-hand
10
 
// side being zero for the Neumann problem).
11
 
//
12
 
// To solve the linear system, we must therefore either make the
13
 
// system nonsingular by adding a constraint, like zero average for x,
14
 
// or modify the right-hand side to satisfy the compatibility
15
 
// condition.
 
7
// Singular solver converges much better
16
8
 
17
9
#include <dolfin.h>
18
10
#include "Poisson.h"
62
54
  // Solve linear system using ordinary linear solver
63
55
  LinearSolver s0(gmres, amg);
64
56
  s0.solve(A, x, b);
65
 
  cout << "Residual: " << residual(A, x, b) << endl;
66
57
 
67
58
  // Solve linear system using special singular solver
68
59
  SingularSolver s1(gmres, amg);
69
60
  s1.solve(A, x, b);
70
 
  cout << "Residual: " << residual(A, x, b) << endl;
71
 
 
72
 
  // Solve modified linear system using ordinary linear solver
 
61
 
 
62
 
 
63
  // Compute constant, c = int_\Omega b and subtract it from b  
 
64
  real c = 0.0; 
 
65
  for (unsigned int i=0; i<b.size(); i++) {
 
66
    c += b[i]; 
 
67
  }
 
68
  c /= b.size(); 
 
69
  Vector bhat(b.size()); 
 
70
  bhat = c; 
 
71
  b -= bhat; 
 
72
 
 
73
  // Solve linear system using ordinary linear solver
73
74
  LinearSolver s2(gmres, amg);
74
 
  normalize(b, average);
75
75
  s2.solve(A, x, b);
76
 
  cout << "Residual: " << residual(A, x, b) << endl;
 
76
 
77
77
 
78
78
  return 0;
79
79
}