CIMPL: Blas3.cpp Source File

00001 //Copyright (c) 2004-2005, Baris Sumengen
00002 //All rights reserved.
00003 //
00004 // CIMPL Matrix Performance Library
00005 //
00006 //Redistribution and use in source and binary
00007 //forms, with or without modification, are
00008 //permitted provided that the following
00009 //conditions are met:
00010 //
00011 //    * No commercial use is allowed. 
00012 //    This software can only be used
00013 //    for non-commercial purposes. This 
00014 //    distribution is mainly intended for
00015 //    academic research and teaching.
00016 //    * Redistributions of source code must
00017 //    retain the above copyright notice, this
00018 //    list of conditions and the following
00019 //    disclaimer.
00020 //    * Redistributions of binary form must
00021 //    mention the above copyright notice, this
00022 //    list of conditions and the following
00023 //    disclaimer in a clearly visible part 
00024 //    in associated product manual, 
00025 //    readme, and web site of the redistributed 
00026 //    software.
00027 //    * Redistributions in binary form must
00028 //    reproduce the above copyright notice,
00029 //    this list of conditions and the
00030 //    following disclaimer in the
00031 //    documentation and/or other materials
00032 //    provided with the distribution.
00033 //    * The name of Baris Sumengen may not be
00034 //    used to endorse or promote products
00035 //    derived from this software without
00036 //    specific prior written permission.
00037 //
00038 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT
00039 //HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
00040 //EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
00041 //NOT LIMITED TO, THE IMPLIED WARRANTIES OF
00042 //MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00043 //PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
00044 //CONTRIBUTORS BE LIABLE FOR ANY
00045 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00046 //EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00047 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
00048 //OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
00049 //DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00050 //HOWEVER CAUSED AND ON ANY THEORY OF
00051 //LIABILITY, WHETHER IN CONTRACT, STRICT
00052 //LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
00053 //OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
00054 //OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00055 //POSSIBILITY OF SUCH DAMAGE.
00056 
00057 
00058 
00059 #include "./Blas.h"
00060 #include "mkl.h"
00061 
00062 
00063 
00064 
00065 namespace Blas
00066 {
00067 
00068 
00072         Matrix<float> Gemm(Matrix<float>& A, Matrix<float>& B)
00073         {
00074                 if(A.Columns() != B.Rows())
00075                 {
00076                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00077                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00078                 }
00079                 Matrix<float> C(A.Rows(), B.Columns(), 0);
00080                 cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00081                 return C;
00082         }
00083 
00084         Matrix<double> Gemm(Matrix<double>& A, Matrix<double>& B)
00085         {
00086                 if(A.Columns() != B.Rows())
00087                 {
00088                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00089                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00090                 }
00091                 Matrix<double> C(A.Rows(), B.Columns(), 0);
00092                 cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00093                 return C;
00094         }
00095 
00096         Matrix<ComplexFloat> Gemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00097         {
00098                 if(A.Columns() != B.Rows())
00099                 {
00100                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00101                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00102                 }
00103                 ComplexFloat temp(1,0);
00104                 Matrix<ComplexFloat> C(A.Rows(), B.Columns(), ComplexFloat(0,0));
00105                 cblas_cgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00106                 return C;
00107         }
00108 
00109         Matrix<ComplexDouble> Gemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00110         {
00111                 if(A.Columns() != B.Rows())
00112                 {
00113                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00114                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00115                 }
00116                 ComplexDouble temp(1,0);
00117                 Matrix<ComplexDouble> C(A.Rows(), B.Columns(), ComplexDouble(0,0));
00118                 cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, A.Rows(), B.Columns(), A.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00119                 return C;
00120         }
00121 
00122 
00123 
00126         Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B)
00127         {
00128                 if(A.Columns() != A.Rows())
00129                 {
00130                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00131                         Utility::RunTimeError("Matrix is not square!");
00132                 }
00133                 if(A.Columns() != B.Rows())
00134                 {
00135                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00136                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00137                 }
00138                 Matrix<float> C(B.Rows(), B.Columns(), 0);
00139                 cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00140                 return C;
00141         }
00142 
00143         Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B)
00144         {
00145                 if(A.Columns() != A.Rows())
00146                 {
00147                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00148                         Utility::RunTimeError("Matrix is not square!");
00149                 }
00150                 if(A.Columns() != B.Rows())
00151                 {
00152                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00153                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00154                 }
00155                 Matrix<double> C(B.Rows(), B.Columns(), 0);
00156                 cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00157                 return C;
00158         }
00159 
00160         Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00161         {
00162                 if(A.Columns() != A.Rows())
00163                 {
00164                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00165                         Utility::RunTimeError("Matrix is not square!");
00166                 }
00167                 if(A.Columns() != B.Rows())
00168                 {
00169                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00170                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00171                 }
00172                 ComplexFloat temp(1,0);
00173                 Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00174                 cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00175                 return C;
00176         }
00177 
00178         Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00179         {
00180                 if(A.Columns() != A.Rows())
00181                 {
00182                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00183                         Utility::RunTimeError("Matrix is not square!");
00184                 }
00185                 if(A.Columns() != B.Rows())
00186                 {
00187                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00188                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00189                 }
00190                 ComplexDouble temp(1,0);
00191                 Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00192                 cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00193                 return C;
00194         }
00195 
00196         Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00197         {
00198                 if(A.Columns() != A.Rows())
00199                 {
00200                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00201                         Utility::RunTimeError("Matrix is not square!");
00202                 }
00203                 if(A.Columns() != B.Rows())
00204                 {
00205                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00206                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00207                 }
00208                 ComplexFloat temp(1,0);
00209                 Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00210                 cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00211                 return C;
00212         }
00213 
00214         Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00215         {
00216                 if(A.Columns() != A.Rows())
00217                 {
00218                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00219                         Utility::RunTimeError("Matrix is not square!");
00220                 }
00221                 if(A.Columns() != B.Rows())
00222                 {
00223                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00224                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00225                 }
00226                 ComplexDouble temp(1,0);
00227                 Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00228                 cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00229                 return C;
00230         }
00231 
00232 
00233         Matrix<float> Symm(Matrix<float>& A, Matrix<float>& B, bool orderReversed)
00234         {
00235                 if(A.Columns() != A.Rows())
00236                 {
00237                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00238                         Utility::RunTimeError("Matrix is not square!");
00239                 }
00240                 if(orderReversed)
00241                 {
00242                         if(A.Rows() != B.Columns())
00243                         {
00244                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00245                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00246                         }
00247                         Matrix<float> C(B.Rows(), B.Columns(), 0);
00248                         cblas_ssymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00249                         return C;
00250                 }
00251                 else
00252                 {
00253                         if(A.Columns() != B.Rows())
00254                         {
00255                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00256                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00257                         }
00258                         Matrix<float> C(B.Rows(), B.Columns(), 0);
00259                         cblas_ssymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00260                         return C;
00261                 }
00262         }
00263 
00264         Matrix<double> Symm(Matrix<double>& A, Matrix<double>& B, bool orderReversed)
00265         {
00266                 if(A.Columns() != A.Rows())
00267                 {
00268                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00269                         Utility::RunTimeError("Matrix is not square!");
00270                 }
00271                 if(orderReversed)
00272                 {
00273                         if(A.Rows() != B.Columns())
00274                         {
00275                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00276                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00277                         }
00278                         Matrix<double> C(B.Rows(), B.Columns(), 0);
00279                         cblas_dsymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00280                         return C;
00281                 }
00282                 else
00283                 {
00284                         if(A.Columns() != B.Rows())
00285                         {
00286                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00287                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00288                         }
00289                         Matrix<double> C(B.Rows(), B.Columns(), 0);
00290                         cblas_dsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), 1, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00291                         return C;
00292                 }
00293         }
00294 
00295         Matrix<ComplexFloat> Symm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed)
00296         {
00297                 if(A.Columns() != A.Rows())
00298                 {
00299                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00300                         Utility::RunTimeError("Matrix is not square!");
00301                 }
00302                 if(orderReversed)
00303                 {
00304                         if(A.Rows() != B.Columns())
00305                         {
00306                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00307                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00308                         }
00309                         ComplexFloat temp(1,0);
00310                         Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00311                         cblas_csymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00312                         return C;
00313                 }
00314                 else
00315                 {
00316                         if(A.Columns() != B.Rows())
00317                         {
00318                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00319                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00320                         }
00321                         ComplexFloat temp(1,0);
00322                         Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00323                         cblas_csymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00324                         return C;
00325                 }
00326         }
00327 
00328         Matrix<ComplexDouble> Symm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed)
00329         {
00330                 if(A.Columns() != A.Rows())
00331                 {
00332                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00333                         Utility::RunTimeError("Matrix is not square!");
00334                 }
00335                 if(orderReversed)
00336                 {
00337                         if(A.Rows() != B.Columns())
00338                         {
00339                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00340                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00341                         }
00342                         ComplexDouble temp(1,0);
00343                         Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00344                         cblas_zsymm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00345                         return C;
00346                 }
00347                 else
00348                 {
00349                         if(A.Columns() != B.Rows())
00350                         {
00351                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00352                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00353                         }
00354                         ComplexDouble temp(1,0);
00355                         Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00356                         cblas_zsymm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00357                         return C;
00358                 }
00359         }
00360 
00361         Matrix<ComplexFloat> Hemm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, bool orderReversed)
00362         {
00363                 if(A.Columns() != A.Rows())
00364                 {
00365                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00366                         Utility::RunTimeError("Matrix is not square!");
00367                 }
00368                 if(orderReversed)
00369                 {
00370                         if(A.Rows() != B.Columns())
00371                         {
00372                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00373                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00374                         }
00375                         ComplexFloat temp(1,0);
00376                         Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00377                         cblas_chemm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00378                         return C;
00379                 }
00380                 else
00381                 {
00382                         if(A.Columns() != B.Rows())
00383                         {
00384                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00385                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00386                         }
00387                         ComplexFloat temp(1,0);
00388                         Matrix<ComplexFloat> C(B.Rows(), B.Columns(), ComplexFloat(0,0));
00389                         cblas_chemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00390                         return C;
00391                 }
00392         }
00393 
00394         Matrix<ComplexDouble> Hemm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, bool orderReversed)
00395         {
00396                 if(A.Columns() != A.Rows())
00397                 {
00398                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00399                         Utility::RunTimeError("Matrix is not square!");
00400                 }
00401                 if(orderReversed)
00402                 {
00403                         if(A.Rows() != B.Columns())
00404                         {
00405                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00406                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00407                         }
00408                         ComplexDouble temp(1,0);
00409                         Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00410                         cblas_zhemm(CblasColMajor, CblasRight, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00411                         return C;
00412                 }
00413                 else
00414                 {
00415                         if(A.Columns() != B.Rows())
00416                         {
00417                                 cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00418                                 Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00419                         }
00420                         ComplexDouble temp(1,0);
00421                         Matrix<ComplexDouble> C(B.Rows(), B.Columns(), ComplexDouble(0,0));
00422                         cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, B.Rows(), B.Columns(), &temp, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00423                         return C;
00424                 }
00425         }
00426 
00427 
00428 
00432 
00433         void Syrk(float alpha, Matrix<float>& A, Matrix<float>& C)
00434         {
00435                 if(C.Columns() != C.Rows())
00436                 {
00437                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00438                         Utility::RunTimeError("Matrix is not square!");
00439                 }
00440                 if(C.Rows() != A.Rows())
00441                 {
00442                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00443                         Utility::RunTimeError("Matrix dimensions does not match!");
00444                 }
00445                 cblas_ssyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
00446         }
00447 
00448         void Syrk(double alpha, Matrix<double>& A, Matrix<double>& C)
00449         {
00450                 if(C.Columns() != C.Rows())
00451                 {
00452                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00453                         Utility::RunTimeError("Matrix is not square!");
00454                 }
00455                 if(C.Rows() != A.Rows())
00456                 {
00457                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00458                         Utility::RunTimeError("Matrix dimensions does not match!");
00459                 }
00460                 cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
00461         }
00462 
00463         void Syrk(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& C)
00464         {
00465                 if(C.Columns() != C.Rows())
00466                 {
00467                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00468                         Utility::RunTimeError("Matrix is not square!");
00469                 }
00470                 if(C.Rows() != A.Rows())
00471                 {
00472                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00473                         Utility::RunTimeError("Matrix dimensions does not match!");
00474                 }
00475                 ComplexFloat temp(1,0);
00476                 cblas_csyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), &temp, C.Data(), C.Rows());
00477         }
00478 
00479         void Syrk(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C)
00480         {
00481                 if(C.Columns() != C.Rows())
00482                 {
00483                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00484                         Utility::RunTimeError("Matrix is not square!");
00485                 }
00486                 if(C.Rows() != A.Rows())
00487                 {
00488                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00489                         Utility::RunTimeError("Matrix dimensions does not match!");
00490                 }
00491                 ComplexDouble temp(1,0);
00492                 cblas_zsyrk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), &temp, C.Data(), C.Rows());
00493         }
00494 
00495 
00496 
00500 
00501         void Syr2k(float alpha, Matrix<float>& A, Matrix<float>& B, Matrix<float>& C)
00502         {
00503                 if(C.Columns() != C.Rows())
00504                 {
00505                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00506                         Utility::RunTimeError("Matrix is not square!");
00507                 }
00508                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00509                 {
00510                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00511                         Utility::RunTimeError("Matrix dimensions does not match!");
00512                 }
00513                 cblas_ssyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00514         }
00515 
00516         void Syr2k(double alpha, Matrix<double>& A, Matrix<double>& B, Matrix<double>& C)
00517         {
00518                 if(C.Columns() != C.Rows())
00519                 {
00520                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00521                         Utility::RunTimeError("Matrix is not square!");
00522                 }
00523                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00524                 {
00525                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00526                         Utility::RunTimeError("Matrix dimensions does not match!");
00527                 }
00528                 cblas_dsyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00529         }
00530 
00531         void Syr2k(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& C)
00532         {
00533                 if(C.Columns() != C.Rows())
00534                 {
00535                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00536                         Utility::RunTimeError("Matrix is not square!");
00537                 }
00538                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00539                 {
00540                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00541                         Utility::RunTimeError("Matrix dimensions does not match!");
00542                 }
00543                 ComplexFloat temp(1,0);
00544                 cblas_csyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00545         }
00546 
00547         void Syr2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C)
00548         {
00549                 if(C.Columns() != C.Rows())
00550                 {
00551                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00552                         Utility::RunTimeError("Matrix is not square!");
00553                 }
00554                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00555                 {
00556                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00557                         Utility::RunTimeError("Matrix dimensions does not match!");
00558                 }
00559                 ComplexDouble temp(1,0);
00560                 cblas_zsyr2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), &temp, C.Data(), C.Rows());
00561         }
00562 
00563 
00564 
00565         // Herk: rank k update of an hermitian matrix
00568 
00569         void Herk(float alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& C)
00570         {
00571                 if(C.Columns() != C.Rows())
00572                 {
00573                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00574                         Utility::RunTimeError("Matrix is not square!");
00575                 }
00576                 if(C.Rows() != A.Rows())
00577                 {
00578                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00579                         Utility::RunTimeError("Matrix dimensions does not match!");
00580                 }
00581                 cblas_cherk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
00582         }
00583 
00584         void Herk(double alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& C)
00585         {
00586                 if(C.Columns() != C.Rows())
00587                 {
00588                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00589                         Utility::RunTimeError("Matrix is not square!");
00590                 }
00591                 if(C.Rows() != A.Rows())
00592                 {
00593                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00594                         Utility::RunTimeError("Matrix dimensions does not match!");
00595                 }
00596                 cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), alpha, A.Data(), A.Rows(), 1, C.Data(), C.Rows());
00597         }
00598 
00599 
00600         // Her2k: rank 2k update of an hermitian matrix
00603 
00604         void Her2k(ComplexFloat alpha, Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B, Matrix<ComplexFloat>& C)
00605         {
00606                 if(C.Columns() != C.Rows())
00607                 {
00608                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00609                         Utility::RunTimeError("Matrix is not square!");
00610                 }
00611                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00612                 {
00613                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00614                         Utility::RunTimeError("Matrix dimensions does not match!");
00615                 }
00616                 cblas_cher2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00617         }
00618 
00619         void Her2k(ComplexDouble alpha, Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B, Matrix<ComplexDouble>& C)
00620         {
00621                 if(C.Columns() != C.Rows())
00622                 {
00623                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00624                         Utility::RunTimeError("Matrix is not square!");
00625                 }
00626                 if(C.Rows() != A.Rows() || C.Rows() != B.Rows() || A.Columns() != B.Columns())
00627                 {
00628                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00629                         Utility::RunTimeError("Matrix dimensions does not match!");
00630                 }
00631                 cblas_zher2k(CblasColMajor, CblasUpper, CblasNoTrans, A.Rows(), A.Columns(), &alpha, A.Data(), A.Rows(), B.Data(), B.Rows(), 1, C.Data(), C.Rows());
00632         }
00633 
00634 
00635 
00638         Matrix<float> Trmm(Matrix<float>& A, Matrix<float>& B)
00639         {
00640                 if(A.Columns() != A.Rows())
00641                 {
00642                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00643                         Utility::RunTimeError("Matrix is not square!");
00644                 }
00645                 if(A.Columns() != B.Rows())
00646                 {
00647                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00648                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00649                 }
00650                 Matrix<float> C = B.Clone();
00651                 cblas_strmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), 1, A.Data(), A.Rows(), C.Data(), C.Rows());
00652                 return C;
00653         }
00654 
00655         Matrix<double> Trmm(Matrix<double>& A, Matrix<double>& B)
00656         {
00657                 if(A.Columns() != A.Rows())
00658                 {
00659                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00660                         Utility::RunTimeError("Matrix is not square!");
00661                 }
00662                 if(A.Columns() != B.Rows())
00663                 {
00664                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00665                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00666                 }
00667                 Matrix<double> C = B.Clone();
00668                 cblas_dtrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), 1, A.Data(), A.Rows(), C.Data(), C.Rows());
00669                 return C;
00670         }
00671 
00672         Matrix<ComplexFloat> Trmm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00673         {
00674                 if(A.Columns() != A.Rows())
00675                 {
00676                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00677                         Utility::RunTimeError("Matrix is not square!");
00678                 }
00679                 if(A.Columns() != B.Rows())
00680                 {
00681                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00682                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00683                 }
00684                 ComplexFloat temp(1,0);
00685                 Matrix<ComplexFloat> C = B.Clone();
00686                 cblas_ctrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), &temp, A.Data(), A.Rows(), C.Data(), C.Rows());
00687                 return C;
00688         }
00689 
00690         Matrix<ComplexDouble> Trmm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00691         {
00692                 if(A.Columns() != A.Rows())
00693                 {
00694                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00695                         Utility::RunTimeError("Matrix is not square!");
00696                 }
00697                 if(A.Columns() != B.Rows())
00698                 {
00699                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00700                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00701                 }
00702                 ComplexDouble temp(1,0);
00703                 Matrix<ComplexDouble> C = B.Clone();
00704                 cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, C.Rows(), C.Columns(), &temp, A.Data(), A.Rows(), C.Data(), C.Rows());
00705                 return C;
00706         }
00707 
00708 
00709 
00710 
00711 
00712         // Trsm: Solve a triangular system of equations with a triangular coefficient matrix
00714         Matrix<float> Trsm(Matrix<float>& A, Matrix<float>& B)
00715         {
00716                 if(A.Columns() != A.Rows())
00717                 {
00718                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00719                         Utility::RunTimeError("Matrix is not square!");
00720                 }
00721                 if(A.Columns() != B.Rows())
00722                 {
00723                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00724                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00725                 }
00726                 Matrix<float> X = B.Clone();
00727                 cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), 1, A.Data(), A.Rows(), X.Data(), X.Rows());
00728                 return X;
00729         }
00730 
00731         Matrix<double> Trsm(Matrix<double>& A, Matrix<double>& B)
00732         {
00733                 if(A.Columns() != A.Rows())
00734                 {
00735                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00736                         Utility::RunTimeError("Matrix is not square!");
00737                 }
00738                 if(A.Columns() != B.Rows())
00739                 {
00740                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00741                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00742                 }
00743                 Matrix<double> X = B.Clone();
00744                 cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), 1, A.Data(), A.Rows(), X.Data(), X.Rows());
00745                 return X;
00746         }
00747 
00748         Matrix<ComplexFloat> Trsm(Matrix<ComplexFloat>& A, Matrix<ComplexFloat>& B)
00749         {
00750                 if(A.Columns() != A.Rows())
00751                 {
00752                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00753                         Utility::RunTimeError("Matrix is not square!");
00754                 }
00755                 if(A.Columns() != B.Rows())
00756                 {
00757                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00758                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00759                 }
00760                 ComplexFloat temp(1,0);
00761                 Matrix<ComplexFloat> X = B.Clone();
00762                 cblas_ctrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), &temp, A.Data(), A.Rows(), X.Data(), X.Rows());
00763                 return X;
00764         }
00765 
00766         Matrix<ComplexDouble> Trsm(Matrix<ComplexDouble>& A, Matrix<ComplexDouble>& B)
00767         {
00768                 if(A.Columns() != A.Rows())
00769                 {
00770                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00771                         Utility::RunTimeError("Matrix is not square!");
00772                 }
00773                 if(A.Columns() != B.Rows())
00774                 {
00775                         cerr << "Line: " << __LINE__ << " File: " << __FILE__ << endl;
00776                         Utility::RunTimeError("Matrix dimensions does not match for matrix multiplication!");
00777                 }
00778                 ComplexDouble temp(1,0);
00779                 Matrix<ComplexDouble> X = B.Clone();
00780                 cblas_ztrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, X.Rows(), X.Columns(), &temp, A.Data(), A.Rows(), X.Data(), X.Rows());
00781                 return X;
00782         }
00783 
00784 
00785 
00786 
00787 };
00788 
00789 
00790 
00791 
00792 
00793 
00794 
00795 
00796 
00797 
00798 
00799 
00800 
00801 
00802 
00803 
00804 
00805 
00806 
00807 
00808 
00809