/[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 5064 by caltinay, Mon May 19 04:47:10 2014 UTC revision 5065 by caltinay, Fri Jun 20 05:49:52 2014 UTC
# Line 19  Line 19 
19    
20  #include <esysUtils/index.h>  #include <esysUtils/index.h>
21  #include <escript/Data.h>  #include <escript/Data.h>
22    #include <escript/SolverOptions.h>
23    
24    #include <cusp/multiply.h>
25    #include <cusp/krylov/bicgstab.h>
26    #include <cusp/krylov/cg.h>
27    #include <cusp/krylov/gmres.h>
28    
29  #define BLOCK_SIZE 128  #define BLOCK_SIZE 128
30    
31  namespace ripley {  namespace ripley {
32    
33    double gettime()
34    {
35    #ifdef _OPENMP
36        return omp_get_wtime();
37    #else
38        struct timeval tv;
39        gettimeofday(&tv, NULL);
40        suseconds_t ret = tv.tv_usec + tv.tv_sec*1e6;
41        return 1e-6*(double)ret;
42    #endif
43    }
44    
45    void list_devices(void)
46    {
47    #ifdef USE_CUDA
48        int deviceCount;
49        CUDA_SAFE_CALL(cudaGetDeviceCount(&deviceCount));
50        if (deviceCount == 0)
51            std::cout << "There is no device supporting CUDA" << std::endl;
52    
53        for (int dev = 0; dev < deviceCount; ++dev) {
54            cudaDeviceProp deviceProp;
55            CUDA_SAFE_CALL(cudaGetDeviceProperties(&deviceProp, dev));
56    
57            if (dev == 0) {
58                if (deviceProp.major == 9999 && deviceProp.minor == 9999)
59                    std::cout << "There is no device supporting CUDA." << std::endl;
60                else if (deviceCount == 1)
61                    std::cout << "There is 1 device supporting CUDA" << std:: endl;
62                else
63                    std::cout << "There are " << deviceCount <<  " devices supporting CUDA" << std:: endl;
64            }
65    
66            std::cout << "\nDevice " << dev << ": \"" << deviceProp.name << "\"" << std::endl;
67            std::cout << "  Major revision number:                         " << deviceProp.major << std::endl;
68            std::cout << "  Minor revision number:                         " << deviceProp.minor << std::endl;
69            std::cout << "  Total amount of global memory:                 " << deviceProp.totalGlobalMem << " bytes" << std::endl;
70        }
71        std::cout << std::endl;
72    #endif
73    }
74    
75  SystemMatrix::SystemMatrix()  SystemMatrix::SystemMatrix()
76  {  {
77      throw RipleyException("Illegal to instantiate SystemMatrix without arguments.");      throw RipleyException("Illegal to instantiate SystemMatrix without arguments.");
# Line 31  SystemMatrix::SystemMatrix() Line 79  SystemMatrix::SystemMatrix()
79    
80  SystemMatrix::SystemMatrix(int blocksize, const escript::FunctionSpace& fs,  SystemMatrix::SystemMatrix(int blocksize, const escript::FunctionSpace& fs,
81                             int nRows, const IndexVector& diagonalOffsets) :                             int nRows, const IndexVector& diagonalOffsets) :
82      AbstractSystemMatrix(blocksize, fs, blocksize, fs),      AbstractSystemMatrix(blocksize, fs, blocksize, fs)
     numRows(nRows*blocksize),  
     offsets(diagonalOffsets)  
83  {  {
84      values.resize(numRows*offsets.size()*blocksize, 0.);      // count nonzero entries
85        int numEntries = 0;
86        for (size_t i = 0; i < diagonalOffsets.size(); i++) {
87            numEntries += blocksize*blocksize *
88                (nRows-std::abs(diagonalOffsets[i])*blocksize) / blocksize;
89        }
90    
91        list_devices();
92    #ifdef USE_CUDA
93        //TODO: give users options...
94        cudaSetDevice(0);
95    #endif
96    
97        mat.resize(nRows*blocksize, nRows*blocksize, numEntries, diagonalOffsets.size()*blocksize);
98        mat.diagonal_offsets.assign(diagonalOffsets.begin(), diagonalOffsets.end());
99  }  }
100    
101  void SystemMatrix::add(const IndexVector& rowIdx,  void SystemMatrix::add(const IndexVector& rowIdx,
# Line 43  void SystemMatrix::add(const IndexVector Line 103  void SystemMatrix::add(const IndexVector
103  {  {
104      const int blockSize = getBlockSize();      const int blockSize = getBlockSize();
105      const int emSize = rowIdx.size();      const int emSize = rowIdx.size();
 #if 0  
     static bool here=false;  
     if (here) return;  
     here=true;  
     double* arr = const_cast<double*>(&array[0]);  
   
     for (int i=0; i<4*blockSize; i++) {  
         for (int j=0; j<4*blockSize; j++) {  
             arr[i*4*blockSize+j] = i*4*blockSize+j;  
         }  
     }  
 #endif  
   
     //for (k in numEq) for (m in numComp)  
     //array[k + numEq*m + numEq*numComp*i + numEq*numComp*emSize*j]  
     //std::cerr << "SystemMatrix::add" << std::endl;  
106      for (int i=0; i<emSize; i++) {      for (int i=0; i<emSize; i++) {
107          for (int j=0; j<emSize; j++) {          for (int j=0; j<emSize; j++) {
108              const int revi = emSize-1-i;              const int revi = emSize-1-i;
109              const int diag = j%2 + revi%2 + 3*(j/2+revi/2 + j/4+revi/4);              const int diag = j%2 + revi%2 + 3*(j/2+revi/2 + j/4+revi/4);
110              for (int k=0; k<blockSize; k++) {              for (int k=0; k<blockSize; k++) {
111                  for (int m=0; m<blockSize; m++) {                  for (int m=0; m<blockSize; m++) {
112                      const int destIdx = rowIdx[i]*blockSize + k + numRows*(diag*blockSize+m);                      const int row = rowIdx[i]*blockSize + k;
113                        const int d = diag*blockSize + m;
114                      const int srcIdx = INDEX4(k, m, i, j, blockSize, blockSize, emSize);                      const int srcIdx = INDEX4(k, m, i, j, blockSize, blockSize, emSize);
115                      values[destIdx] += array[srcIdx];                      mat.values(row, d) += array[srcIdx];
116                  }                  }
117              }              }
118          }          }
# Line 93  void SystemMatrix::ypAx(escript::Data& y Line 138  void SystemMatrix::ypAx(escript::Data& y
138      y.requireWrite();      y.requireWrite();
139      const double* x_dp = x.getSampleDataRO(0);      const double* x_dp = x.getSampleDataRO(0);
140      double* y_dp = y.getSampleDataRW(0);      double* y_dp = y.getSampleDataRW(0);
141      matrixVector(x_dp, 1., y_dp);      double T0 = gettime();
142        HostVectorType xx(x_dp, x_dp+mat.num_rows);
143        HostVectorType yy(mat.num_rows, 0.);
144        cusp::multiply(mat, xx, yy);
145        thrust::copy(yy.begin(), yy.end(), y_dp);
146        std::cerr << "ypAx: " << gettime()-T0 << " seconds." << std::endl;
147  }  }
148    
149  // y = beta y + Ax  // y = beta y + Ax
150  void SystemMatrix::matrixVector(const double* x, double beta, double* y) const  void SystemMatrix::matrixVector(const double* x, double beta, double* y) const
151  {  {
152        //cusp::multiply(mat, x, y);
153    
154      const int blockSize = getBlockSize();      const int blockSize = getBlockSize();
155      if (blockSize == 1) {      if (blockSize == 1) {
156  #pragma omp parallel for  #pragma omp parallel for
157          for (int ch=0; ch<numRows; ch+=BLOCK_SIZE) {          for (size_t ch=0; ch<mat.num_rows; ch+=BLOCK_SIZE) {
158              // initialize chunk              // initialize chunk
159              if (std::abs(beta) > 0.) {              if (std::abs(beta) > 0.) {
160                  if (beta != 1.) {                  if (beta != 1.) {
161                      for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                      for (int row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
162                          y[row] *= beta;                          y[row] *= beta;
163                      }                      }
164                  }                  }
165              } else {              } else {
166                  for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                  for (size_t row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
167                      y[row] = 0.;                      y[row] = 0.;
168                  }                  }
169              }              }
170    
171              // for each diagonal              // for each diagonal
172              for (size_t d=0; d<offsets.size(); d++) {              for (size_t d=0; d<mat.diagonal_offsets.size(); d++) {
173                  for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                  for (int row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
174                      const int col = row + offsets[d];                      const int col = row + mat.diagonal_offsets[d];
175                      if (col >= 0 && col < numRows) {                      if (col >= 0 && col < mat.num_rows) {
176                          y[row] += values[row+d*numRows] * x[col];                          y[row] += mat.values(row, d) * x[col];
177                      }                      }
178                  }                  }
179              }              }
180          }          }
181      } else {      } else {
182  #pragma omp parallel for  #pragma omp parallel for
183          for (int ch=0; ch<numRows; ch+=BLOCK_SIZE) {          for (size_t ch=0; ch < mat.num_rows; ch+=BLOCK_SIZE) {
184              // initialize chunk              // initialize chunk
185              if (std::abs(beta) > 0.) {              if (std::abs(beta) > 0.) {
186                  if (beta != 1.) {                  if (beta != 1.) {
187                      for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                      for (int row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
188                          y[row] *= beta;                          y[row] *= beta;
189                      }                      }
190                  }                  }
191              } else {              } else {
192                  for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                  for (int row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
193                      y[row] = 0.;                      y[row] = 0.;
194                  }                  }
195              }              }
196    
197              // for each diagonal block              // for each diagonal block
198              for (size_t d=0; d<offsets.size(); d++) {              for (size_t d=0; d < mat.diagonal_offsets.size(); d++) {
199                  const int k = offsets[d]*blockSize;                  const int k = mat.diagonal_offsets[d]*blockSize;
200                  for (int row=ch; row<std::min(ch+BLOCK_SIZE,numRows); row++) {                  for (int row=ch; row<std::min(ch+BLOCK_SIZE, mat.num_rows); row++) {
201                      const int col = blockSize*(row/blockSize) + k;                      const int col = blockSize*(row/blockSize) + k;
202                      if (col >= 0 && col <= numRows-blockSize) {                      if (col >= 0 && col <= mat.num_rows-blockSize) {
203                          // for each column in block                          // for each column in block
204                          for (int i=0; i<blockSize; i++) {                          for (int i=0; i<blockSize; i++) {
205                              const double Aij = values[row+(d*blockSize+i)*numRows];                              const double Aij = mat.values(row, d*blockSize+i);
206                              y[row] += Aij * x[col+i];                              y[row] += Aij * x[col+i];
207                          }                          }
208                      }                      }
# Line 160  void SystemMatrix::matrixVector(const do Line 212  void SystemMatrix::matrixVector(const do
212      }      }
213  }  }
214    
215  inline void axpby(size_t N, const double* x, const double* y, double* z,  template<class LinearOperator,
216                    double alpha, double beta)           class Vector>
217  {  void SystemMatrix::runSolver(LinearOperator& A, Vector& x, Vector& b,
218  #pragma omp parallel for                               escript::SolverBuddy& sb) const
219      for (size_t i=0; i<N; i++) {  {
220          z[i] = alpha*x[i] + beta*y[i];      typedef typename LinearOperator::memory_space MemorySpace;
221      }      //sb.isVerbose()
222  }      //sb.getPreconditioner()
223        cusp::default_monitor<double> monitor(b, sb.getIterMax(), sb.getTolerance(), sb.getAbsoluteTolerance());
224  inline void axpy(size_t N, const double* y, double* x, double alpha)      cusp::identity_operator<double, MemorySpace> M(mat.num_rows, mat.num_rows);
225  {  
226  #pragma omp parallel for      double T0 = gettime();
227      for (size_t i=0; i<N; i++) {      switch (sb.getSolverMethod()) {
228          x[i] += alpha*y[i];          case escript::ESCRIPT_DEFAULT:
229      }          case escript::ESCRIPT_PCG:
230  }              cusp::krylov::cg(A, x, b, monitor, M);
231                break;
232  inline double dotc(const std::vector<double>& x, const std::vector<double>& y)          case escript::ESCRIPT_BICGSTAB:
233  {              cusp::krylov::bicgstab(A, x, b, monitor, M);
234      double ret = 0.;              break;
235      for (size_t i=0; i<x.size(); i++)          case escript::ESCRIPT_GMRES:
236          ret += x[i]*y[i];              {
237      return ret;                  const int restart = sb._getRestartForC();
238  }                  if (restart < 1)
239                        throw RipleyException("Invalid restart parameter for GMRES");
240  inline double nrm2(size_t N, const double* x)                  cusp::krylov::gmres(A, x, b, restart, monitor, M);
241  {              }
242      double ret = 0.;              break;
243      for (size_t i=0; i<N; i++)          case escript::ESCRIPT_PRES20:
244          ret += x[i]*x[i];              {
245      return std::sqrt(ret);                  const int restart = 20;
246  }                  cusp::krylov::gmres(A, x, b, restart, monitor, M);
247                }
248  void SystemMatrix::cg(double* x, const double* b) const              break;
249  {          default:
250      const size_t N = numRows;              throw RipleyException("Unsupported solver.");
251        }
252      // allocate workspace      double solvertime = gettime()-T0;
253      std::vector<double> y(N);  
254      std::vector<double> z(N);      if (monitor.converged()) {
255      std::vector<double> r(N);          std::cout << "Solver converged to " << monitor.relative_tolerance()
256      std::vector<double> p(N);              << " relative tolerance after " << monitor.iteration_count()
257                << " iterations and " << solvertime << " seconds."<< std::endl;
258  #pragma omp parallel for      } else {
259      for (size_t i=0; i<N; i++)          std::cout << "Solver reached iteration limit "
260          x[i] = b[i];              << monitor.iteration_limit() << " before converging"
261                << " to " << monitor.relative_tolerance() << " rel. tolerance."
262      // y <- Ax              << std::endl;
     matrixVector(x, 0., &y[0]);  
   
     // r <- b - A*x  
     axpby(N, b, &y[0], &r[0], 1., -1.);  
     
     // z <- M*r  
     //cusp::multiply(M, r, z);  
     z=r;  
   
     // p <- z  
     //blas::copy(z, p);  
 #pragma omp parallel for  
     for (size_t i=0; i<N; i++)  
         p[i] = z[i];  
   
     // rz = <r^H, z>  
     double rz = dotc(r, z);  
   
     const double b_norm = nrm2(N, b);  
   
     int iteration=0;  
   
     while (nrm2(N, &r[0]) > 1e-9*b_norm && iteration<10000)  
     {  
         //std::cerr << nrm2(N, &r[0])<< " " << nrm2(N, x)<<std::endl;  
         iteration++;  
         // y <- Ap  
         matrixVector(&p[0], 0., &y[0]);  
   
         // alpha <- <r,z>/<y,p>  
         const double alpha =  rz / dotc(y, p);  
   
         // x <- x + alpha * p  
         axpy(N, &p[0], x, alpha);  
   
         // r <- r - alpha * y        
         axpy(N, &y[0], &r[0], -alpha);  
   
         // z <- M*r  
         //cusp::multiply(M, r, z);  
 #pragma omp parallel for  
         for (size_t i=0; i<N; i++)  
             z[i] = r[i];  
           
         double rz_old = rz;  
   
         // rz = <r^H, z>  
         rz = dotc(r, z);  
   
         // beta <- <r_{i+1},r_{i+1}>/<r,r>  
         double beta = rz / rz_old;  
           
         // p <- r + beta*p  
         axpby(N, &z[0], &p[0], &p[0], 1., beta);  
263      }      }
264  }  }
265    
# Line 279  void SystemMatrix::setToSolution(escript Line 277  void SystemMatrix::setToSolution(escript
277      }      }
278    
279      std::cerr << "SystemMatrix::setToSolution" << std::endl;      std::cerr << "SystemMatrix::setToSolution" << std::endl;
280        options.attr("resetDiagnostics")();
281        escript::SolverBuddy sb = boost::python::extract<escript::SolverBuddy>(options);
282      out.expand();      out.expand();
283      in.expand();      in.expand();
284    
285      double* out_dp = out.getSampleDataRW(0);      double* out_dp = out.getSampleDataRW(0);
286      const double* in_dp = in.getSampleDataRO(0);      const double* in_dp = in.getSampleDataRO(0);
287      //solve(out_dp, in_dp, &paso_options);      double T0;
288      cg(out_dp, in_dp);  
289      //std::cerr << out.toString() << std::endl;      if (sb.getSolverTarget() == escript::ESCRIPT_TARGET_GPU) {
290    #ifdef USE_CUDA
291            T0 = gettime();
292            DeviceVectorType b(in_dp, in_dp+mat.num_rows);
293            double host2dev = gettime()-T0;
294            std::cerr << "Copy of b: " << host2dev << " seconds." << std::endl;
295            T0 = gettime();
296            DeviceMatrixType dmat = mat;
297            host2dev = gettime()-T0;
298            std::cerr << "Copy of A: " << host2dev << " seconds." << std::endl;
299            DeviceVectorType x(mat.num_rows, 0.);
300            std::cerr << "Solving on CUDA device" << std::endl;
301            runSolver(dmat, x, b, sb);
302            T0 = gettime();
303            thrust::copy(x.begin(), x.end(), out_dp);
304            const double copyTime = gettime()-T0;
305            std::cerr << "Copy of x: " << copyTime << " seconds." << std::endl;
306    #else
307            throw RipleyException("solve: GPU-based solver requested but escript "
308                                  "not compiled with CUDA.");
309    #endif
310        } else { // CPU
311            T0 = gettime();
312            HostVectorType b(in_dp, in_dp+mat.num_rows);
313            double copytime = gettime()-T0;
314            std::cerr << "Copy of b: " << copytime << " seconds." << std::endl;
315            HostVectorType x(mat.num_rows, 0.);
316            std::cerr << "Solving on host" << std::endl;
317            runSolver(mat, x, b, sb);
318            T0 = gettime();
319            thrust::copy(x.begin(), x.end(), out_dp);
320            const double copyTime = gettime()-T0;
321            std::cerr << "Copy of x: " << copyTime << " seconds." << std::endl;
322        }
323    
324  }  }
325    
326  void SystemMatrix::nullifyRowsAndCols(escript::Data& row_q,  void SystemMatrix::nullifyRowsAndCols(escript::Data& row_q,
327                                        escript::Data& col_q,                                        escript::Data& col_q,
328                                        double mdv)                                        double mdv)
329  {  {
330        double T0 = gettime();
331      if (col_q.getDataPointSize() != getColumnBlockSize()) {      if (col_q.getDataPointSize() != getColumnBlockSize()) {
332          throw RipleyException("nullifyRowsAndCols: column block size does not match the number of components of column mask.");          throw RipleyException("nullifyRowsAndCols: column block size does not match the number of components of column mask.");
333      } else if (row_q.getDataPointSize() != getRowBlockSize()) {      } else if (row_q.getDataPointSize() != getRowBlockSize()) {
# Line 309  void SystemMatrix::nullifyRowsAndCols(es Line 345  void SystemMatrix::nullifyRowsAndCols(es
345      const double* colMask = col_q.getSampleDataRO(0);      const double* colMask = col_q.getSampleDataRO(0);
346      const int blockSize = getBlockSize();      const int blockSize = getBlockSize();
347  #pragma omp parallel for  #pragma omp parallel for
348      for (int row=0; row < numRows; row++) {      for (int row=0; row < mat.num_rows; row++) {
349          for (int diag=0; diag < offsets.size(); diag++) {          for (int diag=0; diag < mat.diagonal_offsets.size(); diag++) {
350              const int col = blockSize*(row/blockSize)+offsets[diag]*blockSize;              const int col = blockSize*(row/blockSize)+mat.diagonal_offsets[diag]*blockSize;
351              if (col >= 0 && col <= numRows-blockSize) {              if (col >= 0 && col <= mat.num_rows-blockSize) {
352                  for (int i=0; i<blockSize; i++) {                  for (int i=0; i<blockSize; i++) {
353                      if (rowMask[row] > 0. || colMask[col+i] > 0.) {                      if (rowMask[row] > 0. || colMask[col+i] > 0.) {
354                          values[row+(diag*blockSize+i)*numRows] =                          mat.values(row, diag*blockSize+i) =
355                                                          (row==col+i ? mdv : 0);                                                          (row==col+i ? mdv : 0);
356                      }                      }
357                  }                  }
358              }              }
359          }          }
360      }      }
361        std::cerr << "nullifyRowsAndCols: " << gettime()-T0 << " seconds." << std::endl;
362  }  }
363    
364  void SystemMatrix::saveMM(const std::string& filename) const  void SystemMatrix::saveMM(const std::string& filename) const
365  {  {
366      const int blockSize = getBlockSize();      const int blockSize = getBlockSize();
367    
     // count nonzero entries  
     int numNZ = 0;  
     for (size_t i = 0; i < offsets.size(); i++) {  
         numNZ += blockSize*blockSize*(numRows-std::abs(offsets[i])*blockSize) / blockSize;  
     }  
       
368      std::ofstream f(filename.c_str());      std::ofstream f(filename.c_str());
369      f << "%%%%MatrixMarket matrix coordinate real general" << std::endl;      f << "%%%%MatrixMarket matrix coordinate real general" << std::endl;
370      f << numRows << " " << numRows << " " << numNZ << std::endl;      f << mat.num_rows << " " << mat.num_rows << " " << mat.num_entries << std::endl;
371      for (int row=0; row<numRows; row++) {      for (int row=0; row < mat.num_rows; row++) {
372          for (int diag=0; diag<offsets.size(); diag++) {          for (int diag=0; diag < mat.diagonal_offsets.size(); diag++) {
373              const int col = blockSize*(row/blockSize)+offsets[diag]*blockSize;              const int col = blockSize*(row/blockSize)+mat.diagonal_offsets[diag]*blockSize;
374              if (col >= 0 && col <= numRows-blockSize) {              if (col >= 0 && col <= mat.num_rows-blockSize) {
375                  for (int i=0; i<blockSize; i++) {                  for (int i=0; i<blockSize; i++) {
376                      f << row+1 << " " << col+i+1 << " "                      f << row+1 << " " << col+i+1 << " "
377                            << values[(diag*blockSize+i)*numRows+row] << std::endl;                            << mat.values(row, diag*blockSize+i) << std::endl;
378                  }                  }
379              }              }
380          }          }
# Line 357  void SystemMatrix::saveHB(const std::str Line 388  void SystemMatrix::saveHB(const std::str
388    
389  void SystemMatrix::resetValues()  void SystemMatrix::resetValues()
390  {  {
391      values.assign(values.size(), 0.);      mat.values.values.assign(mat.values.values.size(), 0.);
392  }  }
393    
394  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.5064  
changed lines
  Added in v.5065

  ViewVC Help
Powered by ViewVC 1.1.26