source: Sophya/trunk/SophyaLib/TArray/basarr.cc@ 860

Last change on this file since 860 was 850, checked in by ansari, 26 years ago

Corrections divers + RansomSequence - Reza 10/4/2000
File size: 13.5 KB
Line 
1// Base class for numerical arrays
2// R. Ansari, C.Magneville 03/2000
3
4#include "machdefs.h"
5#include <stdio.h>
6#include <stdlib.h>
7#include "pexceptions.h"
8#include "basarr.h"
9
10
11
12// Variables statiques globales
13char * BaseArray::ck_op_msg_[6] = {"???", "Size(int )", "IsPacked(int )",
14 "Stride(int )", "ElemCheckBound()", "operator()" };
15uint_4 BaseArray::max_nprt_ = 50;
16uint_4 BaseArray::prt_lev_ = 0;
17short BaseArray::default_memory_mapping = CMemoryMapping;
18short BaseArray::default_vector_type = ColumnVector;
19uint_8 BaseArray::openmp_size_threshold = 200000;
20
21// ------ Methodes statiques globales --------
22
23// 1/ Gestion des impressions
24void BaseArray::SetMaxPrint(uint_4 nprt, uint_4 lev)
25{
26 max_nprt_ = nprt;
27 prt_lev_ = (lev < 3) ? lev : 3;
28}
29
30void BaseArray::SetOpenMPSizeThreshold(uint_8 thr)
31{
32 openmp_size_threshold = thr;
33}
34
35
36uint_8 BaseArray::ComputeTotalSize(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset)
37{
38 uint_8 rs = step;
39 for(int k=0; k<ndim; k++) rs *= siz[k];
40 return(rs+offset);
41}
42
43short BaseArray::SetDefaultMemoryMapping(short mm)
44{
45 default_memory_mapping = (mm != CMemoryMapping) ? FortranMemoryMapping : CMemoryMapping;
46 return default_memory_mapping;
47}
48
49short BaseArray::SetDefaultVectorType(short vt)
50{
51 default_vector_type = (vt != ColumnVector) ? RowVector : ColumnVector ;
52 return default_vector_type;
53}
54
55// Memory organisation
56short BaseArray::SelectMemoryMapping(short mm)
57// si mm == CMemoryMapping, elements d'une ligne suivant X (consecutif)
58// sinon (mm == FortranMemoryMapping) elements d'une colonne suivant X
59{
60 if ( (mm == CMemoryMapping) || (mm == FortranMemoryMapping) ) return (mm) ;
61 else return (default_memory_mapping);
62}
63short BaseArray::SelectVectorType(short vt)
64{
65 if ((vt == ColumnVector) || (vt == RowVector)) return(vt);
66 else return(default_vector_type);
67}
68
69void BaseArray::UpdateMemoryMapping(short mm)
70{
71 short vt = default_vector_type;
72 if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
73 if (mm == CMemoryMapping) {
74 marowi_ = 1; macoli_ = 0;
75 }
76 else {
77 marowi_ = 0; macoli_ = 1;
78 }
79
80 if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
81 // Choix automatique Vecteur ligne ou colonne
82 if ( size_[macoli_] == 1) veceli_ = marowi_;
83 else veceli_ = macoli_;
84 }
85 else veceli_ = (vt == ColumnVector ) ? marowi_ : macoli_;
86 ck_memo_vt_ = true; // Check MemMapping and VectorType for CompareSize
87}
88
89void BaseArray::UpdateMemoryMapping(BaseArray const & a, short mm)
90{
91 short vt = default_vector_type;
92 if (mm == SameMemoryMapping) {
93 mm = ((a.marowi_ == 1) ? CMemoryMapping : FortranMemoryMapping);
94 vt = (a.marowi_ == a.veceli_) ? ColumnVector : RowVector;
95 }
96 else if (mm == AutoMemoryMapping) mm = default_memory_mapping;
97
98 if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) mm = default_memory_mapping;
99 if (mm == CMemoryMapping) {
100 marowi_ = 1; macoli_ = 0;
101 }
102 else {
103 marowi_ = 0; macoli_ = 1;
104 }
105 if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
106 // Choix automatique Vecteur ligne ou colonne
107 if ( size_[macoli_] == 1) veceli_ = marowi_;
108 else veceli_ = marowi_;
109 }
110 else veceli_ = (vt == ColumnVector ) ? marowi_ : macoli_;
111 ck_memo_vt_ = true; // Check MemMapping and VectorType for CompareSize
112}
113
114void BaseArray::SetMemoryMapping(short mm)
115{
116 if (mm == SameMemoryMapping) return;
117 if (mm == AutoMemoryMapping) mm = default_memory_mapping;
118 if ( (mm != CMemoryMapping) && (mm != FortranMemoryMapping) ) return;
119 short vt = (marowi_ == veceli_) ? ColumnVector : RowVector;
120 if (mm == CMemoryMapping) {
121 marowi_ = 1; macoli_ = 0;
122 }
123 else {
124 marowi_ = 0; macoli_ = 1;
125 }
126 if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
127 // Choix automatique Vecteur ligne ou colonne
128 if ( size_[macoli_] == 1) veceli_ = marowi_;
129 else veceli_ = macoli_;
130 }
131 else veceli_ = (vt == ColumnVector ) ? marowi_ : macoli_;
132 ck_memo_vt_ = true; // Check MemMapping and VectorType for CompareSize
133}
134
135void BaseArray::SetVectorType(short vt)
136{
137 if (vt == SameVectorType) return;
138 if (vt == AutoVectorType) vt = default_vector_type;
139 if ( (ndim_ == 2) && ((size_[0] == 1) || (size_[1] == 1)) ) {
140 // Choix automatique Vecteur ligne ou colonne
141 if ( size_[macoli_] == 1) veceli_ = marowi_;
142 else veceli_ = macoli_;
143 }
144 else veceli_ = (vt == ColumnVector ) ? marowi_ : macoli_;
145 ck_memo_vt_ = true; // Check MemMapping and VectorType for CompareSize
146}
147
148// -------------------------------------------------------
149// Methodes de la classe
150// -------------------------------------------------------
151
152// Les constructeurs
153BaseArray::BaseArray()
154 : mInfo(NULL)
155// Default constructor
156{
157 ndim_ = 0;
158 for(int k=0; k<BASEARRAY_MAXNDIMS; k++) step_[k] = size_[k] = 0;
159 totsize_ = 0;
160 minstep_ = 0;
161 moystep_ = 0;
162 offset_ = 0;
163 // Default for matrices : Memory organisation and Vector type
164 if (default_memory_mapping == CMemoryMapping) {
165 marowi_ = 1; macoli_ = 0;
166 }
167 else {
168 marowi_ = 0; macoli_ = 1;
169 }
170 veceli_ = (default_vector_type == ColumnVector ) ? marowi_ : macoli_;
171 ck_memo_vt_ = false; // Default : Don't Check MemMapping and VectorType for CompareSize
172}
173
174// Destructeur
175BaseArray::~BaseArray()
176{
177}
178
179
180bool BaseArray::CompareSizes(const BaseArray& a)
181{
182 if (ndim_ != a.ndim_) return(false);
183 for(int k=0; k<ndim_; k++)
184 if (size_[k] != a.size_[k]) return(false);
185 // $CHECK$ Reza doit-on verifier ca
186 if (ck_memo_vt_ && a.ck_memo_vt_)
187 if ( (macoli_ != a.macoli_) || (marowi_ != a.marowi_) ||
188 (veceli_ != a.veceli_) ) return(false);
189 return(true);
190}
191
192void BaseArray::CompactAllDim()
193{
194 if (ndim_ < 2) return;
195 uint_4 ndim = 0;
196 uint_4 size[BASEARRAY_MAXNDIMS];
197 uint_4 step[BASEARRAY_MAXNDIMS];
198 for(int k=0; k<ndim_; k++) {
199 if (size_[k] < 2) continue;
200 size[ndim] = size_[k];
201 step[ndim] = step_[k];
202 ndim++;
203 }
204 if (ndim == 0) {
205 size[0] = size_[0];
206 step[0] = step_[0];
207 ndim = 1;
208 }
209 string exmsg = "BaseArray::CompactAllDim() ";
210 if (!UpdateSizes(ndim, size, step, offset_, exmsg)) throw( ParmError(exmsg) );
211 return;
212}
213
214void BaseArray::CompactTrailingDim()
215{
216 if (ndim_ < 2) return;
217 uint_4 ndim = 0;
218 uint_4 size[BASEARRAY_MAXNDIMS];
219 uint_4 step[BASEARRAY_MAXNDIMS];
220 for(int k=0; k<ndim_; k++) {
221 size[ndim] = size_[k];
222 step[ndim] = step_[k];
223 if (size_[k] > 1) ndim=k;
224 }
225 if (ndim == 0) ndim = 1;
226 string exmsg = "BaseArray::CompactTrailingDim() ";
227 if (!UpdateSizes(ndim, size, step, offset_, exmsg)) throw( ParmError(exmsg) );
228 return;
229}
230
231
232uint_4 BaseArray::MinStepKA() const
233{
234 for(int ka=0; ka<ndim_; ka++)
235 if (step_[ka] == minstep_) return(ka);
236 return(0);
237}
238
239uint_4 BaseArray::MaxSizeKA() const
240{
241 uint_4 ka = 0;
242 uint_4 mx = size_[0];
243 for(int k=0; k<ndim_; k++)
244 if (size_[k] > mx) { ka = k; mx = size_[k]; }
245 return(ka);
246}
247
248
249// Acces lineaire aux elements .... Calcul d'offset
250// --------------------------------------------------
251// Position de l'element 0 du vecteur i selon l'axe ka
252// --------------------------------------------------
253uint_8 BaseArray::Offset(uint_4 ka, uint_8 i) const
254{
255
256 if ( (ndim_ < 1) || (i == 0) ) return(offset_);
257 //#ifdef SO_BOUNDCHECKING
258 if (ka >= ndim_)
259 throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Axe KA Error");
260 if ( i*size_[ka] >= totsize_ )
261 throw RangeCheckError("BaseArray::Offset(uint_4 ka, uint_8 i) Index Error");
262 //#endif
263 uint_4 idx[BASEARRAY_MAXNDIMS];
264 int k;
265 uint_8 rest = i;
266 idx[ka] = 0;
267 for(k=0; k<ndim_; k++) {
268 if (k == ka) continue;
269 idx[k] = rest%size_[k]; rest /= size_[k];
270 }
271 uint_8 off = offset_;
272 for(k=0; k<ndim_; k++) off += idx[k]*step_[k];
273 return (off);
274}
275
276uint_8 BaseArray::Offset(uint_8 ip) const
277{
278 if ( (ndim_ < 1) || (ip == 0) ) return(offset_);
279 //#ifdef SO_BOUNDCHECKING
280 if (ip >= totsize_)
281 throw RangeCheckError("BaseArray::Offset(uint_8 ip) Out of range index ip");
282 //#endif
283
284 uint_4 idx[BASEARRAY_MAXNDIMS];
285 int k;
286 uint_8 rest = ip;
287 for(k=0; k<ndim_; k++) {
288 idx[k] = rest%size_[k]; rest /= size_[k];
289 }
290 //#ifdef SO_BOUNDCHECKING
291 if (rest != 0)
292 throw PError("BaseArray::Offset(uint_8 ip) GUG !!! rest != 0");
293 //#endif
294// if (rest != 0) cerr << " BUG ---- BaseArray::Offset( " << ip << " )" << rest << endl;
295// cerr << " DBG-Offset( " << ip << ")" ;
296// for(k=0; k<ndim_; k++) cerr << idx[k] << "," ;
297// cerr << " ZZZZ " << endl;
298 uint_8 off = offset_;
299 for(k=0; k<ndim_; k++) off += idx[k]*step_[k];
300 return (off);
301}
302
303
304// ----------------------------------------------------
305// Impression, etc ...
306// ----------------------------------------------------
307
308void BaseArray::Show(ostream& os, bool si) const
309{
310 if (ndim_ < 1) {
311 os << "\n--- " << BaseArray::InfoString() << " Unallocated Array ! " << endl;
312 return;
313 }
314 os << "\n--- " << InfoString() ;
315 os << " ND=" << ndim_ << " SizeX*Y*...= " ;
316 for(int k=0; k<ndim_; k++) {
317 os << size_[k];
318 if (k<ndim_-1) os << "x";
319 }
320 os << " ---" << endl;
321 if (prt_lev_ > 0) {
322 os << " TotSize= " << totsize_ << " Step(X Y ...)=" ;
323 for(int k=0; k<ndim_; k++) os << step_[k] << " " ;
324 os << " Offset= " << offset_ << endl;
325 }
326 if (prt_lev_ > 1) {
327 os << " MemoryMapping=" << GetMemoryMapping() << " VecType= " << GetVectorType()
328 << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA()
329 << " VectKA=" << VectKA() << endl;
330 }
331 if (!si && (prt_lev_ < 2)) return;
332 if (mInfo != NULL) os << (*mInfo) << endl;
333
334}
335
336string BaseArray::InfoString() const
337{
338 string rs = "BaseArray Type= ";
339 rs += typeid(*this).name() ;
340 return rs;
341}
342
343DVList& BaseArray::Info()
344{
345if (mInfo == NULL) mInfo = new DVList;
346return(*mInfo);
347}
348
349
350bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
351{
352 if (ndim >= BASEARRAY_MAXNDIMS) {
353 exmsg += " NDim Error"; return false;
354 }
355 if (step < 1) {
356 exmsg += " Step(=0) Error"; return false;
357 }
358
359 minstep_ = moystep_ = step;
360
361 // Flagging bad updates ...
362 ndim_ = 0;
363
364 totsize_ = 1;
365 int k;
366 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
367 size_[k] = 1;
368 step_[k] = 0;
369 }
370 for(k=0; k<ndim; k++) {
371 size_[k] = siz[k] ;
372 step_[k] = totsize_*step;
373 totsize_ *= size_[k];
374 }
375 if (totsize_ < 1) {
376 exmsg += " Size Error"; return false;
377 }
378 offset_ = offset;
379 // Default for matrices : Memory organisation and Vector type
380 if (default_memory_mapping == CMemoryMapping) {
381 marowi_ = 1; macoli_ = 0;
382 }
383 else {
384 marowi_ = 0; macoli_ = 1;
385 }
386 veceli_ = (default_vector_type == ColumnVector ) ? marowi_ : macoli_;
387 ck_memo_vt_ = false; // Default : Don't Check MemMapping and VectorType for CompareSize
388 // Update OK
389 ndim_ = ndim;
390 return true;
391}
392
393bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
394{
395 if (ndim >= BASEARRAY_MAXNDIMS) {
396 exmsg += " NDim Error"; return false;
397 }
398
399 // Flagging bad updates ...
400 ndim_ = 0;
401
402 totsize_ = 1;
403 int k;
404 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
405 size_[k] = 1;
406 step_[k] = 0;
407 }
408 uint_4 minstep = step[0];
409 for(k=0; k<ndim; k++) {
410 size_[k] = siz[k] ;
411 step_[k] = step[k];
412 totsize_ *= size_[k];
413 if (step_[k] < minstep) minstep = step_[k];
414 }
415 if (minstep < 1) {
416 exmsg += " Step(=0) Error"; return false;
417 }
418 if (totsize_ < 1) {
419 exmsg += " Size Error"; return false;
420 }
421 uint_8 plast = 0;
422 for(k=0; k<ndim; k++) plast += (siz[k]-1)*step[k];
423 if (plast == minstep*totsize_ ) moystep_ = minstep;
424 else moystep_ = 0;
425 minstep_ = minstep;
426 offset_ = offset;
427 // Default for matrices : Memory organisation and Vector type
428 if (default_memory_mapping == CMemoryMapping) {
429 marowi_ = 1; macoli_ = 0;
430 }
431 else {
432 marowi_ = 0; macoli_ = 1;
433 }
434 veceli_ = (default_vector_type == ColumnVector ) ? marowi_ : macoli_;
435 ck_memo_vt_ = false; // Default : Don't Check MemMapping and VectorType for CompareSize
436 // Update OK
437 ndim_ = ndim;
438 return true;
439}
440
441bool BaseArray::UpdateSizes(const BaseArray& a, string & exmsg)
442{
443 if (a.ndim_ >= BASEARRAY_MAXNDIMS) {
444 exmsg += " NDim Error"; return false;
445 }
446
447 // Flagging bad updates ...
448 ndim_ = 0;
449
450 totsize_ = 1;
451 int k;
452 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
453 size_[k] = 1;
454 step_[k] = 0;
455 }
456 uint_4 minstep = a.step_[0];
457 for(k=0; k<a.ndim_; k++) {
458 size_[k] = a.size_[k] ;
459 step_[k] = a.step_[k];
460 totsize_ *= size_[k];
461 if (step_[k] < minstep) minstep = step_[k];
462 }
463 if (minstep < 1) {
464 exmsg += " Step(=0) Error"; return false;
465 }
466 if (totsize_ < 1) {
467 exmsg += " Size Error"; return false;
468 }
469
470 minstep_ = a.minstep_;
471 moystep_ = a.moystep_;
472 offset_ = a.offset_;
473 macoli_ = a.macoli_;
474 marowi_ = a.marowi_;
475 veceli_ = a.veceli_;
476 ck_memo_vt_ = a.ck_memo_vt_;
477 // Update OK
478 ndim_ = a.ndim_;
479 return true;
480}
481
482
483void BaseArray::UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const
484{
485 if ( (ndim > ndim_) || (ndim < 1) )
486 throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) NDim Error") );
487 int k;
488 for(k=0; k<ndim; k++)
489 if ( (siz[k]*step[k]+pos[k]) > size_[k] )
490 throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) Size/Pos Error") );
491 uint_8 offset = offset_;
492 for(k=0; k<ndim_; k++) {
493 offset += pos[k]*step_[k];
494 step[k] *= step_[k];
495 }
496 string exm = "BaseArray::UpdateSubArraySizes() ";
497 if (!ra.UpdateSizes(ndim, siz, step, offset, exm))
498 throw( ParmError(exm) );
499 return;
500}
501
502
Note: See TracBrowser for help on using the repository browser.