computation/interfaces/Blas_1.cxx

00001 // Copyright (C) 2001-2009 Vivien Mallet
00002 //
00003 // This file is part of the linear-algebra library Seldon,
00004 // http://seldon.sourceforge.net/.
00005 //
00006 // Seldon is free software; you can redistribute it and/or modify it under the
00007 // terms of the GNU Lesser General Public License as published by the Free
00008 // Software Foundation; either version 2.1 of the License, or (at your option)
00009 // any later version.
00010 //
00011 // Seldon is distributed in the hope that it will be useful, but WITHOUT ANY
00012 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00013 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
00014 // more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public License
00017 // along with Seldon. If not, see http://www.gnu.org/licenses/.
00018 
00019 
00020 #ifndef SELDON_FILE_BLAS_1_CXX
00021 
00022 /*
00023   Functions included in this file:
00024 
00025   xROTG   (GenRot)
00026   xROTMG  (GenModifRot)
00027   xROT    (ApplyRot)
00028   xROTM   (ApplyModifRot)
00029   xSWAP   (Swap)
00030   xSCAL   (Mlt)
00031   xCOPY   (Copy)
00032   xAXPY   (Add)
00033   xDOT    (DotProd)
00034   xDOTU   (DotProd)
00035   SDSDOT  (ScaledDotProd)
00036   xDOTC   (DotProdConj)
00037   xASUM   (Norm1)
00038   xNRM2   (Norm2)
00039   IxAMAX  (GetMaxAbsIndex)
00040 */
00041 
00042 extern "C"
00043 {
00044 #include "cblas.h"
00045 }
00046 
00047 namespace Seldon
00048 {
00049 
00050 
00052   // GenRot //
00053 
00054 
00055   void GenRot(float& a, float& b, float& c, float& d)
00056   {
00057     cblas_srotg(&a, &b, &c, &d);
00058   }
00059 
00060 
00061   void GenRot(double& a, double& b, double& c, double& d)
00062   {
00063     cblas_drotg(&a, &b, &c, &d);
00064   }
00065 
00066 
00067   // GenRot //
00069 
00070 
00071 
00073   // GenModifRot //
00074 
00075 
00076   void GenModifRot(float& d1, float& d2,
00077                    float& x1, const float& y1,
00078                    float* param)
00079   {
00080     cblas_srotmg(&d1, &d2, &x1, y1, param);
00081   }
00082 
00083 
00084   void GenModifRot(double& d1, double& d2,
00085                    double& x1, const double& y1,
00086                    double* param)
00087   {
00088     cblas_drotmg(&d1, &d2, &x1, y1, param);
00089   }
00090 
00091 
00092   // GenModifRot //
00094 
00095 
00096 
00098   // ApplyRot //
00099 
00100 
00101   template <class Allocator>
00102   void ApplyRot(Vector<float, VectFull, Allocator>& X,
00103                 Vector<float, VectFull, Allocator>& Y,
00104                 const float c, const float s)
00105   {
00106     cblas_srot(X.GetLength(), X.GetData(), 1,
00107                Y.GetData(), 1, c, s);
00108   }
00109 
00110 
00111   template <class Allocator>
00112   void ApplyRot(Vector<double, VectFull, Allocator>& X,
00113                 Vector<double, VectFull, Allocator>& Y,
00114                 const double c, const double s)
00115   {
00116     cblas_drot(X.GetLength(), X.GetData(), 1,
00117                Y.GetData(), 1, c, s);
00118   }
00119 
00120 
00121   // ApplyRot //
00123 
00124 
00125 
00127   // ApplyModifRot //
00128 
00129 
00130   template <class Allocator>
00131   void ApplyModifRot(Vector<float, VectFull, Allocator>& X,
00132                      Vector<float, VectFull, Allocator>& Y,
00133                      const float* param)
00134   {
00135     cblas_srotm(X.GetLength(), X.GetData(), 1,
00136                 Y.GetData(), 1, param);
00137   }
00138 
00139 
00140   template <class Allocator>
00141   void ApplyModifRot(Vector<double, VectFull, Allocator>& X,
00142                      Vector<double, VectFull, Allocator>& Y,
00143                      const double* param)
00144   {
00145     cblas_drotm(X.GetLength(), X.GetData(), 1,
00146                 Y.GetData(), 1, param);
00147   }
00148 
00149 
00150   // ApplyModifRot //
00152 
00153 
00154 
00156   // Swap //
00157 
00158 
00159   template <class Allocator>
00160   void Swap(Vector<float, VectFull, Allocator>& X,
00161             Vector<float, VectFull, Allocator>& Y)
00162   {
00163 
00164 #ifdef SELDON_CHECK_DIMENSIONS
00165     CheckDim(X, Y, "Swap(X, Y)", "X <-> Y");
00166 #endif
00167 
00168     cblas_sswap(X.GetLength(), X.GetData(), 1,
00169                 Y.GetData(), 1);
00170   }
00171 
00172 
00173   template <class Allocator>
00174   void Swap(Vector<double, VectFull, Allocator>& X,
00175             Vector<double, VectFull, Allocator>& Y)
00176   {
00177 
00178 #ifdef SELDON_CHECK_DIMENSIONS
00179     CheckDim(X, Y, "Swap(X, Y)", "X <-> Y");
00180 #endif
00181 
00182     cblas_dswap(X.GetLength(), X.GetData(), 1,
00183                 Y.GetData(), 1);
00184   }
00185 
00186 
00187   template <class Allocator>
00188   void Swap(Vector<complex<float>, VectFull, Allocator>& X,
00189             Vector<complex<float>, VectFull, Allocator>& Y)
00190   {
00191 
00192 #ifdef SELDON_CHECK_DIMENSIONS
00193     CheckDim(X, Y, "Swap(X, Y)", "X <-> Y");
00194 #endif
00195 
00196     cblas_cswap(X.GetLength(), reinterpret_cast<void*>(X.GetData()), 1,
00197                 reinterpret_cast<void*>(Y.GetData()), 1);
00198   }
00199 
00200 
00201   template <class Allocator>
00202   void Swap(Vector<complex<double>, VectFull, Allocator>& X,
00203             Vector<complex<double>, VectFull, Allocator>& Y)
00204   {
00205 
00206 #ifdef SELDON_CHECK_DIMENSIONS
00207     CheckDim(X, Y, "Swap(X, Y)", "X <-> Y");
00208 #endif
00209 
00210     cblas_zswap(X.GetLength(), reinterpret_cast<void*>(X.GetData()), 1,
00211                 reinterpret_cast<void*>(Y.GetData()), 1);
00212   }
00213 
00214 
00215   // Swap //
00217 
00218 
00219 
00221   // Mlt //
00222 
00223 
00224   template <class Allocator>
00225   void Mlt(const float alpha,
00226            Vector<float, VectFull, Allocator>& X)
00227   {
00228     cblas_sscal(X.GetLength(), alpha, X.GetData(), 1);
00229   }
00230 
00231 
00232   template <class Allocator>
00233   void Mlt(const double alpha,
00234            Vector<double, VectFull, Allocator>& X)
00235   {
00236     cblas_dscal(X.GetLength(), alpha, X.GetData(), 1);
00237   }
00238 
00239 
00240   template <class Allocator>
00241   void Mlt(const float alpha,
00242            Vector<complex<float>, VectFull, Allocator>& X)
00243   {
00244     cblas_csscal(X.GetLength(), alpha,
00245                  reinterpret_cast<void*>(X.GetData()), 1);
00246   }
00247 
00248 
00249   template <class Allocator>
00250   void Mlt(const double alpha,
00251            Vector<complex<double>, VectFull, Allocator>& X)
00252   {
00253     cblas_zdscal(X.GetLength(), alpha,
00254                  reinterpret_cast<void*>(X.GetData()), 1);
00255   }
00256 
00257 
00258   template <class Allocator>
00259   void Mlt(const complex<float> alpha,
00260            Vector<complex<float>, VectFull, Allocator>& X)
00261   {
00262     cblas_cscal(X.GetLength(),
00263                 reinterpret_cast<const void*>(&alpha),
00264                 reinterpret_cast<void*>(X.GetData()), 1);
00265   }
00266 
00267 
00268   template <class Allocator>
00269   void Mlt(const complex<double> alpha,
00270            Vector<complex<double>, VectFull, Allocator>& X)
00271   {
00272     cblas_zscal(X.GetLength(),
00273                 reinterpret_cast<const void*>(&alpha),
00274                 reinterpret_cast<void*>(X.GetData()), 1);
00275   }
00276 
00277 
00278   // Mlt //
00280 
00281 
00282 
00284   // Copy //
00285 
00286 
00287   template <class Allocator0, class Allocator1>
00288   void Copy(const Vector<float, VectFull, Allocator0>& X,
00289             Vector<float, VectFull, Allocator1>& Y)
00290   {
00291 
00292 #ifdef SELDON_CHECK_DIMENSIONS
00293     CheckDim(X, Y, "Copy(X, Y)", "X -> Y");
00294 #endif
00295 
00296     cblas_scopy(Y.GetLength(),
00297                 reinterpret_cast<const float*>(X.GetData()), 1,
00298                 reinterpret_cast<float*>(Y.GetData()), 1);
00299   }
00300 
00301 
00302   template <class Allocator0, class Allocator1>
00303   void Copy(const Vector<double, VectFull, Allocator0>& X,
00304             Vector<double, VectFull, Allocator1>& Y)
00305   {
00306 
00307 #ifdef SELDON_CHECK_DIMENSIONS
00308     CheckDim(X, Y, "Copy(X, Y)", "X -> Y");
00309 #endif
00310 
00311     cblas_dcopy(Y.GetLength(),
00312                 reinterpret_cast<const double*>(X.GetData()), 1,
00313                 reinterpret_cast<double*>(Y.GetData()), 1);
00314   }
00315 
00316 
00317   template <class Allocator0, class Allocator1>
00318   void Copy(const Vector<complex<float>, VectFull, Allocator0>& X,
00319             Vector<complex<float>, VectFull, Allocator1>& Y)
00320   {
00321 
00322 #ifdef SELDON_CHECK_DIMENSIONS
00323     CheckDim(X, Y, "Copy(X, Y)", "X -> Y");
00324 #endif
00325 
00326     cblas_ccopy(Y.GetLength(),
00327                 reinterpret_cast<const void*>(X.GetData()), 1,
00328                 reinterpret_cast<void*>(Y.GetData()), 1);
00329   }
00330 
00331 
00332   template <class Allocator0, class Allocator1>
00333   void Copy(const Vector<complex<double>, VectFull, Allocator0>& X,
00334             Vector<complex<double>, VectFull, Allocator1>& Y)
00335   {
00336 
00337 #ifdef SELDON_CHECK_DIMENSIONS
00338     CheckDim(X, Y, "Copy(X, Y)", "X -> Y");
00339 #endif
00340 
00341     cblas_zcopy(Y.GetLength(),
00342                 reinterpret_cast<const void*>(X.GetData()), 1,
00343                 reinterpret_cast<void*>(Y.GetData()), 1);
00344   }
00345 
00346 
00347   // Copy //
00349 
00350 
00351 
00353   // Add //
00354 
00355 
00356   template <class Allocator0, class Allocator1>
00357   void Add(const float alpha,
00358            const Vector<float, VectFull, Allocator0>& X,
00359            Vector<float, VectFull, Allocator1>& Y)
00360   {
00361 
00362 #ifdef SELDON_CHECK_DIMENSIONS
00363     CheckDim(X, Y, "Add(alpha, X, Y)");
00364 #endif
00365 
00366     cblas_saxpy(Y.GetLength(),
00367                 alpha,
00368                 reinterpret_cast<const float*>(X.GetData()), 1,
00369                 reinterpret_cast<float*>(Y.GetData()), 1);
00370   }
00371 
00372 
00373   template <class Allocator0, class Allocator1>
00374   void Add(const double alpha,
00375            const Vector<double, VectFull, Allocator0>& X,
00376            Vector<double, VectFull, Allocator1>& Y)
00377   {
00378 
00379 #ifdef SELDON_CHECK_DIMENSIONS
00380     CheckDim(X, Y, "Add(alpha, X, Y)");
00381 #endif
00382 
00383     cblas_daxpy(Y.GetLength(),
00384                 alpha,
00385                 reinterpret_cast<const double*>(X.GetData()), 1,
00386                 reinterpret_cast<double*>(Y.GetData()), 1);
00387   }
00388 
00389 
00390   template <class Allocator0, class Allocator1>
00391   void Add(const complex<float> alpha,
00392            const Vector<complex<float>, VectFull, Allocator0>& X,
00393            Vector<complex<float>, VectFull, Allocator1>& Y)
00394   {
00395 
00396 #ifdef SELDON_CHECK_DIMENSIONS
00397     CheckDim(X, Y, "Add(alpha, X, Y)");
00398 #endif
00399 
00400     cblas_caxpy(Y.GetLength(),
00401                 reinterpret_cast<const void*>(&alpha),
00402                 reinterpret_cast<const void*>(X.GetData()), 1,
00403                 reinterpret_cast<float*>(Y.GetData()), 1);
00404   }
00405 
00406 
00407   template <class Allocator0, class Allocator1>
00408   void Add(const complex<double> alpha,
00409            const Vector<complex<double>, VectFull, Allocator0>& X,
00410            Vector<complex<double>, VectFull, Allocator1>& Y)
00411   {
00412 
00413 #ifdef SELDON_CHECK_DIMENSIONS
00414     CheckDim(X, Y, "Add(alpha, X, Y)");
00415 #endif
00416 
00417     cblas_zaxpy(Y.GetLength(),
00418                 reinterpret_cast<const void*>(&alpha),
00419                 reinterpret_cast<const void*>(X.GetData()), 1,
00420                 reinterpret_cast<double*>(Y.GetData()), 1);
00421   }
00422 
00423 
00424   // Add //
00426 
00427 
00428 
00430   // DotProd //
00431 
00432 
00433   template <class Allocator0, class Allocator1>
00434   float DotProd(const Vector<float, VectFull, Allocator0>& X,
00435                 const Vector<float, VectFull, Allocator1>& Y)
00436   {
00437 
00438 #ifdef SELDON_CHECK_DIMENSIONS
00439     CheckDim(X, Y, "DotProd(X, Y)", "dot(X, Y)");
00440 #endif
00441 
00442     return cblas_sdot(Y.GetLength(),
00443                       reinterpret_cast<const float*>(X.GetData()), 1,
00444                       reinterpret_cast<const float*>(Y.GetData()), 1);
00445   }
00446 
00447 
00448   template <class Allocator0, class Allocator1>
00449   double DotProd(const Vector<double, VectFull, Allocator0>& X,
00450                  const Vector<double, VectFull, Allocator1>& Y)
00451   {
00452 
00453 #ifdef SELDON_CHECK_DIMENSIONS
00454     CheckDim(X, Y, "DotProd(X, Y)", "dot(X, Y)");
00455 #endif
00456 
00457     return cblas_ddot(Y.GetLength(),
00458                       reinterpret_cast<const double*>(X.GetData()), 1,
00459                       reinterpret_cast<const double*>(Y.GetData()), 1);
00460   }
00461 
00462 
00463   template <class Allocator0, class Allocator1>
00464   complex<float>
00465   DotProd(const Vector<complex<float>, VectFull, Allocator0>& X,
00466           const Vector<complex<float>, VectFull, Allocator1>& Y)
00467   {
00468 
00469 #ifdef SELDON_CHECK_DIMENSIONS
00470     CheckDim(X, Y, "DotProd(X, Y)", "dot(X, Y)");
00471 #endif
00472 
00473     complex<float> dotu;
00474     cblas_cdotu_sub(Y.GetLength(),
00475                     reinterpret_cast<const void*>(X.GetData()), 1,
00476                     reinterpret_cast<const void*>(Y.GetData()), 1,
00477                     reinterpret_cast<void*>(&dotu));
00478     return dotu;
00479   }
00480 
00481 
00482   template <class Allocator0, class Allocator1>
00483   complex<double>
00484   DotProd(const Vector<complex<double>, VectFull, Allocator0>& X,
00485           const Vector<complex<double>, VectFull, Allocator1>& Y)
00486   {
00487 
00488 #ifdef SELDON_CHECK_DIMENSIONS
00489     CheckDim(X, Y, "DotProd(X, Y)", "dot(X, Y)");
00490 #endif
00491 
00492     complex<double> dotu;
00493     cblas_zdotu_sub(Y.GetLength(),
00494                     reinterpret_cast<const void*>(X.GetData()), 1,
00495                     reinterpret_cast<const void*>(Y.GetData()), 1,
00496                     reinterpret_cast<void*>(&dotu));
00497     return dotu;
00498   }
00499 
00500 
00501   // DotProd //
00503 
00504 
00505 
00507   // SCALEDDOTPROD //
00508 
00509 
00510   template <class Allocator0, class Allocator1>
00511   float ScaledDotProd(const float alpha,
00512                       const Vector<float, VectFull, Allocator0>& X,
00513                       const Vector<float, VectFull, Allocator1>& Y)
00514   {
00515 
00516 #ifdef SELDON_CHECK_DIMENSIONS
00517     CheckDim(X, Y, "ScaledDotProd(X, Y)", "dot(X, Y)");
00518 #endif
00519 
00520     return cblas_sdsdot(Y.GetLength(), alpha,
00521                         reinterpret_cast<const float*>(X.GetData()), 1,
00522                         reinterpret_cast<const float*>(Y.GetData()), 1);
00523   }
00524 
00525 
00526   // SCALEDDOTPROD //
00528 
00529 
00530 
00532   // DotProjConj //
00533 
00534 
00535   template <class Allocator0, class Allocator1>
00536   complex<float>
00537   DotProdConj(const Vector<complex<float>, VectFull, Allocator0>& X,
00538               const Vector<complex<float>, VectFull, Allocator1>& Y)
00539   {
00540 
00541 #ifdef SELDON_CHECK_DIMENSIONS
00542     CheckDim(X, Y, "DotProdConj(X, Y)", "dot(X, Y)");
00543 #endif
00544 
00545     complex<float> dotc;
00546     cblas_cdotc_sub(Y.GetLength(),
00547                     reinterpret_cast<const void*>(X.GetData()), 1,
00548                     reinterpret_cast<const void*>(Y.GetData()), 1,
00549                     reinterpret_cast<void*>(&dotc));
00550     return dotc;
00551   }
00552 
00553 
00554   template <class Allocator0, class Allocator1>
00555   complex<double>
00556   DotProdConj(const Vector<complex<double>, VectFull, Allocator0>& X,
00557               const Vector<complex<double>, VectFull, Allocator1>& Y)
00558   {
00559 
00560 #ifdef SELDON_CHECK_DIMENSIONS
00561     CheckDim(X, Y, "DotProdConj(X, Y)", "dot(X, Y)");
00562 #endif
00563 
00564     complex<double> dotc;
00565     cblas_zdotc_sub(Y.GetLength(),
00566                     reinterpret_cast<const void*>(X.GetData()), 1,
00567                     reinterpret_cast<const void*>(Y.GetData()), 1,
00568                     reinterpret_cast<void*>(&dotc));
00569     return dotc;
00570   }
00571 
00572 
00573   // DotProdConj //
00575 
00576 
00577 
00579   // Norm1 //
00580 
00581 
00582   template <class Allocator>
00583   float Norm1(const Vector<float, VectFull, Allocator>& X)
00584   {
00585     return cblas_sasum(X.GetLength(),
00586                        reinterpret_cast<const float*>(X.GetData()), 1);
00587   }
00588 
00589 
00590   template <class Allocator>
00591   double Norm1(const Vector<double, VectFull, Allocator>& X)
00592   {
00593     return cblas_dasum(X.GetLength(),
00594                        reinterpret_cast<const double*>(X.GetData()), 1);
00595   }
00596 
00597 
00598   template <class Allocator>
00599   float Norm1(const Vector<complex<float>, VectFull, Allocator>& X)
00600   {
00601     return cblas_scasum(X.GetLength(),
00602                         reinterpret_cast<const void*>(X.GetData()), 1);
00603   }
00604 
00605 
00606   template <class Allocator>
00607   double Norm1(const Vector<complex<double>, VectFull, Allocator>& X)
00608   {
00609     return cblas_dzasum(X.GetLength(),
00610                         reinterpret_cast<const void*>(X.GetData()), 1);
00611   }
00612 
00613 
00614   // Norm1 //
00616 
00617 
00618 
00620   // Norm2 //
00621 
00622 
00623   template <class Allocator>
00624   float Norm2(const Vector<float, VectFull, Allocator>& X)
00625   {
00626     return cblas_snrm2(X.GetLength(),
00627                        reinterpret_cast<const float*>(X.GetData()), 1);
00628   }
00629 
00630 
00631   template <class Allocator>
00632   double Norm2(const Vector<double, VectFull, Allocator>& X)
00633   {
00634     return cblas_dnrm2(X.GetLength(),
00635                        reinterpret_cast<const double*>(X.GetData()), 1);
00636   }
00637 
00638 
00639   template <class Allocator>
00640   float Norm2(const Vector<complex<float>, VectFull, Allocator>& X)
00641   {
00642     return cblas_scnrm2(X.GetLength(),
00643                         reinterpret_cast<const void*>(X.GetData()), 1);
00644   }
00645 
00646 
00647   template <class Allocator>
00648   double Norm2(const Vector<complex<double>, VectFull, Allocator>& X)
00649   {
00650     return cblas_dznrm2(X.GetLength(),
00651                         reinterpret_cast<const void*>(X.GetData()), 1);
00652   }
00653 
00654 
00655   // Norm2 //
00657 
00658 
00659 
00661   // GetMaxAbsIndex //
00662 
00663 
00664   template <class Allocator>
00665   size_t GetMaxAbsIndex(const Vector<float, VectFull, Allocator>& X)
00666   {
00667     return cblas_isamax(X.GetLength(),
00668                         reinterpret_cast<const float*>(X.GetData()), 1);
00669   }
00670 
00671 
00672   template <class Allocator>
00673   size_t GetMaxAbsIndex(const Vector<double, VectFull, Allocator>& X)
00674   {
00675     return cblas_idamax(X.GetLength(),
00676                         reinterpret_cast<const double*>(X.GetData()), 1);
00677   }
00678 
00679 
00680   template <class Allocator>
00681   size_t GetMaxAbsIndex(const Vector<complex<float>, VectFull, Allocator>& X)
00682   {
00683     return cblas_icamax(X.GetLength(),
00684                         reinterpret_cast<const void*>(X.GetData()), 1);
00685   }
00686 
00687 
00688   template <class Allocator>
00689   size_t
00690   GetMaxAbsIndex(const Vector<complex<double>, VectFull, Allocator>& X)
00691   {
00692     return cblas_izamax(X.GetLength(),
00693                         reinterpret_cast<const void*>(X.GetData()), 1);
00694   }
00695 
00696 
00697   // GetMaxAbsIndex //
00699 
00700 
00701 } // namespace Seldon.
00702 
00703 #define SELDON_FILE_BLAS_1_CXX
00704 #endif
Email Contact