← Back to branch summary

~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to src/base/cs_benchmark.c

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-24 00:00:08 UTC
  • mfrom: (6.1.9 sid)
  • Revision ID: package-import@ubuntu.com-20111124000008-2vo99e38267942q5
Tags: 2.1.0-3
Install a missing file

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/base/cs_benchmark.c
1
 
/*============================================================================
2
 
 *
3
 
 *     This file is part of the Code_Saturne Kernel, element of the
4
 
 *     Code_Saturne CFD tool.
5
 
 *
6
 
 *     Copyright (C) 1998-2009 EDF S.A., France
7
 
 *
8
 
 *     contact: saturne-support@edf.fr
9
 
 *
10
 
 *     The Code_Saturne Kernel is free software; you can redistribute it
11
 
 *     and/or modify it under the terms of the GNU General Public License
12
 
 *     as published by the Free Software Foundation; either version 2 of
13
 
 *     the License, or (at your option) any later version.
14
 
 *
15
 
 *     The Code_Saturne Kernel is distributed in the hope that it will be
16
 
 *     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
 
 *     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18
 
 *     GNU General Public License for more details.
19
 
 *
20
 
 *     You should have received a copy of the GNU General Public License
21
 
 *     along with the Code_Saturne Kernel; if not, write to the
22
 
 *     Free Software Foundation, Inc.,
23
 
 *     51 Franklin St, Fifth Floor,
24
 
 *     Boston, MA  02110-1301  USA
25
 
 *
26
 
 *============================================================================*/
27
 
 
28
 
/*============================================================================
29
 
 * Low-level operator benchmarking
30
 
 *============================================================================*/
31
 
 
32
 
#if defined(HAVE_CONFIG_H)
33
 
#include "cs_config.h"
34
 
#endif
35
 
 
36
 
/*----------------------------------------------------------------------------
37
 
 * Standard C library headers
38
 
 *----------------------------------------------------------------------------*/
39
 
 
40
 
#include <stdarg.h>
41
 
#include <stdio.h>
42
 
#include <stdlib.h>
43
 
#include <string.h>
44
 
#include <assert.h>
45
 
#include <math.h>
46
 
 
47
 
#if defined(__STDC_VERSION__)      /* size_t */
48
 
#if (__STDC_VERSION__ == 199901L)
49
 
#    include <stddef.h>
50
 
#  else
51
 
#    include <stdlib.h>
52
 
#  endif
53
 
#else
54
 
#include <stdlib.h>
55
 
#endif
56
 
 
57
 
#if defined(HAVE_MPI)
58
 
#include <mpi.h>
59
 
#endif
60
 
 
61
 
/*----------------------------------------------------------------------------
62
 
 * BFT library headers
63
 
 *----------------------------------------------------------------------------*/
64
 
 
65
 
#include <bft_mem.h>
66
 
#include <bft_error.h>
67
 
#include <bft_printf.h>
68
 
#include <bft_timer.h>
69
 
 
70
 
/*----------------------------------------------------------------------------
71
 
 * FVM library headers
72
 
 *----------------------------------------------------------------------------*/
73
 
 
74
 
/*----------------------------------------------------------------------------
75
 
 * Local headers
76
 
 *----------------------------------------------------------------------------*/
77
 
 
78
 
#include "cs_base.h"
79
 
#include "cs_blas.h"
80
 
#include "cs_halo.h"
81
 
#include "cs_mesh.h"
82
 
#include "cs_mesh_quantities.h"
83
 
#include "cs_matrix.h"
84
 
#include "cs_perio.h"
85
 
 
86
 
/*----------------------------------------------------------------------------
87
 
 *  Header for the current file
88
 
 *----------------------------------------------------------------------------*/
89
 
 
90
 
#include "cs_benchmark.h"
91
 
 
92
 
/*----------------------------------------------------------------------------*/
93
 
 
94
 
BEGIN_C_DECLS
95
 
 
96
 
/*=============================================================================
97
 
 * Local Macro Definitions
98
 
 *============================================================================*/
99
 
 
100
 
/*=============================================================================
101
 
 * Local Structure Definitions
102
 
 *============================================================================*/
103
 
 
104
 
/*============================================================================
105
 
 *  Global variables
106
 
 *============================================================================*/
107
 
 
108
 
/*============================================================================
109
 
 * Private function definitions
110
 
 *============================================================================*/
111
 
 
112
 
/*----------------------------------------------------------------------------
113
 
 * Start timer.
114
 
 *
115
 
 * parameters:
116
 
 *   wt       --> wall-clock time (start in, stop - start out)
117
 
 *   cpu      --> CPU time (start in, stop - start out)
118
 
 *----------------------------------------------------------------------------*/
119
 
 
120
 
static void
121
 
_timer_start(double  *wt,
122
 
             double  *cpu)
123
 
{
124
 
  *wt = bft_timer_wtime();
125
 
  *cpu = bft_timer_cpu_time();
126
 
}
127
 
 
128
 
/*----------------------------------------------------------------------------
129
 
 * Stop timer.
130
 
 *
131
 
 * parameters:
132
 
 *   n_runs     <-- Number of timing runs
133
 
 *   wt         <-> wall-clock time (start in, stop - start out)
134
 
 *   cpu        <-> CPU time (start in, stop - start out)
135
 
 *----------------------------------------------------------------------------*/
136
 
 
137
 
static void
138
 
_timer_stop(int      n_runs,
139
 
            double  *wt,
140
 
            double  *cpu)
141
 
{
142
 
  double wt_s, cpu_s;
143
 
 
144
 
  wt_s = bft_timer_wtime();
145
 
  cpu_s = bft_timer_cpu_time();
146
 
 
147
 
  *wt = (wt_s - *wt) / (double)n_runs;
148
 
  *cpu = (cpu_s - *cpu) / (double)n_runs;
149
 
 }
150
 
 
151
 
/*----------------------------------------------------------------------------
152
 
 * Print overhead.
153
 
 *
154
 
 * parameters:
155
 
 *   wt       <-- wall-clock time
156
 
 *   cpu      <-- CPU time
157
 
 *----------------------------------------------------------------------------*/
158
 
 
159
 
static void
160
 
_print_overhead(double  wt,
161
 
                double  cpu)
162
 
{
163
 
  if (cs_glob_n_ranks == 1)
164
 
    bft_printf(_("  Wall clock : %12.5e\n"
165
 
                 "  CPU :        %12.5e\n"),
166
 
               wt);
167
 
 
168
 
  else {
169
 
 
170
 
    double loc_count[2], glob_min[2], glob_max[2], cpu_tot;
171
 
 
172
 
    loc_count[0] = wt;
173
 
    loc_count[1] = cpu;
174
 
 
175
 
#if defined(HAVE_MPI)
176
 
    MPI_Allreduce(loc_count, glob_min, 2, MPI_DOUBLE, MPI_MIN,
177
 
                  cs_glob_mpi_comm);
178
 
    MPI_Allreduce(loc_count, glob_max, 2, MPI_DOUBLE, MPI_MAX,
179
 
                  cs_glob_mpi_comm);
180
 
    MPI_Allreduce(&cpu, &cpu_tot, 1, MPI_DOUBLE, MPI_SUM,
181
 
                  cs_glob_mpi_comm);
182
 
#else
183
 
    { /* We should never enter here unless we have an alternative to MPI */
184
 
      int i;
185
 
      for (i = 0; i < 3; i++) {
186
 
        glob_min[i] = loc_count[i];
187
 
        glob_max[i] = loc_count[i];
188
 
      }
189
 
      cpu_tot = cpu;
190
 
    }
191
 
#endif
192
 
 
193
 
    bft_printf(_("               Min          Max          Total\n"
194
 
                 "  Wall clock : %12.5e %12.5e\n"
195
 
                 "  CPU :        %12.5e %12.5e %12.5e\n"),
196
 
               glob_min[0], glob_max[0],
197
 
               glob_min[1], glob_max[1], cpu_tot);
198
 
  }
199
 
}
200
 
 
201
 
/*----------------------------------------------------------------------------
202
 
 * Count number of operations.
203
 
 *
204
 
 * parameters:
205
 
 *   n_ops        <-- Local number of operations
206
 
 *   n_ops_single <-- Single-processor equivalent number of operations
207
 
 *                    (without ghosts); ignored if 0
208
 
 *   wt           <-- wall-clock time
209
 
 *   cpu          <-- CPU time
210
 
 *----------------------------------------------------------------------------*/
211
 
 
212
 
static void
213
 
_print_stats(long    n_ops,
214
 
             long    n_ops_single,
215
 
             double  wt,
216
 
             double  cpu)
217
 
{
218
 
  double fm = 1.0 / (1.e9 * wt);
219
 
 
220
 
  if (cs_glob_n_ranks == 1)
221
 
    bft_printf(_("  N ops :      %12ld\n"
222
 
                 "  Wall clock : %12.5e\n"
223
 
                 "  CPU :        %12.5e\n"
224
 
                 "  GFLOPS :     %12.5e\n"),
225
 
               n_ops, wt, cpu, n_ops*fm);
226
 
 
227
 
  else {
228
 
 
229
 
    long n_ops_min, n_ops_max, n_ops_tot;
230
 
    double loc_count[3], glob_min[3], glob_max[3], cpu_tot, fmg;
231
 
 
232
 
    loc_count[0] = wt;
233
 
    loc_count[1] = cpu;
234
 
    loc_count[2] = n_ops*fm;
235
 
 
236
 
#if defined(HAVE_MPI)
237
 
 
238
 
    MPI_Allreduce(&n_ops, &n_ops_min, 1, MPI_LONG, MPI_MIN,
239
 
                  cs_glob_mpi_comm);
240
 
    MPI_Allreduce(&n_ops, &n_ops_max, 1, MPI_LONG, MPI_MAX,
241
 
                  cs_glob_mpi_comm);
242
 
    MPI_Allreduce(&n_ops, &n_ops_tot, 1, MPI_LONG, MPI_SUM,
243
 
                  cs_glob_mpi_comm);
244
 
 
245
 
    MPI_Allreduce(loc_count, glob_min, 3, MPI_DOUBLE, MPI_MIN,
246
 
                  cs_glob_mpi_comm);
247
 
    MPI_Allreduce(loc_count, glob_max, 3, MPI_DOUBLE, MPI_MAX,
248
 
                  cs_glob_mpi_comm);
249
 
    MPI_Allreduce(&cpu, &cpu_tot, 1, MPI_DOUBLE, MPI_SUM,
250
 
                  cs_glob_mpi_comm);
251
 
 
252
 
#else
253
 
    { /* We should never enter here unless we have an alternative to MPI */
254
 
      int i;
255
 
      n_ops_min = n_ops; n_ops_max = n_ops; n_ops_tot = n_ops;
256
 
      for (i = 0; i < 3; i++) {
257
 
        glob_min[i] = loc_count[i];
258
 
        glob_max[i] = loc_count[i];
259
 
      }
260
 
      cpu_tot = cpu;
261
 
    }
262
 
#endif
263
 
 
264
 
    fmg = 1.0 / (1.e9 * glob_max[0]); /* global flops multiplier */
265
 
 
266
 
    if (n_ops_single == 0)
267
 
      bft_printf
268
 
        (_("               Min          Max          Total\n"
269
 
           "  N ops :      %12ld %12ld %12ld\n"
270
 
           "  Wall clock : %12.5e %12.5e\n"
271
 
           "  CPU :        %12.5e %12.5e %12.5e\n"
272
 
           "  GFLOPS :     %12.5e %12.5e %12.5e\n"),
273
 
         n_ops_min, n_ops_max, n_ops_tot,
274
 
         glob_min[0], glob_max[0],
275
 
         glob_min[1], glob_max[1], cpu_tot,
276
 
         glob_min[2], glob_max[2], n_ops_tot*fmg);
277
 
 
278
 
    else
279
 
      bft_printf
280
 
        (_("               Min          Max          Total        Single\n"
281
 
           "  N ops :      %12ld %12ld %12ld %12ld\n"
282
 
           "  Wall clock : %12.5e %12.5e\n"
283
 
           "  CPU :        %12.5e %12.5e %12.5e\n"
284
 
           "  GFLOPS :     %12.5e %12.5e %12.5e %12.5e\n"),
285
 
         n_ops_min, n_ops_max, n_ops_tot, n_ops_single,
286
 
         glob_min[0], glob_max[0],
287
 
         glob_min[1], glob_max[1], cpu_tot,
288
 
         glob_min[2], glob_max[2], n_ops_tot*fmg, n_ops_single*fmg);
289
 
  }
290
 
}
291
 
 
292
 
/*----------------------------------------------------------------------------
293
 
 * Simple dot product.
294
 
 *
295
 
 * parameters:
296
 
 *   global          <-- 0 for local use, 1 for MPI sum
297
 
 *   n_runs          <-- number of operation runs
298
 
 *   n_cells         <-- number of cells (array size)
299
 
 *   x               <-- Vector
300
 
 *   y               <-- Vector
301
 
 *----------------------------------------------------------------------------*/
302
 
 
303
 
static void
304
 
_dot_product_1(int                  global,
305
 
               int                  n_runs,
306
 
               size_t               n_cells,
307
 
               const cs_real_t     *x,
308
 
               const cs_real_t     *y)
309
 
{
310
 
  double wt, cpu;
311
 
  int    run_id;
312
 
  long   n_ops;
313
 
  size_t ii;
314
 
 
315
 
  double test_sum_mult = 1.0/n_runs;
316
 
  double test_sum = 0.0;
317
 
  int _global = global;
318
 
 
319
 
  char type_name[] = "X.Y";
320
 
 
321
 
  if (n_runs < 1)
322
 
    return;
323
 
 
324
 
  if (x == y)
325
 
    type_name[2] = 'X';
326
 
 
327
 
  if (cs_glob_n_ranks == 1)
328
 
    _global = 0;
329
 
 
330
 
  n_ops = n_cells;
331
 
 
332
 
  /* First simple local x.x version */
333
 
 
334
 
  _timer_start(&wt, &cpu);
335
 
 
336
 
#if defined(HAVE_BLAS)
337
 
 
338
 
  test_sum = 0.0;
339
 
 
340
 
  for (run_id = 0; run_id < n_runs; run_id++) {
341
 
    double s1 = cblas_ddot(n_cells, x, 1, y, 1);
342
 
#if defined(HAVE_MPI)
343
 
    if (_global) {
344
 
      double s1_glob = 0.0;
345
 
      MPI_Allreduce(&s1, &s1_glob, 1, MPI_DOUBLE, MPI_SUM,
346
 
                    cs_glob_mpi_comm);
347
 
      s1 = s1_glob;
348
 
    }
349
 
#endif
350
 
    test_sum += test_sum_mult*s1;
351
 
  }
352
 
 
353
 
  _timer_stop(n_runs, &wt, &cpu);
354
 
 
355
 
  if (_global == 0)
356
 
    bft_printf(_("\n"
357
 
                 "Simple local dot product %s with BLAS\n"
358
 
                 "-------------------------------------\n"),
359
 
               type_name);
360
 
  else
361
 
    bft_printf(_("\n"
362
 
                 "Simple global dot product %s with BLAS\n"
363
 
                 "---------------------------------------\n"),
364
 
               type_name);
365
 
 
366
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
367
 
             n_runs, test_sum);
368
 
 
369
 
  _print_stats(n_ops, 0, wt, cpu);
370
 
 
371
 
#endif /* defined(HAVE_BLAS) */
372
 
 
373
 
  test_sum = 0.0;
374
 
 
375
 
  for (run_id = 0; run_id < n_runs; run_id++) {
376
 
    double s1 = 0.0;
377
 
    for (ii = 0; ii < n_cells; ii++)
378
 
      s1 += x[ii] * y[ii];
379
 
#if defined(HAVE_MPI)
380
 
    if (_global) {
381
 
      double s1_glob = 0.0;
382
 
      MPI_Allreduce(&s1, &s1_glob, 1, MPI_DOUBLE, MPI_SUM,
383
 
                    cs_glob_mpi_comm);
384
 
      s1 = s1_glob;
385
 
    }
386
 
#endif
387
 
    test_sum += test_sum_mult*s1;
388
 
  }
389
 
 
390
 
  _timer_stop(n_runs, &wt, &cpu);
391
 
 
392
 
  if (_global == 0)
393
 
    bft_printf(_("\n"
394
 
                 "Simple local dot product %s\n"
395
 
                 "---------------------------\n"),
396
 
               type_name);
397
 
  else
398
 
    bft_printf(_("\n"
399
 
                 "Simple global dot product %s\n"
400
 
                 "----------------------------\n"),
401
 
               type_name);
402
 
 
403
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
404
 
             n_runs, test_sum);
405
 
 
406
 
  _print_stats(n_ops, 0, wt, cpu);
407
 
 
408
 
 }
409
 
 
410
 
/*----------------------------------------------------------------------------
411
 
 * Double local dot product.
412
 
 *
413
 
 * parameters:
414
 
 *   n_runs          <-- number of operation runs
415
 
 *   n_cells         <-- number of cells (array size)
416
 
 *   x               <-- Vector
417
 
 *   y               <-- Vector
418
 
 *----------------------------------------------------------------------------*/
419
 
 
420
 
static void
421
 
_dot_product_2(int                  n_runs,
422
 
               size_t               n_cells,
423
 
               const cs_real_t     *x,
424
 
               const cs_real_t     *y)
425
 
{
426
 
  double wt, cpu;
427
 
  int    run_id;
428
 
  long   n_ops;
429
 
  size_t ii;
430
 
 
431
 
  double test_sum_mult = 1.0/n_runs;
432
 
  double test_sum = 0.0;
433
 
 
434
 
  if (n_runs < 1)
435
 
    return;
436
 
 
437
 
  n_ops = n_cells * 2;
438
 
 
439
 
  /* First simple local x.x version */
440
 
 
441
 
  _timer_start(&wt, &cpu);
442
 
 
443
 
#if defined(HAVE_BLAS)
444
 
 
445
 
  test_sum = 0.0;
446
 
 
447
 
  for (run_id = 0; run_id < n_runs; run_id++) {
448
 
    double s1 = cblas_ddot(n_cells, x, 1, x, 1);
449
 
    double s2 = cblas_ddot(n_cells, x, 1, y, 1);
450
 
    test_sum += test_sum_mult*(s1+s2);
451
 
  }
452
 
 
453
 
  _timer_stop(n_runs, &wt, &cpu);
454
 
 
455
 
  bft_printf(_("\n"
456
 
               "Double local dot product X.X, Y.Y with BLAS\n"
457
 
               "-------------------------------------------\n"));
458
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
459
 
             n_runs, test_sum);
460
 
 
461
 
  _print_stats(n_ops, 0, wt, cpu);
462
 
 
463
 
#endif /* defined(HAVE_BLAS) */
464
 
 
465
 
  test_sum = 0.0;
466
 
 
467
 
  for (run_id = 0; run_id < n_runs; run_id++) {
468
 
    double s1 = 0.0;
469
 
    double s2 = 0.0;
470
 
    for (ii = 0; ii < n_cells; ii++) {
471
 
      s1 += x[ii] * x[ii];
472
 
      s2 += x[ii] * y[ii];
473
 
    }
474
 
    test_sum += test_sum_mult*(s1+s2);
475
 
  }
476
 
 
477
 
  _timer_stop(n_runs, &wt, &cpu);
478
 
 
479
 
  bft_printf(_("\n"
480
 
               "Double local dot product X.X, Y.Y\n"
481
 
               "---------------------------------\n"));
482
 
 
483
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
484
 
             n_runs, test_sum);
485
 
 
486
 
  _print_stats(n_ops, 0, wt, cpu);
487
 
 
488
 
 }
489
 
 
490
 
/*----------------------------------------------------------------------------
491
 
 * y -> ax + y test
492
 
 *
493
 
 * parameters:
494
 
 *   n_runs        <-- number of operation runs
495
 
 *   n_cells       <-- number of cells (array size)
496
 
 *   x             <-- Vector
497
 
 *   y             <-> Vector
498
 
 *----------------------------------------------------------------------------*/
499
 
 
500
 
static void
501
 
_axpy_(int                n_runs,
502
 
       size_t             n_cells,
503
 
       const cs_real_t   *restrict x,
504
 
       cs_real_t         *restrict y)
505
 
{
506
 
  double wt, cpu;
507
 
  int    run_id;
508
 
  long   n_ops;
509
 
  size_t ii;
510
 
 
511
 
  double test_sum_mult = 1.0/n_runs;
512
 
  double test_sum = 0.0;
513
 
 
514
 
  if (n_runs < 1)
515
 
    return;
516
 
 
517
 
  n_ops = n_cells * 2;
518
 
 
519
 
  /* First simple local x.x version */
520
 
 
521
 
  for (ii = 0; ii < n_cells; ii++)
522
 
    y[ii] = 0.0;
523
 
 
524
 
  _timer_start(&wt, &cpu);
525
 
 
526
 
#if defined(HAVE_BLAS)
527
 
 
528
 
  test_sum = 0.0;
529
 
 
530
 
  for (run_id = 0; run_id < n_runs; run_id++) {
531
 
 
532
 
    cblas_daxpy(n_cells, test_sum_mult, x, 1, y, 1);
533
 
 
534
 
    test_sum += test_sum_mult*y[run_id%n_cells];
535
 
 
536
 
  }
537
 
 
538
 
  _timer_stop(n_runs, &wt, &cpu);
539
 
 
540
 
  bft_printf(_("\n"
541
 
               "Y <- aX + Y with BLAS\n"
542
 
               "---------------------\n"));
543
 
 
544
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
545
 
             n_runs, test_sum);
546
 
 
547
 
  _print_stats(n_ops, 0, wt, cpu);
548
 
 
549
 
#endif /* defined(HAVE_BLAS) */
550
 
 
551
 
  test_sum = 0.0;
552
 
 
553
 
  for (run_id = 0; run_id < n_runs; run_id++) {
554
 
 
555
 
    for (ii = 0; ii < n_cells; ii++)
556
 
      y[ii] += test_sum_mult * x[ii];
557
 
 
558
 
    test_sum += test_sum_mult*y[run_id%n_cells];
559
 
 
560
 
  }
561
 
 
562
 
  _timer_stop(n_runs, &wt, &cpu);
563
 
 
564
 
  bft_printf(_("\n"
565
 
               "Y <- aX + Y\n"
566
 
               "-----------\n"));
567
 
 
568
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
569
 
             n_runs, test_sum);
570
 
 
571
 
  _print_stats(n_ops, 0, wt, cpu);
572
 
 
573
 
}
574
 
 
575
 
/*----------------------------------------------------------------------------
576
 
 * Simple divisions on a vector.
577
 
 *
578
 
 * parameters:
579
 
 *   n_runs          <-- number of operation runs
580
 
 *   n_cells         <-- number of cells (array size)
581
 
 *----------------------------------------------------------------------------*/
582
 
 
583
 
static void
584
 
_division_test(int     n_runs,
585
 
               size_t  n_cells)
586
 
{
587
 
  double wt, cpu;
588
 
  int    run_id;
589
 
  long   n_ops;
590
 
  size_t ii;
591
 
 
592
 
  cs_real_t  *x = NULL, *y = NULL, *z = NULL;
593
 
 
594
 
  if (n_runs < 1)
595
 
    return;
596
 
 
597
 
  BFT_MALLOC(x, n_cells, cs_real_t);
598
 
  BFT_MALLOC(y, n_cells, cs_real_t);
599
 
  BFT_MALLOC(z, n_cells, cs_real_t);
600
 
 
601
 
  n_ops = n_cells;
602
 
 
603
 
  /* Division */
604
 
  /*----------*/
605
 
 
606
 
  for (ii = 0; ii < n_cells; ii++) {
607
 
    x[ii] = 2.0 + ii%3;
608
 
    y[ii] = 2.0 + (ii+1)%3;
609
 
  }
610
 
 
611
 
  /* Division of two vectors */
612
 
 
613
 
  _timer_start(&wt, &cpu);
614
 
 
615
 
  for (run_id = 0; run_id < n_runs; run_id++) {
616
 
 
617
 
    for (ii = 0; ii < n_cells; ii++)
618
 
      z[ii] = x[ii] / y[ii];
619
 
 
620
 
  }
621
 
 
622
 
  _timer_stop(n_runs, &wt, &cpu);
623
 
 
624
 
  bft_printf(_("\n"
625
 
               "Division z = x/y\n"
626
 
               "----------------\n"));
627
 
 
628
 
  _print_stats(n_ops, 0, wt, cpu);
629
 
 
630
 
  BFT_FREE(z);
631
 
 
632
 
  /* Copy inverse of a vector */
633
 
 
634
 
  _timer_start(&wt, &cpu);
635
 
 
636
 
  for (run_id = 0; run_id < n_runs; run_id++) {
637
 
 
638
 
    for (ii = 0; ii < n_cells; ii++)
639
 
      y[ii] = 1.0 / x[ii];
640
 
 
641
 
  }
642
 
 
643
 
  _timer_stop(n_runs, &wt, &cpu);
644
 
 
645
 
  bft_printf(_("\n"
646
 
               "Division y = 1/x\n"
647
 
               "----------------\n"));
648
 
 
649
 
  _print_stats(n_ops, 0, wt, cpu);
650
 
 
651
 
  BFT_FREE(y);
652
 
 
653
 
  /* Inverse of a vector */
654
 
 
655
 
  _timer_start(&wt, &cpu);
656
 
 
657
 
  for (run_id = 0; run_id < n_runs; run_id++) {
658
 
 
659
 
    for (ii = 0; ii < n_cells; ii++)
660
 
      x[ii] = 1.0 / x[ii];
661
 
 
662
 
  }
663
 
 
664
 
  _timer_stop(n_runs, &wt, &cpu);
665
 
 
666
 
  bft_printf(_("\n"
667
 
               "Division x <- 1/x\n"
668
 
               "-----------------\n"));
669
 
 
670
 
  _print_stats(n_ops, 0, wt, cpu);
671
 
 
672
 
  BFT_FREE(x);
673
 
 
674
 
 }
675
 
 
676
 
/*----------------------------------------------------------------------------
677
 
 * Simple square root on a vector.
678
 
 *
679
 
 * parameters:
680
 
 *   n_runs          <-- number of operation runs
681
 
 *   n_cells         <-- number of cells (array size)
682
 
 *----------------------------------------------------------------------------*/
683
 
 
684
 
static void
685
 
_sqrt_test(int     n_runs,
686
 
           size_t  n_cells)
687
 
{
688
 
  double wt, cpu;
689
 
  int    run_id;
690
 
  long   n_ops;
691
 
  size_t ii;
692
 
 
693
 
  cs_real_t  *x = NULL, *y = NULL;
694
 
 
695
 
  if (n_runs < 1)
696
 
    return;
697
 
 
698
 
  BFT_MALLOC(x, n_cells, cs_real_t);
699
 
  BFT_MALLOC(y, n_cells, cs_real_t);
700
 
 
701
 
  n_ops = n_cells;
702
 
 
703
 
  for (ii = 0; ii < n_cells; ii++)
704
 
    x[ii] = 2.0 + ii%3;
705
 
 
706
 
  /* Copy square root of a vector */
707
 
 
708
 
  _timer_start(&wt, &cpu);
709
 
 
710
 
  for (run_id = 0; run_id < n_runs; run_id++) {
711
 
 
712
 
    for (ii = 0; ii < n_cells; ii++)
713
 
      y[ii] = sqrt(x[ii]);
714
 
 
715
 
  }
716
 
 
717
 
  _timer_stop(n_runs, &wt, &cpu);
718
 
 
719
 
  bft_printf(_("\n"
720
 
               "y = sqrt(x)\n"
721
 
               "-----------\n"));
722
 
 
723
 
  _print_stats(n_ops, 0, wt, cpu);
724
 
 
725
 
  BFT_FREE(y);
726
 
 
727
 
  /* In place square root of a vector */
728
 
 
729
 
  _timer_start(&wt, &cpu);
730
 
 
731
 
  for (run_id = 0; run_id < n_runs; run_id++) {
732
 
 
733
 
    for (ii = 0; ii < n_cells; ii++)
734
 
      x[ii] = sqrt(x[ii]);
735
 
 
736
 
  }
737
 
 
738
 
  _timer_stop(n_runs, &wt, &cpu);
739
 
 
740
 
  bft_printf(_("\n"
741
 
               "x = sqrt(x)\n"
742
 
               "-----------\n"));
743
 
 
744
 
  _print_stats(n_ops, 0, wt, cpu);
745
 
 
746
 
  BFT_FREE(x);
747
 
 
748
 
}
749
 
 
750
 
/*----------------------------------------------------------------------------
751
 
 * Measure matrix creation + destruction related performance.
752
 
 *
753
 
 * parameters:
754
 
 *   n_runs      <-- number of operation runs
755
 
 *   type_name   <-- type name
756
 
 *   type        <-- matrix type
757
 
 *   symmetric   <-- symmetric structure (if available)
758
 
 *   n_cells     <-- number of local cells
759
 
 *   n_cells_ext <-- number of cells including ghost cells (array size)
760
 
 *   n_faces     <-- local number of internal faces
761
 
 *   cell_num    <-- global cell numbers (1 to n)
762
 
 *   face_cell   <-- face -> cells connectivity (1 to n)
763
 
 *   halo        <-- cell halo structure
764
 
 *   numbering   <-- vectorization or thread-related numbering info, or NULL
765
 
 *----------------------------------------------------------------------------*/
766
 
 
767
 
static void
768
 
_matrix_creation_test(int                    n_runs,
769
 
                      const char            *type_name,
770
 
                      cs_matrix_type_t       type,
771
 
                      cs_bool_t              symmetric,
772
 
                      cs_int_t               n_cells,
773
 
                      cs_int_t               n_cells_ext,
774
 
                      cs_int_t               n_faces,
775
 
                      const fvm_gnum_t      *cell_num,
776
 
                      const cs_int_t        *face_cell,
777
 
                      const cs_halo_t       *halo,
778
 
                      const cs_numbering_t  *numbering)
779
 
{
780
 
  double wt, cpu;
781
 
  int    run_id;
782
 
 
783
 
  cs_matrix_t  *m = NULL;
784
 
 
785
 
  if (n_runs < 1)
786
 
    return;
787
 
 
788
 
  /* First count creation/destruction overhead */
789
 
 
790
 
  _timer_start(&wt, &cpu);
791
 
 
792
 
  for (run_id = 0; run_id < n_runs; run_id++) {
793
 
    m = cs_matrix_create(type,
794
 
                         symmetric,
795
 
                         true,
796
 
                         false,
797
 
                         n_cells,
798
 
                         n_cells_ext,
799
 
                         n_faces,
800
 
                         cell_num,
801
 
                         face_cell,
802
 
                         halo,
803
 
                         numbering);
804
 
    cs_matrix_destroy(&m);
805
 
  }
806
 
 
807
 
  _timer_stop(n_runs, &wt, &cpu);
808
 
 
809
 
  bft_printf(_("\n"
810
 
               "Matrix construction / destruction (%s)\n"
811
 
               "---------------------------------\n"), _(type_name));
812
 
 
813
 
  bft_printf(_("  (calls: %d)\n"), n_runs);
814
 
 
815
 
  _print_overhead(wt, cpu);
816
 
 
817
 
}
818
 
 
819
 
/*----------------------------------------------------------------------------
820
 
 * Measure matrix assignment performance.
821
 
 *
822
 
 * parameters:
823
 
 *   n_runs      <-- number of operation runs
824
 
 *   type_name   <-- type name
825
 
 *   type        <-- matrix type
826
 
 *   sym_struct  <-- symmetric structure (if available)
827
 
 *   sym_coeffs  <-- symmetric coefficients
828
 
 *   n_cells     <-- number of local cells
829
 
 *   n_cells_ext <-- number of cells including ghost cells (array size)
830
 
 *   n_faces     <-- local number of internal faces
831
 
 *   cell_num    <-- global cell numbers (1 to n)
832
 
 *   face_cell   <-- face -> cells connectivity (1 to n)
833
 
 *   halo        <-- cell halo structure
834
 
 *   numbering   <-- vectorization or thread-related numbering info, or NULL
835
 
 *   da          <-- diagonal values
836
 
 *   xa          <-- extradiagonal values
837
 
 *----------------------------------------------------------------------------*/
838
 
 
839
 
static void
840
 
_matrix_assignment_test(int                    n_runs,
841
 
                        const char            *type_name,
842
 
                        cs_matrix_type_t       type,
843
 
                        cs_bool_t              sym_struct,
844
 
                        cs_bool_t              sym_coeffs,
845
 
                        cs_int_t               n_cells,
846
 
                        cs_int_t               n_cells_ext,
847
 
                        cs_int_t               n_faces,
848
 
                        const fvm_gnum_t      *cell_num,
849
 
                        const cs_int_t        *face_cell,
850
 
                        const cs_halo_t       *halo,
851
 
                        const cs_numbering_t  *numbering,
852
 
                        const cs_real_t       *restrict da,
853
 
                        const cs_real_t       *restrict xa)
854
 
{
855
 
  double wt, cpu;
856
 
  int    run_id;
857
 
 
858
 
  cs_matrix_t *m = NULL;
859
 
 
860
 
  if (n_runs < 1)
861
 
    return;
862
 
 
863
 
  /* Count assignment overhead */
864
 
 
865
 
  m = cs_matrix_create(type,
866
 
                       sym_struct,
867
 
                       true,
868
 
                       false,
869
 
                       n_cells,
870
 
                       n_cells_ext,
871
 
                       n_faces,
872
 
                       cell_num,
873
 
                       face_cell,
874
 
                       halo,
875
 
                       numbering);
876
 
 
877
 
  _timer_start(&wt, &cpu);
878
 
 
879
 
  for (run_id = 0; run_id < n_runs; run_id++) {
880
 
    cs_matrix_set_coefficients(m,
881
 
                               sym_coeffs,
882
 
                               da,
883
 
                               xa);
884
 
  }
885
 
 
886
 
  _timer_stop(n_runs, &wt, &cpu);
887
 
 
888
 
  cs_matrix_destroy(&m);
889
 
 
890
 
  bft_printf(_("\n"
891
 
               "Matrix value assignment (%s)\n"
892
 
               "-----------------------\n"), _(type_name));
893
 
 
894
 
  bft_printf(_("  (calls: %d)\n"), n_runs);
895
 
 
896
 
  _print_overhead(wt, cpu);
897
 
 
898
 
}
899
 
 
900
 
/*----------------------------------------------------------------------------
901
 
 * Measure matrix.vector product related performance.
902
 
 *
903
 
 * parameters:
904
 
 *   n_runs      <-- number of operation runs
905
 
 *   type_name   <-- type name
906
 
 *   type        <-- matrix type
907
 
 *   sym_struct  <-- symmetric structure (if available)
908
 
 *   sym_coeffs  <-- symmetric coefficients
909
 
 *   n_cells     <-- number of local cells
910
 
 *   n_cells_ext <-- number of cells including ghost cells (array size)
911
 
 *   n_faces     <-- local number of internal faces
912
 
 *   cell_num    <-- global cell numbers (1 to n)
913
 
 *   face_cell   <-- face -> cells connectivity (1 to n)
914
 
 *   halo        <-- cell halo structure
915
 
 *   numbering   <-- vectorization or thread-related numbering info, or NULL
916
 
 *   da          <-- diagonal values
917
 
 *   xa          <-- extradiagonal values
918
 
 *   x           <-> vector
919
 
 *   y           --> vector
920
 
 *----------------------------------------------------------------------------*/
921
 
 
922
 
static void
923
 
_matrix_vector_test(int                    n_runs,
924
 
                    const char            *type_name,
925
 
                    cs_matrix_type_t       type,
926
 
                    cs_bool_t              sym_struct,
927
 
                    cs_bool_t              sym_coeffs,
928
 
                    cs_int_t               n_cells,
929
 
                    cs_int_t               n_cells_ext,
930
 
                    cs_int_t               n_faces,
931
 
                    const fvm_gnum_t      *cell_num,
932
 
                    const cs_int_t        *face_cell,
933
 
                    const cs_halo_t       *halo,
934
 
                    const cs_numbering_t  *numbering,
935
 
                    const cs_real_t       *restrict da,
936
 
                    const cs_real_t       *restrict xa,
937
 
                    cs_real_t             *restrict x,
938
 
                    cs_real_t             *restrict y)
939
 
{
940
 
  cs_int_t ii;
941
 
  double wt, cpu;
942
 
  int    run_id;
943
 
  long   n_ops, n_ops_glob;
944
 
 
945
 
  double test_sum = 0.0;
946
 
  cs_matrix_t *m = NULL;
947
 
 
948
 
  if (n_runs < 1)
949
 
    return;
950
 
 
951
 
  n_ops = n_cells + n_faces*2;
952
 
 
953
 
  if (cs_glob_n_ranks == 1)
954
 
    n_ops_glob = n_ops;
955
 
  else
956
 
    n_ops_glob = (  cs_glob_mesh->n_g_cells
957
 
                  + cs_glob_mesh->n_g_i_faces*2);
958
 
 
959
 
  m = cs_matrix_create(type,
960
 
                       sym_struct,
961
 
                       true,
962
 
                       false,
963
 
                       n_cells,
964
 
                       n_cells_ext,
965
 
                       n_faces,
966
 
                       cell_num,
967
 
                       face_cell,
968
 
                       halo,
969
 
                       numbering);
970
 
 
971
 
  cs_matrix_set_coefficients(m,
972
 
                             sym_coeffs,
973
 
                             da,
974
 
                             xa);
975
 
 
976
 
  /* Matrix.vector product */
977
 
 
978
 
  _timer_start(&wt, &cpu);
979
 
 
980
 
  for (run_id = 0; run_id < n_runs; run_id++) {
981
 
    cs_matrix_vector_multiply(CS_PERIO_ROTA_COPY,
982
 
                              m,
983
 
                              x,
984
 
                              y);
985
 
    test_sum += y[n_cells-1];
986
 
#if 0
987
 
    for (ii = 0; ii < n_cells; ii++)
988
 
      bft_printf("y[%d] = %12.4f\n", ii, y[ii]);
989
 
#endif
990
 
  }
991
 
 
992
 
  _timer_stop(n_runs, &wt, &cpu);
993
 
 
994
 
  bft_printf(_("\n"
995
 
               "Matrix.vector product (%s)\n"
996
 
               "---------------------\n"), _(type_name));
997
 
 
998
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
999
 
             n_runs, test_sum);
1000
 
 
1001
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1002
 
 
1003
 
  /* Local timing in parallel mode */
1004
 
 
1005
 
  if (cs_glob_n_ranks > 1) {
1006
 
 
1007
 
    test_sum = 0.0;
1008
 
 
1009
 
    _timer_start(&wt, &cpu);
1010
 
 
1011
 
    for (run_id = 0; run_id < n_runs; run_id++) {
1012
 
      cs_matrix_vector_multiply_nosync(m,
1013
 
                                       x,
1014
 
                                       y);
1015
 
      test_sum += y[n_cells-1];
1016
 
    }
1017
 
 
1018
 
    _timer_stop(n_runs, &wt, &cpu);
1019
 
 
1020
 
    bft_printf(_("\n"
1021
 
                 "Local matrix.vector product (%s)\n"
1022
 
                 "---------------------------\n"), _(type_name));
1023
 
 
1024
 
    bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1025
 
               n_runs, test_sum);
1026
 
 
1027
 
    _print_stats(n_ops, n_ops_glob, wt, cpu);
1028
 
 
1029
 
  }
1030
 
 
1031
 
  /* Combined matrix.vector product: alpha.A.x + beta.y */
1032
 
 
1033
 
  test_sum = 0.0;
1034
 
  for (ii = 0; ii < n_cells_ext; y[ii++] = 0.0);
1035
 
 
1036
 
  _timer_start(&wt, &cpu);
1037
 
 
1038
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1039
 
    cs_matrix_alpha_a_x_p_beta_y(CS_PERIO_ROTA_COPY,
1040
 
                                 0.5,
1041
 
                                 0.0,
1042
 
                                 m,
1043
 
                                 x,
1044
 
                                 y);
1045
 
    test_sum += y[n_cells-1];
1046
 
  }
1047
 
 
1048
 
  _timer_stop(n_runs, &wt, &cpu);
1049
 
 
1050
 
  bft_printf(_("\n"
1051
 
               "Matrix.vector product alpha.A.x + beta.y (%s)\n"
1052
 
               "----------------------------------------\n"), _(type_name));
1053
 
 
1054
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1055
 
             n_runs, test_sum);
1056
 
 
1057
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1058
 
 
1059
 
  /* (Matrix - diagonal).vector product (with diagonal structure) */
1060
 
 
1061
 
  cs_matrix_set_coefficients(m,
1062
 
                             sym_coeffs,
1063
 
                             NULL,
1064
 
                             xa);
1065
 
 
1066
 
  test_sum = 0.0;
1067
 
 
1068
 
  _timer_start(&wt, &cpu);
1069
 
 
1070
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1071
 
    cs_matrix_vector_multiply(CS_PERIO_ROTA_COPY,
1072
 
                              m,
1073
 
                              x,
1074
 
                              y);
1075
 
    test_sum += y[n_cells-1];
1076
 
  }
1077
 
 
1078
 
  _timer_stop(n_runs, &wt, &cpu);
1079
 
 
1080
 
  bft_printf(_("\n"
1081
 
               "(Matrix-diagonal).vector product (%s)\n"
1082
 
               "--------------------------------\n"), _(type_name));
1083
 
 
1084
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1085
 
             n_runs, test_sum);
1086
 
 
1087
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1088
 
 
1089
 
  /* (Matrix - diagonal).vector product */
1090
 
 
1091
 
  cs_matrix_destroy(&m);
1092
 
 
1093
 
  n_ops = n_faces*2;
1094
 
 
1095
 
  if (cs_glob_n_ranks == 1)
1096
 
    n_ops_glob = n_ops;
1097
 
  else
1098
 
    n_ops_glob = (  cs_glob_mesh->n_g_cells
1099
 
                  + cs_glob_mesh->n_g_i_faces*2);
1100
 
 
1101
 
  m = cs_matrix_create(type,
1102
 
                       sym_struct,
1103
 
                       false,
1104
 
                       false,
1105
 
                       n_cells,
1106
 
                       n_cells_ext,
1107
 
                       n_faces,
1108
 
                       cell_num,
1109
 
                       face_cell,
1110
 
                       halo,
1111
 
                       numbering);
1112
 
 
1113
 
  cs_matrix_set_coefficients(m,
1114
 
                             sym_coeffs,
1115
 
                             NULL,
1116
 
                             xa);
1117
 
 
1118
 
  test_sum = 0.0;
1119
 
 
1120
 
  _timer_start(&wt, &cpu);
1121
 
 
1122
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1123
 
    cs_matrix_vector_multiply(CS_PERIO_ROTA_COPY,
1124
 
                              m,
1125
 
                              x,
1126
 
                              y);
1127
 
    test_sum += y[n_cells-1];
1128
 
  }
1129
 
 
1130
 
  _timer_stop(n_runs, &wt, &cpu);
1131
 
 
1132
 
  bft_printf(_("\n"
1133
 
               "(Matrix without diagonal).vector product (%s)\n"
1134
 
               "----------------------------------------\n"), _(type_name));
1135
 
 
1136
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1137
 
             n_runs, test_sum);
1138
 
 
1139
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1140
 
 
1141
 
  cs_matrix_destroy(&m);
1142
 
 
1143
 
}
1144
 
 
1145
 
/*----------------------------------------------------------------------------
1146
 
 * Measure matrix.vector product extradiagonal terms related performance
1147
 
 * (symmetric matrix case).
1148
 
 *
1149
 
 * parameters:
1150
 
 *   n_faces         <-- local number of internal faces
1151
 
 *   face_cell       <-- face -> cells connectivity (1 to n)
1152
 
 *   xa              <-- extradiagonal values
1153
 
 *   x               <-- vector
1154
 
 *   y               <-> vector
1155
 
 *----------------------------------------------------------------------------*/
1156
 
 
1157
 
static void
1158
 
_mat_vec_exdiag_native(cs_int_t             n_faces,
1159
 
                       const cs_int_t      *face_cell,
1160
 
                       const cs_real_t     *restrict xa,
1161
 
                       cs_real_t           *restrict x,
1162
 
                       cs_real_t           *restrict y)
1163
 
{
1164
 
  cs_int_t  ii, jj, face_id;
1165
 
 
1166
 
  /* Tell IBM compiler not to alias */
1167
 
#if defined(__xlc__)
1168
 
#pragma disjoint(*x, *y, *xa)
1169
 
#endif
1170
 
 
1171
 
  const cs_int_t *restrict face_cel_p = face_cell;
1172
 
 
1173
 
  for (face_id = 0; face_id < n_faces; face_id++) {
1174
 
    ii = *face_cel_p++ - 1;
1175
 
    jj = *face_cel_p++ - 1;
1176
 
    y[ii] += xa[face_id] * x[jj];
1177
 
    y[jj] += xa[face_id] * x[ii];
1178
 
  }
1179
 
}
1180
 
 
1181
 
/*----------------------------------------------------------------------------
1182
 
 * Measure matrix.vector product extradiagonal terms related performance
1183
 
 * (symmetric matrix case, variant 1).
1184
 
 *
1185
 
 * parameters:
1186
 
 *   n_faces         <-- local number of internal faces
1187
 
 *   face_cell       <-- face -> cells connectivity (1 to n)
1188
 
 *   xa              <-- extradiagonal values
1189
 
 *   x               <-- vector
1190
 
 *   y               <-> vector
1191
 
 *----------------------------------------------------------------------------*/
1192
 
 
1193
 
static void
1194
 
_mat_vec_exdiag_native_v1(cs_int_t             n_faces,
1195
 
                          const cs_int_t      *face_cell,
1196
 
                          const cs_real_t     *restrict xa,
1197
 
                          cs_real_t           *restrict x,
1198
 
                          cs_real_t           *restrict y)
1199
 
{
1200
 
  cs_int_t  ii, ii_prev, kk, face_id, kk_max;
1201
 
  cs_real_t y_it, y_it_prev;
1202
 
 
1203
 
  const int l1_cache_size = 508;
1204
 
 
1205
 
  /*
1206
 
   * 1/ Split y[ii] and y[jj] computation into 2 loops to remove compiler
1207
 
   *    data dependency assertion between y[ii] and y[jj].
1208
 
   * 2/ keep index (*face_cel_p) in L1 cache from y[ii] loop to y[jj] loop
1209
 
   *    and xa in L2 cache.
1210
 
   * 3/ break high frequency occurence of data dependency from one iteration
1211
 
   *    to another in y[ii] loop (nonzero matrix value on the same line ii).
1212
 
   */
1213
 
 
1214
 
  const cs_int_t *restrict face_cel_p = face_cell;
1215
 
 
1216
 
  for (face_id = 0;
1217
 
       face_id < n_faces;
1218
 
       face_id += l1_cache_size) {
1219
 
 
1220
 
    kk_max = CS_MIN((n_faces - face_id), l1_cache_size);
1221
 
 
1222
 
    /* sub-loop to compute y[ii] += xa[face_id] * x[jj] */
1223
 
 
1224
 
    ii = face_cel_p[0] - 1;
1225
 
    ii_prev = ii;
1226
 
    y_it_prev = y[ii_prev] + xa[face_id] * x[face_cel_p[1] - 1];
1227
 
 
1228
 
    for (kk = 1; kk < kk_max; ++kk) {
1229
 
      ii = face_cel_p[2*kk] - 1;
1230
 
      /* y[ii] += xa[face_id+kk] * x[jj]; */
1231
 
      if(ii == ii_prev) {
1232
 
        y_it = y_it_prev;
1233
 
      }
1234
 
      else {
1235
 
        y_it = y[ii];
1236
 
        y[ii_prev] = y_it_prev;
1237
 
      }
1238
 
      ii_prev = ii;
1239
 
      y_it_prev = y_it + xa[face_id+kk] * x[face_cel_p[2*kk+1] - 1];
1240
 
    }
1241
 
    y[ii] = y_it_prev;
1242
 
 
1243
 
    /* sub-loop to compute y[ii] += xa[face_id] * x[jj] */
1244
 
 
1245
 
    for (kk = 0; kk < kk_max; ++kk) {
1246
 
      y[face_cel_p[2*kk+1] - 1]
1247
 
        += xa[face_id+kk] * x[face_cel_p[2*kk] - 1];
1248
 
    }
1249
 
    face_cel_p += 2 * l1_cache_size;
1250
 
  }
1251
 
}
1252
 
 
1253
 
/*----------------------------------------------------------------------------
1254
 
 * Measure matrix.vector product extradiagonal terms related performance
1255
 
 * with contribution to face-based array instead of cell-based array
1256
 
 * (symmetric matrix case).
1257
 
 *
1258
 
 * parameters:
1259
 
 *   n_faces         <-- local number of internal faces
1260
 
 *   face_cell       <-- face -> cells connectivity (1 to n)
1261
 
 *   xa              <-- extradiagonal values
1262
 
 *   x               <-- vector
1263
 
 *   ya              <-> vector
1264
 
 *----------------------------------------------------------------------------*/
1265
 
 
1266
 
static void
1267
 
_mat_vec_exdiag_part_p1(cs_int_t             n_faces,
1268
 
                        const cs_int_t      *face_cell,
1269
 
                        const cs_real_t     *restrict xa,
1270
 
                        cs_real_t           *restrict x,
1271
 
                        cs_real_t           *restrict ya)
1272
 
{
1273
 
  cs_int_t  ii, jj, face_id;
1274
 
 
1275
 
  /* Tell IBM compiler not to alias */
1276
 
#if defined(__xlc__)
1277
 
#pragma disjoint(*x, *xa, *ya)
1278
 
#endif
1279
 
 
1280
 
  const cs_int_t *restrict face_cel_p = face_cell;
1281
 
 
1282
 
  for (face_id = 0; face_id < n_faces; face_id++) {
1283
 
    ii = *face_cel_p++ - 1;
1284
 
    jj = *face_cel_p++ - 1;
1285
 
    ya[face_id] += xa[face_id] * x[ii];
1286
 
    ya[face_id] += xa[face_id] * x[jj];
1287
 
  }
1288
 
}
1289
 
 
1290
 
/*----------------------------------------------------------------------------
1291
 
 * Measure matrix.vector product local extradiagonal part related performance.
1292
 
 *
1293
 
 * parameters:
1294
 
 *   n_runs      <-- number of operation runs
1295
 
 *   n_cells     <-- number of cells
1296
 
 *   n_cells_ext <-- number of cells including ghost cells (array size)
1297
 
 *   n_faces     <-- local number of internal faces
1298
 
 *   face_cell   <-- face -> cells connectivity (1 to n)
1299
 
 *   xa          <-- extradiagonal values
1300
 
 *   x           <-> vector
1301
 
 *   y           --> vector
1302
 
 *----------------------------------------------------------------------------*/
1303
 
 
1304
 
static void
1305
 
_sub_matrix_vector_test(int                  n_runs,
1306
 
                        cs_int_t             n_cells,
1307
 
                        cs_int_t             n_cells_ext,
1308
 
                        cs_int_t             n_faces,
1309
 
                        const cs_int_t      *face_cell,
1310
 
                        const cs_real_t     *restrict xa,
1311
 
                        cs_real_t           *restrict x,
1312
 
                        cs_real_t           *restrict y)
1313
 
{
1314
 
  cs_int_t  jj;
1315
 
  double wt, cpu;
1316
 
  int    run_id;
1317
 
  long   n_ops, n_ops_glob;
1318
 
  double *ya = NULL;
1319
 
 
1320
 
  double test_sum = 0.0;
1321
 
 
1322
 
  if (n_runs < 1)
1323
 
    return;
1324
 
 
1325
 
  n_ops = n_faces*2;
1326
 
 
1327
 
  if (cs_glob_n_ranks == 1)
1328
 
    n_ops_glob = n_ops;
1329
 
  else
1330
 
    n_ops_glob = (  cs_glob_mesh->n_g_cells
1331
 
                  + cs_glob_mesh->n_g_i_faces*2);
1332
 
 
1333
 
  for (jj = 0; jj < n_cells_ext; jj++)
1334
 
    y[jj] = 0.0;
1335
 
 
1336
 
  /* Matrix.vector product, variant 0 */
1337
 
 
1338
 
  _timer_start(&wt, &cpu);
1339
 
 
1340
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1341
 
    _mat_vec_exdiag_native(n_faces, face_cell, xa, x, y);
1342
 
    test_sum += y[n_cells-1];
1343
 
  }
1344
 
 
1345
 
  _timer_stop(n_runs, &wt, &cpu);
1346
 
 
1347
 
  bft_printf(_("\n"
1348
 
               "Matrix.vector product, extradiagonal part, variant 0\n"
1349
 
               "---------------------\n"));
1350
 
 
1351
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1352
 
             n_runs, test_sum);
1353
 
 
1354
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1355
 
 
1356
 
  for (jj = 0; jj < n_cells_ext; jj++)
1357
 
    y[jj] = 0.0;
1358
 
 
1359
 
  test_sum = 0.0;
1360
 
 
1361
 
  /* Matrix.vector product, variant 1 */
1362
 
 
1363
 
  _timer_start(&wt, &cpu);
1364
 
 
1365
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1366
 
    _mat_vec_exdiag_native_v1(n_faces, face_cell, xa, x, y);
1367
 
    test_sum += y[n_cells-1];
1368
 
  }
1369
 
 
1370
 
  _timer_stop(n_runs, &wt, &cpu);
1371
 
 
1372
 
  bft_printf(_("\n"
1373
 
               "Matrix.vector product, extradiagonal part, variant 1\n"
1374
 
               "---------------------\n"));
1375
 
 
1376
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1377
 
             n_runs, test_sum);
1378
 
 
1379
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1380
 
 
1381
 
  /* Matrix.vector product, contribute to faces only */
1382
 
 
1383
 
  BFT_MALLOC(ya, n_faces, cs_real_t);
1384
 
  for (jj = 0; jj < n_faces; jj++)
1385
 
    ya[jj] = 0.0;
1386
 
 
1387
 
  test_sum = 0.0;
1388
 
 
1389
 
  _timer_start(&wt, &cpu);
1390
 
 
1391
 
  for (run_id = 0; run_id < n_runs; run_id++) {
1392
 
    _mat_vec_exdiag_part_p1(n_faces, face_cell, xa, x, ya);
1393
 
    test_sum += y[n_cells-1];
1394
 
  }
1395
 
 
1396
 
  _timer_stop(n_runs, &wt, &cpu);
1397
 
 
1398
 
  BFT_FREE(ya);
1399
 
 
1400
 
  bft_printf(_("\n"
1401
 
               "Matrix.vector product, face values only\n"
1402
 
               "---------------------\n"));
1403
 
 
1404
 
  bft_printf(_("  (calls: %d;  test sum: %12.5f)\n"),
1405
 
             n_runs, test_sum);
1406
 
 
1407
 
  _print_stats(n_ops, n_ops_glob, wt, cpu);
1408
 
 
1409
 
}
1410
 
 
1411
 
/*============================================================================
1412
 
 * Public function definitions
1413
 
 *============================================================================*/
1414
 
 
1415
 
/*----------------------------------------------------------------------------
1416
 
 * Run simple benchmarks.
1417
 
 *
1418
 
 * parameters:
1419
 
 *   mpi_trace_mode <-- indicates if timing mode (0) or MPI trace-friendly
1420
 
 *                      mode (1) is to be used
1421
 
 *----------------------------------------------------------------------------*/
1422
 
 
1423
 
void
1424
 
cs_benchmark(int  mpi_trace_mode)
1425
 
{
1426
 
  /* Local variable definitions */
1427
 
  /*----------------------------*/
1428
 
 
1429
 
  int n_runs;
1430
 
  size_t ii;
1431
 
 
1432
 
  cs_real_t *x1 = NULL, *x2 = NULL;
1433
 
  cs_real_t *y1 = NULL, *y2 = NULL;
1434
 
  cs_real_t *da = NULL, *xa = NULL;
1435
 
 
1436
 
  const cs_mesh_t *mesh = cs_glob_mesh;
1437
 
  const cs_mesh_quantities_t *mesh_v = cs_glob_mesh_quantities;
1438
 
 
1439
 
  size_t n_cells = mesh->n_cells;
1440
 
  size_t n_cells_ext = mesh->n_cells_with_ghosts;
1441
 
  size_t n_faces = mesh->n_i_faces;
1442
 
 
1443
 
  /* Allocate and initialize  working arrays */
1444
 
  /*-----------------------------------------*/
1445
 
 
1446
 
  BFT_MALLOC(x1, n_cells_ext, cs_real_t);
1447
 
  BFT_MALLOC(x2, n_cells_ext, cs_real_t);
1448
 
 
1449
 
  for (ii = 0; ii < n_cells_ext; ii++) {
1450
 
    x1[ii] = mesh_v->cell_cen[ii*3];
1451
 
    x2[ii] = mesh_v->cell_cen[ii*3 + 1];
1452
 
  }
1453
 
 
1454
 
  if (CS_MEM_ALIGN > 0) {
1455
 
 
1456
 
    BFT_MEMALIGN(y1, CS_MEM_ALIGN, n_cells_ext, cs_real_t);
1457
 
    BFT_MEMALIGN(y2, CS_MEM_ALIGN, n_cells_ext, cs_real_t);
1458
 
 
1459
 
  }
1460
 
  else {
1461
 
 
1462
 
    BFT_MALLOC(y1, n_cells_ext, cs_real_t);
1463
 
    BFT_MALLOC(y2, n_cells_ext, cs_real_t);
1464
 
 
1465
 
  }
1466
 
 
1467
 
  BFT_MALLOC(da, n_cells_ext, cs_real_t);
1468
 
  BFT_MALLOC(xa, n_faces*2, cs_real_t);
1469
 
 
1470
 
  for (ii = 0; ii < n_cells_ext; ii++) {
1471
 
    da[ii] = 1.0 + ii/n_cells_ext;
1472
 
    xa[ii] = mesh_v->cell_cen[ii*3 + 1];
1473
 
  }
1474
 
 
1475
 
  for (ii = 0; ii < n_faces; ii++) {
1476
 
    xa[ii] = 0.5*(1.0 + ii/n_faces);
1477
 
    xa[ii + n_faces] = -0.5*(1.0 + ii/n_faces);
1478
 
  }
1479
 
 
1480
 
  /* Run tests */
1481
 
  /*-----------*/
1482
 
 
1483
 
  bft_printf(_("\n"
1484
 
               "Benchmark mode activated\n"
1485
 
               "========================\n"));
1486
 
 
1487
 
  /* Dot product test */
1488
 
  /*------------------*/
1489
 
 
1490
 
  n_runs = (mpi_trace_mode) ? 1 : 10000;
1491
 
 
1492
 
  _dot_product_1(0, n_runs, n_cells, x1, x1);
1493
 
  _dot_product_1(0, n_runs, n_cells, x1, x2);
1494
 
  _dot_product_2(n_runs, n_cells, x1, x2);
1495
 
  _axpy_(n_runs, n_cells, x1, y1);
1496
 
 
1497
 
  _division_test(n_runs/5, n_cells);
1498
 
  _sqrt_test(n_runs/10, n_cells);
1499
 
 
1500
 
#if defined(HAVE_MPI)
1501
 
 
1502
 
  if (cs_glob_n_ranks > 1) {
1503
 
 
1504
 
    n_runs = (mpi_trace_mode) ? 1 : 10000;
1505
 
 
1506
 
    _dot_product_1(1, n_runs, n_cells, x1, x1);
1507
 
    _dot_product_1(1, n_runs, n_cells, x1, x2);
1508
 
 
1509
 
  }
1510
 
 
1511
 
#endif /* HAVE_MPI */
1512
 
 
1513
 
  /* Matrix test */
1514
 
  /*-------------*/
1515
 
 
1516
 
  n_runs = (mpi_trace_mode) ? 0 : 500;
1517
 
 
1518
 
  if (!mpi_trace_mode) {
1519
 
 
1520
 
    /* Creation tests */
1521
 
 
1522
 
    n_runs = 2000;
1523
 
 
1524
 
    _matrix_creation_test(n_runs,
1525
 
                          _("native"),
1526
 
                          CS_MATRIX_NATIVE, false,
1527
 
                          n_cells, n_cells_ext, n_faces,
1528
 
                          mesh->global_cell_num, mesh->i_face_cells,
1529
 
                          mesh->halo,
1530
 
                          mesh->i_face_numbering);
1531
 
 
1532
 
    n_runs = 300;
1533
 
 
1534
 
    _matrix_creation_test(n_runs,
1535
 
                          _("CSR"),
1536
 
                          CS_MATRIX_CSR, false,
1537
 
                          n_cells, n_cells_ext, n_faces,
1538
 
                          mesh->global_cell_num, mesh->i_face_cells,
1539
 
                          mesh->halo,
1540
 
                          mesh->i_face_numbering);
1541
 
 
1542
 
    _matrix_creation_test(n_runs,
1543
 
                          _("CSR sym"),
1544
 
                          CS_MATRIX_CSR, true,
1545
 
                          n_cells, n_cells_ext, n_faces,
1546
 
                          mesh->global_cell_num, mesh->i_face_cells,
1547
 
                          mesh->halo,
1548
 
                          mesh->i_face_numbering);
1549
 
 
1550
 
    /* Assignment tests */
1551
 
 
1552
 
    n_runs = 2000;
1553
 
 
1554
 
    _matrix_assignment_test(n_runs,
1555
 
                            _("native"),
1556
 
                            CS_MATRIX_NATIVE, false, false,
1557
 
                            n_cells, n_cells_ext, n_faces,
1558
 
                            mesh->global_cell_num, mesh->i_face_cells,
1559
 
                            mesh->halo, mesh->i_face_numbering, da, xa);
1560
 
 
1561
 
    _matrix_assignment_test(n_runs,
1562
 
                            _("native, sym coeffs"),
1563
 
                            CS_MATRIX_NATIVE, false, true,
1564
 
                            n_cells, n_cells_ext, n_faces,
1565
 
                            mesh->global_cell_num, mesh->i_face_cells,
1566
 
                            mesh->halo, mesh->i_face_numbering, da, xa);
1567
 
 
1568
 
    n_runs = 300;
1569
 
 
1570
 
    _matrix_assignment_test(n_runs,
1571
 
                            _("CSR"),
1572
 
                            CS_MATRIX_CSR, false, false,
1573
 
                            n_cells, n_cells_ext, n_faces,
1574
 
                            mesh->global_cell_num, mesh->i_face_cells,
1575
 
                            mesh->halo, mesh->i_face_numbering, da, xa);
1576
 
 
1577
 
    _matrix_assignment_test(n_runs,
1578
 
                            _("CSR, sym coeffs"),
1579
 
                            CS_MATRIX_CSR, false, true,
1580
 
                            n_cells, n_cells_ext, n_faces,
1581
 
                            mesh->global_cell_num, mesh->i_face_cells,
1582
 
                            mesh->halo, mesh->i_face_numbering, da, xa);
1583
 
 
1584
 
    _matrix_assignment_test(n_runs,
1585
 
                            _("CSR_sym"),
1586
 
                            CS_MATRIX_CSR, true, true,
1587
 
                            n_cells, n_cells_ext, n_faces,
1588
 
                            mesh->global_cell_num, mesh->i_face_cells,
1589
 
                            mesh->halo, mesh->i_face_numbering, da, xa);
1590
 
 
1591
 
  }
1592
 
 
1593
 
  /* Matrix.vector tests */
1594
 
 
1595
 
  n_runs = (mpi_trace_mode) ? 1 : 1000;
1596
 
 
1597
 
  _matrix_vector_test(n_runs,
1598
 
                      _("native"),
1599
 
                      CS_MATRIX_NATIVE, false, false,
1600
 
                      n_cells, n_cells_ext, n_faces,
1601
 
                      mesh->global_cell_num, mesh->i_face_cells, mesh->halo,
1602
 
                      mesh->i_face_numbering, da, xa, x1, y1);
1603
 
 
1604
 
  _matrix_vector_test(n_runs,
1605
 
                      _("native, sym coeffs"),
1606
 
                      CS_MATRIX_NATIVE, false, true,
1607
 
                      n_cells, n_cells_ext, n_faces,
1608
 
                      mesh->global_cell_num, mesh->i_face_cells, mesh->halo,
1609
 
                      mesh->i_face_numbering, da, xa, x1, y1);
1610
 
 
1611
 
  _matrix_vector_test(n_runs,
1612
 
                      _("CSR"),
1613
 
                      CS_MATRIX_CSR, false, false,
1614
 
                      n_cells, n_cells_ext, n_faces,
1615
 
                      mesh->global_cell_num, mesh->i_face_cells, mesh->halo,
1616
 
                      mesh->i_face_numbering, da, xa, x1, y1);
1617
 
 
1618
 
  _matrix_vector_test(n_runs,
1619
 
                      _("CSR_sym"),
1620
 
                      CS_MATRIX_CSR, true, true,
1621
 
                      n_cells, n_cells_ext, n_faces,
1622
 
                      mesh->global_cell_num, mesh->i_face_cells, mesh->halo,
1623
 
                      mesh->i_face_numbering, da, xa, x1, y1);
1624
 
 
1625
 
  _sub_matrix_vector_test(n_runs,
1626
 
                          n_cells,
1627
 
                          n_cells_ext,
1628
 
                          n_faces,
1629
 
                          mesh->i_face_cells,
1630
 
                          xa,
1631
 
                          x1,
1632
 
                          y1);
1633
 
 
1634
 
  /* Free working arrays */
1635
 
  /*---------------------*/
1636
 
 
1637
 
  BFT_FREE(x1);
1638
 
  BFT_FREE(x2);
1639
 
 
1640
 
  BFT_FREE(y1);
1641
 
  BFT_FREE(y2);
1642
 
 
1643
 
  BFT_FREE(da);
1644
 
  BFT_FREE(xa);
1645
 
 
1646
 
}
1647
 
 
1648
 
/*----------------------------------------------------------------------------*/
1649
 
 
1650
 
END_C_DECLS