Pull request OpenMathLib#132: Add [z/c]gemm conj and no trans matrix …

…tests

Merge in PL/openblas from dev/k.zaytseva/LM-268 to dev-riscv

utest/CMakeLists.txt

@@ -77,6 +77,8 @@ ${DIR_EXT}/test_ztrmv.c
 ${DIR_EXT}/test_ctrmv.c
 $(DIR_EXT)/test_ztrsv.c
 $(DIR_EXT)/test_ctrsv.c
+$(DIR_EXT)/test_zgemm.c
+$(DIR_EXT)/test_cgemm.c
 )
 
 # crashing on travis cl with an error code suggesting resource not found

utest/Makefile

@@ -27,7 +27,7 @@ OBJS_EXT+=$(DIR_EXT)/test_zaxpyc.o $(DIR_EXT)/test_caxpyc.o $(DIR_EXT)/test_ztrs
 OBJS_EXT+=$(DIR_EXT)/test_cgemv_t.o $(DIR_EXT)/test_zgemv_t.o $(DIR_EXT)/test_cgemv_n.o $(DIR_EXT)/test_zgemv_n.o
 OBJS_EXT+=$(DIR_EXT)/test_crot.o $(DIR_EXT)/test_zrot.o $(DIR_EXT)/test_zgbmv.c $(DIR_EXT)/test_cgbmv.c
 OBJS_EXT+=$(DIR_EXT)/test_sgemmt.o $(DIR_EXT)/test_dgemmt.o $(DIR_EXT)/test_cgemmt.o $(DIR_EXT)/test_zgemmt.o
-OBJS_EXT+=$(DIR_EXT)/test_ztrmv.o $(DIR_EXT)/test_ctrmv.o
+OBJS_EXT+=$(DIR_EXT)/test_ztrmv.o $(DIR_EXT)/test_ctrmv.o $(DIR_EXT)/test_zgemm.o $(DIR_EXT)/test_cgemm.o
 
 ifneq ($(NO_LAPACK), 1)
 OBJS += test_potrs.o

utest/test_extensions/test_cgemm.c

@@ -0,0 +1,284 @@
+/*
+ * Copyright (C), 2023, KNS Group LLC (YADRO).
+ * All Rights Reserved.
+ *
+ * This software contains the intellectual property of YADRO
+ * or is licensed to YADRO from third parties. Use of this
+ * software and the intellectual property contained therein is expressly
+ * limited to the terms and conditions of the License Agreement under which
+ * it is provided by YADRO.
+ */
+
+#include "utest/openblas_utest.h"
+#include <cblas.h>
+
+#define DATASIZE 100
+#define INCREMENT 2
+
+struct DATA_CGEMM
+{
+	float a_test[DATASIZE * DATASIZE * 2];
+    float a_verify[DATASIZE * DATASIZE * 2];
+	float b_test[DATASIZE * DATASIZE * 2];
+    float b_verify[DATASIZE * DATASIZE * 2];
+    float c_test[DATASIZE * DATASIZE * 2];
+	float c_verify[DATASIZE * DATASIZE * 2];
+};
+#ifdef BUILD_COMPLEX
+static struct DATA_CGEMM data_cgemm;
+
+static void rand_generate(float *a, blasint n)
+{
+	blasint i;
+	for (i = 0; i < n; i++)
+		a[i] = (float)rand() / (float)RAND_MAX * 5.0f;
+}
+
+/**
+ * Find difference between two rectangle matrix
+ * return norm of differences
+ */
+float smatrix_difference(float *a, float *b, blasint cols, blasint rows, blasint ld)
+{
+    blasint i = 0;
+    blasint j = 0;
+    blasint inc = 1;
+    float norm = 0.0f;
+
+    float *a_ptr = a;
+    float *b_ptr = b;
+
+    for(i = 0; i < rows; i++)
+    {
+        for (j = 0; j < cols; j++) {
+            a_ptr[j] -= b_ptr[j];
+        }
+        norm += cblas_snrm2(cols, a_ptr, inc);
+
+        a_ptr += ld;
+        b_ptr += ld;
+    }
+    return norm/(float)(rows);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices
+ * and comparing it with the non-conjugate cgemm.
+ *
+ * param transa specifies op(A), the transposition (conjugation) operation applied to A
+ * param transb specifies op(B), the transposition (conjugation) operation applied to B
+ * param m specifies the number of rows of the matrix op(A) and of the matrix C
+ * param n specifies the number of columns of the matrix op(B) and the number of columns of the matrix C
+ * param k specifies the number of columns of the matrix op(A) and the number of rows of the matrix op(B)
+ * param alpha - scaling factor for the matrix-matrix product
+ * param lda - leading dimension of matrix A
+ * param ldb - leading dimension of matrix B
+ * param beta - scaling factor for matrix C
+ * param ldc - leading dimension of matrix C
+ * return norm of difference
+ */
+static float check_cgemm(char transa, char transb, blasint m, blasint n, blasint k, 
+                         float *alpha, blasint lda, blasint ldb, float *beta, blasint ldc)
+{
+	blasint i;
+	float alpha_conj[] = {1.0f, 0.0f}; 
+	char transa_verify = transa;
+    char transb_verify = transb;
+
+    int arows = k, acols = m;
+    int brows = n, bcols = k;
+
+    if (transa == 'T' || transa == 'C'){
+        arows = m; acols = k;
+    }
+
+    if (transb == 'T' || transb == 'C'){
+        brows = k; bcols = n;
+    }
+
+	rand_generate(data_cgemm.a_test, arows * lda * 2);
+	rand_generate(data_cgemm.b_test, brows * ldb * 2);
+    rand_generate(data_cgemm.c_test, n * ldc * 2);
+
+	for (i = 0; i < arows * lda * 2; i++)
+		data_cgemm.a_verify[i] = data_cgemm.a_test[i];
+
+	for (i = 0; i < brows * ldb * 2; i++)
+		data_cgemm.b_verify[i] = data_cgemm.b_test[i];
+
+    for (i = 0; i < n * ldc * 2; i++)
+		data_cgemm.c_verify[i] = data_cgemm.c_test[i];
+
+	if (transa == 'R'){
+		cblas_cimatcopy(CblasColMajor, CblasConjNoTrans, arows, acols, alpha_conj, data_cgemm.a_verify, lda, lda);
+		transa_verify = 'N';
+	}
+
+    if (transb == 'R'){
+		cblas_cimatcopy(CblasColMajor, CblasConjNoTrans, brows, bcols, alpha_conj, data_cgemm.b_verify, ldb, ldb);
+		transb_verify = 'N';
+	}
+
+	BLASFUNC(cgemm)(&transa_verify, &transb_verify, &m, &n, &k, alpha, data_cgemm.a_verify, &lda,
+	 				data_cgemm.b_verify, &ldb, beta, data_cgemm.c_verify, &ldc);
+
+	BLASFUNC(cgemm)(&transa, &transb, &m, &n, &k, alpha, data_cgemm.a_test, &lda,
+	 				data_cgemm.b_test, &ldb, beta, data_cgemm.c_test, &ldc);
+
+	return smatrix_difference(data_cgemm.c_test, data_cgemm.c_verify, m, n, ldc*2);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is conjugate and transposed
+ * matrix B is conjugate and not transposed
+ */
+CTEST(cgemm, conjtransa_conjnotransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'C';
+	char transb = 'R';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is not conjugate and not transposed
+ * matrix B is conjugate and not transposed
+ */
+CTEST(cgemm, notransa_conjnotransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'N';
+	char transb = 'R';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is conjugate and not transposed
+ * matrix B is conjugate and transposed
+ */
+CTEST(cgemm, conjnotransa_conjtransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'R';
+	char transb = 'C';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is conjugate and not transposed
+ * matrix B is not conjugate and not transposed
+ */
+CTEST(cgemm, conjnotransa_notransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'R';
+	char transb = 'N';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is conjugate and not transposed
+ * matrix B is conjugate and not transposed
+ */
+CTEST(cgemm, conjnotransa_conjnotransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'R';
+	char transb = 'R';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is conjugate and not transposed
+ * matrix B is transposed
+ */
+CTEST(cgemm, conjnotransa_transb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'R';
+	char transb = 'T';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+
+/**
+ * Test cgemm with the conjugate matrices by conjugating and not transposed matrices 
+ * and comparing it with the non-conjugate cgemm.
+ * Test with the following options:
+ *
+ * matrix A is transposed
+ * matrix B is conjugate and not transposed
+ */
+CTEST(cgemm, transa_conjnotransb)
+{
+	blasint n = DATASIZE, m = DATASIZE, k = DATASIZE;
+    blasint lda = DATASIZE, ldb = DATASIZE, ldc = DATASIZE;
+	char transa = 'T';
+	char transb = 'R';
+	float alpha[] = {-2.0, 1.0f};
+    float beta[] = {1.0f, -1.0f};
+
+	float norm = check_cgemm(transa, transb, m, n, k, alpha, lda, ldb, beta, ldc);
+
+	ASSERT_DBL_NEAR_TOL(0.0f, norm, SINGLE_EPS);
+}
+#endif