1 /////////////////////////////////////////////////////////////////////////////// 2 // // 3 // The Template Matrix/Vector Library for C++ was created by Mike Jarvis // 4 // Copyright (C) 1998 - 2016 // 5 // All rights reserved // 6 // // 7 // The project is hosted at https://code.google.com/p/tmv-cpp/ // 8 // where you can find the current version and current documention. // 9 // // 10 // For concerns or problems with the software, Mike may be contacted at // 11 // mike_jarvis17 [at] gmail. // 12 // // 13 // This software is licensed under a FreeBSD license. The file // 14 // TMV_LICENSE should have bee included with this distribution. // 15 // It not, you can get a copy from https://code.google.com/p/tmv-cpp/. // 16 // // 17 // Essentially, you can use this software however you want provided that // 18 // you include the TMV_LICENSE file in any distribution that uses it. // 19 // // 20 /////////////////////////////////////////////////////////////////////////////// 21 22 23 //#define XDEBUG 24 25 26 #include "TMV_Blas.h" 27 #include "tmv/TMV_SymBandMatrixArithFunc.h" 28 #include "tmv/TMV_SymBandMatrix.h" 29 #include "tmv/TMV_Vector.h" 30 #include "tmv/TMV_VectorArith.h" 31 #include "tmv/TMV_BandMatrix.h" 32 #include "tmv/TMV_BandMatrixArith.h" 33 #ifdef BLAS 34 #include "tmv/TMV_SymBandMatrixArith.h" 35 #endif 36 #include <iostream> 37 38 #ifdef XDEBUG 39 #include "tmv/TMV_MatrixArith.h" 40 #include <iostream> 41 using std::cout; 42 using std::cerr; 43 using std::endl; 44 #endif 45 46 namespace tmv { 47 48 template <class T> cptr() const49 const T* SymBandMatrixComposite<T>::cptr() const 50 { 51 if (!itsm1.get()) { 52 ptrdiff_t s = this->size(); 53 ptrdiff_t lo = this->nlo(); 54 ptrdiff_t len = BandStorageLength(ColMajor,s,s,lo,0); 55 itsm1.resize(len); 56 itsm = itsm1.get(); 57 this->assignTosB(SymBandMatrixView<T>( 58 itsm,s,lo,stepi(),stepj(),diagstep(), 59 Sym,this->uplo(),NonConj 60 TMV_FIRSTLAST1(itsm1.get(),itsm1.get()+len) )); 61 } 62 return itsm; 63 } 64 65 template <class T> stepi() const66 ptrdiff_t SymBandMatrixComposite<T>::stepi() const 67 { return 1; } 68 69 template <class T> stepj() const70 ptrdiff_t SymBandMatrixComposite<T>::stepj() const 71 { return this->nlo(); } 72 73 template <class T> diagstep() const74 ptrdiff_t SymBandMatrixComposite<T>::diagstep() const 75 { return this->nlo() + 1; } 76 77 78 // 79 // MultMV 80 // 81 82 template <bool add, class T, class Ta, class Tx> UnitAMultMV(const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<T> y)83 static void UnitAMultMV( 84 const GenSymBandMatrix<Ta>& A, const GenVector<Tx>& x, 85 VectorView<T> y) 86 { 87 if (add) y += A.lowerBand() * x; 88 else y = A.lowerBand() * x; 89 90 const ptrdiff_t N = A.size(); 91 if (N > 1 && A.nlo() > 0) 92 y.subVector(0,N-1) += A.upperBandOff() * x.subVector(1,N); 93 } 94 95 template <bool add, class T, class Ta, class Tx> NonBlasMultMV(const T alpha,const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<T> y)96 static void NonBlasMultMV( 97 const T alpha, const GenSymBandMatrix<Ta>& A, const GenVector<Tx>& x, 98 VectorView<T> y) 99 // y (+)= alpha * A * x 100 { 101 TMVAssert(A.size() == x.size()); 102 TMVAssert(A.size() == y.size()); 103 TMVAssert(alpha != T(0)); 104 TMVAssert(y.size() > 0); 105 TMVAssert(!SameStorage(x,y)); 106 107 if (A.uplo() == Upper) 108 if (A.isherm()) NonBlasMultMV<add>(alpha,A.adjoint(),x,y); 109 else NonBlasMultMV<add>(alpha,A.transpose(),x,y); 110 else if (y.isconj()) 111 NonBlasMultMV<add>( 112 TMV_CONJ(alpha),A.conjugate(),x.conjugate(),y.conjugate()); 113 else { 114 if (x.step() != 1) { 115 if (TMV_IMAG(alpha) == TMV_RealType(T)(0)) { 116 Vector<Tx> xx = TMV_REAL(alpha)*x; 117 if (y.step() != 1) { 118 Vector<T> yy(y.size()); 119 UnitAMultMV<false>(A,xx,yy.view()); 120 if (!add) y = yy; 121 else y += yy; 122 } else { 123 UnitAMultMV<add>(A,xx,y); 124 } 125 } else { 126 Vector<T> xx = alpha*x; 127 if (y.step()!=1) { 128 Vector<T> yy(y.size()); 129 UnitAMultMV<false>(A,xx,yy.view()); 130 if (add) y += yy; 131 else y = yy; 132 } else { 133 UnitAMultMV<add>(A,xx,y); 134 } 135 } 136 } else if (y.step()!=1 || alpha!=T(1)) { 137 Vector<T> yy(y.size()); 138 UnitAMultMV<false>(A,x,yy.view()); 139 if (add) y += alpha * yy; 140 else y = alpha * yy; 141 } else { 142 TMVAssert(alpha == T(1)); 143 TMVAssert(y.step() == 1); 144 TMVAssert(x.step() == 1); 145 UnitAMultMV<add>(A,x,y); 146 } 147 } 148 } 149 150 #ifdef BLAS 151 template <class T, class Ta, class Tx> BlasMultMV(const T alpha,const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,const int beta,VectorView<T> y)152 static inline void BlasMultMV( 153 const T alpha, const GenSymBandMatrix<Ta>& A, 154 const GenVector<Tx>& x, const int beta, VectorView<T> y) 155 { 156 if (beta == 1) NonBlasMultMV<true>(alpha,A,x,y); 157 else NonBlasMultMV<false>(alpha,A,x,y); 158 } 159 #ifdef INST_DOUBLE 160 template <> BlasMultMV(const double alpha,const GenSymBandMatrix<double> & A,const GenVector<double> & x,const int beta,VectorView<double> y)161 void BlasMultMV( 162 const double alpha, 163 const GenSymBandMatrix<double>& A, const GenVector<double>& x, 164 const int beta, VectorView<double> y) 165 { 166 int n = A.size(); 167 int k = A.nlo(); 168 int lda = A.diagstep(); 169 int xs = x.step(); 170 int ys = y.step(); 171 const double* xp = x.cptr(); 172 if (xs < 0) xp += (n-1)*xs; 173 double* yp = y.ptr(); 174 if (ys < 0) yp += (n-1)*ys; 175 if (beta == 0) y.setZero(); 176 double xbeta(1); 177 const double* Aptr = A.cptr(); 178 if (A.uplo() == Upper) Aptr -= A.nlo(); 179 BLASNAME(dsbmv) ( 180 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 181 BLASV(n),BLASV(k),BLASV(alpha),BLASP(Aptr),BLASV(lda), 182 BLASP(xp),BLASV(xs),BLASV(xbeta), 183 BLASP(yp),BLASV(ys) BLAS1); 184 } 185 template <> BlasMultMV(const std::complex<double> alpha,const GenSymBandMatrix<std::complex<double>> & A,const GenVector<std::complex<double>> & x,const int beta,VectorView<std::complex<double>> y)186 void BlasMultMV( 187 const std::complex<double> alpha, 188 const GenSymBandMatrix<std::complex<double> >& A, 189 const GenVector<std::complex<double> >& x, 190 const int beta, VectorView<std::complex<double> > y) 191 { 192 if (A.isherm()) { 193 int n = A.size(); 194 int k = A.nlo(); 195 int lda = A.diagstep(); 196 int xs = x.step(); 197 int ys = y.step(); 198 const std::complex<double>* xp = x.cptr(); 199 if (xs < 0) xp += (n-1)*xs; 200 std::complex<double>* yp = y.ptr(); 201 if (ys < 0) yp += (n-1)*ys; 202 if (beta == 0) y.setZero(); 203 std::complex<double> xbeta(1); 204 const std::complex<double>* Aptr = A.cptr(); 205 if (A.uplo() == Upper) Aptr -= A.nlo(); 206 BLASNAME(zhbmv) ( 207 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 208 BLASV(n),BLASV(k),BLASP(&alpha),BLASP(Aptr),BLASV(lda), 209 BLASP(xp),BLASV(xs),BLASP(&xbeta), 210 BLASP(yp),BLASV(ys) BLAS1); 211 } else { 212 if (beta == 1) NonBlasMultMV<true>(alpha,A,x,y); 213 else NonBlasMultMV<false>(alpha,A,x,y); 214 } 215 } 216 template <> BlasMultMV(const std::complex<double> alpha,const GenSymBandMatrix<std::complex<double>> & A,const GenVector<double> & x,const int beta,VectorView<std::complex<double>> y)217 void BlasMultMV( 218 const std::complex<double> alpha, 219 const GenSymBandMatrix<std::complex<double> >& A, 220 const GenVector<double>& x, 221 const int beta, VectorView<std::complex<double> > y) 222 { BlasMultMV(alpha,A,Vector<std::complex<double> >(x),beta,y); } 223 template <> BlasMultMV(const std::complex<double> alpha,const GenSymBandMatrix<double> & A,const GenVector<std::complex<double>> & x,const int beta,VectorView<std::complex<double>> y)224 void BlasMultMV( 225 const std::complex<double> alpha, 226 const GenSymBandMatrix<double>& A, 227 const GenVector<std::complex<double> >& x, 228 const int beta, VectorView<std::complex<double> > y) 229 { 230 if (beta == 0) { 231 int n = A.size(); 232 int k = A.nlo(); 233 int lda = A.diagstep(); 234 int xs = 2*x.step(); 235 int ys = 2*y.step(); 236 const double* xp = (const double*) x.cptr(); 237 if (xs < 0) xp += (n-1)*xs; 238 double* yp = (double*) y.ptr(); 239 if (ys < 0) yp += (n-1)*ys; 240 double xalpha(1); 241 y.setZero(); 242 double xbeta(1); 243 const double* Aptr = A.cptr(); 244 if (A.uplo() == Upper) Aptr -= A.nlo(); 245 BLASNAME(dsbmv) ( 246 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 247 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 248 BLASP(xp),BLASV(xs),BLASV(xbeta), 249 BLASP(yp),BLASV(ys) BLAS1); 250 BLASNAME(dsbmv) ( 251 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 252 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 253 BLASP(xp+1),BLASV(xs),BLASV(xbeta), 254 BLASP(yp+1),BLASV(ys) BLAS1); 255 if (x.isconj()) y.conjugateSelf(); 256 y *= alpha; 257 } else if (TMV_IMAG(alpha) == 0. && !x.isconj()) { 258 int n = A.size(); 259 int k = A.nlo(); 260 int lda = A.diagstep(); 261 int xs = 2*x.step(); 262 int ys = 2*y.step(); 263 const double* xp = (const double*) x.cptr(); 264 if (xs < 0) xp += (n-1)*xs; 265 double* yp = (double*) y.ptr(); 266 if (ys < 0) yp += (n-1)*ys; 267 double xalpha(TMV_REAL(alpha)); 268 double xbeta(1); 269 const double* Aptr = A.cptr(); 270 if (A.uplo() == Upper) Aptr -= A.nlo(); 271 BLASNAME(dsbmv) ( 272 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 273 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 274 BLASP(xp),BLASV(xs),BLASV(xbeta), 275 BLASP(yp),BLASV(ys) BLAS1); 276 BLASNAME(dsbmv) ( 277 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 278 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 279 BLASP(xp+1),BLASV(xs),BLASV(xbeta), 280 BLASP(yp+1),BLASV(ys) BLAS1); 281 } else { 282 Vector<std::complex<double> > xx = alpha*x; 283 BlasMultMV(std::complex<double>(1),A,xx,1,y); 284 } 285 } 286 template <> BlasMultMV(const std::complex<double> alpha,const GenSymBandMatrix<double> & A,const GenVector<double> & x,const int beta,VectorView<std::complex<double>> y)287 void BlasMultMV( 288 const std::complex<double> alpha, 289 const GenSymBandMatrix<double>& A, 290 const GenVector<double>& x, 291 const int beta, VectorView<std::complex<double> > y) 292 { 293 int n = A.size(); 294 int k = A.nlo(); 295 int lda = A.diagstep(); 296 int xs = x.step(); 297 int ys = 2*y.step(); 298 const double* xp = x.cptr(); 299 if (xs < 0) xp += (n-1)*xs; 300 double* yp = (double*) y.ptr(); 301 if (ys < 0) yp += (n-1)*ys; 302 double ar(TMV_REAL(alpha)); 303 double ai(TMV_IMAG(alpha)); 304 if (beta == 0) y.setZero(); 305 double xbeta(1); 306 const double* Aptr = A.cptr(); 307 if (A.uplo() == Upper) Aptr -= A.nlo(); 308 if (ar != 0.) { 309 BLASNAME(dsbmv) ( 310 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 311 BLASV(n),BLASV(k),BLASV(ar),BLASP(Aptr),BLASV(lda), 312 BLASP(xp),BLASV(xs),BLASV(xbeta), 313 BLASP(yp),BLASV(ys) BLAS1); 314 } 315 if (ai != 0.) { 316 BLASNAME(dsbmv) ( 317 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 318 BLASV(n),BLASV(k),BLASV(ai),BLASP(Aptr),BLASV(lda), 319 BLASP(xp),BLASV(xs),BLASV(xbeta), 320 BLASP(yp+1),BLASV(ys) BLAS1); 321 } 322 } 323 #endif 324 #ifdef INST_FLOAT 325 template <> BlasMultMV(const float alpha,const GenSymBandMatrix<float> & A,const GenVector<float> & x,const int beta,VectorView<float> y)326 void BlasMultMV( 327 const float alpha, 328 const GenSymBandMatrix<float>& A, const GenVector<float>& x, 329 const int beta, VectorView<float> y) 330 { 331 int n = A.size(); 332 int k = A.nlo(); 333 int lda = A.diagstep(); 334 int xs = x.step(); 335 int ys = y.step(); 336 const float* xp = x.cptr(); 337 if (xs < 0) xp += (n-1)*xs; 338 float* yp = y.ptr(); 339 if (ys < 0) yp += (n-1)*ys; 340 if (beta == 0) y.setZero(); 341 float xbeta(1); 342 const float* Aptr = A.cptr(); 343 if (A.uplo() == Upper) Aptr -= A.nlo(); 344 BLASNAME(ssbmv) ( 345 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 346 BLASV(n),BLASV(k),BLASV(alpha),BLASP(Aptr),BLASV(lda), 347 BLASP(xp),BLASV(xs),BLASV(xbeta), 348 BLASP(yp),BLASV(ys) BLAS1); 349 } 350 template <> BlasMultMV(const std::complex<float> alpha,const GenSymBandMatrix<std::complex<float>> & A,const GenVector<std::complex<float>> & x,const int beta,VectorView<std::complex<float>> y)351 void BlasMultMV( 352 const std::complex<float> alpha, 353 const GenSymBandMatrix<std::complex<float> >& A, 354 const GenVector<std::complex<float> >& x, 355 const int beta, VectorView<std::complex<float> > y) 356 { 357 if (A.isherm()) { 358 int n = A.size(); 359 int k = A.nlo(); 360 int lda = A.diagstep(); 361 int xs = x.step(); 362 int ys = y.step(); 363 const std::complex<float>* xp = x.cptr(); 364 if (xs < 0) xp += (n-1)*xs; 365 std::complex<float>* yp = y.ptr(); 366 if (ys < 0) yp += (n-1)*ys; 367 if (beta == 0) y.setZero(); 368 std::complex<float> xbeta(1); 369 const std::complex<float>* Aptr = A.cptr(); 370 if (A.uplo() == Upper) Aptr -= A.nlo(); 371 BLASNAME(chbmv) ( 372 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 373 BLASV(n),BLASV(k),BLASP(&alpha),BLASP(Aptr),BLASV(lda), 374 BLASP(xp),BLASV(xs),BLASP(&xbeta), 375 BLASP(yp),BLASV(ys) BLAS1); 376 } else { 377 if (beta == 1) NonBlasMultMV<true>(alpha,A,x,y); 378 else NonBlasMultMV<false>(alpha,A,x,y); 379 } 380 } 381 template <> BlasMultMV(const std::complex<float> alpha,const GenSymBandMatrix<std::complex<float>> & A,const GenVector<float> & x,const int beta,VectorView<std::complex<float>> y)382 void BlasMultMV( 383 const std::complex<float> alpha, 384 const GenSymBandMatrix<std::complex<float> >& A, 385 const GenVector<float>& x, 386 const int beta, VectorView<std::complex<float> > y) 387 { BlasMultMV(alpha,A,Vector<std::complex<float> >(x),beta,y); } 388 template <> BlasMultMV(const std::complex<float> alpha,const GenSymBandMatrix<float> & A,const GenVector<std::complex<float>> & x,const int beta,VectorView<std::complex<float>> y)389 void BlasMultMV( 390 const std::complex<float> alpha, 391 const GenSymBandMatrix<float>& A, 392 const GenVector<std::complex<float> >& x, 393 const int beta, VectorView<std::complex<float> > y) 394 { 395 if (beta == 0) { 396 int n = A.size(); 397 int k = A.nlo(); 398 int lda = A.diagstep(); 399 int xs = 2*x.step(); 400 int ys = 2*y.step(); 401 const float* xp = (const float*) x.cptr(); 402 if (xs < 0) xp += (n-1)*xs; 403 float* yp = (float*) y.ptr(); 404 if (ys < 0) yp += (n-1)*ys; 405 float xalpha(1); 406 y.setZero(); 407 float xbeta(1); 408 const float* Aptr = A.cptr(); 409 if (A.uplo() == Upper) Aptr -= A.nlo(); 410 BLASNAME(ssbmv) ( 411 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 412 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 413 BLASP(xp),BLASV(xs),BLASV(xbeta), 414 BLASP(yp),BLASV(ys) BLAS1); 415 BLASNAME(ssbmv) ( 416 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 417 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 418 BLASP(xp+1),BLASV(xs),BLASV(xbeta), 419 BLASP(yp+1),BLASV(ys) BLAS1); 420 if (x.isconj()) y.conjugateSelf(); 421 y *= alpha; 422 } else if (TMV_IMAG(alpha) == 0.F && !x.isconj()) { 423 int n = A.size(); 424 int k = A.nlo(); 425 int lda = A.diagstep(); 426 int xs = 2*x.step(); 427 int ys = 2*y.step(); 428 const float* xp = (const float*) x.cptr(); 429 if (xs < 0) xp += (n-1)*xs; 430 float* yp = (float*) y.ptr(); 431 if (ys < 0) yp += (n-1)*ys; 432 float xalpha(TMV_REAL(alpha)); 433 float xbeta(1); 434 const float* Aptr = A.cptr(); 435 if (A.uplo() == Upper) Aptr -= A.nlo(); 436 BLASNAME(ssbmv) ( 437 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 438 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 439 BLASP(xp),BLASV(xs),BLASV(xbeta), 440 BLASP(yp),BLASV(ys) BLAS1); 441 BLASNAME(ssbmv) ( 442 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 443 BLASV(n),BLASV(k),BLASV(xalpha),BLASP(Aptr),BLASV(lda), 444 BLASP(xp+1),BLASV(xs),BLASV(xbeta), 445 BLASP(yp+1),BLASV(ys) BLAS1); 446 } else { 447 Vector<std::complex<float> > xx = alpha*x; 448 BlasMultMV(std::complex<float>(1),A,xx,1,y); 449 } 450 } 451 template <> BlasMultMV(const std::complex<float> alpha,const GenSymBandMatrix<float> & A,const GenVector<float> & x,const int beta,VectorView<std::complex<float>> y)452 void BlasMultMV( 453 const std::complex<float> alpha, 454 const GenSymBandMatrix<float>& A, 455 const GenVector<float>& x, 456 const int beta, VectorView<std::complex<float> > y) 457 { 458 int n = A.size(); 459 int k = A.nlo(); 460 int lda = A.diagstep(); 461 int xs = x.step(); 462 int ys = 2*y.step(); 463 const float* xp = x.cptr(); 464 if (xs < 0) xp += (n-1)*xs; 465 float* yp = (float*) y.ptr(); 466 if (ys < 0) yp += (n-1)*ys; 467 float ar(TMV_REAL(alpha)); 468 float ai(TMV_IMAG(alpha)); 469 if (beta == 0) y.setZero(); 470 float xbeta(1); 471 const float* Aptr = A.cptr(); 472 if (A.uplo() == Upper) Aptr -= A.nlo(); 473 if (ar != 0.F) { 474 BLASNAME(ssbmv) ( 475 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 476 BLASV(n),BLASV(k),BLASV(ar),BLASP(Aptr),BLASV(lda), 477 BLASP(xp),BLASV(xs),BLASV(xbeta), 478 BLASP(yp),BLASV(ys) BLAS1); 479 } 480 if (ai != 0.F) { 481 BLASNAME(ssbmv) ( 482 BLASCM A.uplo() == Upper?BLASCH_UP:BLASCH_LO, 483 BLASV(n),BLASV(k),BLASV(ai),BLASP(Aptr),BLASV(lda), 484 BLASP(xp),BLASV(xs),BLASV(xbeta), 485 BLASP(yp+1),BLASV(ys) BLAS1); 486 } 487 } 488 #endif 489 #endif // BLAS 490 491 template <bool add, class T, class Ta, class Tx> DoMultMV(const T alpha,const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<T> y)492 static void DoMultMV( 493 const T alpha, const GenSymBandMatrix<Ta>& A, 494 const GenVector<Tx>& x, VectorView<T> y) 495 { 496 TMVAssert(A.rowsize() == x.size()); 497 TMVAssert(A.colsize() == y.size()); 498 TMVAssert(alpha != T(0)); 499 TMVAssert(x.size() > 0); 500 TMVAssert(y.size() > 0); 501 TMVAssert(!SameStorage(x,y)); 502 503 #ifdef BLAS 504 if (!A.iscm() && A.isrm()) 505 if (A.isherm()) DoMultMV<add>(alpha,A.adjoint(),x,y); 506 else DoMultMV<add>(alpha,A.transpose(),x,y); 507 else if (A.isconj()) 508 DoMultMV<add>( 509 TMV_CONJ(alpha),A.conjugate(),x.conjugate(),y.conjugate()); 510 else if (x.step() == 0) { 511 if (x.size() <= 1) 512 DoMultMV<add>( 513 alpha,A,ConstVectorView<Tx>(x.cptr(),x.size(),1,x.ct()),y); 514 else 515 DoMultMV<add>(alpha,A,Vector<Tx>(x),y); 516 } else if (y.step() == 0) { 517 TMVAssert(y.size() <= 1); 518 DoMultMV<add>(alpha,A,x,VectorView<T>(y.ptr(),y.size(),1,y.ct())); 519 } else { 520 #ifdef XDEBUG 521 cout<<"BLAS section: alpha = "<<alpha<<endl; 522 cout<<"A = "<<TMV_Text(A)<<" "<<A<<endl; 523 cout<<"x = "<<TMV_Text(x)<<" step "<<x.step()<<" "<<x<<endl; 524 cout<<"y = "<<TMV_Text(y)<<" step "<<y.step()<<" "<<y<<endl; 525 #endif 526 if (A.iscm()&&(A.nlo()==0 || A.stepj()>0)) { 527 if (!y.isconj() && y.step() != 1) { 528 if (!x.isconj() && x.step() != 1) { 529 BlasMultMV(alpha,A,x,add?1:0,y); 530 } else { 531 Vector<T> xx = alpha*x; 532 BlasMultMV(T(1),A,xx,add?1:0,y); 533 } 534 } else { 535 Vector<T> yy(y.size()); 536 if (!x.isconj() && x.step() != 1) { 537 BlasMultMV(T(1),A,x,0,yy.view()); 538 if (add) y += alpha*yy; 539 else y = alpha*yy; 540 } else { 541 Vector<T> xx = alpha*x; 542 BlasMultMV(T(1),A,xx,0,yy.view()); 543 if (add) y += yy; 544 else y = yy; 545 } 546 } 547 } else { 548 if (TMV_IMAG(alpha) == TMV_RealType(T)(0)) { 549 if (A.isherm()) { 550 if (A.uplo() == Upper) { 551 HermBandMatrix<Ta,Upper|ColMajor> A2 = 552 TMV_REAL(alpha)*A; 553 DoMultMV<add>(T(1),A2,x,y); 554 } else { 555 HermBandMatrix<Ta,Lower|ColMajor> A2 = 556 TMV_REAL(alpha)*A; 557 DoMultMV<add>(T(1),A2,x,y); 558 } 559 } else { 560 if (A.uplo() == Upper) { 561 SymBandMatrix<Ta,Upper|ColMajor> A2 = 562 TMV_REAL(alpha)*A; 563 DoMultMV<add>(T(1),A2,x,y); 564 } else { 565 SymBandMatrix<Ta,Lower|ColMajor> A2 = 566 TMV_REAL(alpha)*A; 567 DoMultMV<add>(T(1),A2,x,y); 568 } 569 } 570 } else { 571 if (!A.issym()) { 572 if (A.uplo() == Upper) { 573 HermBandMatrix<Ta,Upper|ColMajor> A2 = A; 574 DoMultMV<add>(alpha,A2,x,y); 575 } else { 576 HermBandMatrix<Ta,Lower|ColMajor> A2 = A; 577 DoMultMV<add>(alpha,A2,x,y); 578 } 579 } else { 580 if (A.uplo() == Upper) { 581 SymBandMatrix<T,Upper|ColMajor> A2 = alpha*A; 582 DoMultMV<add>(T(1),A2,x,y); 583 } else { 584 SymBandMatrix<T,Lower|ColMajor> A2 = alpha*A; 585 DoMultMV<add>(T(1),A2,x,y); 586 } 587 } 588 } 589 } 590 } 591 #else 592 NonBlasMultMV<add>(alpha,A,x,y); 593 #endif 594 } 595 596 template <bool add, class T, class Ta, class Tx> MultMV(const T alpha,const GenSymBandMatrix<Ta> & A,const GenVector<Tx> & x,VectorView<T> y)597 void MultMV( 598 const T alpha, const GenSymBandMatrix<Ta>& A, const GenVector<Tx>& x, 599 VectorView<T> y) 600 // y (+)= alpha * A * x 601 { 602 TMVAssert(A.rowsize() == x.size()); 603 TMVAssert(A.colsize() == y.size()); 604 #ifdef XDEBUG 605 cout<<"Start MultMV: alpha = "<<alpha<<endl; 606 cout<<"A = "<<TMV_Text(A)<<" "<<A<<endl; 607 cout<<"x = "<<TMV_Text(x)<<" step "<<x.step()<<" "<<x<<endl; 608 cout<<"y = "<<TMV_Text(y)<<" step "<<y.step()<<" "<<y<<endl; 609 Matrix<Ta> A0 = A; 610 Vector<Tx> x0 = x; 611 Vector<T> y0 = y; 612 Vector<T> y2 = alpha*A0*x0; 613 if (add) y2 += y0; 614 #endif 615 616 if (y.size() > 0) { 617 if (x.size()==0 || alpha==T(0)) { 618 if (!add) y.setZero(); 619 } else if (SameStorage(x,y)) { 620 Vector<T> yy(y.size()); 621 DoMultMV<false>(T(1),A,x,yy.view()); 622 if (add) y += alpha*yy; 623 else y = alpha*yy; 624 } else { 625 DoMultMV<add>(alpha,A,x,y); 626 } 627 } 628 #ifdef XDEBUG 629 cout<<"--> y = "<<y<<endl; 630 cout<<"y2 = "<<y2<<endl; 631 if (!(Norm(y-y2) <= 632 0.001*(TMV_ABS(alpha)*Norm(A0)*Norm(x0)+ 633 (add?Norm(y0):TMV_RealType(T)(0))))) { 634 cerr<<"MultMV: alpha = "<<alpha<<endl; 635 cerr<<"A = "<<TMV_Text(A)<<" "<<A0<<endl; 636 cerr<<"x = "<<TMV_Text(x)<<" step "<<x.step()<<" "<<x0<<endl; 637 cerr<<"y = "<<TMV_Text(y)<<" step "<<y.step()<<" "<<y0<<endl; 638 cerr<<"--> y = "<<y<<endl; 639 cerr<<"y2 = "<<y2<<endl; 640 abort(); 641 } 642 #endif 643 } 644 645 #define InstFile "TMV_MultsBV.inst" 646 #include "TMV_Inst.h" 647 #undef InstFile 648 649 } // namespace tmv 650 651 652