Sierra Toolkit  Version of the Day
AggregateLinearSystem.cpp
1 /*------------------------------------------------------------------------*/
2 /* Copyright 2010 Sandia Corporation. */
3 /* Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive */
4 /* license for use of this work by or on behalf of the U.S. Government. */
5 /* Export of this program may require a license from the */
6 /* United States Government. */
7 /*------------------------------------------------------------------------*/
8 
9 
10 #include <stk_linsys/AggregateLinearSystem.hpp>
11 #include <stk_linsys/LinearSystem.hpp>
12 #include <stk_linsys/ImplDetails.hpp>
13 #include <stk_mesh/base/GetBuckets.hpp>
14 
15 #include <stk_linsys/LinsysFunctions.hpp>
16 
17 namespace stk_classic {
18 namespace linsys {
19 
20 AggregateLinearSystem::AggregateLinearSystem(MPI_Comm comm, fei::SharedPtr<fei::Factory> factory, size_t num_matrices, size_t num_rhsvecs)
21  : m_fei_factory(factory),
22  m_linear_system(comm, factory),
23  m_matrices(num_matrices),
24  m_rhsvecs(num_rhsvecs)
25 {
26 }
27 
29 {
30 }
31 
32 void
33 AggregateLinearSystem::set_parameters(Teuchos::ParameterList& paramlist)
34 {
35  m_linear_system.set_parameters(paramlist);
36 }
37 
38 void
39 AggregateLinearSystem::set_num_matrices_rhsvecs(size_t num_matrices, size_t num_rhsvecs)
40 {
41  m_matrices.resize(num_matrices);
42  m_rhsvecs.resize(num_rhsvecs);
43 }
44 
45 void
47 {
48  m_linear_system.synchronize_mappings_and_structure();
49 }
50 
51 void
53 {
54  m_linear_system.create_fei_LinearSystem();
55 
56  fei::SharedPtr<fei::MatrixGraph> mgraph = m_linear_system.get_fei_MatrixGraph();
57 
58  for(size_t i=0; i<m_matrices.size(); ++i) {
59  m_matrices[i] = m_fei_factory->createMatrix(mgraph);
60  }
61 
62  bool is_soln_vec = false;
63  for(size_t i=0; i<m_rhsvecs.size(); ++i) {
64  m_rhsvecs[i] = m_fei_factory->createVector(mgraph,is_soln_vec);
65  }
66 }
67 
68 fei::SharedPtr<fei::Matrix>
70 {
71  if (index >= m_matrices.size()) {
72  throw std::runtime_error("stk_classic::linsys::AggregateLinearSystem::get_matrix ERROR, index out of range.");
73  }
74 
75  return m_matrices[index];
76 }
77 
78 fei::SharedPtr<fei::Vector>
80 {
81  if (index >= m_rhsvecs.size()) {
82  throw std::runtime_error("stk_classic::linsys::AggregateLinearSystem::get_rhsvec ERROR, index out of range.");
83  }
84 
85  return m_rhsvecs[index];
86 }
87 
88 void
89 AggregateLinearSystem::aggregate_system(const std::vector<double>& mat_scalars,
90  const std::vector<double>& rhs_scalars)
91 {
92  if (mat_scalars.size() != m_matrices.size()) {
93  throw std::runtime_error("stk_classic::linsys::AggregateLinearSystem::aggregate_system ERROR, mat_scalars.size() != m_matrices.size().");
94  }
95 
96  if (rhs_scalars.size() != m_rhsvecs.size()) {
97  throw std::runtime_error("stk_classic::linsys::AggregateLinearSystem::aggregate_system ERROR, rhs_scalars.size() != m_rhsvecs.size().");
98  }
99 
100  fei::SharedPtr<fei::LinearSystem> fei_linsys = m_linear_system.get_fei_LinearSystem();
101  fei::SharedPtr<fei::Matrix> matrix = fei_linsys->getMatrix();
102  fei::SharedPtr<fei::Vector> rhsvec = fei_linsys->getRHS();
103 
104  matrix->gatherFromOverlap();
105  matrix->putScalar(0.0);
106 
107  for(size_t i=0; i<m_matrices.size(); ++i) {
108  stk_classic::linsys::add_matrix_to_matrix(mat_scalars[i], *m_matrices[i], *matrix);
109  }
110 
111  rhsvec->gatherFromOverlap();
112  rhsvec->putScalar(0.0);
113 
114  for(size_t i=0; i<m_rhsvecs.size(); ++i) {
115  stk_classic::linsys::add_vector_to_vector(rhs_scalars[i], *m_rhsvecs[i], *rhsvec);
116  }
117 }
118 
119 void
121 {
122  m_linear_system.finalize_assembly();
123 }
124 
125 const DofMapper&
127 {
128  return m_linear_system.get_DofMapper();
129 }
130 
131 DofMapper&
133 {
134  return m_linear_system.get_DofMapper();
135 }
136 
137 void
138 AggregateLinearSystem::reset_to_zero()
139 {
140  for(size_t i=0; i<m_matrices.size(); ++i) {
141  if (m_matrices[i].get() != NULL) m_matrices[i]->putScalar(0);
142  }
143  for(size_t i=0; i<m_rhsvecs.size(); ++i) {
144  if (m_rhsvecs[i].get() != NULL) m_rhsvecs[i]->putScalar(0);
145  }
146 }
147 
148 const fei::SharedPtr<fei::MatrixGraph>
150 {
151  return m_linear_system.get_fei_MatrixGraph();
152 }
153 
154 fei::SharedPtr<fei::MatrixGraph>
156 {
157  return m_linear_system.get_fei_MatrixGraph();
158 }
159 
160 const fei::SharedPtr<fei::LinearSystem>
162 {
163  return m_linear_system.get_fei_LinearSystem();
164 }
165 
166 fei::SharedPtr<fei::LinearSystem>
168 {
169  return m_linear_system.get_fei_LinearSystem();
170 }
171 
172 void
173 AggregateLinearSystem::write_files(const std::string& base_name) const
174 {
175  for(size_t i=0; i<m_matrices.size(); ++i) {
176  if (m_matrices[i].get() == NULL) continue;
177  std::ostringstream ossA;
178  ossA << "A_" << base_name << ".mat"<<i<<".mtx";
179  std::string Aname = ossA.str();
180  m_matrices[i]->writeToFile(Aname.c_str());
181  }
182  for(size_t i=0; i<m_rhsvecs.size(); ++i) {
183  if (m_rhsvecs[i].get() == NULL) continue;
184  std::ostringstream ossb;
185  ossb << "b_" << base_name << ".vec"<<i<<".vec";
186  std::string bname = ossb.str();
187  m_rhsvecs[i]->writeToFile(bname.c_str());
188  }
189 }
190 
191 int
192 AggregateLinearSystem::solve(int &status, const Teuchos::ParameterList & params )
193 {
194  return m_linear_system.solve(status, params);
195 }
196 
197 }//namespace linsys
198 }//namespace stk_classic
199 
const fei::SharedPtr< fei::MatrixGraph > get_fei_MatrixGraph() const
fei::SharedPtr< fei::Vector > get_rhsvec(size_t index)
int solve(int &status, const Teuchos::ParameterList &params)
void aggregate_system(const std::vector< double > &mat_scalars, const std::vector< double > &rhs_scalars)
const fei::SharedPtr< fei::LinearSystem > get_fei_LinearSystem() const
fei::SharedPtr< fei::Matrix > get_matrix(size_t index)
const fei::SharedPtr< fei::LinearSystem > get_fei_LinearSystem() const
const DofMapper & get_DofMapper() const
AggregateLinearSystem(MPI_Comm comm, fei::SharedPtr< fei::Factory > factory, size_t num_matrices=1, size_t num_rhsvecs=1)
void add_matrix_to_matrix(double scalar, const fei::Matrix &src_matrix, fei::Matrix &dest_matrix)
const fei::SharedPtr< fei::MatrixGraph > get_fei_MatrixGraph() const
int solve(int &status, const Teuchos::ParameterList &params)
Sierra Toolkit.
void add_vector_to_vector(double scalar, const fei::Vector &src_vector, fei::Vector &dest_vector)
void set_num_matrices_rhsvecs(size_t num_matrices, size_t num_rhsvecs)