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
34namespace 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
94inline 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
104inline 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)
123inline 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
163inline 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
208inline 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
222inline 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
236inline 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
246inline 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
256inline 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
267inline 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
300inline 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!
381inline 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!
434inline 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
484inline 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
514inline 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
544inline 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
580inline 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
621inline 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
661inline 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
714inline 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
761inline 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
809inline 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
864inline 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
925inline 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
946inline 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
973inline 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
1010inline 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
1048inline 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
1086inline 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
1100inline 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
1114inline 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
1128inline 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
1142inline 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
1149inline 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
1159inline 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
1170inline 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
1188inline 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
1202inline 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
1216inline 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
1230inline 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
1244inline 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
1267inline 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
1293inline 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
1318inline 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
1338inline 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
1378inline 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
1494inline 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
1602inline 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
1714inline 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
1731inline 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
1748inline 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
1770inline 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
1784inline 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
1799inline 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
1814inline 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
1829inline 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
1844inline 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
1859inline 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
1875inline 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
1932inline 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
1987inline 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
2043inline 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
2055inline 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
2071inline 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
2079inline 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
2092inline 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
2100inline 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
2114inline 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
2124inline 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
2139inline 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
2156inline 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
2173inline 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
2192inline 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
2210inline 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
2218inline 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
2231inline 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
2239inline 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
2253inline 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
2263inline 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
2278inline 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
2295inline 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
2312inline 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
2333inline 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
2341inline 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
2354inline 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
2362inline 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
2376inline 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
2386inline 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
2401inline 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
2418inline 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
2435inline 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
2454inline 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
2472inline 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
2480inline 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
2493inline 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
2501inline 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
2515inline 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
2525inline 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
2540inline 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
2557inline 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
2574inline 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}
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.
int Ub(const cimatrix &rm, const int &i) noexcept
Returns the upper bound index.
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
void Resize(cimatrix &A) noexcept
Resizes the matrix.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.