~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to dolfin/la/BlockMatrix.cpp

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2014-09-22 14:35:34 UTC
  • mfrom: (1.1.17) (19.1.23 sid)
  • Revision ID: package-import@ubuntu.com-20140922143534-0yi89jyuqbgdxwm9
Tags: 1.4.0+dfsg-4
* debian/control: Disable libcgal-dev on i386, mipsel and sparc.
* debian/rules: Remove bad directives in pkg-config file dolfin.pc
  (closes: #760658).
* Remove debian/libdolfin-dev.lintian-overrides.

Show diffs side-by-side

added added

removed removed

Lines of Context:
22
22
// Last changed: 2012-03-15
23
23
 
24
24
#include <iostream>
25
 
#include <boost/shared_ptr.hpp>
 
25
#include <memory>
26
26
#include <dolfin/common/Timer.h>
27
27
#include <dolfin/common/NoDeleter.h>
28
28
#include "dolfin/common/utils.h"
35
35
using namespace dolfin;
36
36
 
37
37
//-----------------------------------------------------------------------------
38
 
BlockMatrix::BlockMatrix(std::size_t m, std::size_t n) : matrices(boost::extents[m][n])
 
38
BlockMatrix::BlockMatrix(std::size_t m, std::size_t n)
 
39
  : matrices(boost::extents[m][n])
39
40
{
40
41
  for (std::size_t i = 0; i < m; i++)
41
42
    for (std::size_t j = 0; j < n; j++)
47
48
  // Do nothing
48
49
}
49
50
//-----------------------------------------------------------------------------
50
 
void BlockMatrix::set_block(std::size_t i, std::size_t j, boost::shared_ptr<GenericMatrix> m)
 
51
void BlockMatrix::set_block(std::size_t i, std::size_t j,
 
52
                            std::shared_ptr<GenericMatrix> m)
51
53
{
52
54
  dolfin_assert(i < matrices.shape()[0]);
53
55
  dolfin_assert(j < matrices.shape()[1]);
54
56
  matrices[i][j] = m;
55
57
}
56
58
//-----------------------------------------------------------------------------
57
 
const boost::shared_ptr<GenericMatrix> BlockMatrix::get_block(std::size_t i, std::size_t j) const
 
59
std::shared_ptr<const GenericMatrix>
 
60
BlockMatrix::get_block(std::size_t i, std::size_t j) const
58
61
{
59
62
  dolfin_assert(i < matrices.shape()[0]);
60
63
  dolfin_assert(j < matrices.shape()[1]);
61
64
  return matrices[i][j];
62
65
}
63
66
//-----------------------------------------------------------------------------
64
 
boost::shared_ptr<GenericMatrix> BlockMatrix::get_block(std::size_t i, std::size_t j)
 
67
std::shared_ptr<GenericMatrix> BlockMatrix::get_block(std::size_t i,
 
68
                                                        std::size_t j)
65
69
{
66
70
  dolfin_assert(i < matrices.shape()[0]);
67
71
  dolfin_assert(j < matrices.shape()[1]);
76
80
//-----------------------------------------------------------------------------
77
81
void BlockMatrix::zero()
78
82
{
79
 
  for(std::size_t i = 0; i < matrices.shape()[0]; i++)
80
 
    for(std::size_t j = 0; j < matrices.shape()[1]; j++)
 
83
  for (std::size_t i = 0; i < matrices.shape()[0]; i++)
 
84
    for (std::size_t j = 0; j < matrices.shape()[1]; j++)
81
85
      matrices[i][j]->zero();
82
86
}
83
87
//-----------------------------------------------------------------------------
100
104
    {
101
105
      for (std::size_t j = 0; i < matrices.shape()[1]; j++)
102
106
      {
103
 
        s << "  BlockMatrix (" << i << ", " << j << ")" << std::endl << std::endl;
 
107
        s << "  BlockMatrix (" << i << ", " << j << ")" << std::endl
 
108
          << std::endl;
104
109
        s << indent(indent(matrices[i][j]->str(true))) << std::endl;
105
110
      }
106
111
    }
107
112
  }
108
113
  else
109
 
    s << "<BlockMatrix containing " << matrices.shape()[0] << " x " << matrices.shape()[1] << " blocks>";
 
114
  {
 
115
    s << "<BlockMatrix containing " << matrices.shape()[0] << " x "
 
116
      << matrices.shape()[1] << " blocks>";
 
117
  }
110
118
 
111
119
  return s.str();
112
120
}
123
131
 
124
132
  // Create tempory vector
125
133
  dolfin_assert(matrices[0][0]);
126
 
  boost::shared_ptr<GenericVector>
127
 
    z_tmp = matrices[0][0]->factory().create_vector();
128
134
 
129
135
  // Loop over block rows
130
136
  for(std::size_t row = 0; row < matrices.shape()[0]; row++)
132
138
    // RHS sub-vector
133
139
    GenericVector& _y = *(y.get_block(row));
134
140
 
 
141
    const GenericMatrix& _matA = *matrices[row][0];
 
142
 
135
143
    // Resize y and zero
136
 
    dolfin_assert(matrices[row][0]);
137
 
    const GenericMatrix& _A = *matrices[row][0];
138
 
    _A.resize(_y, 0);
 
144
    if (_y.empty())
 
145
      _matA.init_vector(_y, 0);
139
146
    _y.zero();
140
147
 
141
 
    // Resize z_tmp and zero
142
 
    z_tmp->resize(_y.size());
143
 
    _A.resize(*z_tmp, 0);
144
 
 
145
148
    // Loop over block columns
 
149
    std::shared_ptr<GenericVector>
 
150
      z_tmp = _matA.factory().create_vector();
146
151
    for(std::size_t col = 0; col < matrices.shape()[1]; ++col)
147
152
    {
148
153
      const GenericVector& _x = *(x.get_block(col));
153
158
  }
154
159
}
155
160
//-----------------------------------------------------------------------------
156
 
boost::shared_ptr<GenericMatrix> BlockMatrix::schur_approximation(bool symmetry) const
 
161
std::shared_ptr<GenericMatrix> BlockMatrix::schur_approximation(bool symmetry) const
157
162
{
158
163
  // Currently returns [diag(C * diag(A)^-1 * B) - D]
159
164
  if (!symmetry)
169
174
  GenericMatrix &C = *matrices[1][0];
170
175
  GenericMatrix &D = *matrices[1][1];
171
176
 
172
 
  boost::shared_ptr<GenericMatrix> S(D.copy());
 
177
  std::shared_ptr<GenericMatrix> S(D.copy());
173
178
 
174
179
  std::vector<std::size_t> cols_i;
175
180
  std::vector<double> vals_i;