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