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

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

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