/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /trunk/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5147 by caltinay, Wed Aug 20 06:19:20 2014 UTC revision 5148 by caltinay, Mon Sep 15 01:25:23 2014 UTC
# Line 20  Line 20 
20  #include <pasowrap/SystemMatrixAdapter.h>  #include <pasowrap/SystemMatrixAdapter.h>
21  #include <pasowrap/TransportProblemAdapter.h>  #include <pasowrap/TransportProblemAdapter.h>
22  #include <ripley/domainhelpers.h>  #include <ripley/domainhelpers.h>
23    #include <ripley/RipleySystemMatrix.h>
24    
25  #include <iomanip>  #include <iomanip>
26    
27  namespace bp = boost::python;  namespace bp = boost::python;
28    
29  using namespace std;  using namespace std;
 using paso::SystemMatrixAdapter;  
30  using paso::TransportProblemAdapter;  using paso::TransportProblemAdapter;
31    
32  namespace ripley {  namespace ripley {
# Line 36  void tupleListToMap(DataMap& mapping, co Line 36  void tupleListToMap(DataMap& mapping, co
36      for (int i = 0; i < len(list); i++) {      for (int i = 0; i < len(list); i++) {
37          if (!bp::extract<bp::tuple>(list[i]).check())          if (!bp::extract<bp::tuple>(list[i]).check())
38              throw RipleyException("Passed in list contains objects"              throw RipleyException("Passed in list contains objects"
39                                      " other than tuples");                                    " other than tuples");
40          bp::tuple t = bp::extract<bp::tuple>(list[i]);          bp::tuple t = bp::extract<bp::tuple>(list[i]);
41          if (len(t) != 2 || !bp::extract<string>(t[0]).check() ||          if (len(t) != 2 || !bp::extract<string>(t[0]).check() ||
42                  !bp::extract<escript::Data>(t[1]).check())                  !bp::extract<escript::Data>(t[1]).check())
43              throw RipleyException("The passed in list must contain tuples"              throw RipleyException("The passed in list must contain tuples"
44                  " of the form (string, escript::data)");                  " of the form (string, escript::Data)");
45          mapping[bp::extract<string>(t[0])] = bp::extract<escript::Data>(t[1]);          mapping[bp::extract<string>(t[0])] = bp::extract<escript::Data>(t[1]);
46      }      }
47  }  }
# Line 763  void RipleyDomain::Print_Mesh_Info(bool Line 763  void RipleyDomain::Print_Mesh_Info(bool
763      }      }
764  }  }
765    
766  int RipleyDomain::getSystemMatrixTypeId(int solver, int preconditioner,  int RipleyDomain::getSystemMatrixTypeId(const bp::object& options) const
                                         int package, bool symmetry) const  
767  {  {
768      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,      const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
769              package, symmetry, m_mpiInfo);      int package = sb.getPackage();
770    
771        // use CUSP for single rank and supported solvers+preconditioners,
772        // PASO otherwise
773        if (package == escript::SO_DEFAULT) {
774            if (m_mpiInfo->size == 1) {
775                switch (sb.getSolverMethod()) {
776                    case escript::SO_DEFAULT:
777                    case escript::SO_METHOD_BICGSTAB:
778                    case escript::SO_METHOD_CGLS:
779                    case escript::SO_METHOD_GMRES:
780                    case escript::SO_METHOD_LSQR:
781                    case escript::SO_METHOD_PCG:
782                    case escript::SO_METHOD_PRES20:
783                        package = escript::SO_PACKAGE_CUSP;
784                        break;
785                    default:
786                        package = escript::SO_PACKAGE_PASO;
787                }
788                if (package == escript::SO_PACKAGE_CUSP) {
789                    if (sb.getPreconditioner() != escript::SO_PRECONDITIONER_NONE &&
790                            sb.getPreconditioner() != escript::SO_PRECONDITIONER_JACOBI) {
791                        package = escript::SO_PACKAGE_PASO;
792                    }
793                }
794            } else {
795                package = escript::SO_PACKAGE_PASO;
796            }
797        }
798    
799        if (package == escript::SO_PACKAGE_CUSP) {
800            if (m_mpiInfo->size > 1) {
801                throw RipleyException("CUSP matrices are not supported with more than one rank");
802            }
803            return (int)SMT_CUSP;
804        }
805    
806        // in all other cases we use PASO
807        return (int)SMT_PASO | paso::SystemMatrixAdapter::getSystemMatrixTypeId(
808                sb.getSolverMethod(), sb.getPreconditioner(), sb.getPackage(),
809                sb.isSymmetric(), m_mpiInfo);
810  }  }
811    
812  int RipleyDomain::getTransportTypeId(int solver, int preconditioner,  int RipleyDomain::getTransportTypeId(int solver, int preconditioner,
# Line 799  escript::ASM_ptr RipleyDomain::newSystem Line 838  escript::ASM_ptr RipleyDomain::newSystem
838          reduceColOrder=true;          reduceColOrder=true;
839      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
840          throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
841        // are block sizes identical?
842      if (type & MATRIX_FORMAT_TRILINOS_CRS)      if (row_blocksize != column_blocksize)
843          throw RipleyException("newSystemMatrix: Ripley does not support matrix format TRILINOS_CRS");          throw RipleyException("newSystemMatrix: row/column block sizes must be equal");
844        // are function spaces equal
845        if (reduceRowOrder != reduceColOrder)
846            throw RipleyException("newSystemMatrix: row/column function spaces must be equal");
847    
848      // generate matrix      // generate matrix
849      paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(reduceRowOrder,      //if (reduceRowOrder || reduceColOrder)
850                                                                reduceColOrder));      //    throw RipleyException("newSystemMatrix: reduced order not supported");
851      paso::SystemMatrix_ptr matrix(new paso::SystemMatrix(type, pattern,  
852              row_blocksize, column_blocksize, false));      if (type == (int)SMT_CUSP) {
853      paso::checkPasoError();          const int numMatrixRows = getNumDOF();
854      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,          escript::ASM_ptr sm(new SystemMatrix(m_mpiInfo, row_blocksize,
855                  row_functionspace, column_blocksize, column_functionspace));                      row_functionspace, numMatrixRows, getDiagonalIndices()));
856      return sma;          return sm;
857        } else if (type & (int)SMT_PASO) {
858            paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(
859                                                reduceRowOrder, reduceColOrder));
860            type -= (int)SMT_PASO;
861            paso::SystemMatrix_ptr matrix(new paso::SystemMatrix(type, pattern,
862                    row_blocksize, column_blocksize, false));
863            escript::ASM_ptr sm(new paso::SystemMatrixAdapter(matrix,
864                        row_blocksize, row_functionspace, column_blocksize,
865                        column_functionspace));
866            return sm;
867        } else {
868            throw RipleyException("newSystemMatrix: unknown matrix type ID");
869        }
870  }  }
871    
872  void RipleyDomain::addToSystem(escript::AbstractSystemMatrix& mat,  void RipleyDomain::addToSystem(escript::AbstractSystemMatrix& mat,
# Line 1087  void RipleyDomain::addToSystemMatrix(esc Line 1142  void RipleyDomain::addToSystemMatrix(esc
1142          paso::SystemMatrix_ptr S(sma->getPaso_SystemMatrix());          paso::SystemMatrix_ptr S(sma->getPaso_SystemMatrix());
1143          addToSystemMatrix(S, nodes, numEq, array);          addToSystemMatrix(S, nodes, numEq, array);
1144      } else {      } else {
1145          throw RipleyException("addToSystemMatrix: unknown system matrix type");          SystemMatrix* sm = dynamic_cast<SystemMatrix*>(mat);
1146            if (sm) {
1147                sm->add(nodes, array);
1148            } else {
1149                throw RipleyException("addToSystemMatrix: unknown system matrix type");
1150            }
1151      }      }
1152  }  }
1153    

Legend:
Removed from v.5147  
changed lines
  Added in v.5148

  ViewVC Help
Powered by ViewVC 1.1.26