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 
34 namespace 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 
324  void accumulate_approx(cdotprecision &dp, const rvector_slice & sl1, const rmatrix_subv &rv2)
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 
396  void accumulate_approx(cdotprecision &dp, const rmatrix_subv & rv1, const rvector_slice &sl2)
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
cxsc::idotprecision::get_k
int get_k() const
Get currently set precision for computation of dot products.
Definition: idot.hpp:86
cxsc::CompMat
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
Definition: cimatrix.cpp:45
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::rmatrix
The Data Type rmatrix.
Definition: rmatrix.hpp:471
cxsc::abs
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::rvector
The Data Type rvector.
Definition: rvector.hpp:58
cxsc::idotprecision
The Data Type idotprecision.
Definition: idot.hpp:48
cxsc::Col
cimatrix_subv Col(cimatrix &m, const int &i) noexcept
Returns one column of the matrix as a vector.
Definition: cimatrix.inl:242
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::accumulate_approx
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
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:112
cxsc::rmatrix_subv
The Data Type rmatrix_subv.
Definition: rmatrix.hpp:54
cxsc::cidotprecision
The Data Type cidotprecision.
Definition: cidot.hpp:58
cxsc::dotprecision::set_k
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition: dot.hpp:131
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::DoubleSize
void DoubleSize(cimatrix &A)
Doubles the size of the matrix.
Definition: cimatrix.cpp:83
cxsc::transp
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
Definition: cimatrix.cpp:74
cxsc::cdotprecision
The Data Type cdotprecision.
Definition: cdot.hpp:61
cxsc::Id
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
Definition: cimatrix.cpp:61
cxsc::Resize
void Resize(cimatrix &A) noexcept
Resizes the matrix.
Definition: cimatrix.inl:1211
cxsc::Row
cimatrix_subv Row(cimatrix &m, const int &i) noexcept
Returns one row of the matrix as a vector.
Definition: cimatrix.inl:231
cxsc::VecLen
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966