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 (at) 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 /** ZHEMV performs the matrix-vector operation 13 * 14 * y := alpha*A*x + beta*y, 15 * 16 * where alpha and beta are scalars, x and y are n element vectors and 17 * A is an n by n hermitian matrix. 18 */ 19 int EIGEN_BLAS_FUNC(hemv)(char *uplo, int *n, RealScalar *palpha, RealScalar *pa, int *lda, RealScalar *px, int *incx, RealScalar *pbeta, RealScalar *py, int *incy) 20 { 21 Scalar* a = reinterpret_cast<Scalar*>(pa); 22 Scalar* x = reinterpret_cast<Scalar*>(px); 23 Scalar* y = reinterpret_cast<Scalar*>(py); 24 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 25 Scalar beta = *reinterpret_cast<Scalar*>(pbeta); 26 27 // check arguments 28 int info = 0; 29 if(UPLO(*uplo)==INVALID) info = 1; 30 else if(*n<0) info = 2; 31 else if(*lda<std::max(1,*n)) info = 5; 32 else if(*incx==0) info = 7; 33 else if(*incy==0) info = 10; 34 if(info) 35 return xerbla_(SCALAR_SUFFIX_UP"HEMV ",&info,6); 36 37 if(*n==0) 38 return 1; 39 40 Scalar* actual_x = get_compact_vector(x,*n,*incx); 41 Scalar* actual_y = get_compact_vector(y,*n,*incy); 42 43 if(beta!=Scalar(1)) 44 { 45 if(beta==Scalar(0)) vector(actual_y, *n).setZero(); 46 else vector(actual_y, *n) *= beta; 47 } 48 49 if(alpha!=Scalar(0)) 50 { 51 // TODO performs a direct call to the underlying implementation function 52 if(UPLO(*uplo)==UP) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Upper>() * (alpha * vector(actual_x,*n)); 53 else if(UPLO(*uplo)==LO) vector(actual_y,*n).noalias() += matrix(a,*n,*n,*lda).selfadjointView<Lower>() * (alpha * vector(actual_x,*n)); 54 } 55 56 if(actual_x!=x) delete[] actual_x; 57 if(actual_y!=y) delete[] copy_back(actual_y,y,*n,*incy); 58 59 return 1; 60 } 61 62 /** ZHBMV performs the matrix-vector operation 63 * 64 * y := alpha*A*x + beta*y, 65 * 66 * where alpha and beta are scalars, x and y are n element vectors and 67 * A is an n by n hermitian band matrix, with k super-diagonals. 68 */ 69 // int EIGEN_BLAS_FUNC(hbmv)(char *uplo, int *n, int *k, RealScalar *alpha, RealScalar *a, int *lda, 70 // RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 71 // { 72 // return 1; 73 // } 74 75 /** ZHPMV performs the matrix-vector operation 76 * 77 * y := alpha*A*x + beta*y, 78 * 79 * where alpha and beta are scalars, x and y are n element vectors and 80 * A is an n by n hermitian matrix, supplied in packed form. 81 */ 82 // int EIGEN_BLAS_FUNC(hpmv)(char *uplo, int *n, RealScalar *alpha, RealScalar *ap, RealScalar *x, int *incx, RealScalar *beta, RealScalar *y, int *incy) 83 // { 84 // return 1; 85 // } 86 87 /** ZHPR performs the hermitian rank 1 operation 88 * 89 * A := alpha*x*conjg( x' ) + A, 90 * 91 * where alpha is a real scalar, x is an n element vector and A is an 92 * n by n hermitian matrix, supplied in packed form. 93 */ 94 // int EIGEN_BLAS_FUNC(hpr)(char *uplo, int *n, RealScalar *alpha, RealScalar *x, int *incx, RealScalar *ap) 95 // { 96 // return 1; 97 // } 98 99 /** ZHPR2 performs the hermitian rank 2 operation 100 * 101 * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, 102 * 103 * where alpha is a scalar, x and y are n element vectors and A is an 104 * n by n hermitian matrix, supplied in packed form. 105 */ 106 // int EIGEN_BLAS_FUNC(hpr2)(char *uplo, int *n, RealScalar *palpha, RealScalar *x, int *incx, RealScalar *y, int *incy, RealScalar *ap) 107 // { 108 // return 1; 109 // } 110 111 /** ZHER performs the hermitian rank 1 operation 112 * 113 * A := alpha*x*conjg( x' ) + A, 114 * 115 * where alpha is a real scalar, x is an n element vector and A is an 116 * n by n hermitian matrix. 117 */ 118 int EIGEN_BLAS_FUNC(her)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *pa, int *lda) 119 { 120 Scalar* x = reinterpret_cast<Scalar*>(px); 121 Scalar* a = reinterpret_cast<Scalar*>(pa); 122 RealScalar alpha = *reinterpret_cast<RealScalar*>(palpha); 123 124 int info = 0; 125 if(UPLO(*uplo)==INVALID) info = 1; 126 else if(*n<0) info = 2; 127 else if(*incx==0) info = 5; 128 else if(*lda<std::max(1,*n)) info = 7; 129 if(info) 130 return xerbla_(SCALAR_SUFFIX_UP"HER ",&info,6); 131 132 if(alpha==RealScalar(0)) 133 return 1; 134 135 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 136 137 // TODO perform direct calls to underlying implementation 138 // if(UPLO(*uplo)==LO) matrix(a,*n,*n,*lda).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n), alpha); 139 // else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n), alpha); 140 141 if(UPLO(*uplo)==LO) 142 for(int j=0;j<*n;++j) 143 matrix(a,*n,*n,*lda).col(j).tail(*n-j) += alpha * internal::conj(x_cpy[j]) * vector(x_cpy+j,*n-j); 144 else 145 for(int j=0;j<*n;++j) 146 matrix(a,*n,*n,*lda).col(j).head(j+1) += alpha * internal::conj(x_cpy[j]) * vector(x_cpy,j+1); 147 148 matrix(a,*n,*n,*lda).diagonal().imag().setZero(); 149 150 if(x_cpy!=x) delete[] x_cpy; 151 152 return 1; 153 } 154 155 /** ZHER2 performs the hermitian rank 2 operation 156 * 157 * A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, 158 * 159 * where alpha is a scalar, x and y are n element vectors and A is an n 160 * by n hermitian matrix. 161 */ 162 int EIGEN_BLAS_FUNC(her2)(char *uplo, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 163 { 164 Scalar* x = reinterpret_cast<Scalar*>(px); 165 Scalar* y = reinterpret_cast<Scalar*>(py); 166 Scalar* a = reinterpret_cast<Scalar*>(pa); 167 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 168 169 int info = 0; 170 if(UPLO(*uplo)==INVALID) info = 1; 171 else if(*n<0) info = 2; 172 else if(*incx==0) info = 5; 173 else if(*incy==0) info = 7; 174 else if(*lda<std::max(1,*n)) info = 9; 175 if(info) 176 return xerbla_(SCALAR_SUFFIX_UP"HER2 ",&info,6); 177 178 if(alpha==Scalar(0)) 179 return 1; 180 181 Scalar* x_cpy = get_compact_vector(x, *n, *incx); 182 Scalar* y_cpy = get_compact_vector(y, *n, *incy); 183 184 // TODO perform direct calls to underlying implementation 185 if(UPLO(*uplo)==LO) matrix(a,*n,*n,*lda).selfadjointView<Lower>().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); 186 else if(UPLO(*uplo)==UP) matrix(a,*n,*n,*lda).selfadjointView<Upper>().rankUpdate(vector(x_cpy,*n),vector(y_cpy,*n),alpha); 187 188 matrix(a,*n,*n,*lda).diagonal().imag().setZero(); 189 190 if(x_cpy!=x) delete[] x_cpy; 191 if(y_cpy!=y) delete[] y_cpy; 192 193 return 1; 194 } 195 196 /** ZGERU performs the rank 1 operation 197 * 198 * A := alpha*x*y' + A, 199 * 200 * where alpha is a scalar, x is an m element vector, y is an n element 201 * vector and A is an m by n matrix. 202 */ 203 int EIGEN_BLAS_FUNC(geru)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 204 { 205 Scalar* x = reinterpret_cast<Scalar*>(px); 206 Scalar* y = reinterpret_cast<Scalar*>(py); 207 Scalar* a = reinterpret_cast<Scalar*>(pa); 208 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 209 210 int info = 0; 211 if(*m<0) info = 1; 212 else if(*n<0) info = 2; 213 else if(*incx==0) info = 5; 214 else if(*incy==0) info = 7; 215 else if(*lda<std::max(1,*m)) info = 9; 216 if(info) 217 return xerbla_(SCALAR_SUFFIX_UP"GERU ",&info,6); 218 219 if(alpha==Scalar(0)) 220 return 1; 221 222 Scalar* x_cpy = get_compact_vector(x,*m,*incx); 223 Scalar* y_cpy = get_compact_vector(y,*n,*incy); 224 225 // TODO perform direct calls to underlying implementation 226 matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).transpose(); 227 228 if(x_cpy!=x) delete[] x_cpy; 229 if(y_cpy!=y) delete[] y_cpy; 230 231 return 1; 232 } 233 234 /** ZGERC performs the rank 1 operation 235 * 236 * A := alpha*x*conjg( y' ) + A, 237 * 238 * where alpha is a scalar, x is an m element vector, y is an n element 239 * vector and A is an m by n matrix. 240 */ 241 int EIGEN_BLAS_FUNC(gerc)(int *m, int *n, RealScalar *palpha, RealScalar *px, int *incx, RealScalar *py, int *incy, RealScalar *pa, int *lda) 242 { 243 Scalar* x = reinterpret_cast<Scalar*>(px); 244 Scalar* y = reinterpret_cast<Scalar*>(py); 245 Scalar* a = reinterpret_cast<Scalar*>(pa); 246 Scalar alpha = *reinterpret_cast<Scalar*>(palpha); 247 248 int info = 0; 249 if(*m<0) info = 1; 250 else if(*n<0) info = 2; 251 else if(*incx==0) info = 5; 252 else if(*incy==0) info = 7; 253 else if(*lda<std::max(1,*m)) info = 9; 254 if(info) 255 return xerbla_(SCALAR_SUFFIX_UP"GERC ",&info,6); 256 257 if(alpha==Scalar(0)) 258 return 1; 259 260 Scalar* x_cpy = get_compact_vector(x,*m,*incx); 261 Scalar* y_cpy = get_compact_vector(y,*n,*incy); 262 263 // TODO perform direct calls to underlying implementation 264 matrix(a,*m,*n,*lda) += alpha * vector(x_cpy,*m) * vector(y_cpy,*n).adjoint(); 265 266 if(x_cpy!=x) delete[] x_cpy; 267 if(y_cpy!=y) delete[] y_cpy; 268 269 return 1; 270 } 271
