blis/testsuite/src/test_her.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                    = "her";
  41 static char*     o_types                   = "vm";  // x a
  42 static char*     p_types                   = "uc";  // uploa 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_her_deps( test_params_t* params,
  50                             test_op_t*     op );
  51
  52 void libblis_test_her_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_her_impl( iface_t   iface,
  63                             obj_t*    alpha,
  64                             obj_t*    x,
  65                             obj_t*    a );
  66
  67 void libblis_test_her_check( obj_t*  alpha,
  68                              obj_t*  x,
  69                              obj_t*  a,
  70                              obj_t*  a_orig,
  71                              double* resid );
  72
  73
  74
  75 void libblis_test_her_deps( test_params_t* params, test_op_t* op )
  76 {
  77         libblis_test_randv( params, &(op->ops->randv) );
  78         libblis_test_randm( params, &(op->ops->randm) );
  79         libblis_test_normfv( params, &(op->ops->normfv) );
  80         libblis_test_subv( params, &(op->ops->subv) );
  81         libblis_test_copym( params, &(op->ops->copym) );
  82         libblis_test_scal2v( params, &(op->ops->scal2v) );
  83         libblis_test_dotv( params, &(op->ops->dotv) );
  84         libblis_test_gemv( params, &(op->ops->gemv) );
  85 }
  86
  87
  88
  89 void libblis_test_her( test_params_t* params, test_op_t* op )
  90 {
  91
  92         // Return early if this test has already been done.
  93         if ( op->test_done == TRUE ) return;
  94
  95         // Return early if operation is disabled.
  96         if ( op->op_switch == DISABLE_ALL ||
  97              op->ops->l2_over == DISABLE_ALL ) return;
  98
  99         // Call dependencies first.
 100         if ( TRUE ) libblis_test_her_deps( params, op );
 101
 102         // Execute the test driver for each implementation requested.
 103         if ( op->front_seq == ENABLE )
 104         {
 105                 libblis_test_op_driver( params,
 106                                         op,
 107                                         BLIS_TEST_SEQ_FRONT_END,
 108                                         op_str,
 109                                         p_types,
 110                                         o_types,
 111                                         thresh,
 112                                         libblis_test_her_experiment );
 113         }
 114 }
 115
 116
 117
 118 void libblis_test_her_experiment( test_params_t* params,
 119                                   test_op_t*     op,
 120                                   iface_t        iface,
 121                                   num_t          datatype,
 122                                   char*          pc_str,
 123                                   char*          sc_str,
 124                                   unsigned int   p_cur,
 125                                   perf_t*        perf,
 126                                   double*        resid )
 127 {
 128         unsigned int n_repeats = params->n_repeats;
 129         unsigned int i;
 130
 131         double       time_min  = 1e9;
 132         double       time;
 133
 134         dim_t        m;
 135
 136         uplo_t       uploa;
 137         conj_t       conjx;
 138
 139         obj_t        alpha, x;
 140 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 141         dim_t            test_way = bli_read_nway_from_env( "BLIS_LB_NT" );
 142         obj_t        a[test_way];
 143         obj_t        a_save[test_way];
 144         double       resid_local[test_way];
 145         unsigned int j;
 146 #else
 147         obj_t        a, a_save;
 148 #endif
 149
 150         // Map the dimension specifier to an actual dimension.
 151         m = libblis_test_get_dim_from_prob_size( op->dim_spec[0], p_cur );
 152
 153         // Map parameter characters to BLIS constants.
 154         bli_param_map_char_to_blis_uplo( pc_str[0], &uploa );
 155         bli_param_map_char_to_blis_conj( pc_str[1], &conjx );
 156
 157         // Create test scalars.
 158         bli_obj_scalar_init_detached( datatype, &alpha );
 159
 160         // Create test operands (vectors and/or matrices).
 161         libblis_test_vobj_create( params, datatype,
 162                                   sc_str[0], m,    &x );
 163 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 164         for(i = 0; i < test_way; i++)
 165         {
 166                 libblis_test_mobj_create( params, datatype, BLIS_NO_TRANSPOSE,
 167                                           sc_str[1], m, m, &a[i] );
 168                 libblis_test_mobj_create( params, datatype, BLIS_NO_TRANSPOSE,
 169                                           sc_str[1], m, m, &a_save[i] );
 170         }
 171 #else
 172         libblis_test_mobj_create( params, datatype, BLIS_NO_TRANSPOSE,
 173                                   sc_str[1], m, m, &a );
 174         libblis_test_mobj_create( params, datatype, BLIS_NO_TRANSPOSE,
 175                                   sc_str[1], m, m, &a_save );
 176 #endif
 177
 178         // Set alpha.
 179         //bli_copysc( &BLIS_MINUS_ONE, &alpha );
 180         bli_setsc( -1.0, 0.0, &alpha );
 181
 182         // Randomize x.
 183         bli_randv( &x );
 184
 185 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 186         for(i = 0; i < test_way; i++)
 187         {
 188                 // Set the structure and uplo properties of A.
 189                 bli_obj_set_struc( BLIS_HERMITIAN, a[i] );
 190                 bli_obj_set_uplo( uploa, a[i] );
 191
 192                 // Randomize A, make it densely Hermitian, and zero the unstored triangle
 193                 // to ensure the implementation is reads only from the stored region.
 194                 bli_randm( &a[i] );
 195                 bli_mkherm( &a[i] );
 196                 bli_mktrim( &a[i] );
 197
 198                 // Save A and set its structure and uplo properties.
 199                 bli_obj_set_struc( BLIS_HERMITIAN, a_save[i] );
 200                 bli_obj_set_uplo( uploa, a_save[i] );
 201                 bli_copym( &a[i], &a_save[i] );
 202
 203         }
 204 #else
 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 is reads only from the stored region.
 211         bli_randm( &a );
 212         bli_mkherm( &a );
 213         bli_mktrim( &a );
 214
 215         // Save A and set its structure and uplo properties.
 216         bli_obj_set_struc( BLIS_HERMITIAN, a_save );
 217         bli_obj_set_uplo( uploa, a_save );
 218         bli_copym( &a, &a_save );
 219 #endif
 220
 221         // Apply the remaining parameters.
 222         bli_obj_set_conj( conjx, x );
 223
 224         // Repeat the experiment n_repeats times and record results.
 225         for ( i = 0; i < n_repeats; ++i )
 226         {
 227 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 228                 // Need only one call to initialize the CBLAS OpenCL kernel
 229                 bli_copym( &a_save[0], &a[0] );
 230                 libblis_test_her_impl( iface, &alpha, &x, &a[0] );
 231
 232                 //but need to re-initialize C for each of iteration of n_repeats
 233                 for(j = 0; j < test_way; j++)
 234                 {
 235                         bli_copym( &a_save[i], &a[i] );
 236                 }
 237 #else
 238                 bli_copym( &a_save, &a );
 239                 libblis_test_her_impl( iface, &alpha, &x,  &a );
 240                 bli_copym( &a_save, &a );
 241 #endif
 242
 243                 time = bli_clock();
 244 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 245 #pragma omp parallel num_threads(test_way)
 246                 {
 247                         #pragma omp for
 248                         for(j = 0; j < test_way; j++)
 249                         {
 250                                 libblis_test_her_impl( iface, &alpha, &x, &a[j] );
 251                         }
 252                 }
 253 #else
 254                 libblis_test_her_impl( iface, &alpha, &x, &a );
 255 #endif
 256
 257                 time_min = bli_clock_min_diff( time_min, time );
 258         }
 259 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 260         // Estimate the performance of the best experiment repeat.
 261         perf->gflops = ( 1.0 * m * m ) * test_way / time_min / FLOPS_PER_UNIT_PERF;
 262         if ( bli_obj_is_complex( a[0]) ) perf->gflops *= 4.0;
 263 #else
 264         // Estimate the performance of the best experiment repeat.
 265         perf->gflops = ( 1.0 * m * m ) / time_min / FLOPS_PER_UNIT_PERF;
 266         if ( bli_obj_is_complex( a ) ) perf->gflops *= 4.0;
 267 #endif
 268         perf->time = time_min;
 269
 270 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 271         // Check output of each thread, and send max residue to main
 272         for(i = 0; i < test_way; i++)
 273         {
 274                 // Perform checks.
 275                 libblis_test_her_check( &alpha, &x, &a[i], &a_save[i], &resid_local[i] );
 276                 // Zero out performance and residual if output matrix is empty.
 277                 libblis_test_check_empty_problem( &a[i], perf, &resid_local[i] );
 278
 279                 if(i == 0)
 280                 {
 281                         *resid = resid_local[i];
 282                 }
 283                 else if (resid_local[i] > *resid)
 284                 {
 285                         *resid = resid_local[i];
 286                 }
 287         }
 288
 289 #else
 290         // Perform checks.
 291         libblis_test_her_check( &alpha, &x, &a, &a_save, resid );
 292
 293         // Zero out performance and residual if output matrix is empty.
 294         libblis_test_check_empty_problem( &a, perf, resid );
 295 #endif
 296
 297
 298
 299         // Free the test objects.
 300         bli_obj_free( &x );
 301 #ifdef BLIS_ENABLE_MULTITHREAD_TEST
 302         for(i = 0; i < test_way; i++)
 303         {
 304                 bli_obj_free( &a[i] );
 305                 bli_obj_free( &a_save[i] );
 306         }
 307 #else
 308         bli_obj_free( &a );
 309         bli_obj_free( &a_save );
 310 #endif
 311 }
 312
 313
 314
 315 void libblis_test_her_impl( iface_t   iface,
 316                             obj_t*    alpha,
 317                             obj_t*    x,
 318                             obj_t*    a )
 319 {
 320         switch ( iface )
 321         {
 322         case BLIS_TEST_SEQ_FRONT_END:
 323                 ;
 324 #ifdef CBLAS
 325                 enum CBLAS_ORDER cblas_order;
 326                 enum CBLAS_UPLO cblas_uplo;
 327                 int n;
 328                 int lda, incx;
 329
 330                 void *cblas_a;
 331                 void *cblas_x;
 332
 333                 cblas_a     = bli_obj_buffer( *a );
 334                 cblas_x     = bli_obj_buffer( *x );
 335
 336                 //m     = bli_obj_length( *a );
 337                 n     = bli_obj_width( *a );
 338
 339
 340                 if(bli_obj_is_col_stored(*a))
 341                 {
 342                         cblas_order = CblasColMajor;
 343                         lda = bli_obj_col_stride(*a);
 344                 }
 345                 else if(bli_obj_is_row_stored( *a ) )
 346                 {
 347                         cblas_order = CblasRowMajor;
 348                         lda = bli_obj_row_stride(*a);
 349                 }
 350                 else
 351                 {
 352                         bli_check_error_code( BLIS_INVALID_DIM_STRIDE_COMBINATION );
 353                         return;
 354                 }
 355
 356                 incx = bli_obj_vector_inc(*x);
 357
 358                 if(bli_obj_is_upper( *a ))
 359                 {
 360                         cblas_uplo = CblasUpper;
 361                 }
 362                 else if (bli_obj_is_lower( *a ))
 363                 {
 364                         cblas_uplo = CblasLower;
 365                 }
 366                 else
 367                 {
 368                         bli_check_error_code( BLIS_INVALID_UPLO );
 369                         return;
 370                 }
 371
 372                 if (bli_obj_is_scomplex( *a ))
 373                 {
 374                         float *cblas_alpha;
 375                         cblas_alpha = (float *) bli_obj_buffer( *alpha );
 376
 377
 378                         cblas_cher(cblas_order, cblas_uplo, n, *cblas_alpha, cblas_x, incx, cblas_a, lda);
 379                 }
 380                 else if (bli_obj_is_dcomplex( *a ))
 381                 {
 382                         double *cblas_alpha;
 383                         cblas_alpha = (double *) bli_obj_buffer( *alpha );
 384
 385                         cblas_zher(cblas_order, cblas_uplo, n, *cblas_alpha, cblas_x, incx, cblas_a, lda);
 386                 }
 387                 else
 388                 {
 389                         printf("CBLAS not supported\n");
 390                         libblis_test_printf_error( "Invalid interface type.\n" );
 391                         return;
 392                 }
 393
 394 #else
 395                 bli_her( alpha, x, a );
 396 #endif
 397                 break;
 398
 399                 default:
 400                 libblis_test_printf_error( "Invalid interface type.\n" );
 401         }
 402 }
 403
 404
 405
 406 void libblis_test_her_check( obj_t*  alpha,
 407                              obj_t*  x,
 408                              obj_t*  a,
 409                              obj_t*  a_orig,
 410                              double* resid )
 411 {
 412         num_t  dt      = bli_obj_datatype( *a );
 413         num_t  dt_real = bli_obj_datatype_proj_to_real( *a );
 414
 415         dim_t  m_a     = bli_obj_length( *a );
 416
 417         obj_t  xh, t, v, w;
 418         obj_t  tau, rho, norm;
 419
 420         double junk;
 421
 422         //
 423         // Pre-conditions:
 424         // - x is randomized.
 425         // - a is randomized and Hermitian.
 426         // Note:
 427         // - alpha must be real-valued.
 428         //
 429         // Under these conditions, we assume that the implementation for
 430         //
 431         //   A := A_orig + alpha * conjx(x) * conjx(x)^H
 432         //
 433         // is functioning correctly if
 434         //
 435         //   normf( v - w )
 436         //
 437         // is negligible, where
 438         //
 439         //   v = A * t
 440         //   w = ( A_orig + alpha * conjx(x) * conjx(x)^H ) * t
 441         //     =   A_orig * t + alpha * conjx(x) * conjx(x)^H * t
 442         //     =   A_orig * t + alpha * conjx(x) * rho
 443         //     =   A_orig * t + w
 444         //
 445
 446         bli_mkherm( a );
 447         bli_mkherm( a_orig );
 448         bli_obj_set_struc( BLIS_GENERAL, *a );
 449         bli_obj_set_struc( BLIS_GENERAL, *a_orig );
 450         bli_obj_set_uplo( BLIS_DENSE, *a );
 451         bli_obj_set_uplo( BLIS_DENSE, *a_orig );
 452
 453         bli_obj_scalar_init_detached( dt,      &tau );
 454         bli_obj_scalar_init_detached( dt,      &rho );
 455         bli_obj_scalar_init_detached( dt_real, &norm );
 456
 457         bli_obj_create( dt, m_a, 1, 0, 0, &t );
 458         bli_obj_create( dt, m_a, 1, 0, 0, &v );
 459         bli_obj_create( dt, m_a, 1, 0, 0, &w );
 460
 461         bli_obj_alias_with_conj( BLIS_CONJUGATE, *x, xh );
 462
 463         bli_setsc( 1.0/( double )m_a, -1.0/( double )m_a, &tau );
 464         bli_setv( &tau, &t );
 465
 466         bli_gemv( &BLIS_ONE, a, &t, &BLIS_ZERO, &v );
 467
 468         bli_dotv( &xh, &t, &rho );
 469         bli_mulsc( alpha, &rho );
 470         bli_scal2v( &rho, x, &w );
 471         bli_gemv( &BLIS_ONE, a_orig, &t, &BLIS_ONE, &w );
 472
 473         bli_subv( &w, &v );
 474         bli_normfv( &v, &norm );
 475         bli_getsc( &norm, resid, &junk );
 476
 477         bli_obj_free( &t );
 478         bli_obj_free( &v );
 479         bli_obj_free( &w );
 480 }
 481