C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
imatrix.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: imatrix.cpp,v 1.29 2014/01/30 17:23:45 cxsc Exp $ */
25
26#define _CXSC_CPP
27
28#include "imatrix.hpp"
29#include "matrix.inl"
30#include "imatrix.inl"
31#include "ivecrmat.inl"
32
33#include "idotk.inl"
34
35namespace cxsc {
36
37//Ostrowskis comparison matrix
38rmatrix CompMat ( const imatrix& A) {
39 rmatrix M(Lb(A,1), Ub(A,1), Lb(A,2), Ub(A,2));
40
41 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
42 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
43 if(i-Lb(A,1) == j-Lb(A,2))
44 M[i][j] = AbsMin(A[i][j]);
45 else
46 M[i][j] = -AbsMax(A[i][j]);
47 }
48 }
49
50 return M;
51}
52
53
54imatrix Id ( const imatrix& A ) // Interval identity matrix
55{ //-------------------------
56 int i,j;
57 int lbi = Lb(A,1), ubi = Ub(A,1);
58 int lbj = Lb(A,2), ubj = Ub(A,2);
59 imatrix B(lbi,ubi,lbj,ubj);
60
61 for (i = lbi; i <= ubi; i++)
62 for (j = lbj; j <= ubj; j++)
63 B[i][j] = interval( (i==j) ? 1.0 : 0.0 );
64 return B;
65}
66
67imatrix transp ( const imatrix& A ) // Transposed matrix
68{ //------------------
69 int n;
70 imatrix res(Lb(A,2),Ub(A,2),Lb(A,1),Ub(A,1));
71
72 for (n = Lb(A,1); n <= Ub(A,1); n++) Col(res,n) = Row(A,n);
73 return res;
74}
75
76real MaxRelDiam ( const imatrix_subv& v ) // Maximum relative diameter
77{ //--------------------------
78 real r;
79 int i, l=Lb(v), u=Ub(v);
80
81 r = 0.0;
82 for (i = l; i <= u; i++)
83 if (RelDiam(v[i]) > r) r = RelDiam(v[i]);
84 return r;
85}
86
87// The 'DoubleSize' functions double the number of rows of a matrix
88// or double the length of a vector preserving existing components.
89//------------------------------------------------------------------
90
92{
93 int n = Lb(A,1);
94 Resize(A,n,2*Ub(A,1)-n+1,Lb(A,2),Ub(A,2));
95}
96
97
98 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const imatrix_subv &rv2)
99#if(CXSC_INDEX_CHECK)
100
101#else
102 noexcept
103#endif
104 {
105#if(CXSC_INDEX_CHECK)
106 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, imatrix_subv &)"));
107#endif
108 addDot(dp,rv1,rv2);
109 }
110
111 void accumulate(idotprecision &dp, const ivector & rv1, const imatrix_subv &rv2)
112#if(CXSC_INDEX_CHECK)
113
114#else
115 noexcept
116#endif
117 {
118#if(CXSC_INDEX_CHECK)
119 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, imatrix_subv &)"));
120#endif
121 addDot(dp,rv1,rv2);
122 }
123
124 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const ivector &rv2)
125#if(CXSC_INDEX_CHECK)
126
127#else
128 noexcept
129#endif
130 {
131#if(CXSC_INDEX_CHECK)
132 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, ivector &)"));
133#endif
134 addDot(dp,rv1,rv2);
135 }
136
137
138 void accumulate(idotprecision &dp, const ivector_slice & sl1, const imatrix_subv &rv2)
139#if(CXSC_INDEX_CHECK)
140
141#else
142 noexcept
143#endif
144 {
145#if(CXSC_INDEX_CHECK)
146 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, imatrix_subv &)"));
147#endif
148 addDot(dp,sl1,rv2);
149 }
150
151 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const ivector_slice &sl2)
152#if(CXSC_INDEX_CHECK)
153
154#else
155 noexcept
156#endif
157 {
158#if(CXSC_INDEX_CHECK)
159 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, ivector_slice &)"));
160#endif
161 addDot(dp,rv1,sl2);
162 }
163
164 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const imatrix_subv &rv2)
165#if(CXSC_INDEX_CHECK)
166
167#else
168 noexcept
169#endif
170 {
171#if(CXSC_INDEX_CHECK)
172 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const imatrix_subv &)"));
173#endif
174 idotprecision tmp(0.0);
175 tmp.set_k(dp.get_k());
176 addDot(tmp,rv1,rv2);
177 dp += tmp;
178 }
179
180 void accumulate(cidotprecision &dp, const ivector & rv1, const imatrix_subv &rv2)
181#if(CXSC_INDEX_CHECK)
182
183#else
184 noexcept
185#endif
186 {
187#if(CXSC_INDEX_CHECK)
188 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector &, const imatrix_subv &)"));
189#endif
190 idotprecision tmp(0.0);
191 tmp.set_k(dp.get_k());
192 addDot(tmp,rv1,rv2);
193 dp += tmp;
194 }
195
196 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const ivector &rv2)
197#if(CXSC_INDEX_CHECK)
198
199#else
200 noexcept
201#endif
202 {
203#if(CXSC_INDEX_CHECK)
204 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const ivector &)"));
205#endif
206 idotprecision tmp(0.0);
207 tmp.set_k(dp.get_k());
208 addDot(tmp,rv1,rv2);
209 dp += tmp;
210 }
211
212 void accumulate(cidotprecision &dp, const ivector_slice & sl1, const imatrix_subv &rv2)
213#if(CXSC_INDEX_CHECK)
214
215#else
216 noexcept
217#endif
218 {
219#if(CXSC_INDEX_CHECK)
220 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector_slice &, const imatrix_subv &)"));
221#endif
222 idotprecision tmp(0.0);
223 tmp.set_k(dp.get_k());
224 addDot(tmp,sl1,rv2);
225 dp += tmp;
226 }
227
228 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const ivector_slice &sl2)
229#if(CXSC_INDEX_CHECK)
230
231#else
232 noexcept
233#endif
234 {
235#if(CXSC_INDEX_CHECK)
236 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const ivector_slice &)"));
237#endif
238 idotprecision tmp(0.0);
239 tmp.set_k(dp.get_k());
240 addDot(tmp,rv1,sl2);
241 dp += tmp;
242 }
243
244 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const imatrix_subv &rv2)
245#if(CXSC_INDEX_CHECK)
246
247#else
248 noexcept
249#endif
250 {
251#if(CXSC_INDEX_CHECK)
252 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rmatrix_subv &, imatrix_subv &)"));
253#endif
254 addDot(dp,rv1,rv2);
255 }
256
257 void accumulate(idotprecision &dp, const rvector & rv1, const imatrix_subv &rv2)
258#if(CXSC_INDEX_CHECK)
259
260#else
261 noexcept
262#endif
263 {
264#if(CXSC_INDEX_CHECK)
265 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector &, imatrix_subv &)"));
266#endif
267 addDot(dp,rv1,rv2);
268 }
269
270 void accumulate(idotprecision &dp, const rvector_slice & sl1, const imatrix_subv &rv2)
271#if(CXSC_INDEX_CHECK)
272
273#else
274 noexcept
275#endif
276 {
277#if(CXSC_INDEX_CHECK)
278 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rvector_slice &, imatrix_subv &)"));
279#endif
280 addDot(dp,sl1,rv2);
281 }
282
283 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rmatrix_subv &rv2)
284#if(CXSC_INDEX_CHECK)
285
286#else
287 noexcept
288#endif
289 {
290#if(CXSC_INDEX_CHECK)
291 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rmatrix_subv &)"));
292#endif
293 addDot(dp,rv1,rv2);
294 }
295
296 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rvector &rv2)
297#if(CXSC_INDEX_CHECK)
298
299#else
300 noexcept
301#endif
302 {
303#if(CXSC_INDEX_CHECK)
304 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rvector &)"));
305#endif
306 addDot(dp,rv1,rv2);
307 }
308
309 void accumulate(idotprecision &dp, const imatrix_subv & rv1, const rvector_slice &sl2)
310#if(CXSC_INDEX_CHECK)
311
312#else
313 noexcept
314#endif
315 {
316#if(CXSC_INDEX_CHECK)
317 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const imatrix_subv &, rvector_slice &)"));
318#endif
319 addDot(dp,rv1,sl2);
320 }
321
322 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const imatrix_subv &rv2)
323#if(CXSC_INDEX_CHECK)
324
325#else
326 noexcept
327#endif
328 {
329#if(CXSC_INDEX_CHECK)
330 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const imatrix_subv &)"));
331#endif
332 idotprecision tmp(0.0);
333 tmp.set_k(dp.get_k());
334 addDot(tmp,rv1,rv2);
335 dp += tmp;
336 }
337
338 void accumulate(cidotprecision &dp, const rvector & rv1, const imatrix_subv &rv2)
339#if(CXSC_INDEX_CHECK)
340
341#else
342 noexcept
343#endif
344 {
345#if(CXSC_INDEX_CHECK)
346 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector &, const imatrix_subv &)"));
347#endif
348 idotprecision tmp(0.0);
349 tmp.set_k(dp.get_k());
350 addDot(tmp,rv1,rv2);
351 dp += tmp;
352 }
353
354 void accumulate(cidotprecision &dp, const rvector_slice & sl1, const imatrix_subv &rv2)
355#if(CXSC_INDEX_CHECK)
356
357#else
358 noexcept
359#endif
360 {
361#if(CXSC_INDEX_CHECK)
362 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector_slice&, const imatrix_subv &)"));
363#endif
364 idotprecision tmp(0.0);
365 tmp.set_k(dp.get_k());
366 addDot(tmp,sl1,rv2);
367 dp += tmp;
368 }
369
370 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const rmatrix_subv &rv2)
371#if(CXSC_INDEX_CHECK)
372
373#else
374 noexcept
375#endif
376 {
377#if(CXSC_INDEX_CHECK)
378 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const rmatrix_subv &)"));
379#endif
380 idotprecision tmp(0.0);
381 tmp.set_k(dp.get_k());
382 addDot(tmp,rv1,rv2);
383 dp += tmp;
384 }
385
386 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const rvector &rv2)
387#if(CXSC_INDEX_CHECK)
388
389#else
390 noexcept
391#endif
392 {
393#if(CXSC_INDEX_CHECK)
394 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const rvector &)"));
395#endif
396 idotprecision tmp(0.0);
397 tmp.set_k(dp.get_k());
398 addDot(tmp,rv1,rv2);
399 dp += tmp;
400 }
401
402 void accumulate(cidotprecision &dp, const imatrix_subv & rv1, const rvector_slice &sl2)
403#if(CXSC_INDEX_CHECK)
404
405#else
406 noexcept
407#endif
408 {
409#if(CXSC_INDEX_CHECK)
410 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const imatrix_subv &, const rvector_slice &)"));
411#endif
412 idotprecision tmp(0.0);
413 tmp.set_k(dp.get_k());
414 addDot(tmp,rv1,sl2);
415 dp += tmp;
416 }
417
418 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const ivector &rv2)
419#if(CXSC_INDEX_CHECK)
420
421#else
422 noexcept
423#endif
424 {
425#if(CXSC_INDEX_CHECK)
426 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rmatrix_subv &, ivector &)"));
427#endif
428 addDot(dp,rv1,rv2);
429 }
430
431 void accumulate(idotprecision &dp, const ivector & rv1, const rmatrix_subv &rv2)
432#if(CXSC_INDEX_CHECK)
433
434#else
435 noexcept
436#endif
437 {
438#if(CXSC_INDEX_CHECK)
439 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector &, rmatrix_subv &)"));
440#endif
441 addDot(dp,rv1,rv2);
442 }
443
444 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const ivector &rv2)
445#if(CXSC_INDEX_CHECK)
446
447#else
448 noexcept
449#endif
450 {
451#if(CXSC_INDEX_CHECK)
452 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const ivector &)"));
453#endif
454 idotprecision tmp(0.0);
455 tmp.set_k(dp.get_k());
456 addDot(tmp,rv1,rv2);
457 dp += tmp;
458 }
459
460 void accumulate(cidotprecision &dp, const ivector & rv1, const rmatrix_subv &rv2)
461#if(CXSC_INDEX_CHECK)
462
463#else
464 noexcept
465#endif
466 {
467#if(CXSC_INDEX_CHECK)
468 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector &, const rmatrix_subv &)"));
469#endif
470 idotprecision tmp(0.0);
471 tmp.set_k(dp.get_k());
472 addDot(tmp,rv1,rv2);
473 dp += tmp;
474 }
475
476 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const ivector_slice &rv2)
477#if(CXSC_INDEX_CHECK)
478
479#else
480 noexcept
481#endif
482 {
483#if(CXSC_INDEX_CHECK)
484 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const rmatrix_subv &, ivector_slice &)"));
485#endif
486 addDot(dp,rv1,rv2);
487 }
488
489 void accumulate(idotprecision &dp, const ivector_slice & rv1, const rmatrix_subv &rv2)
490#if(CXSC_INDEX_CHECK)
491
492#else
493 noexcept
494#endif
495 {
496#if(CXSC_INDEX_CHECK)
497 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision &, const ivector_slice &, rmatrix_subv &)"));
498#endif
499 addDot(dp,rv1,rv2);
500 }
501
502 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const ivector_slice &rv2)
503#if(CXSC_INDEX_CHECK)
504
505#else
506 noexcept
507#endif
508 {
509#if(CXSC_INDEX_CHECK)
510 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const ivector_slice &)"));
511#endif
512 idotprecision tmp(0.0);
513 tmp.set_k(dp.get_k());
514 addDot(tmp,rv1,rv2);
515 dp += tmp;
516 }
517
518 void accumulate(cidotprecision &dp, const ivector_slice & rv1, const rmatrix_subv &rv2)
519#if(CXSC_INDEX_CHECK)
520
521#else
522 noexcept
523#endif
524 {
525#if(CXSC_INDEX_CHECK)
526 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const ivector_slice &, const rmatrix_subv &)"));
527#endif
528 idotprecision tmp(0.0);
529 tmp.set_k(dp.get_k());
530 addDot(tmp,rv1,rv2);
531 dp += tmp;
532 }
533
534
535} // namespace cxsc
536
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 imatrix_subv.
Definition imatrix.hpp:56
The Data Type imatrix.
Definition imatrix.hpp:660
The Scalar Type interval.
Definition interval.hpp:55
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 rmatrix_subv.
Definition rmatrix.hpp:54
The Data Type rmatrix.
Definition rmatrix.hpp:471
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.
cimatrix_subv Col(cimatrix &m, const int &i) noexcept
Returns one column of the matrix as a vector.
Definition cimatrix.inl:242
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
Definition cimatrix.cpp:45
real MaxRelDiam(const imatrix_subv &v)
Computes the relative diameter .
Definition imatrix.cpp:76
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
Definition cimatrix.cpp:74
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.
cimatrix_subv Row(cimatrix &m, const int &i) noexcept
Returns one row of the matrix as a vector.
Definition cimatrix.inl:231
void DoubleSize(cimatrix &A)
Doubles the size of the matrix.
Definition cimatrix.cpp:83
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
Definition cimatrix.cpp:61
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition interval.cpp:303
void Resize(cimatrix &A) noexcept
Resizes the matrix.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.
real AbsMin(const interval &x)
Computes the smallest absolute value .
Definition interval.cpp:293