Mentions légales du service

Skip to content
Snippets Groups Projects
Commit caab7a08 authored by Abel Calluaud's avatar Abel Calluaud
Browse files

testing: Update testings to completely test the new lrmm cases

parent 3dc2547d
No related branches found
No related tags found
No related merge requests found
Pipeline #1159792 failed
......@@ -80,7 +80,7 @@ set(TESTS_SOURCES
testrpk_zge2lr_performance.c
testrpk_zge2lr_stability.c
testrpk_zge2lr_tests.c
testrpk_zgemm_tests.c
testrpk_zgemm.c
testrpk_zrradd_tests.c
testrpk_zlaset.c
)
......
......@@ -37,13 +37,16 @@ int testrpk_zcheck_rradd( rpk_ctx_t *lowrank,
const test_matrix_t *B,
test_matrix_t *C );
int testrpk_zcheck_gemm( rpk_ctx_t *lowrank,
rpk_compmeth_t compmeth,
rpk_int_t offx, rpk_int_t offy,
rpk_complex64_t alpha,
int testrpk_zcheck_gemm( rpk_ctx_t *lowrank,
rpk_compmeth_t compmeth,
rpk_int_t offx,
rpk_int_t offy,
rpk_complex64_t alpha,
rpk_trans_t transA,
const test_matrix_t *A,
rpk_trans_t transB,
const test_matrix_t *B,
rpk_complex64_t beta,
const test_matrix_t *C );
rpk_complex64_t beta,
const test_matrix_t *C );
#endif /* _testrpk_z_h_ */
......@@ -322,6 +322,30 @@ testrpk_zgenmat( int mode,
return 0;
}
int testrpk_zgenmat_hermitian(test_matrix_t *matrix) {
int rc = 0;
assert(matrix->m == matrix->n);
matrix->fr = malloc(matrix->ld * matrix->n * sizeof(rpk_complex64_t));
char dist = 'U';
char sym = 'H';
// Lower triangular
int kl = 0;
int ku = 0;
double *S = malloc(matrix->n * sizeof(double));
rpk_complex64_t *work = malloc(3 * matrix->n * sizeof(rpk_complex64_t));
rc = LAPACKE_zlatms_work( LAPACK_COL_MAJOR, matrix->m, matrix->n, dist, ISEED, sym, S, 1, 1.0, 1.0, kl, ku, 'N', matrix->fr, matrix->ld, work );
assert(rc == 0);
free(S);
free(work);
return rc;
}
/**
*******************************************************************************
*
......@@ -375,9 +399,22 @@ testrpk_zgenmat_comp( const rpk_ctx_t *ctx,
test_matrix_t *A )
{
A->fr = malloc( A->ld * A->n * sizeof(rpk_complex64_t) );
testrpk_zgenmat( mode, ctx->tolerance, threshold, A );
ctx->rpk_ge2lr( ctx, rpk_imin( A->m, A->n ),
A->m, A->n, A->fr, A->ld, &(A->lr) );
int rc = testrpk_zgenmat( mode, ctx->tolerance, threshold, A );
assert(rc == 0);
if (A->lrcomp && A->rk > 0) {
ctx->rpk_ge2lr(
ctx, rpk_imin( A->m, A->n ),
A->m, A->n, A->fr, A->ld, &(A->lr)
);
}
else {
A->lr.rk = -1;
A->lr.rkmax = A->m;
A->lr.u = malloc( A->ld * A->n * sizeof(rpk_complex64_t) );
A->lr.v = NULL;
memcpy( A->lr.u, A->fr, A->ld * A->n * sizeof(rpk_complex64_t) );
}
}
/**
......@@ -734,45 +771,53 @@ testrpk_zcheck_gemm( rpk_ctx_t *ctx,
rpk_int_t offx,
rpk_int_t offy,
rpk_complex64_t alpha,
rpk_trans_t transA,
const test_matrix_t *A,
rpk_trans_t transB,
const test_matrix_t *B,
rpk_complex64_t beta,
const test_matrix_t *C )
{
rpk_matrix_t lrC2;
rpk_complex64_t *Clr;
double norm_diff, res;
Clock timer;
int rc = 0;
rpk_matrix_t Cra;
rpk_complex64_t *Clr;
double norm_diff, res;
Clock timer;
int rc = 0;
/* Init lrC */
memset( &lrC2, 0, sizeof( rpk_matrix_t ) );
rpk_zgemm_t zgemm_params;
zgemm_params.transA = transA;
zgemm_params.transB = transB;
zgemm_params.M = (transA == RapackNoTrans)? A->m : A->n;
zgemm_params.N = (transB == RapackNoTrans)? B->n : B->m;
zgemm_params.K = (transA == RapackNoTrans)? A->n : A->m;
zgemm_params.Cm = C->m;
zgemm_params.Cn = C->n;
zgemm_params.offx = offx;
zgemm_params.offy = offy;
zgemm_params.alpha = alpha;
zgemm_params.A = &(A->lr);
zgemm_params.B = &(B->lr);
zgemm_params.beta = beta;
zgemm_params.C = &Cra;
zgemm_params.work = NULL;
zgemm_params.lwork = -1;
zgemm_params.lwused = -1;
zgemm_params.lock = NULL;
/* Init Cra */
memset( &Cra, 0, sizeof( rpk_matrix_t ) );
Clr = malloc( C->m * C->n * sizeof( rpk_complex64_t ) );
/* Backup C into Cra */
rpk_zlrcpy( RapackNoTrans, 1., C->m, C->n, &(C->lr), C->m, C->n, &Cra, 0, 0 );
/* Backup C into C2 */
rpk_zlrcpy( RapackNoTrans, 1., C->m, C->n, &(C->lr), C->m, C->n, &lrC2, 0, 0 );
int nopA = (transA == RapackNoTrans) ? A->n : A->m;
int mopB = (transB == RapackNoTrans) ? B->m : B->n;
assert(nopA == mopB);
/* Compute the low-rank matrix-matrix */
{
rpk_zgemm_t zgemm_params;
zgemm_params.transA = RapackNoTrans;
zgemm_params.transB = RapackConjTrans;
zgemm_params.M = A->m;
zgemm_params.N = B->m;
zgemm_params.K = A->n;
zgemm_params.Cm = C->m;
zgemm_params.Cn = C->n;
zgemm_params.offx = offx;
zgemm_params.offy = offy;
zgemm_params.alpha = alpha;
zgemm_params.A = &(A->lr);
zgemm_params.B = &(B->lr);
zgemm_params.beta = beta;
zgemm_params.C = &lrC2;
zgemm_params.work = NULL;
zgemm_params.lwork = -1;
zgemm_params.lwused = -1;
zgemm_params.lock = NULL;
timer = clockGetLocal();
rpkx_zgemm( ctx, &zgemm_params );
timer = clockGetLocal() - timer;
......@@ -782,22 +827,19 @@ testrpk_zcheck_gemm( rpk_ctx_t *ctx,
* Check || (A+B) - (c(A)+c(B)) || < tol * || A+B ||
* What we really want is || (A+B) - c(A+B) || < tol * || A+B ||
*/
Clr = malloc( C->m * C->n * sizeof( rpk_complex64_t ) );
/* Uncompress the low-rank sum */
rpk_zlr2ge( RapackNoTrans, C->m, C->n,
&lrC2, Clr, C->ld );
rpk_zlr2ge( RapackNoTrans, C->m, C->n, &Cra, Clr, C->ld );
/* Compute the diff */
rpk_zgeadd( RapackNoTrans, C->m, C->n,
-1., C->fr, C->ld,
1., Clr, C->ld );
-1., C->fr, C->ld, 1., Clr, C->ld );
norm_diff = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'f', C->m, C->n,
Clr, C->ld, NULL );
if ( ( C->lr.rk != 0.0 ) && ( ( A->lr.rk + B->lr.rk ) != 0 ) ) {
double errbound = ctx->use_reltol ? C->norm : 1.;
if ( ( Cra.rk != 0.0 ) && ( ( A->lr.rk + B->lr.rk ) != 0 ) ) {
double errbound = (ctx->use_reltol && C->norm > 0) ? C->norm : 1.;
res = norm_diff / ( ctx->tolerance * errbound );
}
else {
......@@ -806,10 +848,7 @@ testrpk_zcheck_gemm( rpk_ctx_t *ctx,
fprintf( stdout, "%7s %4d %e %e %e %e ",
rpk_compmeth_shnames[compmeth],
(int)lrC2.rk, clockVal(timer), C->norm, norm_diff, res );
free( Clr );
rpk_zlrfree( &lrC2 );
(int)Cra.rk, clockVal(timer), C->norm, norm_diff, res );
/* Check the correctness of the result */
if ( res > 10.0 ) {
......@@ -820,8 +859,12 @@ testrpk_zcheck_gemm( rpk_ctx_t *ctx,
fprintf( stdout, "SUCCESS\n" );
}
else {
fprintf( stdout, "FAILED(%d)\n", rc );
fprintf( stdout, "FAILED(%d)\n", rc);
exit(1);
}
free( Clr );
rpk_zlrfree( &Cra );
return rc;
}
/**
*
* @file testrpk_zgemm.c
*
* Tests and validate the rpk_zgemm() routine.
*
* @copyright 2015-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
* Univ. Bordeaux. All rights reserved.
*
* @version 1.0.0
* @author Gregoire Pichon
* @author Mathieu Faverge
* @author Abel Calluaud
* @date 2024-03-11
*
* @precisions normal z -> c d s
*
**/
#include "tests.h"
#include "testrpk_z.h"
#include <math.h>
#define PRINT_RES(_ret_) \
if ( _ret_ == -1 ) { \
printf( "UNDEFINED\n" ); \
} \
else if( _ret_ > 0 ) { \
printf( "FAILED(%d)\n", _ret_ ); \
err++; \
} \
else { \
printf( "SUCCESS\n" ); \
}
int main( int argc, char **argv )
{
rpk_int_t m, n, k, Cm, Cn, offx, offy, use_reltol;
double eps = LAPACKE_dlamch_work( 'e' );
double tolerance = sqrt( eps );
double threshold = tolerance * tolerance;
test_matrix_t A, B, C;
rpk_complex64_t *Cfr;
int nalphas = 3;
int nbetas = 3;
rpk_complex64_t alphas[3] = { 0.0, 1.0, 3.14159 };
rpk_complex64_t betas[3] = { 0.0, 1.0, 2.71828 };
rpk_ctx_t rpkctx = rpk_ctx_default;
rpkctx.use_reltol = 0;
rpkctx.tolerance = tolerance;
rpkctx.rpk_ge2lr = rpkx_zge2lr_svd;
rpkctx.rpk_rradd = rpkx_zrradd_svd;
int ranks[3], r, s, i, j, l, meth;
int err = 0;
int ret = 0;
int mode = 0;
int compressA;
int compressB;
int compressC;
int test_case = 0; // Test case number
rpk_complex64_t alpha;
rpk_complex64_t beta;
CBLAS_TRANSPOSE transA;
CBLAS_TRANSPOSE transB;
for (int ia=0; ia<nalphas; ia++) {
alpha = alphas[ia];
for (int ib=0; ib<nbetas; ib++) {
beta = betas[ib];
for(use_reltol=0; use_reltol<2; use_reltol++) {
rpkctx.use_reltol = use_reltol;
for (s=25; s<=50; s = 2*s) {
ranks[0] = s + 1;
ranks[1] = 16;
ranks[2] = 2;
m = s / 2;
n = s / 4;
k = s / 6;
offx = 0;
offy = 0;
Cm = 2 * m;
Cn = 2 * n;
for (transA = CblasNoTrans; transA <= CblasConjTrans; transA +=1 ) {
for (transB = CblasNoTrans; transB <= CblasConjTrans; transB +=1 ) {
for (compressA = 0; compressA < 2; ++compressA) {
for (compressB = 0; compressB < 2; ++compressB) {
for (compressC = 0; compressC < 2; ++compressC) {
/* Matrix A */
for (i=0; i<3; i++) {
A.m = (transA == CblasNoTrans)? m : k;
A.n = (transA == CblasNoTrans)? k : m;
A.ld = A.m;
A.rk = rpk_imin( A.m, A.n ) / ranks[i];
A.lrcomp = compressA;
testrpk_zgenmat_comp( &rpkctx, mode, threshold, &A );
/* Matrix B */
for (j=0; j<3; j++) {
B.m = (transB == CblasNoTrans)? k : n;
B.n = (transB == CblasNoTrans)? n : k;
B.ld = B.m;
B.rk = rpk_imin( B.m, B.n ) / ranks[j];
B.lrcomp = compressB;
testrpk_zgenmat_comp( &rpkctx, mode, threshold, &B );
/* Matrix C */
for (l=0; l<3; l++) {
C.m = Cm;
C.n = Cn;
C.ld = Cm;
C.rk = rpk_imin( C.m, C.n ) / ranks[l];
C.lrcomp = compressC;
testrpk_zgenmat_comp( &rpkctx, mode, threshold, &C );
++ test_case;
printf( " -- Test gemm %d Cm=%ld, Cn=%ld, m=%ld, n=%ld, k=%ld, alpha=%f, beta=%f, rA=%ld, rB=%ld, rC=%ld, cA=%d, cB=%d, cC=%d tA=%d, tB=%d (%s)\n",
test_case, (long)Cm, (long)Cn, (long)m, (long)n, (long)k,
(double)alpha, (double)beta,
(long)(A.rk), (long)(B.rk), (long)(C.rk),
(int)compressA, (int)compressB, (int)compressC,
(int)transA, (int)transB, use_reltol ? "Relative": "Absolute" );
/* Compute the full rank GEMM */
Cfr = C.fr;
cblas_zgemm( CblasColMajor, transA, transB,
m, n, k,
CBLAS_SADDR( alpha ), A.fr, A.ld,
B.fr, B.ld,
CBLAS_SADDR( beta ), Cfr + C.ld * offy + offx, C.ld );
/* Update the norm of C with the norm of the result */
C.norm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'f', Cm, Cn,
Cfr, C.ld, NULL );
fprintf( stdout, "%7s %4s %12s %12s %12s %12s\n",
"Method", "Rank", "Time", "||C||_f", "||c(C)-C||_f",
"||c(C)-C||_f/(||C||_f * eps)" );
ret = 0;
/* Let's test all methods we have */
for(meth=0; meth<RapackCompressMethodNbr; meth++)
{
rpkctx.rpk_ge2lr = rpk_ge2lr_functions[meth][RapackComplex64-2];
rpkctx.rpk_rradd = rpk_rradd_functions[meth][RapackComplex64-2];
r = testrpk_zcheck_gemm( &rpkctx, meth, offx, offy,
alpha, (rpk_trans_t)transA, &A, (rpk_trans_t)transB, &B, beta, &C );
ret += r * (1 << meth);
}
rpkctx.rpk_ge2lr = rpk_ge2lr_functions[RapackCompressMethodSVD][RapackComplex64-2];
rpkctx.rpk_rradd = rpk_rradd_functions[RapackCompressMethodSVD][RapackComplex64-2];
rpk_zlrfree( &(C.lr) );
free( C.fr );
PRINT_RES( ret );
}
rpk_zlrfree( &(B.lr) );
free( B.fr );
}
rpk_zlrfree( &(A.lr) );
free( A.fr );
}
}
}
}
}
}
}
}
}
}
(void) argc;
(void) argv;
if( err == 0 ) {
printf(" -- All tests PASSED --\n");
return EXIT_SUCCESS;
}
else {
printf(" -- %d tests FAILED --\n", err);
return EXIT_FAILURE;
}
}
/**
*
* @file testrpk_zgemm_tests.c
*
* Tests and validate the rpk_zgemm() routine.
*
* @copyright 2015-2024 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
* Univ. Bordeaux. All rights reserved.
*
* @version 1.0.0
* @author Gregoire Pichon
* @author Mathieu Faverge
* @author Abel Calluaud
* @date 2024-03-11
*
* @precisions normal z -> c d s
*
**/
#include "tests.h"
#include "testrpk_z.h"
#include <math.h>
#define PRINT_RES(_ret_) \
if ( _ret_ == -1 ) { \
printf( "UNDEFINED\n" ); \
} \
else if( _ret_ > 0 ) { \
printf( "FAILED(%d)\n", _ret_ ); \
err++; \
} \
else { \
printf( "SUCCESS\n" ); \
}
int main( int argc, char **argv )
{
rpk_int_t m, n, k, Cm, Cn, offx, offy, use_reltol;
double eps = LAPACKE_dlamch_work( 'e' );
double tolerance = sqrt( eps );
double threshold = tolerance * tolerance;
test_matrix_t A, B, C;
rpk_complex64_t *Cfr;
rpk_complex64_t alpha = -3.48;
rpk_complex64_t beta = 1.;
rpk_ctx_t ctx = rpk_ctx_default;
int ranks[3], r, s, i, j, l, meth;
int err = 0;
int ret = 0;
int mode = 0;
ctx.use_reltol = 0;
ctx.tolerance = tolerance;
ctx.rpk_ge2lr = rpkx_zge2lr_svd;
ctx.rpk_rradd = rpkx_zrradd_svd;
for ( use_reltol = 0; use_reltol < 2; use_reltol++ ) {
ctx.use_reltol = use_reltol;
for ( s = 100; s <= 200; s = 2 * s ) {
ranks[0] = s + 1;
ranks[1] = 16;
ranks[2] = 2;
m = s / 2;
n = s / 4;
k = s / 6;
offx = 1;
offy = 2;
Cm = s;
Cn = s;
/* Matrix A */
for ( i = 0; i < 3; i++ ) {
A.m = m;
A.n = k;
A.ld = m;
A.rk = rpk_imin( m, k ) / ranks[i];
testrpk_zgenmat_comp( &ctx, mode, threshold, &A );
/* Matrix B */
for ( j = 0; j < 3; j++ ) {
B.m = n;
B.n = k;
B.ld = n;
B.rk = rpk_imin( n, k ) / ranks[j];
testrpk_zgenmat_comp( &ctx, mode, threshold, &B );
/* Matrix C */
for ( l = 0; l < 3; l++ ) {
C.m = Cm;
C.n = Cn;
C.ld = Cm;
C.rk = rpk_imin( Cm, Cn ) / ranks[l];
testrpk_zgenmat_comp( &ctx, mode, threshold, &C );
printf( " -- Test LRMM Cm=%ld, Cn=%ld, m=%ld, n=%ld, k=%ld, rA=%ld, rB=%ld, rC=%ld (%s)\n",
(long)Cm, (long)Cn, (long)m, (long)n, (long)k,
(long)(A.lr.rk), (long)(B.lr.rk), (long)(C.lr.rk),
use_reltol ? "Relative": "Absolute" );
/* Compute the full rank GEMM */
Cfr = C.fr;
cblas_zgemm( CblasColMajor, CblasNoTrans, CblasConjTrans,
m, n, k,
CBLAS_SADDR( alpha ), A.fr, A.ld,
B.fr, B.ld,
CBLAS_SADDR( beta ), Cfr + C.ld * offy + offx, C.ld );
/* Update the norm of C with the norm of the result */
C.norm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'f', Cm, Cn,
Cfr, C.ld, NULL );
fprintf( stdout, "%7s %4s %12s %12s %12s %12s\n",
"Method", "Rank", "Time", "||C||_f", "||c(C)-C||_f",
"||c(C)-C||_f/(||C||_f * eps)" );
ret = 0;
/* Let's test all methods we have */
for ( meth = 0; meth < RapackCompressMethodNbr; meth++ ) {
ctx.rpk_ge2lr = rpk_ge2lr_functions[meth][RapackComplex64 - 2];
ctx.rpk_rradd = rpk_rradd_functions[meth][RapackComplex64 - 2];
r = testrpk_zcheck_gemm( &ctx, meth, offx, offy,
alpha, &A, &B, beta, &C );
ret += r * (1 << meth);
}
ctx.rpk_ge2lr = rpk_ge2lr_functions[RapackCompressMethodSVD][RapackComplex64-2];
ctx.rpk_rradd = rpk_rradd_functions[RapackCompressMethodSVD][RapackComplex64-2];
rpk_zlrfree( &(C.lr) );
free( C.fr );
PRINT_RES( ret );
}
rpk_zlrfree( &(B.lr) );
free( B.fr );
}
rpk_zlrfree( &(A.lr) );
free( A.fr );
}
}
}
(void)argc;
(void)argv;
if ( err == 0 ) {
printf( " -- All tests PASSED --\n" );
return EXIT_SUCCESS;
}
else {
printf( " -- %d tests FAILED --\n", err );
return EXIT_FAILURE;
}
}
......@@ -138,6 +138,7 @@ int z_rradd_long()
A.m = all_ma[ima];
A.n = all_na[ima];
A.ld = all_ma[ima];
A.lrcomp = 1;
for (ira=rpk_imax(0,ima-2); ira<nb_ra; ira++) {
A.rk = all_ra[ira];
......@@ -154,6 +155,7 @@ int z_rradd_long()
B.m = all_mb[imb];
B.n = all_nb[imb];
B.ld = all_mb[imb];
B.lrcomp = 1;
for (irb=rpk_imax(0,imb-2); irb<nb_rb; irb++) {
B.rk = all_rb[irb];
......@@ -254,10 +256,12 @@ int z_rradd_short()
A.n = all_na[i];
A.rk = all_ra[i];
A.ld = all_ma[i];
A.lrcomp = 1;
B.m = all_mb[i];
B.n = all_nb[i];
B.rk = all_rb[i];
B.ld = all_mb[i];
B.lrcomp = 1;
for (use_reltol=0; use_reltol < 2; use_reltol++ ) {
printf( " -- Test RRADD MA=LDA=%d, NA=%d, RA=%d, MB=LDB=%d, NB=%d, RB=%d, rkmax=%ld, %s\n",
......
......@@ -41,13 +41,14 @@ struct test_param_s {
typedef struct test_param_s test_param_t;
typedef struct test_matrix_s {
int m; /* Number of rows of the test matrix */
int n; /* Number of columns of the test matrix */
int ld; /* Leading dimension of the test matrix */
int rk; /* Required rank of the test matrix */
double norm; /* Norm of the full rank matrix */
void *fr; /* Full rank matrix */
rpk_matrix_t lr; /* Low rank matrix */
int m; /* Number of rows of the test matrix */
int n; /* Number of columns of the test matrix */
int ld; /* Leading dimension of the test matrix */
int rk; /* Required rank of the test matrix */
double norm; /* Norm of the full rank matrix */
void *fr; /* Full rank matrix */
int lrcomp; /* Indicates low rank matrix storage */
rpk_matrix_t lr; /* Low rank matrix */
} test_matrix_t;
void testGetOptions( int argc, char **argv, test_param_t *params, double eps );
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment