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