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