Newer
Older
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
int EIGEN_BLAS_FUNC(gemm)(const char *opa, const char *opb, const int *m, const int *n, const int *k, const RealScalar *palpha,
const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
// std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*);
static const functype func[12] = {
// array index: NOTR | (NOTR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run),
// array index: TR | (NOTR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run),
// array index: ADJ | (NOTR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run),
0,
// array index: NOTR | (TR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run),
// array index: TR | (TR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run),
// array index: ADJ | (TR << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run),
0,
// array index: NOTR | (ADJ << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run),
// array index: TR | (ADJ << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run),
// array index: ADJ | (ADJ << 2)
(internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run),
0
};
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
const Scalar* b = reinterpret_cast<const Scalar*>(pb);
Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
int info = 0;
if(OP(*opa)==INVALID) info = 1;
else if(OP(*opb)==INVALID) info = 2;
else if(*m<0) info = 3;
else if(*n<0) info = 4;
else if(*k<0) info = 5;
else if(*lda<std::max(1,(OP(*opa)==NOTR)?*m:*k)) info = 8;
else if(*ldb<std::max(1,(OP(*opb)==NOTR)?*k:*n)) info = 10;
else if(*ldc<std::max(1,*m)) info = 13;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"GEMM ",&info,6);
if(beta!=Scalar(1))
{
if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
else matrix(c, *m, *n, *ldc) *= beta;
}
if(*k == 0)
return 0;
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k,1,true);
int code = OP(*opa) | (OP(*opb) << 2);
func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0);
return 0;
}
int EIGEN_BLAS_FUNC(trsm)(const char *side, const char *uplo, const char *opa, const char *diag, const int *m, const int *n,
const RealScalar *palpha, const RealScalar *pa, const int *lda, RealScalar *pb, const int *ldb)
{
// std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking<Scalar,Scalar>&);
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
static const functype func[32] = {
// array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,ColMajor,ColMajor>::run),
// array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,RowMajor,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, Conj, RowMajor,ColMajor>::run),\
0,
// array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,ColMajor,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,RowMajor,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,ColMajor,ColMajor>::run),
// array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,RowMajor,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,ColMajor,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,RowMajor,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run),
// array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run),
// array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run),
0
};
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
int info = 0;
if(SIDE(*side)==INVALID) info = 1;
else if(UPLO(*uplo)==INVALID) info = 2;
else if(OP(*opa)==INVALID) info = 3;
else if(DIAG(*diag)==INVALID) info = 4;
else if(*m<0) info = 5;
else if(*n<0) info = 6;
else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9;
else if(*ldb<std::max(1,*m)) info = 11;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6);
int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m,1,false);
func[code](*m, *n, a, *lda, b, *ldb, blocking);
}
else
{
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n,1,false);
func[code](*n, *m, a, *lda, b, *ldb, blocking);
}
if(alpha!=Scalar(1))
matrix(b,*m,*n,*ldb) *= alpha;
return 0;
}
// b = alpha*op(a)*b for side = 'L'or'l'
// b = alpha*b*op(a) for side = 'R'or'r'
int EIGEN_BLAS_FUNC(trmm)(const char *side, const char *uplo, const char *opa, const char *diag, const int *m, const int *n,
const RealScalar *palpha, const RealScalar *pa, const int *lda, RealScalar *pb, const int *ldb)
{
// std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n";
typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
static const functype func[32] = {
// array index: NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,false,ColMajor,false,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,false,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,false,ColMajor,false,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,false,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
0,
// array index: NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run),
// array index: ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run),
0,
// array index: NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run),
// array index: TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run),
// array index: ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)
(internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run),
0
};
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
int info = 0;
if(SIDE(*side)==INVALID) info = 1;
else if(UPLO(*uplo)==INVALID) info = 2;
else if(OP(*opa)==INVALID) info = 3;
else if(DIAG(*diag)==INVALID) info = 4;
else if(*m<0) info = 5;
else if(*n<0) info = 6;
else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9;
else if(*ldb<std::max(1,*m)) info = 11;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"TRMM ",&info,6);
int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
if(*m==0 || *n==0)
return 1;
// FIXME find a way to avoid this copy
Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb);
matrix(b,*m,*n,*ldb).setZero();
if(SIDE(*side)==LEFT)
{
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m,1,false);
func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha, blocking);
}
else
{
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n,1,false);
func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha, blocking);
}
return 1;
}
// c = alpha*a*b + beta*c for side = 'L'or'l'
// c = alpha*b*a + beta*c for side = 'R'or'r
int EIGEN_BLAS_FUNC(symm)(const char *side, const char *uplo, const int *m, const int *n, const RealScalar *palpha,
const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
// std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n";
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
const Scalar* b = reinterpret_cast<const Scalar*>(pb);
Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
int info = 0;
if(SIDE(*side)==INVALID) info = 1;
else if(UPLO(*uplo)==INVALID) info = 2;
else if(*m<0) info = 3;
else if(*n<0) info = 4;
else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7;
else if(*ldb<std::max(1,*m)) info = 9;
else if(*ldc<std::max(1,*m)) info = 12;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"SYMM ",&info,6);
if(beta!=Scalar(1))
{
if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
else matrix(c, *m, *n, *ldc) *= beta;
}
if(*m==0 || *n==0)
{
return 1;
}
#if ISCOMPLEX
// FIXME add support for symmetric complex matrix
Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size);
if(UPLO(*uplo)==UP)
{
matA.triangularView<Upper>() = matrix(a,size,size,*lda);
matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose();
}
else if(UPLO(*uplo)==LO)
{
matA.triangularView<Lower>() = matrix(a,size,size,*lda);
matA.triangularView<Upper>() = matrix(a,size,size,*lda).transpose();
}
if(SIDE(*side)==LEFT)
matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb);
else if(SIDE(*side)==RIGHT)
matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA;
#else
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,size,1,false);
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, blocking);
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, blocking);
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, blocking);
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, blocking);
else return 0;
else
return 0;
#endif
return 0;
}
// c = alpha*a*a' + beta*c for op = 'N'or'n'
// c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c'
int EIGEN_BLAS_FUNC(syrk)(const char *uplo, const char *op, const int *n, const int *k,
const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
{
// std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
#if !ISCOMPLEX
typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
static const functype func[8] = {
// array index: NOTR | (UP << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run),
// array index: TR | (UP << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run),
// array index: ADJ | (UP << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run),
0,
// array index: NOTR | (LO << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run),
// array index: TR | (LO << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run),
// array index: ADJ | (LO << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run),
0
};
Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
else if(OP(*op)==INVALID || (ISCOMPLEX && OP(*op)==ADJ) ) info = 2;
else if(*n<0) info = 3;
else if(*k<0) info = 4;
else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
else if(*ldc<std::max(1,*n)) info = 10;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);
if(beta!=Scalar(1))
{
if(UPLO(*uplo)==UP)
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
else
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
}
#if ISCOMPLEX
// FIXME add support for symmetric complex matrix
if(UPLO(*uplo)==UP)
{
if(OP(*op)==NOTR)
matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
else
matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
}
else
{
if(OP(*op)==NOTR)
matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
else
matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
}
#else
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*n,*n,*k,1,false);
func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha, blocking);
#endif
return 0;
}
// c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n'
// c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't'
int EIGEN_BLAS_FUNC(syr2k)(const char *uplo, const char *op, const int *n, const int *k, const RealScalar *palpha,
const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
const Scalar* b = reinterpret_cast<const Scalar*>(pb);
Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
// std::cerr << "in syr2k " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << *ldb << " " << beta << " " << *ldc << "\n";
else if(OP(*op)==INVALID || (ISCOMPLEX && OP(*op)==ADJ) ) info = 2;
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
else if(*n<0) info = 3;
else if(*k<0) info = 4;
else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 9;
else if(*ldc<std::max(1,*n)) info = 12;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6);
if(beta!=Scalar(1))
{
if(UPLO(*uplo)==UP)
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
else
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
}
if(*k==0)
return 1;
if(OP(*op)==NOTR)
{
if(UPLO(*uplo)==UP)
{
matrix(c, *n, *n, *ldc).triangularView<Upper>()
+= alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
+ alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
}
else if(UPLO(*uplo)==LO)
matrix(c, *n, *n, *ldc).triangularView<Lower>()
+= alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
+ alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
}
else if(OP(*op)==TR || OP(*op)==ADJ)
{
if(UPLO(*uplo)==UP)
matrix(c, *n, *n, *ldc).triangularView<Upper>()
+= alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
+ alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
else if(UPLO(*uplo)==LO)
matrix(c, *n, *n, *ldc).triangularView<Lower>()
+= alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
+ alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
}
return 0;
}
#if ISCOMPLEX
// c = alpha*a*b + beta*c for side = 'L'or'l'
// c = alpha*b*a + beta*c for side = 'R'or'r
int EIGEN_BLAS_FUNC(hemm)(const char *side, const char *uplo, const int *m, const int *n, const RealScalar *palpha,
const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
const Scalar* b = reinterpret_cast<const Scalar*>(pb);
Scalar alpha = *reinterpret_cast<const Scalar*>(palpha);
Scalar beta = *reinterpret_cast<const Scalar*>(pbeta);
// std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
int info = 0;
if(SIDE(*side)==INVALID) info = 1;
else if(UPLO(*uplo)==INVALID) info = 2;
else if(*m<0) info = 3;
else if(*n<0) info = 4;
else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7;
else if(*ldb<std::max(1,*m)) info = 9;
else if(*ldc<std::max(1,*m)) info = 12;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6);
if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
else if(beta!=Scalar(1)) matrix(c, *m, *n, *ldc) *= beta;
if(*m==0 || *n==0)
{
return 1;
}
int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,size,1,false);
if(SIDE(*side)==LEFT)
{
if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj, ColMajor,false,false, ColMajor>
::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha, blocking);
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, blocking);
else return 0;
}
else if(SIDE(*side)==RIGHT)
{
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>
::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha, blocking);*/
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, blocking);
else return 0;
}
else
{
return 0;
}
return 0;
}
// c = alpha*a*conj(a') + beta*c for op = 'N'or'n'
// c = alpha*conj(a')*a + beta*c for op = 'C'or'c'
int EIGEN_BLAS_FUNC(herk)(const char *uplo, const char *op, const int *n, const int *k,
const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
// std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
static const functype func[8] = {
// array index: NOTR | (UP << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run),
0,
// array index: ADJ | (UP << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run),
0,
// array index: NOTR | (LO << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run),
0,
// array index: ADJ | (LO << 2)
(internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run),
0
};
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
Scalar* c = reinterpret_cast<Scalar*>(pc);
RealScalar alpha = *palpha;
RealScalar beta = *pbeta;
// std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
int info = 0;
if(UPLO(*uplo)==INVALID) info = 1;
else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2;
else if(*n<0) info = 3;
else if(*k<0) info = 4;
else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
else if(*ldc<std::max(1,*n)) info = 10;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6);
int code = OP(*op) | (UPLO(*uplo) << 2);
if(beta!=RealScalar(1))
{
if(UPLO(*uplo)==UP)
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
else
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
if(beta!=Scalar(0))
{
matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
}
}
if(*k>0 && alpha!=RealScalar(0))
{
internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*n,*n,*k,1,false);
func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha, blocking);
matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
}
return 0;
}
// c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n'
// c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c, for op = 'C'or'c'
int EIGEN_BLAS_FUNC(her2k)(const char *uplo, const char *op, const int *n, const int *k,
const RealScalar *palpha, const RealScalar *pa, const int *lda, const RealScalar *pb, const int *ldb, const RealScalar *pbeta, RealScalar *pc, const int *ldc)
const Scalar* a = reinterpret_cast<const Scalar*>(pa);
const Scalar* b = reinterpret_cast<const Scalar*>(pb);
// std::cerr << "in her2k " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << *ldb << " " << beta << " " << *ldc << "\n";
int info = 0;
if(UPLO(*uplo)==INVALID) info = 1;
else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2;
else if(*n<0) info = 3;
else if(*k<0) info = 4;
else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
else if(*ldc<std::max(1,*n)) info = 12;
if(info)
return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6);
if(beta!=RealScalar(1))
{
if(UPLO(*uplo)==UP)
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
else
if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
if(beta!=Scalar(0))
{
matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
}
}
else if(*k>0 && alpha!=Scalar(0))
matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
if(*k==0)
return 1;
if(OP(*op)==NOTR)
{
if(UPLO(*uplo)==UP)
{
matrix(c, *n, *n, *ldc).triangularView<Upper>()
+= alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
+ numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
}
else if(UPLO(*uplo)==LO)
matrix(c, *n, *n, *ldc).triangularView<Lower>()
+= alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
+ numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
}
else if(OP(*op)==ADJ)
{
if(UPLO(*uplo)==UP)
matrix(c, *n, *n, *ldc).triangularView<Upper>()
+= alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
+ numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
else if(UPLO(*uplo)==LO)
matrix(c, *n, *n, *ldc).triangularView<Lower>()
+= alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
+ numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
}
return 1;
}
#endif // ISCOMPLEX