← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to extern/eltopo/common/sparse/sparsematrix.h

  • Committer: Package Import Robot
  • Author(s): Matteo F. Vescovi
  • Date: 2012-07-23 08:54:18 UTC
  • mfrom: (14.2.16 sid)
  • mto: (14.2.19 sid)
  • mto: This revision was merged to the branch mainline in revision 42.
  • Revision ID: package-import@ubuntu.com-20120723085418-9foz30v6afaf5ffs
Tags: 2.63a-2
http://bugs.debian.org/658075
* debian/: Cycles support added (Closes: #658075)
  For now, this top feature has been enabled only
  on [any-amd64 any-i386] architectures because
  of OpenImageIO failing on all others
* debian/: scripts installation path changed
  from /usr/lib to /usr/share:
  + debian/patches/: patchset re-worked for path changing
  + debian/control: "Breaks" field added on yafaray-exporter

Show diffs side-by-side

added added

removed removed

Lines of Context:
extern/eltopo/common/sparse/sparsematrix.h
 
1
#ifndef SPARSEMATRIX_H
 
2
#define SPARSEMATRIX_H
 
3
 
 
4
// using a variant on CSR format
 
5
 
 
6
#include <iostream>
 
7
#include <vector>
 
8
#include "util.h"
 
9
#include <assert.h>
 
10
 
 
11
template<class T>
 
12
struct SparseMatrix
 
13
{
 
14
   unsigned int m, n; // dimensions
 
15
   std::vector<std::vector<unsigned int> > index; // for each row, a list of all column indices (sorted)
 
16
   std::vector<std::vector<T> > value; // values corresponding to index
 
17
 
 
18
   explicit SparseMatrix(unsigned int m_=0)
 
19
      : m(m_), n(m_), index(m_), value(m_)
 
20
   {
 
21
      const unsigned int expected_nonzeros_per_row=7;
 
22
      for(unsigned int i=0; i<m; ++i){
 
23
         index[i].reserve(expected_nonzeros_per_row);
 
24
         value[i].reserve(expected_nonzeros_per_row);
 
25
      }
 
26
   }
 
27
 
 
28
   explicit SparseMatrix(unsigned int m_, unsigned int n_, unsigned int expected_nonzeros_per_row)
 
29
      : m(m_), n(n_), index(m_), value(m_)
 
30
   {
 
31
      for(unsigned int i=0; i<m; ++i){
 
32
         index[i].reserve(expected_nonzeros_per_row);
 
33
         value[i].reserve(expected_nonzeros_per_row);
 
34
      }
 
35
   }
 
36
 
 
37
   void clear(void)
 
38
   {
 
39
      m=n=0;
 
40
      index.clear();
 
41
      value.clear();
 
42
   }
 
43
 
 
44
   void zero(void)
 
45
   {
 
46
      for(unsigned int i=0; i<m; ++i){
 
47
         index[i].resize(0);
 
48
         value[i].resize(0);
 
49
      }
 
50
   }
 
51
 
 
52
   void resize(unsigned int m_, unsigned int n_)
 
53
   {
 
54
      m=m_;
 
55
      n=n_;
 
56
      index.resize(m);
 
57
      value.resize(m);
 
58
   }
 
59
 
 
60
   T operator()(unsigned int i, unsigned int j) const
 
61
   {
 
62
      assert(i<m && j<n);
 
63
      for(unsigned int k=0; k<index[i].size(); ++k){
 
64
         if(index[i][k]==j) return value[i][k];
 
65
         else if(index[i][k]>j) return 0;
 
66
      }
 
67
      return 0;
 
68
   }
 
69
 
 
70
   void set_element(unsigned int i, unsigned int j, T new_value)
 
71
   {
 
72
      assert(i<m && j<n);
 
73
      unsigned int k=0;
 
74
      for(; k<index[i].size(); ++k){
 
75
         if(index[i][k]==j){
 
76
            value[i][k]=new_value;
 
77
            return;
 
78
         }else if(index[i][k]>j){
 
79
            insert(index[i], k, j);
 
80
            insert(value[i], k, new_value);
 
81
            return;
 
82
         }
 
83
      }
 
84
      index[i].push_back(j);
 
85
      value[i].push_back(new_value);
 
86
   }
 
87
 
 
88
   void add_to_element(unsigned int i, unsigned int j, T increment_value)
 
89
   {
 
90
      assert(i<m && j<n);
 
91
      unsigned int k=0;
 
92
      for(; k<index[i].size(); ++k){
 
93
         if(index[i][k]==j){
 
94
            value[i][k]+=increment_value;
 
95
            return;
 
96
         }else if(index[i][k]>j){
 
97
            insert(index[i], k, j);
 
98
            insert(value[i], k, increment_value);
 
99
            return;
 
100
         }
 
101
      }
 
102
      index[i].push_back(j);
 
103
      value[i].push_back(increment_value);
 
104
   }
 
105
 
 
106
   // assumes indices is already sorted
 
107
   void add_sparse_row(unsigned int i, const std::vector<unsigned int> &indices, const std::vector<T> &values)
 
108
   {
 
109
      assert(i<m);
 
110
      unsigned int j=0, k=0;
 
111
      while(j<indices.size() && k<index[i].size()){
 
112
         if(index[i][k]<indices[j]){
 
113
            ++k;
 
114
         }else if(index[i][k]>indices[j]){
 
115
            insert(index[i], k, indices[j]);
 
116
            insert(value[i], k, values[j]);
 
117
            ++j;
 
118
         }else{
 
119
            value[i][k]+=values[j];
 
120
            ++j;
 
121
            ++k;
 
122
         }
 
123
      }
 
124
      for(; j<indices.size(); ++j){
 
125
         index[i].push_back(indices[j]);
 
126
         value[i].push_back(values[j]);
 
127
      }
 
128
   }
 
129
 
 
130
   // assumes indices is already sorted
 
131
   void add_sparse_row(unsigned int i, const std::vector<unsigned int> &indices, T multiplier, const std::vector<T> &values)
 
132
   {
 
133
      assert(i<m);
 
134
      unsigned int j=0, k=0;
 
135
      while(j<indices.size() && k<index[i].size()){
 
136
         if(index[i][k]<indices[j]){
 
137
            ++k;
 
138
         }else if(index[i][k]>indices[j]){
 
139
            insert(index[i], k, indices[j]);
 
140
            insert(value[i], k, multiplier*values[j]);
 
141
            ++j;
 
142
         }else{
 
143
            value[i][k]+=multiplier*values[j];
 
144
            ++j;
 
145
            ++k;
 
146
         }
 
147
      }
 
148
      for(;j<indices.size(); ++j){
 
149
         index[i].push_back(indices[j]);
 
150
         value[i].push_back(multiplier*values[j]);
 
151
      }
 
152
   }
 
153
 
 
154
   // assumes matrix has symmetric structure - so the indices in row i tell us which columns to delete i from
 
155
   void symmetric_remove_row_and_column(unsigned int i)
 
156
   {
 
157
      assert(m==n && i<m);
 
158
      for(unsigned int a=0; a<index[i].size(); ++a){
 
159
         unsigned int j=index[i][a]; // 
 
160
         for(unsigned int b=0; b<index[j].size(); ++b){
 
161
            if(index[j][b]==i){
 
162
               erase(index[j], b);
 
163
               erase(value[j], b);
 
164
               break;
 
165
            }
 
166
         }
 
167
      }
 
168
      index[i].resize(0);
 
169
      value[i].resize(0);
 
170
   }
 
171
 
 
172
   void write_matlab(std::ostream &output, const char *variable_name)
 
173
   {
 
174
      output<<variable_name<<"=sparse([";
 
175
      for(unsigned int i=0; i<m; ++i){
 
176
         for(unsigned int j=0; j<index[i].size(); ++j){
 
177
            output<<i+1<<" ";
 
178
         }
 
179
      }
 
180
      output<<"],...\n  [";
 
181
      for(unsigned int i=0; i<m; ++i){
 
182
         for(unsigned int j=0; j<index[i].size(); ++j){
 
183
            output<<index[i][j]+1<<" ";
 
184
         }
 
185
      }
 
186
      output<<"],...\n  [";
 
187
      for(unsigned int i=0; i<m; ++i){
 
188
         for(unsigned int j=0; j<value[i].size(); ++j){
 
189
            output<<value[i][j]<<" ";
 
190
         }
 
191
      }
 
192
      output<<"], "<<m<<", "<<n<<");"<<std::endl;
 
193
   }
 
194
};
 
195
 
 
196
typedef SparseMatrix<float> SparseMatrixf;
 
197
typedef SparseMatrix<double> SparseMatrixd;
 
198
 
 
199
// perform result=matrix*x
 
200
template<class T>
 
201
void multiply(const SparseMatrix<T> &matrix, const std::vector<T> &x, std::vector<T> &result)
 
202
{
 
203
   assert(matrix.n==x.size());
 
204
   result.resize(matrix.m);
 
205
   for(unsigned int i=0; i<matrix.m; ++i){
 
206
      result[i]=0;
 
207
      for(unsigned int j=0; j<matrix.index[i].size(); ++j){
 
208
         result[i]+=matrix.value[i][j]*x[matrix.index[i][j]];
 
209
      }
 
210
   }
 
211
}
 
212
 
 
213
// perform result=result-matrix*x
 
214
template<class T>
 
215
void multiply_and_subtract(const SparseMatrix<T> &matrix, const std::vector<T> &x, std::vector<T> &result)
 
216
{
 
217
   assert(matrix.n==x.size());
 
218
   result.resize(matrix.m);
 
219
   for(unsigned int i=0; i<matrix.m; ++i){
 
220
      for(unsigned int j=0; j<matrix.index[i].size(); ++j){
 
221
         result[i]-=matrix.value[i][j]*x[matrix.index[i][j]];
 
222
      }
 
223
   }
 
224
}
 
225
 
 
226
// compute C=A*B (B might equal A)
 
227
template<class T>
 
228
void multiply_sparse_matrices(const SparseMatrix<T> &A, const SparseMatrix<T> &B, SparseMatrix<T> &C)
 
229
{
 
230
   assert(A.n==B.m);
 
231
   C.resize(A.m, B.n);
 
232
   C.zero();
 
233
   for(unsigned int i=0; i<A.m; ++i){
 
234
      for(unsigned int p=0; p<A.index[i].size(); ++p){
 
235
         unsigned int k=A.index[i][p];
 
236
         C.add_sparse_row(i, B.index[k], A.value[i][p], B.value[k]);
 
237
      }
 
238
   }
 
239
}
 
240
 
 
241
// compute C=A^T*B (B might equal A)
 
242
template<class T>
 
243
void multiply_AtB(const SparseMatrix<T> &A, const SparseMatrix<T> &B, SparseMatrix<T> &C)
 
244
{
 
245
   assert(A.m==B.m);
 
246
   C.resize(A.n, B.n);
 
247
   C.zero();
 
248
   for(unsigned int k=0; k<A.m; ++k){
 
249
      for(unsigned int p=0; p<A.index[k].size(); ++p){
 
250
         unsigned int i=A.index[k][p];
 
251
         C.add_sparse_row(i, B.index[k], A.value[k][p], B.value[k]);
 
252
      }
 
253
   }
 
254
}
 
255
 
 
256
// compute C=A*D*B (B might equal A)
 
257
template<class T>
 
258
void multiply_sparse_matrices_with_diagonal_weighting(const SparseMatrix<T> &A, const std::vector<T> &diagD,
 
259
                                                      const SparseMatrix<T> &B, SparseMatrix<T> &C)
 
260
{
 
261
   assert(A.n==B.m);
 
262
   assert(diagD.size()==A.n);
 
263
   C.resize(A.m, B.n);
 
264
   C.zero();
 
265
   for(unsigned int i=0; i<A.m; ++i){
 
266
      for(unsigned int p=0; p<A.index[i].size(); ++p){
 
267
         unsigned int k=A.index[i][p];
 
268
         C.add_sparse_row(i, B.index[k], A.value[i][p]*diagD[k], B.value[k]);
 
269
      }
 
270
   }
 
271
}
 
272
 
 
273
// compute C=A^T*D*B (B might equal A)
 
274
template<class T>
 
275
void multiply_AtDB(const SparseMatrix<T> &A, const std::vector<T>& diagD, const SparseMatrix<T> &B, SparseMatrix<T> &C)
 
276
{
 
277
   assert(A.m==B.m);
 
278
   C.resize(A.n, B.n);
 
279
   C.zero();
 
280
   for(unsigned int k=0; k<A.m; ++k){
 
281
      for(unsigned int p=0; p<A.index[k].size(); ++p){
 
282
         unsigned int i=A.index[k][p];
 
283
         C.add_sparse_row(i, B.index[k], A.value[k][p]*diagD[k], B.value[k]);
 
284
      }
 
285
   }
 
286
}
 
287
 
 
288
// set T to the transpose of A
 
289
template<class T>
 
290
void compute_transpose(const SparseMatrix<T> &A, SparseMatrix<T> &C)
 
291
{
 
292
   C.resize(A.n, A.m);
 
293
   C.zero();
 
294
   for(unsigned int i=0; i<A.m; ++i){
 
295
      for(unsigned int p=0; p<A.index[i].size(); ++p){
 
296
         unsigned int k=A.index[i][p];
 
297
         C.index[k].push_back(i);
 
298
         C.value[k].push_back(A.value[i][p]);
 
299
      }
 
300
   }
 
301
}
 
302
 
 
303
template<class T>
 
304
struct FixedSparseMatrix
 
305
{
 
306
   unsigned int m, n; // dimension
 
307
   std::vector<T> value; // nonzero values row by row
 
308
   std::vector<unsigned int> colindex; // corresponding column indices
 
309
   std::vector<unsigned int> rowstart; // where each row starts in value and colindex (and last entry is one past the end, the number of nonzeros)
 
310
 
 
311
   explicit FixedSparseMatrix(unsigned int m_=0)
 
312
      : m(m_), n(m_), value(0), colindex(0), rowstart(m_+1)
 
313
   {}
 
314
 
 
315
   explicit FixedSparseMatrix(unsigned int m_, unsigned int n_)
 
316
      : m(m_), n(n_), value(0), colindex(0), rowstart(m_+1)
 
317
   {}
 
318
 
 
319
   void clear(void)
 
320
   {
 
321
      m=n=0;
 
322
      value.clear();
 
323
      colindex.clear();
 
324
      rowstart.clear();
 
325
   }
 
326
 
 
327
   void resize(unsigned int m_, unsigned int n_)
 
328
   {
 
329
      m=m_;
 
330
      n=n_;
 
331
      rowstart.resize(m+1);
 
332
   }
 
333
 
 
334
   void construct_from_matrix(const SparseMatrix<T> &matrix)
 
335
   {
 
336
      resize(matrix.m, matrix.n);
 
337
      rowstart[0]=0;
 
338
      for(unsigned int i=0; i<m; ++i){
 
339
         rowstart[i+1]=rowstart[i]+matrix.index[i].size();
 
340
      }
 
341
      value.resize(rowstart[m]);
 
342
      colindex.resize(rowstart[m]);
 
343
      unsigned int j=0;
 
344
      for(unsigned int i=0; i<m; ++i){
 
345
         for(unsigned int k=0; k<matrix.index[i].size(); ++k){
 
346
            value[j]=matrix.value[i][k];
 
347
            colindex[j]=matrix.index[i][k];
 
348
            ++j;
 
349
         }
 
350
      }
 
351
   }
 
352
 
 
353
   void write_matlab(std::ostream &output, const char *variable_name)
 
354
   {
 
355
      output<<variable_name<<"=sparse([";
 
356
      for(unsigned int i=0; i<m; ++i){
 
357
         for(unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j){
 
358
            output<<i+1<<" ";
 
359
         }
 
360
      }
 
361
      output<<"],...\n  [";
 
362
      for(unsigned int i=0; i<m; ++i){
 
363
         for(unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j){
 
364
            output<<colindex[j]+1<<" ";
 
365
         }
 
366
      }
 
367
      output<<"],...\n  [";
 
368
      for(unsigned int i=0; i<m; ++i){
 
369
         for(unsigned int j=rowstart[i]; j<rowstart[i+1]; ++j){
 
370
            output<<value[j]<<" ";
 
371
         }
 
372
      }
 
373
      output<<"], "<<m<<", "<<n<<");"<<std::endl;
 
374
   }
 
375
};
 
376
 
 
377
typedef FixedSparseMatrix<float> FixedSparseMatrixf;
 
378
typedef FixedSparseMatrix<double> FixedSparseMatrixd;
 
379
 
 
380
// perform result=matrix*x
 
381
template<class T>
 
382
void multiply(const FixedSparseMatrix<T> &matrix, const std::vector<T> &x, std::vector<T> &result)
 
383
{
 
384
   assert(matrix.n==x.size());
 
385
   result.resize(matrix.m);
 
386
   for(unsigned int i=0; i<matrix.m; ++i){
 
387
      result[i]=0;
 
388
      for(unsigned int j=matrix.rowstart[i]; j<matrix.rowstart[i+1]; ++j){
 
389
         result[i]+=matrix.value[j]*x[matrix.colindex[j]];
 
390
      }
 
391
   }
 
392
}
 
393
 
 
394
// perform result=result-matrix*x
 
395
template<class T>
 
396
void multiply_and_subtract(const FixedSparseMatrix<T> &matrix, const std::vector<T> &x, std::vector<T> &result)
 
397
{
 
398
   assert(matrix.n==x.size());
 
399
   result.resize(matrix.m);
 
400
   for(unsigned int i=0; i<matrix.m; ++i){
 
401
      for(unsigned int j=matrix.rowstart[i]; j<matrix.rowstart[i+1]; ++j){
 
402
         result[i]-=matrix.value[j]*x[matrix.colindex[j]];
 
403
      }
 
404
   }
 
405
}
 
406
 
 
407
#endif