Sierra Toolkit Version of the Day
AggregateLinearSystem.cpp
00001 /*------------------------------------------------------------------------*/
00002 /*                 Copyright 2010 Sandia Corporation.                     */
00003 /*  Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive   */
00004 /*  license for use of this work by or on behalf of the U.S. Government.  */
00005 /*  Export of this program may require a license from the                 */
00006 /*  United States Government.                                             */
00007 /*------------------------------------------------------------------------*/
00008 
00009 
00010 #include <stk_linsys/AggregateLinearSystem.hpp>
00011 #include <stk_linsys/LinearSystem.hpp>
00012 #include <stk_linsys/ImplDetails.hpp>
00013 #include <stk_mesh/base/GetBuckets.hpp>
00014 
00015 #include <stk_linsys/LinsysFunctions.hpp>
00016 
00017 namespace stk {
00018 namespace linsys {
00019 
00020 AggregateLinearSystem::AggregateLinearSystem(MPI_Comm comm, fei::SharedPtr<fei::Factory> factory, size_t num_matrices, size_t num_rhsvecs)
00021  : m_fei_factory(factory),
00022    m_linear_system(comm, factory),
00023    m_matrices(num_matrices),
00024    m_rhsvecs(num_rhsvecs)
00025 {
00026 }
00027 
00028 AggregateLinearSystem::~AggregateLinearSystem()
00029 {
00030 }
00031 
00032 void
00033 AggregateLinearSystem::set_parameters(Teuchos::ParameterList& paramlist)
00034 {
00035   m_linear_system.set_parameters(paramlist);
00036 }
00037 
00038 void
00039 AggregateLinearSystem::set_num_matrices_rhsvecs(size_t num_matrices, size_t num_rhsvecs)
00040 {
00041   m_matrices.resize(num_matrices);
00042   m_rhsvecs.resize(num_rhsvecs);
00043 }
00044 
00045 void
00046 AggregateLinearSystem::synchronize_mappings_and_structure()
00047 {
00048   m_linear_system.synchronize_mappings_and_structure();
00049 }
00050 
00051 void
00052 AggregateLinearSystem::create_fei_LinearSystem()
00053 {
00054   m_linear_system.create_fei_LinearSystem();
00055 
00056   fei::SharedPtr<fei::MatrixGraph> mgraph = m_linear_system.get_fei_MatrixGraph();
00057 
00058   for(size_t i=0; i<m_matrices.size(); ++i) {
00059     m_matrices[i] = m_fei_factory->createMatrix(mgraph);
00060   }
00061 
00062   bool is_soln_vec = false;
00063   for(size_t i=0; i<m_rhsvecs.size(); ++i) {
00064     m_rhsvecs[i] = m_fei_factory->createVector(mgraph,is_soln_vec);
00065   }
00066 }
00067 
00068 fei::SharedPtr<fei::Matrix>
00069 AggregateLinearSystem::get_matrix(size_t index)
00070 {
00071   if (index >= m_matrices.size()) {
00072     throw std::runtime_error("stk::linsys::AggregateLinearSystem::get_matrix ERROR, index out of range.");
00073   }
00074 
00075   return m_matrices[index];
00076 }
00077 
00078 fei::SharedPtr<fei::Vector>
00079 AggregateLinearSystem::get_rhsvec(size_t index)
00080 {
00081   if (index >= m_rhsvecs.size()) {
00082     throw std::runtime_error("stk::linsys::AggregateLinearSystem::get_rhsvec ERROR, index out of range.");
00083   }
00084 
00085   return m_rhsvecs[index];
00086 }
00087 
00088 void
00089 AggregateLinearSystem::aggregate_system(const std::vector<double>& mat_scalars,
00090                                         const std::vector<double>& rhs_scalars)
00091 {
00092   if (mat_scalars.size() != m_matrices.size()) {
00093     throw std::runtime_error("stk::linsys::AggregateLinearSystem::aggregate_system ERROR, mat_scalars.size() != m_matrices.size().");
00094   }
00095 
00096   if (rhs_scalars.size() != m_rhsvecs.size()) {
00097     throw std::runtime_error("stk::linsys::AggregateLinearSystem::aggregate_system ERROR, rhs_scalars.size() != m_rhsvecs.size().");
00098   }
00099 
00100   fei::SharedPtr<fei::LinearSystem> fei_linsys = m_linear_system.get_fei_LinearSystem();
00101   fei::SharedPtr<fei::Matrix> matrix = fei_linsys->getMatrix();
00102   fei::SharedPtr<fei::Vector> rhsvec = fei_linsys->getRHS();
00103 
00104   matrix->gatherFromOverlap();
00105   matrix->putScalar(0.0);
00106 
00107   for(size_t i=0; i<m_matrices.size(); ++i) {
00108     stk::linsys::add_matrix_to_matrix(mat_scalars[i], *m_matrices[i], *matrix);
00109   }
00110 
00111   rhsvec->gatherFromOverlap();
00112   rhsvec->putScalar(0.0);
00113 
00114   for(size_t i=0; i<m_rhsvecs.size(); ++i) {
00115     stk::linsys::add_vector_to_vector(rhs_scalars[i], *m_rhsvecs[i], *rhsvec);
00116   }
00117 }
00118 
00119 void
00120 AggregateLinearSystem::finalize_assembly()
00121 {
00122   m_linear_system.finalize_assembly();
00123 }
00124 
00125 const DofMapper&
00126 AggregateLinearSystem::get_DofMapper() const
00127 {
00128   return m_linear_system.get_DofMapper();
00129 }
00130 
00131 DofMapper&
00132 AggregateLinearSystem::get_DofMapper()
00133 {
00134   return m_linear_system.get_DofMapper();
00135 }
00136 
00137 void
00138 AggregateLinearSystem::reset_to_zero()
00139 {
00140   for(size_t i=0; i<m_matrices.size(); ++i) {
00141     if (m_matrices[i].get() != NULL) m_matrices[i]->putScalar(0);
00142   }
00143   for(size_t i=0; i<m_rhsvecs.size(); ++i) {
00144     if (m_rhsvecs[i].get() != NULL) m_rhsvecs[i]->putScalar(0);
00145   }
00146 }
00147 
00148 const fei::SharedPtr<fei::MatrixGraph>
00149 AggregateLinearSystem::get_fei_MatrixGraph() const
00150 {
00151   return m_linear_system.get_fei_MatrixGraph();
00152 }
00153 
00154 fei::SharedPtr<fei::MatrixGraph>
00155 AggregateLinearSystem::get_fei_MatrixGraph()
00156 {
00157   return m_linear_system.get_fei_MatrixGraph();
00158 }
00159 
00160 const fei::SharedPtr<fei::LinearSystem>
00161 AggregateLinearSystem::get_fei_LinearSystem() const
00162 {
00163   return m_linear_system.get_fei_LinearSystem();
00164 }
00165 
00166 fei::SharedPtr<fei::LinearSystem>
00167 AggregateLinearSystem::get_fei_LinearSystem()
00168 {
00169   return m_linear_system.get_fei_LinearSystem();
00170 }
00171 
00172 void
00173 AggregateLinearSystem::write_files(const std::string& base_name) const
00174 {
00175   for(size_t i=0; i<m_matrices.size(); ++i) {
00176     if (m_matrices[i].get() == NULL) continue;
00177     std::ostringstream ossA;
00178     ossA << "A_" << base_name << ".mat"<<i<<".mtx";
00179     std::string Aname = ossA.str();
00180     m_matrices[i]->writeToFile(Aname.c_str());
00181   }
00182   for(size_t i=0; i<m_rhsvecs.size(); ++i) {
00183     if (m_rhsvecs[i].get() == NULL) continue;
00184     std::ostringstream ossb;
00185     ossb << "b_" << base_name << ".vec"<<i<<".vec";
00186     std::string bname = ossb.str();
00187     m_rhsvecs[i]->writeToFile(bname.c_str());
00188   }
00189 }
00190 
00191 int
00192 AggregateLinearSystem::solve(int &status, const Teuchos::ParameterList & params )
00193 {
00194   return m_linear_system.solve(status, params);
00195 }
00196 
00197 }//namespace linsys
00198 }//namespace stk
00199 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends