C-XSC - A C++ Class Library for Extended Scientific Computing: cxsc

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