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