blis/testsuite/src/test_hemv.c

   1 /*
   2
   3    BLIS
   4    An object-based framework for developing high-performance BLAS-like
   5    libraries.
   6
   7    Copyright (C) 2014, The University of Texas at Austin
   8
   9    Redistribution and use in source and binary forms, with or without
  10    modification, are permitted provided that the following conditions are
  11    met:
  12     - Redistributions of source code must retain the above copyright
  13       notice, this list of conditions and the following disclaimer.
  14     - Redistributions in binary form must reproduce the above copyright
  15       notice, this list of conditions and the following disclaimer in the
  16       documentation and/or other materials provided with the distribution.
  17     - Neither the name of The University of Texas at Austin nor the names
  18       of its contributors may be used to endorse or promote products
  19       derived from this software without specific prior written permission.
  20
  21    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  22    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  23    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  24    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  25    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  26    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  27    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  28    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  29    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  30    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  31    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  32
  33 */
  34
  35 #include "blis.h"
  36 #include "test_libblis.h"
  37
  38
  39 // Static variables.
  40 static char*     op_str                    = "hemv";
  41 static char*     o_types                   = "mvv";  // a x y
  42 static char*     p_types                   = "ucc";  // uploa conja conjx
  43 static thresh_t  thresh[BLIS_NUM_FP_TYPES] = { { 1e-04, 1e-05 },   // warn, pass for s
  44                                                { 1e-04, 1e-05 },   // warn, pass for c
  45                                                { 1e-13, 1e-14 },   // warn, pass for d
  46                                                { 1e-13, 1e-14 } }; // warn, pass for z
  47
  48 // Local prototypes.
  49 void libblis_test_hemv_deps( test_params_t* params,
  50                              test_op_t*     op );
  51
  52 void libblis_test_hemv_experiment( test_params_t* params,
  53                                    test_op_t*     op,
  54                                    iface_t        iface,
  55                                    num_t          datatype,
  56                                    char*          pc_str,
  57                                    char*          sc_str,
  58                                    unsigned int   p_cur,
  59                                    perf_t*        perf,
  60                                    double*        resid );
  61
  62 void libblis_test_hemv_impl( iface_t   iface,
  63                              obj_t*    alpha,
  64                              obj_t*    a,
  65                              obj_t*    x,
  66                              obj_t*    beta,
  67                              obj_t*    y );
  68
  69 void libblis_test_hemv_check( obj_t*  alpha,
  70                               obj_t*  a,
  71                               obj_t*  x,
  72                               obj_t*  beta,
  73                               obj_t*  y,
  74                               obj_t*  y_orig,
  75                               double* resid );
  76
  77
  78
  79 void libblis_test_hemv_deps( test_params_t* params, test_op_t* op )
  80 {
  81         libblis_test_randv( params, &(op->ops->randv) );
  82         libblis_test_randm( params, &(op->ops->randm) );
  83         libblis_test_normfv( params, &(op->ops->normfv) );
  84         libblis_test_subv( params, &(op->ops->subv) );
  85         libblis_test_copyv( params, &(op->ops->copyv) );
  86         libblis_test_scalv( params, &(op->ops->scalv) );
  87         libblis_test_gemv( params, &(op->ops->gemv) );
  88 }
  89
  90
  91
  92 void libblis_test_hemv( test_params_t* params, test_op_t* op )
  93 {
  94
  95         // Return early if this test has already been done.
  96         if ( op->test_done == TRUE ) return;
  97
  98         // Return early if operation is disabled.
  99         if ( op->op_switch == DISABLE_ALL ||
 100              op->ops->l2_over == DISABLE_ALL ) return;
 101
 102         // Call dependencies first.
 103         if ( TRUE ) libblis_test_hemv_deps( params, op );
 104
 105         // Execute the test driver for each implementation requested.
 106         if ( op->front_seq == ENABLE )
 107         {
 108                 libblis_test_op_driver( params,
 109                                         op,
 110                                         BLIS_TEST_SEQ_FRONT_END,
 111                                         op_str,
 112                                         p_types,
 113                                         o_types,
 114                                         thresh,
 115                                         libblis_test_hemv_experiment );
 116         }
 117 }
 118
 119
 120
 121 void libblis_test_hemv_experiment( test_params_t* params,
 122                                    test_op_t*     op,
 123                                    iface_t        iface,
 124                                    num_t          datatype,
 125                                    char*          pc_str,
 126                                    char*          sc_str,
 127                                    unsigned int   p_cur,
 128                                    perf_t*        perf,
 129                                    double*        resid )
 130 {
 131         unsigned int n_repeats = params->n_repeats;
 132         unsigned int i;
 133
 134         double       time_min  = 1e9;
 135         double       time;
 136
 137         dim_t        m;
 138
 139         uplo_t       uploa;
 140         conj_t       conja;
 141         conj_t       conjx;
 142
 143         obj_t        kappa;
 144         obj_t        alpha, a, x, beta;
 145 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 146         dim_t            test_way = bli_read_nway_from_env( "BLIS_LB_NT" );
 147         obj_t        y[test_way];
 148         obj_t        y_save[test_way];
 149         double       resid_local[test_way];
 150         unsigned int j;
 151 #else
 152         obj_t        y, y_save;
 153 #endif
 154
 155
 156         // Map the dimension specifier to an actual dimension.
 157         m = libblis_test_get_dim_from_prob_size( op->dim_spec[0], p_cur );
 158
 159         // Map parameter characters to BLIS constants.
 160         bli_param_map_char_to_blis_uplo( pc_str[0], &uploa );
 161         bli_param_map_char_to_blis_conj( pc_str[1], &conja );
 162         bli_param_map_char_to_blis_conj( pc_str[2], &conjx );
 163
 164         // Create test scalars.
 165         bli_obj_scalar_init_detached( datatype, &alpha );
 166         bli_obj_scalar_init_detached( datatype, &beta );
 167         bli_obj_scalar_init_detached( datatype, &kappa );
 168
 169         // Create test operands (vectors and/or matrices).
 170         libblis_test_mobj_create( params, datatype, BLIS_NO_TRANSPOSE,
 171                                   sc_str[0], m, m, &a );
 172         libblis_test_vobj_create( params, datatype,
 173                                   sc_str[1], m,    &x );
 174 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 175         for(i = 0; i < test_way; i++)
 176         {
 177                 libblis_test_vobj_create( params, datatype,
 178                                           sc_str[2], m,    &y[i] );
 179                 libblis_test_vobj_create( params, datatype,
 180                                           sc_str[2], m,    &y_save[i] );
 181         }
 182 #else
 183         libblis_test_vobj_create( params, datatype,
 184                                       sc_str[2], m,    &y );
 185         libblis_test_vobj_create( params, datatype,
 186                                       sc_str[2], m,    &y_save );
 187 #endif
 188
 189         // Set alpha and beta.
 190 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 191         if ( bli_obj_is_real( y[0] ) )
 192 #else
 193         if ( bli_obj_is_real( y ) )
 194 #endif
 195         {
 196                 bli_setsc(  1.0,  0.0, &alpha );
 197                 bli_setsc( -1.0,  0.0, &beta );
 198         }
 199         else
 200         {
 201                 bli_setsc(  0.0,  1.0, &alpha );
 202                 bli_setsc(  0.0, -1.0, &beta );
 203         }
 204
 205         // Set the structure and uplo properties of A.
 206         bli_obj_set_struc( BLIS_HERMITIAN, a );
 207         bli_obj_set_uplo( uploa, a );
 208
 209         // Randomize A, make it densely Hermitian, and zero the unstored triangle
 210         // to ensure the implementation reads only from the stored region.
 211         bli_randm( &a );
 212         bli_mkherm( &a );
 213         bli_mktrim( &a );
 214
 215         // Randomize x and y, and save y.
 216         bli_randv( &x );
 217 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 218         for(i = 0; i < test_way; i++)
 219         {
 220                 bli_randv( &y[i] );
 221                 bli_copyv( &y[i], &y_save[i] );
 222         }
 223 #else
 224         bli_randv( &y );
 225         bli_copyv( &y, &y_save );
 226 #endif
 227
 228         // Normalize vectors by m.
 229         bli_setsc( 1.0/( double )m, 0.0, &kappa );
 230         bli_scalv( &kappa, &x );
 231 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 232         for(i = 0; i < test_way; i++)
 233         {
 234                 bli_scalv( &kappa, &y[i] );
 235                 bli_scalv( &kappa, &y_save[i] );
 236
 237         }
 238 #else
 239         bli_scalv( &kappa, &y );
 240         bli_scalv( &kappa, &y_save );
 241 #endif
 242
 243         // Apply the remaining parameters.
 244         bli_obj_set_conj( conja, a );
 245         bli_obj_set_conj( conjx, x );
 246
 247         // Repeat the experiment n_repeats times and record results.
 248         for ( i = 0; i < n_repeats; ++i )
 249         {
 250 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 251                 // Need only one call to initialize the CBLAS OpenCL kernel
 252                 bli_copym( &y_save[0], &y[0] );
 253                 libblis_test_hemv_impl( iface, &alpha, &a, &x, &beta, &y[0] );
 254                 //but need to re-initialize C for each of iteration of n_repeats
 255                 for(j = 0; j < test_way; j++)
 256                 {
 257                         bli_copym( &y_save[j], &y[j] );
 258                 }
 259 #else
 260                 bli_copym( &y_save, &y );
 261                 libblis_test_hemv_impl( iface, &alpha, &a, &x, &beta, &y );
 262                 bli_copym( &y_save, &y );
 263 #endif
 264
 265                 time = bli_clock();
 266 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 267 #pragma omp parallel num_threads(test_way)
 268                 {
 269                         #pragma omp for
 270                         for(j = 0; j < test_way; j++)
 271                         {
 272                                 libblis_test_hemv_impl( iface, &alpha, &a, &x, &beta, &y[j] );
 273                         }
 274                 }
 275 #else
 276                 libblis_test_hemv_impl( iface, &alpha, &a, &x, &beta, &y );
 277 #endif
 278
 279                 time_min = bli_clock_min_diff( time_min, time );
 280         }
 281 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 282         // Estimate the performance of the best experiment repeat.
 283         perf->gflops = ( 1.0 * m * m ) *test_way / time_min / FLOPS_PER_UNIT_PERF;
 284         if ( bli_obj_is_complex( y[0] ) ) perf->gflops *= 4.0;
 285 #else
 286         // Estimate the performance of the best experiment repeat.
 287         perf->gflops = ( 1.0 * m * m ) / time_min / FLOPS_PER_UNIT_PERF;
 288         if ( bli_obj_is_complex( y ) ) perf->gflops *= 4.0;
 289 #endif
 290         perf->time = time_min;
 291
 292 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 293         // Check output of each thread, and send max residue to main
 294         for(i = 0; i < test_way; i++)
 295         {
 296                 // Perform checks.
 297                 libblis_test_hemv_check( &alpha, &a, &x, &beta, &y[i], &y_save[i], &resid_local[i] );
 298
 299                 // Set the structure and uplo properties of A. Check manipulates these values
 300                 bli_obj_set_struc( BLIS_HERMITIAN, a );
 301                 bli_obj_set_uplo( uploa, a );
 302                 bli_obj_set_conj( conja, a );
 303
 304
 305                 // Zero out performance and residual if output vector is empty.
 306                 libblis_test_check_empty_problem( &y[i], perf, &resid_local[i] );
 307
 308                 if(i == 0)
 309                 {
 310                         *resid = resid_local[i];
 311                 }
 312                 else if (resid_local[i] > *resid)
 313                 {
 314                         *resid = resid_local[i];
 315                 }
 316         }
 317
 318 #else
 319         // Perform checks.
 320         libblis_test_hemv_check( &alpha, &a, &x, &beta, &y, &y_save, resid );
 321
 322         // Zero out performance and residual if output vector is empty.
 323         libblis_test_check_empty_problem( &y, perf, resid );
 324 #endif
 325
 326
 327         // Free the test objects.
 328         bli_obj_free( &a );
 329         bli_obj_free( &x );
 330 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 331         for(i = 0; i < test_way; i++)
 332         {
 333                 bli_obj_free( &y[i] );
 334                 bli_obj_free( &y_save[i] );
 335         }
 336 #else
 337         bli_obj_free( &y );
 338         bli_obj_free( &y_save );
 339 #endif
 340 }
 341
 342
 343
 344 void libblis_test_hemv_impl( iface_t   iface,
 345                              obj_t*    alpha,
 346                              obj_t*    a,
 347                              obj_t*    x,
 348                              obj_t*    beta,
 349                              obj_t*    y )
 350 {
 351         switch ( iface )
 352         {
 353                 case BLIS_TEST_SEQ_FRONT_END:
 354                         ;
 355 #ifdef CBLAS
 356                         enum CBLAS_ORDER cblas_order;
 357                         enum CBLAS_UPLO cblas_uplo;
 358                         int n;
 359                         int lda, incx, incy;
 360
 361                         void *cblas_alpha, *cblas_beta;
 362                         void *cblas_a, *cblas_x, *cblas_y;
 363
 364                         cblas_alpha = bli_obj_buffer( *alpha );
 365                         cblas_beta  = bli_obj_buffer( *beta );
 366                         cblas_a     = bli_obj_buffer( *a );
 367                         cblas_x     = bli_obj_buffer( *x );
 368                         cblas_y     = bli_obj_buffer( *y );
 369
 370                         //m     = bli_obj_length( *a );
 371                         n     = bli_obj_width( *a );
 372
 373
 374                         if(bli_obj_is_col_stored( *a ))
 375                         {
 376                                 cblas_order = CblasColMajor;
 377                                 lda = bli_obj_col_stride(*a);
 378                         }
 379                         else if(bli_obj_is_row_stored( *a ) )
 380                         {
 381                                 cblas_order = CblasRowMajor;
 382                                 lda = bli_obj_row_stride(*a);
 383                         }
 384                         else
 385                         {
 386                                 bli_check_error_code( BLIS_INVALID_DIM_STRIDE_COMBINATION );
 387                                 return;
 388                         }
 389
 390                         incx = bli_obj_vector_inc(*x);
 391                         incy = bli_obj_vector_inc(*y);
 392
 393                         if(bli_obj_is_upper( *a ))
 394                         {
 395                                 cblas_uplo = CblasUpper;
 396                         }
 397                         else if (bli_obj_is_lower( *a ))
 398                         {
 399                                 cblas_uplo = CblasLower;
 400                         }
 401                         else
 402                         {
 403                                 bli_check_error_code( BLIS_INVALID_UPLO );
 404                                 return;
 405                         }
 406
 407                         if (bli_obj_is_scomplex( *a ))
 408                         {
 409                                 cblas_chemv(cblas_order, cblas_uplo, n, cblas_alpha, cblas_a, lda, cblas_x, incx, cblas_beta, cblas_y, incy);
 410                         }
 411                         else if (bli_obj_is_dcomplex( *a ))
 412                         {
 413                                 cblas_zhemv(cblas_order, cblas_uplo, n, cblas_alpha, cblas_a, lda, cblas_x, incx, cblas_beta, cblas_y, incy);
 414                         }
 415                         else
 416                         {
 417                                 printf("CBLAS not supported\n");
 418                                 return;
 419                         }
 420
 421
 422
 423 #else
 424                 bli_hemv( alpha, a, x, beta, y );
 425 #endif
 426                 break;
 427
 428                 default:
 429                 libblis_test_printf_error( "Invalid interface type.\n" );
 430         }
 431 }
 432
 433
 434
 435 void libblis_test_hemv_check( obj_t*  alpha,
 436                               obj_t*  a,
 437                               obj_t*  x,
 438                               obj_t*  beta,
 439                               obj_t*  y,
 440                               obj_t*  y_orig,
 441                               double* resid )
 442 {
 443         num_t  dt      = bli_obj_datatype( *y );
 444         num_t  dt_real = bli_obj_datatype_proj_to_real( *y );
 445
 446         dim_t  m       = bli_obj_vector_dim( *y );
 447
 448         obj_t  v;
 449         obj_t  norm;
 450
 451         double junk;
 452
 453         //
 454         // Pre-conditions:
 455         // - a is randomized and Hermitian.
 456         // - x is randomized.
 457         // - y_orig is randomized.
 458         // Note:
 459         // - alpha and beta should have non-zero imaginary components in the
 460         //   complex cases in order to more fully exercise the implementation.
 461         //
 462         // Under these conditions, we assume that the implementation for
 463         //
 464         //   y := beta * y_orig + alpha * conja(A) * conjx(x)
 465         //
 466         // is functioning correctly if
 467         //
 468         //   normf( y - v )
 469         //
 470         // is negligible, where
 471         //
 472         //   v = beta * y_orig + alpha * conja(A_dense) * x
 473         //
 474
 475         bli_obj_scalar_init_detached( dt_real, &norm );
 476
 477         bli_obj_create( dt, m, 1, 0, 0, &v );
 478
 479         bli_copyv( y_orig, &v );
 480
 481         bli_mkherm( a );
 482         bli_obj_set_struc( BLIS_GENERAL, *a );
 483         bli_obj_set_uplo( BLIS_DENSE, *a );
 484
 485         bli_gemv( alpha, a, x, beta, &v );
 486
 487         bli_subv( &v, y );
 488         bli_normfv( y, &norm );
 489         bli_getsc( &norm, resid, &junk );
 490
 491         bli_obj_free( &v );
 492 }
 493