← Back to branch summary

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

« back to all changes in this revision

Viewing changes to extern/eltopo/common/sparse/incomplete_qr.cpp

  • 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/incomplete_qr.cpp
1
 
#include <cmath>
2
 
#include <list>
3
 
#include <vector>
4
 
#include "incomplete_qr.h"
5
 
 
6
 
/* The MATLAB code for this algorithm:
7
 
for i=1:n
8
 
   q=A(:,i);
9
 
   for j=find(Rpattern(1:i-1,i))'
10
 
      R(j,i)=Q(:,j)'*q;
11
 
      q=q-Q(:,j)*R(j,i);
12
 
   end
13
 
   R(i,i)=norm(q);
14
 
   q=q.*(abs(q)>Qdroptol*R(i,i));
15
 
   Q(:,i)=q/(norm(q)+1e-300);
16
 
end
17
 
*/
18
 
 
19
 
static inline double sqr(double x)
20
 
{
21
 
   return x*x;
22
 
}
23
 
 
24
 
struct SparseEntry
25
 
{
26
 
   int index;
27
 
   double value;
28
 
 
29
 
   SparseEntry(void) {}
30
 
   SparseEntry(int index_, double value_) : index(index_), value(value_) {}
31
 
   bool operator<(const SparseEntry &e) const { return index<e.index; }
32
 
};
33
 
 
34
 
typedef std::list<SparseEntry> DynamicSparseVector;
35
 
typedef std::vector<SparseEntry> StaticSparseVector;
36
 
 
37
 
static void copy_column(int colstart, int nextcolstart, const int *rowindex, const double *value,
38
 
                        DynamicSparseVector &v)
39
 
{
40
 
   for(int i=colstart; i<nextcolstart; ++i)
41
 
      v.push_back(SparseEntry(rowindex[i], value[i]));
42
 
}
43
 
 
44
 
static double dot_product(const StaticSparseVector &u, const DynamicSparseVector &v)
45
 
{
46
 
   StaticSparseVector::const_iterator p=u.begin();
47
 
   if(p==u.end()) return 0;
48
 
   DynamicSparseVector::const_iterator q=v.begin();
49
 
   if(q==v.end()) return 0;
50
 
   double result=0;
51
 
   for(;;){
52
 
      if(p->index < q->index){
53
 
         ++p; if(p==u.end()) break;
54
 
      }else if(p->index > q->index){
55
 
         ++q; if(q==v.end()) break;
56
 
      }else{ // p->index == q->index
57
 
         result+=p->value*q->value;
58
 
         ++p; if(p==u.end()) break;
59
 
         ++q; if(q==v.end()) break;
60
 
      }
61
 
   }
62
 
   return result;
63
 
}
64
 
 
65
 
static void add_to_vector(double multiplier, const StaticSparseVector &u, DynamicSparseVector &v)
66
 
{
67
 
   StaticSparseVector::const_iterator p=u.begin();
68
 
   if(p==u.end()) return;
69
 
   DynamicSparseVector::iterator q=v.begin();
70
 
   for(;;){
71
 
      if(q==v.end()) break;
72
 
      if(p->index < q->index){
73
 
         q=v.insert(q, SparseEntry(p->index, multiplier*p->value));
74
 
         ++p; if(p==u.end()) return;
75
 
      }else if(p->index > q->index){
76
 
         ++q; if(q==v.end()) break;
77
 
      }else{ // p->index == q->index){
78
 
         q->value+=multiplier*p->value;
79
 
         ++p; if(p==u.end()) return;
80
 
         ++q; if(q==v.end()) break;
81
 
      }
82
 
   }
83
 
   for(;;){
84
 
      v.push_back(SparseEntry(p->index, multiplier*p->value));
85
 
      ++p; if(p==u.end()) return;
86
 
   }
87
 
}
88
 
 
89
 
static double vector_norm(const DynamicSparseVector &v)
90
 
{
91
 
   double norm_squared=0;
92
 
   DynamicSparseVector::const_iterator p;
93
 
   for(p=v.begin(); p!=v.end(); ++p) norm_squared+=sqr(p->value);
94
 
   return std::sqrt(norm_squared);
95
 
}
96
 
 
97
 
static void copy_large_entries(const DynamicSparseVector &u, double threshhold, StaticSparseVector &v)
98
 
{
99
 
   DynamicSparseVector::const_iterator p;
100
 
   for(p=u.begin(); p!=u.end(); ++p){
101
 
      if(std::fabs(p->value)>threshhold) v.push_back(*p);
102
 
   }
103
 
}
104
 
 
105
 
static void normalize(StaticSparseVector &v)
106
 
{
107
 
   double norm_squared=0;
108
 
   StaticSparseVector::iterator p;
109
 
   for(p=v.begin(); p!=v.end(); ++p) norm_squared+=sqr(p->value);
110
 
   if(norm_squared>0){
111
 
      double scale=1/std::sqrt(norm_squared);
112
 
      for(p=v.begin(); p!=v.end(); ++p) p->value*=scale;
113
 
   }
114
 
}
115
 
 
116
 
void incomplete_qr(int m, int n, const int *Acolstart, const int *Arowindex, const double *Avalue,
117
 
                   double Qdroptol,
118
 
                   const int *Rcolstart, const int *Rrowindex, double *Rvalue, double *Rdiag)
119
 
{
120
 
   // column storage for intermediate Q
121
 
   std::vector<StaticSparseVector> Q(n);
122
 
 
123
 
   // for each column of Q, how many more times it will be used in making R
124
 
   // (once this gets decremented to zero, we can free the column)
125
 
   std::vector<int> dependency_count(n,0);
126
 
   for(int i=0; i<n; ++i){
127
 
      for(int a=Rcolstart[i]; a<Rcolstart[i+1]; ++a){
128
 
         int j=Rrowindex[a];
129
 
         if(j<i) ++dependency_count[j];
130
 
      }
131
 
   }
132
 
 
133
 
   // construct R a column at a time
134
 
   for(int i=0; i<n; ++i){
135
 
      DynamicSparseVector q;
136
 
      copy_column(Acolstart[i], Acolstart[i+1], Arowindex, Avalue, q);
137
 
      for(int a=Rcolstart[i]; a<Rcolstart[i+1]; ++a){
138
 
         int j=Rrowindex[a];
139
 
         if(j<i){
140
 
            Rvalue[a]=dot_product(Q[j], q);
141
 
            add_to_vector(-Rvalue[a], Q[j], q);
142
 
            --dependency_count[j];
143
 
            if(dependency_count[j]==0) Q[j].clear();
144
 
         }else // j>i
145
 
            Rvalue[a]=0; // ideally structure of R is upper triangular, so this waste of space won't happen
146
 
      }
147
 
      Rdiag[i]=vector_norm(q);
148
 
      copy_large_entries(q, Qdroptol*Rdiag[i], Q[i]);
149
 
      normalize(Q[i]);
150
 
   }
151
 
}
152