level3_impl.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "common.h"
00026 
00027 int EIGEN_BLAS_FUNC(gemm)(char *opa, char *opb, int *m, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
00028 {
00029 //   std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
00030   typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*);
00031   static functype func[12];
00032 
00033   static bool init = false;
00034   if(!init)
00035   {
00036     for(int k=0; k<12; ++k)
00037       func[k] = 0;
00038     func[NOTR  | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run);
00039     func[TR    | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run);
00040     func[ADJ   | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run);
00041     func[NOTR  | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run);
00042     func[TR    | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run);
00043     func[ADJ   | (TR   << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run);
00044     func[NOTR  | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
00045     func[TR    | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
00046     func[ADJ   | (ADJ  << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run);
00047     init = true;
00048   }
00049 
00050   Scalar* a = reinterpret_cast<Scalar*>(pa);
00051   Scalar* b = reinterpret_cast<Scalar*>(pb);
00052   Scalar* c = reinterpret_cast<Scalar*>(pc);
00053   Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
00054   Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);
00055 
00056   int info = 0;
00057   if(OP(*opa)==INVALID)                                               info = 1;
00058   else if(OP(*opb)==INVALID)                                          info = 2;
00059   else if(*m<0)                                                       info = 3;
00060   else if(*n<0)                                                       info = 4;
00061   else if(*k<0)                                                       info = 5;
00062   else if(*lda<std::max(1,(OP(*opa)==NOTR)?*m:*k))                    info = 8;
00063   else if(*ldb<std::max(1,(OP(*opb)==NOTR)?*k:*n))                    info = 10;
00064   else if(*ldc<std::max(1,*m))                                        info = 13;
00065   if(info)
00066     return xerbla_(SCALAR_SUFFIX_UP"GEMM ",&info,6);
00067 
00068   if(beta!=Scalar(1))
00069   {
00070     if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
00071     else                matrix(c, *m, *n, *ldc) *= beta;
00072   }
00073 
00074   internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k);
00075 
00076   int code = OP(*opa) | (OP(*opb) << 2);
00077   func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0);
00078   return 0;
00079 }
00080 
00081 int EIGEN_BLAS_FUNC(trsm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
00082 {
00083 //   std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
00084   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex);
00085   static functype func[32];
00086 
00087   static bool init = false;
00088   if(!init)
00089   {
00090     for(int k=0; k<32; ++k)
00091       func[k] = 0;
00092 
00093     func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,ColMajor,ColMajor>::run);
00094     func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,RowMajor,ColMajor>::run);
00095     func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          Conj, RowMajor,ColMajor>::run);
00096 
00097     func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,ColMajor,ColMajor>::run);
00098     func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,RowMajor,ColMajor>::run);
00099     func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          Conj, RowMajor,ColMajor>::run);
00100 
00101     func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0,          false,ColMajor,ColMajor>::run);
00102     func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          false,RowMajor,ColMajor>::run);
00103     func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0,          Conj, RowMajor,ColMajor>::run);
00104 
00105     func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0,          false,ColMajor,ColMajor>::run);
00106     func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          false,RowMajor,ColMajor>::run);
00107     func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0,          Conj, RowMajor,ColMajor>::run);
00108 
00109 
00110     func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run);
00111     func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run);
00112     func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
00113 
00114     func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run);
00115     func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run);
00116     func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
00117 
00118     func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run);
00119     func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run);
00120     func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
00121 
00122     func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run);
00123     func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run);
00124     func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
00125 
00126     init = true;
00127   }
00128 
00129   Scalar* a = reinterpret_cast<Scalar*>(pa);
00130   Scalar* b = reinterpret_cast<Scalar*>(pb);
00131   Scalar  alpha = *reinterpret_cast<Scalar*>(palpha);
00132 
00133   int info = 0;
00134   if(SIDE(*side)==INVALID)                                            info = 1;
00135   else if(UPLO(*uplo)==INVALID)                                       info = 2;
00136   else if(OP(*opa)==INVALID)                                          info = 3;
00137   else if(DIAG(*diag)==INVALID)                                       info = 4;
00138   else if(*m<0)                                                       info = 5;
00139   else if(*n<0)                                                       info = 6;
00140   else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 9;
00141   else if(*ldb<std::max(1,*m))                                        info = 11;
00142   if(info)
00143     return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6);
00144 
00145   int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
00146 
00147   if(SIDE(*side)==LEFT)
00148     func[code](*m, *n, a, *lda, b, *ldb);
00149   else
00150     func[code](*n, *m, a, *lda, b, *ldb);
00151 
00152   if(alpha!=Scalar(1))
00153     matrix(b,*m,*n,*ldb) *= alpha;
00154 
00155   return 0;
00156 }
00157 
00158 
00159 // b = alpha*op(a)*b  for side = 'L'or'l'
00160 // b = alpha*b*op(a)  for side = 'R'or'r'
00161 int EIGEN_BLAS_FUNC(trmm)(char *side, char *uplo, char *opa, char *diag, int *m, int *n, RealScalar *palpha,  RealScalar *pa, int *lda, RealScalar *pb, int *ldb)
00162 {
00163 //   std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n";
00164   typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
00165   static functype func[32];
00166   static bool init = false;
00167   if(!init)
00168   {
00169     for(int k=0; k<32; ++k)
00170       func[k] = 0;
00171 
00172     func[NOTR  | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
00173     func[TR    | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
00174     func[ADJ   | (LEFT  << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
00175 
00176     func[NOTR  | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
00177     func[TR    | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
00178     func[ADJ   | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
00179 
00180     func[NOTR  | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          true, ColMajor,false,ColMajor,false,ColMajor>::run);
00181     func[TR    | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,false,ColMajor,false,ColMajor>::run);
00182     func[ADJ   | (LEFT  << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
00183 
00184     func[NOTR  | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0,          false,ColMajor,false,ColMajor,false,ColMajor>::run);
00185     func[TR    | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,false,ColMajor>::run);
00186     func[ADJ   | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0,          false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
00187 
00188     func[NOTR  | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
00189     func[TR    | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
00190     func[ADJ   | (LEFT  << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
00191 
00192     func[NOTR  | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
00193     func[TR    | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
00194     func[ADJ   | (RIGHT << 2) | (UP << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
00195 
00196     func[NOTR  | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
00197     func[TR    | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
00198     func[ADJ   | (LEFT  << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
00199 
00200     func[NOTR  | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
00201     func[TR    | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
00202     func[ADJ   | (RIGHT << 2) | (LO << 3) | (UNIT  << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
00203 
00204     init = true;
00205   }
00206 
00207   Scalar* a = reinterpret_cast<Scalar*>(pa);
00208   Scalar* b = reinterpret_cast<Scalar*>(pb);
00209   Scalar  alpha = *reinterpret_cast<Scalar*>(palpha);
00210 
00211   int info = 0;
00212   if(SIDE(*side)==INVALID)                                            info = 1;
00213   else if(UPLO(*uplo)==INVALID)                                       info = 2;
00214   else if(OP(*opa)==INVALID)                                          info = 3;
00215   else if(DIAG(*diag)==INVALID)                                       info = 4;
00216   else if(*m<0)                                                       info = 5;
00217   else if(*n<0)                                                       info = 6;
00218   else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 9;
00219   else if(*ldb<std::max(1,*m))                                        info = 11;
00220   if(info)
00221     return xerbla_(SCALAR_SUFFIX_UP"TRMM ",&info,6);
00222 
00223   int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
00224 
00225   if(*m==0 || *n==0)
00226     return 1;
00227 
00228   // FIXME find a way to avoid this copy
00229   Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb);
00230   matrix(b,*m,*n,*ldb).setZero();
00231 
00232   if(SIDE(*side)==LEFT)
00233     func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha);
00234   else
00235     func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha);
00236   return 1;
00237 }
00238 
00239 // c = alpha*a*b + beta*c  for side = 'L'or'l'
00240 // c = alpha*b*a + beta*c  for side = 'R'or'r
00241 int EIGEN_BLAS_FUNC(symm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
00242 {
00243 //   std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n";
00244   Scalar* a = reinterpret_cast<Scalar*>(pa);
00245   Scalar* b = reinterpret_cast<Scalar*>(pb);
00246   Scalar* c = reinterpret_cast<Scalar*>(pc);
00247   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00248   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
00249 
00250   int info = 0;
00251   if(SIDE(*side)==INVALID)                                            info = 1;
00252   else if(UPLO(*uplo)==INVALID)                                       info = 2;
00253   else if(*m<0)                                                       info = 3;
00254   else if(*n<0)                                                       info = 4;
00255   else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 7;
00256   else if(*ldb<std::max(1,*m))                                        info = 9;
00257   else if(*ldc<std::max(1,*m))                                        info = 12;
00258   if(info)
00259     return xerbla_(SCALAR_SUFFIX_UP"SYMM ",&info,6);
00260 
00261   if(beta!=Scalar(1))
00262   {
00263     if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
00264     else                matrix(c, *m, *n, *ldc) *= beta;
00265   }
00266 
00267   if(*m==0 || *n==0)
00268   {
00269     return 1;
00270   }
00271 
00272   #if ISCOMPLEX
00273   // FIXME add support for symmetric complex matrix
00274   int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
00275   Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size);
00276   if(UPLO(*uplo)==UP)
00277   {
00278     matA.triangularView<Upper>() = matrix(a,size,size,*lda);
00279     matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose();
00280   }
00281   else if(UPLO(*uplo)==LO)
00282   {
00283     matA.triangularView<Lower>() = matrix(a,size,size,*lda);
00284     matA.triangularView<Upper>() = matrix(a,size,size,*lda).transpose();
00285   }
00286   if(SIDE(*side)==LEFT)
00287     matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb);
00288   else if(SIDE(*side)==RIGHT)
00289     matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA;
00290   #else
00291   if(SIDE(*side)==LEFT)
00292     if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar, DenseIndex, RowMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
00293     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,true,false, ColMajor,false,false, ColMajor>::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
00294     else                      return 0;
00295   else if(SIDE(*side)==RIGHT)
00296     if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, RowMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
00297     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar, DenseIndex, ColMajor,false,false, ColMajor,true,false, ColMajor>::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
00298     else                      return 0;
00299   else
00300     return 0;
00301   #endif
00302 
00303   return 0;
00304 }
00305 
00306 // c = alpha*a*a' + beta*c  for op = 'N'or'n'
00307 // c = alpha*a'*a + beta*c  for op = 'T'or't','C'or'c'
00308 int EIGEN_BLAS_FUNC(syrk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
00309 {
00310 //   std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
00311   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
00312   static functype func[8];
00313 
00314   static bool init = false;
00315   if(!init)
00316   {
00317     for(int k=0; k<8; ++k)
00318       func[k] = 0;
00319 
00320     func[NOTR  | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run);
00321     func[TR    | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run);
00322     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run);
00323 
00324     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run);
00325     func[TR    | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run);
00326     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run);
00327 
00328     init = true;
00329   }
00330 
00331   Scalar* a = reinterpret_cast<Scalar*>(pa);
00332   Scalar* c = reinterpret_cast<Scalar*>(pc);
00333   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00334   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
00335 
00336   int info = 0;
00337   if(UPLO(*uplo)==INVALID)                                            info = 1;
00338   else if(OP(*op)==INVALID)                                           info = 2;
00339   else if(*n<0)                                                       info = 3;
00340   else if(*k<0)                                                       info = 4;
00341   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
00342   else if(*ldc<std::max(1,*n))                                        info = 10;
00343   if(info)
00344     return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);
00345 
00346   if(beta!=Scalar(1))
00347   {
00348     if(UPLO(*uplo)==UP)
00349       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
00350       else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
00351     else
00352       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
00353       else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
00354   }
00355 
00356   #if ISCOMPLEX
00357   // FIXME add support for symmetric complex matrix
00358   if(UPLO(*uplo)==UP)
00359   {
00360     if(OP(*op)==NOTR)
00361       matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
00362     else
00363       matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
00364   }
00365   else
00366   {
00367     if(OP(*op)==NOTR)
00368       matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
00369     else
00370       matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
00371   }
00372   #else
00373   int code = OP(*op) | (UPLO(*uplo) << 2);
00374   func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
00375   #endif
00376 
00377   return 0;
00378 }
00379 
00380 // c = alpha*a*b' + alpha*b*a' + beta*c  for op = 'N'or'n'
00381 // c = alpha*a'*b + alpha*b'*a + beta*c  for op = 'T'or't'
00382 int EIGEN_BLAS_FUNC(syr2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
00383 {
00384   Scalar* a = reinterpret_cast<Scalar*>(pa);
00385   Scalar* b = reinterpret_cast<Scalar*>(pb);
00386   Scalar* c = reinterpret_cast<Scalar*>(pc);
00387   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00388   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
00389 
00390   int info = 0;
00391   if(UPLO(*uplo)==INVALID)                                            info = 1;
00392   else if(OP(*op)==INVALID)                                           info = 2;
00393   else if(*n<0)                                                       info = 3;
00394   else if(*k<0)                                                       info = 4;
00395   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
00396   else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
00397   else if(*ldc<std::max(1,*n))                                        info = 12;
00398   if(info)
00399     return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6);
00400 
00401   if(beta!=Scalar(1))
00402   {
00403     if(UPLO(*uplo)==UP)
00404       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
00405       else                matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
00406     else
00407       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
00408       else                matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
00409   }
00410 
00411   if(*k==0)
00412     return 1;
00413 
00414   if(OP(*op)==NOTR)
00415   {
00416     if(UPLO(*uplo)==UP)
00417     {
00418       matrix(c, *n, *n, *ldc).triangularView<Upper>()
00419         += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
00420         +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
00421     }
00422     else if(UPLO(*uplo)==LO)
00423       matrix(c, *n, *n, *ldc).triangularView<Lower>()
00424         += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
00425         +  alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
00426   }
00427   else if(OP(*op)==TR || OP(*op)==ADJ)
00428   {
00429     if(UPLO(*uplo)==UP)
00430       matrix(c, *n, *n, *ldc).triangularView<Upper>()
00431         += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
00432         +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
00433     else if(UPLO(*uplo)==LO)
00434       matrix(c, *n, *n, *ldc).triangularView<Lower>()
00435         += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
00436         +  alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
00437   }
00438 
00439   return 0;
00440 }
00441 
00442 
00443 #if ISCOMPLEX
00444 
00445 // c = alpha*a*b + beta*c  for side = 'L'or'l'
00446 // c = alpha*b*a + beta*c  for side = 'R'or'r
00447 int EIGEN_BLAS_FUNC(hemm)(char *side, char *uplo, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
00448 {
00449   Scalar* a = reinterpret_cast<Scalar*>(pa);
00450   Scalar* b = reinterpret_cast<Scalar*>(pb);
00451   Scalar* c = reinterpret_cast<Scalar*>(pc);
00452   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00453   Scalar beta  = *reinterpret_cast<Scalar*>(pbeta);
00454 
00455 //   std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
00456 
00457   int info = 0;
00458   if(SIDE(*side)==INVALID)                                            info = 1;
00459   else if(UPLO(*uplo)==INVALID)                                       info = 2;
00460   else if(*m<0)                                                       info = 3;
00461   else if(*n<0)                                                       info = 4;
00462   else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n))                 info = 7;
00463   else if(*ldb<std::max(1,*m))                                        info = 9;
00464   else if(*ldc<std::max(1,*m))                                        info = 12;
00465   if(info)
00466     return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6);
00467 
00468   if(beta==Scalar(0))       matrix(c, *m, *n, *ldc).setZero();
00469   else if(beta!=Scalar(1))  matrix(c, *m, *n, *ldc) *= beta;
00470 
00471   if(*m==0 || *n==0)
00472   {
00473     return 1;
00474   }
00475 
00476   if(SIDE(*side)==LEFT)
00477   {
00478     if(UPLO(*uplo)==UP)       internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj,  ColMajor,false,false, ColMajor>
00479                                 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
00480     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor>
00481                                 ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
00482     else                      return 0;
00483   }
00484   else if(SIDE(*side)==RIGHT)
00485   {
00486     if(UPLO(*uplo)==UP)       matrix(c,*m,*n,*ldc) += alpha * matrix(b,*m,*n,*ldb) * matrix(a,*n,*n,*lda).selfadjointView<Upper>();/*internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, RowMajor,true,Conj,  ColMajor>
00487                                 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/
00488     else if(UPLO(*uplo)==LO)  internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor>
00489                                 ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
00490     else                      return 0;
00491   }
00492   else
00493   {
00494     return 0;
00495   }
00496 
00497   return 0;
00498 }
00499 
00500 // c = alpha*a*conj(a') + beta*c  for op = 'N'or'n'
00501 // c = alpha*conj(a')*a + beta*c  for op  = 'C'or'c'
00502 int EIGEN_BLAS_FUNC(herk)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pbeta, RealScalar *pc, int *ldc)
00503 {
00504   typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar);
00505   static functype func[8];
00506 
00507   static bool init = false;
00508   if(!init)
00509   {
00510     for(int k=0; k<8; ++k)
00511       func[k] = 0;
00512 
00513     func[NOTR  | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run);
00514     func[ADJ   | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run);
00515 
00516     func[NOTR  | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run);
00517     func[ADJ   | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run);
00518 
00519     init = true;
00520   }
00521 
00522   Scalar* a = reinterpret_cast<Scalar*>(pa);
00523   Scalar* c = reinterpret_cast<Scalar*>(pc);
00524   RealScalar alpha = *palpha;
00525   RealScalar beta  = *pbeta;
00526 
00527 //   std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
00528 
00529   int info = 0;
00530   if(UPLO(*uplo)==INVALID)                                            info = 1;
00531   else if((OP(*op)==INVALID) || (OP(*op)==TR))                        info = 2;
00532   else if(*n<0)                                                       info = 3;
00533   else if(*k<0)                                                       info = 4;
00534   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
00535   else if(*ldc<std::max(1,*n))                                        info = 10;
00536   if(info)
00537     return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6);
00538 
00539   int code = OP(*op) | (UPLO(*uplo) << 2);
00540 
00541   if(beta!=RealScalar(1))
00542   {
00543     if(UPLO(*uplo)==UP)
00544       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
00545       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
00546     else
00547       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
00548       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
00549   
00550     if(beta!=Scalar(0))
00551     {
00552       matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
00553       matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
00554     }
00555   }
00556 
00557   if(*k>0 && alpha!=RealScalar(0))
00558   {
00559     func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
00560     matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
00561   }
00562   return 0;
00563 }
00564 
00565 // c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c,  for op = 'N'or'n'
00566 // c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c,  for op = 'C'or'c'
00567 int EIGEN_BLAS_FUNC(her2k)(char *uplo, char *op, int *n, int *k, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, RealScalar *pbeta, RealScalar *pc, int *ldc)
00568 {
00569   Scalar* a = reinterpret_cast<Scalar*>(pa);
00570   Scalar* b = reinterpret_cast<Scalar*>(pb);
00571   Scalar* c = reinterpret_cast<Scalar*>(pc);
00572   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00573   RealScalar beta  = *pbeta;
00574 
00575   int info = 0;
00576   if(UPLO(*uplo)==INVALID)                                            info = 1;
00577   else if((OP(*op)==INVALID) || (OP(*op)==TR))                        info = 2;
00578   else if(*n<0)                                                       info = 3;
00579   else if(*k<0)                                                       info = 4;
00580   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 7;
00581   else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k))                     info = 9;
00582   else if(*ldc<std::max(1,*n))                                        info = 12;
00583   if(info)
00584     return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6);
00585 
00586   if(beta!=RealScalar(1))
00587   {
00588     if(UPLO(*uplo)==UP)
00589       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
00590       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
00591     else
00592       if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
00593       else                matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
00594 
00595     if(beta!=Scalar(0))
00596     {
00597       matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
00598       matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
00599     }
00600   }
00601   else if(*k>0 && alpha!=Scalar(0))
00602     matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
00603 
00604   if(*k==0)
00605     return 1;
00606 
00607   if(OP(*op)==NOTR)
00608   {
00609     if(UPLO(*uplo)==UP)
00610     {
00611       matrix(c, *n, *n, *ldc).triangularView<Upper>()
00612         +=         alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
00613         +  internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
00614     }
00615     else if(UPLO(*uplo)==LO)
00616       matrix(c, *n, *n, *ldc).triangularView<Lower>()
00617         += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
00618         +  internal::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
00619   }
00620   else if(OP(*op)==ADJ)
00621   {
00622     if(UPLO(*uplo)==UP)
00623       matrix(c, *n, *n, *ldc).triangularView<Upper>()
00624         += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
00625         +  internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
00626     else if(UPLO(*uplo)==LO)
00627       matrix(c, *n, *n, *ldc).triangularView<Lower>()
00628         += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
00629         +  internal::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
00630   }
00631 
00632   return 1;
00633 }
00634 
00635 #endif // ISCOMPLEX


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:56