105  bicsb_gespmv<PTDD>(matmul.op1, matmul.op2.arr, arr);

106  }

107  else

108  {

109  cout<< "Detected noncompliant matvec..." << endl;

110  }

111  return *this;

112 }

113 

114 // populate the vector with random entries

115 // currently, only works for T "double" and "float"

116 template <class T, class ITYPE>

117 void Spvec<T,ITYPE>::fillrandom()

118 {

119 #if (__GNUC__ == 4 && (__GNUC_MINOR__ < 7) )

120  RandGen G;

121  for(ITYPE i=0; i< n; ++i)

122  {

123  arr[i] = G.RandReal();

124  }

125 #else

126  std::uniform_real_distribution<T> distribution(0.0f, 1.0f); //Values between 0 and 1

127  std::mt19937 engine; // Mersenne twister MT19937

128  auto generator = std::bind(distribution, engine);

129  std::generate_n(arr, n, generator);

130 #endif

131 }

132 

133 // populate the vector with zeros

134 template <class T, class ITYPE>

135 void Spvec<T,ITYPE>::fillzero()

136 {

137  for(ITYPE i=0; i< n; ++i)

138  {

139  arr[i] = 0;

140  }

141 }

142 

143 template <typename NT, typename IT>

144 void Verify(Spvec<NT, IT> & control, Spvec<NT, IT> & test, string name, IT m)

145 {

146  vector<NT>error(m);

147  std::transform(&control[0], (&control[0])+m, &test[0], error.begin(), absdiff<NT>());

148  auto maxelement = std::max_element(error.begin(), error.end());

149  cout << "Max error is: " << *maxelement << " on y[" << maxelement-error.begin()<<"]=" << test[maxelement-error.begin()] << endl;

150  NT machEps = machineEpsilon<NT>();

151  cout << "Absolute machine epsilon is: " << machEps <<" and y[" << maxelement-error.begin() << "]*EPSILON becomes "

152  << machEps * test[maxelement-error.begin()] << endl;

153 

154  NT sqrtm = sqrt(static_cast<NT>(m));

155  cout << "sqrt(n) * relative error is: " << abs(machEps * test[maxelement-error.begin()]) * sqrtm << endl;

156  if ( (abs(machEps * test[maxelement-error.begin()]) * sqrtm) < abs(*maxelement))

157  {

158  cout << "*** ATTENTION ***: error is more than sqrt(n) times the relative machine epsilon" << endl;

159  }

160 

161 #ifdef DEBUG

162  cout << "<index, control, test>: \n";

163  for(IT i=0; i<m; ++i)

164  {

165  if(error[i] > abs(sqrtm * machEps * test[i]))

166  {

167  cout << i << "\t" << control[i] << " " << test[i] << "\n";

168  }

169  }

170 #endif

171 }

172 

173