[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2009-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
| 5 | //
|
---|
| 6 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
| 7 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
| 8 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
| 9 |
|
---|
| 10 | #include "common.h"
|
---|
| 11 |
|
---|
| 12 | 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)
|
---|
| 13 | {
|
---|
| 14 | // std::cerr << "in gemm " << *opa << " " << *opb << " " << *m << " " << *n << " " << *k << " " << *lda << " " << *ldb << " " << *ldc << " " << *palpha << " " << *pbeta << "\n";
|
---|
| 15 | typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, Scalar, internal::level3_blocking<Scalar,Scalar>&, Eigen::internal::GemmParallelInfo<DenseIndex>*);
|
---|
| 16 | static functype func[12];
|
---|
| 17 |
|
---|
| 18 | static bool init = false;
|
---|
| 19 | if(!init)
|
---|
| 20 | {
|
---|
| 21 | for(int k=0; k<12; ++k)
|
---|
| 22 | func[k] = 0;
|
---|
| 23 | func[NOTR | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,ColMajor,false,ColMajor>::run);
|
---|
| 24 | func[TR | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,false,ColMajor>::run);
|
---|
| 25 | func[ADJ | (NOTR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor>::run);
|
---|
| 26 | func[NOTR | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,false,ColMajor>::run);
|
---|
| 27 | func[TR | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,false,ColMajor>::run);
|
---|
| 28 | func[ADJ | (TR << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,false,ColMajor>::run);
|
---|
| 29 | func[NOTR | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
|
---|
| 30 | func[TR | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,false,Scalar,RowMajor,Conj, ColMajor>::run);
|
---|
| 31 | func[ADJ | (ADJ << 2)] = (internal::general_matrix_matrix_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,RowMajor,Conj, ColMajor>::run);
|
---|
| 32 | init = true;
|
---|
| 33 | }
|
---|
| 34 |
|
---|
| 35 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 36 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 37 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 38 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 39 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
| 40 |
|
---|
| 41 | int info = 0;
|
---|
| 42 | if(OP(*opa)==INVALID) info = 1;
|
---|
| 43 | else if(OP(*opb)==INVALID) info = 2;
|
---|
| 44 | else if(*m<0) info = 3;
|
---|
| 45 | else if(*n<0) info = 4;
|
---|
| 46 | else if(*k<0) info = 5;
|
---|
| 47 | else if(*lda<std::max(1,(OP(*opa)==NOTR)?*m:*k)) info = 8;
|
---|
| 48 | else if(*ldb<std::max(1,(OP(*opb)==NOTR)?*k:*n)) info = 10;
|
---|
| 49 | else if(*ldc<std::max(1,*m)) info = 13;
|
---|
| 50 | if(info)
|
---|
| 51 | return xerbla_(SCALAR_SUFFIX_UP"GEMM ",&info,6);
|
---|
| 52 |
|
---|
| 53 | if(beta!=Scalar(1))
|
---|
| 54 | {
|
---|
| 55 | if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
|
---|
| 56 | else matrix(c, *m, *n, *ldc) *= beta;
|
---|
| 57 | }
|
---|
| 58 |
|
---|
| 59 | internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic> blocking(*m,*n,*k);
|
---|
| 60 |
|
---|
| 61 | int code = OP(*opa) | (OP(*opb) << 2);
|
---|
| 62 | func[code](*m, *n, *k, a, *lda, b, *ldb, c, *ldc, alpha, blocking, 0);
|
---|
| 63 | return 0;
|
---|
| 64 | }
|
---|
| 65 |
|
---|
| 66 | 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)
|
---|
| 67 | {
|
---|
| 68 | // std::cerr << "in trsm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << "," << *n << " " << *palpha << " " << *lda << " " << *ldb<< "\n";
|
---|
| 69 | typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, internal::level3_blocking<Scalar,Scalar>&);
|
---|
| 70 | static functype func[32];
|
---|
| 71 |
|
---|
| 72 | static bool init = false;
|
---|
| 73 | if(!init)
|
---|
| 74 | {
|
---|
| 75 | for(int k=0; k<32; ++k)
|
---|
| 76 | func[k] = 0;
|
---|
| 77 |
|
---|
| 78 | func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,ColMajor,ColMajor>::run);
|
---|
| 79 | func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,RowMajor,ColMajor>::run);
|
---|
| 80 | func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, Conj, RowMajor,ColMajor>::run);
|
---|
| 81 |
|
---|
| 82 | func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,ColMajor,ColMajor>::run);
|
---|
| 83 | func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,RowMajor,ColMajor>::run);
|
---|
| 84 | func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, Conj, RowMajor,ColMajor>::run);
|
---|
| 85 |
|
---|
| 86 | func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|0, false,ColMajor,ColMajor>::run);
|
---|
| 87 | func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, false,RowMajor,ColMajor>::run);
|
---|
| 88 | func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|0, Conj, RowMajor,ColMajor>::run);
|
---|
| 89 |
|
---|
| 90 | func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|0, false,ColMajor,ColMajor>::run);
|
---|
| 91 | func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, false,RowMajor,ColMajor>::run);
|
---|
| 92 | func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|0, Conj, RowMajor,ColMajor>::run);
|
---|
| 93 |
|
---|
| 94 |
|
---|
| 95 | func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,ColMajor,ColMajor>::run);
|
---|
| 96 | func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,RowMajor,ColMajor>::run);
|
---|
| 97 | func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
|
---|
| 98 |
|
---|
| 99 | func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,ColMajor,ColMajor>::run);
|
---|
| 100 | func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,RowMajor,ColMajor>::run);
|
---|
| 101 | func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,Conj, RowMajor,ColMajor>::run);
|
---|
| 102 |
|
---|
| 103 | func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Lower|UnitDiag,false,ColMajor,ColMajor>::run);
|
---|
| 104 | func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,false,RowMajor,ColMajor>::run);
|
---|
| 105 | func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheLeft, Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
|
---|
| 106 |
|
---|
| 107 | func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Lower|UnitDiag,false,ColMajor,ColMajor>::run);
|
---|
| 108 | func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,false,RowMajor,ColMajor>::run);
|
---|
| 109 | func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::triangular_solve_matrix<Scalar,DenseIndex,OnTheRight,Upper|UnitDiag,Conj, RowMajor,ColMajor>::run);
|
---|
| 110 |
|
---|
| 111 | init = true;
|
---|
| 112 | }
|
---|
| 113 |
|
---|
| 114 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 115 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 116 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 117 |
|
---|
| 118 | int info = 0;
|
---|
| 119 | if(SIDE(*side)==INVALID) info = 1;
|
---|
| 120 | else if(UPLO(*uplo)==INVALID) info = 2;
|
---|
| 121 | else if(OP(*opa)==INVALID) info = 3;
|
---|
| 122 | else if(DIAG(*diag)==INVALID) info = 4;
|
---|
| 123 | else if(*m<0) info = 5;
|
---|
| 124 | else if(*n<0) info = 6;
|
---|
| 125 | else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9;
|
---|
| 126 | else if(*ldb<std::max(1,*m)) info = 11;
|
---|
| 127 | if(info)
|
---|
| 128 | return xerbla_(SCALAR_SUFFIX_UP"TRSM ",&info,6);
|
---|
| 129 |
|
---|
| 130 | int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
|
---|
| 131 |
|
---|
| 132 | if(SIDE(*side)==LEFT)
|
---|
| 133 | {
|
---|
| 134 | internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m);
|
---|
| 135 | func[code](*m, *n, a, *lda, b, *ldb, blocking);
|
---|
| 136 | }
|
---|
| 137 | else
|
---|
| 138 | {
|
---|
| 139 | internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n);
|
---|
| 140 | func[code](*n, *m, a, *lda, b, *ldb, blocking);
|
---|
| 141 | }
|
---|
| 142 |
|
---|
| 143 | if(alpha!=Scalar(1))
|
---|
| 144 | matrix(b,*m,*n,*ldb) *= alpha;
|
---|
| 145 |
|
---|
| 146 | return 0;
|
---|
| 147 | }
|
---|
| 148 |
|
---|
| 149 |
|
---|
| 150 | // b = alpha*op(a)*b for side = 'L'or'l'
|
---|
| 151 | // b = alpha*b*op(a) for side = 'R'or'r'
|
---|
| 152 | 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)
|
---|
| 153 | {
|
---|
| 154 | // std::cerr << "in trmm " << *side << " " << *uplo << " " << *opa << " " << *diag << " " << *m << " " << *n << " " << *lda << " " << *ldb << " " << *palpha << "\n";
|
---|
| 155 | typedef void (*functype)(DenseIndex, DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&, internal::level3_blocking<Scalar,Scalar>&);
|
---|
| 156 | static functype func[32];
|
---|
| 157 | static bool init = false;
|
---|
| 158 | if(!init)
|
---|
| 159 | {
|
---|
| 160 | for(int k=0; k<32; ++k)
|
---|
| 161 | func[k] = 0;
|
---|
| 162 |
|
---|
| 163 | func[NOTR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 164 | func[TR | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 165 | func[ADJ | (LEFT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
|
---|
| 166 |
|
---|
| 167 | func[NOTR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 168 | func[TR | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,false,ColMajor>::run);
|
---|
| 169 | func[ADJ | (RIGHT << 2) | (UP << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
|
---|
| 170 |
|
---|
| 171 | func[NOTR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, true, ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 172 | func[TR | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 173 | func[ADJ | (LEFT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
|
---|
| 174 |
|
---|
| 175 | func[NOTR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|0, false,ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 176 | func[TR | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,false,ColMajor>::run);
|
---|
| 177 | func[ADJ | (RIGHT << 2) | (LO << 3) | (NUNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|0, false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
|
---|
| 178 |
|
---|
| 179 | func[NOTR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 180 | func[TR | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 181 | func[ADJ | (LEFT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
|
---|
| 182 |
|
---|
| 183 | func[NOTR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 184 | func[TR | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
|
---|
| 185 | func[ADJ | (RIGHT << 2) | (UP << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
|
---|
| 186 |
|
---|
| 187 | func[NOTR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,true, ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 188 | func[TR | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 189 | func[ADJ | (LEFT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,true, RowMajor,Conj, ColMajor,false,ColMajor>::run);
|
---|
| 190 |
|
---|
| 191 | func[NOTR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Lower|UnitDiag,false,ColMajor,false,ColMajor,false,ColMajor>::run);
|
---|
| 192 | func[TR | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,false,ColMajor>::run);
|
---|
| 193 | func[ADJ | (RIGHT << 2) | (LO << 3) | (UNIT << 4)] = (internal::product_triangular_matrix_matrix<Scalar,DenseIndex,Upper|UnitDiag,false,ColMajor,false,RowMajor,Conj, ColMajor>::run);
|
---|
| 194 |
|
---|
| 195 | init = true;
|
---|
| 196 | }
|
---|
| 197 |
|
---|
| 198 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 199 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 200 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 201 |
|
---|
| 202 | int info = 0;
|
---|
| 203 | if(SIDE(*side)==INVALID) info = 1;
|
---|
| 204 | else if(UPLO(*uplo)==INVALID) info = 2;
|
---|
| 205 | else if(OP(*opa)==INVALID) info = 3;
|
---|
| 206 | else if(DIAG(*diag)==INVALID) info = 4;
|
---|
| 207 | else if(*m<0) info = 5;
|
---|
| 208 | else if(*n<0) info = 6;
|
---|
| 209 | else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 9;
|
---|
| 210 | else if(*ldb<std::max(1,*m)) info = 11;
|
---|
| 211 | if(info)
|
---|
| 212 | return xerbla_(SCALAR_SUFFIX_UP"TRMM ",&info,6);
|
---|
| 213 |
|
---|
| 214 | int code = OP(*opa) | (SIDE(*side) << 2) | (UPLO(*uplo) << 3) | (DIAG(*diag) << 4);
|
---|
| 215 |
|
---|
| 216 | if(*m==0 || *n==0)
|
---|
| 217 | return 1;
|
---|
| 218 |
|
---|
| 219 | // FIXME find a way to avoid this copy
|
---|
| 220 | Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp = matrix(b,*m,*n,*ldb);
|
---|
| 221 | matrix(b,*m,*n,*ldb).setZero();
|
---|
| 222 |
|
---|
| 223 | if(SIDE(*side)==LEFT)
|
---|
| 224 | {
|
---|
| 225 | internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*m);
|
---|
| 226 | func[code](*m, *n, *m, a, *lda, tmp.data(), tmp.outerStride(), b, *ldb, alpha, blocking);
|
---|
| 227 | }
|
---|
| 228 | else
|
---|
| 229 | {
|
---|
| 230 | internal::gemm_blocking_space<ColMajor,Scalar,Scalar,Dynamic,Dynamic,Dynamic,4> blocking(*m,*n,*n);
|
---|
| 231 | func[code](*m, *n, *n, tmp.data(), tmp.outerStride(), a, *lda, b, *ldb, alpha, blocking);
|
---|
| 232 | }
|
---|
| 233 | return 1;
|
---|
| 234 | }
|
---|
| 235 |
|
---|
| 236 | // c = alpha*a*b + beta*c for side = 'L'or'l'
|
---|
| 237 | // c = alpha*b*a + beta*c for side = 'R'or'r
|
---|
| 238 | 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)
|
---|
| 239 | {
|
---|
| 240 | // std::cerr << "in symm " << *side << " " << *uplo << " " << *m << "x" << *n << " lda:" << *lda << " ldb:" << *ldb << " ldc:" << *ldc << " alpha:" << *palpha << " beta:" << *pbeta << "\n";
|
---|
| 241 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 242 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 243 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 244 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 245 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
| 246 |
|
---|
| 247 | int info = 0;
|
---|
| 248 | if(SIDE(*side)==INVALID) info = 1;
|
---|
| 249 | else if(UPLO(*uplo)==INVALID) info = 2;
|
---|
| 250 | else if(*m<0) info = 3;
|
---|
| 251 | else if(*n<0) info = 4;
|
---|
| 252 | else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7;
|
---|
| 253 | else if(*ldb<std::max(1,*m)) info = 9;
|
---|
| 254 | else if(*ldc<std::max(1,*m)) info = 12;
|
---|
| 255 | if(info)
|
---|
| 256 | return xerbla_(SCALAR_SUFFIX_UP"SYMM ",&info,6);
|
---|
| 257 |
|
---|
| 258 | if(beta!=Scalar(1))
|
---|
| 259 | {
|
---|
| 260 | if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
|
---|
| 261 | else matrix(c, *m, *n, *ldc) *= beta;
|
---|
| 262 | }
|
---|
| 263 |
|
---|
| 264 | if(*m==0 || *n==0)
|
---|
| 265 | {
|
---|
| 266 | return 1;
|
---|
| 267 | }
|
---|
| 268 |
|
---|
| 269 | #if ISCOMPLEX
|
---|
| 270 | // FIXME add support for symmetric complex matrix
|
---|
| 271 | int size = (SIDE(*side)==LEFT) ? (*m) : (*n);
|
---|
| 272 | Matrix<Scalar,Dynamic,Dynamic,ColMajor> matA(size,size);
|
---|
| 273 | if(UPLO(*uplo)==UP)
|
---|
| 274 | {
|
---|
| 275 | matA.triangularView<Upper>() = matrix(a,size,size,*lda);
|
---|
| 276 | matA.triangularView<Lower>() = matrix(a,size,size,*lda).transpose();
|
---|
| 277 | }
|
---|
| 278 | else if(UPLO(*uplo)==LO)
|
---|
| 279 | {
|
---|
| 280 | matA.triangularView<Lower>() = matrix(a,size,size,*lda);
|
---|
| 281 | matA.triangularView<Upper>() = matrix(a,size,size,*lda).transpose();
|
---|
| 282 | }
|
---|
| 283 | if(SIDE(*side)==LEFT)
|
---|
| 284 | matrix(c, *m, *n, *ldc) += alpha * matA * matrix(b, *m, *n, *ldb);
|
---|
| 285 | else if(SIDE(*side)==RIGHT)
|
---|
| 286 | matrix(c, *m, *n, *ldc) += alpha * matrix(b, *m, *n, *ldb) * matA;
|
---|
| 287 | #else
|
---|
| 288 | if(SIDE(*side)==LEFT)
|
---|
| 289 | 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);
|
---|
| 290 | 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);
|
---|
| 291 | else return 0;
|
---|
| 292 | else if(SIDE(*side)==RIGHT)
|
---|
| 293 | 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);
|
---|
| 294 | 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);
|
---|
| 295 | else return 0;
|
---|
| 296 | else
|
---|
| 297 | return 0;
|
---|
| 298 | #endif
|
---|
| 299 |
|
---|
| 300 | return 0;
|
---|
| 301 | }
|
---|
| 302 |
|
---|
| 303 | // c = alpha*a*a' + beta*c for op = 'N'or'n'
|
---|
| 304 | // c = alpha*a'*a + beta*c for op = 'T'or't','C'or'c'
|
---|
| 305 | 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)
|
---|
| 306 | {
|
---|
| 307 | // std::cerr << "in syrk " << *uplo << " " << *op << " " << *n << " " << *k << " " << *palpha << " " << *lda << " " << *pbeta << " " << *ldc << "\n";
|
---|
| 308 | #if !ISCOMPLEX
|
---|
| 309 | typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&);
|
---|
| 310 | static functype func[8];
|
---|
| 311 |
|
---|
| 312 | static bool init = false;
|
---|
| 313 | if(!init)
|
---|
| 314 | {
|
---|
| 315 | for(int k=0; k<8; ++k)
|
---|
| 316 | func[k] = 0;
|
---|
| 317 |
|
---|
| 318 | func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Upper>::run);
|
---|
| 319 | func[TR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Upper>::run);
|
---|
| 320 | func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Upper>::run);
|
---|
| 321 |
|
---|
| 322 | func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,ColMajor,Conj, Lower>::run);
|
---|
| 323 | func[TR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,false,Scalar,ColMajor,ColMajor,Conj, Lower>::run);
|
---|
| 324 | func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,ColMajor,false,Lower>::run);
|
---|
| 325 |
|
---|
| 326 | init = true;
|
---|
| 327 | }
|
---|
| 328 | #endif
|
---|
| 329 |
|
---|
| 330 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 331 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 332 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 333 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
| 334 |
|
---|
| 335 | int info = 0;
|
---|
| 336 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
| 337 | else if(OP(*op)==INVALID) info = 2;
|
---|
| 338 | else if(*n<0) info = 3;
|
---|
| 339 | else if(*k<0) info = 4;
|
---|
| 340 | else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
|
---|
| 341 | else if(*ldc<std::max(1,*n)) info = 10;
|
---|
| 342 | if(info)
|
---|
| 343 | return xerbla_(SCALAR_SUFFIX_UP"SYRK ",&info,6);
|
---|
| 344 |
|
---|
| 345 | if(beta!=Scalar(1))
|
---|
| 346 | {
|
---|
| 347 | if(UPLO(*uplo)==UP)
|
---|
| 348 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
|
---|
| 349 | else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
|
---|
| 350 | else
|
---|
| 351 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
|
---|
| 352 | else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
|
---|
| 353 | }
|
---|
| 354 |
|
---|
| 355 | #if ISCOMPLEX
|
---|
| 356 | // FIXME add support for symmetric complex matrix
|
---|
| 357 | if(UPLO(*uplo)==UP)
|
---|
| 358 | {
|
---|
| 359 | if(OP(*op)==NOTR)
|
---|
| 360 | matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
|
---|
| 361 | else
|
---|
| 362 | matrix(c, *n, *n, *ldc).triangularView<Upper>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
|
---|
| 363 | }
|
---|
| 364 | else
|
---|
| 365 | {
|
---|
| 366 | if(OP(*op)==NOTR)
|
---|
| 367 | matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*n,*k,*lda) * matrix(a,*n,*k,*lda).transpose();
|
---|
| 368 | else
|
---|
| 369 | matrix(c, *n, *n, *ldc).triangularView<Lower>() += alpha * matrix(a,*k,*n,*lda).transpose() * matrix(a,*k,*n,*lda);
|
---|
| 370 | }
|
---|
| 371 | #else
|
---|
| 372 | int code = OP(*op) | (UPLO(*uplo) << 2);
|
---|
| 373 | func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
|
---|
| 374 | #endif
|
---|
| 375 |
|
---|
| 376 | return 0;
|
---|
| 377 | }
|
---|
| 378 |
|
---|
| 379 | // c = alpha*a*b' + alpha*b*a' + beta*c for op = 'N'or'n'
|
---|
| 380 | // c = alpha*a'*b + alpha*b'*a + beta*c for op = 'T'or't'
|
---|
| 381 | 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)
|
---|
| 382 | {
|
---|
| 383 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 384 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 385 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 386 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 387 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
| 388 |
|
---|
| 389 | int info = 0;
|
---|
| 390 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
| 391 | else if(OP(*op)==INVALID) info = 2;
|
---|
| 392 | else if(*n<0) info = 3;
|
---|
| 393 | else if(*k<0) info = 4;
|
---|
| 394 | else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
|
---|
| 395 | else if(*ldb<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 9;
|
---|
| 396 | else if(*ldc<std::max(1,*n)) info = 12;
|
---|
| 397 | if(info)
|
---|
| 398 | return xerbla_(SCALAR_SUFFIX_UP"SYR2K",&info,6);
|
---|
| 399 |
|
---|
| 400 | if(beta!=Scalar(1))
|
---|
| 401 | {
|
---|
| 402 | if(UPLO(*uplo)==UP)
|
---|
| 403 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
|
---|
| 404 | else matrix(c, *n, *n, *ldc).triangularView<Upper>() *= beta;
|
---|
| 405 | else
|
---|
| 406 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
|
---|
| 407 | else matrix(c, *n, *n, *ldc).triangularView<Lower>() *= beta;
|
---|
| 408 | }
|
---|
| 409 |
|
---|
| 410 | if(*k==0)
|
---|
| 411 | return 1;
|
---|
| 412 |
|
---|
| 413 | if(OP(*op)==NOTR)
|
---|
| 414 | {
|
---|
| 415 | if(UPLO(*uplo)==UP)
|
---|
| 416 | {
|
---|
| 417 | matrix(c, *n, *n, *ldc).triangularView<Upper>()
|
---|
| 418 | += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
|
---|
| 419 | + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
|
---|
| 420 | }
|
---|
| 421 | else if(UPLO(*uplo)==LO)
|
---|
| 422 | matrix(c, *n, *n, *ldc).triangularView<Lower>()
|
---|
| 423 | += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).transpose()
|
---|
| 424 | + alpha*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).transpose();
|
---|
| 425 | }
|
---|
| 426 | else if(OP(*op)==TR || OP(*op)==ADJ)
|
---|
| 427 | {
|
---|
| 428 | if(UPLO(*uplo)==UP)
|
---|
| 429 | matrix(c, *n, *n, *ldc).triangularView<Upper>()
|
---|
| 430 | += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
|
---|
| 431 | + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
|
---|
| 432 | else if(UPLO(*uplo)==LO)
|
---|
| 433 | matrix(c, *n, *n, *ldc).triangularView<Lower>()
|
---|
| 434 | += alpha*matrix(a, *k, *n, *lda).transpose()*matrix(b, *k, *n, *ldb)
|
---|
| 435 | + alpha*matrix(b, *k, *n, *ldb).transpose()*matrix(a, *k, *n, *lda);
|
---|
| 436 | }
|
---|
| 437 |
|
---|
| 438 | return 0;
|
---|
| 439 | }
|
---|
| 440 |
|
---|
| 441 |
|
---|
| 442 | #if ISCOMPLEX
|
---|
| 443 |
|
---|
| 444 | // c = alpha*a*b + beta*c for side = 'L'or'l'
|
---|
| 445 | // c = alpha*b*a + beta*c for side = 'R'or'r
|
---|
| 446 | 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)
|
---|
| 447 | {
|
---|
| 448 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 449 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 450 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 451 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 452 | Scalar beta = *reinterpret_cast<Scalar*>(pbeta);
|
---|
| 453 |
|
---|
| 454 | // std::cerr << "in hemm " << *side << " " << *uplo << " " << *m << " " << *n << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
|
---|
| 455 |
|
---|
| 456 | int info = 0;
|
---|
| 457 | if(SIDE(*side)==INVALID) info = 1;
|
---|
| 458 | else if(UPLO(*uplo)==INVALID) info = 2;
|
---|
| 459 | else if(*m<0) info = 3;
|
---|
| 460 | else if(*n<0) info = 4;
|
---|
| 461 | else if(*lda<std::max(1,(SIDE(*side)==LEFT)?*m:*n)) info = 7;
|
---|
| 462 | else if(*ldb<std::max(1,*m)) info = 9;
|
---|
| 463 | else if(*ldc<std::max(1,*m)) info = 12;
|
---|
| 464 | if(info)
|
---|
| 465 | return xerbla_(SCALAR_SUFFIX_UP"HEMM ",&info,6);
|
---|
| 466 |
|
---|
| 467 | if(beta==Scalar(0)) matrix(c, *m, *n, *ldc).setZero();
|
---|
| 468 | else if(beta!=Scalar(1)) matrix(c, *m, *n, *ldc) *= beta;
|
---|
| 469 |
|
---|
| 470 | if(*m==0 || *n==0)
|
---|
| 471 | {
|
---|
| 472 | return 1;
|
---|
| 473 | }
|
---|
| 474 |
|
---|
| 475 | if(SIDE(*side)==LEFT)
|
---|
| 476 | {
|
---|
| 477 | if(UPLO(*uplo)==UP) internal::product_selfadjoint_matrix<Scalar,DenseIndex,RowMajor,true,Conj, ColMajor,false,false, ColMajor>
|
---|
| 478 | ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
|
---|
| 479 | else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,true,false, ColMajor,false,false, ColMajor>
|
---|
| 480 | ::run(*m, *n, a, *lda, b, *ldb, c, *ldc, alpha);
|
---|
| 481 | else return 0;
|
---|
| 482 | }
|
---|
| 483 | else if(SIDE(*side)==RIGHT)
|
---|
| 484 | {
|
---|
| 485 | 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>
|
---|
| 486 | ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);*/
|
---|
| 487 | else if(UPLO(*uplo)==LO) internal::product_selfadjoint_matrix<Scalar,DenseIndex,ColMajor,false,false, ColMajor,true,false, ColMajor>
|
---|
| 488 | ::run(*m, *n, b, *ldb, a, *lda, c, *ldc, alpha);
|
---|
| 489 | else return 0;
|
---|
| 490 | }
|
---|
| 491 | else
|
---|
| 492 | {
|
---|
| 493 | return 0;
|
---|
| 494 | }
|
---|
| 495 |
|
---|
| 496 | return 0;
|
---|
| 497 | }
|
---|
| 498 |
|
---|
| 499 | // c = alpha*a*conj(a') + beta*c for op = 'N'or'n'
|
---|
| 500 | // c = alpha*conj(a')*a + beta*c for op = 'C'or'c'
|
---|
| 501 | 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)
|
---|
| 502 | {
|
---|
| 503 | typedef void (*functype)(DenseIndex, DenseIndex, const Scalar *, DenseIndex, const Scalar *, DenseIndex, Scalar *, DenseIndex, const Scalar&);
|
---|
| 504 | static functype func[8];
|
---|
| 505 |
|
---|
| 506 | static bool init = false;
|
---|
| 507 | if(!init)
|
---|
| 508 | {
|
---|
| 509 | for(int k=0; k<8; ++k)
|
---|
| 510 | func[k] = 0;
|
---|
| 511 |
|
---|
| 512 | func[NOTR | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Upper>::run);
|
---|
| 513 | func[ADJ | (UP << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Upper>::run);
|
---|
| 514 |
|
---|
| 515 | func[NOTR | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,ColMajor,false,Scalar,RowMajor,Conj, ColMajor,Lower>::run);
|
---|
| 516 | func[ADJ | (LO << 2)] = (internal::general_matrix_matrix_triangular_product<DenseIndex,Scalar,RowMajor,Conj, Scalar,ColMajor,false,ColMajor,Lower>::run);
|
---|
| 517 |
|
---|
| 518 | init = true;
|
---|
| 519 | }
|
---|
| 520 |
|
---|
| 521 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 522 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 523 | RealScalar alpha = *palpha;
|
---|
| 524 | RealScalar beta = *pbeta;
|
---|
| 525 |
|
---|
| 526 | // std::cerr << "in herk " << *uplo << " " << *op << " " << *n << " " << *k << " " << alpha << " " << *lda << " " << beta << " " << *ldc << "\n";
|
---|
| 527 |
|
---|
| 528 | int info = 0;
|
---|
| 529 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
| 530 | else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2;
|
---|
| 531 | else if(*n<0) info = 3;
|
---|
| 532 | else if(*k<0) info = 4;
|
---|
| 533 | else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
|
---|
| 534 | else if(*ldc<std::max(1,*n)) info = 10;
|
---|
| 535 | if(info)
|
---|
| 536 | return xerbla_(SCALAR_SUFFIX_UP"HERK ",&info,6);
|
---|
| 537 |
|
---|
| 538 | int code = OP(*op) | (UPLO(*uplo) << 2);
|
---|
| 539 |
|
---|
| 540 | if(beta!=RealScalar(1))
|
---|
| 541 | {
|
---|
| 542 | if(UPLO(*uplo)==UP)
|
---|
| 543 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
|
---|
| 544 | else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
|
---|
| 545 | else
|
---|
| 546 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
|
---|
| 547 | else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
|
---|
| 548 |
|
---|
| 549 | if(beta!=Scalar(0))
|
---|
| 550 | {
|
---|
| 551 | matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
|
---|
| 552 | matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
|
---|
| 553 | }
|
---|
| 554 | }
|
---|
| 555 |
|
---|
| 556 | if(*k>0 && alpha!=RealScalar(0))
|
---|
| 557 | {
|
---|
| 558 | func[code](*n, *k, a, *lda, a, *lda, c, *ldc, alpha);
|
---|
| 559 | matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
|
---|
| 560 | }
|
---|
| 561 | return 0;
|
---|
| 562 | }
|
---|
| 563 |
|
---|
| 564 | // c = alpha*a*conj(b') + conj(alpha)*b*conj(a') + beta*c, for op = 'N'or'n'
|
---|
| 565 | // c = alpha*conj(a')*b + conj(alpha)*conj(b')*a + beta*c, for op = 'C'or'c'
|
---|
| 566 | 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)
|
---|
| 567 | {
|
---|
| 568 | Scalar* a = reinterpret_cast<Scalar*>(pa);
|
---|
| 569 | Scalar* b = reinterpret_cast<Scalar*>(pb);
|
---|
| 570 | Scalar* c = reinterpret_cast<Scalar*>(pc);
|
---|
| 571 | Scalar alpha = *reinterpret_cast<Scalar*>(palpha);
|
---|
| 572 | RealScalar beta = *pbeta;
|
---|
| 573 |
|
---|
| 574 | int info = 0;
|
---|
| 575 | if(UPLO(*uplo)==INVALID) info = 1;
|
---|
| 576 | else if((OP(*op)==INVALID) || (OP(*op)==TR)) info = 2;
|
---|
| 577 | else if(*n<0) info = 3;
|
---|
| 578 | else if(*k<0) info = 4;
|
---|
| 579 | else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 7;
|
---|
| 580 | else if(*lda<std::max(1,(OP(*op)==NOTR)?*n:*k)) info = 9;
|
---|
| 581 | else if(*ldc<std::max(1,*n)) info = 12;
|
---|
| 582 | if(info)
|
---|
| 583 | return xerbla_(SCALAR_SUFFIX_UP"HER2K",&info,6);
|
---|
| 584 |
|
---|
| 585 | if(beta!=RealScalar(1))
|
---|
| 586 | {
|
---|
| 587 | if(UPLO(*uplo)==UP)
|
---|
| 588 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Upper>().setZero();
|
---|
| 589 | else matrix(c, *n, *n, *ldc).triangularView<StrictlyUpper>() *= beta;
|
---|
| 590 | else
|
---|
| 591 | if(beta==Scalar(0)) matrix(c, *n, *n, *ldc).triangularView<Lower>().setZero();
|
---|
| 592 | else matrix(c, *n, *n, *ldc).triangularView<StrictlyLower>() *= beta;
|
---|
| 593 |
|
---|
| 594 | if(beta!=Scalar(0))
|
---|
| 595 | {
|
---|
| 596 | matrix(c, *n, *n, *ldc).diagonal().real() *= beta;
|
---|
| 597 | matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
|
---|
| 598 | }
|
---|
| 599 | }
|
---|
| 600 | else if(*k>0 && alpha!=Scalar(0))
|
---|
| 601 | matrix(c, *n, *n, *ldc).diagonal().imag().setZero();
|
---|
| 602 |
|
---|
| 603 | if(*k==0)
|
---|
| 604 | return 1;
|
---|
| 605 |
|
---|
| 606 | if(OP(*op)==NOTR)
|
---|
| 607 | {
|
---|
| 608 | if(UPLO(*uplo)==UP)
|
---|
| 609 | {
|
---|
| 610 | matrix(c, *n, *n, *ldc).triangularView<Upper>()
|
---|
| 611 | += alpha *matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
|
---|
| 612 | + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
|
---|
| 613 | }
|
---|
| 614 | else if(UPLO(*uplo)==LO)
|
---|
| 615 | matrix(c, *n, *n, *ldc).triangularView<Lower>()
|
---|
| 616 | += alpha*matrix(a, *n, *k, *lda)*matrix(b, *n, *k, *ldb).adjoint()
|
---|
| 617 | + numext::conj(alpha)*matrix(b, *n, *k, *ldb)*matrix(a, *n, *k, *lda).adjoint();
|
---|
| 618 | }
|
---|
| 619 | else if(OP(*op)==ADJ)
|
---|
| 620 | {
|
---|
| 621 | if(UPLO(*uplo)==UP)
|
---|
| 622 | matrix(c, *n, *n, *ldc).triangularView<Upper>()
|
---|
| 623 | += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
|
---|
| 624 | + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
|
---|
| 625 | else if(UPLO(*uplo)==LO)
|
---|
| 626 | matrix(c, *n, *n, *ldc).triangularView<Lower>()
|
---|
| 627 | += alpha*matrix(a, *k, *n, *lda).adjoint()*matrix(b, *k, *n, *ldb)
|
---|
| 628 | + numext::conj(alpha)*matrix(b, *k, *n, *ldb).adjoint()*matrix(a, *k, *n, *lda);
|
---|
| 629 | }
|
---|
| 630 |
|
---|
| 631 | return 1;
|
---|
| 632 | }
|
---|
| 633 |
|
---|
| 634 | #endif // ISCOMPLEX
|
---|