00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
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
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
00160
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
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
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
00240
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
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
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
00307
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
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
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
00381
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
00446
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
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>();
00487
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
00501
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
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
00566
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