32#include "cxsc_blas.hpp"
36#ifndef _CXSC_BLAS_SETROUND
37#define _CXSC_BLAS_SETROUND
45 static inline void setround(
int r) {
48 fesetround(FE_DOWNWARD);
51 fesetround(FE_TONEAREST);
54 fesetround(FE_UPWARD);
57 fesetround(FE_TOWARDZERO);
60 fesetround(FE_TONEAREST);
71 static inline int getround() {
72 switch(fegetround()) {
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) {
98 double* xd = (
double*) x.to_blas_array();
99 double* yd = (
double*) y.to_blas_array();
101 res = cblas_ddot(n, xd, inc, yd, inc);
104inline void blasdot(
const rvector& x,
const rvector& y, interval& res) {
107 int rnd = getround();
109 double* xd = (
double*) x.to_blas_array();
110 double* yd = (
double*) y.to_blas_array();
113 SetInf(res, cblas_ddot(n, xd, inc, yd, inc));
115 SetSup(res, cblas_ddot(n, xd, inc, yd, inc));
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) {
126 const int lb1 =
Lb(x);
127 const int lb2 =
Lb(y);
128 int rnd = getround();
130 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
140 bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
142 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
146 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
148 SetRe(res, interval(res_inf,res_sup));
150 bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
152 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
156 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
158 SetIm(res, interval(res_inf,res_sup));
163inline void blasdot(
const ivector& x,
const cvector& y, cinterval& res) {
166 const int lb1 =
Lb(x);
167 const int lb2 =
Lb(y);
168 int rnd = getround();
170 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
180 bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
182 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
186 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
188 SetRe(res, interval(res_inf,res_sup));
190 bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
192 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
196 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
198 SetIm(res, interval(res_inf,res_sup));
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) {
211 double res_re, res_im;
213 double* xd = (
double*) x.to_blas_array();
214 double* yd = (
double*) y.to_blas_array();
216 res_re = cblas_ddot(n, xd, inc, yd, inc2);
217 res_im = cblas_ddot(n, xd, inc, yd+1, inc2);
219 res = complex(res_re,res_im);
222inline void blasdot(
const cvector& x,
const rvector& y, complex& res) {
225 double res_re, res_im;
227 double* xd = (
double*) x.to_blas_array();
228 double* yd = (
double*) y.to_blas_array();
230 res_re = cblas_ddot(n, xd, inc2, yd, inc);
231 res_im = cblas_ddot(n, xd+1, inc2, yd, inc);
233 res = complex(res_re,res_im);
236inline void blasdot(
const cvector& x,
const cvector& y, complex& res) {
240 double* xd = (
double*) x.to_blas_array();
241 double* yd = (
double*) y.to_blas_array();
243 cblas_zdotu_sub(n, xd, inc, yd, inc, (
void*)&res);
246inline void blasdot(
const rvector& x,
const cvector& y, cinterval& res) {
256inline void blasdot(
const cvector& x,
const rvector& y, cinterval& res) {
267inline void blasdot(
const cvector& x,
const cvector& y, cinterval& res) {
268 int rnd = getround();
271 double res_re, res_im;
273 double* xd = (
double*) x.to_blas_array();
274 double* yd = (
double*) y.to_blas_array();
277 res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc);
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));
285 res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc);
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));
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) {
302 for(
int i=0 ; i<n ; i++) {
304 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) {
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]);
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]);
323 }
else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) {
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) {
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]);
336 x_inf[i+1] = -Sup(x[i+lb1]);
337 y_inf[i+1] = Inf(y[i+lb2]);
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]);
344 x_sup[i+1] = Sup(x[i+lb1]);
345 y_sup[i+1] = Sup(y[i+lb2]);
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]);
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]);
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]);
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) {
383 for(
int i=0 ; i<n ; i++) {
385 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 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];
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];
399 }
else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 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];
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];
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];
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];
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++) {
437 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 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];
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];
451 }
else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 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];
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];
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];
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];
484inline void blasdot(
const ivector& x,
const ivector& y, interval& res) {
487 const int lb1 =
Lb(x);
488 const int lb2 =
Lb(y);
489 int rnd = getround();
491 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
501 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
503 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
507 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
511 res = interval(res_inf,res_sup);
514inline void blasdot(
const ivector& x,
const rvector& y, interval& res) {
517 const int lb1 =
Lb(x);
518 const int lb2 =
Lb(y);
519 int rnd = getround();
521 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
531 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
533 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
537 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
541 res = interval(res_inf,res_sup);
544inline void blasdot(
const rvector& x,
const ivector& y, interval& res) {
547 const int lb1 =
Lb(x);
548 const int lb2 =
Lb(y);
549 int rnd = getround();
551 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
561 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
563 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
567 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
571 res = interval(res_inf,res_sup);
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) {
583 const int lb1 =
Lb(x);
584 const int lb2 =
Lb(y);
585 int rnd = getround();
587 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
597 bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
599 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
603 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
605 SetRe(res, interval(res_inf,res_sup));
607 bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
609 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
613 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
615 SetIm(res, interval(res_inf,res_sup));
621inline void blasdot(
const civector& x,
const rvector& y, cinterval& res) {
624 const int lb1 =
Lb(x);
625 const int lb2 =
Lb(y);
626 int rnd = getround();
628 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
638 bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
640 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
644 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
646 SetRe(res, interval(res_inf,res_sup));
648 bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2);
650 res_sup = cblas_ddot(n, dxs, inc, dys, inc);
654 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc);
656 SetIm(res, interval(res_inf,res_sup));
661inline void blasdot(
const civector& x,
const civector& y, cinterval& res) {
664 const int lb1 =
Lb(x);
665 const int lb2 =
Lb(y);
666 int rnd = getround();
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;
673 double* dxmid = xmid.to_blas_array();
674 double* dymid = ymid.to_blas_array();
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]);
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]));
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]);
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]);
696 cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_sup);
698 res_sup += complex(crad_re,crad_im);
702 cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_inf);
704 res_inf -= complex(crad_re,crad_im);
708 UncheckedSetInf(res,res_inf);
709 UncheckedSetSup(res,res_sup);
714inline void blasdot(
const civector& x,
const cvector& y, cinterval& res) {
717 const int lb1 =
Lb(x);
718 int rnd = getround();
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;
725 double* dxmid = xmid.to_blas_array();
726 double* dy = y.to_blas_array();
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]);
734 tmp1 =
abs(Re(y[i+1]));
735 tmp2 =
abs(Im(y[i+1]));
737 crad_re += Re(xrad[i+1]) * tmp1;
738 crad_re += Im(xrad[i+1]) * tmp2;
740 crad_im += Re(xrad[i+1]) * tmp2;
741 crad_im += Im(xrad[i+1]) * tmp1;
744 cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_sup);
746 res_sup += complex(crad_re,crad_im);
750 cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_inf);
752 res_inf -= complex(crad_re,crad_im);
756 UncheckedSetInf(res,res_inf);
757 UncheckedSetSup(res,res_sup);
761inline void blasdot(
const cvector& x,
const civector& y, cinterval& res) {
764 const int lb1 =
Lb(x);
765 int rnd = getround();
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;
772 double* dx = x.to_blas_array();
773 double* dymid = ymid.to_blas_array();
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]);
781 tmp1 =
abs(Re(x[i+1]));
782 tmp2 =
abs(Im(x[i+1]));
784 crad_re += Re(yrad[i+1]) * tmp1;
785 crad_re += Im(yrad[i+1]) * tmp2;
787 crad_im += Re(yrad[i+1]) * tmp2;
788 crad_im += Im(yrad[i+1]) * tmp1;
791 cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_sup);
793 res_sup += complex(crad_re,crad_im);
797 cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_inf);
799 res_inf -= complex(crad_re,crad_im);
803 UncheckedSetInf(res,res_inf);
804 UncheckedSetSup(res,res_sup);
809inline void blasdot(
const civector& x,
const ivector& y, cinterval& res) {
811 const int inc=1, inc2=2;
812 const int lb1 =
Lb(x);
813 const int lb2 =
Lb(y);
814 int rnd = getround();
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;
822 double* dxmid = xmid.to_blas_array();
823 double* dymid = ymid.to_blas_array();
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]);
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]));
838 crad_re += Re(xrad[i+1]) * tmp1 + tmp3 * yrad[i+1];
839 crad_re += Im(xrad[i+1]) * tmp2 + tmp4 * yrad[i+1];
841 crad_im += Re(xrad[i+1]) * tmp2 + tmp3 * yrad[i+1];
842 crad_im += Im(xrad[i+1]) * tmp1 + tmp4 * yrad[i+1];
845 SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc));
846 SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc));
848 res_sup += complex(crad_re,crad_im);
852 SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc));
853 SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc));
855 res_inf -= complex(crad_re,crad_im);
859 UncheckedSetInf(res,res_inf);
860 UncheckedSetSup(res,res_sup);
864inline void blasdot(
const ivector& x,
const civector& y, cinterval& res) {
866 const int inc=1, inc2=2;
867 const int lb1 =
Lb(x);
868 const int lb2 =
Lb(y);
869 int rnd = getround();
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;
877 double* dxmid = xmid.to_blas_array();
878 double* dymid = ymid.to_blas_array();
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]);
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]);
893 crad_re += xrad[i+1] * tmp1 + tmp3 * Re(yrad[i+1]);
894 crad_re += xrad[i+1] * tmp2 + tmp4 * Im(yrad[i+1]);
896 crad_im += xrad[i+1] * tmp2 + tmp3 * Im(yrad[i+1]);
897 crad_im += xrad[i+1] * tmp1 + tmp4 * Re(yrad[i+1]);
900 SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2));
901 SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2));
903 res_sup += complex(crad_re,crad_im);
907 SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2));
908 SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2));
910 res_inf -= complex(crad_re,crad_im);
914 UncheckedSetInf(res,res_inf);
915 UncheckedSetSup(res,res_sup);
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) {
927 double* DA = (
double*) A.to_blas_array();
928 double* dx = (
double*) x.to_blas_array();
929 double* dr = (
double*) r.to_blas_array();
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;
936 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
937 alpha, DA, n, dx, inc, beta, dr, inc);
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();
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;
955 double* dr =
new double[n];
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]);
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]);
973inline void blasmvmul(
const rmatrix& A,
const ivector& x, ivector& r) {
975 const int m =
Ub(A,1)-
Lb(A,1)+1;
977 const int lb1 =
Lb(A,2);
978 const int lb2 =
Lb(x);
979 int rnd = getround();
981 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
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);
996 res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
997 UncheckedSetInf(r[i], -res_inf);
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) {
1012 const int m =
Ub(A,1)-
Lb(A,1)+1;
1014 const int lb1 =
Lb(A,2);
1015 const int lb2 =
Lb(x);
1016 int rnd = getround();
1018 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
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);
1033 res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
1034 UncheckedSetInf(r[i], -res_inf);
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) {
1050 const int m =
Ub(A,1)-
Lb(A,1)+1;
1052 const int lb1 =
Lb(A,2);
1053 const int lb2 =
Lb(x);
1054 int rnd = getround();
1056 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n);
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;
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);
1071 res_inf = cblas_ddot(n, dxi, inc, dyi, inc);
1072 UncheckedSetInf(r[i], -res_inf);
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;
1089 for(
int i=0 ; i<m ; i++)
1090 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1103 for(
int i=0 ; i<m ; i++)
1104 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1117 for(
int i=0 ; i<m ; i++)
1118 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1131 for(
int i=0 ; i<m ; i++)
1132 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1145 for(
int i=0 ; i<m ; i++)
1146 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
1149inline void blasmvmul(
const cmatrix& A,
const cvector& x, civector& r) {
1151 blasmvmul(Re(A),Re(x),tmp1);
1152 blasmvmul(Im(A),Im(x),tmp2);
1154 blasmvmul(Re(A),Im(x),tmp1);
1155 blasmvmul(Im(A),Re(x),tmp2);
1159inline void blasmvmul(
const cmatrix& A,
const rvector& x, civector& r) {
1161 ivector re(n),im(n);
1163 blasmvmul(Re(A),x,re);
1164 blasmvmul(Im(A),x,im);
1170inline void blasmvmul(
const rmatrix& A,
const cvector& x, civector& r) {
1172 ivector re(n),im(n);
1174 blasmvmul(A,Re(x),re);
1175 blasmvmul(A,Im(x),im);
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;
1191 for(
int i=0 ; i<m ; i++)
1192 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1205 for(
int i=0 ; i<m ; i++)
1206 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1219 for(
int i=0 ; i<m ; i++)
1220 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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;
1233 for(
int i=0 ; i<m ; i++)
1234 blasdot(A[i+
Lb(A,1)], x, r[i+
Lb(r)]);
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) {
1246 double* DA = (
double*) A.to_blas_array();
1247 double* dx = (
double*) x.to_blas_array();
1248 double* dr = (
double*) r.to_blas_array();
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);
1256 cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n,
1257 &alpha, DA, n, dx, inc, &beta, dr, inc);
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) {
1269 double* DA = (
double*) A.to_blas_array();
1270 double* dx = (
double*) x.to_blas_array();
1271 double* dr = (
double*) r.to_blas_array();
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;
1279 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
1280 alpha, DA, m, dx, inc, beta, dr, inc);
1282 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n,
1283 alpha, DA, n, dx+1, inc, beta, dr+1, inc);
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) {
1296 double* DA = (
double*) A.to_blas_array();
1297 double* dx = (
double*) tmp.to_blas_array();
1298 double* dr = (
double*) r.to_blas_array();
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);
1306 cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n,
1307 &alpha, DA, n, dx, inc, &beta, dr, inc);
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) {
1320 double* DA = (
double*) A.to_blas_array();
1321 double* DB = (
double*) B.to_blas_array();
1322 double* DC = (
double*) C.to_blas_array();
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;
1329 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1330 n, alpha, DA, n, DB, o, beta, DC, o);
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();
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;
1351 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1352 n, alpha, DA, n, DB, o, beta, DC, o);
1354 for(
int i=0 ; i<m ; i++) {
1355 for(
int j=0 ; j<o ; j++) {
1357 UncheckedSetSup(C[i+
Lb(C,1)][j+
Lb(C,2)], DC[ind]);
1362 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1363 n, alpha, DA, n, DB, o, beta, DC, o);
1365 for(
int i=0 ; i<m ; i++) {
1366 for(
int j=0 ; j<o ; j++) {
1368 UncheckedSetInf(C[i+
Lb(C,1)][j+
Lb(C,2)], DC[ind]);
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;
1383 int lbC_1 =
Lb(C,1);
int lbC_2 =
Lb(C,2);
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();
1396 const double alpha = 1.0, beta = 0.0;
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]);
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]);
1424 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1425 n, alpha, Amid, n, Bmid, o, beta, C1, o);
1429 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1430 n, alpha, Amid, n, Bmid, o, beta, C2, o);
1435 double* Cmid =
new double[m*o];
1436 double* Crad =
new double[m*o];
1438 for(
int i=0 ; i<m ; i++) {
1439 for(
int j=0 ; j<o ; j++) {
1441 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1442 Crad[ind] = Cmid[ind] - C1[ind];
1446 for(
int i=0 ; i<n ; i++) {
1447 for(
int j=0 ; j<o ; j++) {
1449 Babs[ind] += Brad[ind];
1453 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1454 n, alpha, Arad, n, Babs, o, beta, C1, o);
1456 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1457 n, alpha, Aabs, n, Brad, o, beta, C2, o);
1464 Resize(C, lbC_1, lbC_1+m-1, lbC_2, lbC_2+o-1);
1467 for(
int i=0 ; i<m ; i++) {
1468 for(
int j=0 ; j<o ; j++) {
1470 Crad[ind] += C1[ind] + C2[ind];
1471 UncheckedSetSup(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] + Crad[ind]);
1477 for(
int i=0 ; i<m ; i++) {
1478 for(
int j=0 ; j<o ; j++) {
1480 UncheckedSetInf(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] - Crad[ind]);
1494inline void blasmatmul(
const rmatrix &A,
const imatrix &B, imatrix &C) {
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;
1500 int lb1_C =
Lb(C,1);
1501 int lb2_C =
Lb(C,2);
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();
1511 const double alpha = 1.0;
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]));
1529 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1530 n, alpha, DA, n, Bmid, o, beta, C1, o);
1534 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1535 n, alpha, DA, n, Bmid, o, beta, C2, o);
1539 double* Cmid =
new double[m*o];
1540 for(
int i=0 ; i<m ; i++) {
1541 for(
int j=0 ; j<o ; j++) {
1543 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1549 double* Crad =
new double[m*o];
1550 for(
int i=0 ; i<m ; i++) {
1551 for(
int j=0 ; j<o ; j++) {
1553 Crad[ind] = Cmid[ind] - C1[ind];
1559 double* Aabs =
new double[m*n];
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]);
1570 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1571 n, alpha, Aabs, n, Brad, o, beta, Crad, o);
1576 Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1579 for(
int i=0 ; i<m ; i++) {
1580 for(
int j=0 ; j<o ; j++) {
1582 UncheckedSetSup(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] + Crad[ind]);
1588 for(
int i=0 ; i<m ; i++) {
1589 for(
int j=0 ; j<o ; j++) {
1591 UncheckedSetInf(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] - Crad[ind]);
1602inline void blasmatmul(
const imatrix &A,
const rmatrix &B, imatrix &C) {
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;
1608 int lb1_C =
Lb(C,1);
1609 int lb2_C =
Lb(C,2);
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();
1619 const double alpha = 1.0;
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]));
1637 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1638 n, alpha, Amid, n, DB, o, beta, C1, o);
1642 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1643 n, alpha, Amid, n, DB, o, beta, C2, o);
1647 double* Cmid =
new double[m*o];
1648 for(
int i=0 ; i<m ; i++) {
1649 for(
int j=0 ; j<o ; j++) {
1651 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]);
1657 double* Crad =
new double[m*o];
1658 for(
int i=0 ; i<m ; i++) {
1659 for(
int j=0 ; j<o ; j++) {
1661 Crad[ind] = Cmid[ind] - C1[ind];
1667 double* Babs =
new double[n*o];
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]);
1678 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1679 n, alpha, Arad, n, Babs, o, beta, Crad, o);
1684 Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1);
1687 for(
int i=0 ; i<m ; i++) {
1688 for(
int j=0 ; j<o ; j++) {
1690 UncheckedSetSup(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] + Crad[ind]);
1696 for(
int i=0 ; i<m ; i++) {
1697 for(
int j=0 ; j<o ; j++) {
1699 UncheckedSetInf(C[i+
Lb(C,1)][j+
Lb(C,2)], Cmid[ind] - Crad[ind]);
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) {
1716 double* DA = (
double*) A.to_blas_array();
1717 double* DB = (
double*) B.to_blas_array();
1718 double* DC = (
double*) C.to_blas_array();
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);
1726 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1727 n, &alpha, DA, n, DB, o, &beta, DC, o);
1731inline void blasmatmul(
const cmatrix &A,
const rmatrix &B, cmatrix &C) {
1734 double* DA = (
double*) A.to_blas_array();
1735 double* DB = (
double*) tmp.to_blas_array();
1736 double* DC = (
double*) C.to_blas_array();
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);
1744 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1745 n, &alpha, DA, n, DB, o, &beta, DC, o);
1748inline void blasmatmul(
const rmatrix &A,
const cmatrix &B, cmatrix &C) {
1751 double* DA = (
double*) tmp.to_blas_array();
1752 double* DB = (
double*) B.to_blas_array();
1753 double* DC = (
double*) C.to_blas_array();
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);
1761 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o,
1762 n, &alpha, DA, n, DB, o, &beta, DC, o);
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);
1775 blasmatmul(Re(A),Re(B),tmp1);
1776 blasmatmul(Im(A),Im(B),tmp2);
1778 blasmatmul(Re(A),Im(B),tmp1);
1779 blasmatmul(Im(A),Re(B),tmp2);
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;
1790 blasmatmul(Re(A),B,T);
1794 blasmatmul(Im(A),B,T);
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;
1805 blasmatmul(A,Re(B),T);
1809 blasmatmul(A,Im(B),T);
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;
1820 blasmatmul(A,Re(B),T);
1824 blasmatmul(A,Im(B),T);
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;
1835 blasmatmul(Re(A),B,T);
1839 blasmatmul(Im(A),B,T);
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;
1850 blasmatmul(A,Re(B),T);
1854 blasmatmul(A,Im(B),T);
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;
1865 blasmatmul(Re(A),B,T);
1869 blasmatmul(Im(A),B,T);
1875inline void blasmatmul(
const cimatrix &A,
const cmatrix &B, cimatrix &C) {
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();
1883 blasmatmul(Re(A),Re(B),T);
1887 blasmatmul(-Im(A),Im(B),T);
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]);
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]);
1907 blasmatmul(Re(A),Im(B),T);
1911 blasmatmul(Im(A),Re(B),T);
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]);
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]);
1932inline void blasmatmul(
const cmatrix &A,
const cimatrix &B, cimatrix &C) {
1933 int rnd = getround();
1935 imatrix T(
Lb(C,1),
Ub(C,1),
Lb(C,2),
Ub(C,2));
1937 blasmatmul(Re(A),Re(B),T);
1941 blasmatmul(-Im(A),Im(B),T);
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]);
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]);
1961 blasmatmul(Re(A),Im(B),T);
1965 blasmatmul(Im(A),Re(B),T);
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]);
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]);
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();
1994 blasmatmul(Re(A),Re(B),T);
1998 blasmatmul(-Im(A),Im(B),T);
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]);
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]);
2018 blasmatmul(Re(A),Im(B),T);
2022 blasmatmul(Im(A),Re(B),T);
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]);
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]);
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);
2048 blasmatmul(A,Re(B),re);
2049 blasmatmul(A,Im(B),im);
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);
2060 blasmatmul(Re(A),B,re);
2061 blasmatmul(Im(A),B,im);
The namespace cxsc, providing all functionality of the class library C-XSC.
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.
void Resize(cimatrix &A) noexcept
Resizes the matrix.
int Lb(const cimatrix &rm, const int &i) noexcept
Returns the lower bound index.