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

Last change on this file since 813 was 813, checked in by ansari, 25 years ago

Correction bug/amelioarions TArray,TMatrix,TVector - Reza 5/4/2000
File size: 13.4 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 os << "\n--- " << InfoString() ;
311 os << " ND=" << ndim_ << " SizeX*Y*...= " ;
312 for(int k=0; k<ndim_; k++) {
313 os << size_[k];
314 if (k<ndim_-1) os << "x";
315 }
316 os << " ---" << endl;
317 if (prt_lev_ > 0) {
318 os << " TotSize= " << totsize_ << " Step(X Y ...)=" ;
319 for(int k=0; k<ndim_; k++) os << step_[k] << " " ;
320 os << " Offset= " << offset_ << endl;
321 }
322 if (prt_lev_ > 1) {
323 os << " MemoryMapping=" << GetMemoryMapping() << " VecType= " << GetVectorType()
324 << " RowsKA= " << RowsKA() << " ColsKA= " << ColsKA()
325 << " VectKA=" << VectKA() << endl;
326 }
327 if (!si && (prt_lev_ < 2)) return;
328 if (mInfo != NULL) os << (*mInfo) << endl;
329
330}
331
332string BaseArray::InfoString() const
333{
334 string rs = "BaseArray Type= ";
335 rs += typeid(*this).name() ;
336 return rs;
337}
338
339DVList& BaseArray::Info()
340{
341if (mInfo == NULL) mInfo = new DVList;
342return(*mInfo);
343}
344
345
346bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, uint_4 step, uint_8 offset, string & exmsg)
347{
348 if (ndim >= BASEARRAY_MAXNDIMS) {
349 exmsg += " NDim Error"; return false;
350 }
351 if (step < 1) {
352 exmsg += " Step(=0) Error"; return false;
353 }
354
355 minstep_ = moystep_ = step;
356
357 // Flagging bad updates ...
358 ndim_ = 0;
359
360 totsize_ = 1;
361 int k;
362 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
363 size_[k] = 1;
364 step_[k] = 0;
365 }
366 for(k=0; k<ndim; k++) {
367 size_[k] = siz[k] ;
368 step_[k] = totsize_*step;
369 totsize_ *= size_[k];
370 }
371 if (totsize_ < 1) {
372 exmsg += " Size Error"; return false;
373 }
374 offset_ = offset;
375 // Default for matrices : Memory organisation and Vector type
376 if (default_memory_mapping == CMemoryMapping) {
377 marowi_ = 1; macoli_ = 0;
378 }
379 else {
380 marowi_ = 0; macoli_ = 1;
381 }
382 veceli_ = (default_vector_type == ColumnVector ) ? marowi_ : macoli_;
383 ck_memo_vt_ = false; // Default : Don't Check MemMapping and VectorType for CompareSize
384 // Update OK
385 ndim_ = ndim;
386 return true;
387}
388
389bool BaseArray::UpdateSizes(uint_4 ndim, const uint_4 * siz, const uint_4 * step, uint_8 offset, string & exmsg)
390{
391 if (ndim >= BASEARRAY_MAXNDIMS) {
392 exmsg += " NDim Error"; return false;
393 }
394
395 // Flagging bad updates ...
396 ndim_ = 0;
397
398 totsize_ = 1;
399 int k;
400 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
401 size_[k] = 1;
402 step_[k] = 0;
403 }
404 uint_4 minstep = step[0];
405 for(k=0; k<ndim; k++) {
406 size_[k] = siz[k] ;
407 step_[k] = step[k];
408 totsize_ *= size_[k];
409 if (step_[k] < minstep) minstep = step_[k];
410 }
411 if (minstep < 1) {
412 exmsg += " Step(=0) Error"; return false;
413 }
414 if (totsize_ < 1) {
415 exmsg += " Size Error"; return false;
416 }
417 uint_8 plast = 0;
418 for(k=0; k<ndim; k++) plast += (siz[k]-1)*step[k];
419 if (plast == minstep*totsize_ ) moystep_ = minstep;
420 else moystep_ = 0;
421 minstep_ = minstep;
422 offset_ = offset;
423 // Default for matrices : Memory organisation and Vector type
424 if (default_memory_mapping == CMemoryMapping) {
425 marowi_ = 1; macoli_ = 0;
426 }
427 else {
428 marowi_ = 0; macoli_ = 1;
429 }
430 veceli_ = (default_vector_type == ColumnVector ) ? marowi_ : macoli_;
431 ck_memo_vt_ = false; // Default : Don't Check MemMapping and VectorType for CompareSize
432 // Update OK
433 ndim_ = ndim;
434 return true;
435}
436
437bool BaseArray::UpdateSizes(const BaseArray& a, string & exmsg)
438{
439 if (a.ndim_ >= BASEARRAY_MAXNDIMS) {
440 exmsg += " NDim Error"; return false;
441 }
442
443 // Flagging bad updates ...
444 ndim_ = 0;
445
446 totsize_ = 1;
447 int k;
448 for(k=0; k<BASEARRAY_MAXNDIMS; k++) {
449 size_[k] = 1;
450 step_[k] = 0;
451 }
452 uint_4 minstep = a.step_[0];
453 for(k=0; k<a.ndim_; k++) {
454 size_[k] = a.size_[k] ;
455 step_[k] = a.step_[k];
456 totsize_ *= size_[k];
457 if (step_[k] < minstep) minstep = step_[k];
458 }
459 if (minstep < 1) {
460 exmsg += " Step(=0) Error"; return false;
461 }
462 if (totsize_ < 1) {
463 exmsg += " Size Error"; return false;
464 }
465
466 minstep_ = a.minstep_;
467 moystep_ = a.moystep_;
468 offset_ = a.offset_;
469 macoli_ = a.macoli_;
470 marowi_ = a.marowi_;
471 veceli_ = a.veceli_;
472 ck_memo_vt_ = a.ck_memo_vt_;
473 // Update OK
474 ndim_ = a.ndim_;
475 return true;
476}
477
478
479void BaseArray::UpdateSubArraySizes(BaseArray & ra, uint_4 ndim, uint_4 * siz, uint_4 * pos, uint_4 * step) const
480{
481 if ( (ndim > ndim_) || (ndim < 1) )
482 throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) NDim Error") );
483 int k;
484 for(k=0; k<ndim; k++)
485 if ( (siz[k]*step[k]+pos[k]) > size_[k] )
486 throw(SzMismatchError("BaseArray::UpdateSubArraySizes( ... ) Size/Pos Error") );
487 uint_8 offset = offset_;
488 for(k=0; k<ndim_; k++) {
489 offset += pos[k]*step_[k];
490 step[k] *= step_[k];
491 }
492 string exm = "BaseArray::UpdateSubArraySizes() ";
493 if (!ra.UpdateSizes(ndim, siz, step, offset, exm))
494 throw( ParmError(exm) );
495 return;
496}
497
498
Note: See TracBrowser for help on using the repository browser.