Skip to content

Commit 29767e5

Browse files
committed
feat: add c implementation for blas/base/dgemm
--- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: passed - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: passed - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: missing_dependencies - task: lint_c_examples status: missing_dependencies - task: lint_c_benchmarks status: missing_dependencies - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent e09a028 commit 29767e5

22 files changed

Lines changed: 5154 additions & 6 deletions

lib/node_modules/@stdlib/blas/base/dgemm/README.md

Lines changed: 112 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -190,18 +190,79 @@ console.log( C );
190190
#include "stdlib/blas/base/dgemm.h"
191191
```
192192

193-
#### TODO
193+
#### c_dgemm( layout, transA, transB, M, N, K, alpha, \*A, LDA, \*B, LDB, beta, \*C, LDC )
194194

195-
TODO.
195+
Performs the matrix-matrix operation `C = α*op(A)*op(B) + β*C` where `op(X)` is either `op(X) = X` or `op(X) = X^T`, `α` and `β` are scalars, `A`, `B`, and `C` are matrices, with `op(A)` an `M` by `K` matrix, `op(B)` a `K` by `N` matrix, and `C` an `M` by `N` matrix.
196+
197+
```c
198+
#include "stdlib/blas/base/shared.h"
199+
200+
const double A[] = { 1.0, 2.0, 3.0, 4.0 };
201+
const double B[] = { 1.0, 1.0, 0.0, 1.0 };
202+
double C[] = { 1.0, 2.0, 3.0, 4.0 };
203+
204+
c_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, 2, 2, 2, 1.0, A, 2, B, 2, 1.0, C, 2 );
205+
```
206+
207+
The function accepts the following arguments:
208+
209+
- **layout**: `[in] CBLAS_LAYOUT` storage layout.
210+
- **transA**: `[in] CBLAS_TRANSPOSE` specifies whether `A` should be transposed, conjugate-transposed, or not transposed.
211+
- **transB**: `[in] CBLAS_TRANSPOSE` specifies whether `B` should be transposed, conjugate-transposed, or not transposed.
212+
- **M**: `[in] CBLAS_INT` number of rows in the matrix `op(A)` and in the matrix `C`.
213+
- **N**: `[in] CBLAS_INT` number of columns in the matrix `op(B)` and in the matrix `C`.
214+
- **K**: `[in] CBLAS_INT` number of columns in the matrix `op(A)` and number of rows in the matrix `op(B)`.
215+
- **alpha**: `[in] double` scalar constant.
216+
- **A**: `[in] double*` first input matrix.
217+
- **LDA**: `[in] CBLAS_INT` stride of the first dimension of `A` (a.k.a., leading dimension of the matrix `A`).
218+
- **B**: `[in] double*` second input matrix.
219+
- **LDB**: `[in] CBLAS_INT` stride of the first dimension of `B` (a.k.a., leading dimension of the matrix `B`).
220+
- **beta**: `[in] double` scalar constant.
221+
- **C**: `[inout] double*` third input matrix.
222+
- **LDC**: `[in] CBLAS_INT` stride of the first dimension of `C` (a.k.a., leading dimension of the matrix `C`).
223+
224+
```c
225+
void c_dgemm( const CBLAS_LAYOUT layout, const CBLAS_TRANSPOSE transA, const CBLAS_TRANSPOSE transB, const CBLAS_INT M, const CBLAS_INT N, const CBLAS_INT K, const double alpha, const double *A, const CBLAS_INT LDA, const double *B, const CBLAS_INT LDB, const double beta, double *C, const CBLAS_INT LDC )
226+
```
227+
228+
#### c_dgemm_ndarray( transA, transB, M, N, K, alpha, \*A, sa1, sa2, oa, \*B, sb1, sb2, ob, beta, \*C, sc1, sc2, oc )
229+
230+
Performs the matrix-matrix operation `C = α*op(A)*op(B) + β*C`, using alternative indexing semantics and where `op(X)` is either `op(X) = X` or `op(X) = X^T`, `α` and `β` are scalars, `A`, `B`, and `C` are matrices, with `op(A)` an `M` by `K` matrix, `op(B)` a `K` by `N` matrix, and `C` an `M` by `N` matrix.
196231

197232
```c
198-
TODO
233+
#include "stdlib/blas/base/shared.h"
234+
235+
const double A[] = { 1.0, 2.0, 3.0, 4.0 };
236+
const double B[] = { 1.0, 1.0, 0.0, 1.0 };
237+
double C[] = { 1.0, 2.0, 3.0, 4.0 };
238+
239+
c_dgemm_ndarray( CblasNoTrans, CblasNoTrans, 2, 2, 2, 1.0, A, 2, 1, 0, B, 2, 1, 0, 1.0, C, 2, 1, 0 );
199240
```
200241
201-
TODO
242+
The function accepts the following arguments:
243+
244+
- **transA**: `[in] CBLAS_TRANSPOSE` specifies whether `A` should be transposed, conjugate-transposed, or not transposed.
245+
- **transB**: `[in] CBLAS_TRANSPOSE` specifies whether `B` should be transposed, conjugate-transposed, or not transposed.
246+
- **M**: `[in] CBLAS_INT` number of rows in the matrix `op(A)` and in the matrix `C`.
247+
- **N**: `[in] CBLAS_INT` number of columns in the matrix `op(B)` and in the matrix `C`.
248+
- **K**: `[in] CBLAS_INT` number of columns in the matrix `op(A)` and number of rows in the matrix `op(B)`.
249+
- **alpha**: `[in] double` scalar constant.
250+
- **A**: `[in] double*` first input matrix.
251+
- **sa1**: `[in] CBLAS_INT` stride of the first dimension of `A`.
252+
- **sa2**: `[in] CBLAS_INT` stride of the second dimension of `A`.
253+
- **oa**: `[in] CBLAS_INT` starting index for `A`.
254+
- **B**: `[in] double*` second input matrix.
255+
- **sb1**: `[in] CBLAS_INT` stride of the first dimension of `B`.
256+
- **sb2**: `[in] CBLAS_INT` stride of the second dimension of `B`.
257+
- **ob**: `[in] CBLAS_INT` starting index for `B`.
258+
- **beta**: `[in] double` scalar constant.
259+
- **C**: `[inout] double*` third input matrix.
260+
- **sc1**: `[in] CBLAS_INT` stride of the first dimension of `C`.
261+
- **sc2**: `[in] CBLAS_INT` stride of the second dimension of `C`.
262+
- **oc**: `[in] CBLAS_INT` starting index for `C`.
202263
203264
```c
204-
TODO
265+
void c_dgemm_ndarray( const CBLAS_TRANSPOSE transA, const CBLAS_TRANSPOSE transB, const CBLAS_INT M, const CBLAS_INT N, const CBLAS_INT K, const double alpha, const double *A, const CBLAS_INT strideA1, const CBLAS_INT strideA2, const CBLAS_INT offsetA, const double *B, const CBLAS_INT strideB1, const CBLAS_INT strideB2, const CBLAS_INT offsetB, const double beta, double *C, const CBLAS_INT strideC1, const CBLAS_INT strideC2, const CBLAS_INT offsetC )
205266
```
206267

207268
</section>
@@ -223,7 +284,52 @@ TODO
223284
### Examples
224285

225286
```c
226-
TODO
287+
#include "stdlib/blas/base/dgemm.h"
288+
#include "stdlib/blas/base/shared.h"
289+
#include <stdio.h>
290+
291+
int main( void ) {
292+
// Define matrices stored in row-major order:
293+
const double A[ 3*2 ] = {
294+
1.0, 2.0,
295+
3.0, 4.0,
296+
5.0, 6.0
297+
};
298+
const double B[ 2*4 ] = {
299+
1.0, 2.0, 3.0, 4.0,
300+
5.0, 6.0, 7.0, 8.0
301+
};
302+
double C[ 3*4 ] = {
303+
1.0, 2.0, 3.0, 4.0,
304+
5.0, 6.0, 7.0, 8.0,
305+
9.0, 10.0, 11.0, 12.0
306+
};
307+
308+
// Specify matrix dimensions:
309+
const int M = 3;
310+
const int N = 4;
311+
const int K = 2;
312+
313+
// Perform the matrix-matrix operation `C = α*A*B + β*C`:
314+
c_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans, M, N, K, 1.0, A, K, B, N, 1.0, C, N );
315+
316+
// Print the result:
317+
for ( int i = 0; i < M; i++ ) {
318+
for ( int j = 0; j < N; j++ ) {
319+
printf( "C[ %i ][ %i ] = %lf\n", i, j, C[ (i*N)+j ] );
320+
}
321+
}
322+
323+
// Perform the matrix-matrix operation `C = α*A*B + β*C` using alternative indexing semantics:
324+
c_dgemm_ndarray( CblasNoTrans, CblasNoTrans, M, N, K, 1.0, A, K, 1, 0, B, N, 1, 0, 1.0, C, N, 1, 0 );
325+
326+
// Print the result:
327+
for ( int i = 0; i < M; i++ ) {
328+
for ( int j = 0; j < N; j++ ) {
329+
printf( "C[ %i ][ %i ] = %lf\n", i, j, C[ (i*N)+j ] );
330+
}
331+
}
332+
}
227333
```
228334
229335
</section>
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
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+
'use strict';
20+
21+
// MODULES //
22+
23+
var resolve = require( 'path' ).resolve;
24+
var bench = require( '@stdlib/bench' );
25+
var uniform = require( '@stdlib/random/array/uniform' );
26+
var format = require( '@stdlib/string/format' );
27+
var isnan = require( '@stdlib/math/base/assert/is-nan' );
28+
var pow = require( '@stdlib/math/base/special/pow' );
29+
var floor = require( '@stdlib/math/base/special/floor' );
30+
var tryRequire = require( '@stdlib/utils/try-require' );
31+
var pkg = require( './../package.json' ).name;
32+
33+
34+
// VARIABLES //
35+
36+
var dgemm = tryRequire( resolve( __dirname, './../lib/dgemm.native.js' ) );
37+
var opts = {
38+
'skip': ( dgemm instanceof Error )
39+
};
40+
var options = {
41+
'dtype': 'float64'
42+
};
43+
44+
45+
// FUNCTIONS //
46+
47+
/**
48+
* Creates a benchmark function.
49+
*
50+
* @private
51+
* @param {PositiveInteger} N - array dimension size
52+
* @returns {Function} benchmark function
53+
*/
54+
function createBenchmark( N ) {
55+
var A = uniform( N*N, -10.0, 10.0, options );
56+
var B = uniform( N*N, -10.0, 10.0, options );
57+
var C = uniform( N*N, -10.0, 10.0, options );
58+
return benchmark;
59+
60+
/**
61+
* Benchmark function.
62+
*
63+
* @private
64+
* @param {Benchmark} b - benchmark instance
65+
*/
66+
function benchmark( b ) {
67+
var z;
68+
var i;
69+
70+
b.tic();
71+
for ( i = 0; i < b.iterations; i++ ) {
72+
z = dgemm( 'row-major', 'no-transpose', 'no-transpose', N, N, N, 1.0, A, N, B, N, 1.0, C, N );
73+
if ( isnan( z[ i%z.length ] ) ) {
74+
b.fail( 'should not return NaN' );
75+
}
76+
}
77+
b.toc();
78+
if ( isnan( z[ i%z.length ] ) ) {
79+
b.fail( 'should not return NaN' );
80+
}
81+
b.pass( 'benchmark finished' );
82+
b.end();
83+
}
84+
}
85+
86+
87+
// MAIN //
88+
89+
/**
90+
* Main execution sequence.
91+
*
92+
* @private
93+
*/
94+
function main() {
95+
var min;
96+
var max;
97+
var N;
98+
var f;
99+
var i;
100+
101+
min = 1; // 10^min
102+
max = 5; // 10^max
103+
104+
for ( i = min; i <= max; i++ ) {
105+
N = floor( pow( pow( 10, i ), 1.0/2.0 ) );
106+
f = createBenchmark( N );
107+
bench( format( '%s::native:size=%d', pkg, N*N ), opts, f );
108+
}
109+
}
110+
111+
main();
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
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+
'use strict';
20+
21+
// MODULES //
22+
23+
var resolve = require( 'path' ).resolve;
24+
var bench = require( '@stdlib/bench' );
25+
var uniform = require( '@stdlib/random/array/uniform' );
26+
var format = require( '@stdlib/string/format' );
27+
var isnan = require( '@stdlib/math/base/assert/is-nan' );
28+
var pow = require( '@stdlib/math/base/special/pow' );
29+
var floor = require( '@stdlib/math/base/special/floor' );
30+
var tryRequire = require( '@stdlib/utils/try-require' );
31+
var pkg = require( './../package.json' ).name;
32+
33+
34+
// VARIABLES //
35+
36+
var dgemm = tryRequire( resolve( __dirname, './../lib/ndarray.native.js' ) );
37+
var opts = {
38+
'skip': ( dgemm instanceof Error )
39+
};
40+
var options = {
41+
'dtype': 'float64'
42+
};
43+
44+
45+
// FUNCTIONS //
46+
47+
/**
48+
* Creates a benchmark function.
49+
*
50+
* @private
51+
* @param {PositiveInteger} N - array dimension size
52+
* @returns {Function} benchmark function
53+
*/
54+
function createBenchmark( N ) {
55+
var A = uniform( N*N, -10.0, 10.0, options );
56+
var B = uniform( N*N, -10.0, 10.0, options );
57+
var C = uniform( N*N, -10.0, 10.0, options );
58+
return benchmark;
59+
60+
/**
61+
* Benchmark function.
62+
*
63+
* @private
64+
* @param {Benchmark} b - benchmark instance
65+
*/
66+
function benchmark( b ) {
67+
var z;
68+
var i;
69+
70+
b.tic();
71+
for ( i = 0; i < b.iterations; i++ ) {
72+
z = dgemm( 'no-transpose', 'no-transpose', N, N, N, 1.0, A, 1, N, 0, B, 1, N, 0, 1.0, C, 1, N, 0 );
73+
if ( isnan( z[ i%z.length ] ) ) {
74+
b.fail( 'should not return NaN' );
75+
}
76+
}
77+
b.toc();
78+
if ( isnan( z[ i%z.length ] ) ) {
79+
b.fail( 'should not return NaN' );
80+
}
81+
b.pass( 'benchmark finished' );
82+
b.end();
83+
}
84+
}
85+
86+
87+
// MAIN //
88+
89+
/**
90+
* Main execution sequence.
91+
*
92+
* @private
93+
*/
94+
function main() {
95+
var min;
96+
var max;
97+
var N;
98+
var f;
99+
var i;
100+
101+
min = 1; // 10^min
102+
max = 5; // 10^max
103+
104+
for ( i = min; i <= max; i++ ) {
105+
N = floor( pow( pow( 10, i ), 1.0/2.0 ) );
106+
f = createBenchmark( N );
107+
bench( format( '%s::native:ndarray:size=%d', pkg, N*N ), opts, f );
108+
}
109+
}
110+
111+
main();

0 commit comments

Comments
 (0)