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
34namespace cxsc {
35
36int 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
44int 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
53ivector 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
62int 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
80int 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
87rvector 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
97real 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
108real 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//----------------------------------------------------------------------------
123int 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//------------------------------------------------------------------
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
The Data Type cidotprecision.
Definition cidot.hpp:58
int get_k() const
Get currently set precision for computation of dot products.
Definition cidot.hpp:89
The Data Type idotprecision.
Definition idot.hpp:48
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition idot.hpp:88
The Data Type ivector_slice.
Definition ivector.hpp:963
The Data Type ivector.
Definition ivector.hpp:55
The Scalar Type real.
Definition real.hpp:114
The Data Type rvector_slice.
Definition rvector.hpp:1064
The Data Type rvector.
Definition rvector.hpp:58
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition interval.cpp:288
real MaxRelDiam(const imatrix_subv &v)
Computes the relative diameter .
Definition imatrix.cpp:76
int UlpAcc(const interval &x, int n)
Checks if the diameter of the interval is ulps.
Definition interval.cpp:335
real RelDiam(const interval &x)
Computes the relative diameter .
Definition interval.cpp:316
int Ub(const cimatrix &rm, const int &i) noexcept
Returns the upper bound index.
int in(const cinterval &x, const cinterval &y)
Checks if first argument is part of second argument.
void DoubleSize(cimatrix &A)
Doubles the size of the matrix.
Definition cimatrix.cpp:83
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739
void Resize(cimatrix &A) noexcept
Resizes the matrix.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.
int Zero(ivector &x)
Checks if vector is zero vector.
Definition ivector.cpp:80
cinterval Blow(cinterval x, const real &eps)
Performs an epsilon inflation.