re_vision: level2_impl.h Source File

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(gemv)(char *opa, int *m, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *pb, int *incb, RealScalar *pbeta, RealScalar *pc, int *incc)
00028 {
00029   typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int , Scalar *, int, Scalar);
00030   static functype func[4];
00031 
00032   static bool init = false;
00033   if(!init)
00034   {
00035     for(int k=0; k<4; ++k)
00036       func[k] = 0;
00037 
00038     func[NOTR] = (internal::general_matrix_vector_product<int,Scalar,ColMajor,false,Scalar,false>::run);
00039     func[TR  ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,false,Scalar,false>::run);
00040     func[ADJ ] = (internal::general_matrix_vector_product<int,Scalar,RowMajor,Conj, Scalar,false>::run);
00041 
00042     init = true;
00043   }
00044 
00045   Scalar* a = reinterpret_cast<Scalar*>(pa);
00046   Scalar* b = reinterpret_cast<Scalar*>(pb);
00047   Scalar* c = reinterpret_cast<Scalar*>(pc);
00048   Scalar alpha  = *reinterpret_cast<Scalar*>(palpha);
00049   Scalar beta   = *reinterpret_cast<Scalar*>(pbeta);
00050 
00051   // check arguments
00052   int info = 0;
00053   if(OP(*opa)==INVALID)           info = 1;
00054   else if(*m<0)                   info = 2;
00055   else if(*n<0)                   info = 3;
00056   else if(*lda<std::max(1,*m))    info = 6;
00057   else if(*incb==0)               info = 8;
00058   else if(*incc==0)               info = 11;
00059   if(info)
00060     return xerbla_(SCALAR_SUFFIX_UP"GEMV ",&info,6);
00061 
00062   if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
00063     return 0;
00064 
00065   int actual_m = *m;
00066   int actual_n = *n;
00067   if(OP(*opa)!=NOTR)
00068     std::swap(actual_m,actual_n);
00069 
00070   Scalar* actual_b = get_compact_vector(b,actual_n,*incb);
00071   Scalar* actual_c = get_compact_vector(c,actual_m,*incc);
00072 
00073   if(beta!=Scalar(1))
00074   {
00075     if(beta==Scalar(0)) vector(actual_c, actual_m).setZero();
00076     else                vector(actual_c, actual_m) *= beta;
00077   }
00078 
00079   int code = OP(*opa);
00080   func[code](actual_m, actual_n, a, *lda, actual_b, 1, actual_c, 1, alpha);
00081 
00082   if(actual_b!=b) delete[] actual_b;
00083   if(actual_c!=c) delete[] copy_back(actual_c,c,actual_m,*incc);
00084 
00085   return 1;
00086 }
00087 
00088 int EIGEN_BLAS_FUNC(trsv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
00089 {
00090   typedef void (*functype)(int, const Scalar *, int, Scalar *);
00091   static functype func[16];
00092 
00093   static bool init = false;
00094   if(!init)
00095   {
00096     for(int k=0; k<16; ++k)
00097       func[k] = 0;
00098 
00099     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,ColMajor>::run);
00100     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,RowMajor>::run);
00101     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       Conj, RowMajor>::run);
00102 
00103     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|0,       false,ColMajor>::run);
00104     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       false,RowMajor>::run);
00105     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|0,       Conj, RowMajor>::run);
00106 
00107     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,ColMajor>::run);
00108     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,RowMajor>::run);
00109     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,Conj, RowMajor>::run);
00110 
00111     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Lower|UnitDiag,false,ColMajor>::run);
00112     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,false,RowMajor>::run);
00113     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::triangular_solve_vector<Scalar,Scalar,int,OnTheLeft, Upper|UnitDiag,Conj, RowMajor>::run);
00114 
00115     init = true;
00116   }
00117 
00118   Scalar* a = reinterpret_cast<Scalar*>(pa);
00119   Scalar* b = reinterpret_cast<Scalar*>(pb);
00120 
00121   int info = 0;
00122   if(UPLO(*uplo)==INVALID)                                            info = 1;
00123   else if(OP(*opa)==INVALID)                                          info = 2;
00124   else if(DIAG(*diag)==INVALID)                                       info = 3;
00125   else if(*n<0)                                                       info = 4;
00126   else if(*lda<std::max(1,*n))                                        info = 6;
00127   else if(*incb==0)                                                   info = 8;
00128   if(info)
00129     return xerbla_(SCALAR_SUFFIX_UP"TRSV ",&info,6);
00130 
00131   Scalar* actual_b = get_compact_vector(b,*n,*incb);
00132 
00133   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
00134   func[code](*n, a, *lda, actual_b);
00135 
00136   if(actual_b!=b) delete[] copy_back(actual_b,b,*n,*incb);
00137 
00138   return 0;
00139 }
00140 
00141 
00142 
00143 int EIGEN_BLAS_FUNC(trmv)(char *uplo, char *opa, char *diag, int *n, RealScalar *pa, int *lda, RealScalar *pb, int *incb)
00144 {
00145   typedef void (*functype)(int, int, const Scalar *, int, const Scalar *, int, Scalar *, int, Scalar);
00146   static functype func[16];
00147 
00148   static bool init = false;
00149   if(!init)
00150   {
00151     for(int k=0; k<16; ++k)
00152       func[k] = 0;
00153 
00154     func[NOTR  | (UP << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Upper|0,       Scalar,false,Scalar,false,ColMajor>::run);
00155     func[TR    | (UP << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Lower|0,       Scalar,false,Scalar,false,RowMajor>::run);
00156     func[ADJ   | (UP << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Lower|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
00157 
00158     func[NOTR  | (LO << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Lower|0,       Scalar,false,Scalar,false,ColMajor>::run);
00159     func[TR    | (LO << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Upper|0,       Scalar,false,Scalar,false,RowMajor>::run);
00160     func[ADJ   | (LO << 2) | (NUNIT << 3)] = (internal::product_triangular_matrix_vector<int,Upper|0,       Scalar,Conj, Scalar,false,RowMajor>::run);
00161 
00162     func[NOTR  | (UP << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Upper|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
00163     func[TR    | (UP << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Lower|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
00164     func[ADJ   | (UP << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Lower|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
00165 
00166     func[NOTR  | (LO << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Lower|UnitDiag,Scalar,false,Scalar,false,ColMajor>::run);
00167     func[TR    | (LO << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Upper|UnitDiag,Scalar,false,Scalar,false,RowMajor>::run);
00168     func[ADJ   | (LO << 2) | (UNIT  << 3)] = (internal::product_triangular_matrix_vector<int,Upper|UnitDiag,Scalar,Conj, Scalar,false,RowMajor>::run);
00169 
00170     init = true;
00171   }
00172 
00173   Scalar* a = reinterpret_cast<Scalar*>(pa);
00174   Scalar* b = reinterpret_cast<Scalar*>(pb);
00175 
00176   int info = 0;
00177   if(UPLO(*uplo)==INVALID)                                            info = 1;
00178   else if(OP(*opa)==INVALID)                                          info = 2;
00179   else if(DIAG(*diag)==INVALID)                                       info = 3;
00180   else if(*n<0)                                                       info = 4;
00181   else if(*lda<std::max(1,*n))                                        info = 6;
00182   else if(*incb==0)                                                   info = 8;
00183   if(info)
00184     return xerbla_(SCALAR_SUFFIX_UP"TRMV ",&info,6);
00185 
00186   if(*n==0)
00187     return 1;
00188 
00189   Scalar* actual_b = get_compact_vector(b,*n,*incb);
00190   Matrix<Scalar,Dynamic,1> res(*n);
00191   res.setZero();
00192 
00193   int code = OP(*opa) | (UPLO(*uplo) << 2) | (DIAG(*diag) << 3);
00194   if(code>=16 || func[code]==0)
00195     return 0;
00196 
00197   func[code](*n, *n, a, *lda, actual_b, 1, res.data(), 1, Scalar(1));
00198 
00199   copy_back(res.data(),b,*n,*incb);
00200   if(actual_b!=b) delete[] actual_b;
00201 
00202   return 0;
00203 }
00204 
00212 int EIGEN_BLAS_FUNC(gbmv)(char *trans, int *m, int *n, int *kl, int *ku, RealScalar *palpha, RealScalar *pa, int *lda,
00213                           RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy)
00214 {
00215   Scalar* a = reinterpret_cast<Scalar*>(pa);
00216   Scalar* x = reinterpret_cast<Scalar*>(px);
00217   Scalar* y = reinterpret_cast<Scalar*>(py);
00218   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00219   Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
00220   int coeff_rows = *kl+*ku+1;
00221   
00222   int info = 0;
00223        if(OP(*trans)==INVALID)                                        info = 1;
00224   else if(*m<0)                                                       info = 2;
00225   else if(*n<0)                                                       info = 3;
00226   else if(*kl<0)                                                      info = 4;
00227   else if(*ku<0)                                                      info = 5;
00228   else if(*lda<coeff_rows)                                            info = 8;
00229   else if(*incx==0)                                                   info = 10;
00230   else if(*incy==0)                                                   info = 13;
00231   if(info)
00232     return xerbla_(SCALAR_SUFFIX_UP"GBMV ",&info,6);
00233   
00234   if(*m==0 || *n==0 || (alpha==Scalar(0) && beta==Scalar(1)))
00235     return 0;
00236   
00237   int actual_m = *m;
00238   int actual_n = *n;
00239   if(OP(*trans)!=NOTR)
00240     std::swap(actual_m,actual_n);
00241   
00242   Scalar* actual_x = get_compact_vector(x,actual_n,*incx);
00243   Scalar* actual_y = get_compact_vector(y,actual_m,*incy);
00244   
00245   if(beta!=Scalar(1))
00246   {
00247     if(beta==Scalar(0)) vector(actual_y, actual_m).setZero();
00248     else                vector(actual_y, actual_m) *= beta;
00249   }
00250   
00251   MatrixType mat_coeffs(a,coeff_rows,*n,*lda);
00252   
00253   int nb = std::min(*n,(*m)+(*ku));
00254   for(int j=0; j<nb; ++j)
00255   {
00256     int start = std::max(0,j - *ku);
00257     int end = std::min((*m)-1,j + *kl);
00258     int len = end - start + 1;
00259     int offset = (*ku) - j + start;
00260     if(OP(*trans)==NOTR)
00261       vector(actual_y+start,len) += (alpha*actual_x[j]) * mat_coeffs.col(j).segment(offset,len);
00262     else if(OP(*trans)==TR)
00263       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).transpose() * vector(actual_x+start,len) ).value();
00264     else
00265       actual_y[j] += alpha * ( mat_coeffs.col(j).segment(offset,len).adjoint()   * vector(actual_x+start,len) ).value();
00266   }    
00267   
00268   if(actual_x!=x) delete[] actual_x;
00269   if(actual_y!=y) delete[] copy_back(actual_y,y,actual_m,*incy);
00270   
00271   return 0;
00272 }
00273 
00281 // int EIGEN_BLAS_FUNC(tbmv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx)
00282 // {
00283 //   return 1;
00284 // }
00285 
00297 // int EIGEN_BLAS_FUNC(tbsv)(char *uplo, char *trans, char *diag, int *n, int *k, RealScalar *a, int *lda, RealScalar *x, int *incx)
00298 // {
00299 //   return 1;
00300 // }
00301 
00309 // int EIGEN_BLAS_FUNC(tpmv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx)
00310 // {
00311 //   return 1;
00312 // }
00313 
00324 // int EIGEN_BLAS_FUNC(tpsv)(char *uplo, char *trans, char *diag, int *n, RealScalar *ap, RealScalar *x, int *incx)
00325 // {
00326 //   return 1;
00327 // }
00328 
00336 int EIGEN_BLAS_FUNC(ger)(int *m, int *n, Scalar *palpha, Scalar *px, int *incx, Scalar *py, int *incy, Scalar *pa, int *lda)
00337 {
00338   Scalar* x = reinterpret_cast<Scalar*>(px);
00339   Scalar* y = reinterpret_cast<Scalar*>(py);
00340   Scalar* a = reinterpret_cast<Scalar*>(pa);
00341   Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
00342 
00343   int info = 0;
00344        if(*m<0)                                                       info = 1;
00345   else if(*n<0)                                                       info = 2;
00346   else if(*incx==0)                                                   info = 5;
00347   else if(*incy==0)                                                   info = 7;
00348   else if(*lda<std::max(1,*m))                                        info = 9;
00349   if(info)
00350     return xerbla_(SCALAR_SUFFIX_UP"GER  ",&info,6);
00351 
00352   if(alpha==Scalar(0))
00353     return 1;
00354 
00355   Scalar* x_cpy = get_compact_vector(x,*m,*incx);
00356   Scalar* y_cpy = get_compact_vector(y,*n,*incy);
00357 
00358   // TODO perform direct calls to underlying implementation
00359   matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).adjoint();
00360 
00361   if(x_cpy!=x)  delete[] x_cpy;
00362   if(y_cpy!=y)  delete[] y_cpy;
00363 
00364   return 1;
00365 }
00366 
00367