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) throw() {
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  throw(OP_WITH_WRONG_DIM)
84 #else
85  throw()
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  throw(OP_WITH_WRONG_DIM)
102 #else
103  throw()
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  throw(OP_WITH_WRONG_DIM)
120 #else
121  throw()
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  throw(OP_WITH_WRONG_DIM)
138 #else
139  throw()
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  throw(OP_WITH_WRONG_DIM)
156 #else
157  throw()
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  throw(OP_WITH_WRONG_DIM)
174 #else
175  throw()
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  throw(OP_WITH_WRONG_DIM)
192 #else
193  throw()
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  throw(OP_WITH_WRONG_DIM)
211 #else
212  throw()
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  throw(OP_WITH_WRONG_DIM)
229 #else
230  throw()
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  throw(OP_WITH_WRONG_DIM)
247 #else
248  throw()
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  throw(OP_WITH_WRONG_DIM)
263 #else
264  throw()
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  throw(OP_WITH_WRONG_DIM)
279 #else
280  throw()
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  throw(OP_WITH_WRONG_DIM)
296 #else
297  throw()
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  throw(OP_WITH_WRONG_DIM)
314 #else
315  throw()
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  throw(OP_WITH_WRONG_DIM)
333 #else
334  throw()
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  throw(OP_WITH_WRONG_DIM)
350 #else
351  throw()
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  throw(OP_WITH_WRONG_DIM)
367 #else
368  throw()
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  throw(OP_WITH_WRONG_DIM)
386 #else
387  throw()
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  throw(OP_WITH_WRONG_DIM)
405 #else
406  throw()
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  throw(OP_WITH_WRONG_DIM)
422 #else
423  throw()
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
void set_k(unsigned int i)
Set precision for computation of dot products.
Definition: dot.hpp:131
The Data Type idotprecision.
Definition: idot.hpp:47
The Data Type dotprecision.
Definition: dot.hpp:111
int Lb(const cimatrix &rm, const int &i)
Returns the lower bound index.
Definition: cimatrix.inl:1156
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cimatrix_subv Col(cimatrix &m, const int &i)
Returns one column of the matrix as a vector.
Definition: cimatrix.inl:242
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966
int get_k() const
Get currently set precision for computation of dot products.
Definition: idot.hpp:86
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
void DoubleSize(cimatrix &A)
Doubles the size of the matrix.
Definition: cimatrix.cpp:83
The Data Type cidotprecision.
Definition: cidot.hpp:57
cimatrix Id(const cimatrix &A)
Returns the Identity matrix.
Definition: cimatrix.cpp:61
The Data Type rvector_slice.
Definition: rvector.hpp:1063
The Data Type rmatrix_subv.
Definition: rmatrix.hpp:53
The Data Type rvector.
Definition: rvector.hpp:57
The Data Type cdotprecision.
Definition: cdot.hpp:60
void Resize(cimatrix &A)
Resizes the matrix.
Definition: cimatrix.inl:1211
The Data Type rmatrix.
Definition: rmatrix.hpp:470
cimatrix_subv Row(cimatrix &m, const int &i)
Returns one row of the matrix as a vector.
Definition: cimatrix.inl:231
int get_k() const
Get currently set precision for computation of dot products.
Definition: cidot.hpp:89
rmatrix CompMat(const cimatrix &A)
Returns Ostrowski's comparison matrix.
Definition: cimatrix.cpp:45
int Ub(const cimatrix &rm, const int &i)
Returns the upper bound index.
Definition: cimatrix.inl:1163
cimatrix transp(const cimatrix &A)
Returns the transposed matrix.
Definition: cimatrix.cpp:74
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737