← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

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

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

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