← Back to branch summary

~chaffra/+junk/trilinos

« back to all changes in this revision

Viewing changes to packages/teuchos/src/Teuchos_BLAS.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme, Johannes Ring
  • Date: 2009-12-13 12:53:22 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20091213125322-in0nrdjc55deqsw9
Tags: 10.0.3.dfsg-1
[Christophe Prud'homme]
* New upstream release

[Johannes Ring]
* debian/patches/libname.patch: Add prefix 'libtrilinos_' to all
  libraries. 
* debian/patches/soname.patch: Add soversion to libraries.
* debian/watch: Update download URL.
* debian/control:
  - Remove python-numeric from Build-Depends (virtual package).
  - Remove automake and autotools from Build-Depends and add cmake to
    reflect switch to CMake.
  - Add python-support to Build-Depends.
* debian/rules: 
  - Cleanup and updates for switch to CMake.

Show diffs side-by-side

added added

removed removed

Lines of Context:
packages/teuchos/src/Teuchos_BLAS.hpp
74
74
#include "Teuchos_ScalarTraits.hpp"
75
75
#include "Teuchos_OrdinalTraits.hpp"
76
76
#include "Teuchos_BLAS_types.hpp"
 
77
#include "Teuchos_TestForException.hpp"
77
78
 
78
79
/*! \class Teuchos::BLAS
79
80
    \brief The Templated BLAS Wrapper Class.
115
116
  extern const char EUploChar[];
116
117
  extern const char EDiagChar[];
117
118
 
 
119
  //! Default implementation for BLAS routines
 
120
  /*!
 
121
   * This class provides the default implementation for the BLAS routines.  It
 
122
   * is put in a separate class so that specializations of BLAS for other types
 
123
   * still have this implementation available.
 
124
   */
118
125
  template<typename OrdinalType, typename ScalarType>
119
 
  class BLAS
 
126
  class DefaultBLASImpl
120
127
  {    
121
128
 
122
129
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
126
133
    //@{ 
127
134
    
128
135
    //! Default constructor.
129
 
    inline BLAS(void) {}
 
136
    inline DefaultBLASImpl(void) {}
130
137
 
131
138
    //! Copy constructor.
132
 
    inline BLAS(const BLAS<OrdinalType, ScalarType>& /*BLAS_source*/) {}
 
139
    inline DefaultBLASImpl(const DefaultBLASImpl<OrdinalType, ScalarType>& /*BLAS_source*/) {}
133
140
 
134
141
    //! Destructor.
135
 
    inline virtual ~BLAS(void) {}
 
142
    inline virtual ~DefaultBLASImpl(void) {}
136
143
    //@}
137
144
 
138
145
    //! @name Level 1 BLAS Routines.
151
158
    void COPY(const OrdinalType n, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const;
152
159
 
153
160
    //! Perform the operation: \c y \c <- \c y+alpha*x.
154
 
    void AXPY(const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const;
 
161
    template <typename alpha_type, typename x_type>
 
162
    void AXPY(const OrdinalType n, const alpha_type alpha, const x_type* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const;
155
163
 
156
164
    //! Sum the absolute values of the entries of \c x.
157
165
    typename ScalarTraits<ScalarType>::magnitudeType ASUM(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const;
158
166
 
159
167
    //! Form the dot product of the vectors \c x and \c y.
160
 
    ScalarType DOT(const OrdinalType n, const ScalarType* x, const OrdinalType incx, const ScalarType* y, const OrdinalType incy) const;
 
168
    template <typename x_type, typename y_type>
 
169
    ScalarType DOT(const OrdinalType n, const x_type* x, const OrdinalType incx, const y_type* y, const OrdinalType incy) const;
161
170
 
162
171
    //! Compute the 2-norm of the std::vector \c x.
163
172
    typename ScalarTraits<ScalarType>::magnitudeType NRM2(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const;
171
180
    //@{ 
172
181
 
173
182
    //! Performs the matrix-std::vector operation:  \c y \c <- \c alpha*A*x+beta*y or \c y \c <- \c alpha*A'*x+beta*y where \c A is a general \c m by \c n matrix.
174
 
    void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, 
175
 
              const OrdinalType lda, const ScalarType* x, const OrdinalType incx, const ScalarType beta, ScalarType* y, const OrdinalType incy) const;
 
183
    template <typename alpha_type, typename A_type, typename x_type, typename beta_type>
 
184
    void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, 
 
185
              const OrdinalType lda, const x_type* x, const OrdinalType incx, const beta_type beta, ScalarType* y, const OrdinalType incy) const;
176
186
 
177
187
    //! Performs the matrix-std::vector operation:  \c x \c <- \c A*x or \c x \c <- \c A'*x where \c A is a unit/non-unit \c n by \c n upper/lower triangular matrix.
178
 
    void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const ScalarType* A, 
 
188
    template <typename A_type>
 
189
    void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const A_type* A, 
179
190
              const OrdinalType lda, ScalarType* x, const OrdinalType incx) const;
180
191
 
181
 
    //! Performs the rank 1 operation:  \c A \c <- \c alpha*x*y'+A.
182
 
    void GER(const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, 
183
 
             const ScalarType* y, const OrdinalType incy, ScalarType* A, const OrdinalType lda) const;
 
192
    //! \brief Performs the rank 1 operation:  \c A \c <- \c alpha*x*y'+A. 
 
193
    template <typename alpha_type, typename x_type, typename y_type>
 
194
    void GER(const OrdinalType m, const OrdinalType n, const alpha_type alpha, const x_type* x, const OrdinalType incx, 
 
195
             const y_type* y, const OrdinalType incy, ScalarType* A, const OrdinalType lda) const;
184
196
    //@}
185
197
    
186
198
    //! @name Level 3 BLAS Routines. 
187
199
    //@{ 
188
200
 
189
201
    //! Performs the matrix-matrix operation: \c C \c <- \c alpha*op(A)*op(B)+beta*C where \c op(A) is either \c A or \c A', \c op(B) is either \c B or \c B', and C is an \c m by \c k matrix.
190
 
    void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, const ScalarType* B, const OrdinalType ldb, const ScalarType beta, ScalarType* C, const OrdinalType ldc) const;
 
202
    template <typename alpha_type, typename A_type, typename B_type, typename beta_type>
 
203
    void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type* A, const OrdinalType lda, const B_type* B, const OrdinalType ldb, const beta_type beta, ScalarType* C, const OrdinalType ldc) const;
191
204
 
192
205
    //! Performs the matrix-matrix operation: \c C \c <- \c alpha*A*B+beta*C or \c C \c <- \c alpha*B*A+beta*C where \c A is an \c m by \c m or \c n by \c n symmetric matrix and \c B is a general matrix.
193
 
    void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, const ScalarType* B, const OrdinalType ldb, const ScalarType beta, ScalarType* C, const OrdinalType ldc) const;
 
206
    template <typename alpha_type, typename A_type, typename B_type, typename beta_type>
 
207
    void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, const OrdinalType lda, const B_type* B, const OrdinalType ldb, const beta_type beta, ScalarType* C, const OrdinalType ldc) const;
194
208
 
195
209
    //! Performs the matrix-matrix operation: \c C \c <- \c alpha*op(A)*B+beta*C or \c C \c <- \c alpha*B*op(A)+beta*C where \c op(A) is an unit/non-unit, upper/lower triangular matrix and \c B is a general matrix.
 
210
    template <typename alpha_type, typename A_type>
196
211
    void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n,
197
 
                const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
 
212
                const alpha_type alpha, const A_type* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
198
213
 
199
214
    //! Solves the matrix equations:  \c op(A)*X=alpha*B or \c X*op(A)=alpha*B where \c X and \c B are \c m by \c n matrices, \c A is a unit/non-unit, upper/lower triangular matrix and \c op(A) is \c A or \c A'.  The matrix \c X is overwritten on \c B.
 
215
    template <typename alpha_type, typename A_type>
200
216
    void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n,
201
 
                const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
 
217
                const alpha_type alpha, const A_type* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
 
218
    //@}
 
219
  };
 
220
 
 
221
  template<typename OrdinalType, typename ScalarType>
 
222
  class BLAS : public DefaultBLASImpl<OrdinalType,ScalarType>
 
223
  {    
 
224
 
 
225
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
 
226
    
 
227
  public:
 
228
    //! @name Constructor/Destructor.
 
229
    //@{ 
 
230
    
 
231
    //! Default constructor.
 
232
    inline BLAS(void) {}
 
233
 
 
234
    //! Copy constructor.
 
235
    inline BLAS(const BLAS<OrdinalType, ScalarType>& /*BLAS_source*/) {}
 
236
 
 
237
    //! Destructor.
 
238
    inline virtual ~BLAS(void) {}
202
239
    //@}
203
240
  };
204
241
 
207
244
//------------------------------------------------------------------------------------------
208
245
    
209
246
  template<typename OrdinalType, typename ScalarType>
210
 
  void BLAS<OrdinalType, ScalarType>::ROTG(ScalarType* da, ScalarType* db, MagnitudeType* c, ScalarType* s) const
 
247
  void DefaultBLASImpl<OrdinalType, ScalarType>::ROTG(ScalarType* da, ScalarType* db, MagnitudeType* c, ScalarType* s) const
211
248
  {
212
 
    ScalarType scale, r;
213
 
    ScalarType roe = ScalarTraits<ScalarType>::zero();
 
249
    ScalarType roe, alpha;
214
250
    ScalarType zero = ScalarTraits<ScalarType>::zero();
215
251
    ScalarType one = ScalarTraits<ScalarType>::one();
 
252
    MagnitudeType scale, norm;
 
253
    MagnitudeType m_one = ScalarTraits<MagnitudeType>::one();
 
254
    MagnitudeType m_zero = ScalarTraits<MagnitudeType>::zero();
216
255
 
 
256
    roe = *db;
217
257
    if ( ScalarTraits<ScalarType>::magnitude( *da ) > ScalarTraits<ScalarType>::magnitude( *db ) ) { roe = *da; }
218
258
    scale = ScalarTraits<ScalarType>::magnitude( *da ) + ScalarTraits<ScalarType>::magnitude( *db );
219
 
    if ( scale == zero ) // There is nothing to do.
 
259
    if ( scale == m_zero ) // There is nothing to do.
220
260
    {
221
 
      *c = one;
 
261
      *c = m_one;
222
262
      *s = zero;
223
263
      *da = zero; *db = zero;
224
 
    } else { // Compute the Givens rotation.
225
 
      r = scale*ScalarTraits<ScalarType>::squareroot( ( *da/scale)*(*da/scale) + (*db/scale)*(*db/scale) );
226
 
      if ( roe < zero ) { r *= -one; }
227
 
      *c = *da / r;
228
 
      *s = *db / r;
229
 
      *db = ScalarTraits<ScalarType>::one();
 
264
    } 
 
265
    else if ( *da == zero ) // Still nothing to do.
 
266
    { 
 
267
      *c = m_zero;
 
268
      *s = one;
 
269
      *da = *db; *db = zero;
 
270
    } 
 
271
    else 
 
272
    { // Compute the Givens rotation.
 
273
      norm = scale*ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::squareroot( ( *da/scale)*(*da/scale) + (*db/scale)*(*db/scale) ) );
 
274
      alpha = roe / ScalarTraits<ScalarType>::magnitude(roe);
 
275
      *c = ScalarTraits<ScalarType>::magnitude(*da) / norm;
 
276
      *s = alpha * ScalarTraits<ScalarType>::conjugate(*db) / norm;
 
277
      *db = one;
230
278
      if( ScalarTraits<ScalarType>::magnitude( *da ) > ScalarTraits<ScalarType>::magnitude( *db ) ){ *db = *s; }
231
279
      if( ScalarTraits<ScalarType>::magnitude( *db ) >= ScalarTraits<ScalarType>::magnitude( *da ) &&
232
 
           *c != ScalarTraits<ScalarType>::zero() ) { *db = one / *c; }
233
 
      *da = r;
 
280
           *c != ScalarTraits<MagnitudeType>::zero() ) { *db = one / *c; }
 
281
      *da = norm * alpha;
234
282
    }
235
283
  } /* end ROTG */
236
284
      
237
285
  template<typename OrdinalType, typename ScalarType>
238
 
  void BLAS<OrdinalType,ScalarType>::ROT(const OrdinalType n, ScalarType* dx, const OrdinalType incx, ScalarType* dy, const OrdinalType incy, MagnitudeType* c, ScalarType* s) const
 
286
  void DefaultBLASImpl<OrdinalType,ScalarType>::ROT(const OrdinalType n, ScalarType* dx, const OrdinalType incx, ScalarType* dy, const OrdinalType incy, MagnitudeType* c, ScalarType* s) const
239
287
  {
240
288
    // ToDo:  Implement this.
241
289
  }
242
290
 
243
291
  template<typename OrdinalType, typename ScalarType>
244
 
  void BLAS<OrdinalType, ScalarType>::SCAL(const OrdinalType n, const ScalarType alpha, ScalarType* x, const OrdinalType incx) const
 
292
  void DefaultBLASImpl<OrdinalType, ScalarType>::SCAL(const OrdinalType n, const ScalarType alpha, ScalarType* x, const OrdinalType incx) const
245
293
  {
246
294
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
247
295
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
259
307
  } /* end SCAL */
260
308
 
261
309
  template<typename OrdinalType, typename ScalarType>
262
 
  void BLAS<OrdinalType, ScalarType>::COPY(const OrdinalType n, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
 
310
  void DefaultBLASImpl<OrdinalType, ScalarType>::COPY(const OrdinalType n, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
263
311
  {
264
312
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
265
313
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
279
327
  } /* end COPY */
280
328
 
281
329
  template<typename OrdinalType, typename ScalarType>
282
 
  void BLAS<OrdinalType, ScalarType>::AXPY(const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
 
330
  template <typename alpha_type, typename x_type>
 
331
  void DefaultBLASImpl<OrdinalType, ScalarType>::AXPY(const OrdinalType n, const alpha_type alpha, const x_type* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
283
332
  {
284
333
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
285
334
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
286
335
    OrdinalType i, ix = izero, iy = izero;
287
 
    if( n > izero && alpha != ScalarTraits<ScalarType>::zero())
 
336
    if( n > izero && alpha != ScalarTraits<alpha_type>::zero())
288
337
      {
289
338
        // Set the initial indices (ix, iy).
290
339
        if (incx < izero) { ix = (-n+ione)*incx; }
300
349
  } /* end AXPY */
301
350
 
302
351
  template<typename OrdinalType, typename ScalarType>
303
 
  typename ScalarTraits<ScalarType>::magnitudeType BLAS<OrdinalType, ScalarType>::ASUM(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
 
352
  typename ScalarTraits<ScalarType>::magnitudeType DefaultBLASImpl<OrdinalType, ScalarType>::ASUM(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
304
353
  {
305
354
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
306
355
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
321
370
  } /* end ASUM */
322
371
  
323
372
  template<typename OrdinalType, typename ScalarType>
324
 
  ScalarType BLAS<OrdinalType, ScalarType>::DOT(const OrdinalType n, const ScalarType* x, const OrdinalType incx, const ScalarType* y, const OrdinalType incy) const
 
373
  template <typename x_type, typename y_type>
 
374
  ScalarType DefaultBLASImpl<OrdinalType, ScalarType>::DOT(const OrdinalType n, const x_type* x, const OrdinalType incx, const y_type* y, const OrdinalType incy) const
325
375
  {
326
376
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
327
377
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
335
385
 
336
386
        for(i = izero; i < n; i++)
337
387
          {
338
 
            result += ScalarTraits<ScalarType>::conjugate(x[ix]) * y[iy];
 
388
            result += ScalarTraits<x_type>::conjugate(x[ix]) * y[iy];
339
389
            ix += incx;
340
390
            iy += incy;
341
391
          }
344
394
  } /* end DOT */
345
395
  
346
396
  template<typename OrdinalType, typename ScalarType>
347
 
  typename ScalarTraits<ScalarType>::magnitudeType BLAS<OrdinalType, ScalarType>::NRM2(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
 
397
  typename ScalarTraits<ScalarType>::magnitudeType DefaultBLASImpl<OrdinalType, ScalarType>::NRM2(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
348
398
  {
349
399
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
350
400
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
358
408
    
359
409
        for(i = izero; i < n; i++)
360
410
          {
361
 
            result += ScalarTraits<ScalarType>::conjugate(x[ix]) * x[ix];
 
411
            result += ScalarTraits<ScalarType>::magnitude(ScalarTraits<ScalarType>::conjugate(x[ix]) * x[ix]);
362
412
            ix += incx;
363
413
          }
364
 
        result = ScalarTraits<ScalarType>::squareroot(result);
 
414
        result = ScalarTraits<typename ScalarTraits<ScalarType>::magnitudeType>::squareroot(result);
365
415
      } 
366
416
    return result;
367
417
  } /* end NRM2 */
368
418
  
369
419
  template<typename OrdinalType, typename ScalarType>
370
 
  OrdinalType BLAS<OrdinalType, ScalarType>::IAMAX(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
 
420
  OrdinalType DefaultBLASImpl<OrdinalType, ScalarType>::IAMAX(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
371
421
  {
372
422
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
373
423
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
374
424
    OrdinalType result = izero, ix = izero, i;
375
 
    ScalarType maxval;
 
425
    typename ScalarTraits<ScalarType>::magnitudeType maxval =
 
426
      ScalarTraits<typename ScalarTraits<ScalarType>::magnitudeType>::zero();
376
427
 
377
428
    if ( n > izero ) 
378
429
      {
395
446
//------------------------------------------------------------------------------------------
396
447
//      LEVEL 2 BLAS ROUTINES
397
448
//------------------------------------------------------------------------------------------
398
 
 
399
449
  template<typename OrdinalType, typename ScalarType>
400
 
  void BLAS<OrdinalType, ScalarType>::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, const ScalarType* x, const OrdinalType incx, const ScalarType beta, ScalarType* y, const OrdinalType incy) const
 
450
  template <typename alpha_type, typename A_type, typename x_type, typename beta_type>
 
451
  void DefaultBLASImpl<OrdinalType, ScalarType>::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, const OrdinalType lda, const x_type* x, const OrdinalType incx, const beta_type beta, ScalarType* y, const OrdinalType incy) const
401
452
  {
402
453
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
403
454
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
404
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
405
 
    ScalarType one = ScalarTraits<ScalarType>::one();
 
455
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
456
    beta_type beta_zero = ScalarTraits<beta_type>::zero();
 
457
    x_type x_zero = ScalarTraits<x_type>::zero();
 
458
    ScalarType y_zero = ScalarTraits<ScalarType>::zero();
 
459
    beta_type beta_one = ScalarTraits<beta_type>::one();
406
460
    bool BadArgument = false;
407
461
 
 
462
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex && trans == CONJ_TRANS, std::logic_error,
 
463
        "Teuchos::BLAS::GEMV() does not currently support CONJ_TRANS for complex data types.");
 
464
 
408
465
    // Quick return if there is nothing to do!
409
 
    if( m == izero || n == izero || ( alpha == zero && beta == one ) ){ return; }
 
466
    if( m == izero || n == izero || ( alpha == alpha_zero && beta == beta_one ) ){ return; }
410
467
    
411
468
    // Otherwise, we need to check the argument list.
412
469
    if( m < izero ) { 
450
507
 
451
508
      // Form y = beta*y
452
509
      ix = kx; iy = ky;
453
 
      if(beta != one) {
 
510
      if(beta != beta_one) {
454
511
        if (incy == ione) {
455
 
          if (beta == zero) {
456
 
            for(i = izero; i < leny; i++) { y[i] = zero; }
 
512
          if (beta == beta_zero) {
 
513
            for(i = izero; i < leny; i++) { y[i] = y_zero; }
457
514
          } else {
458
515
            for(i = izero; i < leny; i++) { y[i] *= beta; }
459
516
          }
460
517
        } else {
461
 
          if (beta == zero) {
 
518
          if (beta == beta_zero) {
462
519
            for(i = izero; i < leny; i++) {
463
 
              y[iy] = zero;
 
520
              y[iy] = y_zero;
464
521
              iy += incy;
465
522
            }
466
523
          } else {
473
530
      }
474
531
        
475
532
      // Return if we don't have to do anything more.
476
 
      if(alpha == zero) { return; }
 
533
      if(alpha == alpha_zero) { return; }
477
534
 
478
535
      if( ETranspChar[trans] == 'N' ) {
479
536
        // Form y = alpha*A*y
480
537
        jx = kx;
481
538
        if (incy == ione) {
482
539
          for(j = izero; j < n; j++) {
483
 
            if (x[jx] != zero) {
 
540
            if (x[jx] != x_zero) {
484
541
              temp = alpha*x[jx];
485
542
              for(i = izero; i < m; i++) {
486
543
                y[i] += temp*A[j*lda + i];
490
547
          }
491
548
        } else {
492
549
          for(j = izero; j < n; j++) {
493
 
            if (x[jx] != zero) {
 
550
            if (x[jx] != x_zero) {
494
551
              temp = alpha*x[jx];
495
552
              iy = ky;
496
553
              for(i = izero; i < m; i++) {
505
562
        jy = ky;
506
563
        if (incx == ione) {
507
564
          for(j = izero; j < n; j++) {
508
 
            temp = zero;
 
565
            temp = y_zero;
509
566
            for(i = izero; i < m; i++) {
510
567
              temp += A[j*lda + i]*x[i];
511
568
            }
514
571
          }
515
572
        } else {
516
573
          for(j = izero; j < n; j++) {
517
 
            temp = zero;
 
574
            temp = y_zero;
518
575
            ix = kx;
519
576
            for (i = izero; i < m; i++) {
520
577
              temp += A[j*lda + i]*x[ix];
529
586
  } /* end GEMV */
530
587
 
531
588
 template<typename OrdinalType, typename ScalarType>
532
 
 void BLAS<OrdinalType, ScalarType>::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const ScalarType* A, const OrdinalType lda, ScalarType* x, const OrdinalType incx) const
 
589
 template <typename A_type>
 
590
 void DefaultBLASImpl<OrdinalType, ScalarType>::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const A_type* A, const OrdinalType lda, ScalarType* x, const OrdinalType incx) const
533
591
  {
534
592
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
535
593
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
536
594
    ScalarType zero = ScalarTraits<ScalarType>::zero();
537
595
    bool BadArgument = false;
538
596
 
 
597
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex && trans == CONJ_TRANS, std::logic_error,
 
598
            "Teuchos::BLAS::TRMV() does not currently support CONJ_TRANS for complex data types.");
 
599
 
539
600
    // Quick return if there is nothing to do!
540
601
    if( n == izero ){ return; }
541
602
    
685
746
  } /* end TRMV */
686
747
        
687
748
  template<typename OrdinalType, typename ScalarType>
688
 
  void BLAS<OrdinalType, ScalarType>::GER(const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, const ScalarType* y, const OrdinalType incy, ScalarType* A, const OrdinalType lda) const
 
749
  template <typename alpha_type, typename x_type, typename y_type>
 
750
  void DefaultBLASImpl<OrdinalType, ScalarType>::GER(const OrdinalType m, const OrdinalType n, const alpha_type alpha, const x_type* x, const OrdinalType incx, const y_type* y, const OrdinalType incy, ScalarType* A, const OrdinalType lda) const
689
751
  {
690
752
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
691
753
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
692
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
 
754
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
755
    y_type y_zero = ScalarTraits<y_type>::zero();
693
756
    bool BadArgument = false;
694
757
 
 
758
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex, std::logic_error,
 
759
            "Teuchos::BLAS::GER() does not currently support complex data types.");
 
760
 
695
761
    // Quick return if there is nothing to do!
696
 
    if( m == izero || n == izero || alpha == zero ){ return; }
 
762
    if( m == izero || n == izero || alpha == alpha_zero ){ return; }
697
763
    
698
764
    // Otherwise, we need to check the argument list.
699
765
    if( m < izero ) { 
728
794
      // Start the operations for incx == 1
729
795
      if( incx == ione ) {
730
796
        for( j=izero; j<n; j++ ) {
731
 
          if ( y[jy] != zero ) {
 
797
          if ( y[jy] != y_zero ) {
732
798
            temp = alpha*y[jy];
733
799
            for ( i=izero; i<m; i++ ) {
734
800
              A[j*lda + i] += x[i]*temp;
739
805
      } 
740
806
      else { // Start the operations for incx != 1
741
807
        for( j=izero; j<n; j++ ) {
742
 
          if ( y[jy] != zero ) {
 
808
          if ( y[jy] != y_zero ) {
743
809
            temp = alpha*y[jy];
744
810
            ix = kx;
745
811
            for( i=izero; i<m; i++ ) {
758
824
//------------------------------------------------------------------------------------------
759
825
        
760
826
  template<typename OrdinalType, typename ScalarType>
761
 
  void BLAS<OrdinalType, ScalarType>::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, const ScalarType* B, const OrdinalType ldb, const ScalarType beta, ScalarType* C, const OrdinalType ldc) const
 
827
  template <typename alpha_type, typename A_type, typename B_type, typename beta_type>
 
828
  void DefaultBLASImpl<OrdinalType, ScalarType>::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const alpha_type alpha, const A_type* A, const OrdinalType lda, const B_type* B, const OrdinalType ldb, const beta_type beta, ScalarType* C, const OrdinalType ldc) const
762
829
  {
 
830
 
 
831
    typedef TypeNameTraits<OrdinalType> OTNT;
 
832
    typedef TypeNameTraits<ScalarType> STNT;
 
833
 
763
834
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
764
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
765
 
    ScalarType one = ScalarTraits<ScalarType>::one();
 
835
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
836
    beta_type beta_zero = ScalarTraits<beta_type>::zero();
 
837
    B_type B_zero = ScalarTraits<B_type>::zero();
 
838
    ScalarType C_zero = ScalarTraits<ScalarType>::zero();
 
839
    beta_type beta_one = ScalarTraits<beta_type>::one();
766
840
    OrdinalType i, j, p;
767
841
    OrdinalType NRowA = m, NRowB = k;
768
842
    ScalarType temp;
769
843
    bool BadArgument = false;
770
844
 
 
845
    TEST_FOR_EXCEPTION(
 
846
      Teuchos::ScalarTraits<ScalarType>::isComplex
 
847
      && (transa == CONJ_TRANS || transb == CONJ_TRANS),
 
848
      std::logic_error,
 
849
            "Teuchos::BLAS<"<<OTNT::name()<<","<<STNT::name()<<">::GEMM()"
 
850
      " does not currently support CONJ_TRANS for complex data types."
 
851
      );
 
852
 
771
853
    // Change dimensions of matrix if either matrix is transposed
772
854
    if( !(ETranspChar[transa]=='N') ) {
773
855
      NRowA = k;
777
859
    }
778
860
 
779
861
    // Quick return if there is nothing to do!
780
 
    if( (m==izero) || (n==izero) || (((alpha==zero)||(k==izero)) && (beta==one)) ){ return; }
 
862
    if( (m==izero) || (n==izero) || (((alpha==alpha_zero)||(k==izero)) && (beta==beta_one)) ){ return; }
781
863
    if( m < izero ) { 
782
864
      std::cout << "BLAS::GEMM Error: M == " << m << std::endl;     
783
865
      BadArgument = true;
806
888
    if(!BadArgument) {
807
889
 
808
890
      // Only need to scale the resulting matrix C.
809
 
      if( alpha == zero ) {
810
 
        if( beta == zero ) {
 
891
      if( alpha == alpha_zero ) {
 
892
        if( beta == beta_zero ) {
811
893
          for (j=izero; j<n; j++) {
812
894
            for (i=izero; i<m; i++) {
813
 
              C[j*ldc + i] = zero;
 
895
              C[j*ldc + i] = C_zero;
814
896
            }
815
897
          }
816
898
        } else {
829
911
        if ( ETranspChar[transa]=='N' ) {
830
912
          // Compute C = alpha*A*B + beta*C
831
913
          for (j=izero; j<n; j++) {
832
 
            if( beta == zero ) {
 
914
            if( beta == beta_zero ) {
833
915
              for (i=izero; i<m; i++){
834
 
                C[j*ldc + i] = zero;
 
916
                C[j*ldc + i] = C_zero;
835
917
              }
836
 
            } else if( beta != one ) {
 
918
            } else if( beta != beta_one ) {
837
919
              for (i=izero; i<m; i++){
838
920
                C[j*ldc + i] *= beta;
839
921
              }
840
922
            }
841
923
            for (p=izero; p<k; p++){
842
 
              if (B[j*ldb + p] != zero ){
 
924
              if (B[j*ldb + p] != B_zero ){
843
925
                temp = alpha*B[j*ldb + p];
844
926
                for (i=izero; i<m; i++) {
845
927
                  C[j*ldc + i] += temp*A[p*lda + i];
851
933
          // Compute C = alpha*A'*B + beta*C
852
934
          for (j=izero; j<n; j++) {
853
935
            for (i=izero; i<m; i++) {
854
 
              temp = zero;
 
936
              temp = C_zero;
855
937
              for (p=izero; p<k; p++) {
856
938
                temp += A[i*lda + p]*B[j*ldb + p];
857
939
              }
858
 
              if (beta == zero) {
 
940
              if (beta == beta_zero) {
859
941
                C[j*ldc + i] = alpha*temp;
860
942
              } else {
861
943
                C[j*ldc + i] = alpha*temp + beta*C[j*ldc + i];
867
949
        if ( ETranspChar[transa]=='N' ) {
868
950
          // Compute C = alpha*A*B' + beta*C
869
951
          for (j=izero; j<n; j++) {
870
 
            if (beta == zero) {
 
952
            if (beta == beta_zero) {
871
953
              for (i=izero; i<m; i++) {
872
 
                C[j*ldc + i] = zero;
 
954
                C[j*ldc + i] = C_zero;
873
955
              } 
874
 
            } else if ( beta != one ) {
 
956
            } else if ( beta != beta_one ) {
875
957
              for (i=izero; i<m; i++) {
876
958
                C[j*ldc + i] *= beta;
877
959
              }
878
960
            }
879
961
            for (p=izero; p<k; p++) {
880
 
              if (B[p*ldb + j] != zero) {
 
962
              if (B[p*ldb + j] != B_zero) {
881
963
                temp = alpha*B[p*ldb + j];
882
964
                for (i=izero; i<m; i++) {
883
965
                  C[j*ldc + i] += temp*A[p*lda + i];
889
971
          // Compute C += alpha*A'*B' + beta*C
890
972
          for (j=izero; j<n; j++) {
891
973
            for (i=izero; i<m; i++) {
892
 
              temp = zero;
 
974
              temp = C_zero;
893
975
              for (p=izero; p<k; p++) {
894
976
                temp += A[i*lda + p]*B[p*ldb + j];
895
977
              }
896
 
              if (beta == zero) {
 
978
              if (beta == beta_zero) {
897
979
                C[j*ldc + i] = alpha*temp;
898
980
              } else {
899
981
                C[j*ldc + i] = alpha*temp + beta*C[j*ldc + i];
907
989
 
908
990
 
909
991
  template<typename OrdinalType, typename ScalarType>
910
 
  void BLAS<OrdinalType, ScalarType>::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, const ScalarType* B, const OrdinalType ldb, const ScalarType beta, ScalarType* C, const OrdinalType ldc) const
 
992
  template <typename alpha_type, typename A_type, typename B_type, typename beta_type>
 
993
  void DefaultBLASImpl<OrdinalType, ScalarType>::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, const OrdinalType lda, const B_type* B, const OrdinalType ldb, const beta_type beta, ScalarType* C, const OrdinalType ldc) const
911
994
  {
912
995
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
913
996
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
914
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
915
 
    ScalarType one = ScalarTraits<ScalarType>::one();
 
997
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
998
    beta_type beta_zero = ScalarTraits<beta_type>::zero();
 
999
    ScalarType C_zero = ScalarTraits<ScalarType>::zero();
 
1000
    beta_type beta_one = ScalarTraits<beta_type>::one();
916
1001
    OrdinalType i, j, k, NRowA = m;
917
1002
    ScalarType temp1, temp2;
918
1003
    bool BadArgument = false;
919
1004
    bool Upper = (EUploChar[uplo] == 'U');
920
1005
    if (ESideChar[side] == 'R') { NRowA = n; }
 
1006
 
 
1007
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex, std::logic_error,
 
1008
            "Teuchos::BLAS::SYMM() does not currently support complex data types.");
921
1009
    
922
1010
    // Quick return.
923
 
    if ( (m==izero) || (n==izero) || ( (alpha==zero)&&(beta==one) ) ) { return; }
 
1011
    if ( (m==izero) || (n==izero) || ( (alpha==alpha_zero)&&(beta==beta_one) ) ) { return; }
924
1012
    if( m < 0 ) { 
925
1013
      std::cout << "BLAS::SYMM Error: M == "<< m << std::endl;
926
1014
      BadArgument = true; }
940
1028
    if(!BadArgument) {
941
1029
 
942
1030
      // Only need to scale C and return.
943
 
      if (alpha == zero) {
944
 
        if (beta == zero ) {
 
1031
      if (alpha == alpha_zero) {
 
1032
        if (beta == beta_zero ) {
945
1033
          for (j=izero; j<n; j++) {
946
1034
            for (i=izero; i<m; i++) {
947
 
              C[j*ldc + i] = zero;
 
1035
              C[j*ldc + i] = C_zero;
948
1036
            }
949
1037
          }
950
1038
        } else {
965
1053
          for (j=izero; j<n; j++) {
966
1054
            for (i=izero; i<m; i++) {
967
1055
              temp1 = alpha*B[j*ldb + i];
968
 
              temp2 = zero;
 
1056
              temp2 = C_zero;
969
1057
              for (k=izero; k<i; k++) {
970
1058
                C[j*ldc + k] += temp1*A[i*lda + k];
971
1059
                temp2 += B[j*ldb + k]*A[i*lda + k];
972
1060
              }
973
 
              if (beta == zero) {
 
1061
              if (beta == beta_zero) {
974
1062
                C[j*ldc + i] = temp1*A[i*lda + i] + alpha*temp2;
975
1063
              } else {
976
1064
                C[j*ldc + i] = beta*C[j*ldc + i] + temp1*A[i*lda + i] + alpha*temp2;
982
1070
          for (j=izero; j<n; j++) {
983
1071
            for (i=m-ione; i>-ione; i--) {
984
1072
              temp1 = alpha*B[j*ldb + i];
985
 
              temp2 = zero;
 
1073
              temp2 = C_zero;
986
1074
              for (k=i+ione; k<m; k++) {
987
1075
                C[j*ldc + k] += temp1*A[i*lda + k];
988
1076
                temp2 += B[j*ldb + k]*A[i*lda + k];
989
1077
              }
990
 
              if (beta == zero) {
 
1078
              if (beta == beta_zero) {
991
1079
                C[j*ldc + i] = temp1*A[i*lda + i] + alpha*temp2;
992
1080
              } else {
993
1081
                C[j*ldc + i] = beta*C[j*ldc + i] + temp1*A[i*lda + i] + alpha*temp2;
999
1087
        // Compute C = alpha*B*A + beta*C.
1000
1088
        for (j=izero; j<n; j++) {
1001
1089
          temp1 = alpha*A[j*lda + j];
1002
 
          if (beta == zero) {
 
1090
          if (beta == beta_zero) {
1003
1091
            for (i=izero; i<m; i++) {
1004
1092
              C[j*ldc + i] = temp1*B[j*ldb + i];
1005
1093
            }
1034
1122
} // end SYMM
1035
1123
  
1036
1124
  template<typename OrdinalType, typename ScalarType>
1037
 
  void BLAS<OrdinalType, ScalarType>::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const
 
1125
  template <typename alpha_type, typename A_type>
 
1126
  void DefaultBLASImpl<OrdinalType, ScalarType>::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const
1038
1127
  {
1039
1128
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
1040
1129
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
1041
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
 
1130
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
1131
    A_type A_zero = ScalarTraits<A_type>::zero();
 
1132
    ScalarType B_zero = ScalarTraits<ScalarType>::zero();
1042
1133
    ScalarType one = ScalarTraits<ScalarType>::one();
1043
1134
    OrdinalType i, j, k, NRowA = m;
1044
1135
    ScalarType temp;
1047
1138
    bool NoUnit = (EDiagChar[diag] == 'N');
1048
1139
    bool Upper = (EUploChar[uplo] == 'U');
1049
1140
 
 
1141
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex && transa == CONJ_TRANS, std::logic_error,
 
1142
            "Teuchos::BLAS::TRMM() does not currently support CONJ_TRANS for complex data types.");
 
1143
 
1050
1144
    if(!LSide) { NRowA = n; }
1051
1145
 
1052
1146
    // Quick return.
1067
1161
    if(!BadArgument) {
1068
1162
 
1069
1163
      // B only needs to be zeroed out.
1070
 
      if( alpha == zero ) {
 
1164
      if( alpha == alpha_zero ) {
1071
1165
        for( j=izero; j<n; j++ ) {
1072
1166
          for( i=izero; i<m; i++ ) {
1073
 
            B[j*ldb + i] = zero;
 
1167
            B[j*ldb + i] = B_zero;
1074
1168
          }
1075
1169
        }
1076
1170
        return;
1087
1181
            // A is upper triangular
1088
1182
            for( j=izero; j<n; j++ ) {
1089
1183
              for( k=izero; k<m; k++) {
1090
 
                if ( B[j*ldb + k] != zero ) {
 
1184
                if ( B[j*ldb + k] != B_zero ) {
1091
1185
                  temp = alpha*B[j*ldb + k];
1092
1186
                  for( i=izero; i<k; i++ ) {
1093
1187
                    B[j*ldb + i] += temp*A[k*lda + i];
1102
1196
            // A is lower triangular
1103
1197
            for( j=izero; j<n; j++ ) {
1104
1198
              for( k=m-ione; k>-ione; k-- ) {
1105
 
                if( B[j*ldb + k] != zero ) {
 
1199
                if( B[j*ldb + k] != B_zero ) {
1106
1200
                  temp = alpha*B[j*ldb + k];
1107
1201
                  B[j*ldb + k] = temp;
1108
1202
                  if ( NoUnit )
1158
1252
                B[j*ldb + i] *= temp;
1159
1253
              }
1160
1254
              for( k=izero; k<j; k++ ) {
1161
 
                if( A[j*lda + k] != zero ) {
 
1255
                if( A[j*lda + k] != A_zero ) {
1162
1256
                  temp = alpha*A[j*lda + k];
1163
1257
                  for( i=izero; i<m; i++ ) {
1164
1258
                    B[j*ldb + i] += temp*B[k*ldb + i];
1176
1270
                B[j*ldb + i] *= temp;
1177
1271
              }
1178
1272
              for( k=j+ione; k<n; k++ ) {
1179
 
                if( A[j*lda + k] != zero ) {
 
1273
                if( A[j*lda + k] != A_zero ) {
1180
1274
                  temp = alpha*A[j*lda + k];
1181
1275
                  for( i=izero; i<m; i++ ) {
1182
1276
                    B[j*ldb + i] += temp*B[k*ldb + i];
1191
1285
          if( Upper ) {
1192
1286
            for( k=izero; k<n; k++ ) {
1193
1287
              for( j=izero; j<k; j++ ) {
1194
 
                if( A[k*lda + j] != zero ) {
 
1288
                if( A[k*lda + j] != A_zero ) {
1195
1289
                  temp = alpha*A[k*lda + j];
1196
1290
                  for( i=izero; i<m; i++ ) {
1197
1291
                    B[j*ldb + i] += temp*B[k*ldb + i];
1210
1304
          } else {
1211
1305
            for( k=n-ione; k>-ione; k-- ) {
1212
1306
              for( j=k+ione; j<n; j++ ) {
1213
 
                if( A[k*lda + j] != zero ) {
 
1307
                if( A[k*lda + j] != A_zero ) {
1214
1308
                  temp = alpha*A[k*lda + j];
1215
1309
                  for( i=izero; i<m; i++ ) {
1216
1310
                    B[j*ldb + i] += temp*B[k*ldb + i];
1233
1327
  } // end TRMM
1234
1328
  
1235
1329
  template<typename OrdinalType, typename ScalarType>
1236
 
  void BLAS<OrdinalType, ScalarType>::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const
 
1330
  template <typename alpha_type, typename A_type>
 
1331
  void DefaultBLASImpl<OrdinalType, ScalarType>::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const alpha_type alpha, const A_type* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const
1237
1332
  {
1238
1333
    OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
1239
1334
    OrdinalType ione = OrdinalTraits<OrdinalType>::one();
1240
 
    ScalarType zero = ScalarTraits<ScalarType>::zero();
1241
 
    ScalarType one = ScalarTraits<ScalarType>::one();
 
1335
    alpha_type alpha_zero = ScalarTraits<alpha_type>::zero();
 
1336
    A_type A_zero = ScalarTraits<A_type>::zero();
 
1337
    ScalarType B_zero = ScalarTraits<ScalarType>::zero();
 
1338
    alpha_type alpha_one = ScalarTraits<alpha_type>::one();
 
1339
    ScalarType B_one = ScalarTraits<ScalarType>::one();
1242
1340
    ScalarType temp;
1243
1341
    OrdinalType NRowA = m;
1244
1342
    bool BadArgument = false;
1245
1343
    bool NoUnit = (EDiagChar[diag]=='N');
1246
1344
    
 
1345
    TEST_FOR_EXCEPTION(Teuchos::ScalarTraits<ScalarType>::isComplex && transa == CONJ_TRANS, std::logic_error,
 
1346
            "Teuchos::BLAS::TRSM() does not currently support CONJ_TRANS for complex data types.");
 
1347
 
1247
1348
    if (!(ESideChar[side] == 'L')) { NRowA = n; }
1248
1349
 
1249
1350
    // Quick return.
1265
1366
      {
1266
1367
        int i, j, k;
1267
1368
        // Set the solution to the zero std::vector.
1268
 
        if(alpha == zero) {
 
1369
        if(alpha == alpha_zero) {
1269
1370
            for(j = izero; j < n; j++) {
1270
1371
                for( i = izero; i < m; i++) {
1271
 
                    B[j*ldb+i] = zero;
 
1372
                    B[j*ldb+i] = B_zero;
1272
1373
                }
1273
1374
            }
1274
1375
        }
1286
1387
                        // A is upper triangular.
1287
1388
                        for(j = izero; j < n; j++) {
1288
1389
                            // Perform alpha*B if alpha is not 1.
1289
 
                            if(alpha != one) {
 
1390
                          if(alpha != alpha_one) {
1290
1391
                                for( i = izero; i < m; i++) {
1291
1392
                                    B[j*ldb+i] *= alpha;
1292
1393
                                }
1294
1395
                            // Perform a backsolve for column j of B.
1295
1396
                            for(k = (m - ione); k > -ione; k--) {
1296
1397
                                // If this entry is zero, we don't have to do anything.
1297
 
                                if (B[j*ldb + k] != zero) {
 
1398
                                if (B[j*ldb + k] != B_zero) {
1298
1399
                                    if (NoUnit) {
1299
1400
                                        B[j*ldb + k] /= A[k*lda + k];
1300
1401
                                    }
1309
1410
                    { // A is lower triangular.
1310
1411
                        for(j = izero; j < n; j++) {
1311
1412
                            // Perform alpha*B if alpha is not 1.
1312
 
                            if(alpha != one) {
 
1413
                            if(alpha != alpha_one) {
1313
1414
                                for( i = izero; i < m; i++) {
1314
1415
                                    B[j*ldb+i] *= alpha;
1315
1416
                                }
1317
1418
                            // Perform a forward solve for column j of B.
1318
1419
                            for(k = izero; k < m; k++) {
1319
1420
                                // If this entry is zero, we don't have to do anything.
1320
 
                                if (B[j*ldb + k] != zero) {   
 
1421
                                if (B[j*ldb + k] != B_zero) {   
1321
1422
                                    if (NoUnit) {
1322
1423
                                        B[j*ldb + k] /= A[k*lda + k];
1323
1424
                                    }
1379
1480
                        // Perform a backsolve for column j of B.
1380
1481
                        for(j = izero; j < n; j++) {
1381
1482
                            // Perform alpha*B if alpha is not 1.
1382
 
                            if(alpha != one) {
 
1483
                            if(alpha != alpha_one) {
1383
1484
                                for( i = izero; i < m; i++) {
1384
1485
                                    B[j*ldb+i] *= alpha;
1385
1486
                                }
1386
1487
                            }
1387
1488
                            for(k = izero; k < j; k++) {
1388
1489
                                // If this entry is zero, we don't have to do anything.
1389
 
                                if (A[j*lda + k] != zero) {
 
1490
                                if (A[j*lda + k] != A_zero) {
1390
1491
                                    for(i = izero; i < m; i++) {
1391
1492
                                        B[j*ldb + i] -= A[j*lda + k] * B[k*ldb + i];
1392
1493
                                    }
1393
1494
                                }
1394
1495
                            }
1395
1496
                            if (NoUnit) {
1396
 
                                temp = one/A[j*lda + j];
 
1497
                                temp = B_one/A[j*lda + j];
1397
1498
                                for(i = izero; i < m; i++) {
1398
1499
                                    B[j*ldb + i] *= temp;
1399
1500
                                }
1404
1505
                    { // A is lower triangular.
1405
1506
                        for(j = (n - ione); j > -ione; j--) {
1406
1507
                            // Perform alpha*B if alpha is not 1.
1407
 
                            if(alpha != one) {
 
1508
                            if(alpha != alpha_one) {
1408
1509
                                for( i = izero; i < m; i++) {
1409
1510
                                    B[j*ldb+i] *= alpha;
1410
1511
                                }
1412
1513
                            // Perform a forward solve for column j of B.
1413
1514
                            for(k = j+ione; k < n; k++) {
1414
1515
                                // If this entry is zero, we don't have to do anything.
1415
 
                                if (A[j*lda + k] != zero) {
 
1516
                                if (A[j*lda + k] != A_zero) {
1416
1517
                                    for(i = izero; i < m; i++) {
1417
1518
                                        B[j*ldb + i] -= A[j*lda + k] * B[k*ldb + i]; 
1418
1519
                                    }
1419
1520
                                } 
1420
1521
                            }
1421
1522
                            if (NoUnit) {
1422
 
                                temp = one/A[j*lda + j];
 
1523
                                temp = B_one/A[j*lda + j];
1423
1524
                                for(i = izero; i < m; i++) {
1424
1525
                                    B[j*ldb + i] *= temp;
1425
1526
                                }
1435
1536
                        // A is upper triangular.
1436
1537
                        for(k = (n - ione); k > -ione; k--) {
1437
1538
                            if (NoUnit) {
1438
 
                                temp = one/A[k*lda + k];
 
1539
                                temp = B_one/A[k*lda + k];
1439
1540
                                for(i = izero; i < m; i++) {
1440
1541
                                    B[k*ldb + i] *= temp;
1441
1542
                                }
1442
1543
                            }
1443
1544
                            for(j = izero; j < k; j++) {
1444
 
                                if (A[k*lda + j] != zero) {
 
1545
                                if (A[k*lda + j] != A_zero) {
1445
1546
                                    temp = A[k*lda + j];
1446
1547
                                    for(i = izero; i < m; i++) {
1447
1548
                                        B[j*ldb + i] -= temp*B[k*ldb + i];
1448
1549
                                    }
1449
1550
                                }
1450
1551
                            }
1451
 
                            if (alpha != one) {
 
1552
                            if (alpha != alpha_one) {
1452
1553
                                for (i = izero; i < m; i++) {
1453
1554
                                    B[k*ldb + i] *= alpha;
1454
1555
                                }
1459
1560
                    { // A is lower triangular.
1460
1561
                        for(k = izero; k < n; k++) {
1461
1562
                            if (NoUnit) {
1462
 
                                temp = one/A[k*lda + k];
 
1563
                                temp = B_one/A[k*lda + k];
1463
1564
                                for (i = izero; i < m; i++) {
1464
1565
                                    B[k*ldb + i] *= temp;
1465
1566
                                }
1466
1567
                            }
1467
1568
                            for(j = k+ione; j < n; j++) {
1468
 
                                if(A[k*lda + j] != zero) {
 
1569
                                if(A[k*lda + j] != A_zero) {
1469
1570
                                    temp = A[k*lda + j];
1470
1571
                                    for(i = izero; i < m; i++) {
1471
1572
                                        B[j*ldb + i] -= temp*B[k*ldb + i];
1472
1573
                                    }
1473
1574
                                }
1474
1575
                            }
1475
 
                            if (alpha != one) {
 
1576
                            if (alpha != alpha_one) {
1476
1577
                                for (i = izero; i < m; i++) {
1477
1578
                                    B[k*ldb + i] *= alpha;
1478
1579
                                }
1487
1588
 
1488
1589
  // Explicit instantiation for template<int,float>
1489
1590
 
1490
 
#ifdef HAVE_TEUCHOS_BLASFLOAT
1491
 
 
1492
1591
  template <>
1493
1592
  class BLAS<int, float>
1494
1593
  {    
1514
1613
    void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const int m, const int n, const float alpha, const float* A, const int lda, float* B, const int ldb) const;
1515
1614
  };
1516
1615
 
1517
 
#endif // HAVE_TEUCHOS_BLASFLOAT
1518
 
 
1519
1616
  // Explicit instantiation for template<int,double>
1520
1617
 
1521
1618
  template<>
1543
1640
    void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const int m, const int n, const double alpha, const double* A, const int lda, double* B, const int ldb) const;
1544
1641
  };
1545
1642
 
1546
 
#ifdef HAVE_TEUCHOS_BLASFLOAT
1547
 
 
1548
1643
  // Explicit instantiation for template<int,complex<float> >
1549
1644
 
1550
1645
  template<>
1572
1667
    void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const int m, const int n, const std::complex<float> alpha, const std::complex<float>* A, const int lda, std::complex<float>* B, const int ldb) const;
1573
1668
  };
1574
1669
 
1575
 
#endif // HAVE_TEUCHOS_BLASFLOAT
1576
 
 
1577
1670
  // Explicit instantiation for template<int,complex<double> >
1578
1671
 
1579
1672
  template<>