C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: cxsc_blas.inl,v 1.20 2014/01/30 17:23:44 cxsc Exp $ */ 00025 00026 /* 00027 ** FastPLSS: A library of (parallel) verified linear (interval) system 00028 ** solvers using C-XSC (V 0.2) 00029 */ 00030 00031 #include <fenv.h> 00032 #include "cxsc_blas.hpp" 00033 00034 namespace cxsc { 00035 00036 #ifndef _CXSC_BLAS_SETROUND 00037 #define _CXSC_BLAS_SETROUND 00038 /* 00039 Sets the rounding mode according to the parameter r: 00040 r=-1: Round downwards 00041 r=0 : Round to nearest 00042 r=1 : Round upwards 00043 r=2 : Round toward zero 00044 */ 00045 static inline void setround(int r) { 00046 switch(r) { 00047 case -1: 00048 fesetround(FE_DOWNWARD); 00049 break; 00050 case 0: 00051 fesetround(FE_TONEAREST); 00052 break; 00053 case 1: 00054 fesetround(FE_UPWARD); 00055 break; 00056 case 2: 00057 fesetround(FE_TOWARDZERO); 00058 break; 00059 default: 00060 fesetround(FE_TONEAREST); 00061 } 00062 } 00063 00064 /* 00065 Sets the rounding mode according to the parameter r: 00066 r=-1: Round downwards 00067 r=0 : Round to nearest 00068 r=1 : Round upwards 00069 r=2 : Round toward zero 00070 */ 00071 static inline int getround() { 00072 switch(fegetround()) { 00073 case FE_DOWNWARD: 00074 return -1; 00075 break; 00076 case FE_TONEAREST: 00077 return 0; 00078 break; 00079 case FE_UPWARD: 00080 return 1; 00081 break; 00082 case FE_TOWARDZERO: 00083 return 2; 00084 break; 00085 default: 00086 return 0; 00087 } 00088 } 00089 #endif 00090 00091 #ifdef _CXSC_BLAS_RVECTOR 00092 #ifndef _CXSC_BLAS_RVECTOR_INC 00093 #define _CXSC_BLAS_RVECTOR_INC 00094 inline void blasdot(const rvector& x, const rvector& y, real& res) { 00095 int n = VecLen(x); 00096 int inc=1; 00097 00098 double* xd = (double*) x.to_blas_array(); 00099 double* yd = (double*) y.to_blas_array(); 00100 00101 res = cblas_ddot(n, xd, inc, yd, inc); 00102 } 00103 00104 inline void blasdot(const rvector& x, const rvector& y, interval& res) { 00105 int n = VecLen(x); 00106 int inc=1; 00107 int rnd = getround(); 00108 00109 double* xd = (double*) x.to_blas_array(); 00110 double* yd = (double*) y.to_blas_array(); 00111 00112 setround(-1); 00113 SetInf(res, cblas_ddot(n, xd, inc, yd, inc)); 00114 setround(1); 00115 SetSup(res, cblas_ddot(n, xd, inc, yd, inc)); 00116 setround(rnd); 00117 } 00118 #endif 00119 #endif 00120 00121 #if defined(_CXSC_BLAS_CVECTOR) && defined(_CXSC_BLAS_IVECTOR) 00122 #if !defined(_CXSC_BLAS_CVECTOR_INC) || !defined(_CXSC_BLAS_IVECTOR_INC) 00123 inline void blasdot(const cvector& x, const ivector& y, cinterval& res) { 00124 const int n = VecLen(x); 00125 const int inc=1; 00126 const int lb1 = Lb(x); 00127 const int lb2 = Lb(y); 00128 int rnd = getround(); 00129 00130 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00131 00132 double* dxi = x_inf.to_blas_array(); 00133 double* dxs = x_sup.to_blas_array(); 00134 double* dyi = y_inf.to_blas_array(); 00135 double* dys = y_sup.to_blas_array(); 00136 double res_inf,res_sup; 00137 00138 setround(1); 00139 00140 bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00141 00142 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00143 00144 //setround(-1); 00145 00146 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00147 00148 SetRe(res, interval(res_inf,res_sup)); 00149 00150 bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00151 00152 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00153 00154 //setround(-1); 00155 00156 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00157 00158 SetIm(res, interval(res_inf,res_sup)); 00159 00160 setround(rnd); 00161 } 00162 00163 inline void blasdot(const ivector& x, const cvector& y, cinterval& res) { 00164 const int n = VecLen(x); 00165 const int inc=1; 00166 const int lb1 = Lb(x); 00167 const int lb2 = Lb(y); 00168 int rnd = getround(); 00169 00170 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00171 00172 double* dxi = x_inf.to_blas_array(); 00173 double* dxs = x_sup.to_blas_array(); 00174 double* dyi = y_inf.to_blas_array(); 00175 double* dys = y_sup.to_blas_array(); 00176 double res_inf,res_sup; 00177 00178 setround(1); 00179 00180 bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00181 00182 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00183 00184 //setround(-1); 00185 00186 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00187 00188 SetRe(res, interval(res_inf,res_sup)); 00189 00190 bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00191 00192 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00193 00194 //setround(-1); 00195 00196 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00197 00198 SetIm(res, interval(res_inf,res_sup)); 00199 00200 setround(rnd); 00201 } 00202 #endif 00203 #endif 00204 00205 #ifdef _CXSC_BLAS_CVECTOR 00206 #ifndef _CXSC_BLAS_CVECTOR_INC 00207 #define _CXSC_BLAS_CVECTOR_INC 00208 inline void blasdot(const rvector& x, const cvector& y, complex& res) { 00209 int n = VecLen(x); 00210 int inc=1, inc2=2; 00211 double res_re, res_im; 00212 00213 double* xd = (double*) x.to_blas_array(); 00214 double* yd = (double*) y.to_blas_array(); 00215 00216 res_re = cblas_ddot(n, xd, inc, yd, inc2); 00217 res_im = cblas_ddot(n, xd, inc, yd+1, inc2); 00218 00219 res = complex(res_re,res_im); 00220 } 00221 00222 inline void blasdot(const cvector& x, const rvector& y, complex& res) { 00223 int n = VecLen(x); 00224 int inc=1, inc2=2; 00225 double res_re, res_im; 00226 00227 double* xd = (double*) x.to_blas_array(); 00228 double* yd = (double*) y.to_blas_array(); 00229 00230 res_re = cblas_ddot(n, xd, inc2, yd, inc); 00231 res_im = cblas_ddot(n, xd+1, inc2, yd, inc); 00232 00233 res = complex(res_re,res_im); 00234 } 00235 00236 inline void blasdot(const cvector& x, const cvector& y, complex& res) { 00237 int n = VecLen(x); 00238 int inc=1; 00239 00240 double* xd = (double*) x.to_blas_array(); 00241 double* yd = (double*) y.to_blas_array(); 00242 00243 cblas_zdotu_sub(n, xd, inc, yd, inc, (void*)&res); 00244 } 00245 00246 inline void blasdot(const rvector& x, const cvector& y, cinterval& res) { 00247 interval re,im; 00248 00249 blasdot(x,Re(y),re); 00250 blasdot(x,Im(y),im); 00251 00252 SetRe(res,re); 00253 SetIm(res,im); 00254 } 00255 00256 inline void blasdot(const cvector& x, const rvector& y, cinterval& res) { 00257 interval re,im; 00258 00259 blasdot(Re(x),y,re); 00260 blasdot(Im(x),y,im); 00261 00262 SetRe(res,re); 00263 SetIm(res,im); 00264 } 00265 00266 00267 inline void blasdot(const cvector& x, const cvector& y, cinterval& res) { 00268 int rnd = getround(); 00269 int n = VecLen(x); 00270 int inc=2; 00271 double res_re, res_im; 00272 00273 double* xd = (double*) x.to_blas_array(); 00274 double* yd = (double*) y.to_blas_array(); 00275 00276 setround(1); 00277 res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc); 00278 00279 setround(-1); 00280 res_re += cblas_ddot(n, xd, inc, yd, inc); 00281 res_im = cblas_ddot(n, xd, inc, yd+1, inc); 00282 res_im += cblas_ddot(n, xd+1, inc, yd, inc); 00283 UncheckedSetInf(res,complex(res_re,res_im)); 00284 00285 res_re = - cblas_ddot(n, xd+1, inc, yd+1, inc); 00286 setround(1); 00287 res_re += cblas_ddot(n, xd, inc, yd, inc); 00288 res_im = cblas_ddot(n, xd, inc, yd+1, inc); 00289 res_im += cblas_ddot(n, xd+1, inc, yd, inc); 00290 UncheckedSetSup(res,complex(res_re,res_im)); 00291 00292 setround(rnd); 00293 } 00294 #endif 00295 #endif 00296 00297 #ifdef _CXSC_BLAS_IVECTOR 00298 #ifndef _CXSC_BLAS_IVECTOR_INC 00299 #define _CXSC_BLAS_IVECTOR_INC 00300 inline void bsort(const ivector &x, const ivector &y, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb1, int lb2) { 00301 00302 for(int i=0 ; i<n ; i++) { 00303 00304 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) { 00305 00306 if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) { 00307 x_inf[i+1] = -Inf(x[i+lb1]); 00308 y_inf[i+1] = Inf(y[i+lb2]); 00309 x_sup[i+1] = Sup(x[i+lb1]); 00310 y_sup[i+1] = Sup(y[i+lb2]); 00311 } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) { 00312 x_inf[i+1] = -Sup(x[i+lb1]); 00313 y_inf[i+1] = Inf(y[i+lb2]); 00314 x_sup[i+1] = Sup(x[i+lb1]); 00315 y_sup[i+1] = Sup(y[i+lb2]); 00316 } else { 00317 x_inf[i+1] = -Sup(x[i+lb1]); 00318 y_inf[i+1] = Inf(y[i+lb2]); 00319 x_sup[i+1] = Inf(x[i+lb1]); 00320 y_sup[i+1] = Sup(y[i+lb2]); 00321 } 00322 00323 } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) { 00324 00325 if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) { 00326 x_inf[i+1] = -Inf(x[i+lb1]); 00327 y_inf[i+1] = Sup(y[i+lb2]); 00328 x_sup[i+1] = Sup(x[i+lb1]); 00329 y_sup[i+1] = Sup(y[i+lb2]); 00330 } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) { 00331 00332 if((-Inf(x[i+lb1]))*Sup(y[i+lb2]) >= Sup(x[i+lb1])*(-Inf(y[i+lb2]))) { 00333 x_inf[i+1] = -Inf(x[i+lb1]); 00334 y_inf[i+1] = Sup(y[i+lb2]); 00335 } else { 00336 x_inf[i+1] = -Sup(x[i+lb1]); 00337 y_inf[i+1] = Inf(y[i+lb2]); 00338 } 00339 00340 if(Inf(x[i+lb1])*Inf(y[i+lb2]) >= Sup(x[i+lb1])*Sup(y[i+lb2])) { 00341 x_sup[i+1] = Inf(x[i+lb1]); 00342 y_sup[i+1] = Inf(y[i+lb2]); 00343 } else { 00344 x_sup[i+1] = Sup(x[i+lb1]); 00345 y_sup[i+1] = Sup(y[i+lb2]); 00346 } 00347 00348 } else { 00349 x_inf[i+1] = Sup(x[i+lb1]); 00350 y_inf[i+1] = -Inf(y[i+lb2]); 00351 x_sup[i+1] = Inf(x[i+lb1]); 00352 y_sup[i+1] = Inf(y[i+lb2]); 00353 } 00354 00355 } else { 00356 00357 if(Inf(y[i+lb2]) >= 0 && Sup(y[i+lb2]) >= 0) { 00358 x_inf[i+1] = -Inf(x[i+lb1]); 00359 y_inf[i+1] = Sup(y[i+lb2]); 00360 x_sup[i+1] = Sup(x[i+lb1]); 00361 y_sup[i+1] = Inf(y[i+lb2]); 00362 } else if(Inf(y[i+lb2]) < 0 && Sup(y[i+lb2]) >= 0) { 00363 x_inf[i+1] = -Inf(x[i+lb1]); 00364 y_inf[i+1] = Sup(y[i+lb2]); 00365 x_sup[i+1] = Inf(x[i+lb1]); 00366 y_sup[i+1] = Inf(y[i+lb2]); 00367 } else { 00368 x_inf[i+1] = -Sup(x[i+lb1]); 00369 y_inf[i+1] = Sup(y[i+lb2]); 00370 x_sup[i+1] = Inf(x[i+lb1]); 00371 y_sup[i+1] = Inf(y[i+lb2]); 00372 } 00373 00374 } 00375 00376 } 00377 } 00378 00379 //Sorts inf and sup of an interval vector and a real vector into two real vector for the computation of the infimum 00380 //and two real vectors for the computation of the supremum. Rounding mode must be set to upwards! 00381 inline void bsort(const ivector &x, const rvector &y, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb1, int lb2) { 00382 00383 for(int i=0 ; i<n ; i++) { 00384 00385 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) { 00386 00387 if(y[i+lb2] >= 0) { 00388 x_inf[i+1] = -Inf(x[i+lb1]); 00389 y_inf[i+1] = y[i+lb2]; 00390 x_sup[i+1] = Sup(x[i+lb1]); 00391 y_sup[i+1] = y[i+lb2]; 00392 } else { 00393 x_inf[i+1] = -Sup(x[i+lb1]); 00394 y_inf[i+1] = y[i+lb2]; 00395 x_sup[i+1] = Inf(x[i+lb1]); 00396 y_sup[i+1] = y[i+lb2]; 00397 } 00398 00399 } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) { 00400 00401 if(y[i+lb2] >= 0) { 00402 x_inf[i+1] = -Inf(x[i+lb1]); 00403 y_inf[i+1] = y[i+lb2]; 00404 x_sup[i+1] = Sup(x[i+lb1]); 00405 y_sup[i+1] = y[i+lb2]; 00406 } else { 00407 x_inf[i+1] = -Sup(x[i+lb1]); 00408 y_inf[i+1] = y[i+lb2]; 00409 x_sup[i+1] = Inf(x[i+lb1]); 00410 y_sup[i+1] = y[i+lb2]; 00411 } 00412 00413 } else { 00414 00415 if(y[i+lb2] >= 0) { 00416 x_inf[i+1] = -Inf(x[i+lb1]); 00417 y_inf[i+1] = y[i+lb2]; 00418 x_sup[i+1] = Sup(x[i+lb1]); 00419 y_sup[i+1] = y[i+lb2]; 00420 } else { 00421 x_inf[i+1] = -Sup(x[i+lb1]); 00422 y_inf[i+1] = y[i+lb2]; 00423 x_sup[i+1] = Inf(x[i+lb1]); 00424 y_sup[i+1] = y[i+lb2]; 00425 } 00426 00427 } 00428 00429 } 00430 } 00431 00432 //Sorts inf and sup of an interval vector and a real vector into two real vector for the computation of the infimum 00433 //and two real vectors for the computation of the supremum. Rounding mode must be set to upwards! 00434 inline void bsort(const rvector &y, const ivector &x, rvector& x_inf, rvector& y_inf, rvector &x_sup, rvector &y_sup, int n, int lb2, int lb1) { 00435 for(int i=0 ; i<n ; i++) { 00436 00437 if(Inf(x[i+lb1]) >= 0 && Sup(x[i+lb1]) >= 0) { 00438 00439 if(y[i+lb2] >= 0) { 00440 x_inf[i+1] = -Inf(x[i+lb1]); 00441 y_inf[i+1] = y[i+lb2]; 00442 x_sup[i+1] = Sup(x[i+lb1]); 00443 y_sup[i+1] = y[i+lb2]; 00444 } else { 00445 x_inf[i+1] = -Sup(x[i+lb1]); 00446 y_inf[i+1] = y[i+lb2]; 00447 x_sup[i+1] = Inf(x[i+lb1]); 00448 y_sup[i+1] = y[i+lb2]; 00449 } 00450 00451 } else if(Inf(x[i+lb1]) < 0 && Sup(x[i+lb1]) >= 0) { 00452 00453 if(y[i+lb2] >= 0) { 00454 x_inf[i+1] = -Inf(x[i+lb1]); 00455 y_inf[i+1] = y[i+lb2]; 00456 x_sup[i+1] = Sup(x[i+lb1]); 00457 y_sup[i+1] = y[i+lb2]; 00458 } else { 00459 x_inf[i+1] = -Sup(x[i+lb1]); 00460 y_inf[i+1] = y[i+lb2]; 00461 x_sup[i+1] = Inf(x[i+lb1]); 00462 y_sup[i+1] = y[i+lb2]; 00463 } 00464 00465 } else { 00466 00467 if(y[i+lb2] >= 0) { 00468 x_inf[i+1] = -Inf(x[i+lb1]); 00469 y_inf[i+1] = y[i+lb2]; 00470 x_sup[i+1] = Sup(x[i+lb1]); 00471 y_sup[i+1] = y[i+lb2]; 00472 } else { 00473 x_inf[i+1] = -Sup(x[i+lb1]); 00474 y_inf[i+1] = y[i+lb2]; 00475 x_sup[i+1] = Inf(x[i+lb1]); 00476 y_sup[i+1] = y[i+lb2]; 00477 } 00478 00479 } 00480 00481 } 00482 } 00483 00484 inline void blasdot(const ivector& x, const ivector& y, interval& res) { 00485 const int n = VecLen(x); 00486 const int inc=1; 00487 const int lb1 = Lb(x); 00488 const int lb2 = Lb(y); 00489 int rnd = getround(); 00490 00491 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00492 00493 double* dxi = x_inf.to_blas_array(); 00494 double* dxs = x_sup.to_blas_array(); 00495 double* dyi = y_inf.to_blas_array(); 00496 double* dys = y_sup.to_blas_array(); 00497 double res_inf,res_sup; 00498 00499 setround(1); 00500 00501 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00502 00503 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00504 00505 //setround(-1); 00506 00507 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00508 00509 setround(rnd); 00510 00511 res = interval(res_inf,res_sup); 00512 } 00513 00514 inline void blasdot(const ivector& x, const rvector& y, interval& res) { 00515 const int n = VecLen(x); 00516 const int inc=1; 00517 const int lb1 = Lb(x); 00518 const int lb2 = Lb(y); 00519 int rnd = getround(); 00520 00521 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00522 00523 double* dxi = x_inf.to_blas_array(); 00524 double* dxs = x_sup.to_blas_array(); 00525 double* dyi = y_inf.to_blas_array(); 00526 double* dys = y_sup.to_blas_array(); 00527 double res_inf,res_sup; 00528 00529 setround(1); 00530 00531 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00532 00533 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00534 00535 //setround(-1); 00536 00537 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00538 00539 setround(rnd); 00540 00541 res = interval(res_inf,res_sup); 00542 } 00543 00544 inline void blasdot(const rvector& x, const ivector& y, interval& res) { 00545 const int n = VecLen(x); 00546 const int inc=1; 00547 const int lb1 = Lb(x); 00548 const int lb2 = Lb(y); 00549 int rnd = getround(); 00550 00551 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00552 00553 double* dxi = x_inf.to_blas_array(); 00554 double* dxs = x_sup.to_blas_array(); 00555 double* dyi = y_inf.to_blas_array(); 00556 double* dys = y_sup.to_blas_array(); 00557 double res_inf,res_sup; 00558 00559 setround(1); 00560 00561 bsort(x,y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00562 00563 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00564 00565 //setround(-1); 00566 00567 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00568 00569 setround(rnd); 00570 00571 res = interval(res_inf,res_sup); 00572 } 00573 #endif 00574 #endif 00575 00576 00577 #ifdef _CXSC_BLAS_CIVECTOR 00578 #ifndef _CXSC_BLAS_CIVECTOR_INC 00579 #define _CXSC_BLAS_CIVECTOR_INC 00580 inline void blasdot(const rvector& x, const civector& y, cinterval& res) { 00581 const int n = VecLen(x); 00582 const int inc=1; 00583 const int lb1 = Lb(x); 00584 const int lb2 = Lb(y); 00585 int rnd = getround(); 00586 00587 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00588 00589 double* dxi = x_inf.to_blas_array(); 00590 double* dxs = x_sup.to_blas_array(); 00591 double* dyi = y_inf.to_blas_array(); 00592 double* dys = y_sup.to_blas_array(); 00593 double res_inf,res_sup; 00594 00595 setround(1); 00596 00597 bsort(x,Re(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00598 00599 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00600 00601 //setround(-1); 00602 00603 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00604 00605 SetRe(res, interval(res_inf,res_sup)); 00606 00607 bsort(x,Im(y),x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00608 00609 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00610 00611 //setround(-1); 00612 00613 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00614 00615 SetIm(res, interval(res_inf,res_sup)); 00616 00617 setround(rnd); 00618 } 00619 00620 00621 inline void blasdot(const civector& x, const rvector& y, cinterval& res) { 00622 const int n = VecLen(x); 00623 const int inc=1; 00624 const int lb1 = Lb(x); 00625 const int lb2 = Lb(y); 00626 int rnd = getround(); 00627 00628 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00629 00630 double* dxi = x_inf.to_blas_array(); 00631 double* dxs = x_sup.to_blas_array(); 00632 double* dyi = y_inf.to_blas_array(); 00633 double* dys = y_sup.to_blas_array(); 00634 double res_inf,res_sup; 00635 00636 setround(1); 00637 00638 bsort(Re(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00639 00640 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00641 00642 //setround(-1); 00643 00644 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00645 00646 SetRe(res, interval(res_inf,res_sup)); 00647 00648 bsort(Im(x),y,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00649 00650 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00651 00652 //setround(-1); 00653 00654 res_inf = -cblas_ddot(n, dxi, inc, dyi, inc); 00655 00656 SetIm(res, interval(res_inf,res_sup)); 00657 00658 setround(rnd); 00659 } 00660 00661 inline void blasdot(const civector& x, const civector& y, cinterval& res) { 00662 const int n = VecLen(x); 00663 const int inc=1; 00664 const int lb1 = Lb(x); 00665 const int lb2 = Lb(y); 00666 int rnd = getround(); 00667 00668 cvector xmid(n),xrad(n),ymid(n),yrad(n); 00669 real crad_re = 0.0, crad_im = 0.0; 00670 complex res_inf(0.0,0.0),res_sup(0.0,0.0); 00671 real tmp1,tmp2,tmp3,tmp4; 00672 00673 double* dxmid = xmid.to_blas_array(); 00674 double* dymid = ymid.to_blas_array(); 00675 00676 setround(1); 00677 00678 for(int i=0 ; i<n ; i++) { 00679 xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i])); 00680 ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i])); 00681 xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]); 00682 yrad[i+1] = ymid[i+1] - Inf(y[i+lb2]); 00683 00684 tmp1 = abs(Re(ymid[i+1])) + Re(yrad[i+1]); 00685 tmp2 = abs(Im(ymid[i+1])) + Im(yrad[i+1]); 00686 tmp3 = abs(Re(xmid[i+1])); 00687 tmp4 = abs(Im(xmid[i+1])); 00688 00689 crad_re += Re(xrad[i+1]) * tmp1 + tmp3 * Re(yrad[i+1]); 00690 crad_re += Im(xrad[i+1]) * tmp2 + tmp4 * Im(yrad[i+1]); 00691 00692 crad_im += Re(xrad[i+1]) * tmp2 + tmp3 * Im(yrad[i+1]); 00693 crad_im += Im(xrad[i+1]) * tmp1 + tmp4 * Re(yrad[i+1]); 00694 } 00695 00696 cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_sup); 00697 00698 res_sup += complex(crad_re,crad_im); 00699 00700 setround(-1); 00701 00702 cblas_zdotu_sub(n, dxmid, inc, dymid, inc, &res_inf); 00703 00704 res_inf -= complex(crad_re,crad_im); 00705 00706 setround(rnd); 00707 00708 UncheckedSetInf(res,res_inf); 00709 UncheckedSetSup(res,res_sup); 00710 00711 } 00712 00713 00714 inline void blasdot(const civector& x, const cvector& y, cinterval& res) { 00715 const int n = VecLen(x); 00716 const int inc=1; 00717 const int lb1 = Lb(x); 00718 int rnd = getround(); 00719 00720 cvector xmid(n),xrad(n); 00721 real crad_re = 0.0, crad_im = 0.0; 00722 complex res_inf(0.0,0.0),res_sup(0.0,0.0); 00723 real tmp1,tmp2,tmp3,tmp4; 00724 00725 double* dxmid = xmid.to_blas_array(); 00726 double* dy = y.to_blas_array(); 00727 00728 setround(1); 00729 00730 for(int i=0 ; i<n ; i++) { 00731 xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i])); 00732 xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]); 00733 00734 tmp1 = abs(Re(y[i+1])); 00735 tmp2 = abs(Im(y[i+1])); 00736 00737 crad_re += Re(xrad[i+1]) * tmp1; 00738 crad_re += Im(xrad[i+1]) * tmp2; 00739 00740 crad_im += Re(xrad[i+1]) * tmp2; 00741 crad_im += Im(xrad[i+1]) * tmp1; 00742 } 00743 00744 cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_sup); 00745 00746 res_sup += complex(crad_re,crad_im); 00747 00748 setround(-1); 00749 00750 cblas_zdotu_sub(n, dxmid, inc, dy, inc, &res_inf); 00751 00752 res_inf -= complex(crad_re,crad_im); 00753 00754 setround(rnd); 00755 00756 UncheckedSetInf(res,res_inf); 00757 UncheckedSetSup(res,res_sup); 00758 00759 } 00760 00761 inline void blasdot(const cvector& x, const civector& y, cinterval& res) { 00762 const int n = VecLen(x); 00763 const int inc=1; 00764 const int lb1 = Lb(x); 00765 int rnd = getround(); 00766 00767 cvector ymid(n),yrad(n); 00768 real crad_re = 0.0, crad_im = 0.0; 00769 complex res_inf(0.0,0.0),res_sup(0.0,0.0); 00770 real tmp1,tmp2,tmp3,tmp4; 00771 00772 double* dx = x.to_blas_array(); 00773 double* dymid = ymid.to_blas_array(); 00774 00775 setround(1); 00776 00777 for(int i=0 ; i<n ; i++) { 00778 ymid[i+1] = Inf(y[lb1+i]) + 0.5*(Sup(y[lb1+i])-Inf(y[lb1+i])); 00779 yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]); 00780 00781 tmp1 = abs(Re(x[i+1])); 00782 tmp2 = abs(Im(x[i+1])); 00783 00784 crad_re += Re(yrad[i+1]) * tmp1; 00785 crad_re += Im(yrad[i+1]) * tmp2; 00786 00787 crad_im += Re(yrad[i+1]) * tmp2; 00788 crad_im += Im(yrad[i+1]) * tmp1; 00789 } 00790 00791 cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_sup); 00792 00793 res_sup += complex(crad_re,crad_im); 00794 00795 setround(-1); 00796 00797 cblas_zdotu_sub(n, dx, inc, dymid, inc, &res_inf); 00798 00799 res_inf -= complex(crad_re,crad_im); 00800 00801 setround(rnd); 00802 00803 UncheckedSetInf(res,res_inf); 00804 UncheckedSetSup(res,res_sup); 00805 00806 } 00807 00808 00809 inline void blasdot(const civector& x, const ivector& y, cinterval& res) { 00810 const int n = VecLen(x); 00811 const int inc=1, inc2=2; 00812 const int lb1 = Lb(x); 00813 const int lb2 = Lb(y); 00814 int rnd = getround(); 00815 00816 cvector xmid(n),xrad(n); 00817 rvector ymid(n),yrad(n); 00818 real crad_re = 0.0, crad_im = 0.0; 00819 complex res_inf(0.0,0.0),res_sup(0.0,0.0); 00820 real tmp1,tmp2,tmp3,tmp4; 00821 00822 double* dxmid = xmid.to_blas_array(); 00823 double* dymid = ymid.to_blas_array(); 00824 00825 setround(1); 00826 00827 for(int i=0 ; i<n ; i++) { 00828 xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i])); 00829 ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i])); 00830 xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]); 00831 yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]); 00832 00833 tmp1 = abs(ymid[i+1]) + yrad[i+1]; 00834 tmp2 = abs(ymid[i+1]) + yrad[i+1]; 00835 tmp3 = abs(Re(xmid[i+1])); 00836 tmp4 = abs(Im(xmid[i+1])); 00837 00838 crad_re += Re(xrad[i+1]) * tmp1 + tmp3 * yrad[i+1]; 00839 crad_re += Im(xrad[i+1]) * tmp2 + tmp4 * yrad[i+1]; 00840 00841 crad_im += Re(xrad[i+1]) * tmp2 + tmp3 * yrad[i+1]; 00842 crad_im += Im(xrad[i+1]) * tmp1 + tmp4 * yrad[i+1]; 00843 } 00844 00845 SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc)); 00846 SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc)); 00847 00848 res_sup += complex(crad_re,crad_im); 00849 00850 setround(-1); 00851 00852 SetRe(res_sup, cblas_ddot(n, dxmid, inc2, dymid, inc)); 00853 SetIm(res_sup, cblas_ddot(n, dxmid+1, inc2, dymid, inc)); 00854 00855 res_inf -= complex(crad_re,crad_im); 00856 00857 setround(rnd); 00858 00859 UncheckedSetInf(res,res_inf); 00860 UncheckedSetSup(res,res_sup); 00861 } 00862 00863 00864 inline void blasdot(const ivector& x, const civector& y, cinterval& res) { 00865 const int n = VecLen(x); 00866 const int inc=1, inc2=2; 00867 const int lb1 = Lb(x); 00868 const int lb2 = Lb(y); 00869 int rnd = getround(); 00870 00871 rvector xmid(n),xrad(n); 00872 cvector ymid(n),yrad(n); 00873 real crad_re = 0.0, crad_im = 0.0; 00874 complex res_inf(0.0,0.0),res_sup(0.0,0.0); 00875 real tmp1,tmp2,tmp3,tmp4; 00876 00877 double* dxmid = xmid.to_blas_array(); 00878 double* dymid = ymid.to_blas_array(); 00879 00880 setround(1); 00881 00882 for(int i=0 ; i<n ; i++) { 00883 xmid[i+1] = Inf(x[lb1+i]) + 0.5*(Sup(x[lb1+i])-Inf(x[lb1+i])); 00884 ymid[i+1] = Inf(y[lb2+i]) + 0.5*(Sup(y[lb2+i])-Inf(y[lb2+i])); 00885 xrad[i+1] = xmid[i+1] - Inf(x[i+lb1]); 00886 yrad[i+1] = ymid[i+1] - Inf(y[i+lb1]); 00887 00888 tmp1 = abs(Re(ymid[i+1])) + Re(yrad[i+1]); 00889 tmp2 = abs(Im(ymid[i+1])) + Im(yrad[i+1]); 00890 tmp3 = abs(xmid[i+1]); 00891 tmp4 = abs(xmid[i+1]); 00892 00893 crad_re += xrad[i+1] * tmp1 + tmp3 * Re(yrad[i+1]); 00894 crad_re += xrad[i+1] * tmp2 + tmp4 * Im(yrad[i+1]); 00895 00896 crad_im += xrad[i+1] * tmp2 + tmp3 * Im(yrad[i+1]); 00897 crad_im += xrad[i+1] * tmp1 + tmp4 * Re(yrad[i+1]); 00898 } 00899 00900 SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2)); 00901 SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2)); 00902 00903 res_sup += complex(crad_re,crad_im); 00904 00905 setround(-1); 00906 00907 SetRe(res_sup, cblas_ddot(n, dxmid, inc, dymid, inc2)); 00908 SetIm(res_sup, cblas_ddot(n, dxmid, inc, dymid+1, inc2)); 00909 00910 res_inf -= complex(crad_re,crad_im); 00911 00912 setround(rnd); 00913 00914 UncheckedSetInf(res,res_inf); 00915 UncheckedSetSup(res,res_sup); 00916 } 00917 #endif 00918 #endif 00919 /***************************************************************************/ 00920 00921 00922 #ifdef _CXSC_BLAS_RMATRIX 00923 #ifndef _CXSC_BLAS_RMATVEC_INC 00924 #define _CXSC_BLAS_RMATVEC_INC 00925 inline void blasmvmul(const rmatrix& A, const rvector& x, rvector& r) { 00926 00927 double* DA = (double*) A.to_blas_array(); 00928 double* dx = (double*) x.to_blas_array(); 00929 double* dr = (double*) r.to_blas_array(); 00930 00931 const int m = Ub(A,1) - Lb(A,1) + 1; 00932 const int n = Ub(A,2) - Lb(A,2) + 1; 00933 const double alpha = 1.0, beta = 0.0; 00934 const int inc = 1; 00935 00936 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 00937 alpha, DA, n, dx, inc, beta, dr, inc); 00938 } 00939 #endif 00940 #endif 00941 00942 #ifdef _CXSC_BLAS_RMATRIX 00943 #ifdef _CXSC_BLAS_IVECTOR 00944 #ifndef _CXSC_BLAS_RIMATVEC_INC 00945 #define _CXSC_BLAS_RIMATVEC_INC 00946 inline void blasmvmul(const rmatrix& A, const rvector& x, ivector& r) { 00947 int rnd = getround(); 00948 double* DA = (double*) A.to_blas_array(); 00949 double* dx = (double*) x.to_blas_array(); 00950 00951 const int m = Ub(A,1) - Lb(A,1) + 1; 00952 const int n = Ub(A,2) - Lb(A,2) + 1; 00953 const double alpha = 1.0, beta = 0.0; 00954 const int inc = 1; 00955 double* dr = new double[n]; 00956 00957 setround(1); 00958 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 00959 alpha, DA, n, dx, inc, beta, dr, inc); 00960 for(int i=0 ; i<n ; i++) 00961 UncheckedSetSup(r[i+Lb(r)], dr[i]); 00962 00963 setround(-1); 00964 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 00965 alpha, DA, n, dx, inc, beta, dr, inc); 00966 for(int i=0 ; i<n ; i++) 00967 UncheckedSetInf(r[i+Lb(r)], dr[i]); 00968 00969 delete[] dr; 00970 setround(rnd); 00971 } 00972 00973 inline void blasmvmul(const rmatrix& A, const ivector& x, ivector& r) { 00974 const int n = VecLen(x); 00975 const int m = Ub(A,1)-Lb(A,1)+1; 00976 const int inc=1; 00977 const int lb1 = Lb(A,2); 00978 const int lb2 = Lb(x); 00979 int rnd = getround(); 00980 00981 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 00982 00983 double* dxi = x_inf.to_blas_array(); 00984 double* dxs = x_sup.to_blas_array(); 00985 double* dyi = y_inf.to_blas_array(); 00986 double* dys = y_sup.to_blas_array(); 00987 double res_inf,res_sup; 00988 00989 setround(1); 00990 00991 for(int i=1 ; i<=m ; i++) { 00992 bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 00993 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 00994 UncheckedSetSup(r[i], res_sup); 00995 // setround(-1); 00996 res_inf = cblas_ddot(n, dxi, inc, dyi, inc); 00997 UncheckedSetInf(r[i], -res_inf); 00998 } 00999 01000 setround(rnd); 01001 } 01002 #endif 01003 #endif 01004 #endif 01005 01006 #ifdef _CXSC_BLAS_IMATRIX 01007 #ifdef _CXSC_BLAS_IVECTOR 01008 #ifndef _CXSC_BLAS_IRMATVEC_INC 01009 #define _CXSC_BLAS_IRMATVEC_INC 01010 inline void blasmvmul(const imatrix& A, const rvector& x, ivector& r) { 01011 const int n = VecLen(x); 01012 const int m = Ub(A,1)-Lb(A,1)+1; 01013 const int inc=1; 01014 const int lb1 = Lb(A,2); 01015 const int lb2 = Lb(x); 01016 int rnd = getround(); 01017 01018 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 01019 01020 double* dxi = x_inf.to_blas_array(); 01021 double* dxs = x_sup.to_blas_array(); 01022 double* dyi = y_inf.to_blas_array(); 01023 double* dys = y_sup.to_blas_array(); 01024 double res_inf,res_sup; 01025 01026 setround(1); 01027 01028 for(int i=1 ; i<=m ; i++) { 01029 bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 01030 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 01031 UncheckedSetSup(r[i], res_sup); 01032 // setround(-1); 01033 res_inf = cblas_ddot(n, dxi, inc, dyi, inc); 01034 UncheckedSetInf(r[i], -res_inf); 01035 } 01036 01037 setround(rnd); 01038 01039 } 01040 #endif 01041 #endif 01042 #endif 01043 01044 #ifdef _CXSC_BLAS_IMATRIX 01045 #ifdef _CXSC_BLAS_IVECTOR 01046 #ifndef _CXSC_BLAS_IIMATVEC_INC 01047 #define _CXSC_BLAS_IIMATVEC_INC 01048 inline void blasmvmul(const imatrix& A, const ivector& x, ivector& r) { 01049 const int n = VecLen(x); 01050 const int m = Ub(A,1)-Lb(A,1)+1; 01051 const int inc=1; 01052 const int lb1 = Lb(A,2); 01053 const int lb2 = Lb(x); 01054 int rnd = getround(); 01055 01056 rvector x_inf(n),x_sup(n),y_inf(n),y_sup(n); 01057 01058 double* dxi = x_inf.to_blas_array(); 01059 double* dxs = x_sup.to_blas_array(); 01060 double* dyi = y_inf.to_blas_array(); 01061 double* dys = y_sup.to_blas_array(); 01062 double res_inf,res_sup; 01063 01064 setround(1); 01065 01066 for(int i=1 ; i<=m ; i++) { 01067 bsort(A[i+Lb(A,1)-1],x,x_inf,y_inf,x_sup,y_sup,n,lb1,lb2); 01068 res_sup = cblas_ddot(n, dxs, inc, dys, inc); 01069 UncheckedSetSup(r[i], res_sup); 01070 // setround(-1); 01071 res_inf = cblas_ddot(n, dxi, inc, dyi, inc); 01072 UncheckedSetInf(r[i], -res_inf); 01073 } 01074 01075 setround(rnd); 01076 01077 } 01078 #endif 01079 #endif 01080 #endif 01081 01082 #ifdef _CXSC_BLAS_CMATRIX 01083 #ifdef _CXSC_BLAS_CIVECTOR 01084 #ifndef _CXSC_BLAS_CIMATVEC_INC 01085 #define _CXSC_BLAS_CIMATVEC_INC 01086 inline void blasmvmul(const cmatrix& A, const ivector& x, civector& r) { 01087 const int m = Ub(A,1)-Lb(A,1)+1; 01088 01089 for(int i=0 ; i<m ; i++) 01090 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01091 } 01092 #endif 01093 #endif 01094 #endif 01095 01096 #ifdef _CXSC_BLAS_IMATRIX 01097 #ifdef _CXSC_BLAS_CIVECTOR 01098 #ifndef _CXSC_BLAS_ICMATVEC_INC 01099 #define _CXSC_BLAS_ICMATVEC_INC 01100 inline void blasmvmul(const imatrix& A, const cvector& x, civector& r) { 01101 const int m = Ub(A,1)-Lb(A,1)+1; 01102 01103 for(int i=0 ; i<m ; i++) 01104 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01105 } 01106 #endif 01107 #endif 01108 #endif 01109 01110 #ifdef _CXSC_BLAS_RMATRIX 01111 #ifdef _CXSC_BLAS_CIVECTOR 01112 #ifndef _CXSC_BLAS_RCIMATVEC_INC 01113 #define _CXSC_BLAS_RCIMATVEC_INC 01114 inline void blasmvmul(const rmatrix& A, const civector& x, civector& r) { 01115 const int m = Ub(A,1)-Lb(A,1)+1; 01116 01117 for(int i=0 ; i<m ; i++) 01118 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01119 } 01120 #endif 01121 #endif 01122 #endif 01123 01124 #ifdef _CXSC_BLAS_CIMATRIX 01125 #ifdef _CXSC_BLAS_CIVECTOR 01126 #ifndef _CXSC_BLAS_CIRMATVEC_INC 01127 #define _CXSC_BLAS_CIRMATVEC_INC 01128 inline void blasmvmul(const cimatrix& A, const rvector& x, civector& r) { 01129 const int m = Ub(A,1)-Lb(A,1)+1; 01130 01131 for(int i=0 ; i<m ; i++) 01132 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01133 } 01134 #endif 01135 #endif 01136 #endif 01137 01138 #ifdef _CXSC_BLAS_CMATRIX 01139 #ifdef _CXSC_BLAS_CIVECTOR 01140 #ifndef _CXSC_BLAS_CCIMATVEC_INC 01141 #define _CXSC_BLAS_CCIMATVEC_INC 01142 inline void blasmvmul(const cmatrix& A, const civector& x, civector& r) { 01143 const int m = Ub(A,1)-Lb(A,1)+1; 01144 01145 for(int i=0 ; i<m ; i++) 01146 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01147 } 01148 01149 inline void blasmvmul(const cmatrix& A, const cvector& x, civector& r) { 01150 ivector tmp1, tmp2; 01151 blasmvmul(Re(A),Re(x),tmp1); 01152 blasmvmul(Im(A),Im(x),tmp2); 01153 SetRe(r,tmp1-tmp2); 01154 blasmvmul(Re(A),Im(x),tmp1); 01155 blasmvmul(Im(A),Re(x),tmp2); 01156 SetIm(r,tmp1+tmp2); 01157 } 01158 01159 inline void blasmvmul(const cmatrix& A, const rvector& x, civector& r) { 01160 int n = VecLen(r); 01161 ivector re(n),im(n); 01162 01163 blasmvmul(Re(A),x,re); 01164 blasmvmul(Im(A),x,im); 01165 01166 SetRe(r,re); 01167 SetIm(r,im); 01168 } 01169 01170 inline void blasmvmul(const rmatrix& A, const cvector& x, civector& r) { 01171 int n = VecLen(r); 01172 ivector re(n),im(n); 01173 01174 blasmvmul(A,Re(x),re); 01175 blasmvmul(A,Im(x),im); 01176 01177 SetRe(r,re); 01178 SetIm(r,im); 01179 } 01180 #endif 01181 #endif 01182 #endif 01183 01184 #ifdef _CXSC_BLAS_CIMATRIX 01185 #ifdef _CXSC_BLAS_CIVECTOR 01186 #ifndef _CXSC_BLAS_CICMATVEC_INC 01187 #define _CXSC_BLAS_CICMATVEC_INC 01188 inline void blasmvmul(const cimatrix& A, const cvector& x, civector& r) { 01189 const int m = Ub(A,1)-Lb(A,1)+1; 01190 01191 for(int i=0 ; i<m ; i++) 01192 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01193 } 01194 #endif 01195 #endif 01196 #endif 01197 01198 #ifdef _CXSC_BLAS_IMATRIX 01199 #ifdef _CXSC_BLAS_CIVECTOR 01200 #ifndef _CXSC_BLAS_ICIMATVEC_INC 01201 #define _CXSC_BLAS_ICIMATVEC_INC 01202 inline void blasmvmul(const imatrix& A, const civector& x, civector& r) { 01203 const int m = Ub(A,1)-Lb(A,1)+1; 01204 01205 for(int i=0 ; i<m ; i++) 01206 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01207 } 01208 #endif 01209 #endif 01210 #endif 01211 01212 #ifdef _CXSC_BLAS_CIMATRIX 01213 #ifdef _CXSC_BLAS_CIVECTOR 01214 #ifndef _CXSC_BLAS_CIIMATVEC_INC 01215 #define _CXSC_BLAS_CIIMATVEC_INC 01216 inline void blasmvmul(const cimatrix& A, const ivector& x, civector& r) { 01217 const int m = Ub(A,1)-Lb(A,1)+1; 01218 01219 for(int i=0 ; i<m ; i++) 01220 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01221 } 01222 #endif 01223 #endif 01224 #endif 01225 01226 #ifdef _CXSC_BLAS_CIMATRIX 01227 #ifdef _CXSC_BLAS_CIVECTOR 01228 #ifndef _CXSC_BLAS_CICIMATVEC_INC 01229 #define _CXSC_BLAS_CICIMATVEC_INC 01230 inline void blasmvmul(const cimatrix& A, const civector& x, civector& r) { 01231 const int m = Ub(A,1)-Lb(A,1)+1; 01232 01233 for(int i=0 ; i<m ; i++) 01234 blasdot(A[i+Lb(A,1)], x, r[i+Lb(r)]); 01235 } 01236 #endif 01237 #endif 01238 #endif 01239 01240 #ifdef _CXSC_BLAS_CMATRIX 01241 #ifdef _CXSC_BLAS_CVECTOR 01242 #ifndef _CXSC_BLAS_CCMATVEC_INC 01243 #define _CXSC_BLAS_CCMATVEC_INC 01244 inline void blasmvmul(const cmatrix& A, const cvector& x, cvector& r) { 01245 01246 double* DA = (double*) A.to_blas_array(); 01247 double* dx = (double*) x.to_blas_array(); 01248 double* dr = (double*) r.to_blas_array(); 01249 01250 const int m = Ub(A,1) - Lb(A,1) + 1; 01251 const int n = Ub(A,2) - Lb(A,2) + 1; 01252 const complex alpha(1.0,0.0); 01253 const complex beta(0.0,0.0); 01254 const int inc = 1; 01255 01256 cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n, 01257 &alpha, DA, n, dx, inc, &beta, dr, inc); 01258 } 01259 #endif 01260 #endif 01261 #endif 01262 01263 #ifdef _CXSC_BLAS_RMATRIX 01264 #ifdef _CXSC_BLAS_CVECTOR 01265 #ifndef _CXSC_BLAS_RCMATVEC_INC 01266 #define _CXSC_BLAS_RCMATVEC_INC 01267 inline void blasmvmul(const rmatrix& A, const cvector& x, cvector& r) { 01268 01269 double* DA = (double*) A.to_blas_array(); 01270 double* dx = (double*) x.to_blas_array(); 01271 double* dr = (double*) r.to_blas_array(); 01272 01273 const int m = Ub(A,1) - Lb(A,1) + 1; 01274 const int n = Ub(A,2) - Lb(A,2) + 1; 01275 const double alpha = 1.0; 01276 const double beta = 0.0; 01277 const int inc = 2; 01278 01279 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 01280 alpha, DA, m, dx, inc, beta, dr, inc); 01281 01282 cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 01283 alpha, DA, n, dx+1, inc, beta, dr+1, inc); 01284 } 01285 #endif 01286 #endif 01287 #endif 01288 01289 #ifdef _CXSC_BLAS_CMATRIX 01290 #ifdef _CXSC_BLAS_CVECTOR 01291 #ifndef _CXSC_BLAS_CRMATVEC_INC 01292 #define _CXSC_BLAS_CRMATVEC_INC 01293 inline void blasmvmul(const cmatrix& A, const rvector& x, cvector& r) { 01294 cvector tmp(x); 01295 01296 double* DA = (double*) A.to_blas_array(); 01297 double* dx = (double*) tmp.to_blas_array(); 01298 double* dr = (double*) r.to_blas_array(); 01299 01300 const int m = Ub(A,1) - Lb(A,1) + 1; 01301 const int n = Ub(A,2) - Lb(A,2) + 1; 01302 const complex alpha(1.0,0.0); 01303 const complex beta(0.0,0.0); 01304 const int inc = 1; 01305 01306 cblas_zgemv(CblasRowMajor, CblasNoTrans, m, n, 01307 &alpha, DA, n, dx, inc, &beta, dr, inc); 01308 } 01309 #endif 01310 #endif 01311 #endif 01312 01313 /***************************************************************************/ 01314 01315 #ifdef _CXSC_BLAS_RMATRIX 01316 #ifndef _CXSC_BLAS_RMATRIX_INC 01317 #define _CXSC_BLAS_RMATRIX_INC 01318 inline void blasmatmul(const rmatrix &A, const rmatrix &B, rmatrix &C) { 01319 01320 double* DA = (double*) A.to_blas_array(); 01321 double* DB = (double*) B.to_blas_array(); 01322 double* DC = (double*) C.to_blas_array(); 01323 01324 const int m = Ub(A,1) - Lb(A,1) + 1; 01325 const int n = Ub(A,2) - Lb(A,2) + 1; 01326 const int o = Ub(C,2) - Lb(C,2) + 1; 01327 const double alpha = 1.0, beta = 0.0; 01328 01329 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01330 n, alpha, DA, n, DB, o, beta, DC, o); 01331 } 01332 #endif 01333 #endif 01334 01335 #ifdef _CXSC_BLAS_IMATRIX 01336 #ifndef _CXSC_BLAS_IMATRIX_INC 01337 #define _CXSC_BLAS_IMATRIX_INC 01338 inline void blasmatmul(const rmatrix &A, const rmatrix &B, imatrix &C) { 01339 int rnd = getround(); 01340 double* DA = (double*) A.to_blas_array(); 01341 double* DB = (double*) B.to_blas_array(); 01342 int ind; 01343 01344 const int m = Ub(A,1) - Lb(A,1) + 1; 01345 const int n = Ub(A,2) - Lb(A,2) + 1; 01346 const int o = Ub(C,2) - Lb(C,2) + 1; 01347 double* DC = new double[m*o]; 01348 const double alpha = 1.0, beta = 0.0; 01349 01350 setround(1); 01351 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01352 n, alpha, DA, n, DB, o, beta, DC, o); 01353 01354 for(int i=0 ; i<m ; i++) { 01355 for(int j=0 ; j<o ; j++) { 01356 ind = i*o+j; 01357 UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], DC[ind]); 01358 } 01359 } 01360 01361 setround(-1); 01362 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01363 n, alpha, DA, n, DB, o, beta, DC, o); 01364 01365 for(int i=0 ; i<m ; i++) { 01366 for(int j=0 ; j<o ; j++) { 01367 ind = i*o+j; 01368 UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], DC[ind]); 01369 } 01370 } 01371 01372 delete[] DC; 01373 01374 setround(rnd); 01375 } 01376 01377 01378 inline void blasmatmul(const imatrix &A, const imatrix &B, imatrix &C) { 01379 const int m = Ub(A,1) - Lb(A,1) + 1; 01380 const int n = Ub(A,2) - Lb(A,2) + 1; 01381 const int o = Ub(C,2) - Lb(C,2) + 1; 01382 01383 int lbC_1 = Lb(C,1); int lbC_2 = Lb(C,2); 01384 Resize(C); 01385 01386 double* Amid = new double[m*n]; 01387 double* Aabs = new double[m*n]; 01388 double* Arad = new double[m*n]; 01389 double* Bmid = new double[n*o]; 01390 double* Babs = new double[n*o]; 01391 double* Brad = new double[n*o]; 01392 double* C1 = new double[m*o]; 01393 double* C2 = new double[m*o]; 01394 int rnd = getround(); 01395 01396 const double alpha = 1.0, beta = 0.0; 01397 01398 setround(1); 01399 01400 int ind; 01401 01402 //Compute mid and rad of A 01403 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) { 01404 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) { 01405 ind = (i-Lb(A,1))*n+(j-Lb(A,2)); 01406 Amid[ind] = _double(Inf(A[i][j]) + 0.5*(Sup(A[i][j]) - Inf(A[i][j]))); 01407 Arad[ind] = _double(Amid[ind] - Inf(A[i][j])); 01408 Aabs[ind] = fabs(Amid[ind]); 01409 } 01410 } 01411 01412 //Compute mid and rad of B 01413 for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) { 01414 for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) { 01415 ind = (i-Lb(B,1))*o+(j-Lb(B,2)); 01416 Bmid[ind] = _double(Inf(B[i][j]) + 0.5*(Sup(B[i][j]) - Inf(B[i][j]))); 01417 Brad[ind] = _double(Bmid[ind] - Inf(B[i][j])); 01418 Babs[ind] = fabs(Bmid[ind]); 01419 } 01420 } 01421 01422 setround(-1); 01423 01424 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01425 n, alpha, Amid, n, Bmid, o, beta, C1, o); 01426 01427 setround(1); 01428 01429 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01430 n, alpha, Amid, n, Bmid, o, beta, C2, o); 01431 01432 delete[] Amid; 01433 delete[] Bmid; 01434 01435 double* Cmid = new double[m*o]; 01436 double* Crad = new double[m*o]; 01437 01438 for(int i=0 ; i<m ; i++) { 01439 for(int j=0 ; j<o ; j++) { 01440 ind = i*o+j; 01441 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]); 01442 Crad[ind] = Cmid[ind] - C1[ind]; 01443 } 01444 } 01445 01446 for(int i=0 ; i<n ; i++) { 01447 for(int j=0 ; j<o ; j++) { 01448 ind = i*o+j; 01449 Babs[ind] += Brad[ind]; 01450 } 01451 } 01452 01453 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01454 n, alpha, Arad, n, Babs, o, beta, C1, o); 01455 01456 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01457 n, alpha, Aabs, n, Brad, o, beta, C2, o); 01458 01459 delete[] Aabs; 01460 delete[] Arad; 01461 delete[] Babs; 01462 delete[] Brad; 01463 01464 Resize(C, lbC_1, lbC_1+m-1, lbC_2, lbC_2+o-1); 01465 //C = imatrix(lbC_1, lbC_1+m-1, lbC_2, lbC_2+o-1); 01466 01467 for(int i=0 ; i<m ; i++) { 01468 for(int j=0 ; j<o ; j++) { 01469 ind = i*o+j; 01470 Crad[ind] += C1[ind] + C2[ind]; 01471 UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]); 01472 } 01473 } 01474 01475 setround(-1); 01476 01477 for(int i=0 ; i<m ; i++) { 01478 for(int j=0 ; j<o ; j++) { 01479 ind = i*o+j; 01480 UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]); 01481 } 01482 } 01483 01484 setround(rnd); 01485 01486 delete[] C1; 01487 delete[] C2; 01488 delete[] Crad; 01489 delete[] Cmid; 01490 01491 } 01492 01493 01494 inline void blasmatmul(const rmatrix &A, const imatrix &B, imatrix &C) { 01495 01496 const int m = Ub(A,1) - Lb(A,1) + 1; 01497 const int n = Ub(A,2) - Lb(A,2) + 1; 01498 const int o = Ub(B,2) - Lb(B,2) + 1; 01499 01500 int lb1_C = Lb(C,1); 01501 int lb2_C = Lb(C,2); 01502 Resize(C); 01503 01504 double* DA = (double*)A.to_blas_array(); 01505 double* Bmid = new double[n*o]; 01506 double* Brad = new double[n*o]; 01507 double* C1 = new double[m*o]; 01508 double* C2 = new double[m*o]; 01509 int rnd = getround(); 01510 01511 const double alpha = 1.0; 01512 double beta = 0.0; 01513 01514 setround(1); 01515 01516 int ind; 01517 01518 //Compute mid and rad of B 01519 for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) { 01520 for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) { 01521 ind = (i-Lb(B,1))*o+(j-Lb(B,2)); 01522 Bmid[ind] = _double(Inf(B[i][j]) + 0.5*(Sup(B[i][j]) - Inf(B[i][j]))); 01523 Brad[ind] = _double(Bmid[ind] - Inf(B[i][j])); 01524 } 01525 } 01526 01527 setround(-1); 01528 01529 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01530 n, alpha, DA, n, Bmid, o, beta, C1, o); 01531 01532 setround(1); 01533 01534 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01535 n, alpha, DA, n, Bmid, o, beta, C2, o); 01536 01537 delete[] Bmid; 01538 01539 double* Cmid = new double[m*o]; 01540 for(int i=0 ; i<m ; i++) { 01541 for(int j=0 ; j<o ; j++) { 01542 ind = i*o+j; 01543 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]); 01544 } 01545 } 01546 01547 delete[] C2; 01548 01549 double* Crad = new double[m*o]; 01550 for(int i=0 ; i<m ; i++) { 01551 for(int j=0 ; j<o ; j++) { 01552 ind = i*o+j; 01553 Crad[ind] = Cmid[ind] - C1[ind]; 01554 } 01555 } 01556 01557 delete[] C1; 01558 01559 double* Aabs = new double[m*n]; 01560 //Compute abs(A) 01561 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) { 01562 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) { 01563 ind = (i-Lb(A,1))*n+(j-Lb(A,2)); 01564 Aabs[ind] = fabs(DA[ind]); 01565 } 01566 } 01567 01568 beta = 1.0; 01569 01570 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01571 n, alpha, Aabs, n, Brad, o, beta, Crad, o); 01572 01573 delete[] Aabs; 01574 delete[] Brad; 01575 01576 Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1); 01577 //C = imatrix(lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1); 01578 01579 for(int i=0 ; i<m ; i++) { 01580 for(int j=0 ; j<o ; j++) { 01581 ind = i*o+j; 01582 UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]); 01583 } 01584 } 01585 01586 setround(-1); 01587 01588 for(int i=0 ; i<m ; i++) { 01589 for(int j=0 ; j<o ; j++) { 01590 ind = i*o+j; 01591 UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]); 01592 } 01593 } 01594 01595 setround(rnd); 01596 01597 delete[] Crad; 01598 delete[] Cmid; 01599 } 01600 01601 01602 inline void blasmatmul(const imatrix &A, const rmatrix &B, imatrix &C) { 01603 01604 const int m = Ub(A,1) - Lb(A,1) + 1; 01605 const int n = Ub(A,2) - Lb(A,2) + 1; 01606 const int o = Ub(C,2) - Lb(C,2) + 1; 01607 01608 int lb1_C = Lb(C,1); 01609 int lb2_C = Lb(C,2); 01610 Resize(C); 01611 01612 double* DB = (double*) B.to_blas_array(); 01613 double* Amid = new double[m*n]; 01614 double* Arad = new double[m*n]; 01615 double* C1 = new double[m*o]; 01616 double* C2 = new double[m*o]; 01617 int rnd = getround(); 01618 01619 const double alpha = 1.0; 01620 double beta = 0.0; 01621 01622 setround(1); 01623 01624 int ind; 01625 01626 //Compute mid and rad of A 01627 for(int i=Lb(A,1) ; i<=Ub(A,1) ; i++) { 01628 for(int j=Lb(A,2) ; j<=Ub(A,2) ; j++) { 01629 ind = (i-Lb(A,1))*n+(j-Lb(A,2)); 01630 Amid[ind] = _double(Inf(A[i][j]) + 0.5*(Sup(A[i][j]) - Inf(A[i][j]))); 01631 Arad[ind] = _double(Amid[ind] - Inf(A[i][j])); 01632 } 01633 } 01634 01635 setround(-1); 01636 01637 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01638 n, alpha, Amid, n, DB, o, beta, C1, o); 01639 01640 setround(1); 01641 01642 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01643 n, alpha, Amid, n, DB, o, beta, C2, o); 01644 01645 delete[] Amid; 01646 01647 double* Cmid = new double[m*o]; 01648 for(int i=0 ; i<m ; i++) { 01649 for(int j=0 ; j<o ; j++) { 01650 ind = i*o+j; 01651 Cmid[ind] = C1[ind] + 0.5*(C2[ind] - C1[ind]); 01652 } 01653 } 01654 01655 delete[] C2; 01656 01657 double* Crad = new double[m*o]; 01658 for(int i=0 ; i<m ; i++) { 01659 for(int j=0 ; j<o ; j++) { 01660 ind = i*o+j; 01661 Crad[ind] = Cmid[ind] - C1[ind]; 01662 } 01663 } 01664 01665 delete[] C1; 01666 01667 double* Babs = new double[n*o]; 01668 //Compute abs of B 01669 for(int i=Lb(B,1) ; i<=Ub(B,1) ; i++) { 01670 for(int j=Lb(B,2) ; j<=Ub(B,2) ; j++) { 01671 ind = (i-Lb(B,1))*o+(j-Lb(B,2)); 01672 Babs[ind] = fabs(DB[ind]); 01673 } 01674 } 01675 01676 beta = 1.0; 01677 01678 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01679 n, alpha, Arad, n, Babs, o, beta, Crad, o); 01680 01681 delete[] Arad; 01682 delete[] Babs; 01683 01684 Resize(C, lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1); 01685 //C = imatrix(lb1_C, lb1_C+m-1, lb2_C, lb2_C+o-1); 01686 01687 for(int i=0 ; i<m ; i++) { 01688 for(int j=0 ; j<o ; j++) { 01689 ind = i*o+j; 01690 UncheckedSetSup(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] + Crad[ind]); 01691 } 01692 } 01693 01694 setround(-1); 01695 01696 for(int i=0 ; i<m ; i++) { 01697 for(int j=0 ; j<o ; j++) { 01698 ind = i*o+j; 01699 UncheckedSetInf(C[i+Lb(C,1)][j+Lb(C,2)], Cmid[ind] - Crad[ind]); 01700 } 01701 } 01702 01703 setround(rnd); 01704 01705 delete[] Crad; 01706 delete[] Cmid; 01707 } 01708 #endif 01709 #endif 01710 01711 #ifdef _CXSC_BLAS_CMATRIX 01712 #ifndef _CXSC_BLAS_CMATRIX_INC 01713 #define _CXSC_BLAS_CMATRIX_INC 01714 inline void blasmatmul(const cmatrix &A, const cmatrix &B, cmatrix &C) { 01715 01716 double* DA = (double*) A.to_blas_array(); 01717 double* DB = (double*) B.to_blas_array(); 01718 double* DC = (double*) C.to_blas_array(); 01719 01720 const int m = Ub(A,1) - Lb(A,1) + 1; 01721 const int n = Ub(A,2) - Lb(A,2) + 1; 01722 const int o = Ub(C,2) - Lb(C,2) + 1; 01723 const complex alpha(1.0,0.0); 01724 const complex beta(0.0,0.0); 01725 01726 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01727 n, &alpha, DA, n, DB, o, &beta, DC, o); 01728 } 01729 01730 01731 inline void blasmatmul(const cmatrix &A, const rmatrix &B, cmatrix &C) { 01732 cmatrix tmp(B); 01733 01734 double* DA = (double*) A.to_blas_array(); 01735 double* DB = (double*) tmp.to_blas_array(); 01736 double* DC = (double*) C.to_blas_array(); 01737 01738 const int m = Ub(A,1) - Lb(A,1) + 1; 01739 const int n = Ub(A,2) - Lb(A,2) + 1; 01740 const int o = Ub(C,2) - Lb(C,2) + 1; 01741 const complex alpha(1.0,0.0); 01742 const complex beta(0.0,0.0); 01743 01744 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01745 n, &alpha, DA, n, DB, o, &beta, DC, o); 01746 } 01747 01748 inline void blasmatmul(const rmatrix &A, const cmatrix &B, cmatrix &C) { 01749 cmatrix tmp(A); 01750 01751 double* DA = (double*) tmp.to_blas_array(); 01752 double* DB = (double*) B.to_blas_array(); 01753 double* DC = (double*) C.to_blas_array(); 01754 01755 const int m = Ub(A,1) - Lb(A,1) + 1; 01756 const int n = Ub(A,2) - Lb(A,2) + 1; 01757 const int o = Ub(C,2) - Lb(C,2) + 1; 01758 const complex alpha(1.0,0.0); 01759 const complex beta(0.0,0.0); 01760 01761 cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, o, 01762 n, &alpha, DA, n, DB, o, &beta, DC, o); 01763 } 01764 #endif 01765 #endif 01766 01767 #ifdef _CXSC_BLAS_CIMATRIX 01768 #ifndef _CXSC_BLAS_CIMATRIX_INC 01769 #define _CXSC_BLAS_CIMATRIX_INC 01770 inline void blasmatmul(const cmatrix &A, const cmatrix &B, cimatrix &C) { 01771 const int m = Ub(A,1) - Lb(A,1) + 1; 01772 const int o = Ub(C,2) - Lb(C,2) + 1; 01773 imatrix tmp1(m,o), tmp2(m,o); 01774 01775 blasmatmul(Re(A),Re(B),tmp1); 01776 blasmatmul(Im(A),Im(B),tmp2); 01777 SetRe(C,tmp1-tmp2); 01778 blasmatmul(Re(A),Im(B),tmp1); 01779 blasmatmul(Im(A),Re(B),tmp2); 01780 SetIm(C,tmp1+tmp2); 01781 } 01782 01783 01784 inline void blasmatmul(const cmatrix &A, const imatrix &B, cimatrix &C) { 01785 const int m = Ub(A,1) - Lb(A,1) + 1; 01786 const int o = Ub(C,2) - Lb(C,2) + 1; 01787 01788 imatrix T(m,o); 01789 01790 blasmatmul(Re(A),B,T); 01791 01792 SetRe(C, T); 01793 01794 blasmatmul(Im(A),B,T); 01795 01796 SetIm(C, T); 01797 } 01798 01799 inline void blasmatmul(const imatrix &A, const cmatrix &B, cimatrix &C) { 01800 const int m = Ub(A,1) - Lb(A,1) + 1; 01801 const int o = Ub(C,2) - Lb(C,2) + 1; 01802 01803 imatrix T(m,o); 01804 01805 blasmatmul(A,Re(B),T); 01806 01807 SetRe(C, T); 01808 01809 blasmatmul(A,Im(B),T); 01810 01811 SetIm(C, T); 01812 } 01813 01814 inline void blasmatmul(const rmatrix &A, const cimatrix &B, cimatrix &C) { 01815 const int m = Ub(A,1) - Lb(A,1) + 1; 01816 const int o = Ub(C,2) - Lb(C,2) + 1; 01817 01818 imatrix T(m,o); 01819 01820 blasmatmul(A,Re(B),T); 01821 01822 SetRe(C, T); 01823 01824 blasmatmul(A,Im(B),T); 01825 01826 SetIm(C, T); 01827 } 01828 01829 inline void blasmatmul(const cimatrix &A, const rmatrix &B, cimatrix &C) { 01830 const int m = Ub(A,1) - Lb(A,1) + 1; 01831 const int o = Ub(C,2) - Lb(C,2) + 1; 01832 01833 imatrix T(m,o); 01834 01835 blasmatmul(Re(A),B,T); 01836 01837 SetRe(C, T); 01838 01839 blasmatmul(Im(A),B,T); 01840 01841 SetIm(C, T); 01842 } 01843 01844 inline void blasmatmul(const imatrix &A, const cimatrix &B, cimatrix &C) { 01845 const int m = Ub(A,1) - Lb(A,1) + 1; 01846 const int o = Ub(C,2) - Lb(C,2) + 1; 01847 01848 imatrix T(m,o); 01849 01850 blasmatmul(A,Re(B),T); 01851 01852 SetRe(C, T); 01853 01854 blasmatmul(A,Im(B),T); 01855 01856 SetIm(C, T); 01857 } 01858 01859 inline void blasmatmul(const cimatrix &A, const imatrix &B, cimatrix &C) { 01860 const int m = Ub(A,1) - Lb(A,1) + 1; 01861 const int o = Ub(C,2) - Lb(C,2) + 1; 01862 01863 imatrix T(m,o); 01864 01865 blasmatmul(Re(A),B,T); 01866 01867 SetRe(C, T); 01868 01869 blasmatmul(Im(A),B,T); 01870 01871 SetIm(C, T); 01872 } 01873 01874 01875 inline void blasmatmul(const cimatrix &A, const cmatrix &B, cimatrix &C) { 01876 01877 const int m = Ub(A,1) - Lb(A,1) + 1; 01878 const int o = Ub(C,2) - Lb(C,2) + 1; 01879 int rnd = getround(); 01880 01881 imatrix T(m,o); 01882 01883 blasmatmul(Re(A),Re(B),T); 01884 01885 SetRe(C, T); 01886 01887 blasmatmul(-Im(A),Im(B),T); 01888 01889 setround(-1); 01890 01891 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01892 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01893 InfRe(C[i][j]) += Inf(T[i][j]); 01894 } 01895 } 01896 01897 setround(1); 01898 01899 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01900 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01901 SupRe(C[i][j]) += Sup(T[i][j]); 01902 } 01903 } 01904 01905 setround(0); 01906 01907 blasmatmul(Re(A),Im(B),T); 01908 01909 SetIm(C, T); 01910 01911 blasmatmul(Im(A),Re(B),T); 01912 01913 setround(-1); 01914 01915 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01916 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01917 InfIm(C[i][j]) += Inf(T[i][j]); 01918 } 01919 } 01920 01921 setround(1); 01922 01923 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01924 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01925 SupIm(C[i][j]) += Sup(T[i][j]); 01926 } 01927 } 01928 01929 setround(rnd); 01930 } 01931 01932 inline void blasmatmul(const cmatrix &A, const cimatrix &B, cimatrix &C) { 01933 int rnd = getround(); 01934 01935 imatrix T(Lb(C,1),Ub(C,1),Lb(C,2),Ub(C,2)); 01936 01937 blasmatmul(Re(A),Re(B),T); 01938 01939 SetRe(C, T); 01940 01941 blasmatmul(-Im(A),Im(B),T); 01942 01943 setround(-1); 01944 01945 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01946 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01947 InfRe(C[i][j]) += Inf(T[i][j]); 01948 } 01949 } 01950 01951 setround(1); 01952 01953 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01954 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01955 SupRe(C[i][j]) += Sup(T[i][j]); 01956 } 01957 } 01958 01959 setround(0); 01960 01961 blasmatmul(Re(A),Im(B),T); 01962 01963 SetIm(C, T); 01964 01965 blasmatmul(Im(A),Re(B),T); 01966 01967 setround(-1); 01968 01969 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01970 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01971 InfIm(C[i][j]) += Inf(T[i][j]); 01972 } 01973 } 01974 01975 setround(1); 01976 01977 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 01978 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 01979 SupIm(C[i][j]) += Sup(T[i][j]); 01980 } 01981 } 01982 01983 setround(rnd); 01984 } 01985 01986 01987 inline void blasmatmul(const cimatrix &A, const cimatrix &B, cimatrix &C) { 01988 const int m = Ub(A,1) - Lb(A,1) + 1; 01989 const int o = Ub(C,2) - Lb(C,2) + 1; 01990 int rnd = getround(); 01991 01992 imatrix T(m,o); 01993 01994 blasmatmul(Re(A),Re(B),T); 01995 01996 SetRe(C, T); 01997 01998 blasmatmul(-Im(A),Im(B),T); 01999 02000 setround(-1); 02001 02002 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 02003 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 02004 InfRe(C[i][j]) += Inf(T[i][j]); 02005 } 02006 } 02007 02008 setround(1); 02009 02010 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 02011 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 02012 SupRe(C[i][j]) += Sup(T[i][j]); 02013 } 02014 } 02015 02016 setround(0); 02017 02018 blasmatmul(Re(A),Im(B),T); 02019 02020 SetIm(C, T); 02021 02022 blasmatmul(Im(A),Re(B),T); 02023 02024 setround(-1); 02025 02026 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 02027 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 02028 InfIm(C[i][j]) += Inf(T[i][j]); 02029 } 02030 } 02031 02032 setround(1); 02033 02034 for(int i=Lb(C,1) ; i<=Ub(C,1) ; i++) { 02035 for(int j=Lb(C,2) ; j<=Ub(C,2) ; j++) { 02036 SupIm(C[i][j]) += Sup(T[i][j]); 02037 } 02038 } 02039 02040 setround(rnd); 02041 } 02042 02043 inline void blasmatmul(const rmatrix &A, const cmatrix &B, cimatrix &C) { 02044 int m = Ub(A,1)-Lb(A,1)+1; 02045 int n = Ub(A,2)-Lb(A,2)+1; 02046 imatrix re(m,n),im(m,n); 02047 02048 blasmatmul(A,Re(B),re); 02049 blasmatmul(A,Im(B),im); 02050 02051 SetRe(C,re); 02052 SetIm(C,im); 02053 } 02054 02055 inline void blasmatmul(const cmatrix &A, const rmatrix &B, cimatrix &C) { 02056 int m = Ub(A,1)-Lb(A,1)+1; 02057 int n = Ub(A,2)-Lb(A,2)+1; 02058 imatrix re(m,n),im(m,n); 02059 02060 blasmatmul(Re(A),B,re); 02061 blasmatmul(Im(A),B,im); 02062 02063 SetRe(C,re); 02064 SetIm(C,im); 02065 } 02066 02067 //========================================================================= 02068 02069 /* 02070 //Computes B=A+B 02071 inline void blasadd(const rvector& A, rvector& B) { 02072 double* DA = A.to_blas_array(); 02073 double* DB = B.to_blas_array(); 02074 02075 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 1 ); 02076 } 02077 02078 //Computes B=A+B 02079 inline void blasadd(const rvector& A, ivector& B) { 02080 double* DA = A.to_blas_array(); 02081 double* DB = B.to_blas_array(); 02082 int rnd = getround(); 02083 02084 setround(-1); 02085 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 2 ); 02086 setround(1); 02087 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB+1, 2 ); 02088 setround(rnd); 02089 } 02090 02091 //Computes B=A+B 02092 inline void blasadd(const rvector& A, cvector& B) { 02093 double* DA = A.to_blas_array(); 02094 double* DB = B.to_blas_array(); 02095 02096 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 2 ); 02097 } 02098 02099 //Computes B=A+B 02100 inline void blasadd(const rvector& A, civector& B) { 02101 double* DA = A.to_blas_array(); 02102 double* DB = B.to_blas_array(); 02103 int rnd = getround(); 02104 02105 setround(-1); 02106 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB, 4 ); 02107 setround(1); 02108 cblas_daxpy(VecLen(A), 1.0, DA, 1, DB+1, 4 ); 02109 setround(rnd); 02110 } 02111 02112 02113 //Computes B=A+B 02114 inline void blasadd(const cvector& A, cvector& B) { 02115 double* DA = A.to_blas_array(); 02116 double* DB = B.to_blas_array(); 02117 02118 02119 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 2 ); 02120 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 2 ); 02121 } 02122 02123 //Computes B=A+B 02124 inline void blasadd(const cvector& A, civector& B) { 02125 double* DA = A.to_blas_array(); 02126 double* DB = B.to_blas_array(); 02127 int rnd = getround(); 02128 02129 setround(-1); 02130 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 4 ); 02131 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB+2, 4 ); 02132 setround(1); 02133 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 4 ); 02134 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+3, 4 ); 02135 setround(rnd); 02136 } 02137 02138 //Computes B=A+B 02139 inline void blasadd(const ivector& A, ivector& B) { 02140 double* DA = A.to_blas_array(); 02141 double* DB = B.to_blas_array(); 02142 int rnd = getround(); 02143 02144 setround(-1); 02145 02146 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 2 ); 02147 02148 setround(1); 02149 02150 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 2 ); 02151 02152 setround(rnd); 02153 } 02154 02155 //Computes B=A+B 02156 inline void blasadd(const ivector& A, civector& B) { 02157 double* DA = A.to_blas_array(); 02158 double* DB = B.to_blas_array(); 02159 int rnd = getround(); 02160 02161 setround(-1); 02162 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB, 4 ); 02163 cblas_daxpy(VecLen(A), 1.0, DA, 2, DB+2, 4 ); 02164 02165 setround(1); 02166 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+1, 4 ); 02167 cblas_daxpy(VecLen(A), 1.0, DA+1, 2, DB+3, 4 ); 02168 02169 setround(rnd); 02170 } 02171 02172 //Computes B=A+B 02173 inline void blasadd(const civector& A, civector& B) { 02174 double* DA = A.to_blas_array(); 02175 double* DB = B.to_blas_array(); 02176 int rnd = getround(); 02177 02178 setround(-1); 02179 02180 cblas_daxpy(VecLen(A), 1.0, DA, 4, DB, 4 ); 02181 cblas_daxpy(VecLen(A), 1.0, DA+2, 4, DB+2, 4 ); 02182 02183 setround(1); 02184 02185 cblas_daxpy(VecLen(A), 1.0, DA+1, 4, DB+1, 4 ); 02186 cblas_daxpy(VecLen(A), 1.0, DA+3, 4, DB+3, 4 ); 02187 02188 setround(rnd); 02189 } 02190 02191 //Computes C=A+B 02192 inline void blasadd(const cvector& A, const ivector& B, civector& C) { 02193 C = A; 02194 double* DB = B.to_blas_array(); 02195 double* DC = C.to_blas_array(); 02196 int rnd = getround(); 02197 02198 setround(-1); 02199 02200 cblas_daxpy(VecLen(A), 1.0, DB, 4, DC, 4 ); 02201 02202 setround(1); 02203 02204 cblas_daxpy(VecLen(A), 1.0, DB+1, 4, DC+1, 4 ); 02205 02206 setround(rnd); 02207 } 02208 02209 //Computes B=-A+B 02210 inline void blassub(const rvector& A, rvector& B) { 02211 double* DA = A.to_blas_array(); 02212 double* DB = B.to_blas_array(); 02213 02214 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 1 ); 02215 } 02216 02217 //Computes B=-A+B 02218 inline void blassub(const rvector& A, ivector& B) { 02219 double* DA = A.to_blas_array(); 02220 double* DB = B.to_blas_array(); 02221 int rnd = getround(); 02222 02223 setround(-1); 02224 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 2 ); 02225 setround(1); 02226 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB+1, 2 ); 02227 setround(rnd); 02228 } 02229 02230 //Computes B=-A+B 02231 inline void blassub(const rvector& A, cvector& B) { 02232 double* DA = A.to_blas_array(); 02233 double* DB = B.to_blas_array(); 02234 02235 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 2 ); 02236 } 02237 02238 //Computes B=-A+B 02239 inline void blassub(const rvector& A, civector& B) { 02240 double* DA = A.to_blas_array(); 02241 double* DB = B.to_blas_array(); 02242 int rnd = getround(); 02243 02244 setround(-1); 02245 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB, 4 ); 02246 setround(1); 02247 cblas_daxpy(VecLen(A), -1.0, DA, 1, DB+1, 4 ); 02248 setround(rnd); 02249 } 02250 02251 02252 //Computes B=-A+B 02253 inline void blassub(const cvector& A, cvector& B) { 02254 double* DA = A.to_blas_array(); 02255 double* DB = B.to_blas_array(); 02256 02257 02258 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB, 2 ); 02259 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+1, 2 ); 02260 } 02261 02262 //Computes B=-A+B 02263 inline void blassub(const cvector& A, civector& B) { 02264 double* DA = A.to_blas_array(); 02265 double* DB = B.to_blas_array(); 02266 int rnd = getround(); 02267 02268 setround(-1); 02269 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB, 4 ); 02270 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+2, 4 ); 02271 setround(1); 02272 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+1, 4 ); 02273 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+3, 4 ); 02274 setround(rnd); 02275 } 02276 02277 //Computes B=-A+B 02278 inline void blassub(const ivector& A, ivector& B) { 02279 double* DA = A.to_blas_array(); 02280 double* DB = B.to_blas_array(); 02281 int rnd = getround(); 02282 02283 setround(-1); 02284 02285 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB, 2 ); 02286 02287 setround(1); 02288 02289 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+1, 2 ); 02290 02291 setround(rnd); 02292 } 02293 02294 //Computes B=-A+B 02295 inline void blassub(const ivector& A, civector& B) { 02296 double* DA = A.to_blas_array(); 02297 double* DB = B.to_blas_array(); 02298 int rnd = getround(); 02299 02300 setround(-1); 02301 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB, 4 ); 02302 cblas_daxpy(VecLen(A), -1.0, DA+1, 2, DB+2, 4 ); 02303 02304 setround(1); 02305 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+1, 4 ); 02306 cblas_daxpy(VecLen(A), -1.0, DA, 2, DB+3, 4 ); 02307 02308 setround(rnd); 02309 } 02310 02311 //Computes B=-A+B 02312 inline void blassub(const civector& A, civector& B) { 02313 double* DA = A.to_blas_array(); 02314 double* DB = B.to_blas_array(); 02315 int rnd = getround(); 02316 02317 setround(-1); 02318 02319 cblas_daxpy(VecLen(A), -1.0, DA+1, 4, DB, 4 ); 02320 cblas_daxpy(VecLen(A), -1.0, DA+3, 4, DB+2, 4 ); 02321 02322 setround(1); 02323 02324 cblas_daxpy(VecLen(A), -1.0, DA, 4, DB+1, 4 ); 02325 cblas_daxpy(VecLen(A), -1.0, DA+2, 4, DB+3, 4 ); 02326 02327 setround(rnd); 02328 } 02329 02330 //============================================================== 02331 02332 //Computes B=A+B 02333 inline void blasadd(const rmatrix& A, rmatrix& B) { 02334 double* DA = A.to_blas_array(); 02335 double* DB = B.to_blas_array(); 02336 02337 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 1 ); 02338 } 02339 02340 //Computes B=A+B 02341 inline void blasadd(const rmatrix& A, imatrix& B) { 02342 double* DA = A.to_blas_array(); 02343 double* DB = B.to_blas_array(); 02344 int rnd = getround(); 02345 02346 setround(-1); 02347 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 2 ); 02348 setround(1); 02349 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB+1, 2 ); 02350 setround(rnd); 02351 } 02352 02353 //Computes B=A+B 02354 inline void blasadd(const rmatrix& A, cmatrix& B) { 02355 double* DA = A.to_blas_array(); 02356 double* DB = B.to_blas_array(); 02357 02358 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 2 ); 02359 } 02360 02361 //Computes B=A+B 02362 inline void blasadd(const rmatrix& A, cimatrix& B) { 02363 double* DA = A.to_blas_array(); 02364 double* DB = B.to_blas_array(); 02365 int rnd = getround(); 02366 02367 setround(-1); 02368 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB, 4 ); 02369 setround(1); 02370 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 1, DB+1, 4 ); 02371 setround(rnd); 02372 } 02373 02374 02375 //Computes B=A+B 02376 inline void blasadd(const cmatrix& A, cmatrix& B) { 02377 double* DA = A.to_blas_array(); 02378 double* DB = B.to_blas_array(); 02379 02380 02381 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 2 ); 02382 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 2 ); 02383 } 02384 02385 //Computes B=A+B 02386 inline void blasadd(const cmatrix& A, cimatrix& B) { 02387 double* DA = A.to_blas_array(); 02388 double* DB = B.to_blas_array(); 02389 int rnd = getround(); 02390 02391 setround(-1); 02392 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 4 ); 02393 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB+2, 4 ); 02394 setround(1); 02395 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 4 ); 02396 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+3, 4 ); 02397 setround(rnd); 02398 } 02399 02400 //Computes B=A+B 02401 inline void blasadd(const imatrix& A, imatrix& B) { 02402 double* DA = A.to_blas_array(); 02403 double* DB = B.to_blas_array(); 02404 int rnd = getround(); 02405 02406 setround(-1); 02407 02408 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 2 ); 02409 02410 setround(1); 02411 02412 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 2 ); 02413 02414 setround(rnd); 02415 } 02416 02417 //Computes B=A+B 02418 inline void blasadd(const imatrix& A, cimatrix& B) { 02419 double* DA = A.to_blas_array(); 02420 double* DB = B.to_blas_array(); 02421 int rnd = getround(); 02422 02423 setround(-1); 02424 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB, 4 ); 02425 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 2, DB+2, 4 ); 02426 02427 setround(1); 02428 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+1, 4 ); 02429 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 2, DB+3, 4 ); 02430 02431 setround(rnd); 02432 } 02433 02434 //Computes B=A+B 02435 inline void blasadd(const cimatrix& A, cimatrix& B) { 02436 double* DA = A.to_blas_array(); 02437 double* DB = B.to_blas_array(); 02438 int rnd = getround(); 02439 02440 setround(-1); 02441 02442 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA, 4, DB, 4 ); 02443 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+2, 4, DB+2, 4 ); 02444 02445 setround(1); 02446 02447 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+1, 4, DB+1, 4 ); 02448 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DA+3, 4, DB+3, 4 ); 02449 02450 setround(rnd); 02451 } 02452 02453 //Computes C=A+B 02454 inline void blasadd(const cmatrix& A, const imatrix& B, cimatrix& C) { 02455 C = A; 02456 double* DB = B.to_blas_array(); 02457 double* DC = C.to_blas_array(); 02458 int rnd = getround(); 02459 02460 setround(-1); 02461 02462 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DB, 4, DC, 4 ); 02463 02464 setround(1); 02465 02466 cblas_daxpy(RowLen(A)*ColLen(A), 1.0, DB+1, 4, DC+1, 4 ); 02467 02468 setround(rnd); 02469 } 02470 02471 //Computes B=-A+B 02472 inline void blassub(const rmatrix& A, rmatrix& B) { 02473 double* DA = A.to_blas_array(); 02474 double* DB = B.to_blas_array(); 02475 02476 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 1 ); 02477 } 02478 02479 //Computes B=-A+B 02480 inline void blassub(const rmatrix& A, imatrix& B) { 02481 double* DA = A.to_blas_array(); 02482 double* DB = B.to_blas_array(); 02483 int rnd = getround(); 02484 02485 setround(-1); 02486 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 2 ); 02487 setround(1); 02488 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB+1, 2 ); 02489 setround(rnd); 02490 } 02491 02492 //Computes B=-A+B 02493 inline void blassub(const rmatrix& A, cmatrix& B) { 02494 double* DA = A.to_blas_array(); 02495 double* DB = B.to_blas_array(); 02496 02497 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 2 ); 02498 } 02499 02500 //Computes B=-A+B 02501 inline void blassub(const rmatrix& A, cimatrix& B) { 02502 double* DA = A.to_blas_array(); 02503 double* DB = B.to_blas_array(); 02504 int rnd = getround(); 02505 02506 setround(-1); 02507 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB, 4 ); 02508 setround(1); 02509 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 1, DB+1, 4 ); 02510 setround(rnd); 02511 } 02512 02513 02514 //Computes B=-A+B 02515 inline void blassub(const cmatrix& A, cmatrix& B) { 02516 double* DA = A.to_blas_array(); 02517 double* DB = B.to_blas_array(); 02518 02519 02520 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB, 2 ); 02521 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+1, 2 ); 02522 } 02523 02524 //Computes B=-A+B 02525 inline void blassub(const cmatrix& A, cimatrix& B) { 02526 double* DA = A.to_blas_array(); 02527 double* DB = B.to_blas_array(); 02528 int rnd = getround(); 02529 02530 setround(-1); 02531 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB, 4 ); 02532 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+2, 4 ); 02533 setround(1); 02534 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+1, 4 ); 02535 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+3, 4 ); 02536 setround(rnd); 02537 } 02538 02539 //Computes B=-A+B 02540 inline void blassub(const imatrix& A, imatrix& B) { 02541 double* DA = A.to_blas_array(); 02542 double* DB = B.to_blas_array(); 02543 int rnd = getround(); 02544 02545 setround(-1); 02546 02547 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB, 2 ); 02548 02549 setround(1); 02550 02551 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+1, 2 ); 02552 02553 setround(rnd); 02554 } 02555 02556 //Computes B=-A+B 02557 inline void blassub(const imatrix& A, cimatrix& B) { 02558 double* DA = A.to_blas_array(); 02559 double* DB = B.to_blas_array(); 02560 int rnd = getround(); 02561 02562 setround(-1); 02563 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB, 4 ); 02564 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 2, DB+2, 4 ); 02565 02566 setround(1); 02567 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+1, 4 ); 02568 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 2, DB+3, 4 ); 02569 02570 setround(rnd); 02571 } 02572 02573 //Computes B=-A+B 02574 inline void blassub(const cimatrix& A, cimatrix& B) { 02575 double* DA = A.to_blas_array(); 02576 double* DB = B.to_blas_array(); 02577 int rnd = getround(); 02578 02579 setround(-1); 02580 02581 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+1, 4, DB, 4 ); 02582 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+3, 4, DB+2, 4 ); 02583 02584 setround(1); 02585 02586 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA, 4, DB+1, 4 ); 02587 cblas_daxpy(RowLen(A)*ColLen(A), -1.0, DA+2, 4, DB+3, 4 ); 02588 02589 setround(rnd); 02590 } */ 02591 02592 #endif 02593 #endif 02594 02595 }