From 338381e1606472f6127eb472ab6ed2d55393360e Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Fri, 23 Oct 2026 16:34:08 +0200 Subject: [PATCH] fixed some conjugation issues in cblas replacement wrappers --- noncblas.cc | 120 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 111 insertions(+), 9 deletions(-) diff --git a/noncblas.cc b/noncblas.cc index 316594a..c85fbba 100644 --- a/noncblas.cc +++ b/noncblas.cc @@ -244,29 +244,52 @@ float cblas_scnrm2(const int N, const void *X, const int incX){ -//the following ones are f2c-compatible, but is it truly portable??? +//Fortran functions returning complex: gfortran (and most current compilers) return the value as C99 complex, +//the old f2c/g77 convention (and gfortran -ff2c) passes a hidden pointer to the result as the first argument - define F2C_COMPLEX_RETURN for that case +#ifdef F2C_COMPLEX_RETURN extern "C" void FORNAME(zdotu) (void *retval, const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +extern "C" void FORNAME(zdotc) (void *retval, const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +#else +extern "C" std::complex FORNAME(zdotu) (const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +extern "C" std::complex FORNAME(zdotc) (const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +#endif + void cblas_zdotu_sub(const int N, const void *X, const int incX, const void *Y, const int incY, void *dotu){ #ifdef FORINT const FINT ntmp=N; const FINT incxtmp=incX; const FINT incytmp=incY; +#ifdef F2C_COMPLEX_RETURN FORNAME(zdotu) (dotu,&ntmp,X,&incxtmp,Y,&incytmp); #else + *(std::complex *)dotu = FORNAME(zdotu) (&ntmp,X,&incxtmp,Y,&incytmp); +#endif +#else +#ifdef F2C_COMPLEX_RETURN FORNAME(zdotu) (dotu,&N,X,&incX,Y,&incY); +#else + *(std::complex *)dotu = FORNAME(zdotu) (&N,X,&incX,Y,&incY); +#endif #endif } -extern "C" void FORNAME(zdotc) (void *retval, const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); void cblas_zdotc_sub(const int N, const void *X, const int incX, const void *Y, const int incY, void *dotc){ #ifdef FORINT const FINT ntmp=N; const FINT incxtmp=incX; const FINT incytmp=incY; +#ifdef F2C_COMPLEX_RETURN FORNAME(zdotc) (dotc,&ntmp,X,&incxtmp,Y,&incytmp); #else + *(std::complex *)dotc = FORNAME(zdotc) (&ntmp,X,&incxtmp,Y,&incytmp); +#endif +#else +#ifdef F2C_COMPLEX_RETURN FORNAME(zdotc) (dotc,&N,X,&incX,Y,&incY); +#else + *(std::complex *)dotc = FORNAME(zdotc) (&N,X,&incX,Y,&incY); +#endif #endif } @@ -327,14 +350,72 @@ void cblas_zhpmv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); if(Uplo!=CblasLower) laerror("CblasLower uplo asserted"); char U = BLAS_FORTRANCASE('u'); +//row-major lower packed hermitian H is column-major upper packed conj(H), therefore +//conj(y) = conj(beta) conj(y) + conj(alpha) conj(H) conj(x) is computed by zhpmv on conjugated vectors and scalars +std::complex alphac = std::conj(*(const std::complex *)alpha); +std::complex betac = std::conj(*(const std::complex *)beta); +std::complex *p; +p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} +p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} #ifdef FORINT const FINT ntmp=N; const FINT incxtmp=incX; const FINT incytmp=incY; - FORNAME(zhpmv) (&U,&ntmp, alpha, Ap, X, &incxtmp, beta, Y, &incytmp); + FORNAME(zhpmv) (&U,&ntmp, &alphac, Ap, X, &incxtmp, &betac, Y, &incytmp); #else - FORNAME(zhpmv) (&U,&N, alpha, Ap, X, &incX, beta, Y, &incY); + FORNAME(zhpmv) (&U,&N, &alphac, Ap, X, &incX, &betac, Y, &incY); #endif +//conjugate back +p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} +p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} +} + + +extern "C" void FORNAME(dspr2) (const char *uplo, const FINT *n, const double *alpha, const double *x, const FINT *incx, const double *y, const FINT *incy, double *ap); +void cblas_dspr2(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, + const int N, const double alpha, const double *X, + const int incX, const double *Y, const int incY, double *Ap) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); + if(Uplo!=CblasLower) laerror("CblasLower uplo asserted"); + char U = BLAS_FORTRANCASE('u'); +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + FORNAME(dspr2) (&U,&ntmp, &alpha, X, &incxtmp, Y, &incytmp, Ap); +#else + FORNAME(dspr2) (&U,&N, &alpha, X, &incX, Y, &incY, Ap); +#endif +} + + +//row-major lower packed hermitian H is column-major upper packed conj(H), therefore +//conj(H) += conj(a) conj(x) conj(y)^H + a conj(y) conj(x)^H is computed by zhpr2 with alpha=a and conjugated swapped vectors +extern "C" void FORNAME(zhpr2) (const char *uplo, const FINT *n, const void *alpha, const void *x, const FINT *incx, const void *y, const FINT *incy, void *ap); +void cblas_zhpr2(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, const int N, + const void *alpha, const void *X, const int incX, + const void *Y, const int incY, void *Ap) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); + if(Uplo!=CblasLower) laerror("CblasLower uplo asserted"); + char U = BLAS_FORTRANCASE('u'); +bool same = (X==Y && incX==incY); +//conjugate x and y +std::complex *p; +p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} +if(!same) {p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;}} +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + FORNAME(zhpr2) (&U,&ntmp, alpha, Y, &incytmp, X, &incxtmp, Ap); +#else + FORNAME(zhpr2) (&U,&N, alpha, Y, &incY, X, &incX, Ap); +#endif +//conjugate back +p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} +if(!same) {p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;}} } @@ -381,7 +462,8 @@ void cblas_sger(const CBLAS_ORDER Order, const int M, const int N, -extern "C" void FORNAME(zgerc) (const FINT *m, const FINT *n, const void *alpha, const void *x, const FINT *incx, const void *y, const FINT *incy, void *a, const FINT *lda); +extern "C" void FORNAME(zgeru) (const FINT *m, const FINT *n, const void *alpha, const void *x, const FINT *incx, const void *y, const FINT *incy, void *a, const FINT *lda); +//row-major A += alpha x y^H is column-major A^T += alpha conj(y) x^T, i.e. zgeru with conjugated y void cblas_zgerc(const CBLAS_ORDER Order, const int M, const int N, const void *alpha, const void *X, const int incX, const void *Y, const int incY, void *A, const int lda) @@ -396,16 +478,15 @@ p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= const FINT incxtmp=incX; const FINT incytmp=incY; const FINT ldatmp=lda; - FORNAME(zgerc) (&ntmp, &mtmp, alpha, Y, &incytmp, X, &incxtmp, A, &ldatmp); + FORNAME(zgeru) (&ntmp, &mtmp, alpha, Y, &incytmp, X, &incxtmp, A, &ldatmp); #else - FORNAME(zgerc)(&N, &M, alpha, Y, &incY, X, &incX, A, &lda); + FORNAME(zgeru)(&N, &M, alpha, Y, &incY, X, &incX, A, &lda); #endif //conjugate y back p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} } -extern "C" void FORNAME(zgeru) (const FINT *m, const FINT *n, const void *alpha, const void *x, const FINT *incx, const void *y, const FINT *incy, void *a, const FINT *lda); void cblas_zgeru(const CBLAS_ORDER Order, const int M, const int N, const void *alpha, const void *X, const int incX, const void *Y, const int incY, void *A, const int lda) @@ -562,7 +643,6 @@ void cblas_zgemv(const CBLAS_ORDER Order, void *Y, const int incY) { if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); - if(TransA == CblasConjTrans) laerror("zgemv with CblasConjTrans not supportted"); //swap n-m and toggle transposition char transa = BLAS_FORTRANCASE(TransA==CblasNoTrans?'T':'N'); #ifdef FORINT @@ -571,6 +651,28 @@ void cblas_zgemv(const CBLAS_ORDER Order, const FINT ldatmp=lda; const FINT incxtmp=incX; const FINT incytmp=incY; +#endif +if(TransA == CblasConjTrans) + { + //row-major A is column-major A^T, so A^H = conj(A^T) and + //conj(y) = conj(beta) conj(y) + conj(alpha) A^T conj(x) is computed by non-transposed zgemv on conjugated vectors and scalars + const int lx = M, ly = N; //lengths of x and y for op(A)=A^H + std::complex alphac = std::conj(*(const std::complex *)alpha); + std::complex betac = std::conj(*(const std::complex *)beta); + std::complex *p; + p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} + p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} +#ifdef FORINT + FORNAME(zgemv) (&transa, &ntmp, &mtmp, &alphac, A, &ldatmp, X, &incxtmp, &betac, Y, &incytmp ); +#else + FORNAME(zgemv) (&transa, &N, &M, &alphac, A, &lda, X, &incX, &betac, Y, &incY ); +#endif + //conjugate back + p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= incX;} + p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} + return; + } +#ifdef FORINT FORNAME(zgemv) (&transa, &ntmp, &mtmp, alpha, A, &ldatmp, X, &incxtmp, beta, Y, &incytmp ); #else FORNAME(zgemv) (&transa, &N, &M, alpha, A, &lda, X, &incX, beta, Y, &incY );