added missing float support in noncblas

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<float> FORNAME(cdotu) (const FINT *n, const void *x, const FINT *incx, const void *y, const FINT *incy);
+extern "C" std::complex<float> 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<float> *)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<float> *)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<float> *)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<float> *)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<double> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+=
 p= (std::complex<double> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> alphac = std::conj(*(const std::complex<float> *)alpha);
+std::complex<float> betac = std::conj(*(const std::complex<float> *)beta);
+std::complex<float> *p;
+p= (std::complex<float> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+= incX;}
+p= (std::complex<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+= incX;}
+p= (std::complex<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<double> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+=
 if(!same) {p= (std::complex<double> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> *p;
+p= (std::complex<float> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+= incX;}
+if(!same) {p= (std::complex<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> *) X; for(int i=0; i<N; ++i) {p->imag(-p->imag()); p+= incX;}
+if(!same) {p= (std::complex<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> *p;
+p= (std::complex<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> *) Y; for(int i=0; i<N; ++i) {p->imag(-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<float> alphac = std::conj(*(const std::complex<float> *)alpha);
+	std::complex<float> betac = std::conj(*(const std::complex<float> *)beta);
+	std::complex<float> *p;
+	p= (std::complex<float> *) X; for(int i=0; i<lx; ++i) {p->imag(-p->imag()); p+= incX;}
+	p= (std::complex<float> *) Y; for(int i=0; i<ly; ++i) {p->imag(-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<float> *) X; for(int i=0; i<lx; ++i) {p->imag(-p->imag()); p+= incX;}
+	p= (std::complex<float> *) Y; for(int i=0; i<ly; ++i) {p->imag(-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