1// ==========================================================================
2// $Source: /var/lib/cvs/Givaro/src/library/poly1/givpoly1axpy.inl,v $
3// Copyright(c)'1994-2009 by The Givaro group
4// This file is part of Givaro.
5// Givaro is governed by the CeCILL-B license under French law
6// and abiding by the rules of distribution of free software.
7// see the COPYRIGHT file for more details.
8// Authors: J-G. Dumas
9// $Id: givpoly1axpy.inl,v 1.6 2011-02-02 16:23:56 briceboyer Exp $
10// ==========================================================================
11
12#ifndef __GIVARO_poly1axpy_INL
13#define __GIVARO_poly1axpy_INL
14
15namespace Givaro {
16
17    // axpy, axmy, maxpy
18    // J.G.D. 16.11.2006
19    // A lot can be done to optimize those
20    // Except for axpy, axpyin, maxpy with a a Type_t,
21    // all of them use a temporary vector where
22    // a temporary value only would be sufficient.
23
24    template <class Domain>
25    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axpy  (Rep& r, const Rep& a, const Rep& x, const Rep& y) const
26    {
27        return this->addin( this->mul(r,a,x), y );
28    }
29
30
31    template <class Domain>
32    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axpy  (Rep& r, const Type_t& a, const Rep& x, const Rep& y) const
33    {
34        typename Rep::const_iterator ix = x.begin(), iy = y.begin();
35        if (y.size() > x.size()) {
36            r.resize(y.size());
37            typename Rep::iterator ir = r.begin();
38            for( ; ix != x.end(); ++ir, ++ix, ++iy)
39                this->_domain.axpy(*ir, a, *ix, *iy);
40            for( ; ir != r.end(); ++ir, ++iy)
41                this->_domain.assign(*ir, *iy);
42        } else {
43            r.resize(x.size());
44            typename Rep::iterator ir = r.begin();
45            for( ; iy != y.end(); ++ir, ++ix, ++iy)
46                this->_domain.axpy(*ir, a, *ix, *iy);
47            for( ; ir != r.end(); ++ir, ++ix)
48                this->_domain.mul(*ir, a, *ix);
49        }
50        return r;
51    }
52
53    template <class Domain>
54    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axpyin(Rep& r, const Rep& a, const Rep& x) const
55    {
56        Rep tmp; this->init(tmp);
57        this->assign(tmp,r);
58        return this->axpy(r,a,x,tmp);
59    }
60
61    template <class Domain>
62    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axpyin(Rep& r, const Type_t& a, const Rep& x) const{
63        typename Rep::const_iterator ix = x.begin();
64        if (x.size() > r.size()) {
65            for(typename Rep::iterator ir = r.begin() ; ir != r.end(); ++ir, ++ix)
66                this->_domain.axpyin(*ir, a, *ix);
67            Type_t tmp;
68            for( ; ix != x.end(); ++ix)
69                r.push_back( this->_domain.mul(tmp, a, *ix) );
70        } else {
71            for(typename Rep::iterator ir = r.begin() ; ix != x.end(); ++ir, ++ix)
72                this->_domain.axpyin(*ir, a, *ix);
73        }
74        return r;
75    }
76    // -- maxpy: r <- c - a * b
77    template <class Domain>
78    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::maxpy  (Rep& r, const Rep& a, const Rep& b, const Rep& c) const{
79        Rep tmp; this->init(tmp);
80        return this->sub(r,c,this->mul(tmp,a,b));
81    }
82
83    template <class Domain>
84    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::maxpy  (Rep& r, const Type_t& a, const Rep& b, const Rep& c) const{
85        size_t sC = c.size();
86        size_t sB = b.size();
87        size_t sR = r.size();
88        if (sB == 0) { r.copy(c); return r; }
89        if (sC == 0) { return this->negin( this->mul(r,a,b) ); }
90        size_t i, max = sC < sB ? sB : sC;
91        if (sR != max) r.resize(max);
92        if (sC < sB)
93        {
94            for (i=0; i<sC; ++i) _domain.maxpy(r[i], a, b[i], c[i]);
95            for (; i<sB; ++i) _domain.negin( _domain.mul(r[i], a, b[i]) );
96        }
97        else {
98            for (i=0; i<sB; ++i) _domain.maxpy(r[i], a, b[i], c[i]);
99            for (; i<sC; ++i) _domain.assign(r[i], c[i]);
100        }
101        return r;
102    }
103
104    // -- maxpyin: r -= a*b
105    template <class Domain>
106    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::maxpyin(Rep& r, const Rep& a, const Rep& b) const{
107        Rep tmp; this->init(tmp);
108        return this->subin(r, this->mul(tmp,a,b));
109    }
110    template <class Domain>
111    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::maxpyin(Rep& r, const Type_t& a, const Rep& b) const{
112        Rep tmp; this->init(tmp);
113        return this->subin(r, this->mul(tmp,a,b));
114    }
115    // -- axmy: r <- a * x - y
116    template <class Domain>
117    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axmy  (Rep& r, const Rep& a, const Rep& x, const Rep& y) const{
118        return this->subin(this->mul(r, a, x),y);
119    }
120    template <class Domain>
121    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axmy  (Rep& r, const Type_t& a, const Rep& x, const Rep& y) const{
122        return this->subin(this->mul(r, a, x),y);
123    }
124    // -- axmyin: r <- a * x - r
125    template <class Domain>
126    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axmyin (Rep& r, const Rep& a, const Rep& x) const
127    {
128        this->maxpyin(r, a, x);
129        return this->negin(r);
130    }
131    template <class Domain>
132    inline typename Poly1Dom<Domain,Dense>::Rep& Poly1Dom<Domain,Dense>::axmyin (Rep& r, const Type_t& a, const Rep& x) const
133    {
134        this->maxpyin(r, a, x);
135        return this->negin(r);
136    }
137
138} // Givaro
139
140#endif // __GIVARO_poly1axpy_INL
141/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
142// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
143