|
| 1 | +/** |
| 2 | +* @license Apache-2.0 |
| 3 | +* |
| 4 | +* Copyright (c) 2026 The Stdlib Authors. |
| 5 | +* |
| 6 | +* Licensed under the Apache License, Version 2.0 (the "License"); |
| 7 | +* you may not use this file except in compliance with the License. |
| 8 | +* You may obtain a copy of the License at |
| 9 | +* |
| 10 | +* http://www.apache.org/licenses/LICENSE-2.0 |
| 11 | +* |
| 12 | +* Unless required by applicable law or agreed to in writing, software |
| 13 | +* distributed under the License is distributed on an "AS IS" BASIS, |
| 14 | +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| 15 | +* See the License for the specific language governing permissions and |
| 16 | +* limitations under the License. |
| 17 | +*/ |
| 18 | + |
| 19 | +#include "stdlib/blas/base/shared.h" |
| 20 | +#include "stdlib/blas/base/cgemv.h" |
| 21 | +#include <stdlib.h> |
| 22 | +#include <stdio.h> |
| 23 | +#include <math.h> |
| 24 | +#include <time.h> |
| 25 | +#include <sys/time.h> |
| 26 | +#include "stdlib/complex/float32/ctor.h" |
| 27 | + |
| 28 | +#define NAME "cgemv" |
| 29 | +#define ITERATIONS 10000000 |
| 30 | +#define REPEATS 3 |
| 31 | +#define MIN 1 |
| 32 | +#define MAX 6 |
| 33 | + |
| 34 | +/** |
| 35 | +* Prints the TAP version. |
| 36 | +*/ |
| 37 | +static void print_version( void ) { |
| 38 | + printf( "TAP version 13\n" ); |
| 39 | +} |
| 40 | + |
| 41 | +/** |
| 42 | +* Prints the TAP summary. |
| 43 | +* |
| 44 | +* @param total total number of tests |
| 45 | +* @param passing total number of passing tests |
| 46 | +*/ |
| 47 | +static void print_summary( int total, int passing ) { |
| 48 | + printf( "#\n" ); |
| 49 | + printf( "1..%d\n", total ); // TAP plan |
| 50 | + printf( "# total %d\n", total ); |
| 51 | + printf( "# pass %d\n", passing ); |
| 52 | + printf( "#\n" ); |
| 53 | + printf( "# ok\n" ); |
| 54 | +} |
| 55 | + |
| 56 | +/** |
| 57 | +* Prints benchmarks results. |
| 58 | +* |
| 59 | +* @param iterations number of iterations |
| 60 | +* @param elapsed elapsed time in seconds |
| 61 | +*/ |
| 62 | +static void print_results( int iterations, double elapsed ) { |
| 63 | + double rate = (double)iterations / elapsed; |
| 64 | + printf( " ---\n" ); |
| 65 | + printf( " iterations: %d\n", iterations ); |
| 66 | + printf( " elapsed: %0.9f\n", elapsed ); |
| 67 | + printf( " rate: %0.9f\n", rate ); |
| 68 | + printf( " ...\n" ); |
| 69 | +} |
| 70 | + |
| 71 | +/** |
| 72 | +* Returns a clock time. |
| 73 | +* |
| 74 | +* @return clock time |
| 75 | +*/ |
| 76 | +static double tic( void ) { |
| 77 | + struct timeval now; |
| 78 | + gettimeofday( &now, NULL ); |
| 79 | + return (double)now.tv_sec + (double)now.tv_usec/1.0e6; |
| 80 | +} |
| 81 | + |
| 82 | +/** |
| 83 | +* Generates a random number on the interval [0,1). |
| 84 | +* |
| 85 | +* @return random number |
| 86 | +*/ |
| 87 | +static float rand_float( void ) { |
| 88 | + int r = rand(); |
| 89 | + return (float)r / ( (float)RAND_MAX + 1.0f ); |
| 90 | +} |
| 91 | + |
| 92 | +/** |
| 93 | +* Runs a benchmark. |
| 94 | +* |
| 95 | +* @param iterations number of iterations |
| 96 | +* @param N array dimension size |
| 97 | +* @return elapsed time in seconds |
| 98 | +*/ |
| 99 | +static double benchmark1( int iterations, int N ) { |
| 100 | + stdlib_complex64_t alpha; |
| 101 | + stdlib_complex64_t beta; |
| 102 | + double elapsed; |
| 103 | + float *A; |
| 104 | + float *x; |
| 105 | + float *y; |
| 106 | + double t; |
| 107 | + int i; |
| 108 | + int j; |
| 109 | + |
| 110 | + A = (float *)malloc( N * N * 2 * sizeof( float ) ); |
| 111 | + x = (float *)malloc( N * 2 * sizeof( float ) ); |
| 112 | + y = (float *)malloc( N * 2 * sizeof( float ) ); |
| 113 | + alpha = stdlib_complex64( 0.5f, 0.5f ); |
| 114 | + beta = stdlib_complex64( 0.5f, -0.5f ); |
| 115 | + |
| 116 | + for ( i = 0, j = 0; i < N; i++, j += 2 ) { |
| 117 | + x[ i ] = ( rand_float()*2.0f ) - 1.0f; |
| 118 | + x[ i+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 119 | + y[ i ] = ( rand_float()*2.0f ) - 1.0f; |
| 120 | + y[ i+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 121 | + A[ j ] = ( rand_float()*2.0f ) - 1.0f; |
| 122 | + A[ j+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 123 | + } |
| 124 | + t = tic(); |
| 125 | + for ( i = 0; i < iterations; i++ ) { |
| 126 | + // cppcheck-suppress uninitvar |
| 127 | + c_cgemv( CblasRowMajor, CblasNoTrans, N, N, alpha, (void *)A, N, (void *)x, 1, beta, (void *)y, 1 ); |
| 128 | + if ( y[ i%N ] != y[ i%N ] ) { |
| 129 | + printf( "should not return NaN\n" ); |
| 130 | + break; |
| 131 | + } |
| 132 | + } |
| 133 | + elapsed = tic() - t; |
| 134 | + if ( y[ i%N ] != y[ i%N ] ) { |
| 135 | + printf( "should not return NaN\n" ); |
| 136 | + } |
| 137 | + free( A ); |
| 138 | + free( x ); |
| 139 | + free( y ); |
| 140 | + return elapsed; |
| 141 | +} |
| 142 | + |
| 143 | +/** |
| 144 | +* Runs a benchmark. |
| 145 | +* |
| 146 | +* @param iterations number of iterations |
| 147 | +* @param N array dimension size |
| 148 | +* @return elapsed time in seconds |
| 149 | +*/ |
| 150 | +static double benchmark2( int iterations, int N ) { |
| 151 | + stdlib_complex64_t alpha; |
| 152 | + stdlib_complex64_t beta; |
| 153 | + double elapsed; |
| 154 | + float *A; |
| 155 | + float *x; |
| 156 | + float *y; |
| 157 | + double t; |
| 158 | + int i; |
| 159 | + int j; |
| 160 | + |
| 161 | + A = (float *)malloc( N * N * 2 * sizeof( float ) ); |
| 162 | + x = (float *)malloc( N * 2 * sizeof( float ) ); |
| 163 | + y = (float *)malloc( N * 2 * sizeof( float ) ); |
| 164 | + alpha = stdlib_complex64( 0.5f, 0.5f ); |
| 165 | + beta = stdlib_complex64( 0.5f, -0.5f ); |
| 166 | + |
| 167 | + for ( i = 0, j = 0; i < N; i++, j += 2 ) { |
| 168 | + x[ i ] = ( rand_float()*2.0f ) - 1.0f; |
| 169 | + x[ i+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 170 | + y[ i ] = ( rand_float()*2.0f ) - 1.0f; |
| 171 | + y[ i+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 172 | + A[ j ] = ( rand_float()*2.0f ) - 1.0f; |
| 173 | + A[ j+1 ] = ( rand_float()*2.0f ) - 1.0f; |
| 174 | + } |
| 175 | + t = tic(); |
| 176 | + for ( i = 0; i < iterations; i++ ) { |
| 177 | + // cppcheck-suppress uninitvar |
| 178 | + c_cgemv_ndarray( CblasNoTrans, N, N, alpha, (void *)A, N, 1, 0, (void *)x, 1, 0, beta, (void *)y, 1, 0 ); |
| 179 | + if ( y[ i%N ] != y[ i%N ] ) { |
| 180 | + printf( "should not return NaN\n" ); |
| 181 | + break; |
| 182 | + } |
| 183 | + } |
| 184 | + elapsed = tic() - t; |
| 185 | + if ( y[ i%N ] != y[ i%N ] ){ |
| 186 | + printf( "should not return NaN\n" ); |
| 187 | + } |
| 188 | + free( A ); |
| 189 | + free( x ); |
| 190 | + free( y ); |
| 191 | + return elapsed; |
| 192 | +} |
| 193 | + |
| 194 | +/** |
| 195 | +* Main execution sequence. |
| 196 | +*/ |
| 197 | +int main( void ) { |
| 198 | + double elapsed; |
| 199 | + int count; |
| 200 | + int iter; |
| 201 | + int N; |
| 202 | + int i; |
| 203 | + int j; |
| 204 | + |
| 205 | + // Use the current time to seed the random number generator: |
| 206 | + srand( time( NULL ) ); |
| 207 | + |
| 208 | + print_version(); |
| 209 | + count = 0; |
| 210 | + for ( i = MIN; i <= MAX; i++ ) { |
| 211 | + N = floor( pow( pow( 10, i ), 1.0/2.0 ) ); |
| 212 | + iter = ITERATIONS / pow( 10, i-1 ); |
| 213 | + for ( j = 0; j < REPEATS; j++ ) { |
| 214 | + count += 1; |
| 215 | + printf( "# c::%s:size=%d\n", NAME, N*N ); |
| 216 | + elapsed = benchmark1( iter, N ); |
| 217 | + print_results( iter, elapsed ); |
| 218 | + printf( "ok %d benchmark finished\n", count ); |
| 219 | + } |
| 220 | + for ( j = 0; j < REPEATS; j++ ) { |
| 221 | + count += 1; |
| 222 | + printf( "# c::%s:ndarray:size=%d\n", NAME, N*N ); |
| 223 | + elapsed = benchmark2( iter, N ); |
| 224 | + print_results( iter, elapsed ); |
| 225 | + printf( "ok %d benchmark finished\n", count ); |
| 226 | + } |
| 227 | + } |
| 228 | + print_summary( count, count ); |
| 229 | +} |
0 commit comments