Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members | Examples | Search for 
src

Teuchos_BLAS.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //                    Teuchos: Common Tools Package
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 // Kris
00030 // 06.16.03 -- Start over from scratch
00031 // 06.16.03 -- Initial templatization (Tpetra_BLAS.cpp is no longer needed)
00032 // 06.18.03 -- Changed xxxxx_() function calls to XXXXX_F77()
00033 //          -- Added warning messages for generic calls
00034 // 07.08.03 -- Move into Teuchos package/namespace
00035 // 07.24.03 -- The first iteration of BLAS generics is nearing completion. Caveats:
00036 //             * TRSM isn't finished yet; it works for one case at the moment (left side, upper tri., no transpose, no unit diag.)
00037 //             * Many of the generic implementations are quite inefficient, ugly, or both. I wrote these to be easy to debug, not for efficiency or legibility. The next iteration will improve both of these aspects as much as possible.
00038 //             * Very little verification of input parameters is done, save for the character-type arguments (TRANS, etc.) which is quite robust.
00039 //             * All of the routines that make use of both an incx and incy parameter (which includes much of the L1 BLAS) are set up to work iff incx == incy && incx > 0. Allowing for differing/negative values of incx/incy should be relatively trivial.
00040 //             * All of the L2/L3 routines assume that the entire matrix is being used (that is, if A is mxn, lda = m); they don't work on submatrices yet. This *should* be a reasonably trivial thing to fix, as well.
00041 //          -- Removed warning messages for generic calls
00042 // 08.08.03 -- TRSM now works for all cases where SIDE == L and DIAG == N. DIAG == U is implemented but does not work correctly; SIDE == R is not yet implemented.
00043 // 08.14.03 -- TRSM now works for all cases and accepts (and uses) leading-dimension information.
00044 // 09.26.03 -- character input replaced with enumerated input to cause compiling errors and not run-time errors ( suggested by RAB ).
00045 
00046 #ifndef _TEUCHOS_BLAS_HPP_
00047 #define _TEUCHOS_BLAS_HPP_
00048 
00056 /* for INTEL_CXML, the second arg may need to be changed to 'one'.  If so
00057 the appropriate declaration of one will need to be added back into
00058 functions that include the macro:
00059 */
00060 #if defined (INTEL_CXML)
00061         unsigned int one=1;
00062 #endif
00063 
00064 #ifdef CHAR_MACRO
00065 #undef CHAR_MACRO
00066 #endif
00067 #if defined (INTEL_CXML)
00068 #define CHAR_MACRO(char_var) &char_var, one
00069 #else
00070 #define CHAR_MACRO(char_var) &char_var
00071 #endif
00072 
00073 #include "Teuchos_ConfigDefs.hpp"
00074 #include "Teuchos_BLAS_wrappers.hpp"
00075 #include "Teuchos_BLAS_types.hpp"
00076 #include "Teuchos_ScalarTraits.hpp"
00077 #include "Teuchos_OrdinalTraits.hpp"
00078 
00112 namespace Teuchos
00113 {
00114   extern const char ESideChar[];
00115   extern const char ETranspChar[];
00116   extern const char EUploChar[];
00117   extern const char EDiagChar[];
00118 
00119   template<typename OrdinalType, typename ScalarType>
00120   class BLAS
00121   {    
00122 
00123     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00124     
00125   public:
00127 
00128     
00130     inline BLAS(void) {};
00131 
00133     inline BLAS(const BLAS<OrdinalType, ScalarType>& BLAS_source) {};
00134 
00136     inline virtual ~BLAS(void) {};
00138 
00140 
00141 
00143     void ROTG(ScalarType* da, ScalarType* db, MagnitudeType* c, ScalarType* s) const;
00144 
00146     void ROT(const OrdinalType n, ScalarType* dx, const OrdinalType incx, ScalarType* dy, const OrdinalType incy, MagnitudeType* c, ScalarType* s) const;
00147 
00149     void SCAL(const OrdinalType n, const ScalarType alpha, ScalarType* x, const OrdinalType incx) const;
00150 
00152     void COPY(const OrdinalType n, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const;
00153 
00155     void AXPY(const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const;
00156 
00158     typename ScalarTraits<ScalarType>::magnitudeType ASUM(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const;
00159 
00161     ScalarType DOT(const OrdinalType n, const ScalarType* x, const OrdinalType incx, const ScalarType* y, const OrdinalType incy) const;
00162 
00164     typename ScalarTraits<ScalarType>::magnitudeType NRM2(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const;
00165 
00167     OrdinalType IAMAX(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const;
00168 
00170 
00172 
00173 
00175     void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* A, 
00176         const OrdinalType lda, const ScalarType* x, const OrdinalType incx, const ScalarType beta, ScalarType* y, const OrdinalType incy) const;
00177 
00179     void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const ScalarType* A, 
00180         const OrdinalType lda, ScalarType* x, const OrdinalType incx) const;
00181 
00183     void GER(const OrdinalType m, const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, 
00184        const ScalarType* y, const OrdinalType incy, ScalarType* A, const OrdinalType lda) const;
00186     
00188 
00189 
00191     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;
00192 
00194     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;
00195 
00197     void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n,
00198                 const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
00199 
00201     void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n,
00202                 const ScalarType alpha, const ScalarType* A, const OrdinalType lda, ScalarType* B, const OrdinalType ldb) const;
00204   };
00205 
00206 //------------------------------------------------------------------------------------------
00207 //      LEVEL 1 BLAS ROUTINES  
00208 //------------------------------------------------------------------------------------------
00209     
00210   template<typename OrdinalType, typename ScalarType>
00211   void BLAS<OrdinalType, ScalarType>::ROTG(ScalarType* da, ScalarType* db, MagnitudeType* c, ScalarType* s) const
00212   {
00213     ScalarType roe, scale, r;
00214     ScalarType zero = ScalarTraits<ScalarType>::zero();
00215     ScalarType one = ScalarTraits<ScalarType>::one();
00216 
00217     if ( ScalarTraits<ScalarType>::magnitude( *da ) > ScalarTraits<ScalarType>::magnitude( *db ) ) { roe = *da; }
00218     scale = ScalarTraits<ScalarType>::magnitude( *da ) + ScalarTraits<ScalarType>::magnitude( *db );
00219     if ( scale == zero ) // There is nothing to do.
00220     {
00221       *c = one;
00222       *s = zero;
00223       *da = zero; *db = zero;
00224     } else { // Compute the Givens rotation.
00225       r = scale*ScalarTraits<ScalarType>::squareroot( ( *da/scale)*(*da/scale) + (*db/scale)*(*db/scale) );
00226       if ( roe < zero ) { r *= -one; }
00227       *c = *da / r;
00228       *s = *db / r;
00229       *db = ScalarTraits<ScalarType>::one();
00230       if( ScalarTraits<ScalarType>::magnitude( *da ) > ScalarTraits<ScalarType>::magnitude( *db ) ){ *db = *s; }
00231       if( ScalarTraits<ScalarType>::magnitude( *db ) >= ScalarTraits<ScalarType>::magnitude( *da ) &&
00232      *c != ScalarTraits<ScalarType>::zero() ) { *db = one / *c; }
00233       *da = r;
00234     }
00235   } /* end ROTG */
00236       
00237   template<typename OrdinalType, typename ScalarType>
00238   void BLAS<OrdinalType,ScalarType>::ROT(const OrdinalType n, ScalarType* dx, const OrdinalType incx, ScalarType* dy, const OrdinalType incy, MagnitudeType* c, ScalarType* s) const
00239   {
00240     // ToDo:  Implement this.
00241   }
00242 
00243   template<typename OrdinalType, typename ScalarType>
00244   void BLAS<OrdinalType, ScalarType>::SCAL(const OrdinalType n, const ScalarType alpha, ScalarType* x, const OrdinalType incx) const
00245   {
00246     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00247     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00248     OrdinalType i, ix = izero;
00249     if ( n > izero ) {
00250         // Set the initial index (ix).
00251         if (incx < izero) { ix = (-n+ione)*incx; } 
00252         // Scale the vector.
00253         for(i = izero; i < n; i++)
00254         {
00255             x[ix] *= alpha;
00256             ix += incx;
00257         }
00258     }
00259   } /* end SCAL */
00260 
00261   template<typename OrdinalType, typename ScalarType>
00262   void BLAS<OrdinalType, ScalarType>::COPY(const OrdinalType n, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
00263   {
00264     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00265     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00266     OrdinalType i, ix = izero, iy = izero;
00267     if ( n > izero ) {
00268   // Set the initial indices (ix, iy).
00269       if (incx < izero) { ix = (-n+ione)*incx; }
00270       if (incy < izero) { iy = (-n+ione)*incy; }
00271 
00272         for(i = izero; i < n; i++)
00273           {
00274       y[iy] = x[ix];
00275       ix += incx;
00276       iy += incy;
00277           }
00278       }
00279   } /* end COPY */
00280 
00281   template<typename OrdinalType, typename ScalarType>
00282   void BLAS<OrdinalType, ScalarType>::AXPY(const OrdinalType n, const ScalarType alpha, const ScalarType* x, const OrdinalType incx, ScalarType* y, const OrdinalType incy) const
00283   {
00284     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00285     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00286     OrdinalType i, ix = izero, iy = izero;
00287     if( n > izero && alpha != ScalarTraits<ScalarType>::zero())
00288       {
00289   // Set the initial indices (ix, iy).
00290       if (incx < izero) { ix = (-n+ione)*incx; }
00291       if (incy < izero) { iy = (-n+ione)*incy; }
00292 
00293         for(i = izero; i < n; i++)
00294           {
00295       y[iy] += alpha * x[ix];
00296       ix += incx;
00297       iy += incy;
00298           }
00299       }
00300   } /* end AXPY */
00301 
00302   template<typename OrdinalType, typename ScalarType>
00303   typename ScalarTraits<ScalarType>::magnitudeType BLAS<OrdinalType, ScalarType>::ASUM(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
00304   {
00305     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00306     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00307     typename ScalarTraits<ScalarType>::magnitudeType result = 
00308       ScalarTraits<typename ScalarTraits<ScalarType>::magnitudeType>::zero();
00309     OrdinalType i, ix = izero;
00310     if( n > izero ) {
00311   // Set the initial indices
00312   if (incx < izero) { ix = (-n+ione)*incx; }
00313 
00314       for(i = izero; i < n; i++)
00315           {
00316       result += ScalarTraits<ScalarType>::magnitude(x[ix]);
00317       ix += incx;
00318           }
00319     } 
00320    return result;
00321   } /* end ASUM */
00322   
00323   template<typename OrdinalType, typename ScalarType>
00324   ScalarType BLAS<OrdinalType, ScalarType>::DOT(const OrdinalType n, const ScalarType* x, const OrdinalType incx, const ScalarType* y, const OrdinalType incy) const
00325   {
00326     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00327     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00328     ScalarType result = ScalarTraits<ScalarType>::zero();
00329     OrdinalType i, ix = izero, iy = izero;
00330     if( n > izero ) 
00331       {
00332   // Set the initial indices (ix,iy).       
00333   if (incx < izero) { ix = (-n+ione)*incx; }
00334   if (incy < izero) { iy = (-n+ione)*incy; }
00335 
00336   for(i = izero; i < n; i++)
00337     {
00338       result += ScalarTraits<ScalarType>::conjugate(x[ix]) * y[iy];
00339       ix += incx;
00340       iy += incy;
00341     }
00342       }
00343     return result;
00344   } /* end DOT */
00345   
00346   template<typename OrdinalType, typename ScalarType>
00347   typename ScalarTraits<ScalarType>::magnitudeType BLAS<OrdinalType, ScalarType>::NRM2(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
00348   {
00349     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00350     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00351     typename ScalarTraits<ScalarType>::magnitudeType result = 
00352       ScalarTraits<typename ScalarTraits<ScalarType>::magnitudeType>::zero();
00353     OrdinalType i, ix = izero;
00354     if ( n > izero ) 
00355       {
00356   // Set the initial index.
00357   if (incx < izero) { ix = (-n+ione)*incx; }  
00358     
00359   for(i = izero; i < n; i++)
00360           {
00361       result += ScalarTraits<ScalarType>::conjugate(x[ix]) * x[ix];
00362       ix += incx;
00363           }
00364       result = ScalarTraits<ScalarType>::squareroot(result);
00365       } 
00366     return result;
00367   } /* end NRM2 */
00368   
00369   template<typename OrdinalType, typename ScalarType>
00370   OrdinalType BLAS<OrdinalType, ScalarType>::IAMAX(const OrdinalType n, const ScalarType* x, const OrdinalType incx) const
00371   {
00372     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00373     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00374     OrdinalType result = izero, ix = izero, i;
00375     ScalarType maxval;
00376 
00377     if ( n > izero ) 
00378       {
00379   if (incx < izero) { ix = (-n+ione)*incx; }
00380   maxval = ScalarTraits<ScalarType>::magnitude(x[ix]);
00381   ix += incx;
00382       for(i = ione; i < n; i++)
00383           {
00384       if(ScalarTraits<ScalarType>::magnitude(x[ix]) > maxval)
00385         {
00386         result = i;
00387           maxval = ScalarTraits<ScalarType>::magnitude(x[ix]);
00388         }
00389       ix += incx;
00390     }
00391       }
00392     return result + 1; // the BLAS I?AMAX functions return 1-indexed (Fortran-style) values
00393   } /* end IAMAX */
00394 
00395 //------------------------------------------------------------------------------------------
00396 //      LEVEL 2 BLAS ROUTINES
00397 //------------------------------------------------------------------------------------------
00398 
00399   template<typename OrdinalType, typename ScalarType>
00400   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
00401   {
00402     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00403     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00404     ScalarType zero = ScalarTraits<ScalarType>::zero();
00405     ScalarType one = ScalarTraits<ScalarType>::one();
00406     bool BadArgument = false;
00407 
00408     // Quick return if there is nothing to do!
00409     if( m == izero || n == izero || ( alpha == zero && beta == one ) ){ return; }
00410     
00411     // Otherwise, we need to check the argument list.
00412     if( m < izero ) { 
00413   cout << "BLAS::GEMV Error: M == " << m << endl;     
00414   BadArgument = true;
00415     }
00416     if( n < izero ) { 
00417   cout << "BLAS::GEMV Error: N == " << n << endl;     
00418   BadArgument = true;
00419     }
00420     if( lda < m ) { 
00421   cout << "BLAS::GEMV Error: LDA < MAX(1,M)"<< endl;      
00422   BadArgument = true;
00423     }
00424     if( incx == izero ) {
00425   cout << "BLAS::GEMV Error: INCX == 0"<< endl;
00426   BadArgument = true;
00427     }
00428     if( incy == izero ) {
00429   cout << "BLAS::GEMV Error: INCY == 0"<< endl;
00430   BadArgument = true;
00431     }
00432 
00433     if(!BadArgument) {
00434       OrdinalType i, j, lenx, leny, ix, iy, jx, jy; 
00435       OrdinalType kx = izero, ky = izero;
00436       ScalarType temp;
00437 
00438       // Determine the lengths of the vectors x and y.
00439       if(ETranspChar[trans] == 'N') {
00440   lenx = n;
00441   leny = m;
00442       } else {
00443   lenx = m;
00444   leny = n;
00445       }
00446 
00447       // Set the starting pointers for the vectors x and y if incx/y < 0.
00448       if (incx < izero ) { kx =  (ione - lenx)*incx; }
00449       if (incy < izero ) { ky =  (ione - leny)*incy; }
00450 
00451       // Form y = beta*y
00452       ix = kx; iy = ky;
00453       if(beta != one) {
00454   if (incy == ione) {
00455     if (beta == zero) {
00456       for(i = izero; i < leny; i++) { y[i] = zero; }
00457     } else {
00458       for(i = izero; i < leny; i++) { y[i] *= beta; }
00459     }
00460   } else {
00461     if (beta == zero) {
00462       for(i = izero; i < leny; i++) {
00463         y[iy] = zero;
00464         iy += incy;
00465       }
00466     } else {
00467       for(i = izero; i < leny; i++) {
00468         y[iy] *= beta;
00469         iy += incy;
00470       }
00471     }
00472   }
00473       }
00474   
00475       // Return if we don't have to do anything more.
00476       if(alpha == zero) { return; }
00477 
00478       if( ETranspChar[trans] == 'N' ) {
00479   // Form y = alpha*A*y
00480   jx = kx;
00481   if (incy == ione) {
00482     for(j = izero; j < n; j++) {
00483       if (x[jx] != zero) {
00484         temp = alpha*x[jx];
00485         for(i = izero; i < m; i++) {
00486     y[i] += temp*A[j*lda + i];
00487         }
00488       }
00489       jx += incx;
00490     }
00491   } else {
00492     for(j = izero; j < n; j++) {
00493       if (x[jx] != zero) {
00494         temp = alpha*x[jx];
00495         iy = ky;
00496         for(i = izero; i < m; i++) {
00497     y[iy] += temp*A[j*lda + i];
00498     iy += incy;
00499         }
00500       }
00501       jx += incx;
00502     }
00503   }
00504       } else {
00505   jy = ky;
00506   if (incx == ione) {
00507     for(j = izero; j < n; j++) {
00508       temp = zero;
00509       for(i = izero; i < m; i++) {
00510         temp += A[j*lda + i]*x[i];
00511       }
00512       y[jy] += alpha*temp;
00513       jy += incy;
00514     }
00515   } else {
00516     for(j = izero; j < n; j++) {
00517       temp = zero;
00518       ix = kx;
00519       for (i = izero; i < m; i++) {
00520         temp += A[j*lda + i]*x[ix];
00521         ix += incx;
00522       }
00523       y[jy] += alpha*temp;
00524       jy += incy;
00525     }
00526   }
00527       }
00528     } /* if (!BadArgument) */
00529   } /* end GEMV */
00530 
00531  template<typename OrdinalType, typename ScalarType>
00532  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
00533   {
00534     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00535     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00536     ScalarType zero = ScalarTraits<ScalarType>::zero();
00537     ScalarType one = ScalarTraits<ScalarType>::one();
00538     bool BadArgument = false;
00539 
00540     // Quick return if there is nothing to do!
00541     if( n == izero ){ return; }
00542     
00543     // Otherwise, we need to check the argument list.
00544     if( n < izero ) { 
00545       cout << "BLAS::TRMV Error: N == " << n << endl;     
00546       BadArgument = true;
00547     }
00548     if( lda < n ) { 
00549       cout << "BLAS::TRMV Error: LDA < MAX(1,N)"<< endl;      
00550       BadArgument = true;
00551     }
00552     if( incx == izero ) {
00553       cout << "BLAS::TRMV Error: INCX == 0"<< endl;
00554       BadArgument = true;
00555     }
00556 
00557     if(!BadArgument) {
00558       OrdinalType i, j, ix, jx, kx = izero;
00559       ScalarType temp;
00560       bool NoUnit = (EDiagChar[diag] == 'N');
00561 
00562       // Set the starting pointer for the vector x if incx < 0.
00563       if (incx < izero) { kx = (-n+ione)*incx; }
00564 
00565       // Start the operations for a nontransposed triangular matrix 
00566       if (ETranspChar[trans] == 'N') {
00567   /* Compute x = A*x */
00568   if (EUploChar[uplo] == 'U') {
00569     /* A is an upper triangular matrix */
00570     if (incx == ione) {
00571       for (j=izero; j<n; j++) {
00572         if (x[j] != zero) {
00573     temp = x[j];
00574     for (i=izero; i<j; i++) {
00575       x[i] += temp*A[j*lda + i];
00576     }
00577     if (NoUnit) 
00578       x[j] *= A[j*lda + j];
00579         }
00580       }
00581     } else {
00582       jx = kx;
00583       for (j=izero; j<n; j++) {
00584         if (x[jx] != zero) {
00585     temp = x[jx];
00586     ix = kx;
00587     for (i=izero; i<j; i++) {
00588       x[ix] += temp*A[j*lda + i];
00589       ix += incx;
00590     }
00591     if (NoUnit)
00592       x[jx] *= A[j*lda + j];
00593         }
00594         jx += incx;
00595       }
00596     } /* if (incx == ione) */
00597   } else { /* A is a lower triangular matrix */
00598     if (incx == ione) {
00599       for (j=n-ione; j>-ione; j--) {
00600         if (x[j] != zero) {
00601     temp = x[j];
00602     for (i=n-ione; i>j; i--) {
00603       x[i] += temp*A[j*lda + i];
00604     }
00605     if (NoUnit)
00606       x[j] *= A[j*lda + j];
00607         }
00608       }
00609     } else {
00610       kx += (n-ione)*incx;
00611       jx = kx;
00612       for (j=n-ione; j>-ione; j--) {
00613         if (x[jx] != zero) {
00614     temp = x[jx];
00615     ix = kx;
00616     for (i=n-ione; i>j; i--) {
00617       x[ix] += temp*A[j*lda + i];
00618       ix -= incx;
00619     }
00620     if (NoUnit) 
00621       x[jx] *= A[j*lda + j];
00622         }
00623         jx -= incx;
00624       }
00625     }
00626   } /* if (EUploChar[uplo]=='U') */
00627       } else { /* A is transposed/conjugated */
00628   /* Compute x = A'*x */
00629   if (EUploChar[uplo]=='U') {
00630     /* A is an upper triangular matrix */
00631     if (incx == ione) {
00632       for (j=n-ione; j>-ione; j--) {
00633         temp = x[j];
00634         if (NoUnit)
00635     temp *= A[j*lda + j];
00636         for (i=j-ione; i>-ione; i--) {
00637     temp += A[j*lda + i]*x[i];
00638         }
00639         x[j] = temp;
00640       }
00641     } else {
00642       jx = kx + (n-ione)*incx;
00643       for (j=n-ione; j>-ione; j--) {
00644         temp = x[jx];
00645         ix = jx;
00646         if (NoUnit)
00647     temp *= A[j*lda + j];
00648         for (i=j-ione; i>-ione; i--) {
00649     ix -= incx;
00650     temp += A[j*lda + i]*x[ix];
00651         }
00652         x[jx] = temp;
00653         jx -= incx;
00654       }
00655     }
00656   } else {
00657     /* A is a lower triangular matrix */
00658     if (incx == ione) {
00659       for (j=izero; j<n; j++) {
00660         temp = x[j];
00661         if (NoUnit)
00662     temp *= A[j*lda + j];
00663         for (i=j+ione; i<n; i++) {
00664     temp += A[j*lda + i]*x[i];
00665         }
00666         x[j] = temp;
00667       }
00668     } else {
00669       jx = kx;
00670       for (j=izero; j<n; j++) {
00671         temp = x[jx];
00672         ix = jx;
00673         if (NoUnit) 
00674     temp *= A[j*lda + j];
00675         for (i=j+ione; i<n; i++) {
00676     ix += incx;
00677     temp += A[j*lda + i]*x[ix];
00678         }
00679         x[jx] = temp;
00680         jx += incx;       
00681       }
00682     }
00683   } /* if (EUploChar[uplo]=='U') */
00684       } /* if (ETranspChar[trans]=='N') */
00685     } /* if (!BadArgument) */
00686   } /* end TRMV */
00687         
00688   template<typename OrdinalType, typename ScalarType>
00689   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
00690   {
00691     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00692     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00693     ScalarType zero = ScalarTraits<ScalarType>::zero();
00694     ScalarType one = ScalarTraits<ScalarType>::one();
00695     bool BadArgument = false;
00696 
00697     // Quick return if there is nothing to do!
00698     if( m == izero || n == izero || alpha == zero ){ return; }
00699     
00700     // Otherwise, we need to check the argument list.
00701     if( m < izero ) { 
00702   cout << "BLAS::GER Error: M == " << m << endl;      
00703   BadArgument = true;
00704     }
00705     if( n < izero ) { 
00706   cout << "BLAS::GER Error: N == " << n << endl;      
00707   BadArgument = true;
00708     }
00709     if( lda < m ) { 
00710   cout << "BLAS::GER Error: LDA < MAX(1,M)"<< endl;     
00711   BadArgument = true;
00712     }
00713     if( incx == 0 ) {
00714   cout << "BLAS::GER Error: INCX == 0"<< endl;
00715   BadArgument = true;
00716     }
00717     if( incy == 0 ) {
00718   cout << "BLAS::GER Error: INCY == 0"<< endl;
00719   BadArgument = true;
00720     }
00721 
00722     if(!BadArgument) {
00723       OrdinalType i, j, ix, jy = izero, kx = izero;
00724       ScalarType temp;
00725 
00726       // Set the starting pointers for the vectors x and y if incx/y < 0.
00727       if (incx < izero) { kx = (-m+ione)*incx; }
00728       if (incy < izero) { jy = (-n+ione)*incy; }
00729 
00730       // Start the operations for incx == 1
00731       if( incx == ione ) {
00732   for( j=izero; j<n; j++ ) {
00733     if ( y[jy] != zero ) {
00734       temp = alpha*y[jy];
00735       for ( i=izero; i<m; i++ ) {
00736         A[j*lda + i] += x[i]*temp;
00737       }
00738     }
00739     jy += incy;
00740   }
00741       } 
00742       else { // Start the operations for incx != 1
00743   for( j=izero; j<n; j++ ) {
00744     if ( y[jy] != zero ) {
00745       temp = alpha*y[jy];
00746       ix = kx;
00747       for( i=izero; i<m; i++ ) {
00748         A[j*lda + i] += x[ix]*temp;
00749         ix += incx;
00750       }
00751     }
00752     jy += incy;
00753   }
00754       }
00755     } /* if(!BadArgument) */
00756   } /* end GER */
00757   
00758 //------------------------------------------------------------------------------------------
00759 //      LEVEL 3 BLAS ROUTINES
00760 //------------------------------------------------------------------------------------------
00761         
00762   template<typename OrdinalType, typename ScalarType>
00763   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
00764   {
00765     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00766     ScalarType zero = ScalarTraits<ScalarType>::zero();
00767     ScalarType one = ScalarTraits<ScalarType>::one();
00768     OrdinalType i, j, p;
00769     OrdinalType NRowA = m, NRowB = k;
00770     ScalarType temp;
00771     bool BadArgument = false;
00772 
00773     // Change dimensions of matrix if either matrix is transposed
00774     if( !(ETranspChar[transa]=='N') ) {
00775       NRowA = k;
00776     }
00777     if( !(ETranspChar[transb]=='N') ) {
00778       NRowB = n;
00779     }
00780 
00781     // Quick return if there is nothing to do!
00782     if( (m==izero) || (n==izero) || (((alpha==zero)||(k==izero)) && (beta==one)) ){ return; }
00783     if( m < izero ) { 
00784       cout << "BLAS::GEMM Error: M == " << m << endl;     
00785       BadArgument = true;
00786     }
00787     if( n < izero ) { 
00788       cout << "BLAS::GEMM Error: N == " << n << endl;     
00789       BadArgument = true;
00790     }
00791     if( k < izero ) { 
00792       cout << "BLAS::GEMM Error: K == " << k << endl;     
00793       BadArgument = true;
00794     }
00795     if( lda < NRowA ) { 
00796       cout << "BLAS::GEMM Error: LDA < MAX(1,M)"<< endl;      
00797       BadArgument = true;
00798     }
00799     if( ldb < NRowB ) { 
00800       cout << "BLAS::GEMM Error: LDB < MAX(1,K)"<< endl;      
00801       BadArgument = true;
00802     }
00803      if( ldc < m ) { 
00804       cout << "BLAS::GEMM Error: LDC < MAX(1,M)"<< endl;      
00805       BadArgument = true;
00806     }
00807 
00808     if(!BadArgument) {
00809 
00810       // Only need to scale the resulting matrix C.
00811       if( alpha == zero ) {
00812   if( beta == zero ) {
00813     for (j=izero; j<n; j++) {
00814       for (i=izero; i<m; i++) {
00815         C[j*ldc + i] = zero;
00816       }
00817     }
00818   } else {
00819     for (j=izero; j<n; j++) {
00820       for (i=izero; i<m; i++) {
00821         C[j*ldc + i] *= beta;
00822       }
00823     }
00824   }
00825   return;
00826       }
00827       //
00828       // Now start the operations.
00829       //
00830       if ( ETranspChar[transb]=='N' ) {
00831   if ( ETranspChar[transa]=='N' ) {
00832     // Compute C = alpha*A*B + beta*C
00833     for (j=izero; j<n; j++) {
00834       if( beta == zero ) {
00835         for (i=izero; i<m; i++){
00836     C[j*ldc + i] = zero;
00837         }
00838       } else if( beta != one ) {
00839         for (i=izero; i<m; i++){
00840     C[j*ldc + i] *= beta;
00841         }
00842       }
00843       for (p=izero; p<k; p++){
00844         if (B[j*ldb + p] != zero ){
00845     temp = alpha*B[j*ldb + p];
00846     for (i=izero; i<m; i++) {
00847       C[j*ldc + i] += temp*A[p*lda + i];
00848     }
00849         }
00850       }
00851     }
00852   } else {
00853     // Compute C = alpha*A'*B + beta*C
00854     for (j=izero; j<n; j++) {
00855       for (i=izero; i<m; i++) {
00856         temp = zero;
00857         for (p=izero; p<k; p++) {
00858     temp += A[i*lda + p]*B[j*ldb + p];
00859         }
00860         if (beta == zero) {
00861     C[j*ldc + i] = alpha*temp;
00862         } else {
00863     C[j*ldc + i] = alpha*temp + beta*C[j*ldc + i];
00864         }
00865       }
00866     }
00867   }
00868       } else {
00869   if ( ETranspChar[transa]=='N' ) {
00870     // Compute C = alpha*A*B' + beta*C
00871     for (j=izero; j<n; j++) {
00872       if (beta == zero) {
00873         for (i=izero; i<m; i++) {
00874     C[j*ldc + i] = zero;
00875         } 
00876       } else if ( beta != one ) {
00877         for (i=izero; i<m; i++) {
00878     C[j*ldc + i] *= beta;
00879         }
00880       }
00881       for (p=izero; p<k; p++) {
00882         if (B[p*ldb + j] != zero) {
00883     temp = alpha*B[p*ldb + j];
00884     for (i=izero; i<m; i++) {
00885       C[j*ldc + i] += temp*A[p*lda + i];
00886     }
00887         }
00888       }
00889     }
00890   } else {
00891     // Compute C += alpha*A'*B' + beta*C
00892     for (j=izero; j<n; j++) {
00893       for (i=izero; i<m; i++) {
00894         temp = zero;
00895         for (p=izero; p<k; p++) {
00896     temp += A[i*lda + p]*B[p*ldb + j];
00897         }
00898         if (beta == zero) {
00899     C[j*ldc + i] = alpha*temp;
00900         } else {
00901     C[j*ldc + i] = alpha*temp + beta*C[j*ldc + i];
00902         }
00903       }
00904     }
00905   } // end if (ETranspChar[transa]=='N') ...
00906       } // end if (ETranspChar[transb]=='N') ...
00907     } // end if (!BadArgument) ...
00908   } // end of GEMM
00909 
00910 
00911   template<typename OrdinalType, typename ScalarType>
00912   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
00913   {
00914     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
00915     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
00916     ScalarType zero = ScalarTraits<ScalarType>::zero();
00917     ScalarType one = ScalarTraits<ScalarType>::one();
00918     OrdinalType i, j, k, NRowA = m;
00919     ScalarType temp1, temp2;
00920     bool BadArgument = false;
00921     bool Upper = (EUploChar[uplo] == 'U');
00922     if (ESideChar[side] == 'R') { NRowA = n; }
00923     
00924     // Quick return.
00925     if ( (m==izero) || (n==izero) || ( (alpha==zero)&&(beta==one) ) ) { return; }
00926     if( m < 0 ) { 
00927       cout << "BLAS::SYMM Error: M == "<< m << endl;
00928       BadArgument = true; }
00929     if( n < 0 ) {
00930       cout << "BLAS::SYMM Error: N == "<< n << endl;
00931       BadArgument = true; }
00932     if( lda < NRowA ) {
00933       cout << "BLAS::SYMM Error: LDA == "<<lda<<endl;
00934       BadArgument = true; }
00935     if( ldb < m ) {
00936       cout << "BLAS::SYMM Error: LDB == "<<ldb<<endl;
00937       BadArgument = true; }
00938     if( ldc < m ) {
00939       cout << "BLAS::SYMM Error: LDC == "<<ldc<<endl;
00940       BadArgument = true; }
00941 
00942     if(!BadArgument) {
00943 
00944       // Only need to scale C and return.
00945       if (alpha == zero) {
00946   if (beta == zero ) {
00947     for (j=izero; j<n; j++) {
00948       for (i=izero; i<m; i++) {
00949         C[j*ldc + i] = zero;
00950       }
00951     }
00952   } else {
00953     for (j=izero; j<n; j++) {
00954       for (i=izero; i<m; i++) {
00955         C[j*ldc + i] *= beta;
00956       }
00957     }
00958   }
00959   return;
00960       }
00961 
00962       if ( ESideChar[side] == 'L') {
00963   // Compute C = alpha*A*B + beta*C
00964 
00965   if (Upper) {
00966     // The symmetric part of A is stored in the upper triangular part of the matrix.
00967     for (j=izero; j<n; j++) {
00968       for (i=izero; i<m; i++) {
00969         temp1 = alpha*B[j*ldb + i];
00970         temp2 = zero;
00971         for (k=izero; k<i; k++) {
00972     C[j*ldc + k] += temp1*A[i*lda + k];
00973     temp2 += B[j*ldb + k]*A[i*lda + k];
00974         }
00975         if (beta == zero) {
00976     C[j*ldc + i] = temp1*A[i*lda + i] + alpha*temp2;
00977         } else {
00978     C[j*ldc + i] = beta*C[j*ldc + i] + temp1*A[i*lda + i] + alpha*temp2;
00979         }
00980       }
00981     }
00982   } else {
00983     // The symmetric part of A is stored in the lower triangular part of the matrix.
00984     for (j=izero; j<n; j++) {
00985       for (i=m-ione; i>-ione; i--) {
00986         temp1 = alpha*B[j*ldb + i];
00987         temp2 = zero;
00988         for (k=i+ione; k<m; k++) {
00989     C[j*ldc + k] += temp1*A[i*lda + k];
00990     temp2 += B[j*ldb + k]*A[i*lda + k];
00991         }
00992         if (beta == zero) {
00993     C[j*ldc + i] = temp1*A[i*lda + i] + alpha*temp2;
00994         } else {
00995     C[j*ldc + i] = beta*C[j*ldc + i] + temp1*A[i*lda + i] + alpha*temp2;
00996         }
00997       }
00998     }
00999   }
01000       } else {
01001   // Compute C = alpha*B*A + beta*C.
01002   for (j=izero; j<n; j++) {
01003     temp1 = alpha*A[j*lda + j];
01004     if (beta == zero) {
01005       for (i=izero; i<m; i++) {
01006         C[j*ldc + i] = temp1*B[j*ldb + i];
01007       }
01008     } else {
01009       for (i=izero; i<m; i++) {
01010         C[j*ldc + i] = beta*C[j*ldc + i] + temp1*B[j*ldb + i];
01011       }
01012     }
01013     for (k=izero; k<j; k++) {
01014       if (Upper) {
01015         temp1 = alpha*A[j*lda + k];
01016       } else {
01017         temp1 = alpha*A[k*lda + j];
01018       }
01019       for (i=izero; i<m; i++) {
01020         C[j*ldc + i] += temp1*B[k*ldb + i];
01021       }
01022     }
01023     for (k=j+ione; k<n; k++) {
01024       if (Upper) {
01025         temp1 = alpha*A[k*lda + j];
01026       } else {
01027         temp1 = alpha*A[j*lda + k];
01028       }
01029       for (i=izero; i<m; i++) {
01030         C[j*ldc + i] += temp1*B[k*ldb + i];
01031       }
01032     }
01033   }
01034       } // end if (ESideChar[side]=='L') ...
01035     } // end if(!BadArgument) ...
01036 } // end SYMM
01037   
01038   template<typename OrdinalType, typename ScalarType>
01039   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
01040   {
01041     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
01042     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
01043     ScalarType zero = ScalarTraits<ScalarType>::zero();
01044     ScalarType one = ScalarTraits<ScalarType>::one();
01045     OrdinalType i, j, k, NRowA = m;
01046     ScalarType temp;
01047     bool BadArgument = false;
01048     bool LSide = (ESideChar[side] == 'L');
01049     bool NoUnit = (EDiagChar[diag] == 'N');
01050     bool Upper = (EUploChar[uplo] == 'U');
01051 
01052     if(!LSide) { NRowA = n; }
01053 
01054     // Quick return.
01055     if (n==izero || m==izero) { return; }
01056     if( m < 0 ) {
01057       cout << "BLAS::TRMM Error: M == "<< m <<endl;
01058       BadArgument = true; }
01059     if( n < 0 ) {
01060       cout << "BLAS::TRMM Error: N == "<< n <<endl;
01061       BadArgument = true; }
01062     if( lda < NRowA ) {
01063       cout << "BLAS::TRMM Error: LDA == "<< lda << endl;
01064       BadArgument = true; }
01065     if( ldb < m ) {
01066       cout << "BLAS::TRMM Error: M == "<< ldb << endl;
01067       BadArgument = true; }
01068 
01069     if(!BadArgument) {
01070 
01071       // B only needs to be zeroed out.
01072       if( alpha == zero ) {
01073   for( j=izero; j<n; j++ ) {
01074     for( i=izero; i<m; i++ ) {
01075       B[j*ldb + i] = zero;
01076     }
01077   }
01078   return;
01079       }
01080       
01081       //  Start the computations. 
01082       if ( LSide ) {
01083   // A is on the left side of B.
01084   
01085   if ( ETranspChar[transa]=='N' ) {
01086     // Compute B = alpha*A*B
01087 
01088     if ( Upper ) {
01089       // A is upper triangular
01090       for( j=izero; j<n; j++ ) {
01091         for( k=izero; k<m; k++) {
01092     if ( B[j*ldb + k] != zero ) {
01093       temp = alpha*B[j*ldb + k];
01094       for( i=izero; i<k; i++ ) {
01095         B[j*ldb + i] += temp*A[k*lda + i];
01096       }
01097       if ( NoUnit )
01098         temp *=A[k*lda + k];
01099       B[j*ldb + k] = temp;
01100     }
01101         }
01102       }
01103     } else {
01104       // A is lower triangular
01105       for( j=izero; j<n; j++ ) {
01106         for( k=m-ione; k>-ione; k-- ) {
01107     if( B[j*ldb + k] != zero ) {
01108       temp = alpha*B[j*ldb + k];
01109       B[j*ldb + k] = temp;
01110       if ( NoUnit )
01111         B[j*ldb + k] *= A[k*lda + k];
01112       for( i=k+ione; i<m; i++ ) {
01113         B[j*ldb + i] += temp*A[k*lda + i];
01114       }
01115     }
01116         }
01117       }
01118     }
01119   } else {
01120     // Compute B = alpha*A'*B
01121     if( Upper ) {
01122       for( j=izero; j<n; j++ ) {
01123         for( i=m-ione; i>-ione; i-- ) {
01124     temp = B[j*ldb + i];
01125     if( NoUnit )
01126       temp *= A[i*lda + i];
01127     for( k=izero; k<i; k++ ) {
01128       temp += A[i*lda + k]*B[j*ldb + k];
01129     }
01130     B[j*ldb + i] = alpha*temp;
01131         }
01132       }
01133     } else {
01134       for( j=izero; j<n; j++ ) {
01135         for( i=izero; i<m; i++ ) {
01136     temp = B[j*ldb + i];
01137     if( NoUnit ) 
01138       temp *= A[i*lda + i];
01139     for( k=i+ione; k<m; k++ ) {
01140       temp += A[i*lda + k]*B[j*ldb + k];
01141     }
01142     B[j*ldb + i] = alpha*temp;
01143         }
01144       }
01145     }
01146   }
01147       } else {
01148   // A is on the right hand side of B.
01149   
01150   if( ETranspChar[transa] == 'N' ) {
01151     // Compute B = alpha*B*A
01152 
01153     if( Upper ) {
01154       // A is upper triangular.
01155       for( j=n-ione; j>-ione; j-- ) {
01156         temp = alpha;
01157         if( NoUnit )
01158     temp *= A[j*lda + j];
01159         for( i=izero; i<m; i++ ) {
01160     B[j*ldb + i] *= temp;
01161         }
01162         for( k=izero; k<j; k++ ) {
01163     if( A[j*lda + k] != zero ) {
01164       temp = alpha*A[j*lda + k];
01165       for( i=izero; i<m; i++ ) {
01166         B[j*ldb + i] += temp*B[k*ldb + i];
01167       }
01168     }
01169         }
01170       }
01171     } else {
01172       // A is lower triangular.
01173       for( j=izero; j<n; j++ ) {
01174         temp = alpha;
01175         if( NoUnit )
01176     temp *= A[j*lda + j];
01177         for( i=izero; i<m; i++ ) {
01178     B[j*ldb + i] *= temp;
01179         }
01180         for( k=j+ione; k<n; k++ ) {
01181     if( A[j*lda + k] != zero ) {
01182       temp = alpha*A[j*lda + k];
01183       for( i=izero; i<m; i++ ) {
01184         B[j*ldb + i] += temp*B[k*ldb + i];
01185       }
01186     }
01187         }
01188       }
01189     }
01190   } else {
01191     // Compute B = alpha*B*A'
01192 
01193     if( Upper ) {
01194       for( k=izero; k<n; k++ ) {
01195         for( j=izero; j<k; j++ ) {
01196     if( A[k*lda + j] != zero ) {
01197       temp = alpha*A[k*lda + j];
01198       for( i=izero; i<m; i++ ) {
01199         B[j*ldb + i] += temp*B[k*ldb + i];
01200       }
01201     }
01202         }
01203         temp = alpha;
01204         if( NoUnit ) 
01205     temp *= A[k*lda + k];
01206         if( temp != one ) {
01207     for( i=izero; i<m; i++ ) {
01208       B[k*ldb + i] *= temp;
01209     }
01210         }
01211       }
01212     } else {
01213       for( k=n-ione; k>-ione; k-- ) {
01214         for( j=k+ione; j<n; j++ ) {
01215     if( A[k*lda + j] != zero ) {
01216       temp = alpha*A[k*lda + j];
01217       for( i=izero; i<m; i++ ) {
01218         B[j*ldb + i] += temp*B[k*ldb + i];
01219       }
01220     }
01221         }
01222         temp = alpha;
01223         if( NoUnit )
01224     temp *= A[k*lda + k];
01225         if( temp != one ) {
01226     for( i=izero; i<m; i++ ) {
01227       B[k*ldb + i] *= temp;
01228     }
01229         }
01230       }
01231     }
01232   } // end if( ETranspChar[transa] == 'N' ) ...
01233       } // end if ( LSide ) ...
01234     } // end if (!BadArgument)
01235   } // end TRMM
01236   
01237   template<typename OrdinalType, typename ScalarType>
01238   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
01239   {
01240     OrdinalType izero = OrdinalTraits<OrdinalType>::zero();
01241     OrdinalType ione = OrdinalTraits<OrdinalType>::one();
01242     ScalarType zero = ScalarTraits<ScalarType>::zero();
01243     ScalarType one = ScalarTraits<ScalarType>::one();
01244     ScalarType temp;
01245     OrdinalType NRowA = m;
01246     bool BadArgument = false;
01247     bool NoUnit = (EDiagChar[diag]=='N');
01248     
01249     if (!(ESideChar[side] == 'L')) { NRowA = n; }
01250 
01251     // Quick return.
01252     if (n == izero || m == izero) { return; }
01253     if( m < izero ) {
01254       cout << "BLAS::TRSM Error: M == "<<m<<endl;
01255       BadArgument = true; }
01256     if( n < izero ) {
01257       cout << "BLAS::TRSM Error: N == "<<n<<endl;
01258       BadArgument = true; }
01259     if( lda < NRowA ) {
01260       cout << "BLAS::TRSM Error: LDA == "<<lda<<endl;
01261       BadArgument = true; }
01262     if( ldb < m ) {
01263       cout << "BLAS::TRSM Error: LDB == "<<ldb<<endl;
01264       BadArgument = true; }
01265 
01266     if(!BadArgument)
01267       {
01268   int i, j, k;
01269   // Set the solution to the zero vector.
01270   if(alpha == zero) {
01271       for(j = izero; j < n; j++) {
01272         for( i = izero; i < m; i++) {
01273         B[j*ldb+i] = zero;
01274           }
01275       }
01276   }
01277   else 
01278   { // Start the operations.
01279       if(ESideChar[side] == 'L') {
01280     //
01281         // Perform computations for OP(A)*X = alpha*B     
01282     //
01283     if(ETranspChar[transa] == 'N') {
01284         //
01285         //  Compute B = alpha*inv( A )*B
01286         //
01287         if(EUploChar[uplo] == 'U') { 
01288       // A is upper triangular.
01289       for(j = izero; j < n; j++) {
01290               // Perform alpha*B if alpha is not 1.
01291               if(alpha != one) {
01292                 for( i = izero; i < m; i++) {
01293                 B[j*ldb+i] *= alpha;
01294             }
01295           }
01296           // Perform a backsolve for column j of B.
01297           for(k = (m - ione); k > -ione; k--) {
01298         // If this entry is zero, we don't have to do anything.
01299         if (B[j*ldb + k] != zero) {
01300             if (NoUnit) {
01301           B[j*ldb + k] /= A[k*lda + k];
01302             }
01303             for(i = izero; i < k; i++) {
01304           B[j*ldb + i] -= B[j*ldb + k] * A[k*lda + i];
01305             }
01306         }
01307           }
01308       }
01309         }
01310         else 
01311         { // A is lower triangular.
01312                         for(j = izero; j < n; j++) {
01313                             // Perform alpha*B if alpha is not 1.
01314                             if(alpha != one) {
01315                                 for( i = izero; i < m; i++) {
01316                                     B[j*ldb+i] *= alpha;
01317                                 }
01318                             }
01319                             // Perform a forward solve for column j of B.
01320                             for(k = izero; k < m; k++) {
01321                                 // If this entry is zero, we don't have to do anything.
01322                                 if (B[j*ldb + k] != zero) {   
01323                                     if (NoUnit) {
01324                                         B[j*ldb + k] /= A[k*lda + k];
01325                                     }
01326                                     for(i = k+ione; i < m; i++) {
01327                                         B[j*ldb + i] -= B[j*ldb + k] * A[k*lda + i];
01328                                     }
01329                                 }
01330                             }
01331                         }
01332         } // end if (uplo == 'U')
01333     }  // if (transa =='N') 
01334         else { 
01335         //
01336         //  Compute B = alpha*inv( A' )*B
01337         //
01338         if(EUploChar[uplo] == 'U') { 
01339       // A is upper triangular.
01340       for(j = izero; j < n; j++) {
01341                   for( i = izero; i < m; i++) {
01342             temp = alpha*B[j*ldb+i];
01343             for(k = izero; k < i; k++) {
01344             temp -= A[i*lda + k] * B[j*ldb + k];
01345         }
01346         if (NoUnit) {
01347             temp /= A[i*lda + i];
01348         }
01349         B[j*ldb + i] = temp;
01350           }
01351       }
01352         }
01353         else
01354         { // A is lower triangular.
01355                         for(j = izero; j < n; j++) {
01356                             for(i = (m - ione) ; i > -ione; i--) {
01357                                 temp = alpha*B[j*ldb+i];
01358                               for(k = i+ione; k < m; k++) {
01359             temp -= A[i*lda + k] * B[j*ldb + k];
01360         }
01361         if (NoUnit) {
01362             temp /= A[i*lda + i];
01363         }
01364         B[j*ldb + i] = temp;
01365                             }
01366                         }
01367         }
01368     }
01369       }  // if (side == 'L')
01370       else { 
01371          // side == 'R'
01372          //
01373          // Perform computations for X*OP(A) = alpha*B      
01374          //
01375         if (ETranspChar[transa] == 'N') {
01376         //
01377         //  Compute B = alpha*B*inv( A )
01378         //
01379         if(EUploChar[uplo] == 'U') { 
01380       // A is upper triangular.
01381           // Perform a backsolve for column j of B.
01382       for(j = izero; j < n; j++) {
01383               // Perform alpha*B if alpha is not 1.
01384               if(alpha != one) {
01385                 for( i = izero; i < m; i++) {
01386                 B[j*ldb+i] *= alpha;
01387             }
01388           }
01389           for(k = izero; k < j; k++) {
01390         // If this entry is zero, we don't have to do anything.
01391         if (A[j*lda + k] != zero) {
01392             for(i = izero; i < m; i++) {
01393           B[j*ldb + i] -= A[j*lda + k] * B[k*ldb + i];
01394             }
01395         }
01396           }
01397           if (NoUnit) {
01398         temp = one/A[j*lda + j];
01399         for(i = izero; i < m; i++) {
01400             B[j*ldb + i] *= temp;
01401         }
01402           }
01403       }
01404         }
01405         else 
01406         { // A is lower triangular.
01407                         for(j = (n - ione); j > -ione; j--) {
01408                             // Perform alpha*B if alpha is not 1.
01409                             if(alpha != one) {
01410                                 for( i = izero; i < m; i++) {
01411                                     B[j*ldb+i] *= alpha;
01412                                 }
01413                             }
01414                             // Perform a forward solve for column j of B.
01415                             for(k = j+ione; k < n; k++) {
01416                                 // If this entry is zero, we don't have to do anything.
01417         if (A[j*lda + k] != zero) {
01418             for(i = izero; i < m; i++) {
01419                                         B[j*ldb + i] -= A[j*lda + k] * B[k*ldb + i]; 
01420                                     }
01421                                 } 
01422                             }
01423           if (NoUnit) {
01424         temp = one/A[j*lda + j];
01425         for(i = izero; i < m; i++) {
01426             B[j*ldb + i] *= temp;
01427         }
01428           }     
01429                         }
01430         } // end if (uplo == 'U')
01431     }  // if (transa =='N') 
01432         else { 
01433         //
01434         //  Compute B = alpha*B*inv( A' )
01435         //
01436         if(EUploChar[uplo] == 'U') { 
01437       // A is upper triangular.
01438       for(k = (n - ione); k > -ione; k--) {
01439           if (NoUnit) {
01440         temp = one/A[k*lda + k];
01441                     for(i = izero; i < m; i++) {
01442                 B[k*ldb + i] *= temp;
01443         }
01444           }
01445           for(j = izero; j < k; j++) {
01446         if (A[k*lda + j] != zero) {
01447             temp = A[k*lda + j];
01448             for(i = izero; i < m; i++) {
01449           B[j*ldb + i] -= temp*B[k*ldb + i];
01450             }
01451         }
01452           }
01453           if (alpha != one) {
01454         for (i = izero; i < m; i++) {
01455             B[k*ldb + i] *= alpha;
01456         }
01457           }
01458       }
01459         }
01460         else
01461         { // A is lower triangular.
01462       for(k = izero; k < n; k++) {
01463           if (NoUnit) {
01464         temp = one/A[k*lda + k];
01465         for (i = izero; i < m; i++) {
01466             B[k*ldb + i] *= temp;
01467         }
01468           }
01469           for(j = k+ione; j < n; j++) {
01470         if(A[k*lda + j] != zero) {
01471             temp = A[k*lda + j];
01472             for(i = izero; i < m; i++) {
01473           B[j*ldb + i] -= temp*B[k*ldb + i];
01474             }
01475         }
01476           }
01477           if (alpha != one) {
01478         for (i = izero; i < m; i++) {
01479             B[k*ldb + i] *= alpha;
01480         }
01481           }
01482                         }
01483         }
01484     }   
01485       }
01486   }
01487     }
01488   }
01489   
01490 #ifndef DOXYGEN_SHOULD_SKIP_THIS
01491 
01492   template<typename OrdinalType>
01493   class BLAS<OrdinalType, float>
01494   {    
01495   public:
01496     inline BLAS(void) {};
01497     inline BLAS(const BLAS<OrdinalType, float>& BLAS_source) {};
01498     inline virtual ~BLAS(void) {};
01499     void ROTG(float* da, float* db, float* c, float* s) const;
01500     void ROT(const OrdinalType n, float* dx, const OrdinalType incx, float* dy, const OrdinalType incy, float* c, float* s) const;
01501     float ASUM(const OrdinalType n, const float* x, const OrdinalType incx) const;
01502     void AXPY(const OrdinalType n, const float alpha, const float* x, const OrdinalType incx, float* y, const OrdinalType incy) const;
01503     void COPY(const OrdinalType n, const float* x, const OrdinalType incx, float* y, const OrdinalType incy) const;
01504     float DOT(const OrdinalType n, const float* x, const OrdinalType incx, const float* y, const OrdinalType incy) const;
01505     float NRM2(const OrdinalType n, const float* x, const OrdinalType incx) const;
01506     void SCAL(const OrdinalType n, const float alpha, float* x, const OrdinalType incx) const;
01507     OrdinalType IAMAX(const OrdinalType n, const float* x, const OrdinalType incx) const;
01508     void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, const float* x, const OrdinalType incx, const float beta, float* y, const OrdinalType incy) const;
01509     void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const float* A, const OrdinalType lda, float* x, const OrdinalType incx) const;
01510     void GER(const OrdinalType m, const OrdinalType n, const float alpha, const float* x, const OrdinalType incx, const float* y, const OrdinalType incy, float* A, const OrdinalType lda) const;
01511     void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const float alpha, const float* A, const OrdinalType lda, const float* B, const OrdinalType ldb, const float beta, float* C, const OrdinalType ldc) const;
01512     void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, const float *B, const OrdinalType ldb, const float beta, float *C, const OrdinalType ldc) const;
01513     void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, float* B, const OrdinalType ldb) const;
01514     void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, float* B, const OrdinalType ldb) const;
01515   };
01516 
01517   template<typename OrdinalType>
01518   void BLAS<OrdinalType, float>::ROTG(float* da, float* db, float* c, float* s) const
01519   { SROTG_F77(da, db, c, s ); }
01520 
01521   template<typename OrdinalType>
01522   void BLAS<OrdinalType, float>::ROT(const OrdinalType n, float* dx, const OrdinalType incx, float* dy, const OrdinalType incy, float* c, float* s) const
01523   { SROT_F77(&n, dx, &incx, dy, &incy, c, s); }
01524 
01525   template<typename OrdinalType>
01526   float BLAS<OrdinalType, float>::ASUM(const OrdinalType n, const float* x, const OrdinalType incx) const
01527   { return SASUM_F77(&n, x, &incx); }
01528   
01529   template<typename OrdinalType>
01530   void BLAS<OrdinalType, float>::AXPY(const OrdinalType n, const float alpha, const float* x, const OrdinalType incx, float* y, const OrdinalType incy) const
01531   { SAXPY_F77(&n, &alpha, x, &incx, y, &incy); }
01532   
01533   template<typename OrdinalType>
01534   void BLAS<OrdinalType, float>::COPY(const OrdinalType n, const float* x, const OrdinalType incx, float* y, const OrdinalType incy) const 
01535   { SCOPY_F77(&n, x, &incx, y, &incy); }
01536   
01537   template<typename OrdinalType>
01538   float BLAS<OrdinalType, float>::DOT(const OrdinalType n, const float* x, const OrdinalType incx, const float* y, const OrdinalType incy) const
01539   { return SDOT_F77(&n, x, &incx, y, &incy); }
01540   
01541   template<typename OrdinalType>
01542   OrdinalType BLAS<OrdinalType, float>::IAMAX(const OrdinalType n, const float* x, const OrdinalType incx) const
01543   { return ISAMAX_F77(&n, x, &incx); }
01544 
01545   template<typename OrdinalType>
01546   float BLAS<OrdinalType, float>::NRM2(const OrdinalType n, const float* x, const OrdinalType incx) const
01547   { return SNRM2_F77(&n, x, &incx); }
01548   
01549   template<typename OrdinalType>
01550   void BLAS<OrdinalType, float>::SCAL(const OrdinalType n, const float alpha, float* x, const OrdinalType incx) const
01551   { SSCAL_F77(&n, &alpha, x, &incx); }
01552   
01553   template<typename OrdinalType>
01554   void BLAS<OrdinalType, float>::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, const float* x, const OrdinalType incx, const float beta, float* y, const OrdinalType incy) const
01555   { SGEMV_F77(CHAR_MACRO(ETranspChar[trans]), &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }
01556   
01557   template<typename OrdinalType>
01558   void BLAS<OrdinalType, float>::GER(const OrdinalType m, const OrdinalType n, const float alpha, const float* x, const OrdinalType incx, const float* y, const OrdinalType incy, float* A, const OrdinalType lda) const
01559   { SGER_F77(&m, &n, &alpha, x, &incx, y, &incy, A, &lda); }
01560 
01561   template<typename OrdinalType>
01562   void BLAS<OrdinalType, float>::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const float* A, const OrdinalType lda, float* x, const OrdinalType incx) const
01563   { STRMV_F77(CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[trans]), CHAR_MACRO(EDiagChar[diag]), &n, A, &lda, x, &incx); }
01564   
01565   template<typename OrdinalType>
01566   void BLAS<OrdinalType, float>::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const float alpha, const float* A, const OrdinalType lda, const float* B, const OrdinalType ldb, const float beta, float* C, const OrdinalType ldc) const
01567   { SGEMM_F77(CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(ETranspChar[transb]), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01568   
01569   template<typename OrdinalType>
01570   void BLAS<OrdinalType, float>::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, const float* B, const OrdinalType ldb, const float beta, float* C, const OrdinalType ldc) const
01571   { SSYMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), &m, &n, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01572   
01573   template<typename OrdinalType>
01574   void BLAS<OrdinalType, float>::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, float* B, const OrdinalType ldb) const
01575   { STRMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01576   
01577   template<typename OrdinalType>
01578   void BLAS<OrdinalType, float>::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const float alpha, const float* A, const OrdinalType lda, float* B, const OrdinalType ldb) const
01579   { STRSM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01580 
01581 
01582   template<typename OrdinalType>
01583   class BLAS<OrdinalType, double>
01584   {    
01585   public:
01586     inline BLAS(void) {};
01587     inline BLAS(const BLAS<OrdinalType, double>& BLAS_source) {};
01588     inline virtual ~BLAS(void) {};
01589     void ROTG(double* da, double* db, double* c, double* s) const;
01590     void ROT(const OrdinalType n, double* dx, const OrdinalType incx, double* dy, const OrdinalType incy, double* c, double* s) const;
01591     double ASUM(const OrdinalType n, const double* x, const OrdinalType incx) const;
01592     void AXPY(const OrdinalType n, const double alpha, const double* x, const OrdinalType incx, double* y, const OrdinalType incy) const;
01593     void COPY(const OrdinalType n, const double* x, const OrdinalType incx, double* y, const OrdinalType incy) const;
01594     double DOT(const OrdinalType n, const double* x, const OrdinalType incx, const double* y, const OrdinalType incy) const;
01595     double NRM2(const OrdinalType n, const double* x, const OrdinalType incx) const;
01596     void SCAL(const OrdinalType n, const double alpha, double* x, const OrdinalType incx) const;
01597     OrdinalType IAMAX(const OrdinalType n, const double* x, const OrdinalType incx) const;
01598     void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, const double* x, const OrdinalType incx, const double beta, double* y, const OrdinalType incy) const;
01599     void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const double* A, const OrdinalType lda, double* x, const OrdinalType incx) const;
01600     void GER(const OrdinalType m, const OrdinalType n, const double alpha, const double* x, const OrdinalType incx, const double* y, const OrdinalType incy, double* A, const OrdinalType lda) const;
01601     void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const double alpha, const double* A, const OrdinalType lda, const double* B, const OrdinalType ldb, const double beta, double* C, const OrdinalType ldc) const;
01602     void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, const double *B, const OrdinalType ldb, const double beta, double *C, const OrdinalType ldc) const;
01603     void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, double* B, const OrdinalType ldb) const;
01604     void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, double* B, const OrdinalType ldb) const;
01605   };
01606   
01607   template<typename OrdinalType>
01608   void BLAS<OrdinalType, double>::ROTG(double* da, double* db, double* c, double* s) const
01609   { DROTG_F77(da, db, c, s); }
01610 
01611   template<typename OrdinalType>
01612   void BLAS<OrdinalType, double>::ROT(const OrdinalType n, double* dx, const OrdinalType incx, double* dy, const OrdinalType incy, double* c, double* s) const
01613   { DROT_F77(&n, dx, &incx, dy, &incy, c, s); }
01614 
01615   template<typename OrdinalType>
01616   double BLAS<OrdinalType, double>::ASUM(const OrdinalType n, const double* x, const OrdinalType incx) const
01617   { return DASUM_F77(&n, x, &incx); }
01618   
01619   template<typename OrdinalType>
01620   void BLAS<OrdinalType, double>::AXPY(const OrdinalType n, const double alpha, const double* x, const OrdinalType incx, double* y, const OrdinalType incy) const
01621   { DAXPY_F77(&n, &alpha, x, &incx, y, &incy); }
01622   
01623   template<typename OrdinalType>
01624   void BLAS<OrdinalType, double>::COPY(const OrdinalType n, const double* x, const OrdinalType incx, double* y, const OrdinalType incy) const
01625   { DCOPY_F77(&n, x, &incx, y, &incy); }
01626   
01627   template<typename OrdinalType>
01628   double BLAS<OrdinalType, double>::DOT(const OrdinalType n, const double* x, const OrdinalType incx, const double* y, const OrdinalType incy) const
01629   { return DDOT_F77(&n, x, &incx, y, &incy); }
01630   
01631   template<typename OrdinalType>
01632   OrdinalType BLAS<OrdinalType, double>::IAMAX(const OrdinalType n, const double* x, const OrdinalType incx) const
01633   { return IDAMAX_F77(&n, x, &incx); }
01634 
01635   template<typename OrdinalType>
01636   double BLAS<OrdinalType, double>::NRM2(const OrdinalType n, const double* x, const OrdinalType incx) const
01637   { return DNRM2_F77(&n, x, &incx); }
01638   
01639   template<typename OrdinalType>
01640   void BLAS<OrdinalType, double>::SCAL(const OrdinalType n, const double alpha, double* x, const OrdinalType incx) const
01641   { DSCAL_F77(&n, &alpha, x, &incx); }
01642   
01643   template<typename OrdinalType>
01644   void BLAS<OrdinalType, double>::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, const double* x, const OrdinalType incx, const double beta, double* y, const OrdinalType incy) const
01645   { DGEMV_F77(CHAR_MACRO(ETranspChar[trans]), &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }
01646   
01647   template<typename OrdinalType>
01648   void BLAS<OrdinalType, double>::GER(const OrdinalType m, const OrdinalType n, const double alpha, const double* x, const OrdinalType incx, const double* y, const OrdinalType incy, double* A, const OrdinalType lda) const
01649   { DGER_F77(&m, &n, &alpha, x, &incx, y, &incy, A, &lda); }
01650 
01651   template<typename OrdinalType>
01652   void BLAS<OrdinalType, double>::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const double* A, const OrdinalType lda, double* x, const OrdinalType incx) const
01653   { DTRMV_F77(CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[trans]), CHAR_MACRO(EDiagChar[diag]), &n, A, &lda, x, &incx); }
01654   
01655   template<typename OrdinalType>
01656   void BLAS<OrdinalType, double>::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const double alpha, const double* A, const OrdinalType lda, const double* B, const OrdinalType ldb, const double beta, double* C, const OrdinalType ldc) const
01657   { DGEMM_F77(CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(ETranspChar[transb]), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01658   
01659   template<typename OrdinalType>
01660   void BLAS<OrdinalType, double>::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, const double *B, const OrdinalType ldb, const double beta, double *C, const OrdinalType ldc) const
01661   { DSYMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), &m, &n, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01662   
01663   template<typename OrdinalType>
01664   void BLAS<OrdinalType, double>::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, double* B, const OrdinalType ldb) const
01665   { DTRMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01666 
01667   template<typename OrdinalType>
01668   void BLAS<OrdinalType, double>::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const double alpha, const double* A, const OrdinalType lda, double* B, const OrdinalType ldb) const
01669   { DTRSM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01670   
01671 #ifdef HAVE_TEUCHOS_COMPLEX
01672 
01673   template<typename OrdinalType>
01674   class BLAS<OrdinalType, complex<float> >
01675   {    
01676   public:
01677     inline BLAS(void) {};
01678     inline BLAS(const BLAS<OrdinalType, complex<float> >& BLAS_source) {};
01679     inline virtual ~BLAS(void) {};
01680     void ROTG(complex<float>* da, complex<float>* db, float* c, complex<float>* s) const;
01681     void ROT(const OrdinalType n, complex<float>* dx, const OrdinalType incx, complex<float>* dy, const OrdinalType incy, float* c, complex<float>* s) const;
01682     float ASUM(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const;
01683     void AXPY(const OrdinalType n, const complex<float> alpha, const complex<float>* x, const OrdinalType incx, complex<float>* y, const OrdinalType incy) const;
01684     void COPY(const OrdinalType n, const complex<float>* x, const OrdinalType incx, complex<float>* y, const OrdinalType incy) const;
01685     complex<float> DOT(const OrdinalType n, const complex<float>* x, const OrdinalType incx, const complex<float>* y, const OrdinalType incy) const;
01686     float NRM2(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const;
01687     void SCAL(const OrdinalType n, const complex<float> alpha, complex<float>* x, const OrdinalType incx) const;
01688     OrdinalType IAMAX(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const;
01689     void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float>* x, const OrdinalType incx, const complex<float> beta, complex<float>* y, const OrdinalType incy) const;
01690     void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const complex<float>* A, const OrdinalType lda, complex<float>* x, const OrdinalType incx) const;
01691     void GER(const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* x, const OrdinalType incx, const complex<float>* y, const OrdinalType incy, complex<float>* A, const OrdinalType lda) const;
01692     void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float>* B, const OrdinalType ldb, const complex<float> beta, complex<float>* C, const OrdinalType ldc) const;
01693     void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float> *B, const OrdinalType ldb, const complex<float> beta, complex<float> *C, const OrdinalType ldc) const;
01694     void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, complex<float>* B, const OrdinalType ldb) const;
01695     void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, complex<float>* B, const OrdinalType ldb) const;
01696   };
01697 
01698   template<typename OrdinalType>
01699   void BLAS<OrdinalType, complex<float> >::ROTG(complex<float>* da, complex<float>* db, float* c, complex<float>* s) const
01700   { CROTG_F77(da, db, c, s ); }
01701 
01702   template<typename OrdinalType>
01703   void BLAS<OrdinalType, complex<float> >::ROT(const OrdinalType n, complex<float>* dx, const OrdinalType incx, complex<float>* dy, const OrdinalType incy, float* c, complex<float>* s) const
01704   { CROT_F77(&n, dx, &incx, dy, &incy, c, s); }
01705 
01706   template<typename OrdinalType>
01707   float BLAS<OrdinalType, complex<float> >::ASUM(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const
01708   { return CASUM_F77(&n, x, &incx); }
01709   
01710   template<typename OrdinalType>
01711   void BLAS<OrdinalType, complex<float> >::AXPY(const OrdinalType n, const complex<float> alpha, const complex<float>* x, const OrdinalType incx, complex<float>* y, const OrdinalType incy) const
01712   { CAXPY_F77(&n, &alpha, x, &incx, y, &incy); }
01713   
01714   template<typename OrdinalType>
01715   void BLAS<OrdinalType, complex<float> >::COPY(const OrdinalType n, const complex<float>* x, const OrdinalType incx, complex<float>* y, const OrdinalType incy) const
01716   { CCOPY_F77(&n, x, &incx, y, &incy); }
01717   
01718   template<typename OrdinalType>
01719   complex<float> BLAS<OrdinalType, complex<float> >::DOT(const OrdinalType n, const complex<float>* x, const OrdinalType incx, const complex<float>* y, const OrdinalType incy) const
01720   { return CDOT_F77(&n, x, &incx, y, &incy); }
01721   
01722   template<typename OrdinalType>
01723   OrdinalType BLAS<OrdinalType, complex<float> >::IAMAX(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const
01724   { return ICAMAX_F77(&n, x, &incx); }
01725 
01726   template<typename OrdinalType>
01727   float BLAS<OrdinalType, complex<float> >::NRM2(const OrdinalType n, const complex<float>* x, const OrdinalType incx) const
01728   { return CNRM2_F77(&n, x, &incx); }
01729   
01730   template<typename OrdinalType>
01731   void BLAS<OrdinalType, complex<float> >::SCAL(const OrdinalType n, const complex<float> alpha, complex<float>* x, const OrdinalType incx) const
01732   { CSCAL_F77(&n, &alpha, x, &incx); }
01733   
01734   template<typename OrdinalType>
01735   void BLAS<OrdinalType, complex<float> >::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float>* x, const OrdinalType incx, const complex<float> beta, complex<float>* y, const OrdinalType incy) const
01736   { CGEMV_F77(CHAR_MACRO(ETranspChar[trans]), &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }
01737   
01738   template<typename OrdinalType>
01739   void BLAS<OrdinalType, complex<float> >::GER(const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* x, const OrdinalType incx, const complex<float>* y, const OrdinalType incy, complex<float>* A, const OrdinalType lda) const
01740   { CGER_F77(&m, &n, &alpha, x, &incx, y, &incy, A, &lda); }
01741 
01742   template<typename OrdinalType>
01743   void BLAS<OrdinalType, complex<float> >::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const complex<float>* A, const OrdinalType lda, complex<float>* x, const OrdinalType incx) const
01744   { CTRMV_F77(CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[trans]), CHAR_MACRO(EDiagChar[diag]), &n, A, &lda, x, &incx); }
01745   
01746   template<typename OrdinalType>
01747   void BLAS<OrdinalType, complex<float> >::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float>* B, const OrdinalType ldb, const complex<float> beta, complex<float>* C, const OrdinalType ldc) const
01748   { CGEMM_F77(CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(ETranspChar[transb]), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); } 
01749  
01750   template<typename OrdinalType>
01751   void BLAS<OrdinalType, complex<float> >::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, const complex<float>* B, const OrdinalType ldb, const complex<float> beta, complex<float>* C, const OrdinalType ldc) const
01752   { CSYMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), &m, &n, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01753   
01754   template<typename OrdinalType>
01755   void BLAS<OrdinalType, complex<float> >::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, complex<float>* B, const OrdinalType ldb) const
01756   { CTRMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01757   
01758   template<typename OrdinalType>
01759   void BLAS<OrdinalType, complex<float> >::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<float> alpha, const complex<float>* A, const OrdinalType lda, complex<float>* B, const OrdinalType ldb) const
01760   { CTRSM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01761 
01762 
01763   template<typename OrdinalType>
01764   class BLAS<OrdinalType, complex<double> >
01765   {    
01766   public:
01767     inline BLAS(void) {};
01768     inline BLAS(const BLAS<OrdinalType, complex<double> >& BLAS_source) {};
01769     inline virtual ~BLAS(void) {};
01770     void ROTG(complex<double>* da, complex<double>* db, double* c, complex<double>* s) const;
01771     void ROT(const OrdinalType n, complex<double>* dx, const OrdinalType incx, complex<double>* dy, const OrdinalType incy, double* c, complex<double>* s) const;
01772     double ASUM(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const;
01773     void AXPY(const OrdinalType n, const complex<double> alpha, const complex<double>* x, const OrdinalType incx, complex<double>* y, const OrdinalType incy) const;
01774     void COPY(const OrdinalType n, const complex<double>* x, const OrdinalType incx, complex<double>* y, const OrdinalType incy) const;
01775     complex<double> DOT(const OrdinalType n, const complex<double>* x, const OrdinalType incx, const complex<double>* y, const OrdinalType incy) const;
01776     double NRM2(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const;
01777     void SCAL(const OrdinalType n, const complex<double> alpha, complex<double>* x, const OrdinalType incx) const;
01778     OrdinalType IAMAX(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const;
01779     void GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double>* x, const OrdinalType incx, const complex<double> beta, complex<double>* y, const OrdinalType incy) const;
01780     void TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const complex<double>* A, const OrdinalType lda, complex<double>* x, const OrdinalType incx) const;
01781     void GER(const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* x, const OrdinalType incx, const complex<double>* y, const OrdinalType incy, complex<double>* A, const OrdinalType lda) const;
01782     void GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double>* B, const OrdinalType ldb, const complex<double> beta, complex<double>* C, const OrdinalType ldc) const;
01783     void SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double> *B, const OrdinalType ldb, const complex<double> beta, complex<double> *C, const OrdinalType ldc) const;
01784     void TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, complex<double>* B, const OrdinalType ldb) const;
01785     void TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, complex<double>* B, const OrdinalType ldb) const;
01786   };
01787   
01788   template<typename OrdinalType>
01789   void BLAS<OrdinalType, complex<double> >::ROTG(complex<double>* da, complex<double>* db, double* c, complex<double>* s) const
01790   { ZROTG_F77(da, db, c, s); }
01791 
01792   template<typename OrdinalType>
01793   void BLAS<OrdinalType, complex<double> >::ROT(const OrdinalType n, complex<double>* dx, const OrdinalType incx, complex<double>* dy, const OrdinalType incy, double* c, complex<double>* s) const
01794   { ZROT_F77(&n, dx, &incx, dy, &incy, c, s); }
01795 
01796   template<typename OrdinalType>
01797   double BLAS<OrdinalType, complex<double> >::ASUM(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const
01798   { return ZASUM_F77(&n, x, &incx); }
01799   
01800   template<typename OrdinalType>
01801   void BLAS<OrdinalType, complex<double> >::AXPY(const OrdinalType n, const complex<double> alpha, const complex<double>* x, const OrdinalType incx, complex<double>* y, const OrdinalType incy) const
01802   { ZAXPY_F77(&n, &alpha, x, &incx, y, &incy); }
01803   
01804   template<typename OrdinalType>
01805   void BLAS<OrdinalType, complex<double> >::COPY(const OrdinalType n, const complex<double>* x, const OrdinalType incx, complex<double>* y, const OrdinalType incy) const
01806   { ZCOPY_F77(&n, x, &incx, y, &incy); }
01807   
01808   template<typename OrdinalType>
01809   complex<double> BLAS<OrdinalType, complex<double> >::DOT(const OrdinalType n, const complex<double>* x, const OrdinalType incx, const complex<double>* y, const OrdinalType incy) const
01810   { return ZDOT_F77(&n, x, &incx, y, &incy); }
01811   
01812   template<typename OrdinalType>
01813   OrdinalType BLAS<OrdinalType, complex<double> >::IAMAX(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const
01814   { return IZAMAX_F77(&n, x, &incx); }
01815 
01816   template<typename OrdinalType>
01817   double BLAS<OrdinalType, complex<double> >::NRM2(const OrdinalType n, const complex<double>* x, const OrdinalType incx) const
01818   { return ZNRM2_F77(&n, x, &incx); }
01819   
01820   template<typename OrdinalType>
01821   void BLAS<OrdinalType, complex<double> >::SCAL(const OrdinalType n, const complex<double> alpha, complex<double>* x, const OrdinalType incx) const
01822   { ZSCAL_F77(&n, &alpha, x, &incx); }
01823   
01824   template<typename OrdinalType>
01825   void BLAS<OrdinalType, complex<double> >::GEMV(ETransp trans, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double>* x, const OrdinalType incx, const complex<double> beta, complex<double>* y, const OrdinalType incy) const
01826   { ZGEMV_F77(CHAR_MACRO(ETranspChar[trans]), &m, &n, &alpha, A, &lda, x, &incx, &beta, y, &incy); }
01827   
01828   template<typename OrdinalType>
01829   void BLAS<OrdinalType, complex<double> >::GER(const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* x, const OrdinalType incx, const complex<double>* y, const OrdinalType incy, complex<double>* A, const OrdinalType lda) const
01830   { ZGER_F77(&m, &n, &alpha, x, &incx, y, &incy, A, &lda); }
01831 
01832   template<typename OrdinalType>
01833   void BLAS<OrdinalType, complex<double> >::TRMV(EUplo uplo, ETransp trans, EDiag diag, const OrdinalType n, const complex<double>* A, const OrdinalType lda, complex<double>* x, const OrdinalType incx) const
01834   { ZTRMV_F77(CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[trans]), CHAR_MACRO(EDiagChar[diag]), &n, A, &lda, x, &incx); }
01835   
01836   template<typename OrdinalType>
01837   void BLAS<OrdinalType, complex<double> >::GEMM(ETransp transa, ETransp transb, const OrdinalType m, const OrdinalType n, const OrdinalType k, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double>* B, const OrdinalType ldb, const complex<double> beta, complex<double>* C, const OrdinalType ldc) const
01838   { ZGEMM_F77(CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(ETranspChar[transb]), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01839   
01840   template<typename OrdinalType>
01841   void BLAS<OrdinalType, complex<double> >::SYMM(ESide side, EUplo uplo, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, const complex<double> *B, const OrdinalType ldb, const complex<double> beta, complex<double> *C, const OrdinalType ldc) const
01842   { ZSYMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), &m, &n, &alpha, A, &lda, B, &ldb, &beta, C, &ldc); }
01843   
01844   template<typename OrdinalType>
01845   void BLAS<OrdinalType, complex<double> >::TRMM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, complex<double>* B, const OrdinalType ldb) const
01846   { ZTRMM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01847 
01848   template<typename OrdinalType>
01849   void BLAS<OrdinalType, complex<double> >::TRSM(ESide side, EUplo uplo, ETransp transa, EDiag diag, const OrdinalType m, const OrdinalType n, const complex<double> alpha, const complex<double>* A, const OrdinalType lda, complex<double>* B, const OrdinalType ldb) const
01850   { ZTRSM_F77(CHAR_MACRO(ESideChar[side]), CHAR_MACRO(EUploChar[uplo]), CHAR_MACRO(ETranspChar[transa]), CHAR_MACRO(EDiagChar[diag]), &m, &n, &alpha, A, &lda, B, &ldb); }
01851   
01852 #endif // HAVE_TEUCHOS_COMPLEX
01853 
01854 #endif // DOXYGEN_SHOULD_SKIP_THIS
01855 
01856 } // namespace Teuchos
01857 
01858 #endif // _TEUCHOS_BLAS_HPP_
Generated on Thu Sep 18 12:30:29 2008 for Teuchos - Trilinos Tools Package by doxygen 1.3.9.1