Charmxi: Fix reversed arugments passed to generateSignature.
[charm.git] / src / langs / charj / src / charj / libs / Array.h
1 #ifndef ARRAY_H
2 #define ARRAY_H
3
4 #include <cstdarg>
5 #include <cstdio>
6 #include <cstdlib>
7 #include <cassert>
8 #include <iostream>
9 #include <algorithm>
10 #include <charm++.h>
11 #if CMK_HAS_CBLAS
12 #include <cblas.h>
13 #endif
14
15 namespace CharjArray {
16   class Range {
17   public:
18     int size, start, stop;
19     Range() {}
20     Range(int size_) : size(size_), start(0), stop(size) {}
21     Range(int start_, int stop_) :
22     size(stop_ - start_), start(start_), stop(stop_) {
23       assert(stop >= start);
24     }
25   };
26
27   template<int dims>
28   class Domain {
29   public:
30     Range ranges[dims];
31     
32     Domain() {}
33
34     Domain(Range ranges_[]) {
35       for (int i = 0; i < dims; i++) 
36         ranges[i] = ranges_[i];      
37     }
38
39     Domain(Range range) {
40       ranges[0] = range;
41     }
42
43     Domain(Range range1, Range range2) {
44       // TODO: fix Charj generator so it uses the array
45       ranges[0] = range1;
46       ranges[1] = range2;
47     }
48
49     int size() const {
50       int total = 0;
51
52       for (int i = 0; i < dims; i++)
53         if (total == 0)
54           total = ranges[i].size;
55         else
56           total *= ranges[i].size;
57
58       return total;
59     }
60   };
61
62   template<int dims>
63   class RowMajor {
64   public:
65     static int access(const int i, const Domain<dims> &d) {
66       return i - d.ranges[0].start;
67     }
68     static int access(const int i, const int j, const Domain<dims> &d) {
69       return (i - d.ranges[0].start) * d.ranges[1].size + j -
70         d.ranges[1].start;
71     }
72     // Generic access method, not used right now.
73     // static int access(const int *i, const Domain<dims> &d) {
74     //   int off = i[0];
75     //   int dimoff = 1;
76     //   for (int j = ndims-1; j > 0; --j) {
77     //     dimoff *= d.ranges[j].size;
78     //     off += dimoff * (i[j] - d.ranges[j].start);
79     //   }
80     //   return off;
81     // }
82   };
83
84   template<int dims>
85   class ColMajor {
86   public:
87     static int access(const int i, const Domain<dims> &d) {
88       return i - d.ranges[0].start;
89     }
90     static int access(const int i, const int j, const Domain<dims> &d) {
91       return (j - d.ranges[1].start) * d.ranges[1].size + i -
92         d.ranges[0].start;
93     }
94   };
95
96   template<class type, int dims = 1, class atype = RowMajor<dims> >
97   class Array {
98   private:
99     Domain<dims> domain;
100     type *block;
101     int ref_cout;
102     Array* ref_parent;
103
104   public:
105     Array(Domain<dims> domain_) : ref_parent(0) {
106       init(domain_);
107     }
108
109     Array(type **block_) {
110       block = *block_;
111     }
112
113     Array() : ref_parent(0) {
114
115     }
116
117     Array(Array* parent, Domain<dims> domain_) : ref_parent(parent) {
118       domain = domain_;
119       block = parent->block;
120     }
121
122     void init(Domain<dims> &domain_) {
123       domain = domain_;
124       //if (atype == ROW_MAJOR)
125       block = new type[domain.size()];
126     }
127
128     ~Array() {
129       delete[] block;
130     }
131
132     /*type* operator[] (const Domain<dims> &domain) {
133       return block[domain.ranges[0].size];
134       }*/
135
136     type& operator[] (const int i) {
137       return block[atype::access(i, domain)];
138     }
139
140     const type& operator[] (const int i) const {
141       return block[atype::access(i, domain)];
142     }
143
144     type& access(const int i, const int j) {
145       return block[atype::access(i, j, domain)];
146     }
147
148     const type& access(const int i, const int j) const {
149       return block[atype::access(i, j, domain)];
150     }
151
152     Array<type, dims, atype>* operator[] (const Domain<dims> &domain) {
153       return new Array<type, dims, atype>(this, domain);
154     }
155
156     int size() const {
157       return domain.size();
158     }
159
160     int size(int dim) const {
161       return domain.ranges[dim].size;
162     }
163
164     void pup(PUP::er& p) { }
165
166     void fill(const type &t) {
167       for (int i = 0; i < domain.size(); ++i)
168         block[i] = t;
169     }
170
171     /// Do these arrays have the same shape and contents?
172     bool operator==(const Array &rhs) const
173     {
174       for (int i = 0; i < dims; ++i)
175         if (this->size(i) != rhs.size(i))
176           return false;
177
178       for (int i = 0; i < this->size(); ++i)
179         if (this->block[i] != rhs.block[i])
180           return false;
181
182       return true;
183     }
184     bool operator!=(const Array &rhs) const
185     {
186       return !(*this == rhs);
187     }
188   };
189
190   /**
191      A local Matrix class for various sorts of linear-algebra work.
192
193      Indexed from 0, to reflect the C-heritage of Charj.
194    */
195   template <typename V, class atype = RowMajor<2> >
196   class Matrix : public Array<V, 2, atype>
197   {
198   public:
199     Matrix() { }
200     /// A square matrix
201     Matrix(unsigned int n) : Array<V,2,atype>(Domain<2>(n,n)) { }
202
203     /// A identity matrix
204     static Matrix* ident(int n)
205     {
206       Matrix *ret = new Matrix(n);
207       ret->fill(0);
208
209       for (int i = 0; i < n; ++i)
210         ret->access(i,i) = 1;
211
212       return ret;
213     }
214   };
215
216   template <typename T, class atype = RowMajor<1> >
217   class Vector : public Array<T, 1, atype>
218   {
219   public:
220     Vector() { }
221     Vector(unsigned int n) : Array<T, 1, atype>(Range(n)) { }
222   };
223
224   /// Compute the inner (dot) product v1^T * v2
225   // To compute v1^H * v2, call as dot(v1.C(), v2)
226   template<typename T, class atype1, class atype2>
227   T dot(const Vector<T, atype1> *pv1, const Vector<T, atype2> *pv2)
228   {
229     const Vector<T, atype1> &v1 = *pv1, &v2 = *pv2;
230     assert(v1.size() == v2.size());
231     // XXX: This default initialization worries me some, since it
232     // won't necessarily be an additive identity for all T. - Phil
233     T ret = T();
234     int size = v1.size();
235     for (int i = 0; i < size; ++i)
236       ret += v1[i] * v2[i];
237     return ret;
238   }
239 #if CMK_HAS_CBLAS
240   template <>
241   float dot<float, RowMajor<1>, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv1,
242                                               const Vector<float, RowMajor<1> > *pv2)
243   {
244     const Vector<float, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
245     assert(v1.size() == v2.size());
246     return cblas_sdot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
247   }
248   template <>
249   double dot<double, RowMajor<1>, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv1,
250                                                 const Vector<double, RowMajor<1> > *pv2)
251   {
252     const Vector<double, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
253     assert(v1.size() == v2.size());
254     return cblas_ddot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
255   }
256 #endif
257
258   /// Computer the 1-norm of the given vector
259   template<typename T, class atype>
260     T norm1(const Vector<T, atype> *pv)
261   {
262     const Vector<T, atype> &v = *pv;
263     // XXX: See comment about additive identity in dot(), above
264     T ret = T();
265     int size = v.size();
266     for (int i = 0; i < size; ++i)
267       ret += v[i];
268     return ret;
269   }
270
271   /// Compute the Euclidean (2) norm of the given vector
272   template<typename T, class atype>
273     T norm2(const Vector<T, atype> *pv)
274   {
275     const Vector<T, atype> &v = *pv;
276     // XXX: See comment about additive identity in dot(), above
277     T ret = T();
278     int size = v.size();
279     for (int i = 0; i < size; ++i)
280       ret += v[i] * v[i];
281     return sqrt(ret);
282   }
283 #if CMK_HAS_CBLAS
284   template<>
285     float norm2<float, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv)
286   {
287     const Vector<float, RowMajor<1> > &v = *pv;
288     return cblas_snrm2(v.size(), &(v[0]), 1);
289   }
290   template<>
291     double norm2<double, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv)
292   {
293     const Vector<double, RowMajor<1> > &v = *pv;
294     return cblas_dnrm2(v.size(), &(v[0]), 1);
295   }
296 #endif
297
298   /// Compute the infinity (max) norm of the given vector
299   // Will fail on zero-length vectors
300   template<typename T, class atype>
301     T normI(const Vector<T, atype> *pv)
302   {
303     const Vector<T, atype> &v = *pv;
304     T ret = v[0];
305     int size = v.size();
306     for (int i = 1; i < size; ++i)
307       ret = max(ret, v[i]);
308     return ret;
309   }
310
311   /// Scale a vector by some constant
312   template<typename T, typename U, class atype>
313     void scale(const T &t, Vector<U, atype> *pv)
314   {
315     const Vector<T, atype> &v = *pv;
316     int size = v.size();
317     for (int i = 0; i < size; ++i)
318       v[i] = t * v[i];
319   }
320 #if CMK_HAS_CBLAS
321   template<>
322     void scale<float, float, RowMajor<1> >(const float &t,
323                                            Vector<float, RowMajor<1> > *pv)
324   {
325     Vector<float, RowMajor<1> > &v = *pv;
326     cblas_sscal(v.size(), t, &(v[0]), 1);
327   }
328   template<>
329     void scale<double, double, RowMajor<1> >(const double &t,
330                                              Vector<double, RowMajor<1> > *pv)
331   {
332     Vector<double, RowMajor<1> > &v = *pv;
333     cblas_dscal(v.size(), t, &(v[0]), 1);
334   }
335 #endif
336
337   /// Add one vector to a scaled version of another
338   template<typename T, typename U, class atype>
339     void axpy(const T &a, const Vector<U, atype> *px, Vector<U, atype> *py)
340   {
341     Vector<T, atype> &x = *px;
342     const Vector<T, atype> &y = *py;
343     int size = x.size();
344     assert(size == y.size());
345     for (int i = 0; i < size; ++i)
346       x[i] = a * x[i] + y[i];
347   }
348 #if CMK_HAS_CBLAS
349   template<>
350     void axpy<float, float, RowMajor<1> >(const float &a,
351                                           const Vector<float, RowMajor<1> > *px,
352                                           Vector<float, RowMajor<1> > *py)
353   {
354     const Vector<float, RowMajor<1> > &x = *px;
355     Vector<float, RowMajor<1> > &y = *py;
356     int size = x.size();
357     assert(size == y.size());
358     cblas_saxpy(size, a, &(x[0]), 1, &(y[0]), 1);
359   }
360   template<>
361     void axpy<double, double, RowMajor<1> >(const double &a,
362                                           const Vector<double, RowMajor<1> > *px,
363                                           Vector<double, RowMajor<1> > *py)
364   {
365     const Vector<double, RowMajor<1> > &x = *px;
366     Vector<double, RowMajor<1> > &y = *py;
367     int size = x.size();
368     assert(size == y.size());
369     cblas_daxpy(size, a, &(x[0]), 1, &(y[0]), 1);
370   }
371 #endif
372 }
373
374 #endif