/[escript]/branches/diaplayground/ripley/src/RipleySystemMatrix.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/RipleySystemMatrix.cpp

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

revision 5135 by caltinay, Mon Sep 8 07:17:34 2014 UTC revision 5136 by caltinay, Tue Sep 9 07:13:55 2014 UTC
# Line 30  Line 30 
30    
31  #define BLOCK_SIZE 128  #define BLOCK_SIZE 128
32    
33    namespace bp = boost::python;
34    
35  namespace ripley {  namespace ripley {
36    
37  double gettime()  double gettime()
# Line 94  SystemMatrix::SystemMatrix(esysUtils::JM Line 96  SystemMatrix::SystemMatrix(esysUtils::JM
96      for (size_t i = 0; i < diagonalOffsets.size(); i++) {      for (size_t i = 0; i < diagonalOffsets.size(); i++) {
97          numEntries += blocksize*blocksize*(nRows-std::abs(diagonalOffsets[i]));          numEntries += blocksize*blocksize*(nRows-std::abs(diagonalOffsets[i]));
98      }      }
99      std::cerr << "Matrix has " <<numEntries<<" entries."<<std::endl;      //std::cerr << "Matrix has " <<numEntries<<" entries."<<std::endl;
100    
101      if (cudaDevices.size() == 0)      if (cudaDevices.size() == 0)
102          checkCUDA();          checkCUDA();
# Line 149  void SystemMatrix::ypAx(escript::Data& y Line 151  void SystemMatrix::ypAx(escript::Data& y
151      HostVectorType yy(mat.num_rows, 0.);      HostVectorType yy(mat.num_rows, 0.);
152      cusp::multiply(mat, xx, yy);      cusp::multiply(mat, xx, yy);
153      thrust::copy(yy.begin(), yy.end(), y_dp);      thrust::copy(yy.begin(), yy.end(), y_dp);
154      std::cerr << "ypAx: " << gettime()-T0 << " seconds." << std::endl;      //std::cerr << "ypAx: " << gettime()-T0 << " seconds." << std::endl;
155  }  }
156    
157  template<class LinearOperator,  template<class LinearOperator,
# Line 161  void SystemMatrix::runSolver(LinearOpera Line 163  void SystemMatrix::runSolver(LinearOpera
163      typedef typename LinearOperator::memory_space MemorySpace;      typedef typename LinearOperator::memory_space MemorySpace;
164      //cusp::verbose_monitor<double> monitor(b, sb.getIterMax(), sb.getTolerance(), sb.getAbsoluteTolerance());      //cusp::verbose_monitor<double> monitor(b, sb.getIterMax(), sb.getTolerance(), sb.getAbsoluteTolerance());
165      cusp::default_monitor<double> monitor(b, sb.getIterMax(), sb.getTolerance(), sb.getAbsoluteTolerance());      cusp::default_monitor<double> monitor(b, sb.getIterMax(), sb.getTolerance(), sb.getAbsoluteTolerance());
166        int solver = sb.getSolverMethod();
167        if (solver == escript::SO_DEFAULT) {
168            if (sb.isSymmetric()) {
169                solver = escript::SO_METHOD_PCG;
170            } else {
171                solver = escript::SO_METHOD_BICGSTAB;
172            }
173        }
174    
175      double T0 = gettime();      double T0 = gettime();
176      switch (sb.getSolverMethod()) {      switch (solver) {
177          case escript::SO_DEFAULT:          case escript::SO_DEFAULT:
178          case escript::SO_METHOD_PCG:          case escript::SO_METHOD_PCG:
179              cusp::krylov::cg(A, x, b, monitor, M);              cusp::krylov::cg(A, x, b, monitor, M);
# Line 176  void SystemMatrix::runSolver(LinearOpera Line 186  void SystemMatrix::runSolver(LinearOpera
186              break;              break;
187          case escript::SO_METHOD_GMRES:          case escript::SO_METHOD_GMRES:
188              {              {
189                  const int restart = sb._getRestartForC();                  const int restart = (sb.getRestart()==0 ? 1000 : sb.getRestart());
190                  if (restart < 1)                  if (restart < 1)
191                      throw RipleyException("Invalid restart parameter for GMRES");                      throw RipleyException("Invalid restart parameter for GMRES");
192                  cusp::krylov::gmres(A, x, b, restart, monitor, M);                  cusp::krylov::gmres(A, x, b, restart, monitor, M);
# Line 208  void SystemMatrix::runSolver(LinearOpera Line 218  void SystemMatrix::runSolver(LinearOpera
218  }  }
219    
220  void SystemMatrix::setToSolution(escript::Data& out, escript::Data& in,  void SystemMatrix::setToSolution(escript::Data& out, escript::Data& in,
221                                   boost::python::object& options) const                                   bp::object& options) const
222  {  {
223      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
224          throw RipleyException("solve: ripley's block diagonal matrix "          throw RipleyException("solve: ripley's block diagonal matrix "
# Line 225  void SystemMatrix::setToSolution(escript Line 235  void SystemMatrix::setToSolution(escript
235      }      }
236    
237      options.attr("resetDiagnostics")();      options.attr("resetDiagnostics")();
238      escript::SolverBuddy sb = boost::python::extract<escript::SolverBuddy>(options);      escript::SolverBuddy sb = bp::extract<escript::SolverBuddy>(options);
239      out.expand();      out.expand();
240      in.expand();      in.expand();
241    
# Line 351  void SystemMatrix::nullifyRowsAndCols(es Line 361  void SystemMatrix::nullifyRowsAndCols(es
361              }              }
362          }          }
363      }      }
364      std::cerr << "nullifyRowsAndCols: " << gettime()-T0 << " seconds." << std::endl;      //std::cerr << "nullifyRowsAndCols: " << gettime()-T0 << " seconds." << std::endl;
365      matrixAltered = true;      matrixAltered = true;
366  }  }
367    

Legend:
Removed from v.5135  
changed lines
  Added in v.5136

  ViewVC Help
Powered by ViewVC 1.1.26