C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
ivector.cpp
1 /*
2 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3 **
4 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5 ** Universitaet Karlsruhe, Germany
6 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7 ** Universitaet Wuppertal, Germany
8 **
9 ** This library is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU Library General Public
11 ** License as published by the Free Software Foundation; either
12 ** version 2 of the License, or (at your option) any later version.
13 **
14 ** This library is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 ** Library General Public License for more details.
18 **
19 ** You should have received a copy of the GNU Library General Public
20 ** License along with this library; if not, write to the Free
21 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 */
23 
24 /* CVS $Id: ivector.cpp,v 1.26 2014/01/30 17:23:45 cxsc Exp $ */
25 
26 #define _CXSC_CPP
27 
28 #include "ivector.hpp"
29 #include "vector.inl"
30 #include "ivector.inl"
31 
32 #include "idotk.inl"
33 
34 namespace cxsc {
35 
36 int in ( const ivector& x, const ivector& y ) // Inner inclusion for two ivectors
37 { //---------------------------------
38  int i = Lb(x), n = Ub(x), ok = 1;
39 
40  while (ok && i <= n) { ok = in(x[i],y[i]); i++; }
41  return ok;
42 }
43 
44 int in ( int x, ivector& y ) // Inner inclusion of an integer in a ivector
45 { //-------------------------------------------
46  int i = Lb(y), n = Ub(y), ok = 1;
47  double d = x;
48 
49  while (ok && i <= n) { ok = in(d,y[i]); i++; }
50  return ok;
51 }
52 
53 ivector Blow ( const ivector& x, real eps ) // Epsilon deflation
54 { //------------------
55  int i;
56  ivector h(Lb(x),Ub(x));
57 
58  for (i = Lb(x); i <= Ub(x); i++) h[i] = Blow(x[i],eps);
59  return h;
60 }
61 
62 int Disjoint ( ivector& a, ivector& b ) // Test for disjointedness
63 { //------------------------
64  int al = Lb(a), au = Ub(a), bl = Lb(b);
65  int disjointed, i, d;
66 
67  d = bl - al;
68 
69  i = al;
70  disjointed = 0;
71  do {
72  if (Disjoint(a[i],b[i+d]))
73  disjointed = 1;
74  else
75  i++;
76  } while ( !(disjointed || i > au) );
77  return disjointed;
78 }
79 
80 int Zero ( ivector& x ) // Check for zero vector
81 { //----------------------
82  int i, ok;
83  for (i = Lb(x), ok = 1; i <= Ub(x) && ok; i++) ok = (x[i] == 0.0);
84  return ok;
85 }
86 
87 rvector mid ( ivector& v ) // Vector of midpoints
88 { //--------------------
89  int i;
90  int l = Lb(v), u = Ub(v);
91  rvector w(l,u);
92 
93  for (i = l; i <= u; i++) w[i] = mid(v[i]);
94  return w;
95 }
96 
97 real MaxRelDiam ( const ivector& v ) // Maximum relative diameter
98 { //--------------------------
99  real r;
100  int i, l=Lb(v), u=Ub(v);
101 
102  r = 0.0;
103  for (i = l; i <= u; i++)
104  if (RelDiam(v[i]) > r) r = RelDiam(v[i]);
105  return r;
106 }
107 
108 real MaxRelDiam ( const ivector_slice& v ) // Maximum relative diameter
109 { //--------------------------
110  real r;
111  int i, l=Lb(v), u=Ub(v);
112 
113  r = 0.0;
114  for (i = l; i <= u; i++)
115  if (RelDiam(v[i]) > r) r = RelDiam(v[i]);
116  return r;
117 }
118 
119 //----------------------------------------------------------------------------
120 // Checks if the diameter of the interval vector 'v' is less or equal to 'n'
121 // ulps. Ulp is an abbreviation for: units in the last place of the mantissa.
122 //----------------------------------------------------------------------------
123 int UlpAcc ( ivector& v, int n )
124 {
125  int k, upper;
126 
127  for (k = Lb(v), upper = Ub(v); (k < upper) && UlpAcc(v[k],n); k++);
128  return UlpAcc(v[k],n);
129 }
130 
131 // The 'DoubleSize' functions double the number of rows of a matrix
132 // or double the length of a vector preserving existing components.
133 //------------------------------------------------------------------
134 void DoubleSize ( ivector& x )
135 {
136  int n = Lb(x);
137  Resize(x,n,2*Ub(x)-n+1);
138 }
139 
140 
141 //accurate dot product definitions
142 
143 
144  void accumulate(idotprecision &dp, const ivector & rv1, const ivector &rv2)
145 #if(CXSC_INDEX_CHECK)
146 
147 #else
148  noexcept
149 #endif
150  {
151 #if(CXSC_INDEX_CHECK)
152  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, ivector &)"));
153 #endif
154  addDot(dp,rv1,rv2);
155  }
156 
157 
158  void accumulate(idotprecision &dp, const ivector_slice & sl, const ivector &rv)
159 #if(CXSC_INDEX_CHECK)
160 
161 #else
162  noexcept
163 #endif
164  {
165 #if(CXSC_INDEX_CHECK)
166  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, ivector &)"));
167 #endif
168  addDot(dp,sl,rv);
169  }
170 
171 
172  void accumulate(idotprecision &dp, const ivector &rv, const ivector_slice &sl)
173 #if(CXSC_INDEX_CHECK)
174 
175 #else
176  noexcept
177 #endif
178  {
179 #if(CXSC_INDEX_CHECK)
180  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, ivector_slice &)"));
181 #endif
182  addDot(dp,rv,sl);
183  }
184 
185 
186  void accumulate(idotprecision &dp, const ivector_slice & sl1, const ivector_slice &sl2)
187 #if(CXSC_INDEX_CHECK)
188 
189 #else
190  noexcept
191 #endif
192  {
193 #if(CXSC_INDEX_CHECK)
194  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, ivector_slice &)"));
195 #endif
196  addDot(dp,sl1,sl2);
197  }
198 
199 
200  void accumulate(cidotprecision &dp, const ivector & rv1, const ivector &rv2)
201 #if(CXSC_INDEX_CHECK)
202 
203 #else
204  noexcept
205 #endif
206  {
207 #if(CXSC_INDEX_CHECK)
208  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, ivector &)"));
209 #endif
210  idotprecision tmp = Re(dp);
211  tmp.set_k(dp.get_k());
212  addDot(tmp,rv1,rv2);
213  SetRe(dp,tmp);
214  }
215 
216 
217  void accumulate(cidotprecision &dp, const ivector_slice & sl, const ivector &rv)
218 #if(CXSC_INDEX_CHECK)
219 
220 #else
221  noexcept
222 #endif
223  {
224 #if(CXSC_INDEX_CHECK)
225  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, ivector &)"));
226 #endif
227  idotprecision tmp = Re(dp);
228  tmp.set_k(dp.get_k());
229  addDot(tmp,sl,rv);
230  SetRe(dp,tmp);
231  }
232 
233 
234  void accumulate(cidotprecision &dp, const ivector &rv, const ivector_slice &sl)
235 #if(CXSC_INDEX_CHECK)
236 
237 #else
238  noexcept
239 #endif
240  {
241 #if(CXSC_INDEX_CHECK)
242  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, ivector_slice &)"));
243 #endif
244  idotprecision tmp = Re(dp);
245  tmp.set_k(dp.get_k());
246  addDot(tmp,rv,sl);
247  SetRe(dp,tmp);
248  }
249 
250 
251  void accumulate(cidotprecision &dp, const ivector_slice & sl1, const ivector_slice &sl2)
252 #if(CXSC_INDEX_CHECK)
253 
254 #else
255  noexcept
256 #endif
257  {
258 #if(CXSC_INDEX_CHECK)
259  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, ivector_slice &)"));
260 #endif
261  idotprecision tmp = Re(dp);
262  tmp.set_k(dp.get_k());
263  addDot(tmp,sl1,sl2);
264  SetRe(dp,tmp);
265  }
266 
267 
268  void accumulate(idotprecision &dp, const rvector & rv1, const ivector &rv2)
269 #if(CXSC_INDEX_CHECK)
270 
271 #else
272  noexcept
273 #endif
274  {
275 #if(CXSC_INDEX_CHECK)
276  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, ivector &)"));
277 #endif
278  addDot(dp,rv1,rv2);
279  }
280 
281 
282  void accumulate(idotprecision &dp, const ivector & rv1, const rvector &rv2)
283 #if(CXSC_INDEX_CHECK)
284 
285 #else
286  noexcept
287 #endif
288  {
289 #if(CXSC_INDEX_CHECK)
290  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rvector &)"));
291 #endif
292  addDot(dp,rv1,rv2);
293  }
294 
295 
296  void accumulate(idotprecision &dp, const rvector_slice & sl, const ivector &rv)
297 #if(CXSC_INDEX_CHECK)
298 
299 #else
300  noexcept
301 #endif
302  {
303 #if(CXSC_INDEX_CHECK)
304  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, ivector &)"));
305 #endif
306  addDot(dp,sl,rv);
307  }
308 
309 
310  void accumulate(idotprecision &dp,const ivector_slice &sl,const rvector &rv)
311 #if(CXSC_INDEX_CHECK)
312 
313 #else
314  noexcept
315 #endif
316  {
317 #if(CXSC_INDEX_CHECK)
318  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rvector &)"));
319 #endif
320  addDot(dp,sl,rv);
321  }
322 
323 
324  void accumulate(idotprecision &dp, const rvector &rv, const ivector_slice &sl)
325 #if(CXSC_INDEX_CHECK)
326 
327 #else
328  noexcept
329 #endif
330  {
331 #if(CXSC_INDEX_CHECK)
332  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, ivector_slice &)"));
333 #endif
334  addDot(dp,rv,sl);
335  }
336 
337 
338  void accumulate(idotprecision &dp,const ivector &rv,const rvector_slice &sl)
339 #if(CXSC_INDEX_CHECK)
340 
341 #else
342  noexcept
343 #endif
344  {
345 #if(CXSC_INDEX_CHECK)
346  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rvector_slice &)"));
347 #endif
348  addDot(dp,rv,sl);
349  }
350 
351 
352  void accumulate(idotprecision &dp, const ivector_slice & sl1, const rvector_slice &sl2)
353 #if(CXSC_INDEX_CHECK)
354 
355 #else
356  noexcept
357 #endif
358  {
359 #if(CXSC_INDEX_CHECK)
360  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rvector_slice &)"));
361 #endif
362  addDot(dp,sl1,sl2);
363  }
364 
365 
366  void accumulate(idotprecision &dp, const rvector_slice & sl1, const ivector_slice &sl2)
367 #if(CXSC_INDEX_CHECK)
368 
369 #else
370  noexcept
371 #endif
372  {
373 #if(CXSC_INDEX_CHECK)
374  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, ivector_slice &)"));
375 #endif
376  addDot(dp,sl1,sl2);
377  }
378 
379 
380  void accumulate(cidotprecision &dp, const rvector & rv1, const ivector &rv2)
381 #if(CXSC_INDEX_CHECK)
382 
383 #else
384  noexcept
385 #endif
386  {
387 #if(CXSC_INDEX_CHECK)
388  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector &, ivector &)"));
389 #endif
390  idotprecision tmp = Re(dp);
391  tmp.set_k(dp.get_k());
392  addDot(tmp,rv1,rv2);
393  SetRe(dp,tmp);
394  }
395 
396 
397  void accumulate(cidotprecision &dp, const ivector & rv1, const rvector &rv2)
398 #if(CXSC_INDEX_CHECK)
399 
400 #else
401  noexcept
402 #endif
403  {
404 #if(CXSC_INDEX_CHECK)
405  if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, rvector &)"));
406 #endif
407  idotprecision tmp = Re(dp);
408  tmp.set_k(dp.get_k());
409  addDot(tmp,rv1,rv2);
410  SetRe(dp,tmp);
411  }
412 
413 
414  void accumulate(cidotprecision &dp, const rvector_slice & sl, const ivector &rv)
415 #if(CXSC_INDEX_CHECK)
416 
417 #else
418  noexcept
419 #endif
420  {
421 #if(CXSC_INDEX_CHECK)
422  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector_slice &, ivector &)"));
423 #endif
424  idotprecision tmp = Re(dp);
425  tmp.set_k(dp.get_k());
426  addDot(tmp,sl,rv);
427  SetRe(dp,tmp);
428  }
429 
430 
431  void accumulate(cidotprecision &dp,const ivector_slice &sl,const rvector &rv)
432 #if(CXSC_INDEX_CHECK)
433 
434 #else
435  noexcept
436 #endif
437  {
438 #if(CXSC_INDEX_CHECK)
439  if(VecLen(sl)!=VecLen(rv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, rvector &)"));
440 #endif
441  idotprecision tmp = Re(dp);
442  tmp.set_k(dp.get_k());
443  addDot(tmp,sl,rv);
444  SetRe(dp,tmp);
445  }
446 
447 
448  void accumulate(cidotprecision &dp, const rvector &rv, const ivector_slice &sl)
449 #if(CXSC_INDEX_CHECK)
450 
451 #else
452  noexcept
453 #endif
454  {
455 #if(CXSC_INDEX_CHECK)
456  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector &, ivector_slice &)"));
457 #endif
458  idotprecision tmp = Re(dp);
459  tmp.set_k(dp.get_k());
460  addDot(tmp,rv,sl);
461  SetRe(dp,tmp);
462  }
463 
464 
465  void accumulate(cidotprecision &dp,const ivector &rv,const rvector_slice &sl)
466 #if(CXSC_INDEX_CHECK)
467 
468 #else
469  noexcept
470 #endif
471  {
472 #if(CXSC_INDEX_CHECK)
473  if(VecLen(rv)!=VecLen(sl)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector &, rvector_slice &)"));
474 #endif
475  idotprecision tmp = Re(dp);
476  tmp.set_k(dp.get_k());
477  addDot(tmp,rv,sl);
478  SetRe(dp,tmp);
479  }
480 
481 
482  void accumulate(cidotprecision &dp, const ivector_slice & sl1, const rvector_slice &sl2)
483 #if(CXSC_INDEX_CHECK)
484 
485 #else
486  noexcept
487 #endif
488  {
489 #if(CXSC_INDEX_CHECK)
490  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const ivector_slice &, rvector_slice &)"));
491 #endif
492  idotprecision tmp = Re(dp);
493  tmp.set_k(dp.get_k());
494  addDot(tmp,sl1,sl2);
495  SetRe(dp,tmp);
496  }
497 
498 
499  void accumulate(cidotprecision &dp, const rvector_slice & sl1, const ivector_slice &sl2)
500 #if(CXSC_INDEX_CHECK)
501 
502 #else
503  noexcept
504 #endif
505  {
506 #if(CXSC_INDEX_CHECK)
507  if(VecLen(sl1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision &, const rvector_slice &, ivector_slice &)"));
508 #endif
509  idotprecision tmp = Re(dp);
510  tmp.set_k(dp.get_k());
511  addDot(tmp,sl1,sl2);
512  SetRe(dp,tmp);
513  }
514 
515 
516  //Summation accumulates
517  void accumulate(idotprecision &dp, const ivector& v) {
518  addSum(Inf(dp),Inf(v));
519  addSum(Sup(dp),Sup(v));
520  }
521 
522  void accumulate(cidotprecision &dp, const ivector& v) {
523  addSum(InfRe(dp),Inf(v));
524  addSum(SupRe(dp),Sup(v));
525  }
526 
527 } // namespace cxsc
528 
cxsc::mid
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition: cimatrix.inl:739
cxsc::cidotprecision::get_k
int get_k() const
Get currently set precision for computation of dot products.
Definition: cidot.hpp:89
cxsc::Ub
int Ub(const cimatrix &rm, const int &i) noexcept
Returns the upper bound index.
Definition: cimatrix.inl:1163
cxsc::Zero
int Zero(ivector &x)
Checks if vector is zero vector.
Definition: ivector.cpp:80
cxsc::rvector
The Data Type rvector.
Definition: rvector.hpp:58
cxsc::idotprecision
The Data Type idotprecision.
Definition: idot.hpp:48
cxsc::Blow
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.
Definition: cinterval.cpp:665
cxsc::in
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
Definition: cinterval.cpp:654
cxsc::ivector_slice
The Data Type ivector_slice.
Definition: ivector.hpp:963
cxsc::rvector_slice
The Data Type rvector_slice.
Definition: rvector.hpp:1064
cxsc::Lb
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.
Definition: cimatrix.inl:1156
cxsc::ivector
The Data Type ivector.
Definition: ivector.hpp:55
cxsc::idotprecision::set_k
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition: idot.hpp:88
cxsc::cidotprecision
The Data Type cidotprecision.
Definition: cidot.hpp:58
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::UlpAcc
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition: interval.cpp:335
cxsc::DoubleSize
void DoubleSize(cimatrix &A)
Doubles the size of the matrix.
Definition: cimatrix.cpp:83
cxsc::RelDiam
real RelDiam(const interval &x)
Computes the relative diameter .
Definition: interval.cpp:316
cxsc::Resize
void Resize(cimatrix &A) noexcept
Resizes the matrix.
Definition: cimatrix.inl:1211
cxsc::Disjoint
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition: interval.cpp:288
cxsc::VecLen
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966
cxsc::MaxRelDiam
real MaxRelDiam(const imatrix_subv &v)
Computes the relative diameter .
Definition: imatrix.cpp:76
cxsc::real
The Scalar Type real.
Definition: real.hpp:114