1/* linbox/blackbox/nag-sparse.h 2 * Copyright (C) 2002 Rich Seagraves 3 * 4 * Written by Rich Seagraves <seagrave@cis.udel.edu> 5 * Modified by Zhendong Wan, -bds 6 * ------------------------------------ 7 * 8 * 9 * ========LICENCE======== 10 * This file is part of the library LinBox. 11 * 12 * LinBox is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU Lesser General Public 14 * License as published by the Free Software Foundation; either 15 * version 2.1 of the License, or (at your option) any later version. 16 * 17 * This library is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 20 * Lesser General Public License for more details. 21 * 22 * You should have received a copy of the GNU Lesser General Public 23 * License along with this library; if not, write to the Free Software 24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 25 * ========LICENCE======== 26 *. 27 */ 28 29namespace LinBox 30{ 31 32 /* 33 // for GF(3^10), apply to packvec will do adds with periodic normalizations 34 template<> 35 template<Field> 36 Packvec<typename Field::Accumulator>& ZeroOne<Field>::apply<Packvec<typename Field::Accumulator>, Packvec<typename Field::Accumulator> > 37 (Pacvec<typename Field::Accumulator>& y, 38 Pacvec<typename Field::Accumulator>& x) 39 { return y; } 40 41 // more stuff for paley graphs 42 template<class OutVector, class InVector> 43 template<> 44 OutVector & ZeroOne<Special3_10Field>::apply(OutVector & y, const InVector & x) const 45 { 46 typedef Special3_10Field Field; 47 typedef typename Packvec<Field::Accumulator> Vector; 48 Vector xx, yy; 49 field().convert(xx, x); 50 apply(yy, xx); 51 field().convert(y, yy); 52 } 53 */ 54 55 template<class Field> 56 template<class OutVector, class InVector> 57 OutVector & ZeroOne<Field>::apply(OutVector & y, const InVector & x) const 58 //OutVector & ZeroOne<Field>::apply(OutVector & y, const InVector & x) 59 { 60 61 //std::cout << endl; 62 //std::cout << " new apply pass " << endl; 63 //std::cout << this->rowdim() << " " << this->coldim() << endl; 64 //std::cout << " inside apply " << endl; 65 //std::cout << x.size() << " " << y.size() << endl; 66 67 linbox_check((y.size()==rowdim())&&(x.size()==coldim())); 68 69 FieldAXPY<Field> accum (field()); 70 71 typename OutVector::iterator yp; 72 typename InVector::const_iterator xp; 73 PointerVector::const_iterator ip; 74 IndexVector::const_iterator jp; 75 76 77 //std::cout << " before sorting " << endl; 78 79 if( !sorted ) 80 { 81 //std::cout << " -- apply: switch sort (" << _rowdim << ", " << _coldim << ", " << nnz() << ", " << sorted << " " << &(this->sorted) << ") -- " << std::endl; 82 switch_sort(); 83 //std::cout << " -- apply: switch sort (" << _rowdim << ", " << _coldim << ", " << nnz() << ", " << sorted << ") -- " << std::endl; 84 } 85 //std::cout << " -- zo apply: " << this << "(" << _rowdim << ", " << _coldim << ", " << nnz() << ")" << std::endl; 86 87 xp=x.begin(); 88 yp=y.begin(); 89 accum.reset(); 90 91 92 //std::cout << " before for loop " << endl; 93 //std::cout << _indexP.end() - _indexP.begin() << " "; 94 //std::cout << _index.end() - _index.begin() << endl; 95 96 //std::cout << " stupid test here " << (*_indexP.begin()) - _index.begin() << endl; 97 //cout << &_index << " of size "; 98 //cout << sizeof(_index) << " "; 99 //cout << &(*_index.begin()) << endl; 100 101 for(ip = _indexP.begin(); ip !=_indexP.end()-1; ++ip, ++yp) 102 //for(ip = _indexP.begin(); ip < _indexP.end()-2; ++ip, ++yp) // zigzag way 103 { 104 //std::cout << " inside the outer for loop " << ip - _indexP.begin() << endl; 105 for(jp = *ip; jp !=*(ip + 1); ++jp) 106 { 107 //std::cout << jp - _index.begin() << endl; 108 accum.accumulate_special( *(xp + (ptrdiff_t)*jp) ); 109 } 110 //std::cout << " accumulate is done for one iteration " << endl; 111 //std::cout << " before accum.get " << yp - y.begin() << endl; 112 accum.get(*yp); 113 //std::cout << " before accum.reset " << endl; 114 accum.reset(); 115 116 /* // zigzag way 117 ++ip;++yp; 118 119 for(jp = *(ip + 1); jp >*ip; --jp) 120 accum.accumulate_special( *(xp + *(jp-1)) ); 121 accum.get(*yp); 122 accum.reset(); 123 */ 124 125 } 126 127 return y; 128 } 129 130 /* if you want to keep two copies for the matrix, one of which is sorted by row, 131 * the other by column, then you want to use this applyTranspose function. In 132 * this case, un-comment this one and comment out the applyTranspose further down 133 */ 134 /* 135 template<class Field> 136 template<class OutVector, class InVector> 137 OutVector & ZeroOne<Field>::applyTranspose(OutVector & y, const InVector & x) const 138 { 139 linbox_check((y.size()==coldim())&&(x.size()==rowdim())); 140 141 FieldAXPY<Field> accum (field()); 142 143 typename OutVector::iterator yp; 144 typename InVector::const_iterator xp; 145 PointerVector::const_iterator ip; 146 IndexVector::const_iterator jp; 147 148 xp=x.begin(); 149 yp=y.begin(); 150 accum.reset(); 151 152 for(ip = _colP.begin(); ip < _colP.end()-1; ++ip, ++yp) 153 { 154 for(jp = *ip; jp <*(ip + 1); ++jp) 155 accum.accumulate_special( *(xp + *jp) ); 156 accum.get(*yp); 157 accum.reset(); 158 } 159 160 return y; 161 162 } 163 */ 164 165 template<class Field> 166 template<class OutVector, class InVector> 167 OutVector & ZeroOne<Field>::applyTranspose(OutVector & y, const InVector & x) const 168 { 169 linbox_check((y.size()==coldim())&&(x.size()==rowdim())); 170 171 FieldAXPY<Field> accum (field()); 172 173 typename OutVector::iterator yp; 174 typename InVector::const_iterator xp; 175 PointerVector::const_iterator ip; 176 IndexVector::const_iterator jp; 177 178 xp=x.begin(); 179 yp=y.begin(); 180 accum.reset(); 181 182 if( sorted ) 183 //{ 184 //std::cout << " -- in apply transpose, before switch sort -- " << std::endl; 185 switch_sort(); 186 //std::cout << " -- in apply transpose, after switch sort -- " << std::endl; 187 //} 188 189 for(ip = _indexP.begin(); ip < _indexP.end()-1; ++ip, ++yp) 190 { 191 for(jp = *ip; jp <*(ip + 1); ++jp) 192 accum.accumulate_special( *(xp +(ptrdiff_t) *jp) ); 193 accum.get(*yp); 194 accum.reset(); 195 } 196 return y; 197 } 198 199}//End of LinBox 200 201// Local Variables: 202// mode: C++ 203// tab-width: 4 204// indent-tabs-mode: nil 205// c-basic-offset: 4 206// End: 207// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 208