C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
cxsc_blas.inl
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 }