diff --git a/noncblas.cc b/noncblas.cc index e2471f7..a44bfe9 100644 --- a/noncblas.cc +++ b/noncblas.cc @@ -185,6 +185,18 @@ void cblas_zcopy(const int N, const void *X, const int incX, void *Y, const int #endif } +extern "C" void FORNAME(ccopy) (const FINT *n, const void *x, const FINT *incx, void *y, const FINT *incy); +void cblas_ccopy(const int N, const void *X, const int incX, void *Y, const int incY){ +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + FORNAME(ccopy) (&ntmp,X,&incxtmp,Y,&incytmp); +#else + FORNAME(ccopy) (&N,X,&incX,Y,&incY); +#endif +} + extern "C" void FORNAME(zaxpy) (const FINT *n, const void *a, const void *x, const FINT *incx, void *y, const FINT *incy); void cblas_zaxpy(const int N, const void *alpha, const void *X, const int incX, void *Y, const int incY){ #ifdef FORINT @@ -197,6 +209,18 @@ void cblas_zaxpy(const int N, const void *alpha, const void *X, const int incX, #endif } +extern "C" void FORNAME(caxpy) (const FINT *n, const void *a, const void *x, const FINT *incx, void *y, const FINT *incy); +void cblas_caxpy(const int N, const void *alpha, const void *X, const int incX, void *Y, const int incY){ +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + FORNAME(caxpy) (&ntmp,alpha,X,&incxtmp,Y,&incytmp); +#else + FORNAME(caxpy) (&N,alpha,X,&incX,Y,&incY); +#endif +} + extern "C" void FORNAME(zscal) (const FINT *n, const void *a, void *x, const FINT *incx); void cblas_zscal(const int N, const void *alpha, void *X, const int incX){ #ifdef FORINT @@ -208,6 +232,17 @@ void cblas_zscal(const int N, const void *alpha, void *X, const int incX){ #endif } +extern "C" void FORNAME(cscal) (const FINT *n, const void *a, void *x, const FINT *incx); +void cblas_cscal(const int N, const void *alpha, void *X, const int incX){ +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + FORNAME(cscal)(&ntmp,alpha,X,&incxtmp); +#else + FORNAME(cscal)(&N,alpha,X,&incX); +#endif +} + extern "C" void FORNAME(zdscal) (const FINT *n, const double *a, void *x, const FINT *incx); void cblas_zdscal(const int N, const double alpha, void *X, const int incX){ #ifdef FORINT @@ -219,6 +254,17 @@ void cblas_zdscal(const int N, const double alpha, void *X, const int incX){ #endif } +extern "C" void FORNAME(csscal) (const FINT *n, const float *a, void *x, const FINT *incx); +void cblas_csscal(const int N, const float alpha, void *X, const int incX){ +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + FORNAME(csscal)(&ntmp,&alpha,X,&incxtmp); +#else + FORNAME(csscal)(&N,&alpha,X,&incX); +#endif +} + extern "C" double FORNAME(dznrm2) (const FINT *n, const void *x, const FINT *incx); double cblas_dznrm2(const int N, const void *X, const int incX){ @@ -294,6 +340,54 @@ void cblas_zdotc_sub(const int N, const void *X, const int incX, const void *Y, } +#ifdef F2C_COMPLEX_RETURN +extern "C" void FORNAME(cdotu) (void *retval, const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +extern "C" void FORNAME(cdotc) (void *retval, const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +#else +extern "C" std::complex FORNAME(cdotu) (const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +extern "C" std::complex FORNAME(cdotc) (const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy); +#endif + +void cblas_cdotu_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(cdotu) (dotu,&ntmp,X,&incxtmp,Y,&incytmp); +#else + *(std::complex *)dotu = FORNAME(cdotu) (&ntmp,X,&incxtmp,Y,&incytmp); +#endif +#else +#ifdef F2C_COMPLEX_RETURN + FORNAME(cdotu) (dotu,&N,X,&incX,Y,&incY); +#else + *(std::complex *)dotu = FORNAME(cdotu) (&N,X,&incX,Y,&incY); +#endif +#endif +} + + +void cblas_cdotc_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(cdotc) (dotc,&ntmp,X,&incxtmp,Y,&incytmp); +#else + *(std::complex *)dotc = FORNAME(cdotc) (&ntmp,X,&incxtmp,Y,&incytmp); +#endif +#else +#ifdef F2C_COMPLEX_RETURN + FORNAME(cdotc) (dotc,&N,X,&incX,Y,&incY); +#else + *(std::complex *)dotc = FORNAME(cdotc) (&N,X,&incX,Y,&incY); +#endif +#endif +} + + //Level 2 and Level 3 on symmetric/hermitian packed matrices - straightforward //enum CBLAS_UPLO {CblasUpper=121, CblasLower=122}; //enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102 }; @@ -382,6 +476,47 @@ p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} } +extern "C" void FORNAME(chpmv) (const char *uplo, const FINT *n, const void *alpha, const void *ap, const void *x, const FINT *incx, const void *beta, void *y, const FINT *incy); +void cblas_chpmv(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, + const int N, const void *alpha, const void *Ap, + const void *X, const int incX, + const void *beta, void *Y, const int incY) +{ + char U = BLAS_FORTRANCASE('u'); + if(Order==CblasColMajor && Uplo==CblasUpper) //native fortran convention, no conjugation needed + { +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + FORNAME(chpmv) (&U,&ntmp, alpha, Ap, X, &incxtmp, beta, Y, &incytmp); +#else + FORNAME(chpmv) (&U,&N, alpha, Ap, X, &incX, beta, Y, &incY); +#endif + return; + } + if(Order!=CblasRowMajor) laerror("CblasRowMajor or CblasColMajor+CblasUpper order asserted"); + if(Uplo!=CblasLower) laerror("CblasLower uplo asserted"); +//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 chpmv 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(chpmv) (&U,&ntmp, &alphac, Ap, X, &incxtmp, &betac, Y, &incytmp); +#else + FORNAME(chpmv) (&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, @@ -401,6 +536,24 @@ void cblas_dspr2(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, #endif } +extern "C" void FORNAME(sspr2) (const char *uplo, const FINT *n, const float *alpha, const float *x, const FINT *incx, const float *y, const FINT *incy, float *ap); +void cblas_sspr2(const CBLAS_ORDER Order, const CBLAS_UPLO Uplo, + const int N, const float alpha, const float *X, + const int incX, const float *Y, const int incY, float *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(sspr2) (&U,&ntmp, &alpha, X, &incxtmp, Y, &incytmp, Ap); +#else + FORNAME(sspr2) (&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 @@ -430,6 +583,32 @@ p= (std::complex *) X; for(int i=0; iimag(-p->imag()); p+= if(!same) {p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;}} } +extern "C" void FORNAME(chpr2) (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_chpr2(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(chpr2) (&U,&ntmp, alpha, Y, &incytmp, X, &incxtmp, Ap); +#else + FORNAME(chpr2) (&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;}} +} + //Level 2 and Level 3 on general matrices - take into account the transposed storage of matrices in Fortran and C @@ -516,6 +695,47 @@ void cblas_zgeru(const CBLAS_ORDER Order, const int M, const int N, #endif } +extern "C" void FORNAME(cgeru) (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_cgerc(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) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); +//conjugate y +std::complex *p; +p= (std::complex *) Y; for(int i=0; iimag(-p->imag()); p+= incY;} +#ifdef FORINT + const FINT mtmp=M; + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + const FINT ldatmp=lda; + FORNAME(cgeru) (&ntmp, &mtmp, alpha, Y, &incytmp, X, &incxtmp, A, &ldatmp); +#else + FORNAME(cgeru)(&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;} +} + + +void cblas_cgeru(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) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); +#ifdef FORINT + const FINT mtmp=M; + const FINT ntmp=N; + const FINT incxtmp=incX; + const FINT incytmp=incY; + const FINT ldatmp=lda; + FORNAME(cgeru) (&ntmp, &mtmp, alpha, Y, &incytmp, X, &incxtmp, A, &ldatmp); +#else + FORNAME(cgeru)(&N, &M, alpha, Y, &incY, X, &incX, A, &lda); +#endif +} + extern "C" void FORNAME(dgemm) (const char *transa, const char *transb, const FINT *m, const FINT *n, const FINT *k, const double *alpha, const double *a, const FINT *lda, const double *b, const FINT *ldb, const double *beta, double *c, const FINT *ldc); void cblas_dgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA, @@ -595,6 +815,32 @@ void cblas_zgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA, #endif } +extern "C" void FORNAME(cgemm) (const char *transa, const char *transb, const FINT *m, const FINT *n, const FINT *k, const void *alpha, const void *a, const FINT *lda, const void *b, const FINT *ldb, const void *beta, void *c, const FINT *ldc); +void cblas_cgemm(const CBLAS_ORDER Order, const CBLAS_TRANSPOSE TransA, + const CBLAS_TRANSPOSE TransB, const int M, const int N, + const int K, const void *alpha, const void *A, + const int lda, const void *B, const int ldb, + const void *beta, void *C, const int ldc) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); + //swap a-b, m-n + char transb = BLAS_FORTRANCASE(TransB==CblasConjTrans?'C':(TransB==CblasNoTrans?'N':'T')); + char transa = BLAS_FORTRANCASE(TransA==CblasConjTrans?'C':(TransA==CblasNoTrans?'N':'T')); +#ifdef FORINT + const FINT mtmp=M; + const FINT ntmp=N; + const FINT ktmp=K; + const FINT ldatmp=lda; + const FINT ldbtmp=ldb; + const FINT ldctmp=ldc; + FORNAME(cgemm) (&transb,&transa, + &ntmp, &mtmp, &ktmp, alpha, B, &ldbtmp, A, &ldatmp, beta, C, &ldctmp); +#else + FORNAME(cgemm) ( &transb,&transa, + &N, &M, &K, alpha, B, &ldb, A, &lda, beta, C, &ldc); +#endif +} + extern "C" void FORNAME(dgemv) (const char *TRANS, const FINT *M, const FINT *N, const double *ALPHA, const double *A, const FINT *LDA, const double *X, const FINT *INCX, const double *BETA, double *Y, const FINT *INCY); void cblas_dgemv(const CBLAS_ORDER Order, @@ -691,6 +937,50 @@ if(TransA == CblasConjTrans) #endif } +extern "C" void FORNAME(cgemv) (const char *TRANS, const FINT *M, const FINT *N, const void *ALPHA, const void *A, const FINT *LDA, const void *X, const FINT *INCX, const void *BETA, void *Y, const FINT *INCY); +void cblas_cgemv(const CBLAS_ORDER Order, + const CBLAS_TRANSPOSE TransA, const int M, const int N, + const void *alpha, const void *A, const int lda, + const void *X, const int incX, const void *beta, + void *Y, const int incY) +{ + if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); + //swap n-m and toggle transposition + char transa = BLAS_FORTRANCASE(TransA==CblasNoTrans?'T':'N'); +#ifdef FORINT + const FINT mtmp=M; + const FINT ntmp=N; + 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 cgemv 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(cgemv) (&transa, &ntmp, &mtmp, &alphac, A, &ldatmp, X, &incxtmp, &betac, Y, &incytmp ); +#else + FORNAME(cgemv) (&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(cgemv) (&transa, &ntmp, &mtmp, alpha, A, &ldatmp, X, &incxtmp, beta, Y, &incytmp ); +#else + FORNAME(cgemv) (&transa, &N, &M, alpha, A, &lda, X, &incX, beta, Y, &incY ); +#endif +} + extern "C" FINT FORNAME(idamax) (const FINT *N, const double *DX, const FINT *INCX); CBLAS_INDEX cblas_idamax(const int N, const double *X, const int incX) { #ifdef FORINT @@ -713,6 +1003,17 @@ CBLAS_INDEX cblas_izamax(const int N, const void *X, const int incX) { #endif } +extern "C" FINT FORNAME(icamax) (const FINT *N, const void *DX, const FINT *INCX); +CBLAS_INDEX cblas_icamax(const int N, const void *X, const int incX) { +#ifdef FORINT + const FINT ntmp=N; + const FINT incxtmp=incX; + return (CBLAS_INDEX)FORNAME(icamax)(&ntmp, X, &incxtmp)-1; +#else + return (CBLAS_INDEX)FORNAME(icamax)(&N, X, &incX)-1; +#endif +} + extern "C" FINT FORNAME(isamax) (const FINT *N, const float *DX, const FINT *INCX); CBLAS_INDEX cblas_isamax(const int N, const float *X, const int incX) { #ifdef FORINT @@ -788,7 +1089,6 @@ int clapack_sgesv(const CBLAS_ORDER Order, const int N, const int NRHS, float *A, const int lda, int *ipiv, float *B, const int ldb) { -std::cout <<"In my clapack_sgesv\n"; FINT INFO=0; if(Order!=CblasRowMajor) laerror("CblasRowMajor order asserted"); //B should be in the same physical order, just transpose A in place and the LU result on output @@ -848,4 +1148,16 @@ void cblas_zswap(const int N, void *X, const int incX, void *Y, const int incY){ #endif } +extern "C" void FORNAME(cswap)(const FINT *N, void *X, const FINT *incX, void *Y, const FINT *incY); +void cblas_cswap(const int N, void *X, const int incX, void *Y, const int incY){ +#ifdef FORINT + const FINT N_tmp = N; + const FINT incX_tmp = incX; + const FINT incY_tmp = incY; + FORNAME(cswap)(&N_tmp, X, &incX_tmp, Y, &incY_tmp); +#else + FORNAME(cswap)(&N, X, &incX, Y, &incY); +#endif +} + #endif