C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
rmatrix.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: rmatrix.cpp,v 1.28 2014/01/30 17:23:48 cxsc Exp $ */
25
26#define _CXSC_CPP
27
28#include "rmatrix.hpp"
29#include "matrix.inl"
30#include "rmatrix.inl"
31
32#include "dotk.inl"
33
34namespace cxsc {
35
36 //Ostrowskis comparison matrix
37 rmatrix CompMat ( const rmatrix& A) noexcept {
38 rmatrix M(Lb(A,1), Ub(A,1), Lb(A,2), Ub(A,2));
39
40 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
41 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
42 if(i-Lb(A,1) == j-Lb(A,2))
43 M[i][j] = abs(A[i][j]);
44 else
45 M[i][j] = -abs(A[i][j]);
46 }
47 }
48
49 return M;
50}
51
52 rmatrix Id ( const rmatrix& A ) // Real identity matrix
53 { //---------------------
54 int i,j;
55 int lbi = Lb(A,1), ubi = Ub(A,1);
56 int lbj = Lb(A,2), ubj = Ub(A,2);
57 rmatrix B(lbi,ubi,lbj,ubj);
58
59 for (i = lbi; i <= ubi; i++)
60 for (j = lbj; j <= ubj; j++)
61 B[i][j] = (i==j) ? 1.0 : 0.0;
62 return B;
63 }
64
65 rmatrix transp ( const rmatrix& A ) // Transposed matrix
66 { //------------------
67 int n;
68 rmatrix res(Lb(A,2),Ub(A,2),Lb(A,1),Ub(A,1));
69
70 for (n = Lb(A,1); n <= Ub(A,1); n++) Col(res,n) = Row(A,n);
71 return res;
72 }
73
74 void DoubleSize ( rmatrix& A )
75 {
76 int n = Lb(A,1);
77 Resize(A,n,2*Ub(A,1)-n+1,Lb(A,2),Ub(A,2));
78 }
79
80
81 void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
82#if(CXSC_INDEX_CHECK)
83
84#else
85 noexcept
86#endif
87 {
88#if(CXSC_INDEX_CHECK)
89 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
90#endif
91 addDot(dp,rv1,rv2);
92 }
93
94 void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2) {
95 addDot_op(dp,rv1,rv2);
96 }
97
98
99 void accumulate(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
100#if(CXSC_INDEX_CHECK)
101
102#else
103 noexcept
104#endif
105 {
106#if(CXSC_INDEX_CHECK)
107 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector &, const rmatrix_subv &)"));
108#endif
109 addDot(dp,rv1,rv2);
110 }
111
112 void accumulate_approx(dotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2) {
113 addDot_op(dp,rv1,rv2);
114 }
115
116
117 void accumulate(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
118#if(CXSC_INDEX_CHECK)
119
120#else
121 noexcept
122#endif
123 {
124#if(CXSC_INDEX_CHECK)
125 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector &)"));
126#endif
127 addDot(dp,rv1,rv2);
128 }
129
130 void accumulate_approx(dotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2) {
131 addDot_op(dp,rv1,rv2);
132 }
133
134
135 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
136#if(CXSC_INDEX_CHECK)
137
138#else
139 noexcept
140#endif
141 {
142#if(CXSC_INDEX_CHECK)
143 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
144#endif
145 dotprecision tmp(0.0);
146 tmp.set_k(dp.get_k());
147 addDot(tmp,rv1,rv2);
148 dp += tmp;
149 }
150
151
152
153 void accumulate(idotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
154#if(CXSC_INDEX_CHECK)
155
156#else
157 noexcept
158#endif
159 {
160#if(CXSC_INDEX_CHECK)
161 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector&, const rmatrix_subv &)"));
162#endif
163 dotprecision tmp(0.0);
164 tmp.set_k(dp.get_k());
165 addDot(tmp,rv1,rv2);
166 dp += tmp;
167 }
168
169
170
171 void accumulate(idotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
172#if(CXSC_INDEX_CHECK)
173
174#else
175 noexcept
176#endif
177 {
178#if(CXSC_INDEX_CHECK)
179 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rmatrix_subv&, const rvector &)"));
180#endif
181 dotprecision tmp(0.0);
182 tmp.set_k(dp.get_k());
183 addDot(tmp,rv1,rv2);
184 dp += tmp;
185 }
186
187
188
189 void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
190#if(CXSC_INDEX_CHECK)
191
192#else
193 noexcept
194#endif
195 {
196#if(CXSC_INDEX_CHECK)
197 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rmatrix_subv &)"));
198#endif
199 addDot(Re(dp),rv1,rv2);
200 }
201
202 void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rmatrix_subv &rv2)
203 {
204 addDot_op(Re(dp),rv1,rv2);
205 }
206
207
208 void accumulate(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
209#if(CXSC_INDEX_CHECK)
210
211#else
212 noexcept
213#endif
214 {
215#if(CXSC_INDEX_CHECK)
216 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector&, const rmatrix_subv &)"));
217#endif
218 addDot(Re(dp),rv1,rv2);
219 }
220
221 void accumulate_approx(cdotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
222 {
223 addDot_op(Re(dp),rv1,rv2);
224 }
225
226 void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
227#if(CXSC_INDEX_CHECK)
228
229#else
230 noexcept
231#endif
232 {
233#if(CXSC_INDEX_CHECK)
234 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector &)"));
235#endif
236 addDot(Re(dp),rv1,rv2);
237 }
238
239 void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
240 {
241 addDot_op(Re(dp),rv1,rv2);
242 }
243
244 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rmatrix_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(cidotprecision&, const rmatrix_subv &, const rmatrix_subv &)"));
253#endif
254 dotprecision tmp(0.0);
255 tmp.set_k(dp.get_k());
256 addDot(tmp,rv1,rv2);
257 dp += tmp;
258 }
259
260 void accumulate(cidotprecision &dp, const rvector & rv1, const rmatrix_subv &rv2)
261#if(CXSC_INDEX_CHECK)
262
263#else
264 noexcept
265#endif
266 {
267#if(CXSC_INDEX_CHECK)
268 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector &, const rmatrix_subv &)"));
269#endif
270 dotprecision tmp(0.0);
271 tmp.set_k(dp.get_k());
272 addDot(tmp,rv1,rv2);
273 dp += tmp;
274 }
275
276 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rvector &rv2)
277#if(CXSC_INDEX_CHECK)
278
279#else
280 noexcept
281#endif
282 {
283#if(CXSC_INDEX_CHECK)
284 if(VecLen(rv1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rvector &)"));
285#endif
286 dotprecision tmp(0.0);
287 tmp.set_k(dp.get_k());
288 addDot(tmp,rv1,rv2);
289 dp += tmp;
290 }
291
292
293 void accumulate(dotprecision &dp,const rvector_slice &sl,const rmatrix_subv &sv)
294#if(CXSC_INDEX_CHECK)
295
296#else
297 noexcept
298#endif
299 {
300#if(CXSC_INDEX_CHECK)
301 if(VecLen(sl)!=VecLen(sv)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rvector_slice &, const rmatrix_subv &)"));
302#endif
303 addDot(dp,sl,sv);
304 }
305
307 addDot_op(dp,sl,sv);
308 }
309
310
311 void accumulate(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
312#if(CXSC_INDEX_CHECK)
313
314#else
315 noexcept
316#endif
317 {
318#if(CXSC_INDEX_CHECK)
319 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rvector_slice&, const rmatrix_subv &)"));
320#endif
321 addDot(Re(dp),sl1,rv2);
322 }
323
325 {
326 addDot_op(Re(dp),sl1,rv2);
327 }
328
329
330 void accumulate(idotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
331#if(CXSC_INDEX_CHECK)
332
333#else
334 noexcept
335#endif
336 {
337#if(CXSC_INDEX_CHECK)
338 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(idotprecision&, const rvector_slice&, const rmatrix_subv &)"));
339#endif
340 dotprecision tmp(0.0);
341 tmp.set_k(dp.get_k());
342 addDot(tmp,sl1,rv2);
343 dp += tmp;
344 }
345
346
347 void accumulate(cidotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
348#if(CXSC_INDEX_CHECK)
349
350#else
351 noexcept
352#endif
353 {
354#if(CXSC_INDEX_CHECK)
355 if(VecLen(sl1)!=VecLen(rv2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rvector_slice &, const rmatrix_subv &)"));
356#endif
357 dotprecision tmp(0.0);
358 tmp.set_k(dp.get_k());
359 addDot(tmp,sl1,rv2);
360 dp += tmp;
361 }
362
363
364 void accumulate(dotprecision &dp,const rmatrix_subv &mv,const rvector_slice &vs)
365#if(CXSC_INDEX_CHECK)
366
367#else
368 noexcept
369#endif
370 {
371#if(CXSC_INDEX_CHECK)
372 if(VecLen(mv)!=VecLen(vs)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(dotprecision&, const rmatrix_subv &, const rvector_slice &)"));
373#endif
374 addDot(dp,vs,mv);
375 }
376
378 addDot_op(dp,vs,mv);
379 }
380
381
382
383 void accumulate(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
384#if(CXSC_INDEX_CHECK)
385
386#else
387 noexcept
388#endif
389 {
390#if(CXSC_INDEX_CHECK)
391 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cdotprecision&, const rmatrix_subv&, const rvector_slice &)"));
392#endif
393 addDot(Re(dp),rv1,sl2);
394 }
395
397 {
398 addDot_op(Re(dp),rv1,sl2);
399 }
400
401
402 void accumulate(idotprecision &dp, const rmatrix_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(idotprecision&, const rmatrix_subv&, const rvector_slice &)"));
411#endif
412 dotprecision tmp(0.0);
413 tmp.set_k(dp.get_k());
414 addDot(tmp,rv1,sl2);
415 dp += tmp;
416 }
417
418
419 void accumulate(cidotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
420#if(CXSC_INDEX_CHECK)
421
422#else
423 noexcept
424#endif
425 {
426#if(CXSC_INDEX_CHECK)
427 if(VecLen(rv1)!=VecLen(sl2)) cxscthrow(OP_WITH_WRONG_DIM("void accumulate(cidotprecision&, const rmatrix_subv &, const rvector_slice &)"));
428#endif
429 dotprecision tmp(0.0);
430 tmp.set_k(dp.get_k());
431 addDot(tmp,rv1,sl2);
432 dp += tmp;
433 }
434
435
436} // namespace cxsc
The Data Type cdotprecision.
Definition cdot.hpp:61
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 dotprecision.
Definition dot.hpp:112
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition dot.hpp:131
The Data Type idotprecision.
Definition idot.hpp:48
int get_k() const
Get currently set precision for computation of dot products.
Definition idot.hpp:86
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.
void accumulate_approx(cdotprecision &dp, const cmatrix_subv &rv1, const cmatrix_subv &rv2)
The accurate scalar product of the last two arguments added to the value of the first argument (witho...
Definition cmatrix.cpp:99
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
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
Definition cimatrix.cpp:74
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
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
void Resize(cimatrix &A) noexcept
Resizes the matrix.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.