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