1// ========================================================================== 2// $Source: /var/lib/cvs/Givaro/src/library/matrix/givmatsparseops.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: T. Gautier 9// $Id: givmatsparseops.inl,v 1.3 2011-02-02 16:23:56 briceboyer Exp $ 10// ========================================================================== 11 12#error "this looks very much like dead code" 13 14namespace Givaro { 15#pragma message "#warning this file will probably not compile" 16 17 // -- map of a unary operator, with operator()( Type_t& res ) 18 // res and u could be aliases if OP permits it 19 template<class Domain> 20 template<class UNOP> 21 inline void MatrixDom<Domain,Sparse>:: 22 map ( Rep& res, UNOP& op ) const 23 { 24 size_t sz = res._data.size(); 25 for (size_t i=0; i<sz; ++i) op(res._data[i]); 26 } 27 28 29 template<class Domain> 30 template<class UNOP> 31 inline void MatrixDom<Domain,Sparse>:: 32 map ( Rep& res, UNOP& op, const Rep& u ) const 33 { 34 res.resize(nrow(u), ncol(u)); 35 res._index.resize(u._index.size()); 36 res._data.resize(u._data.size()); 37 size_t sz = res._data.size(); 38 for (size_t i=0; i<sz; ++i) { 39 res._index[i] = u._index[i]; 40 op(res._data[i], u._data[i]); 41 } 42 } 43 44 // -- Comparaizon 45 template<class Domain> 46 inline int MatrixDom<Domain,Sparse>::areEqual 47 ( const Rep& P, const Rep& Q) const 48 { 49 if (nrow(P) != nrow(Q)) return 0; 50 if (ncol(P) != ncol(Q)) return 0; 51 size_t sz = P._data.size(); 52 if (sz != Q._data.size()) return 0; 53 for(size_t i=0; i<sz; ++i) { 54 if (P._index[i] != Q._index[i]) return 0; 55 if (_domain.areNEqual(P._data[i] != Q._data[i])) return 0; 56 } 57 return 1; 58 } 59 60 template<class Domain> 61 inline int MatrixDom<Domain,Sparse>::areNEqual 62 ( const Rep& P, const Rep& Q) const 63 { 64 return !areEqual(P,Q); 65 } 66 67 template<class Domain> 68 inline int MatrixDom<Domain,Sparse>::isZero ( const Rep& P ) const 69 { 70 if (nrow(P) == 0) return 1; // -- col souhld be 0 71 if (ncol(P) == 0) { cerr << " Error: inconsistent data structure"; 72 return 1; } // -- row souhld be 0 73 size_t sz = P._data.size(); 74 if (sz == 0) return 1; 75 for(size_t i=0; i<sz; ++i) 76 if (!_domain.isZero(P._data[i])) return 0; 77 return 1; 78 } 79 80 template<class Domain> 81 inline void MatrixDom<Domain,Sparse>::mulin 82 ( Rep& res, const Type_t& u ) const 83 { 84 Curried2<MulOp<Domain> > op(_domain, u); 85 map(res, op); 86 } 87 88 template<class Domain> 89 inline void MatrixDom<Domain,Sparse>::mul 90 ( Rep& res, const Type_t& u, const Rep& v ) const 91 { 92 Curried1<MulOp<Domain> > op(_domain, u); 93 map(res, op, v); 94 } 95 96 template<class Domain> 97 inline void MatrixDom<Domain,Sparse>::mul 98 ( Rep& res, const Rep& u, const Type_t& v ) const 99 { 100 Curried2<MulOp<Domain> > op(_domain, v); 101 map(res, op, u); 102 } 103 104 105 template<class Domain> 106 inline void MatrixDom<Domain,Sparse>::neg ( Rep& R, const Rep& P ) const 107 { 108 size_t sz = P._data.size(); 109 if (R._data.size() != sz) { 110 R.resize(nrow(P), ncol(P)); 111 R._index.copy(P._index); 112 R._data.resize(P._data); 113 } 114 for(size_t i=0; i<sz; ++i) _domain.neg(R._data[i], P._data[i]); 115 } 116 117 // VD is the vector domain for res and u 118 template<class Domain> 119 inline void MatrixDom<Domain,Sparse>::mul 120 ( VectorDom<Domain,Dense>::Rep& R, 121 const Rep& M, 122 const VectorDom<Domain,Dense>& VD, 123 const VectorDom<Domain,Dense>::Rep& U ) const 124 { 125 Indice_t i,j; 126 Indice_t irow, erow; 127 for (i = nrow(M); --i>=0;) 128 { 129 // -- update the i-th row of R 130 _domain.assign(R[i], _domain.zero); 131 irow = M._rows[i]; // - first index of the ith row 132 erow = M._rows[i+1]; // - last index of the ith row 133 for (; irow != erow; ++irow) 134 _domain.axpyin(R[i], M._data[irow], U[M._index[irow]]); 135 } 136 } 137 138 template<class Domain> 139 inline void MatrixDom<Domain,Sparse>::multrans 140 ( typename VectorDom<Domain,Dense>::Rep& R, 141 const Rep& M, 142 const VectorDom<Domain,Dense>& VS, 143 const typename VectorDom<Domain,Dense>::Rep& U ) const 144 { 145 Indice_t i,j; 146 Indice_t irow, erow; 147 for (i = 0; i<ncol(M); ++i) 148 _domain.assign(R[i], _domain.zero); 149 for (i = nrow(M); --i>=0;) 150 { 151 // -- update the i-th row of R 152 irow = M._rows[i]; // - first index of the ith row 153 erow = M._rows[i+1]; // - last index of the ith row 154 for (; irow != erow; ++irow) 155 _domain.axpyin(R[M._index[irow]], U[i], M._data[irow]); 156 } 157 } 158 159 160 161 template<class Domain> 162 void MatrixDom<Domain, Sparse>::compact 163 ( Rep& Ms, 164 const MatrixDom<Domain, Dense>& MD, 165 const MatrixDom<Domain, Dense>::Rep& Md) 166 { 167 // -- Should compare _domain and MD.subdomain(): to be equal! 168 169 // -- Iterate by row M and store non nul entry 170 Indice_t next_rows =0; // next entry to write in _storage._rows 171 Indice_t next_val = 0; // next entry to write in _data & _index 172 Indice_t nrows = MD.nrow(Md); 173 Indice_t ncols = MD.ncol(Md); 174 size_t size = 0; // size of _data and _index 175 Ms.resize( nrows, ncols ); 176 Indice_t i,j; 177 Ms._rows[0] = 0; 178 Type_t val; 179 for (i=0; i<nrows; ++i) 180 { 181 Ms._rows[1+i] = Ms._rows[i]; 182 for (j=0; j<ncols; ++j) 183 { 184 _domain.assign(val, Md(i,j)); 185 if ( !_domain.isZero(val) ) 186 { 187 Ms._data.resize(size+1); 188 Ms._index.resize(size+1); 189 _domain.assign(Ms._data[size], val); 190 Ms._index[size] = j; 191 ++size; Ms._rows[i+1] +=1; 192 } 193 } 194 195 cout << "i:" << i << "----> "; 196 for (Indice_t k=0; k<=nrows; ++k) 197 cout << "," << Ms._rows[k]; 198 cout << endl; 199 } 200 } 201 202 template <class Domain> 203 inline ostream& MatrixDom<Domain,Sparse>::write( ostream& sout ) const 204 { 205 return _domain.write(sout << '(') << ",Sparse)"; 206 } 207 208 template <class Domain> 209 inline istream& MatrixDom<Domain,Sparse>::read( istream& sin ) 210 { 211 char ch; 212 sin >> std::ws >> ch; 213 if (ch != '(') 214 GivError::throw_error( 215 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no '('")); 216 217 // -- read subdomain: 218 _domain.read(sin); 219 220 // -- read , 221 sin >> std::ws >> ch; 222 if (ch != ',') 223 GivError::throw_error( 224 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ','")); 225 226 // -- read dense: 227 char name[7]; 228 sin >> std::ws; sin.read(name, 6); name[6] = (char)0; 229 if (strcmp(name, "Sparse") !=0) 230 GivError::throw_error( 231 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no 'Sparse'")); 232 233 // -- read ) 234 sin >> std::ws >> ch; 235 if (ch != ')') 236 GivError::throw_error( 237 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ')'")); 238 return sin; 239 } 240 241 template <class Domain> 242 ostream& MatrixDom<Domain,Sparse>::write( ostream& sout, const Rep& A) const 243 { 244 sout << '[' << nrow(A) << ',' << ncol(A) << ','; 245 IntDom D; 246 { 247 VectorDom<IntDom,Dense> VDi (D); 248 VDi.write(sout, A._rows) << ','; 249 VDi.write(sout, A._index) << ','; 250 }{ 251 const VectorDom<Domain,Dense> VD (_domain); 252 VD.write(sout, A._data); 253 } 254 return sout << ']'; 255 } 256 257 template <class Domain> 258 istream& MatrixDom<Domain,Sparse>::read( istream& sin, Rep& R) const 259 { 260 long nr,nc; 261 char ch; 262 sin >> std::ws >> ch; 263 if (ch != '[') 264 GivError::throw_error( 265 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no '['")); 266 267 sin >> std::ws >> nr; 268 sin >> std::ws >> ch; 269 if (ch != ',') 270 GivError::throw_error( 271 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ','")); 272 sin >> std::ws >> nc; 273 274 R.resize(nr,nc); 275 276 sin >> std::ws >> ch; 277 if (ch != ',') 278 GivError::throw_error( 279 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ','")); 280 281 VectorDom<IntDom,Dense> VDi = IntDom(); 282 VectorDom<Domain,Dense> VD (_domain); 283 VDi.read(sin, R._rows) >> std::ws >> ch; 284 if (ch != ',') 285 GivError::throw_error( 286 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ','")); 287 VDi.read(sin, R._index) >> std::ws >> ch; 288 if (ch != ',') 289 GivError::throw_error( 290 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ','")); 291 VD.read(sin, R._data) >> std::ws >> ch; 292 if (ch != ']') 293 GivError::throw_error( 294 GivBadFormat("MatrixDom<Domain,Sparse>::read: syntax error no ']'")); 295 return sin; 296 } 297 298} // Givaro 299 300/* -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 301// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 302