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