Add a test for the distributed random generation of RHS to make sure we always...

Add a test for the distributed random generation of RHS to make sure we always work with the same RHS

tests/CMakeLists.txt

@@ -52,6 +52,7 @@ if ( SPM_WITH_MPI )
   list( APPEND TESTS
     spm_scatter_gather_tests.c
     spm_dist_norm_tests.c
+    spm_dist_genrhs_tests.c
     )
 endif()
 
@@ -68,7 +69,8 @@ set( SPM_DOF_TESTS
   spm_dof_expand_tests spm_dof_norm_tests spm_dof_matvec_tests)
 set( SPM_MPI_TESTS
   spm_scatter_gather_tests
-  spm_dist_norm_tests )
+  spm_dist_norm_tests
+  spm_dist_genrhs_tests )
 
 # List of run types
 set( RUNTYPE shm )

tests/spm_dist_genrhs_tests.c

+/**
+ *
+ * @file spm_dist_norm_tests.c
+ *
+ * Tests and validate the spm_genrhs routines in the case of random distributed vectors.
+ *
+ * @copyright 2015-2017 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
+ *                      Univ. Bordeaux. All rights reserved.
+ *
+ * @version 6.0.0
+ * @author Mathieu Faverge
+ * @author Theophile Terraz
+ * @author Tony Delarue
+ * @date 2020-05-19
+ *
+ **/
+#include <stdint.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <assert.h>
+#include <time.h>
+#include <spm_tests.h>
+
+#if !defined(SPM_WITH_MPI)
+#error "This test should not be compiled in non distributed version"
+#endif
+
+int main (int argc, char **argv)
+{
+    char         *filename;
+    spmatrix_t    original, *spm, *spmdist;
+    spm_driver_t  driver;
+    int           clustnbr = 1;
+    int           clustnum = 0;
+    int           baseval, root;
+    int           rc = SPM_SUCCESS;
+    int           err = 0;
+    int           dof, dofmax = 4;
+    int           to_free = 0;
+    int           nrhs = 3;
+    size_t        sizeloc, sizedst;
+    void         *bloc, *bdst;
+
+    MPI_Init( &argc, &argv );
+
+    /**
+     * Get options from command line
+     */
+    spmGetOptions( argc, argv,
+                   &driver, &filename );
+
+    rc = spmReadDriver( driver, filename, &original );
+    free(filename);
+
+    if ( rc != SPM_SUCCESS ) {
+        fprintf(stderr, "ERROR: Could not read the file, stop the test !!!\n");
+        return EXIT_FAILURE;
+    }
+
+    MPI_Comm_size( MPI_COMM_WORLD, &clustnbr );
+    MPI_Comm_rank( MPI_COMM_WORLD, &clustnum );
+
+    if ( original.flttype == SpmPattern ) {
+        spmGenFakeValues( &original );
+    }
+
+    spmPrintInfo( &original, stdout );
+
+    printf(" -- SPM GenRHS Test --\n");
+
+    for( dof=-1; dof<2; dof++ )
+    {
+        if ( dof >= 0 ) {
+            spm = spmDofExtend( &original, dof, dofmax );
+            to_free = 1;
+        }
+        else {
+            spm = malloc(sizeof(spmatrix_t));
+            memcpy( spm, &original, sizeof(spmatrix_t) );
+            to_free = 0;
+        }
+
+        if ( spm == NULL ) {
+            fprintf( stderr, "Issue to extend the matrix\n" );
+            continue;
+        }
+
+        sizeloc = spm_size_of( spm->flttype ) * spm->nexp * nrhs;
+        bloc    = malloc( sizeloc );
+
+        memset( bloc, 0xbeefdead, sizeloc );
+        if ( spmGenRHS( SpmRhsRndB, nrhs, spm,
+                        NULL, spm->nexp, bloc, spm->nexp ) != SPM_SUCCESS ) {
+            fprintf( stderr, "Issue to generate the local rhs\n" );
+            continue;
+        }
+
+        for( root=-1; root<clustnbr; root++ )
+        {
+            spmdist = spmScatter( spm, -1, NULL, 1, -1, MPI_COMM_WORLD );
+            if ( spmdist == NULL ) {
+                fprintf( stderr, "Failed to scatter the spm\n" );
+                err++;
+                continue;
+            }
+
+            sizedst = spm_size_of( spmdist->flttype ) * spmdist->nexp * nrhs;
+            bdst    = malloc( sizedst );
+
+            for( baseval=0; baseval<2; baseval++ )
+            {
+                spmBase( spmdist, baseval );
+
+                if(clustnum == 0) {
+                    printf( " Case: %s - base(%d) - dof(%s) - root(%d): ",
+                            fltnames[spmdist->flttype],
+                            baseval, dofname[dof+1], root );
+                }
+
+                memset( bdst, 0xdeadbeef, sizedst );
+                if ( spmGenRHS( SpmRhsRndB, nrhs, spmdist,
+                                NULL, spmdist->nexp, bdst, spmdist->nexp ) != SPM_SUCCESS ) {
+                    err++;
+                    continue;
+                }
+
+                switch( spmdist->flttype ){
+                case SpmComplex64:
+                    rc = z_spm_dist_genrhs_check( spmdist, nrhs, bloc, bdst, root );
+                    break;
+
+                case SpmComplex32:
+                    rc = c_spm_dist_genrhs_check( spmdist, nrhs, bloc, bdst, root );
+                    break;
+
+                case SpmFloat:
+                    rc = s_spm_dist_genrhs_check( spmdist, nrhs, bloc, bdst, root );
+                    break;
+
+                case SpmDouble:
+                default:
+                    rc = d_spm_dist_genrhs_check( spmdist, nrhs, bloc, bdst, root );
+                }
+
+                if ( clustnum == 0 ) {
+                    if ( rc == 0 ) {
+                        printf( "SUCCESS\n" );
+                    }
+                    else {
+                        printf( "FAILURE\n" );
+                    }
+                }
+                err = (rc != 0) ? err+1 : err;
+            }
+
+            free( bdst );
+            spmExit( spmdist );
+            free( spmdist );
+        }
+
+        if ( spm != &original ) {
+            if( to_free ){
+                spmExit( spm  );
+            }
+            free( spm );
+        }
+        free( bloc );
+    }
+
+    spmExit( &original  );
+    MPI_Finalize();
+
+    if( err == 0 ) {
+        if(clustnum == 0) {
+            printf(" -- All tests PASSED --\n");
+        }
+        return EXIT_SUCCESS;
+    }
+    else
+    {
+        printf(" -- %d tests FAILED --\n", err);
+        return EXIT_FAILURE;
+    }
+}

tests/spm_tests.h

@@ -83,21 +83,29 @@ void z_spm_print_check( char *filename, const spmatrix_t *spm );
 int  z_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  z_spm_norm_check( const spmatrix_t *spm );
 int  z_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
+int  z_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
+                              const spm_complex64_t *bloc, const spm_complex64_t *bdst, int root );
 
 void c_spm_print_check( char *filename, const spmatrix_t *spm );
 int  c_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  c_spm_norm_check( const spmatrix_t *spm );
 int  c_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
+int  c_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
+                              const spm_complex32_t *bloc, const spm_complex32_t *bdst, int root );
 
 void d_spm_print_check( char *filename, const spmatrix_t *spm );
 int  d_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  d_spm_norm_check( const spmatrix_t *spm );
 int  d_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
+int  d_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
+                              const double *bloc, const double *bdst, int root );
 
 void s_spm_print_check( char *filename, const spmatrix_t *spm );
 int  s_spm_matvec_check( spm_trans_t trans, const spmatrix_t *spm );
 int  s_spm_norm_check( const spmatrix_t *spm );
 int  s_spm_dist_norm_check( const spmatrix_t *spm, const spmatrix_t *spmdist );
+int  s_spm_dist_genrhs_check( const spmatrix_t *spm, spm_int_t nrhs,
+                              const float *bloc, const float *bdst, int root );
 
 void p_spm_print_check( char *filename, const spmatrix_t *spm );
 

tests/z_spm_tests.c

@@ -334,4 +334,55 @@ z_spm_dist_norm_check( const spmatrix_t *spm,
 
     return ret;
 }
+
+int
+z_spm_dist_genrhs_check( const spmatrix_t      *spm,
+                         spm_int_t              nrhs,
+                         const spm_complex64_t *bloc,
+                         const spm_complex64_t *bdst,
+                         int                    root )
+{
+    static spm_complex64_t mzone = (spm_complex64_t)-1.;
+    spm_complex64_t       *tmp;
+    double                 norm, normr, eps, result;
+    int                    rc, ret = 0;
+    spm_int_t              M;
+
+    eps  = LAPACKE_dlamch_work('e');
+    M    = spm->gNexp;
+    norm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'M', M, nrhs, bloc, M );
+
+    /* Let's gather the distributed RHS */
+    tmp = z_spmReduceRHS( spm, nrhs, bdst, spm->nexp, root );
+
+    rc = 0;
+    if ( tmp != NULL ) {
+        cblas_zaxpy( M * nrhs, CBLAS_SADDR(mzone), bloc, 1, tmp, 1 );
+        normr = LAPACKE_zlange( LAPACK_COL_MAJOR, 'M', M, nrhs, tmp, M );
+
+        result = fabs(normr) / (norm * eps);
+        /**
+         * By default the rhs is scaled by the frobenius norm of A, thus we need
+         * to take into account the accumulation error un distributed on the
+         * norm to validate the test here.
+         */
+        result /=  spm->gnnzexp;
+        if ( result > 1. ) {
+            fprintf( stderr, "[%2d] || X_global || = %e, || X_global -  X_dist || = %e, error=%e\n",
+                     (int)(spm->clustnum), norm, normr, result );
+            rc = 1;
+        }
+
+        free( tmp );
+    }
+    else {
+        if ( (spm->clustnum == root) || (root == -1) ) {
+            rc = 1;
+        }
+    }
+
+    MPI_Allreduce( &rc, &ret, 1, MPI_INT, MPI_MAX, spm->comm );
+
+    return ret;
+}
 #endif