text: #ifndef CRS_MATRIX_HPP
  #define CRS_MATRIX_HPP
  
  #include <fstream>
  #include <dirent.h>
  #include <sys/stat.h>
  #include <sys/types.h>
  #include <climits>
  #include <vector>
  
  //#include "../core/dvector.hpp"
  //#include "../core/dunordered_map.hpp"
  //#include "../core/mpihelper.hpp"
  //#include "../core/prefix_sum.hpp"
  
  //#include "rowmajor_matrix.hpp"
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <sys/stat.h>
  
  #undef FTRACE_A
  #undef FTRACE_B
  #if defined(_SX)||defined(__ve) && 0
  #include <ftrace.h>
  #define FTRACE_A(name) ftrace_region_begin(name)
  #define FTRACE_B(name) ftrace_region_end(name)
  #else
  #define FTRACE_A(name)
  #define FTRACE_B(name)
  #endif
  
  
  
  #define CRS_SPMM_THR 32
  #define CRS_SPMM_VLEN 256
  #define SPARSE_VECTOR_VLEN 256
  #define TO_SKIP_REBALANCE 256
  
  //This file was copied from Frovedis.
  
  
  template <class T, class I = size_t, class O = size_t>
  struct crs_matrix_local {
      crs_matrix_local() {}
      crs_matrix_local(size_t nrows, size_t ncols) :
          local_num_row(nrows), local_num_col(ncols) {}
      void view_data(){
          printf("crs data\n");
          printf("----------------------------\n");
          printf("col %d, row %d\n",local_num_col,local_num_row);
          printf("first off %d, last off %d\n",off[0],off[local_num_row]);
          printf("----------------------------\n");
      }
  
      T* val; // size is off[local_num_row]
      I* idx; // size is off[local_num_row]
      O* off; // size is local_num_row + 1 ("0" is added)
      T* buf; // size is local_num_row
      size_t local_num_row;
      size_t local_num_col;
  
  
  };
  
  
  
  
  template <class T, class I, class O>
  void crs_matrix_spmv_impl(const crs_matrix_local<T,I,O>& mat, T* retp, const T* vp) {
      const T* valp = &mat.val[0];
      const I* idxp = &mat.idx[0];
      const O* offp = &mat.off[0];
  
      //    printf("crs_matrix_spmv_impl\n");
      //    printf("c %d, r %d\n",mat.local_num_col,mat.local_num_row);
      //    printf("val %d, idx %d, off %d\n",mat.val.size(),mat.idx.size(),mat.off.size());
      //    printf("off: s %d, e %d\n",mat.off[0],mat.off[mat.local_num_row]);
  
      FTRACE_A("crs impl");
  
  #if 0
      size_t size = mat.val.size();
      std::vector<T> crspart(size);
      T* crspartp = crspart.data();
  
  
  
  
      FTRACE_A("crs impl mul");
  
      size_t size8 = (size>>3)<<3;
  
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
      for(size_t i = 0; i < size8; i+=8) {
          size_t pos0 = idxp[i+0];
          size_t pos1 = idxp[i+1];
          size_t pos2 = idxp[i+2];
          size_t pos3 = idxp[i+3];
          size_t pos4 = idxp[i+4];
          size_t pos5 = idxp[i+5];
          size_t pos6 = idxp[i+6];
          size_t pos7 = idxp[i+7];
  
  
          crspartp[i+0] = valp[i+0] * vp[pos0];
          crspartp[i+1] = valp[i+1] * vp[pos1];
          crspartp[i+2] = valp[i+2] * vp[pos2];
          crspartp[i+3] = valp[i+3] * vp[pos3];
          crspartp[i+4] = valp[i+4] * vp[pos4];
          crspartp[i+5] = valp[i+5] * vp[pos5];
          crspartp[i+6] = valp[i+6] * vp[pos6];
          crspartp[i+7] = valp[i+7] * vp[pos7];
      }
  
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
      for(size_t i = size8; i < size; i++) {
          crspartp[i] = valp[i] * vp[idxp[i]];
      }
  
  
      //#pragma _NEC vovertake
      ////#pragma cdir on_adb(vp)
      //    for(size_t i = 0; i < size; i++) {
      //        crspartp[i] = valp[i] * vp[idxp[i]];
      //    }
      FTRACE_B("crs impl mul");
  
      FTRACE_A("crs impl sum");
      //    for(size_t r = 0; r < mat.local_num_row; r++) {
      //        if(offp[r+1]-offp[r]>0){
      //#pragma _NEC vovertake
      ////#pragma cdir on_adb(vp)
      //            for(O c = offp[r]; c < offp[r+1]; c++) {
      //                retp[r] = retp[r] + crspartp[c];
      //            }
      //        }
      //        else if(offp[r+1]-offp[r]<0) printf("err\n");
      //    }
      for(size_t r = 0; r < mat.local_num_row; r++) {
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
          for(O c = offp[r]; c < offp[r+1]; c++) {
              retp[r] = retp[r] + crspartp[c];
          }
      }
      FTRACE_B("crs impl sum");
  
  #elif 0
  
  #define MIN(a,b)        ((a)<(b)?(a):(b))
  
      int row8 = (mat.local_num_row>>3)<<3;
  
      FTRACE_A("crs 8");
      for(size_t r = 0; r < row8; r+=8) {
          int len0 = offp[r+1] - offp[r];
          int len1 = offp[r+2] - offp[r+1];
          int len2 = offp[r+3] - offp[r+2];
          int len3 = offp[r+4] - offp[r+3];
          int len4 = offp[r+5] - offp[r+4];
          int len5 = offp[r+6] - offp[r+5];
          int len6 = offp[r+7] - offp[r+6];
          int len7 = offp[r+8] - offp[r+7];
          int len = MIN(MIN(MIN(len0,len1),MIN(len2,len3)),MIN(MIN(len4,len5),MIN(len6,len7)));
          int off0 = offp[r];
          int off1 = offp[r+1] - offp[r];
          int off2 = offp[r+2] - offp[r];
          int off3 = offp[r+3] - offp[r];
          int off4 = offp[r+4] - offp[r];
          int off5 = offp[r+5] - offp[r];
          int off6 = offp[r+6] - offp[r];
          int off7 = offp[r+7] - offp[r];
          //        printf("r %d, len %d\n",r,len);
  
          if(r>15) FTRACE_A("crs 8 in A");
          else FTRACE_A("crs 8 in B");
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
          for(O c = 0; c < len; c++) {
              const I pos0 = idxp[c+off0];
              const I pos1 = idxp[c+off1];
              const I pos2 = idxp[c+off2];
              const I pos3 = idxp[c+off3];
              const I pos4 = idxp[c+off4];
              const I pos5 = idxp[c+off5];
              const I pos6 = idxp[c+off6];
              const I pos7 = idxp[c+off7];
              retp[r+0] = retp[r+0] + valp[c+off0] * vp[pos0];
              retp[r+1] = retp[r+1] + valp[c+off1] * vp[pos1];
              retp[r+2] = retp[r+2] + valp[c+off2] * vp[pos2];
              retp[r+3] = retp[r+3] + valp[c+off3] * vp[pos3];
              retp[r+4] = retp[r+4] + valp[c+off4] * vp[pos4];
              retp[r+5] = retp[r+5] + valp[c+off5] * vp[pos5];
              retp[r+6] = retp[r+6] + valp[c+off6] * vp[pos6];
              retp[r+7] = retp[r+7] + valp[c+off7] * vp[pos7];
          }
          if(r>15) FTRACE_B("crs 8 in A");
          else FTRACE_B("crs 8 in B");
  
          if(r>15) FTRACE_A("crs in A");
          else FTRACE_A("crs in B");
          for(size_t rr = 0;rr<8;rr++){
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
              for(O c = offp[r+rr]+len; c < offp[r+1+rr]; c++) {
                  retp[r+rr] = retp[r+rr] + valp[c] * vp[idxp[c]];
              }
          }
          if(r>15) FTRACE_B("crs in A");
          else FTRACE_B("crs in B");
      }
      FTRACE_B("crs 8");
      for(size_t r = row8; r < mat.local_num_row; r++) {
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
          for(O c = offp[r]; c < offp[r+1]; c++) {
              retp[r] = retp[r] + valp[c] * vp[idxp[c]];
          }
      }
  
  #else
  
      //    for(size_t r = 0; r < mat.local_num_row; r++) {
      //        if(offp[r+1]-offp[r]>0){
      //#pragma _NEC vovertake
      ////#pragma cdir on_adb(vp)
      //            for(O c = offp[r]; c < offp[r+1]; c++) {
      //                retp[r] += valp[c] * vp[idxp[c]];
      //            }
      //        }
      //        else if(offp[r+1]-offp[r]<0) printf("err\n");
      //    }
  
  
      for(size_t r = 0; r < mat.local_num_row; r++) {
  #pragma _NEC vovertake
  //#pragma cdir on_adb(vp)
          for(O c = offp[r]; c < offp[r+1]; c++) {
              retp[r] += valp[c] * vp[idxp[c]];
          }
      }
  #endif
      FTRACE_B("crs impl");
      //    printf("crs_matrix_spmv_impl end\n");
  }
  
  
  template <class T, class I, class O>
  std::vector<T> operator*(const crs_matrix_local<T,I,O>& mat, const std::vector<T>& v) {
      std::vector<T> ret(mat.local_num_row);
      if(mat.local_num_col != v.size())
          throw std::runtime_error("operator*: size of vector does not match");
      //    printf("c %d, r %d, ret %d, v %d\n",mat.local_num_col,mat.local_num_row,ret.size(),v.size());
      crs_matrix_spmv_impl(mat, ret.data(), v.data());
      return ret;
  }
  
  
  
  
  
  template <class T, class I, class O>
  void crs_matrix_spmm_impl2(const crs_matrix_local<T,I,O>& mat,
                             T* retvalp, const T* vvalp, size_t num_col) {
      const T* valp = &mat.val[0];
      const I* idxp = &mat.idx[0];
      const O* offp = &mat.off[0];
      T current_sum[CRS_SPMM_VLEN];
  #pragma _NEC vreg(current_sum)
      for(size_t i = 0; i < CRS_SPMM_VLEN; i++) {
          current_sum[i] = 0;
      }
      size_t each = num_col / CRS_SPMM_VLEN;
      size_t rest = num_col % CRS_SPMM_VLEN;
      for(size_t r = 0; r < mat.local_num_row; r++) {
          for(size_t e = 0; e < each; e++) {
              for(O c = offp[r]; c < offp[r+1]; c++) {
                  for(size_t mc = 0; mc < CRS_SPMM_VLEN; mc++) {
                      current_sum[mc] +=
                              valp[c] * vvalp[idxp[c] * num_col + CRS_SPMM_VLEN * e + mc];
                  }
              }
              for(size_t mc = 0; mc < CRS_SPMM_VLEN; mc++) {
                  retvalp[r * num_col + CRS_SPMM_VLEN * e + mc] += current_sum[mc];
              }
              for(size_t i = 0; i < CRS_SPMM_VLEN; i++) {
                  current_sum[i] = 0;
              }
          }
          for(O c = offp[r]; c < offp[r+1]; c++) {
              for(size_t mc = 0; mc < rest; mc++) {
                  current_sum[mc] +=
                          valp[c] * vvalp[idxp[c] * num_col + CRS_SPMM_VLEN * each + mc];
              }
          }
          for(size_t mc = 0; mc < rest; mc++) {
              retvalp[r * num_col + CRS_SPMM_VLEN * each + mc] += current_sum[mc];
          }
          for(size_t i = 0; i < rest; i++) {
              current_sum[i] = 0;
          }
      }
  }
  
  template <class T, class I, class O>
  void crs_matrix_spmm_impl(const crs_matrix_local<T,I,O>& mat,
                            T* retvalp, const T* vvalp, size_t num_col) {
      if(num_col < CRS_SPMM_THR) {
          const T* valp = &mat.val[0];
          const I* idxp = &mat.idx[0];
          const O* offp = &mat.off[0];
          for(size_t r = 0; r < mat.local_num_row; r++) {
              size_t mc = 0;
              for(; mc + 15 < num_col; mc += 16) {
                  for(O c = offp[r]; c < offp[r+1]; c++) {
                      retvalp[r * num_col + mc] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc];
                      retvalp[r * num_col + mc + 1] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 1];
                      retvalp[r * num_col + mc + 2] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 2];
                      retvalp[r * num_col + mc + 3] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 3];
                      retvalp[r * num_col + mc + 4] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 4];
                      retvalp[r * num_col + mc + 5] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 5];
                      retvalp[r * num_col + mc + 6] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 6];
                      retvalp[r * num_col + mc + 7] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 7];
                      retvalp[r * num_col + mc + 8] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 8];
                      retvalp[r * num_col + mc + 9] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 9];
                      retvalp[r * num_col + mc + 10] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 10];
                      retvalp[r * num_col + mc + 11] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 11];
                      retvalp[r * num_col + mc + 12] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 12];
                      retvalp[r * num_col + mc + 13] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 13];
                      retvalp[r * num_col + mc + 14] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 14];
                      retvalp[r * num_col + mc + 15] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 15];
                  }
              }
              for(; mc + 7 < num_col; mc += 8) {
                  for(O c = offp[r]; c < offp[r+1]; c++) {
                      retvalp[r * num_col + mc] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc];
                      retvalp[r * num_col + mc + 1] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 1];
                      retvalp[r * num_col + mc + 2] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 2];
                      retvalp[r * num_col + mc + 3] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 3];
                      retvalp[r * num_col + mc + 4] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 4];
                      retvalp[r * num_col + mc + 5] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 5];
                      retvalp[r * num_col + mc + 6] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 6];
                      retvalp[r * num_col + mc + 7] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 7];
                  }
              }
              for(; mc + 3 < num_col; mc += 4) {
                  for(O c = offp[r]; c < offp[r+1]; c++) {
                      retvalp[r * num_col + mc] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc];
                      retvalp[r * num_col + mc + 1] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 1];
                      retvalp[r * num_col + mc + 2] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 2];
                      retvalp[r * num_col + mc + 3] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc + 3];
                  }
              }
              for(; mc < num_col; mc++) {
                  for(O c = offp[r]; c < offp[r+1]; c++) {
                      retvalp[r * num_col + mc] +=
                              valp[c] * vvalp[idxp[c] * num_col + mc];
                  }
              }
          }
          /*
        for(size_t r = 0; r < mat.local_num_row; r++) {
        O c = offp[r];
        for(; c + CRS_VLEN < offp[r+1]; c += CRS_VLEN) {
        for(size_t mc = 0; mc < num_col; mc++) {
        //#pragma cdir on_adb(vvalp)
        for(O i = 0; i < CRS_VLEN; i++) {
        retvalp[r * num_col + mc] +=
        valp[c+i] * vvalp[idxp[c+i] * num_col + mc];
        }
        }
        }
        for(size_t mc = 0; mc < num_col; mc++) {
        //#pragma cdir on_adb(vvalp)
        for(O i = 0; c + i < offp[r+1]; i++) {
        retvalp[r * num_col + mc] +=
        valp[c+i] * vvalp[idxp[c+i] * num_col + mc];
        }
        }
        c = offp[r+1];
        }
      */
      } else {
          crs_matrix_spmm_impl2(mat, retvalp, vvalp, num_col);
      }
  }
  
  
  
  #endif
  