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