18a91c3b5a49639e562744ce501cae5269301913
[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(type& block_) {
114       block = &block_;
115     }
116
117     Array() : ref_parent(0) {
118
119     }
120
121     Array(Array* parent, Domain<dims> domain_) : ref_parent(parent) {
122       domain = domain_;
123       block = parent->block;
124     }
125
126     void init(Domain<dims> &domain_) {
127       domain = domain_;
128       //if (atype == ROW_MAJOR)
129       block = new type[domain.size()];
130       printf("Array: allocating memory, size=%d, base pointer=%p\n",
131              domain.size(), block);
132     }
133
134     type* raw() { return block; }
135
136     ~Array() {
137       delete[] block;
138     }
139
140     /*type* operator[] (const Domain<dims> &domain) {
141       return block[domain.ranges[0].size];
142       }*/
143
144     type& operator[] (const int i) {
145       return block[atype::access(i, domain)];
146     }
147
148     const type& operator[] (const int i) const {
149       return block[atype::access(i, domain)];
150     }
151
152     type& access(const int i) {
153       return this->operator[](i);
154     }
155
156     type& access(const int i, const int j) {
157       //printf("Array: accessing, index (%d,%d), offset=%d, base pointer=%p\n",
158       //i, j, atype::access(i, j, domain), block);
159       return block[atype::access(i, j, domain)];
160     }
161
162     type& access(const int i, const Range r) {
163       Domain<1> d(r);
164       //printf("Array: accessing subrange, size = %d, range (%d,%d), base pointer=%p\n",
165       //d.size(), r.start, r.stop, block);
166       type* buf = new type[d.size()];
167       for (int j = 0; j < d.size(); j++) {
168         //printf("Array: copying element (%d,%d), base pointer=%p\n", i, j, block);
169         buf[j] = block[atype::access(i, j, domain)];
170       }
171       return *buf;
172     }
173
174     const type& access(const int i, const int j) const {
175       return block[atype::access(i, j, domain)];
176     }
177
178     Array<type, dims, atype>* operator[] (const Domain<dims> &domain) {
179       return new Array<type, dims, atype>(this, domain);
180     }
181
182     int size() const {
183       return domain.size();
184     }
185
186     int size(int dim) const {
187       return domain.ranges[dim].size;
188     }
189
190     void pup(PUP::er& p) { }
191
192     void fill(const type &t) {
193       for (int i = 0; i < domain.size(); ++i)
194         block[i] = t;
195     }
196
197     /// Do these arrays have the same shape and contents?
198     bool operator==(const Array &rhs) const
199     {
200       for (int i = 0; i < dims; ++i)
201         if (this->size(i) != rhs.size(i))
202           return false;
203
204       for (int i = 0; i < this->size(); ++i)
205         if (this->block[i] != rhs.block[i])
206           return false;
207
208       return true;
209     }
210     bool operator!=(const Array &rhs) const
211     {
212       return !(*this == rhs);
213     }
214   };
215
216   /**
217      A local Matrix class for various sorts of linear-algebra work.
218
219      Indexed from 0, to reflect the C-heritage of Charj.
220    */
221   template <typename V, class atype = RowMajor<2> >
222   class Matrix : public Array<V, 2, atype>
223   {
224   public:
225     Matrix() { }
226     /// A square matrix
227     Matrix(unsigned int n) : Array<V,2,atype>(Domain<2>(n,n)) { }
228
229     /// A identity matrix
230     static Matrix* ident(int n)
231     {
232       Matrix *ret = new Matrix(n);
233       ret->fill(0);
234
235       for (int i = 0; i < n; ++i)
236         ret->access(i,i) = 1;
237
238       return ret;
239     }
240   };
241
242   template <typename T, class atype = RowMajor<1> >
243   class Vector : public Array<T, 1, atype>
244   {
245   public:
246     Vector() { }
247     Vector(unsigned int n) : Array<T, 1, atype>(Range(n)) { }
248   };
249
250   /// Compute the inner (dot) product v1^T * v2
251   // To compute v1^H * v2, call as dot(v1.C(), v2)
252   template<typename T, class atype1, class atype2>
253   T dot(const Vector<T, atype1> *pv1, const Vector<T, atype2> *pv2)
254   {
255     const Vector<T, atype1> &v1 = *pv1, &v2 = *pv2;
256     assert(v1.size() == v2.size());
257     // XXX: This default initialization worries me some, since it
258     // won't necessarily be an additive identity for all T. - Phil
259     T ret = T();
260     int size = v1.size();
261     for (int i = 0; i < size; ++i)
262       ret += v1[i] * v2[i];
263     return ret;
264   }
265 #if CMK_HAS_CBLAS
266   template <>
267   float dot<float, RowMajor<1>, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv1,
268                                               const Vector<float, RowMajor<1> > *pv2)
269   {
270     const Vector<float, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
271     assert(v1.size() == v2.size());
272     return cblas_sdot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
273   }
274   template <>
275   double dot<double, RowMajor<1>, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv1,
276                                                 const Vector<double, RowMajor<1> > *pv2)
277   {
278     const Vector<double, RowMajor<1> > &v1 = *pv1, &v2 = *pv2;
279     assert(v1.size() == v2.size());
280     return cblas_ddot(v1.size(), &(v1[0]), 1, &(v2[0]), 1);
281   }
282 #endif
283
284   /// Computer the 1-norm of the given vector
285   template<typename T, class atype>
286     T norm1(const Vector<T, atype> *pv)
287   {
288     const Vector<T, atype> &v = *pv;
289     // XXX: See comment about additive identity in dot(), above
290     T ret = T();
291     int size = v.size();
292     for (int i = 0; i < size; ++i)
293       ret += v[i];
294     return ret;
295   }
296
297   /// Compute the Euclidean (2) norm of the given vector
298   template<typename T, class atype>
299     T norm2(const Vector<T, atype> *pv)
300   {
301     const Vector<T, atype> &v = *pv;
302     // XXX: See comment about additive identity in dot(), above
303     T ret = T();
304     int size = v.size();
305     for (int i = 0; i < size; ++i)
306       ret += v[i] * v[i];
307     return sqrt(ret);
308   }
309 #if CMK_HAS_CBLAS
310   template<>
311     float norm2<float, RowMajor<1> >(const Vector<float, RowMajor<1> > *pv)
312   {
313     const Vector<float, RowMajor<1> > &v = *pv;
314     return cblas_snrm2(v.size(), &(v[0]), 1);
315   }
316   template<>
317     double norm2<double, RowMajor<1> >(const Vector<double, RowMajor<1> > *pv)
318   {
319     const Vector<double, RowMajor<1> > &v = *pv;
320     return cblas_dnrm2(v.size(), &(v[0]), 1);
321   }
322 #endif
323
324   /// Compute the infinity (max) norm of the given vector
325   // Will fail on zero-length vectors
326   template<typename T, class atype>
327     T normI(const Vector<T, atype> *pv)
328   {
329     const Vector<T, atype> &v = *pv;
330     T ret = v[0];
331     int size = v.size();
332     for (int i = 1; i < size; ++i)
333       ret = max(ret, v[i]);
334     return ret;
335   }
336
337   /// Scale a vector by some constant
338   template<typename T, typename U, class atype>
339     void scale(const T &t, Vector<U, atype> *pv)
340   {
341     const Vector<T, atype> &v = *pv;
342     int size = v.size();
343     for (int i = 0; i < size; ++i)
344       v[i] = t * v[i];
345   }
346 #if CMK_HAS_CBLAS
347   template<>
348     void scale<float, float, RowMajor<1> >(const float &t,
349                                            Vector<float, RowMajor<1> > *pv)
350   {
351     Vector<float, RowMajor<1> > &v = *pv;
352     cblas_sscal(v.size(), t, &(v[0]), 1);
353   }
354   template<>
355     void scale<double, double, RowMajor<1> >(const double &t,
356                                              Vector<double, RowMajor<1> > *pv)
357   {
358     Vector<double, RowMajor<1> > &v = *pv;
359     cblas_dscal(v.size(), t, &(v[0]), 1);
360   }
361 #endif
362
363   /// Add one vector to a scaled version of another
364   template<typename T, typename U, class atype>
365     void axpy(const T &a, const Vector<U, atype> *px, Vector<U, atype> *py)
366   {
367     Vector<T, atype> &x = *px;
368     const Vector<T, atype> &y = *py;
369     int size = x.size();
370     assert(size == y.size());
371     for (int i = 0; i < size; ++i)
372       x[i] = a * x[i] + y[i];
373   }
374 #if CMK_HAS_CBLAS
375   template<>
376     void axpy<float, float, RowMajor<1> >(const float &a,
377                                           const Vector<float, RowMajor<1> > *px,
378                                           Vector<float, RowMajor<1> > *py)
379   {
380     const Vector<float, RowMajor<1> > &x = *px;
381     Vector<float, RowMajor<1> > &y = *py;
382     int size = x.size();
383     assert(size == y.size());
384     cblas_saxpy(size, a, &(x[0]), 1, &(y[0]), 1);
385   }
386   template<>
387     void axpy<double, double, RowMajor<1> >(const double &a,
388                                           const Vector<double, RowMajor<1> > *px,
389                                           Vector<double, RowMajor<1> > *py)
390   {
391     const Vector<double, RowMajor<1> > &x = *px;
392     Vector<double, RowMajor<1> > &y = *py;
393     int size = x.size();
394     assert(size == y.size());
395     cblas_daxpy(size, a, &(x[0]), 1, &(y[0]), 1);
396   }
397 #endif
398 }
399
400 #endif