apidoc/html/cxsc__blas_8inl_source.html

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 }