1/* 2 * Copyright (C) 2014 the FFLAS-FFPACK group 3 * 4 * Written by Brice Boyer (briceboyer) <boyer.brice@gmail.com> 5 * 6 * 7 * ========LICENCE======== 8 * This file is part of the library FFLAS-FFPACK. 9 * 10 * FFLAS-FFPACK is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 23 * ========LICENCE======== 24 *. 25 */ 26 27/** @file fflas/fflas_level1.h 28 * @brief Vector operations 29 * or anything of \f$n\f$ complexity 30 */ 31 32#ifndef __FFLASFFPACK_fflas_fflas_level1_INL 33#define __FFLASFFPACK_fflas_fflas_level1_INL 34 35namespace FFLAS { 36 37 38 39 //--------------------------------------------------------------------- 40 // Level 1 routines 41 //--------------------------------------------------------------------- 42 43 /** freduce 44 * \f$x \gets x mod F\f$. 45 * @param F field 46 * @param n size of the vectors 47 * \param X vector in \p F 48 * \param incX stride of \p X 49 * @bug use cblas_(d)scal when possible 50 */ 51 template<class Field> 52 void 53 freduce (const Field& F, const size_t n, 54 typename Field::Element_ptr X, const size_t incX); 55 56 /** freduce 57 * \f$x \gets y mod F\f$. 58 * @param F field 59 * @param n size of the vectors 60 * \param Y vector of \p Element 61 * \param incY stride of \p Y 62 * \param X vector in \p F 63 * \param incX stride of \p X 64 * @bug use cblas_(d)scal when possible 65 */ 66 template<class Field> 67 void 68 freduce (const Field& F, const size_t n, 69 typename Field::ConstElement_ptr Y, const size_t incY, 70 typename Field::Element_ptr X, const size_t incX); 71 72 /** finit 73 * \f$x \gets y mod F\f$. 74 * @param F field 75 * @param n size of the vectors 76 * \param Y vector of \p OtherElement 77 * \param incY stride of \p Y 78 * \param X vector in \p F 79 * \param incX stride of \p X 80 * @bug use cblas_(d)scal when possible 81 */ 82 template<class Field, class OtherElement_ptr> 83 void 84 finit (const Field& F, const size_t n, 85 const OtherElement_ptr Y, const size_t incY, 86 typename Field::Element_ptr X, const size_t incX); 87 88 /** finit 89 * Initializes \p X in \p F$. 90 * @param F field 91 * @param n size of the vectors 92 * \param X vector in \p F 93 * \param incX stride of \p X 94 */ 95 template<class Field> 96 void 97 finit (const Field& F, const size_t n, 98 typename Field::Element_ptr X, const size_t incX); 99 100 /** fconvert 101 * \f$x \gets y mod F\f$. 102 * @param F field 103 * @param n size of the vectors 104 * \param Y vector of \p F 105 * \param incY stride of \p Y 106 * \param X vector in \p OtherElement 107 * \param incX stride of \p X 108 * @bug use cblas_(d)scal when possible 109 */ 110 template<class Field, class OtherElement_ptr> 111 void 112 fconvert (const Field& F, const size_t n, 113 OtherElement_ptr X, const size_t incX, 114 typename Field::ConstElement_ptr Y, const size_t incY) 115 { 116 OtherElement_ptr Xi = X ; 117 typename Field::ConstElement_ptr Yi = Y ; 118 for (; Xi < X+n*incX; Xi+=incX, Yi += incY ) 119 F.convert( *Xi , *Yi); 120 } 121 122 /** fnegin 123 * \f$x \gets - x\f$. 124 * @param F field 125 * @param n size of the vectors 126 * \param X vector in \p F 127 * \param incX stride of \p X 128 * @bug use cblas_(d)scal when possible 129 */ 130 template<class Field> 131 void 132 fnegin (const Field& F, const size_t n, 133 typename Field::Element_ptr X, const size_t incX) 134 { 135 typename Field::Element_ptr Xi = X ; 136 for (; Xi < X+n*incX; Xi+=incX ) 137 F.negin( *Xi ); 138 } 139 140 /** fneg 141 * \f$x \gets - y\f$. 142 * @param F field 143 * @param n size of the vectors 144 * \param X vector in \p F 145 * \param incX stride of \p X 146 * \param Y vector in \p F 147 * \param incY stride of \p Y 148 * @bug use cblas_(d)scal when possible 149 */ 150 template<class Field> 151 void 152 fneg (const Field& F, const size_t n, 153 typename Field::ConstElement_ptr Y, const size_t incY, 154 typename Field::Element_ptr X, const size_t incX) 155 { 156 typename Field::Element_ptr Xi = X ; 157 typename Field::ConstElement_ptr Yi = Y ; 158 for (; Xi < X+n*incX; Xi+=incX,Yi+=incY ) 159 F.neg( *Xi, *Yi ); 160 } 161 162 /** \brief fzero : \f$A \gets 0 \f$. 163 * @param F field 164 * @param n number of elements to zero 165 * \param X vector in \p F 166 * \param incX stride of \p X 167 */ 168 template<class Field> 169 void 170 fzero (const Field& F, const size_t n, 171 typename Field::Element_ptr X, const size_t incX) 172 { 173 if (incX == 1) { // contigous data 174 // memset(X,(int)F.zero,n); // might be bogus ? 175 for (size_t i = 0 ; i < n ; ++i) 176 F.assign(*(X+i), F.zero); 177 178 } 179 else { // not contiguous (strided) 180 for (size_t i = 0 ; i < n ; ++i) 181 F.assign(*(X+i*incX), F.zero); 182 } 183 } 184 185 /** \brief frand : \f$A \gets random \f$. 186 * @param F field 187 * @param G randomiterator 188 * @param n number of elements to randomize 189 * \param X vector in \p F 190 * \param incX stride of \p X 191 */ 192 template<class Field, class RandIter> 193 void 194 frand (const Field& F, RandIter& G, const size_t n, 195 typename Field::Element_ptr X, const size_t incX) 196 { 197 if (incX == 1) { // contigous data 198 // memset(X,(int)F.zero,n); // might be bogus ? 199 for (size_t i = 0 ; i < n ; ++i) 200 G.random(*(X+i)); 201 } 202 else { // not contiguous (strided) 203 for (size_t i = 0 ; i < n ; ++i) 204 G.random(*(X+i*incX)); 205 } 206 } 207 208 /** \brief fiszero : test \f$X = 0 \f$. 209 * @param F field 210 * @param n vector dimension 211 * \param X vector in \p F 212 * \param incX increment of \p X 213 */ 214 template<class Field> 215 bool 216 fiszero (const Field& F, const size_t n, 217 typename Field::ConstElement_ptr X, const size_t incX) 218 { 219 bool res=true; 220 typename Field::ConstElement_ptr Xi = X; 221 for (size_t i = 0 ; i < n ; ++i, Xi += incX) 222 res = res && F.isZero (*Xi); 223 return res; 224 } 225 226 /** \brief fequal : test \f$X = Y \f$. 227 * @param F field 228 * @param n vector dimension 229 * \param X vector in \p F 230 * \param incX increment of \p X 231 * \param Y vector in \p F 232 * \param incY increment of \p Y 233 */ 234 template<class Field> 235 bool 236 fequal (const Field& F, const size_t n, 237 typename Field::ConstElement_ptr X, const size_t incX, 238 typename Field::ConstElement_ptr Y, const size_t incY) 239 { 240 bool res=true; 241 for (size_t i = 0 ; i < n ; ++i) 242 res = res && F.areEqual (X [i*incX], Y [i*incY]); 243 return res; 244 } 245 246 /** \brief fassign : \f$x \gets y \f$. 247 * X is preallocated 248 * @todo variant for triagular matrix 249 * @param F field 250 * @param N size of the vectors 251 * \param [out] X vector in \p F 252 * \param incX stride of \p X 253 * \param [in] Y vector in \p F 254 * \param incY stride of \p Y 255 */ 256 template<class Field> 257 void 258 fassign (const Field& F, const size_t N, 259 typename Field::ConstElement_ptr Y, const size_t incY , 260 typename Field::Element_ptr X, const size_t incX); 261 262 263 /** fscalin 264 * \f$x \gets \alpha \cdot x\f$. 265 * @param F field 266 * @param n size of the vectors 267 * @param alpha scalar 268 * \param X vector in \p F 269 * \param incX stride of \p X 270 * @bug use cblas_(d)scal when possible 271 * @internal 272 * @todo check if comparison with +/-1,0 is necessary. 273 */ 274 template<class Field> 275 void 276 fscalin (const Field& F, const size_t n, const typename Field::Element alpha, 277 typename Field::Element_ptr X, const size_t incX); 278 279 280 /** fscal 281 * \f$y \gets \alpha \cdot x\f$. 282 * @param F field 283 * @param n size of the vectors 284 * @param alpha scalar 285 * \param[in] X vector in \p F 286 * \param incX stride of \p X 287 * \param[out] Y vector in \p F 288 * \param incY stride of \p Y 289 * @bug use cblas_(d)scal when possible 290 * @internal 291 * @todo check if comparison with +/-1,0 is necessary. 292 */ 293 template<class Field> 294 void 295 fscal (const Field& F, const size_t n 296 , const typename Field::Element alpha 297 , typename Field::ConstElement_ptr X, const size_t incX 298 , typename Field::Element_ptr Y, const size_t incY); 299 300 301 302 /** \brief faxpy : \f$y \gets \alpha \cdot x + y\f$. 303 * @param F field 304 * @param N size of the vectors 305 * @param alpha scalar 306 * \param[in] X vector in \p F 307 * \param incX stride of \p X 308 * \param[in,out] Y vector in \p F 309 * \param incY stride of \p Y 310 */ 311 template<class Field> 312 void 313 faxpy (const Field& F, const size_t N, 314 const typename Field::Element alpha, 315 typename Field::ConstElement_ptr X, const size_t incX, 316 typename Field::Element_ptr Y, const size_t incY ); 317 318 /** \brief faxpby : \f$y \gets \alpha \cdot x + \beta \cdot y\f$. 319 * @param F field 320 * @param N size of the vectors 321 * @param alpha scalar 322 * \param[in] X vector in \p F 323 * \param incX stride of \p X 324 * \param beta scalar 325 * \param[in,out] Y vector in \p F 326 * \param incY stride of \p Y 327 * \note this is a catlas function 328 */ 329 template<class Field> 330 void 331 faxpby (const Field& F, const size_t N, 332 const typename Field::Element alpha, 333 typename Field::ConstElement_ptr X, const size_t incX, 334 const typename Field::Element beta, 335 typename Field::Element_ptr Y, const size_t incY ); 336 337 338 /** \brief fdot: dot product \f$x^T y\f$. 339 * @param F field 340 * @param N size of the vectors 341 * \param X vector in \p F 342 * \param incX stride of \p X 343 * \param Y vector in \p F 344 * \param incY stride of \p Y 345 */ 346 template<class Field> 347 typename Field::Element 348 fdot (const Field& F, const size_t N, 349 typename Field::ConstElement_ptr X, const size_t incX, 350 typename Field::ConstElement_ptr Y, const size_t incY); 351 352 353 template<typename Field> 354 typename Field::Element 355 fdot (const Field& F, const size_t N, 356 typename Field::ConstElement_ptr X, const size_t incX, 357 typename Field::ConstElement_ptr Y, const size_t incY, 358 const ParSeqHelper::Sequential seq); 359 360 template<typename Field, class Cut, class Param> 361 typename Field::Element 362 fdot (const Field& F, const size_t N, 363 typename Field::ConstElement_ptr X, const size_t incX, 364 typename Field::ConstElement_ptr Y, const size_t incY, 365 const ParSeqHelper::Parallel<Cut,Param> par); 366 367 /** \brief fswap: \f$ X \leftrightarrow Y\f$. 368 * @bug use cblas_dswap when double 369 * @param F field 370 * @param N size of the vectors 371 * \param X vector in \p F 372 * \param incX stride of \p X 373 * \param Y vector in \p F 374 * \param incY stride of \p Y 375 */ 376 template<class Field> 377 void 378 fswap (const Field& F, const size_t N, typename Field::Element_ptr X, const size_t incX, 379 typename Field::Element_ptr Y, const size_t incY ) 380 { 381 382 typename Field::Element tmp; F.init(tmp); 383 typename Field::Element_ptr Xi = X; 384 typename Field::Element_ptr Yi=Y; 385 for (; Xi < X+N*incX; Xi+=incX, Yi+=incY ){ 386 F.assign( tmp, *Xi ); 387 F.assign( *Xi, *Yi ); 388 F.assign( *Yi, tmp ); 389 } 390 } 391 392 template <class Field> 393 void 394 pfadd (const Field & F, const size_t M, const size_t N, 395 typename Field::ConstElement_ptr A, const size_t lda, 396 typename Field::ConstElement_ptr B, const size_t ldb, 397 typename Field::Element_ptr C, const size_t ldc, const size_t numths); 398 399 template <class Field> 400 void 401 pfsub (const Field & F, const size_t M, const size_t N, 402 typename Field::ConstElement_ptr A, const size_t lda, 403 typename Field::ConstElement_ptr B, const size_t ldb, 404 typename Field::Element_ptr C, const size_t ldc, const size_t numths); 405 406 template <class Field> 407 void 408 pfaddin (const Field& F, const size_t M, const size_t N, 409 typename Field::ConstElement_ptr B, const size_t ldb, 410 typename Field::Element_ptr C, const size_t ldc, size_t numths); 411 412 template <class Field> 413 void 414 pfsubin (const Field& F, const size_t M, const size_t N, 415 typename Field::ConstElement_ptr B, const size_t ldb, 416 typename Field::Element_ptr C, const size_t ldc, size_t numths); 417 418 419 template <class Field> 420 void 421 fadd (const Field& F, const size_t N, 422 typename Field::ConstElement_ptr A, const size_t inca, 423 typename Field::ConstElement_ptr B, const size_t incb, 424 typename Field::Element_ptr C, const size_t incc); 425 426 template <class Field> 427 void 428 fsub (const Field& F, const size_t N, 429 typename Field::ConstElement_ptr A, const size_t inca, 430 typename Field::ConstElement_ptr B, const size_t incb, 431 typename Field::Element_ptr C, const size_t incc); 432 433 template <class Field> 434 void 435 faddin (const Field& F, const size_t N, 436 typename Field::ConstElement_ptr B, const size_t incb, 437 typename Field::Element_ptr C, const size_t incc); 438 439 template <class Field> 440 void 441 fsubin (const Field& F, const size_t N, 442 typename Field::ConstElement_ptr B, const size_t incb, 443 typename Field::Element_ptr C, const size_t incc); 444 445 446 template <class Field> 447 void 448 fadd (const Field& F, const size_t N, 449 typename Field::ConstElement_ptr A, const size_t inca, 450 const typename Field::Element alpha, 451 typename Field::ConstElement_ptr B, const size_t incb, 452 typename Field::Element_ptr C, const size_t incc); 453 454} // FFLAS 455 456 457#endif // __FFLASFFPACK_fflas_fflas_level1_INL 458/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 459// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 460