C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
cxsc_blas.inl
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: cxsc_blas.inl,v 1.20 2014/01/30 17:23:44 cxsc Exp $ */
25 
26 /*
27 ** FastPLSS: A library of (parallel) verified linear (interval) system
28 ** solvers using C-XSC (V 0.2)
29 */
30 
31 #include <fenv.h>
32 #include "cxsc_blas.hpp"
33 
34 namespace cxsc {
35 
36 #ifndef _CXSC_BLAS_SETROUND
37 #define _CXSC_BLAS_SETROUND
38  /*
39  Sets the rounding mode according to the parameter r:
40  r=-1: Round downwards
41  r=0 : Round to nearest
42  r=1 : Round upwards
43  r=2 : Round toward zero
44  */
45  static inline void setround(int r) {
46  switch(r) {
47  case -1:
48  fesetround(FE_DOWNWARD);
49  break;
50  case 0:
51  fesetround(FE_TONEAREST);
52  break;
53  case 1:
54  fesetround(FE_UPWARD);
55  break;
56  case 2:
57  fesetround(FE_TOWARDZERO);
58  break;
59  default:
60  fesetround(FE_TONEAREST);
61  }
62  }
63 
64  /*
65  Sets the rounding mode according to the parameter r:
66  r=-1: Round downwards
67  r=0 : Round to nearest
68  r=1 : Round upwards
69  r=2 : Round toward zero
70  */
71  static inline int getround() {
72  switch(fegetround()) {
73  case FE_DOWNWARD:
74  return -1;
75  break;
76  case FE_TONEAREST:
77  return 0;
78  break;
79  case FE_UPWARD:
80  return 1;
81  break;
82  case FE_TOWARDZERO:
83  return 2;
84  break;
85  default:
86  return 0;
87  }
88  }
89 #endif
90 
91 #ifdef _CXSC_BLAS_RVECTOR
92 #ifndef _CXSC_BLAS_RVECTOR_INC
93 #define _CXSC_BLAS_RVECTOR_INC
94 inline void blasdot(const rvector& x, const rvector& y, real& res) {
95  int n = VecLen(x);
96  int inc=1;
97 
98  double* xd = (double*) x.to_blas_array();
99  double* yd = (double*) y.to_blas_array();
100 
101  res = cblas_ddot(n, xd, inc, yd, inc);
102 }
103 
104 inline void blasdot(const rvector& x, const rvector& y, interval& res) {
105  int n = VecLen(x);
106  int inc=1;
107  int rnd = getround();
108 
109  double* xd = (double*) x.to_blas_array();
110  double* yd = (double*) y.to_blas_array();
111 
112  setround(-1);
113  SetInf(res, cblas_ddot(n, xd, inc, yd, inc));
114  setround(1);
115  SetSup(res, cblas_ddot(n, xd, inc, yd, inc));
116  setround(rnd);
117 }
118 #endif
119 #endif
120 
121 #if defined(_CXSC_BLAS_CVECTOR) && defined(_CXSC_BLAS_IVECTOR)
122 #if !defined(_CXSC_BLAS_CVECTOR_INC) || !defined(_CXSC_BLAS_IVECTOR_INC)
123 inline void blasdot(const cvector& x, const ivector& y, cinterval& res) {
124  const int n = VecLen(x);
125  const int inc=1;
126  const int lb1 = Lb(x);
127  const int lb2 = Lb(y);
128  int rnd = getround();
129 
130  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
131 
132  double* dxi = x_inf.to_blas_array();
133  double* dxs = x_sup.to_blas_array();
134  double* dyi = y_inf.to_blas_array();
135  double* dys = y_sup.to_blas_array();
136  double res_inf,res_sup;
137 
138  setround(1);
139 
140  bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
141 
142  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
143 
144  //setround(-1);
145 
146  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
147 
148  SetRe(res, interval(res_inf,res_sup));
149 
150  bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
151 
152  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
153 
154  //setround(-1);
155 
156  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
157 
158  SetIm(res, interval(res_inf,res_sup));
159 
160  setround(rnd);
161 }
162 
163 inline void blasdot(const ivector& x, const cvector& y, cinterval& res) {
164  const int n = VecLen(x);
165  const int inc=1;
166  const int lb1 = Lb(x);
167  const int lb2 = Lb(y);
168  int rnd = getround();
169 
170  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
171 
172  double* dxi = x_inf.to_blas_array();
173  double* dxs = x_sup.to_blas_array();
174  double* dyi = y_inf.to_blas_array();
175  double* dys = y_sup.to_blas_array();
176  double res_inf,res_sup;
177 
178  setround(1);
179 
180  bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
181 
182  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
183 
184  //setround(-1);
185 
186  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
187 
188  SetRe(res, interval(res_inf,res_sup));
189 
190  bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
191 
192  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
193 
194  //setround(-1);
195 
196  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
197 
198  SetIm(res, interval(res_inf,res_sup));
199 
200  setround(rnd);
201 }
202 #endif
203 #endif
204 
205 #ifdef _CXSC_BLAS_CVECTOR
206 #ifndef _CXSC_BLAS_CVECTOR_INC
207 #define _CXSC_BLAS_CVECTOR_INC
208 inline void blasdot(const rvector& x, const cvector& y, complex& res) {
209  int n = VecLen(x);
210  int inc=1, inc2=2;
211  double res_re, res_im;
212 
213  double* xd = (double*) x.to_blas_array();
214  double* yd = (double*) y.to_blas_array();
215 
216  res_re = cblas_ddot(n, xd, inc, yd, inc2);
217  res_im = cblas_ddot(n, xd, inc, yd+1, inc2);
218 
219  res = complex(res_re,res_im);
220 }
221 
222 inline void blasdot(const cvector& x, const rvector& y, complex& res) {
223  int n = VecLen(x);
224  int inc=1, inc2=2;
225  double res_re, res_im;
226 
227  double* xd = (double*) x.to_blas_array();
228  double* yd = (double*) y.to_blas_array();
229 
230  res_re = cblas_ddot(n, xd, inc2, yd, inc);
231  res_im = cblas_ddot(n, xd+1, inc2, yd, inc);
232 
233  res = complex(res_re,res_im);
234 }
235 
236 inline void blasdot(const cvector& x, const cvector& y, complex& res) {
237  int n = VecLen(x);
238  int inc=1;
239 
240  double* xd = (double*) x.to_blas_array();
241  double* yd = (double*) y.to_blas_array();
242 
243  cblas_zdotu_sub(n, xd, inc, yd, inc, (void*)&res);
244 }
245 
246 inline void blasdot(const rvector& x, const cvector& y, cinterval& res) {
247  interval re,im;
248 
249  blasdot(x,Re(y),re);
250  blasdot(x,Im(y),im);
251 
252  SetRe(res,re);
253  SetIm(res,im);
254 }
255 
256 inline void blasdot(const cvector& x, const rvector& y, cinterval& res) {
257  interval re,im;
258 
259  blasdot(Re(x),y,re);
260  blasdot(Im(x),y,im);
261 
262  SetRe(res,re);
263  SetIm(res,im);
264 }
265 
266 
267 inline void blasdot(const cvector& x, const cvector& y, cinterval& res) {
268  int rnd = getround();
269  int n = VecLen(x);
270  int inc=2;
271  double res_re, res_im;
272 
273  double* xd = (double*) x.to_blas_array();
274  double* yd = (double*) y.to_blas_array();
275 
276  setround(1);
277  res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc);
278 
279  setround(-1);
280  res_re += cblas_ddot(n, xd, inc, yd, inc);
281  res_im = cblas_ddot(n, xd, inc, yd+1, inc);
282  res_im += cblas_ddot(n, xd+1, inc, yd, inc);
283  UncheckedSetInf(res,complex(res_re,res_im));
284 
285  res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc);
286  setround(1);
287  res_re += cblas_ddot(n, xd, inc, yd, inc);
288  res_im = cblas_ddot(n, xd, inc, yd+1, inc);
289  res_im += cblas_ddot(n, xd+1, inc, yd, inc);
290  UncheckedSetSup(res,complex(res_re,res_im));
291 
292  setround(rnd);
293 }
294 #endif
295 #endif
296 
297 #ifdef _CXSC_BLAS_IVECTOR
298 #ifndef _CXSC_BLAS_IVECTOR_INC
299 #define _CXSC_BLAS_IVECTOR_INC
300 inline void bsort(const ivector &x, const ivector &y, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb1, int lb2) {
301 
302  for(int i=0 ; i<n ; i++) {
303 
304  if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) {
305 
306  if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) {
307  x_inf[i+1] = -Inf(x[i+lb1]);
308  y_inf[i+1] = Inf(y[i+lb2]);
309  x_sup[i+1] = Sup(x[i+lb1]);
310  y_sup[i+1] = Sup(y[i+lb2]);
311  } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) {
312  x_inf[i+1] = -Sup(x[i+lb1]);
313  y_inf[i+1] = Inf(y[i+lb2]);
314  x_sup[i+1] = Sup(x[i+lb1]);
315  y_sup[i+1] = Sup(y[i+lb2]);
316  } else {
317  x_inf[i+1] = -Sup(x[i+lb1]);
318  y_inf[i+1] = Inf(y[i+lb2]);
319  x_sup[i+1] = Inf(x[i+lb1]);
320  y_sup[i+1] = Sup(y[i+lb2]);
321  }
322 
323  } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) {
324 
325  if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) {
326  x_inf[i+1] = -Inf(x[i+lb1]);
327  y_inf[i+1] = Sup(y[i+lb2]);
328  x_sup[i+1] = Sup(x[i+lb1]);
329  y_sup[i+1] = Sup(y[i+lb2]);
330  } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) {
331 
332  if((-Inf(x[i+lb1]))*Sup(y[i+lb2]) >= Sup(x[i+lb1])*(-Inf(y[i+lb2]))) {
333  x_inf[i+1] = -Inf(x[i+lb1]);
334  y_inf[i+1] = Sup(y[i+lb2]);
335  } else {
336  x_inf[i+1] = -Sup(x[i+lb1]);
337  y_inf[i+1] = Inf(y[i+lb2]);
338  }
339 
340  if(Inf(x[i+lb1])*Inf(y[i+lb2]) >= Sup(x[i+lb1])*Sup(y[i+lb2])) {
341  x_sup[i+1] = Inf(x[i+lb1]);
342  y_sup[i+1] = Inf(y[i+lb2]);
343  } else {
344  x_sup[i+1] = Sup(x[i+lb1]);
345  y_sup[i+1] = Sup(y[i+lb2]);
346  }
347 
348  } else {
349  x_inf[i+1] = Sup(x[i+lb1]);
350  y_inf[i+1] = -Inf(y[i+lb2]);
351  x_sup[i+1] = Inf(x[i+lb1]);
352  y_sup[i+1] = Inf(y[i+lb2]);
353  }
354 
355  } else {
356 
357  if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) {
358  x_inf[i+1] = -Inf(x[i+lb1]);
359  y_inf[i+1] = Sup(y[i+lb2]);
360  x_sup[i+1] = Sup(x[i+lb1]);
361  y_sup[i+1] = Inf(y[i+lb2]);
362  } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) {
363  x_inf[i+1] = -Inf(x[i+lb1]);
364  y_inf[i+1] = Sup(y[i+lb2]);
365  x_sup[i+1] = Inf(x[i+lb1]);
366  y_sup[i+1] = Inf(y[i+lb2]);
367  } else {
368  x_inf[i+1] = -Sup(x[i+lb1]);
369  y_inf[i+1] = Sup(y[i+lb2]);
370  x_sup[i+1] = Inf(x[i+lb1]);
371  y_sup[i+1] = Inf(y[i+lb2]);
372  }
373 
374  }
375 
376  }
377 }
378 
379 //Sorts inf and sup of an interval vector and a real vector into two real vector for the computation of the infimum
380 //and two real vectors for the computation of the supremum. Rounding mode must be set to upwards!
381 inline void bsort(const ivector &x, const rvector &y, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb1, int lb2) {
382 
383  for(int i=0 ; i<n ; i++) {
384 
385  if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) {
386 
387  if(y[i+lb2] >= 0) {
388  x_inf[i+1] = -Inf(x[i+lb1]);
389  y_inf[i+1] = y[i+lb2];
390  x_sup[i+1] = Sup(x[i+lb1]);
391  y_sup[i+1] = y[i+lb2];
392  } else {
393  x_inf[i+1] = -Sup(x[i+lb1]);
394  y_inf[i+1] = y[i+lb2];
395  x_sup[i+1] = Inf(x[i+lb1]);
396  y_sup[i+1] = y[i+lb2];
397  }
398 
399  } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) {
400 
401  if(y[i+lb2] >= 0) {
402  x_inf[i+1] = -Inf(x[i+lb1]);
403  y_inf[i+1] = y[i+lb2];
404  x_sup[i+1] = Sup(x[i+lb1]);
405  y_sup[i+1] = y[i+lb2];
406  } else {
407  x_inf[i+1] = -Sup(x[i+lb1]);
408  y_inf[i+1] = y[i+lb2];
409  x_sup[i+1] = Inf(x[i+lb1]);
410  y_sup[i+1] = y[i+lb2];
411  }
412 
413  } else {
414 
415  if(y[i+lb2] >= 0) {
416  x_inf[i+1] = -Inf(x[i+lb1]);
417  y_inf[i+1] = y[i+lb2];
418  x_sup[i+1] = Sup(x[i+lb1]);
419  y_sup[i+1] = y[i+lb2];
420  } else {
421  x_inf[i+1] = -Sup(x[i+lb1]);
422  y_inf[i+1] = y[i+lb2];
423  x_sup[i+1] = Inf(x[i+lb1]);
424  y_sup[i+1] = y[i+lb2];
425  }
426 
427  }
428 
429  }
430 }
431 
432 //Sorts inf and sup of an interval vector and a real vector into two real vector for the computation of the infimum
433 //and two real vectors for the computation of the supremum. Rounding mode must be set to upwards!
434 inline void bsort(const rvector &y, const ivector &x, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb2, int lb1) {
435  for(int i=0 ; i<n ; i++) {
436 
437  if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) {
438 
439  if(y[i+lb2] >= 0) {
440  x_inf[i+1] = -Inf(x[i+lb1]);
441  y_inf[i+1] = y[i+lb2];
442  x_sup[i+1] = Sup(x[i+lb1]);
443  y_sup[i+1] = y[i+lb2];
444  } else {
445  x_inf[i+1] = -Sup(x[i+lb1]);
446  y_inf[i+1] = y[i+lb2];
447  x_sup[i+1] = Inf(x[i+lb1]);
448  y_sup[i+1] = y[i+lb2];
449  }
450 
451  } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) {
452 
453  if(y[i+lb2] >= 0) {
454  x_inf[i+1] = -Inf(x[i+lb1]);
455  y_inf[i+1] = y[i+lb2];
456  x_sup[i+1] = Sup(x[i+lb1]);
457  y_sup[i+1] = y[i+lb2];
458  } else {
459  x_inf[i+1] = -Sup(x[i+lb1]);
460  y_inf[i+1] = y[i+lb2];
461  x_sup[i+1] = Inf(x[i+lb1]);
462  y_sup[i+1] = y[i+lb2];
463  }
464 
465  } else {
466 
467  if(y[i+lb2] >= 0) {
468  x_inf[i+1] = -Inf(x[i+lb1]);
469  y_inf[i+1] = y[i+lb2];
470  x_sup[i+1] = Sup(x[i+lb1]);
471  y_sup[i+1] = y[i+lb2];
472  } else {
473  x_inf[i+1] = -Sup(x[i+lb1]);
474  y_inf[i+1] = y[i+lb2];
475  x_sup[i+1] = Inf(x[i+lb1]);
476  y_sup[i+1] = y[i+lb2];
477  }
478 
479  }
480 
481  }
482 }
483 
484 inline void blasdot(const ivector& x, const ivector& y, interval& res) {
485  const int n = VecLen(x);
486  const int inc=1;
487  const int lb1 = Lb(x);
488  const int lb2 = Lb(y);
489  int rnd = getround();
490 
491  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
492 
493  double* dxi = x_inf.to_blas_array();
494  double* dxs = x_sup.to_blas_array();
495  double* dyi = y_inf.to_blas_array();
496  double* dys = y_sup.to_blas_array();
497  double res_inf,res_sup;
498 
499  setround(1);
500 
501  bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
502 
503  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
504 
505  //setround(-1);
506 
507  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
508 
509  setround(rnd);
510 
511  res = interval(res_inf,res_sup);
512 }
513 
514 inline void blasdot(const ivector& x, const rvector& y, interval& res) {
515  const int n = VecLen(x);
516  const int inc=1;
517  const int lb1 = Lb(x);
518  const int lb2 = Lb(y);
519  int rnd = getround();
520 
521  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
522 
523  double* dxi = x_inf.to_blas_array();
524  double* dxs = x_sup.to_blas_array();
525  double* dyi = y_inf.to_blas_array();
526  double* dys = y_sup.to_blas_array();
527  double res_inf,res_sup;
528 
529  setround(1);
530 
531  bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
532 
533  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
534 
535  //setround(-1);
536 
537  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
538 
539  setround(rnd);
540 
541  res = interval(res_inf,res_sup);
542 }
543 
544 inline void blasdot(const rvector& x, const ivector& y, interval& res) {
545  const int n = VecLen(x);
546  const int inc=1;
547  const int lb1 = Lb(x);
548  const int lb2 = Lb(y);
549  int rnd = getround();
550 
551  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
552 
553  double* dxi = x_inf.to_blas_array();
554  double* dxs = x_sup.to_blas_array();
555  double* dyi = y_inf.to_blas_array();
556  double* dys = y_sup.to_blas_array();
557  double res_inf,res_sup;
558 
559  setround(1);
560 
561  bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
562 
563  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
564 
565  //setround(-1);
566 
567  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
568 
569  setround(rnd);
570 
571  res = interval(res_inf,res_sup);
572 }
573 #endif
574 #endif
575 
576 
577 #ifdef _CXSC_BLAS_CIVECTOR
578 #ifndef _CXSC_BLAS_CIVECTOR_INC
579 #define _CXSC_BLAS_CIVECTOR_INC
580 inline void blasdot(const rvector& x, const civector& y, cinterval& res) {
581  const int n = VecLen(x);
582  const int inc=1;
583  const int lb1 = Lb(x);
584  const int lb2 = Lb(y);
585  int rnd = getround();
586 
587  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
588 
589  double* dxi = x_inf.to_blas_array();
590  double* dxs = x_sup.to_blas_array();
591  double* dyi = y_inf.to_blas_array();
592  double* dys = y_sup.to_blas_array();
593  double res_inf,res_sup;
594 
595  setround(1);
596 
597  bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
598 
599  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
600 
601  //setround(-1);
602 
603  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
604 
605  SetRe(res, interval(res_inf,res_sup));
606 
607  bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
608 
609  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
610 
611  //setround(-1);
612 
613  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
614 
615  SetIm(res, interval(res_inf,res_sup));
616 
617  setround(rnd);
618 }
619 
620 
621 inline void blasdot(const civector& x, const rvector& y, cinterval& res) {
622  const int n = VecLen(x);
623  const int inc=1;
624  const int lb1 = Lb(x);
625  const int lb2 = Lb(y);
626  int rnd = getround();
627 
628  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
629 
630  double* dxi = x_inf.to_blas_array();
631  double* dxs = x_sup.to_blas_array();
632  double* dyi = y_inf.to_blas_array();
633  double* dys = y_sup.to_blas_array();
634  double res_inf,res_sup;
635 
636  setround(1);
637 
638  bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
639 
640  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
641 
642  //setround(-1);
643 
644  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
645 
646  SetRe(res, interval(res_inf,res_sup));
647 
648  bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
649 
650  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
651 
652  //setround(-1);
653 
654  res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
655 
656  SetIm(res, interval(res_inf,res_sup));
657 
658  setround(rnd);
659 }
660 
661 inline void blasdot(const civector& x, const civector& y, cinterval& res) {
662  const int n = VecLen(x);
663  const int inc=1;
664  const int lb1 = Lb(x);
665  const int lb2 = Lb(y);
666  int rnd = getround();
667 
668  cvector xmid(n),xrad(n),ymid(n),yrad(n);
669  real crad_re = 0.0, crad_im = 0.0;
670  complex res_inf(0.0,0.0),res_sup(0.0,0.0);
671  real tmp1,tmp2,tmp3,tmp4;
672 
673  double* dxmid = xmid.to_blas_array();
674  double* dymid = ymid.to_blas_array();
675 
676  setround(1);
677 
678  for(int i=0 ; i<n ; i++) {
679  xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i]));
680  ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i]));
681  xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]);
682  yrad[i+1] = ymid[i+1] - Inf(y[i+lb2]);
683 
684  tmp1 = abs(Re(ymid[i+1])) + Re(yrad[i+1]);
685  tmp2 = abs(Im(ymid[i+1])) + Im(yrad[i+1]);
686  tmp3 = abs(Re(xmid[i+1]));
687  tmp4 = abs(Im(xmid[i+1]));
688 
689  crad_re += Re(xrad[i+1]) * tmp1 + tmp3 * Re(yrad[i+1]);
690  crad_re += Im(xrad[i+1]) * tmp2 + tmp4 * Im(yrad[i+1]);
691 
692  crad_im += Re(xrad[i+1]) * tmp2 + tmp3 * Im(yrad[i+1]);
693  crad_im += Im(xrad[i+1]) * tmp1 + tmp4 * Re(yrad[i+1]);
694  }
695 
696  cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_sup);
697 
698  res_sup += complex(crad_re,crad_im);
699 
700  setround(-1);
701 
702  cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_inf);
703 
704  res_inf -= complex(crad_re,crad_im);
705 
706  setround(rnd);
707 
708  UncheckedSetInf(res,res_inf);
709  UncheckedSetSup(res,res_sup);
710 
711 }
712 
713 
714 inline void blasdot(const civector& x, const cvector& y, cinterval& res) {
715  const int n = VecLen(x);
716  const int inc=1;
717  const int lb1 = Lb(x);
718  int rnd = getround();
719 
720  cvector xmid(n),xrad(n);
721  real crad_re = 0.0, crad_im = 0.0;
722  complex res_inf(0.0,0.0),res_sup(0.0,0.0);
723  real tmp1,tmp2,tmp3,tmp4;
724 
725  double* dxmid = xmid.to_blas_array();
726  double* dy = y.to_blas_array();
727 
728  setround(1);
729 
730  for(int i=0 ; i<n ; i++) {
731  xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i]));
732  xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]);
733 
734  tmp1 = abs(Re(y[i+1]));
735  tmp2 = abs(Im(y[i+1]));
736 
737  crad_re += Re(xrad[i+1]) * tmp1;
738  crad_re += Im(xrad[i+1]) * tmp2;
739 
740  crad_im += Re(xrad[i+1]) * tmp2;
741  crad_im += Im(xrad[i+1]) * tmp1;
742  }
743 
744  cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_sup);
745 
746  res_sup += complex(crad_re,crad_im);
747 
748  setround(-1);
749 
750  cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_inf);
751 
752  res_inf -= complex(crad_re,crad_im);
753 
754  setround(rnd);
755 
756  UncheckedSetInf(res,res_inf);
757  UncheckedSetSup(res,res_sup);
758 
759 }
760 
761 inline void blasdot(const cvector& x, const civector& y, cinterval& res) {
762  const int n = VecLen(x);
763  const int inc=1;
764  const int lb1 = Lb(x);
765  int rnd = getround();
766 
767  cvector ymid(n),yrad(n);
768  real crad_re = 0.0, crad_im = 0.0;
769  complex res_inf(0.0,0.0),res_sup(0.0,0.0);
770  real tmp1,tmp2,tmp3,tmp4;
771 
772  double* dx = x.to_blas_array();
773  double* dymid = ymid.to_blas_array();
774 
775  setround(1);
776 
777  for(int i=0 ; i<n ; i++) {
778  ymid[i+1] = Inf(y[lb1+i]) + 0.5*(Sup(y[lb1+i])-Inf(y[lb1+i]));
779  yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]);
780 
781  tmp1 = abs(Re(x[i+1]));
782  tmp2 = abs(Im(x[i+1]));
783 
784  crad_re += Re(yrad[i+1]) * tmp1;
785  crad_re += Im(yrad[i+1]) * tmp2;
786 
787  crad_im += Re(yrad[i+1]) * tmp2;
788  crad_im += Im(yrad[i+1]) * tmp1;
789  }
790 
791  cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_sup);
792 
793  res_sup += complex(crad_re,crad_im);
794 
795  setround(-1);
796 
797  cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_inf);
798 
799  res_inf -= complex(crad_re,crad_im);
800 
801  setround(rnd);
802 
803  UncheckedSetInf(res,res_inf);
804  UncheckedSetSup(res,res_sup);
805 
806 }
807 
808 
809 inline void blasdot(const civector& x, const ivector& y, cinterval& res) {
810  const int n = VecLen(x);
811  const int inc=1, inc2=2;
812  const int lb1 = Lb(x);
813  const int lb2 = Lb(y);
814  int rnd = getround();
815 
816  cvector xmid(n),xrad(n);
817  rvector ymid(n),yrad(n);
818  real crad_re = 0.0, crad_im = 0.0;
819  complex res_inf(0.0,0.0),res_sup(0.0,0.0);
820  real tmp1,tmp2,tmp3,tmp4;
821 
822  double* dxmid = xmid.to_blas_array();
823  double* dymid = ymid.to_blas_array();
824 
825  setround(1);
826 
827  for(int i=0 ; i<n ; i++) {
828  xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i]));
829  ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i]));
830  xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]);
831  yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]);
832 
833  tmp1 = abs(ymid[i+1]) + yrad[i+1];
834  tmp2 = abs(ymid[i+1]) + yrad[i+1];
835  tmp3 = abs(Re(xmid[i+1]));
836  tmp4 = abs(Im(xmid[i+1]));
837 
838  crad_re += Re(xrad[i+1]) * tmp1 + tmp3 * yrad[i+1];
839  crad_re += Im(xrad[i+1]) * tmp2 + tmp4 * yrad[i+1];
840 
841  crad_im += Re(xrad[i+1]) * tmp2 + tmp3 * yrad[i+1];
842  crad_im += Im(xrad[i+1]) * tmp1 + tmp4 * yrad[i+1];
843  }
844 
845  SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc));
846  SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc));
847 
848  res_sup += complex(crad_re,crad_im);
849 
850  setround(-1);
851 
852  SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc));
853  SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc));
854 
855  res_inf -= complex(crad_re,crad_im);
856 
857  setround(rnd);
858 
859  UncheckedSetInf(res,res_inf);
860  UncheckedSetSup(res,res_sup);
861 }
862 
863 
864 inline void blasdot(const ivector& x, const civector& y, cinterval& res) {
865  const int n = VecLen(x);
866  const int inc=1, inc2=2;
867  const int lb1 = Lb(x);
868  const int lb2 = Lb(y);
869  int rnd = getround();
870 
871  rvector xmid(n),xrad(n);
872  cvector ymid(n),yrad(n);
873  real crad_re = 0.0, crad_im = 0.0;
874  complex res_inf(0.0,0.0),res_sup(0.0,0.0);
875  real tmp1,tmp2,tmp3,tmp4;
876 
877  double* dxmid = xmid.to_blas_array();
878  double* dymid = ymid.to_blas_array();
879 
880  setround(1);
881 
882  for(int i=0 ; i<n ; i++) {
883  xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i]));
884  ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i]));
885  xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]);
886  yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]);
887 
888  tmp1 = abs(Re(ymid[i+1])) + Re(yrad[i+1]);
889  tmp2 = abs(Im(ymid[i+1])) + Im(yrad[i+1]);
890  tmp3 = abs(xmid[i+1]);
891  tmp4 = abs(xmid[i+1]);
892 
893  crad_re += xrad[i+1] * tmp1 + tmp3 * Re(yrad[i+1]);
894  crad_re += xrad[i+1] * tmp2 + tmp4 * Im(yrad[i+1]);
895 
896  crad_im += xrad[i+1] * tmp2 + tmp3 * Im(yrad[i+1]);
897  crad_im += xrad[i+1] * tmp1 + tmp4 * Re(yrad[i+1]);
898  }
899 
900  SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2));
901  SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2));
902 
903  res_sup += complex(crad_re,crad_im);
904 
905  setround(-1);
906 
907  SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2));
908  SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2));
909 
910  res_inf -= complex(crad_re,crad_im);
911 
912  setround(rnd);
913 
914  UncheckedSetInf(res,res_inf);
915  UncheckedSetSup(res,res_sup);
916 }
917 #endif
918 #endif
919 /***************************************************************************/
920 
921 
922 #ifdef _CXSC_BLAS_RMATRIX
923 #ifndef _CXSC_BLAS_RMATVEC_INC
924 #define _CXSC_BLAS_RMATVEC_INC
925 inline void blasmvmul(const rmatrix& A, const rvector& x, rvector& r) {
926 
927  double* DA = (double*) A.to_blas_array();
928  double* dx = (double*) x.to_blas_array();
929  double* dr = (double*) r.to_blas_array();
930 
931  const int m = Ub(A,1) - Lb(A,1) + 1;
932  const int n = Ub(A,2) - Lb(A,2) + 1;
933  const double alpha = 1.0, beta = 0.0;
934  const int inc = 1;
935 
936  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
937  alpha, DA, n, dx, inc, beta, dr, inc);
938 }
939 #endif
940 #endif
941 
942 #ifdef _CXSC_BLAS_RMATRIX
943 #ifdef _CXSC_BLAS_IVECTOR
944 #ifndef _CXSC_BLAS_RIMATVEC_INC
945 #define _CXSC_BLAS_RIMATVEC_INC
946 inline void blasmvmul(const rmatrix& A, const rvector& x, ivector& r) {
947  int rnd = getround();
948  double* DA = (double*) A.to_blas_array();
949  double* dx = (double*) x.to_blas_array();
950 
951  const int m = Ub(A,1) - Lb(A,1) + 1;
952  const int n = Ub(A,2) - Lb(A,2) + 1;
953  const double alpha = 1.0, beta = 0.0;
954  const int inc = 1;
955  double* dr = new double[n];
956 
957  setround(1);
958  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
959  alpha, DA, n, dx, inc, beta, dr, inc);
960  for(int i=0 ; i<n ; i++)
961  UncheckedSetSup(r[i+Lb(r)], dr[i]);
962 
963  setround(-1);
964  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
965  alpha, DA, n, dx, inc, beta, dr, inc);
966  for(int i=0 ; i<n ; i++)
967  UncheckedSetInf(r[i+Lb(r)], dr[i]);
968 
969  delete[] dr;
970  setround(rnd);
971 }
972 
973 inline void blasmvmul(const rmatrix& A, const ivector& x, ivector& r) {
974  const int n = VecLen(x);
975  const int m = Ub(A,1)-Lb(A,1)+1;
976  const int inc=1;
977  const int lb1 = Lb(A,2);
978  const int lb2 = Lb(x);
979  int rnd = getround();
980 
981  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
982 
983  double* dxi = x_inf.to_blas_array();
984  double* dxs = x_sup.to_blas_array();
985  double* dyi = y_inf.to_blas_array();
986  double* dys = y_sup.to_blas_array();
987  double res_inf,res_sup;
988 
989  setround(1);
990 
991  for(int i=1 ; i<=m ; i++) {
992  bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
993  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
994  UncheckedSetSup(r[i], res_sup);
995  // setround(-1);
996  res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
997  UncheckedSetInf(r[i], -res_inf);
998  }
999 
1000  setround(rnd);
1001 }
1002 #endif
1003 #endif
1004 #endif
1005 
1006 #ifdef _CXSC_BLAS_IMATRIX
1007 #ifdef _CXSC_BLAS_IVECTOR
1008 #ifndef _CXSC_BLAS_IRMATVEC_INC
1009 #define _CXSC_BLAS_IRMATVEC_INC
1010 inline void blasmvmul(const imatrix& A, const rvector& x, ivector& r) {
1011  const int n = VecLen(x);
1012  const int m = Ub(A,1)-Lb(A,1)+1;
1013  const int inc=1;
1014  const int lb1 = Lb(A,2);
1015  const int lb2 = Lb(x);
1016  int rnd = getround();
1017 
1018  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
1019 
1020  double* dxi = x_inf.to_blas_array();
1021  double* dxs = x_sup.to_blas_array();
1022  double* dyi = y_inf.to_blas_array();
1023  double* dys = y_sup.to_blas_array();
1024  double res_inf,res_sup;
1025 
1026  setround(1);
1027 
1028  for(int i=1 ; i<=m ; i++) {
1029  bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
1030  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
1031  UncheckedSetSup(r[i], res_sup);
1032  // setround(-1);
1033  res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
1034  UncheckedSetInf(r[i], -res_inf);
1035  }
1036 
1037  setround(rnd);
1038 
1039 }
1040 #endif
1041 #endif
1042 #endif
1043 
1044 #ifdef _CXSC_BLAS_IMATRIX
1045 #ifdef _CXSC_BLAS_IVECTOR
1046 #ifndef _CXSC_BLAS_IIMATVEC_INC
1047 #define _CXSC_BLAS_IIMATVEC_INC
1048 inline void blasmvmul(const imatrix& A, const ivector& x, ivector& r) {
1049  const int n = VecLen(x);
1050  const int m = Ub(A,1)-Lb(A,1)+1;
1051  const int inc=1;
1052  const int lb1 = Lb(A,2);
1053  const int lb2 = Lb(x);
1054  int rnd = getround();
1055 
1056  rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
1057 
1058  double* dxi = x_inf.to_blas_array();
1059  double* dxs = x_sup.to_blas_array();
1060  double* dyi = y_inf.to_blas_array();
1061  double* dys = y_sup.to_blas_array();
1062  double res_inf,res_sup;
1063 
1064  setround(1);
1065 
1066  for(int i=1 ; i<=m ; i++) {
1067  bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
1068  res_sup = cblas_ddot(n, dxs, inc, dys, inc);
1069  UncheckedSetSup(r[i], res_sup);
1070  // setround(-1);
1071  res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
1072  UncheckedSetInf(r[i], -res_inf);
1073  }
1074 
1075  setround(rnd);
1076 
1077 }
1078 #endif
1079 #endif
1080 #endif
1081 
1082 #ifdef _CXSC_BLAS_CMATRIX
1083 #ifdef _CXSC_BLAS_CIVECTOR
1084 #ifndef _CXSC_BLAS_CIMATVEC_INC
1085 #define _CXSC_BLAS_CIMATVEC_INC
1086 inline void blasmvmul(const cmatrix& A, const ivector& x, civector& r) {
1087  const int m = Ub(A,1)-Lb(A,1)+1;
1088 
1089  for(int i=0 ; i<m ; i++)
1090  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1091 }
1092 #endif
1093 #endif
1094 #endif
1095 
1096 #ifdef _CXSC_BLAS_IMATRIX
1097 #ifdef _CXSC_BLAS_CIVECTOR
1098 #ifndef _CXSC_BLAS_ICMATVEC_INC
1099 #define _CXSC_BLAS_ICMATVEC_INC
1100 inline void blasmvmul(const imatrix& A, const cvector& x, civector& r) {
1101  const int m = Ub(A,1)-Lb(A,1)+1;
1102 
1103  for(int i=0 ; i<m ; i++)
1104  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1105 }
1106 #endif
1107 #endif
1108 #endif
1109 
1110 #ifdef _CXSC_BLAS_RMATRIX
1111 #ifdef _CXSC_BLAS_CIVECTOR
1112 #ifndef _CXSC_BLAS_RCIMATVEC_INC
1113 #define _CXSC_BLAS_RCIMATVEC_INC
1114 inline void blasmvmul(const rmatrix& A, const civector& x, civector& r) {
1115  const int m = Ub(A,1)-Lb(A,1)+1;
1116 
1117  for(int i=0 ; i<m ; i++)
1118  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1119 }
1120 #endif
1121 #endif
1122 #endif
1123 
1124 #ifdef _CXSC_BLAS_CIMATRIX
1125 #ifdef _CXSC_BLAS_CIVECTOR
1126 #ifndef _CXSC_BLAS_CIRMATVEC_INC
1127 #define _CXSC_BLAS_CIRMATVEC_INC
1128 inline void blasmvmul(const cimatrix& A, const rvector& x, civector& r) {
1129  const int m = Ub(A,1)-Lb(A,1)+1;
1130 
1131  for(int i=0 ; i<m ; i++)
1132  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1133 }
1134 #endif
1135 #endif
1136 #endif
1137 
1138 #ifdef _CXSC_BLAS_CMATRIX
1139 #ifdef _CXSC_BLAS_CIVECTOR
1140 #ifndef _CXSC_BLAS_CCIMATVEC_INC
1141 #define _CXSC_BLAS_CCIMATVEC_INC
1142 inline void blasmvmul(const cmatrix& A, const civector& x, civector& r) {
1143  const int m = Ub(A,1)-Lb(A,1)+1;
1144 
1145  for(int i=0 ; i<m ; i++)
1146  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1147 }
1148 
1149 inline void blasmvmul(const cmatrix& A, const cvector& x, civector& r) {
1150  ivector tmp1, tmp2;
1151  blasmvmul(Re(A),Re(x),tmp1);
1152  blasmvmul(Im(A),Im(x),tmp2);
1153  SetRe(r,tmp1-tmp2);
1154  blasmvmul(Re(A),Im(x),tmp1);
1155  blasmvmul(Im(A),Re(x),tmp2);
1156  SetIm(r,tmp1+tmp2);
1157 }
1158 
1159 inline void blasmvmul(const cmatrix& A, const rvector& x, civector& r) {
1160  int n = VecLen(r);
1161  ivector re(n),im(n);
1162 
1163  blasmvmul(Re(A),x,re);
1164  blasmvmul(Im(A),x,im);
1165 
1166  SetRe(r,re);
1167  SetIm(r,im);
1168 }
1169 
1170 inline void blasmvmul(const rmatrix& A, const cvector& x, civector& r) {
1171  int n = VecLen(r);
1172  ivector re(n),im(n);
1173 
1174  blasmvmul(A,Re(x),re);
1175  blasmvmul(A,Im(x),im);
1176 
1177  SetRe(r,re);
1178  SetIm(r,im);
1179 }
1180 #endif
1181 #endif
1182 #endif
1183 
1184 #ifdef _CXSC_BLAS_CIMATRIX
1185 #ifdef _CXSC_BLAS_CIVECTOR
1186 #ifndef _CXSC_BLAS_CICMATVEC_INC
1187 #define _CXSC_BLAS_CICMATVEC_INC
1188 inline void blasmvmul(const cimatrix& A, const cvector& x, civector& r) {
1189  const int m = Ub(A,1)-Lb(A,1)+1;
1190 
1191  for(int i=0 ; i<m ; i++)
1192  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1193 }
1194 #endif
1195 #endif
1196 #endif
1197 
1198 #ifdef _CXSC_BLAS_IMATRIX
1199 #ifdef _CXSC_BLAS_CIVECTOR
1200 #ifndef _CXSC_BLAS_ICIMATVEC_INC
1201 #define _CXSC_BLAS_ICIMATVEC_INC
1202 inline void blasmvmul(const imatrix& A, const civector& x, civector& r) {
1203  const int m = Ub(A,1)-Lb(A,1)+1;
1204 
1205  for(int i=0 ; i<m ; i++)
1206  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1207 }
1208 #endif
1209 #endif
1210 #endif
1211 
1212 #ifdef _CXSC_BLAS_CIMATRIX
1213 #ifdef _CXSC_BLAS_CIVECTOR
1214 #ifndef _CXSC_BLAS_CIIMATVEC_INC
1215 #define _CXSC_BLAS_CIIMATVEC_INC
1216 inline void blasmvmul(const cimatrix& A, const ivector& x, civector& r) {
1217  const int m = Ub(A,1)-Lb(A,1)+1;
1218 
1219  for(int i=0 ; i<m ; i++)
1220  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1221 }
1222 #endif
1223 #endif
1224 #endif
1225 
1226 #ifdef _CXSC_BLAS_CIMATRIX
1227 #ifdef _CXSC_BLAS_CIVECTOR
1228 #ifndef _CXSC_BLAS_CICIMATVEC_INC
1229 #define _CXSC_BLAS_CICIMATVEC_INC
1230 inline void blasmvmul(const cimatrix& A, const civector& x, civector& r) {
1231  const int m = Ub(A,1)-Lb(A,1)+1;
1232 
1233  for(int i=0 ; i<m ; i++)
1234  blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]);
1235 }
1236 #endif
1237 #endif
1238 #endif
1239 
1240 #ifdef _CXSC_BLAS_CMATRIX
1241 #ifdef _CXSC_BLAS_CVECTOR
1242 #ifndef _CXSC_BLAS_CCMATVEC_INC
1243 #define _CXSC_BLAS_CCMATVEC_INC
1244 inline void blasmvmul(const cmatrix& A, const cvector& x, cvector& r) {
1245 
1246  double* DA = (double*) A.to_blas_array();
1247  double* dx = (double*) x.to_blas_array();
1248  double* dr = (double*) r.to_blas_array();
1249 
1250  const int m = Ub(A,1) - Lb(A,1) + 1;
1251  const int n = Ub(A,2) - Lb(A,2) + 1;
1252  const complex alpha(1.0,0.0);
1253  const complex beta(0.0,0.0);
1254  const int inc = 1;
1255 
1256  cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n,
1257  &alpha, DA, n, dx, inc, &beta, dr, inc);
1258 }
1259 #endif
1260 #endif
1261 #endif
1262 
1263 #ifdef _CXSC_BLAS_RMATRIX
1264 #ifdef _CXSC_BLAS_CVECTOR
1265 #ifndef _CXSC_BLAS_RCMATVEC_INC
1266 #define _CXSC_BLAS_RCMATVEC_INC
1267 inline void blasmvmul(const rmatrix& A, const cvector& x, cvector& r) {
1268 
1269  double* DA = (double*) A.to_blas_array();
1270  double* dx = (double*) x.to_blas_array();
1271  double* dr = (double*) r.to_blas_array();
1272 
1273  const int m = Ub(A,1) - Lb(A,1) + 1;
1274  const int n = Ub(A,2) - Lb(A,2) + 1;
1275  const double alpha = 1.0;
1276  const double beta = 0.0;
1277  const int inc = 2;
1278 
1279  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
1280  alpha, DA, m, dx, inc, beta, dr, inc);
1281 
1282  cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
1283  alpha, DA, n, dx+1, inc, beta, dr+1, inc);
1284 }
1285 #endif
1286 #endif
1287 #endif
1288 
1289 #ifdef _CXSC_BLAS_CMATRIX
1290 #ifdef _CXSC_BLAS_CVECTOR
1291 #ifndef _CXSC_BLAS_CRMATVEC_INC
1292 #define _CXSC_BLAS_CRMATVEC_INC
1293 inline void blasmvmul(const cmatrix& A, const rvector& x, cvector& r) {
1294  cvector tmp(x);
1295 
1296  double* DA = (double*) A.to_blas_array();
1297  double* dx = (double*) tmp.to_blas_array();
1298  double* dr = (double*) r.to_blas_array();
1299 
1300  const int m = Ub(A,1) - Lb(A,1) + 1;
1301  const int n = Ub(A,2) - Lb(A,2) + 1;
1302  const complex alpha(1.0,0.0);
1303  const complex beta(0.0,0.0);
1304  const int inc = 1;
1305 
1306  cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n,
1307  &alpha, DA, n, dx, inc, &beta, dr, inc);
1308 }
1309 #endif
1310 #endif
1311 #endif
1312 
1313 /***************************************************************************/
1314 
1315 #ifdef _CXSC_BLAS_RMATRIX
1316 #ifndef _CXSC_BLAS_RMATRIX_INC
1317 #define _CXSC_BLAS_RMATRIX_INC
1318 inline void blasmatmul(const rmatrix &A, const rmatrix &B, rmatrix &C) {
1319 
1320  double* DA = (double*) A.to_blas_array();
1321  double* DB = (double*) B.to_blas_array();
1322  double* DC = (double*) C.to_blas_array();
1323 
1324  const int m = Ub(A,1) - Lb(A,1) + 1;
1325  const int n = Ub(A,2) - Lb(A,2) + 1;
1326  const int o = Ub(C,2) - Lb(C,2) + 1;
1327  const double alpha = 1.0, beta = 0.0;
1328 
1329  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1330  n, alpha, DA, n, DB, o, beta, DC, o);
1331 }
1332 #endif
1333 #endif
1334 
1335 #ifdef _CXSC_BLAS_IMATRIX
1336 #ifndef _CXSC_BLAS_IMATRIX_INC
1337 #define _CXSC_BLAS_IMATRIX_INC
1338 inline void blasmatmul(const rmatrix &A, const rmatrix &B, imatrix &C) {
1339  int rnd = getround();
1340  double* DA = (double*) A.to_blas_array();
1341  double* DB = (double*) B.to_blas_array();
1342  int ind;
1343 
1344  const int m = Ub(A,1) - Lb(A,1) + 1;
1345  const int n = Ub(A,2) - Lb(A,2) + 1;
1346  const int o = Ub(C,2) - Lb(C,2) + 1;
1347  double* DC = new double[m*o];
1348  const double alpha = 1.0, beta = 0.0;
1349 
1350  setround(1);
1351  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1352  n, alpha, DA, n, DB, o, beta, DC, o);
1353 
1354  for(int i=0 ; i<m ; i++) {
1355  for(int j=0 ; j<o ; j++) {
1356  ind = i*o+j;
1357  UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], DC[ind]);
1358  }
1359  }
1360 
1361  setround(-1);
1362  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1363  n, alpha, DA, n, DB, o, beta, DC, o);
1364 
1365  for(int i=0 ; i<m ; i++) {
1366  for(int j=0 ; j<o ; j++) {
1367  ind = i*o+j;
1368  UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], DC[ind]);
1369  }
1370  }
1371 
1372  delete[] DC;
1373 
1374  setround(rnd);
1375 }
1376 
1377 
1378 inline void blasmatmul(const imatrix &A, const imatrix &B, imatrix &C) {
1379  const int m = Ub(A,1) - Lb(A,1) + 1;
1380  const int n = Ub(A,2) - Lb(A,2) + 1;
1381  const int o = Ub(C,2) - Lb(C,2) + 1;
1382 
1383  int lbC_1 = Lb(C,1); int lbC_2 = Lb(C,2);
1384  Resize(C);
1385 
1386  double* Amid = new double[m*n];
1387  double* Aabs = new double[m*n];
1388  double* Arad = new double[m*n];
1389  double* Bmid = new double[n*o];
1390  double* Babs = new double[n*o];
1391  double* Brad = new double[n*o];
1392  double* C1 = new double[m*o];
1393  double* C2 = new double[m*o];
1394  int rnd = getround();
1395 
1396  const double alpha = 1.0, beta = 0.0;
1397 
1398  setround(1);
1399 
1400  int ind;
1401 
1402  //Compute mid and rad of A
1403  for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
1404  for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
1405  ind = (i-Lb(A,1))*n+(j-Lb(A,2));
1406  Amid[ind] = _double(Inf(A[i][j]) + 0.5*(Sup(A[i][j]) - Inf(A[i][j])));
1407  Arad[ind] = _double(Amid[ind] - Inf(A[i][j]));
1408  Aabs[ind] = fabs(Amid[ind]);
1409  }
1410  }
1411 
1412  //Compute mid and rad of B
1413  for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) {
1414  for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) {
1415  ind = (i-Lb(B,1))*o+(j-Lb(B,2));
1416  Bmid[ind] = _double(Inf(B[i][j]) + 0.5*(Sup(B[i][j]) - Inf(B[i][j])));
1417  Brad[ind] = _double(Bmid[ind] - Inf(B[i][j]));
1418  Babs[ind] = fabs(Bmid[ind]);
1419  }
1420  }
1421 
1422  setround(-1);
1423 
1424  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1425  n, alpha, Amid, n, Bmid, o, beta, C1, o);
1426 
1427  setround(1);
1428 
1429  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1430  n, alpha, Amid, n, Bmid, o, beta, C2, o);
1431 
1432  delete[] Amid;
1433  delete[] Bmid;
1434 
1435  double* Cmid = new double[m*o];
1436  double* Crad = new double[m*o];
1437 
1438  for(int i=0 ; i<m ; i++) {
1439  for(int j=0 ; j<o ; j++) {
1440  ind = i*o+j;
1441  Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1442  Crad[ind] = Cmid[ind] - C1[ind];
1443  }
1444  }
1445 
1446  for(int i=0 ; i<n ; i++) {
1447  for(int j=0 ; j<o ; j++) {
1448  ind = i*o+j;
1449  Babs[ind] += Brad[ind];
1450  }
1451  }
1452 
1453  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1454  n, alpha, Arad, n, Babs, o, beta, C1, o);
1455 
1456  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1457  n, alpha, Aabs, n, Brad, o, beta, C2, o);
1458 
1459  delete[] Aabs;
1460  delete[] Arad;
1461  delete[] Babs;
1462  delete[] Brad;
1463 
1464  Resize(C, lbC_1, lbC_1+m-1, lbC_2, lbC_2+o-1);
1465  //C = imatrix(lbC_1, lbC_1+m-1, lbC_2, lbC_2+o-1);
1466 
1467  for(int i=0 ; i<m ; i++) {
1468  for(int j=0 ; j<o ; j++) {
1469  ind = i*o+j;
1470  Crad[ind] += C1[ind] + C2[ind];
1471  UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]);
1472  }
1473  }
1474 
1475  setround(-1);
1476 
1477  for(int i=0 ; i<m ; i++) {
1478  for(int j=0 ; j<o ; j++) {
1479  ind = i*o+j;
1480  UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]);
1481  }
1482  }
1483 
1484  setround(rnd);
1485 
1486  delete[] C1;
1487  delete[] C2;
1488  delete[] Crad;
1489  delete[] Cmid;
1490 
1491 }
1492 
1493 
1494 inline void blasmatmul(const rmatrix &A, const imatrix &B, imatrix &C) {
1495 
1496  const int m = Ub(A,1) - Lb(A,1) + 1;
1497  const int n = Ub(A,2) - Lb(A,2) + 1;
1498  const int o = Ub(B,2) - Lb(B,2) + 1;
1499 
1500  int lb1_C = Lb(C,1);
1501  int lb2_C = Lb(C,2);
1502  Resize(C);
1503 
1504  double* DA = (double*)A.to_blas_array();
1505  double* Bmid = new double[n*o];
1506  double* Brad = new double[n*o];
1507  double* C1 = new double[m*o];
1508  double* C2 = new double[m*o];
1509  int rnd = getround();
1510 
1511  const double alpha = 1.0;
1512  double beta = 0.0;
1513 
1514  setround(1);
1515 
1516  int ind;
1517 
1518  //Compute mid and rad of B
1519  for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) {
1520  for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) {
1521  ind = (i-Lb(B,1))*o+(j-Lb(B,2));
1522  Bmid[ind] = _double(Inf(B[i][j]) + 0.5*(Sup(B[i][j]) - Inf(B[i][j])));
1523  Brad[ind] = _double(Bmid[ind] - Inf(B[i][j]));
1524  }
1525  }
1526 
1527  setround(-1);
1528 
1529  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1530  n, alpha, DA, n, Bmid, o, beta, C1, o);
1531 
1532  setround(1);
1533 
1534  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1535  n, alpha, DA, n, Bmid, o, beta, C2, o);
1536 
1537  delete[] Bmid;
1538 
1539  double* Cmid = new double[m*o];
1540  for(int i=0 ; i<m ; i++) {
1541  for(int j=0 ; j<o ; j++) {
1542  ind = i*o+j;
1543  Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1544  }
1545  }
1546 
1547  delete[] C2;
1548 
1549  double* Crad = new double[m*o];
1550  for(int i=0 ; i<m ; i++) {
1551  for(int j=0 ; j<o ; j++) {
1552  ind = i*o+j;
1553  Crad[ind] = Cmid[ind] - C1[ind];
1554  }
1555  }
1556 
1557  delete[] C1;
1558 
1559  double* Aabs = new double[m*n];
1560  //Compute abs(A)
1561  for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
1562  for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
1563  ind = (i-Lb(A,1))*n+(j-Lb(A,2));
1564  Aabs[ind] = fabs(DA[ind]);
1565  }
1566  }
1567 
1568  beta = 1.0;
1569 
1570  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1571  n, alpha, Aabs, n, Brad, o, beta, Crad, o);
1572 
1573  delete[] Aabs;
1574  delete[] Brad;
1575 
1576  Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1577  //C = imatrix(lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1578 
1579  for(int i=0 ; i<m ; i++) {
1580  for(int j=0 ; j<o ; j++) {
1581  ind = i*o+j;
1582  UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]);
1583  }
1584  }
1585 
1586  setround(-1);
1587 
1588  for(int i=0 ; i<m ; i++) {
1589  for(int j=0 ; j<o ; j++) {
1590  ind = i*o+j;
1591  UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]);
1592  }
1593  }
1594 
1595  setround(rnd);
1596 
1597  delete[] Crad;
1598  delete[] Cmid;
1599 }
1600 
1601 
1602 inline void blasmatmul(const imatrix &A, const rmatrix &B, imatrix &C) {
1603 
1604  const int m = Ub(A,1) - Lb(A,1) + 1;
1605  const int n = Ub(A,2) - Lb(A,2) + 1;
1606  const int o = Ub(C,2) - Lb(C,2) + 1;
1607 
1608  int lb1_C = Lb(C,1);
1609  int lb2_C = Lb(C,2);
1610  Resize(C);
1611 
1612  double* DB = (double*) B.to_blas_array();
1613  double* Amid = new double[m*n];
1614  double* Arad = new double[m*n];
1615  double* C1 = new double[m*o];
1616  double* C2 = new double[m*o];
1617  int rnd = getround();
1618 
1619  const double alpha = 1.0;
1620  double beta = 0.0;
1621 
1622  setround(1);
1623 
1624  int ind;
1625 
1626  //Compute mid and rad of A
1627  for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) {
1628  for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) {
1629  ind = (i-Lb(A,1))*n+(j-Lb(A,2));
1630  Amid[ind] = _double(Inf(A[i][j]) + 0.5*(Sup(A[i][j]) - Inf(A[i][j])));
1631  Arad[ind] = _double(Amid[ind] - Inf(A[i][j]));
1632  }
1633  }
1634 
1635  setround(-1);
1636 
1637  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1638  n, alpha, Amid, n, DB, o, beta, C1, o);
1639 
1640  setround(1);
1641 
1642  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1643  n, alpha, Amid, n, DB, o, beta, C2, o);
1644 
1645  delete[] Amid;
1646 
1647  double* Cmid = new double[m*o];
1648  for(int i=0 ; i<m ; i++) {
1649  for(int j=0 ; j<o ; j++) {
1650  ind = i*o+j;
1651  Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1652  }
1653  }
1654 
1655  delete[] C2;
1656 
1657  double* Crad = new double[m*o];
1658  for(int i=0 ; i<m ; i++) {
1659  for(int j=0 ; j<o ; j++) {
1660  ind = i*o+j;
1661  Crad[ind] = Cmid[ind] - C1[ind];
1662  }
1663  }
1664 
1665  delete[] C1;
1666 
1667  double* Babs = new double[n*o];
1668  //Compute abs of B
1669  for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) {
1670  for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) {
1671  ind = (i-Lb(B,1))*o+(j-Lb(B,2));
1672  Babs[ind] = fabs(DB[ind]);
1673  }
1674  }
1675 
1676  beta = 1.0;
1677 
1678  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1679  n, alpha, Arad, n, Babs, o, beta, Crad, o);
1680 
1681  delete[] Arad;
1682  delete[] Babs;
1683 
1684  Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1685  //C = imatrix(lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1686 
1687  for(int i=0 ; i<m ; i++) {
1688  for(int j=0 ; j<o ; j++) {
1689  ind = i*o+j;
1690  UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]);
1691  }
1692  }
1693 
1694  setround(-1);
1695 
1696  for(int i=0 ; i<m ; i++) {
1697  for(int j=0 ; j<o ; j++) {
1698  ind = i*o+j;
1699  UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]);
1700  }
1701  }
1702 
1703  setround(rnd);
1704 
1705  delete[] Crad;
1706  delete[] Cmid;
1707 }
1708 #endif
1709 #endif
1710 
1711 #ifdef _CXSC_BLAS_CMATRIX
1712 #ifndef _CXSC_BLAS_CMATRIX_INC
1713 #define _CXSC_BLAS_CMATRIX_INC
1714 inline void blasmatmul(const cmatrix &A, const cmatrix &B, cmatrix &C) {
1715 
1716  double* DA = (double*) A.to_blas_array();
1717  double* DB = (double*) B.to_blas_array();
1718  double* DC = (double*) C.to_blas_array();
1719 
1720  const int m = Ub(A,1) - Lb(A,1) + 1;
1721  const int n = Ub(A,2) - Lb(A,2) + 1;
1722  const int o = Ub(C,2) - Lb(C,2) + 1;
1723  const complex alpha(1.0,0.0);
1724  const complex beta(0.0,0.0);
1725 
1726  cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1727  n, &alpha, DA, n, DB, o, &beta, DC, o);
1728 }
1729 
1730 
1731 inline void blasmatmul(const cmatrix &A, const rmatrix &B, cmatrix &C) {
1732  cmatrix tmp(B);
1733 
1734  double* DA = (double*) A.to_blas_array();
1735  double* DB = (double*) tmp.to_blas_array();
1736  double* DC = (double*) C.to_blas_array();
1737 
1738  const int m = Ub(A,1) - Lb(A,1) + 1;
1739  const int n = Ub(A,2) - Lb(A,2) + 1;
1740  const int o = Ub(C,2) - Lb(C,2) + 1;
1741  const complex alpha(1.0,0.0);
1742  const complex beta(0.0,0.0);
1743 
1744  cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1745  n, &alpha, DA, n, DB, o, &beta, DC, o);
1746 }
1747 
1748 inline void blasmatmul(const rmatrix &A, const cmatrix &B, cmatrix &C) {
1749  cmatrix tmp(A);
1750 
1751  double* DA = (double*) tmp.to_blas_array();
1752  double* DB = (double*) B.to_blas_array();
1753  double* DC = (double*) C.to_blas_array();
1754 
1755  const int m = Ub(A,1) - Lb(A,1) + 1;
1756  const int n = Ub(A,2) - Lb(A,2) + 1;
1757  const int o = Ub(C,2) - Lb(C,2) + 1;
1758  const complex alpha(1.0,0.0);
1759  const complex beta(0.0,0.0);
1760 
1761  cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1762  n, &alpha, DA, n, DB, o, &beta, DC, o);
1763 }
1764 #endif
1765 #endif
1766 
1767 #ifdef _CXSC_BLAS_CIMATRIX
1768 #ifndef _CXSC_BLAS_CIMATRIX_INC
1769 #define _CXSC_BLAS_CIMATRIX_INC
1770 inline void blasmatmul(const cmatrix &A, const cmatrix &B, cimatrix &C) {
1771  const int m = Ub(A,1) - Lb(A,1) + 1;
1772  const int o = Ub(C,2) - Lb(C,2) + 1;
1773  imatrix tmp1(m,o), tmp2(m,o);
1774 
1775  blasmatmul(Re(A),Re(B),tmp1);
1776  blasmatmul(Im(A),Im(B),tmp2);
1777  SetRe(C,tmp1-tmp2);
1778  blasmatmul(Re(A),Im(B),tmp1);
1779  blasmatmul(Im(A),Re(B),tmp2);
1780  SetIm(C,tmp1+tmp2);
1781 }
1782 
1783 
1784 inline void blasmatmul(const cmatrix &A, const imatrix &B, cimatrix &C) {
1785  const int m = Ub(A,1) - Lb(A,1) + 1;
1786  const int o = Ub(C,2) - Lb(C,2) + 1;
1787 
1788  imatrix T(m,o);
1789 
1790  blasmatmul(Re(A),B,T);
1791 
1792  SetRe(C, T);
1793 
1794  blasmatmul(Im(A),B,T);
1795 
1796  SetIm(C, T);
1797 }
1798 
1799 inline void blasmatmul(const imatrix &A, const cmatrix &B, cimatrix &C) {
1800  const int m = Ub(A,1) - Lb(A,1) + 1;
1801  const int o = Ub(C,2) - Lb(C,2) + 1;
1802 
1803  imatrix T(m,o);
1804 
1805  blasmatmul(A,Re(B),T);
1806 
1807  SetRe(C, T);
1808 
1809  blasmatmul(A,Im(B),T);
1810 
1811  SetIm(C, T);
1812 }
1813 
1814 inline void blasmatmul(const rmatrix &A, const cimatrix &B, cimatrix &C) {
1815  const int m = Ub(A,1) - Lb(A,1) + 1;
1816  const int o = Ub(C,2) - Lb(C,2) + 1;
1817 
1818  imatrix T(m,o);
1819 
1820  blasmatmul(A,Re(B),T);
1821 
1822  SetRe(C, T);
1823 
1824  blasmatmul(A,Im(B),T);
1825 
1826  SetIm(C, T);
1827 }
1828 
1829 inline void blasmatmul(const cimatrix &A, const rmatrix &B, cimatrix &C) {
1830  const int m = Ub(A,1) - Lb(A,1) + 1;
1831  const int o = Ub(C,2) - Lb(C,2) + 1;
1832 
1833  imatrix T(m,o);
1834 
1835  blasmatmul(Re(A),B,T);
1836 
1837  SetRe(C, T);
1838 
1839  blasmatmul(Im(A),B,T);
1840 
1841  SetIm(C, T);
1842 }
1843 
1844 inline void blasmatmul(const imatrix &A, const cimatrix &B, cimatrix &C) {
1845  const int m = Ub(A,1) - Lb(A,1) + 1;
1846  const int o = Ub(C,2) - Lb(C,2) + 1;
1847 
1848  imatrix T(m,o);
1849 
1850  blasmatmul(A,Re(B),T);
1851 
1852  SetRe(C, T);
1853 
1854  blasmatmul(A,Im(B),T);
1855 
1856  SetIm(C, T);
1857 }
1858 
1859 inline void blasmatmul(const cimatrix &A, const imatrix &B, cimatrix &C) {
1860  const int m = Ub(A,1) - Lb(A,1) + 1;
1861  const int o = Ub(C,2) - Lb(C,2) + 1;
1862 
1863  imatrix T(m,o);
1864 
1865  blasmatmul(Re(A),B,T);
1866 
1867  SetRe(C, T);
1868 
1869  blasmatmul(Im(A),B,T);
1870 
1871  SetIm(C, T);
1872 }
1873 
1874 
1875 inline void blasmatmul(const cimatrix &A, const cmatrix &B, cimatrix &C) {
1876 
1877  const int m = Ub(A,1) - Lb(A,1) + 1;
1878  const int o = Ub(C,2) - Lb(C,2) + 1;
1879  int rnd = getround();
1880 
1881  imatrix T(m,o);
1882 
1883  blasmatmul(Re(A),Re(B),T);
1884 
1885  SetRe(C, T);
1886 
1887  blasmatmul(-Im(A),Im(B),T);
1888 
1889  setround(-1);
1890 
1891  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1892  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1893  InfRe(C[i][j]) += Inf(T[i][j]);
1894  }
1895  }
1896 
1897  setround(1);
1898 
1899  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1900  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1901  SupRe(C[i][j]) += Sup(T[i][j]);
1902  }
1903  }
1904 
1905  setround(0);
1906 
1907  blasmatmul(Re(A),Im(B),T);
1908 
1909  SetIm(C, T);
1910 
1911  blasmatmul(Im(A),Re(B),T);
1912 
1913  setround(-1);
1914 
1915  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1916  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1917  InfIm(C[i][j]) += Inf(T[i][j]);
1918  }
1919  }
1920 
1921  setround(1);
1922 
1923  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1924  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1925  SupIm(C[i][j]) += Sup(T[i][j]);
1926  }
1927  }
1928 
1929  setround(rnd);
1930 }
1931 
1932 inline void blasmatmul(const cmatrix &A, const cimatrix &B, cimatrix &C) {
1933  int rnd = getround();
1934 
1935  imatrix T(Lb(C,1),Ub(C,1),Lb(C,2),Ub(C,2));
1936 
1937  blasmatmul(Re(A),Re(B),T);
1938 
1939  SetRe(C, T);
1940 
1941  blasmatmul(-Im(A),Im(B),T);
1942 
1943  setround(-1);
1944 
1945  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1946  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1947  InfRe(C[i][j]) += Inf(T[i][j]);
1948  }
1949  }
1950 
1951  setround(1);
1952 
1953  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1954  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1955  SupRe(C[i][j]) += Sup(T[i][j]);
1956  }
1957  }
1958 
1959  setround(0);
1960 
1961  blasmatmul(Re(A),Im(B),T);
1962 
1963  SetIm(C, T);
1964 
1965  blasmatmul(Im(A),Re(B),T);
1966 
1967  setround(-1);
1968 
1969  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1970  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1971  InfIm(C[i][j]) += Inf(T[i][j]);
1972  }
1973  }
1974 
1975  setround(1);
1976 
1977  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
1978  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
1979  SupIm(C[i][j]) += Sup(T[i][j]);
1980  }
1981  }
1982 
1983  setround(rnd);
1984 }
1985 
1986 
1987 inline void blasmatmul(const cimatrix &A, const cimatrix &B, cimatrix &C) {
1988  const int m = Ub(A,1) - Lb(A,1) + 1;
1989  const int o = Ub(C,2) - Lb(C,2) + 1;
1990  int rnd = getround();
1991 
1992  imatrix T(m,o);
1993 
1994  blasmatmul(Re(A),Re(B),T);
1995 
1996  SetRe(C, T);
1997 
1998  blasmatmul(-Im(A),Im(B),T);
1999 
2000  setround(-1);
2001 
2002  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
2003  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
2004  InfRe(C[i][j]) += Inf(T[i][j]);
2005  }
2006  }
2007 
2008  setround(1);
2009 
2010  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
2011  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
2012  SupRe(C[i][j]) += Sup(T[i][j]);
2013  }
2014  }
2015 
2016  setround(0);
2017 
2018  blasmatmul(Re(A),Im(B),T);
2019 
2020  SetIm(C, T);
2021 
2022  blasmatmul(Im(A),Re(B),T);
2023 
2024  setround(-1);
2025 
2026  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
2027  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
2028  InfIm(C[i][j]) += Inf(T[i][j]);
2029  }
2030  }
2031 
2032  setround(1);
2033 
2034  for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) {
2035  for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) {
2036  SupIm(C[i][j]) += Sup(T[i][j]);
2037  }
2038  }
2039 
2040  setround(rnd);
2041 }
2042 
2043 inline void blasmatmul(const rmatrix &A, const cmatrix &B, cimatrix &C) {
2044  int m = Ub(A,1)-Lb(A,1)+1;
2045  int n = Ub(A,2)-Lb(A,2)+1;
2046  imatrix re(m,n),im(m,n);
2047 
2048  blasmatmul(A,Re(B),re);
2049  blasmatmul(A,Im(B),im);
2050 
2051  SetRe(C,re);
2052  SetIm(C,im);
2053 }
2054 
2055 inline void blasmatmul(const cmatrix &A, const rmatrix &B, cimatrix &C) {
2056  int m = Ub(A,1)-Lb(A,1)+1;
2057  int n = Ub(A,2)-Lb(A,2)+1;
2058  imatrix re(m,n),im(m,n);
2059 
2060  blasmatmul(Re(A),B,re);
2061  blasmatmul(Im(A),B,im);
2062 
2063  SetRe(C,re);
2064  SetIm(C,im);
2065 }
2066 
2067 //=========================================================================
2068 
2069 /*
2070 //Computes B=A+B
2071 inline void blasadd(const rvector& A, rvector& B) {
2072  double* DA = A.to_blas_array();
2073  double* DB = B.to_blas_array();
2074 
2075  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 1 );
2076 }
2077 
2078 //Computes B=A+B
2079 inline void blasadd(const rvector& A, ivector& B) {
2080  double* DA = A.to_blas_array();
2081  double* DB = B.to_blas_array();
2082  int rnd = getround();
2083 
2084  setround(-1);
2085  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 2 );
2086  setround(1);
2087  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB+1, 2 );
2088  setround(rnd);
2089 }
2090 
2091 //Computes B=A+B
2092 inline void blasadd(const rvector& A, cvector& B) {
2093  double* DA = A.to_blas_array();
2094  double* DB = B.to_blas_array();
2095 
2096  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 2 );
2097 }
2098 
2099 //Computes B=A+B
2100 inline void blasadd(const rvector& A, civector& B) {
2101  double* DA = A.to_blas_array();
2102  double* DB = B.to_blas_array();
2103  int rnd = getround();
2104 
2105  setround(-1);
2106  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 4 );
2107  setround(1);
2108  cblas_daxpy(VecLen(A), 1.0, DA, 1, DB+1, 4 );
2109  setround(rnd);
2110 }
2111 
2112 
2113 //Computes B=A+B
2114 inline void blasadd(const cvector& A, cvector& B) {
2115  double* DA = A.to_blas_array();
2116  double* DB = B.to_blas_array();
2117 
2118 
2119  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 2 );
2120  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 2 );
2121 }
2122 
2123 //Computes B=A+B
2124 inline void blasadd(const cvector& A, civector& B) {
2125  double* DA = A.to_blas_array();
2126  double* DB = B.to_blas_array();
2127  int rnd = getround();
2128 
2129  setround(-1);
2130  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 4 );
2131  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB+2, 4 );
2132  setround(1);
2133  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 4 );
2134  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+3, 4 );
2135  setround(rnd);
2136 }
2137 
2138 //Computes B=A+B
2139 inline void blasadd(const ivector& A, ivector& B) {
2140  double* DA = A.to_blas_array();
2141  double* DB = B.to_blas_array();
2142  int rnd = getround();
2143 
2144  setround(-1);
2145 
2146  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 2 );
2147 
2148  setround(1);
2149 
2150  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 2 );
2151 
2152  setround(rnd);
2153 }
2154 
2155 //Computes B=A+B
2156 inline void blasadd(const ivector& A, civector& B) {
2157  double* DA = A.to_blas_array();
2158  double* DB = B.to_blas_array();
2159  int rnd = getround();
2160 
2161  setround(-1);
2162  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 4 );
2163  cblas_daxpy(VecLen(A), 1.0, DA, 2, DB+2, 4 );
2164 
2165  setround(1);
2166  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 4 );
2167  cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+3, 4 );
2168 
2169  setround(rnd);
2170 }
2171 
2172 //Computes B=A+B
2173 inline void blasadd(const civector& A, civector& B) {
2174  double* DA = A.to_blas_array();
2175  double* DB = B.to_blas_array();
2176  int rnd = getround();
2177 
2178  setround(-1);
2179 
2180  cblas_daxpy(VecLen(A), 1.0, DA, 4, DB, 4 );
2181  cblas_daxpy(VecLen(A), 1.0, DA+2, 4, DB+2, 4 );
2182 
2183  setround(1);
2184 
2185  cblas_daxpy(VecLen(A), 1.0, DA+1, 4, DB+1, 4 );
2186  cblas_daxpy(VecLen(A), 1.0, DA+3, 4, DB+3, 4 );
2187 
2188  setround(rnd);
2189 }
2190 
2191 //Computes C=A+B
2192 inline void blasadd(const cvector& A, const ivector& B, civector& C) {
2193  C = A;
2194  double* DB = B.to_blas_array();
2195  double* DC = C.to_blas_array();
2196  int rnd = getround();
2197 
2198  setround(-1);
2199 
2200  cblas_daxpy(VecLen(A), 1.0, DB, 4, DC, 4 );
2201 
2202  setround(1);
2203 
2204  cblas_daxpy(VecLen(A), 1.0, DB+1, 4, DC+1, 4 );
2205 
2206  setround(rnd);
2207 }
2208 
2209 //Computes B=-A+B
2210 inline void blassub(const rvector& A, rvector& B) {
2211  double* DA = A.to_blas_array();
2212  double* DB = B.to_blas_array();
2213 
2214  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 1 );
2215 }
2216 
2217 //Computes B=-A+B
2218 inline void blassub(const rvector& A, ivector& B) {
2219  double* DA = A.to_blas_array();
2220  double* DB = B.to_blas_array();
2221  int rnd = getround();
2222 
2223  setround(-1);
2224  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 2 );
2225  setround(1);
2226  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB+1, 2 );
2227  setround(rnd);
2228 }
2229 
2230 //Computes B=-A+B
2231 inline void blassub(const rvector& A, cvector& B) {
2232  double* DA = A.to_blas_array();
2233  double* DB = B.to_blas_array();
2234 
2235  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 2 );
2236 }
2237 
2238 //Computes B=-A+B
2239 inline void blassub(const rvector& A, civector& B) {
2240  double* DA = A.to_blas_array();
2241  double* DB = B.to_blas_array();
2242  int rnd = getround();
2243 
2244  setround(-1);
2245  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 4 );
2246  setround(1);
2247  cblas_daxpy(VecLen(A), -1.0, DA, 1, DB+1, 4 );
2248  setround(rnd);
2249 }
2250 
2251 
2252 //Computes B=-A+B
2253 inline void blassub(const cvector& A, cvector& B) {
2254  double* DA = A.to_blas_array();
2255  double* DB = B.to_blas_array();
2256 
2257 
2258  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB, 2 );
2259  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+1, 2 );
2260 }
2261 
2262 //Computes B=-A+B
2263 inline void blassub(const cvector& A, civector& B) {
2264  double* DA = A.to_blas_array();
2265  double* DB = B.to_blas_array();
2266  int rnd = getround();
2267 
2268  setround(-1);
2269  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB, 4 );
2270  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+2, 4 );
2271  setround(1);
2272  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+1, 4 );
2273  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+3, 4 );
2274  setround(rnd);
2275 }
2276 
2277 //Computes B=-A+B
2278 inline void blassub(const ivector& A, ivector& B) {
2279  double* DA = A.to_blas_array();
2280  double* DB = B.to_blas_array();
2281  int rnd = getround();
2282 
2283  setround(-1);
2284 
2285  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB, 2 );
2286 
2287  setround(1);
2288 
2289  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+1, 2 );
2290 
2291  setround(rnd);
2292 }
2293 
2294 //Computes B=-A+B
2295 inline void blassub(const ivector& A, civector& B) {
2296  double* DA = A.to_blas_array();
2297  double* DB = B.to_blas_array();
2298  int rnd = getround();
2299 
2300  setround(-1);
2301  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB, 4 );
2302  cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+2, 4 );
2303 
2304  setround(1);
2305  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+1, 4 );
2306  cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+3, 4 );
2307 
2308  setround(rnd);
2309 }
2310 
2311 //Computes B=-A+B
2312 inline void blassub(const civector& A, civector& B) {
2313  double* DA = A.to_blas_array();
2314  double* DB = B.to_blas_array();
2315  int rnd = getround();
2316 
2317  setround(-1);
2318 
2319  cblas_daxpy(VecLen(A), -1.0, DA+1, 4, DB, 4 );
2320  cblas_daxpy(VecLen(A), -1.0, DA+3, 4, DB+2, 4 );
2321 
2322  setround(1);
2323 
2324  cblas_daxpy(VecLen(A), -1.0, DA, 4, DB+1, 4 );
2325  cblas_daxpy(VecLen(A), -1.0, DA+2, 4, DB+3, 4 );
2326 
2327  setround(rnd);
2328 }
2329 
2330 //==============================================================
2331 
2332 //Computes B=A+B
2333 inline void blasadd(const rmatrix& A, rmatrix& B) {
2334  double* DA = A.to_blas_array();
2335  double* DB = B.to_blas_array();
2336 
2337  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 1 );
2338 }
2339 
2340 //Computes B=A+B
2341 inline void blasadd(const rmatrix& A, imatrix& B) {
2342  double* DA = A.to_blas_array();
2343  double* DB = B.to_blas_array();
2344  int rnd = getround();
2345 
2346  setround(-1);
2347  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 2 );
2348  setround(1);
2349  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB+1, 2 );
2350  setround(rnd);
2351 }
2352 
2353 //Computes B=A+B
2354 inline void blasadd(const rmatrix& A, cmatrix& B) {
2355  double* DA = A.to_blas_array();
2356  double* DB = B.to_blas_array();
2357 
2358  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 2 );
2359 }
2360 
2361 //Computes B=A+B
2362 inline void blasadd(const rmatrix& A, cimatrix& B) {
2363  double* DA = A.to_blas_array();
2364  double* DB = B.to_blas_array();
2365  int rnd = getround();
2366 
2367  setround(-1);
2368  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 4 );
2369  setround(1);
2370  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB+1, 4 );
2371  setround(rnd);
2372 }
2373 
2374 
2375 //Computes B=A+B
2376 inline void blasadd(const cmatrix& A, cmatrix& B) {
2377  double* DA = A.to_blas_array();
2378  double* DB = B.to_blas_array();
2379 
2380 
2381  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 2 );
2382  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 2 );
2383 }
2384 
2385 //Computes B=A+B
2386 inline void blasadd(const cmatrix& A, cimatrix& B) {
2387  double* DA = A.to_blas_array();
2388  double* DB = B.to_blas_array();
2389  int rnd = getround();
2390 
2391  setround(-1);
2392  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 4 );
2393  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB+2, 4 );
2394  setround(1);
2395  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 4 );
2396  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+3, 4 );
2397  setround(rnd);
2398 }
2399 
2400 //Computes B=A+B
2401 inline void blasadd(const imatrix& A, imatrix& B) {
2402  double* DA = A.to_blas_array();
2403  double* DB = B.to_blas_array();
2404  int rnd = getround();
2405 
2406  setround(-1);
2407 
2408  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 2 );
2409 
2410  setround(1);
2411 
2412  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 2 );
2413 
2414  setround(rnd);
2415 }
2416 
2417 //Computes B=A+B
2418 inline void blasadd(const imatrix& A, cimatrix& B) {
2419  double* DA = A.to_blas_array();
2420  double* DB = B.to_blas_array();
2421  int rnd = getround();
2422 
2423  setround(-1);
2424  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 4 );
2425  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB+2, 4 );
2426 
2427  setround(1);
2428  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 4 );
2429  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+3, 4 );
2430 
2431  setround(rnd);
2432 }
2433 
2434 //Computes B=A+B
2435 inline void blasadd(const cimatrix& A, cimatrix& B) {
2436  double* DA = A.to_blas_array();
2437  double* DB = B.to_blas_array();
2438  int rnd = getround();
2439 
2440  setround(-1);
2441 
2442  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 4, DB, 4 );
2443  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+2, 4, DB+2, 4 );
2444 
2445  setround(1);
2446 
2447  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 4, DB+1, 4 );
2448  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+3, 4, DB+3, 4 );
2449 
2450  setround(rnd);
2451 }
2452 
2453 //Computes C=A+B
2454 inline void blasadd(const cmatrix& A, const imatrix& B, cimatrix& C) {
2455  C = A;
2456  double* DB = B.to_blas_array();
2457  double* DC = C.to_blas_array();
2458  int rnd = getround();
2459 
2460  setround(-1);
2461 
2462  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DB, 4, DC, 4 );
2463 
2464  setround(1);
2465 
2466  cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DB+1, 4, DC+1, 4 );
2467 
2468  setround(rnd);
2469 }
2470 
2471 //Computes B=-A+B
2472 inline void blassub(const rmatrix& A, rmatrix& B) {
2473  double* DA = A.to_blas_array();
2474  double* DB = B.to_blas_array();
2475 
2476  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 1 );
2477 }
2478 
2479 //Computes B=-A+B
2480 inline void blassub(const rmatrix& A, imatrix& B) {
2481  double* DA = A.to_blas_array();
2482  double* DB = B.to_blas_array();
2483  int rnd = getround();
2484 
2485  setround(-1);
2486  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 2 );
2487  setround(1);
2488  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB+1, 2 );
2489  setround(rnd);
2490 }
2491 
2492 //Computes B=-A+B
2493 inline void blassub(const rmatrix& A, cmatrix& B) {
2494  double* DA = A.to_blas_array();
2495  double* DB = B.to_blas_array();
2496 
2497  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 2 );
2498 }
2499 
2500 //Computes B=-A+B
2501 inline void blassub(const rmatrix& A, cimatrix& B) {
2502  double* DA = A.to_blas_array();
2503  double* DB = B.to_blas_array();
2504  int rnd = getround();
2505 
2506  setround(-1);
2507  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 4 );
2508  setround(1);
2509  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB+1, 4 );
2510  setround(rnd);
2511 }
2512 
2513 
2514 //Computes B=-A+B
2515 inline void blassub(const cmatrix& A, cmatrix& B) {
2516  double* DA = A.to_blas_array();
2517  double* DB = B.to_blas_array();
2518 
2519 
2520  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB, 2 );
2521  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+1, 2 );
2522 }
2523 
2524 //Computes B=-A+B
2525 inline void blassub(const cmatrix& A, cimatrix& B) {
2526  double* DA = A.to_blas_array();
2527  double* DB = B.to_blas_array();
2528  int rnd = getround();
2529 
2530  setround(-1);
2531  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB, 4 );
2532  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+2, 4 );
2533  setround(1);
2534  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+1, 4 );
2535  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+3, 4 );
2536  setround(rnd);
2537 }
2538 
2539 //Computes B=-A+B
2540 inline void blassub(const imatrix& A, imatrix& B) {
2541  double* DA = A.to_blas_array();
2542  double* DB = B.to_blas_array();
2543  int rnd = getround();
2544 
2545  setround(-1);
2546 
2547  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB, 2 );
2548 
2549  setround(1);
2550 
2551  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+1, 2 );
2552 
2553  setround(rnd);
2554 }
2555 
2556 //Computes B=-A+B
2557 inline void blassub(const imatrix& A, cimatrix& B) {
2558  double* DA = A.to_blas_array();
2559  double* DB = B.to_blas_array();
2560  int rnd = getround();
2561 
2562  setround(-1);
2563  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB, 4 );
2564  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+2, 4 );
2565 
2566  setround(1);
2567  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+1, 4 );
2568  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+3, 4 );
2569 
2570  setround(rnd);
2571 }
2572 
2573 //Computes B=-A+B
2574 inline void blassub(const cimatrix& A, cimatrix& B) {
2575  double* DA = A.to_blas_array();
2576  double* DB = B.to_blas_array();
2577  int rnd = getround();
2578 
2579  setround(-1);
2580 
2581  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 4, DB, 4 );
2582  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+3, 4, DB+2, 4 );
2583 
2584  setround(1);
2585 
2586  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 4, DB+1, 4 );
2587  cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+2, 4, DB+3, 4 );
2588 
2589  setround(rnd);
2590 } */
2591 
2592 #endif
2593 #endif
2594 
2595 }
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
int VecLen(const scimatrix_subv &S)
Returns the length of the subvector.
Definition: scimatrix.hpp:9966
void Resize(cimatrix &A)
Resizes the matrix.
Definition: cimatrix.inl:1211
int Ub(const cimatrix &rm, const int &i)
Returns the upper bound index.
Definition: cimatrix.inl:1163
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737