matrix_algorithms.hpp 2.56 KB
Newer Older
1
// Software License for MTL
2
//
3
4
5
6
7
// Copyright (c) 2007 The Trustees of Indiana University.
//               2008 Dresden University of Technology and the Trustees of Indiana University.
//               2010 SimuNova UG (haftungsbeschränkt), www.simunova.com.
// All rights reserved.
// Authors: Peter Gottschling and Andrew Lumsdaine
8
//
9
// This file is part of the Matrix Template Library
10
//
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
// See also license.mtl.txt in the distribution.
//
// Algorithm inspired by Nick Vannieuwenhoven, written by Cornelius Steinhardt



#ifndef MTL_MATRIX_ALGORITHMS_INCLUDE
#define MTL_MATRIX_ALGORITHMS_INCLUDE

#include <boost/numeric/mtl/interface/vpt.hpp>
#include <boost/numeric/mtl/mtl.hpp>
#include <boost/numeric/mtl/matrix/element.hpp>
#include <boost/numeric/mtl/matrix/element_structure.hpp>

#include <iostream>

namespace mtl {  namespace matrix {

29
30
31
/// Construct the sparse data structure from an elementstructure
template< typename ElementStructure, typename Matrix, typename Vector>
void assemble_compressed(const ElementStructure& es, Matrix& A, Vector& order)
32
33
34
35
36
37
38
39
40
41
42
{

  	typedef typename ElementStructure::element_type::value_type   value_type;
 	typedef typename ElementStructure::element_iterator           iterator;
 	typedef typename ElementStructure::element_type               element_type;
 	typedef typename element_type::index_type                     index_type;
 	typedef typename element_type::matrix_type                    matrix_type;
	typedef typename matrix_type::size_type                       size_type;
	A.change_dim(es.get_total_vars(), es.get_total_vars());
	set_to_zero(A);
	value_type zero(0);
43

44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
	{//start inserterblock
	  mtl::matrix::inserter<Matrix, mtl::operations::update_plus<value_type> >  ins(A);
	  for(iterator it = es.element_begin(); it != es.element_end(); ++it) {
		element_type& element = *it;
		const index_type& idx = element.get_indices();
		matrix_type& values = element.get_values();
		for(int i = 0; i < element.nb_vars(); ++i) {
			for(int j = 0; j < element.nb_vars(); ++j) {
				if(values(i,j) != zero) {
					ins[size_type(order(idx(i)))][size_type(order(idx(j)))] << values(i,j);
				}
			}
		}
	  }
	}//end inserterblock
}


62
63
64
/// Construct the sparse data structure from an elementstructure
template< typename ElementStructure, typename Matrix>
void assemble_compressed(const ElementStructure& es, Matrix& A)
65
66
67
68
69
70
71
72
73
74
{
  typedef typename  ElementStructure::element_type::matrix_type::size_type size_type;
  mtl::dense_vector<size_type> ident(es.get_total_vars());
  iota(ident);
  assemble_compressed(es, A, ident);
}

}}//end namespace mtl

#endif // MTL_MATRIX_ALGORITHMS_INCLUDE