Skip to content

feat: add C implementation to blas/base/ssyr #7127

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 12 commits into from
Jul 17, 2025
Merged

feat: add C implementation to blas/base/ssyr #7127

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
8e6c1fd
refactor: c implementation of blas/base/ssyr
ShabiShett07 May 30, 2025
f5414a1
chore: update packages
ShabiShett07 May 30, 2025
530e616
chore: minor clean-up
ShabiShett07 May 30, 2025
0cb1c6b
chore: minor clean-up
ShabiShett07 May 30, 2025
bae7044
chore: minor clean-up
ShabiShett07 May 30, 2025
2df183e
chore: minor clean-up
ShabiShett07 May 30, 2025
4bbde86
Merge branch 'develop' into feature/c/ssyr
ShabiShett07 Jul 11, 2025
1117111
refactor: replicate dsyr changes
ShabiShett07 Jul 12, 2025
3ed76f5
fix: change dsyr to ssyr
ShabiShett07 Jul 12, 2025
ca04fbe
fix: change parameter naming
ShabiShett07 Jul 12, 2025
62abab2
refactor: c implementation
ShabiShett07 Jul 12, 2025
8d9d9a8
chore: clean-up
kgryte Jul 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
121 changes: 70 additions & 51 deletions lib/node_modules/@stdlib/blas/base/ssyr/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@

@license Apache-2.0

Copyright (c) 2025 The Stdlib Authors.
Copyright (c) 2024 The Stdlib Authors.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
Expand All @@ -20,7 +20,7 @@ limitations under the License.

# ssyr

> Perform the symmetric rank 1 operation `A = α*x*x**T + A`.
> Perform the symmetric rank 1 operation `A = α*x*x^T + A`.

<section class="usage">

Expand All @@ -32,16 +32,16 @@ var ssyr = require( '@stdlib/blas/base/ssyr' );

#### ssyr( order, uplo, N, α, x, sx, A, LDA )

Performs the symmetric rank 1 operation `A = α*x*x**T + A` where `α` is a scalar, `x` is an `N` element vector, and `A` is an `N` by `N` symmetric matrix.
Performs the symmetric rank 1 operation `A = α*x*x^T + A` where `α` is a scalar, `x` is an `N` element vector, and `A` is an `N` by `N` symmetric matrix.

```javascript
var Float32Array = require( '@stdlib/array/float32' );

var A = new Float32Array( [ 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 0.0, 0.0, 1.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 ] );
var x = new Float32Array( [ 1.0, 2.0, 3.0 ] );

ssyr( 'row-major', 'upper', 3, 1.0, x, 1, A, 3 );
// A => <Float32Array>[ 2.0, 4.0, 6.0, 0.0, 5.0, 8.0, 0.0, 0.0, 10.0 ]
// A => <Float32Array>[ 2.0, 4.0, 6.0, 2.0, 5.0, 8.0, 3.0, 2.0, 10.0 ]
```

The function has the following parameters:
Expand All @@ -51,20 +51,20 @@ The function has the following parameters:
- **N**: number of elements along each dimension of `A`.
- **α**: scalar constant.
- **x**: input [`Float32Array`][mdn-float32array].
- **sx**: index increment for `x`.
- **sx**: stride length for `x`.
- **A**: input matrix stored in linear memory as a [`Float32Array`][mdn-float32array].
- **lda**: stride of the first dimension of `A` (a.k.a., leading dimension of the matrix `A`).
- **LDA**: stride of the first dimension of `A` (a.k.a., leading dimension of the matrix `A`).

The stride parameters determine how elements in the input arrays are accessed at runtime. For example, to iterate over every other element of `x` in reverse order,
The stride parameters determine how elements in the input arrays are accessed at runtime. For example, to iterate over the elements of `x` in reverse order,

```javascript
var Float32Array = require( '@stdlib/array/float32' );

var A = new Float32Array( [ 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 0.0, 0.0, 1.0 ] );
var x = new Float32Array( [ 1.0, 2.0, 3.0, 4.0, 5.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 ] );
var x = new Float32Array( [ 3.0, 2.0, 1.0 ] );

ssyr( 'row-major', 'upper', 3, 1.0, x, -2, A, 3 );
// A => <Float32Array>[ 26.0, 17.0, 8.0, 0.0, 10.0, 5.0, 0.0, 0.0, 2.0 ]
ssyr( 'row-major', 'upper', 3, 1.0, x, -1, A, 3 );
// A => <Float32Array>[ 2.0, 4.0, 6.0, 2.0, 5.0, 8.0, 3.0, 2.0, 10.0 ]
```

Note that indexing is relative to the first index. To introduce an offset, use [`typed array`][mdn-typed-array] views.
Expand All @@ -75,14 +75,14 @@ Note that indexing is relative to the first index. To introduce an offset, use [
var Float32Array = require( '@stdlib/array/float32' );

// Initial arrays...
var x0 = new Float32Array( [ 1.0, 1.0, 1.0, 1.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 0.0, 0.0, 1.0 ] );
var x0 = new Float32Array( [ 0.0, 3.0, 2.0, 1.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 ] );

// Create offset views...
var x1 = new Float32Array( x0.buffer, x0.BYTES_PER_ELEMENT*1 ); // start at 2nd element

ssyr( 'row-major', 'upper', 3, 1.0, x1, -1, A, 3 );
// A => <Float32Array>[ 2.0, 3.0, 4.0, 0.0, 2.0, 3.0, 0.0, 0.0, 2.0 ]
// A => <Float32Array>[ 2.0, 4.0, 6.0, 2.0, 5.0, 8.0, 3.0, 2.0, 10.0 ]
```

#### ssyr.ndarray( uplo, N, α, x, sx, ox, A, sa1, sa2, oa )
Expand All @@ -92,11 +92,11 @@ Performs the symmetric rank 1 operation `A = α*x*x^T + A`, using alternative in
```javascript
var Float32Array = require( '@stdlib/array/float32' );

var A = new Float32Array( [ 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 0.0, 0.0, 1.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 ] );
var x = new Float32Array( [ 1.0, 2.0, 3.0 ] );

ssyr.ndarray( 'upper', 3, 1.0, x, 1, 0, A, 3, 1, 0 );
// A => <Float32Array>[ 2.0, 4.0, 6.0, 0.0, 5.0, 8.0, 0.0, 0.0, 10.0 ]
// A => <Float32Array>[ 2.0, 4.0, 6.0, 2.0, 5.0, 8.0, 3.0, 2.0, 10.0 ]
```

The function has the following additional parameters:
Expand All @@ -111,11 +111,11 @@ While [`typed array`][mdn-typed-array] views mandate a view offset based on the
```javascript
var Float32Array = require( '@stdlib/array/float32' );

var A = new Float32Array( [ 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 0.0, 0.0, 1.0 ] );
var A = new Float32Array( [ 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 ] );
var x = new Float32Array( [ 1.0, 2.0, 3.0, 4.0, 5.0 ] );

ssyr.ndarray( 'upper', 3, 1.0, x, -2, 4, A, 3, 1, 0 );
// A => <Float32Array>[ 26.0, 17.0, 8.0, 0.0, 10.0, 5.0, 0.0, 0.0, 2.0 ]
// A => <Float32Array>[ 26.0, 17.0, 8.0, 2.0, 10.0, 5.0, 3.0, 2.0, 2.0 ]
```

</section>
Expand Down Expand Up @@ -149,11 +149,18 @@ var opts = {

var N = 3;

var A = ones( N*N, opts.dtype );
// Create N-by-N symmetric matrices:
var A1 = ones( N*N, opts.dtype );
var A2 = ones( N*N, opts.dtype );

// Create a random vector:
var x = discreteUniform( N, -10.0, 10.0, opts );

ssyr( 'row-major', 'upper', 3, 1.0, x, 1, A, 3 );
console.log( A );
ssyr( 'row-major', 'upper', 3, 1.0, x, 1, A1, 3 );
console.log( A1 );

ssyr.ndarray( 'upper', 3, 1.0, x, 1, 0, A2, 3, 1, 0 );
console.log( A2 );
```

</section>
Expand Down Expand Up @@ -186,62 +193,62 @@ console.log( A );
#include "stdlib/blas/base/ssyr.h"
```

#### c_ssyr( order, uplo, N, alpha, \*X, strideX, \*A, LDA )
#### c_ssyr( layout, uplo, N, alpha, \*X, sx, \*A, LDA )

Performs the symmetric rank 1 operation `A = α*x*x^T + A`.
Performs the symmetric rank 1 operation `A = α*x*x^T + A` where `α` is a scalar, `x` is an `N` element vector, and `A` is an `N` by `N` symmetric matrix.

```c
#include "stdlib/blas/base/shared.h"

float A[] = { 1.0f, 0.0f, 0.0f, 2.0f, 1.0f, 0.0f, 3.0f, 2.0f, 1.0f };
const float x[] = { 1.0f, 2.0f, 3.0f };
float A[] = { 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 };
const float x[] = { 1.0, 2.0, 3.0 };

c_ssyr( CblasColMajor, CblasUpper, 3, 1.0f, x, 1, A, 3 );
c_ssyr( CblasColMajor, CblasUpper, 3, 1.0, x, 1, A, 3 );
```

The function accepts the following arguments:

- **order**: `[in] CBLAS_LAYOUT` storage layout.
- **layout**: `[in] CBLAS_LAYOUT` storage layout.
- **uplo**: `[in] CBLAS_UPLO` specifies whether the upper or lower triangular part of the symmetric matrix `A` should be referenced.
- **N**: `[in] CBLAS_INT` number of elements along each dimension of `A`.
- **alpha**: `[in] float` scalar.
- **alpha**: `[in] float` scalar constant.
- **X**: `[in] float*` input array.
- **strideX**: `[in] CBLAS_INT` index increment for `X`.
- **sx**: `[in] CBLAS_INT` stride length for `X`.
- **A**: `[inout] float*` input matrix.
- **LDA**: `[in] CBLAS_INT` stride of the first dimension of `A` (a.k.a., leading dimension of the matrix `A`).

```c
void c_ssyr( const CBLAS_LAYOUT order, const CBLAS_UPLO uplo, const CBLAS_INT N, const float alpha, const float *X, const CBLAS_INT strideX, float *A, const CBLAS_INT LDA )
void c_ssyr( const CBLAS_LAYOUT layout, const CBLAS_UPLO uplo, const CBLAS_INT N, const float alpha, const float *X, const CBLAS_INT strideX, float *A, const CBLAS_INT LDA )
```

#### c_ssyr_ndarray( uplo, N, alpha, \*X, strideX, offsetX, \*A, sa1, sa2, oa )
#### c_ssyr_ndarray( uplo, N, alpha, \*X, sx, ox, \*A, sa1, sa2, oa )

Performs the symmetric rank 1 operation `A = α*x*x^T + A` using alternative indexing semantics.
Performs the symmetric rank 1 operation `A = α*x*x^T + A`, using alternative indexing semantics and where `α` is a scalar, `x` is an `N` element vector, and `A` is an `N` by `N` symmetric matrix.

```c
#include "stdlib/blas/base/shared.h"

float A[] = { 1.0f, 2.0f, 3.0f, 0.0f, 1.0f, 2.0f, 0.0f, 0.0f, 1.0f };
const float x[] = { 1.0f, 2.0f, 3.0f };
float A[] = { 1.0, 2.0, 3.0, 2.0, 1.0, 2.0, 3.0, 2.0, 1.0 };
const float x[] = { 1.0, 2.0, 3.0 };

c_ssyr_ndarray( CblasUpper, 3, 1.0f, x, 1, 0, A, 3, 1, 0 );
c_ssyr_ndarray( CblasUpper, 3, 1.0, x, 1, 0, A, 3, 1, 0 );
```

The function accepts the following arguments:

- **uplo**: `[in] CBLAS_UPLO` specifies whether the upper or lower triangular part of the symmetric matrix `A` should be referenced.
- **N**: `[in] CBLAS_INT` number of elements along each dimension of `A`.
- **alpha**: `[in] float` scalar.
- **alpha**: `[in] float` scalar constant.
- **X**: `[in] float*` input array.
- **strideX**: `[in] CBLAS_INT` index increment for `X`.
- **offsetX**: `[in] CBLAS_INT` starting index for `X`.
- **sx**: `[in] CBLAS_INT` stride length for `X`.
- **ox**: `[in] CBLAS_INT` starting index for `X`.
- **A**: `[inout] float*` input matrix.
- **sa1**: `[in] CBLAS_INT` stride of the first dimension of `A`.
- **sa2**: `[in] CBLAS_INT` stride of the second dimension of `A`.
- **oa**: `[in] CBLAS_INT` starting index for `A`.

```c
void c_ssyr_ndarray( const CBLAS_UPLO uplo, const CBLAS_INT N, const float alpha, const float *X, const CBLAS_INT strideX, const CBLAS_INT offsetX, float *A, const CBLAS_INT sa1, const CBLAS_INT sa2, const CBLAS_INT oa )
void c_ssyr_ndarray( const CBLAS_UPLO uplo, const CBLAS_INT N, const float alpha, const float *X, const CBLAS_INT strideX, const CBLAS_INT offsetX, float *A, const CBLAS_INT strideA1, const CBLAS_INT strideA2, const CBLAS_INT offsetA )
```

</section>
Expand All @@ -268,27 +275,39 @@ void c_ssyr_ndarray( const CBLAS_UPLO uplo, const CBLAS_INT N, const float alpha
#include <stdio.h>

int main( void ) {
// Create strided arrays:
float A[] = { 1.0f, 0.0f, 0.0f, 2.0f, 1.0f, 0.0f, 3.0f, 2.0f, 1.0f };
const float x[] = { 1.0f, 2.0f, 3.0f };

// Specify the number of elements along each dimension of `A`:
// Define 3x3 symmetric matrices stored in row-major layout:
float A1[ 3*3 ] = {
1.0f, 2.0f, 3.0f,
2.0f, 1.0f, 2.0f,
3.0f, 2.0f, 1.0f
};

float A2[ 3*3 ] = {
1.0f, 2.0f, 3.0f,
2.0f, 1.0f, 2.0f,
3.0f, 2.0f, 1.0f
};

// Define a vector:
const float x[ 3 ] = { 1.0f, 2.0f, 3.0f };

// Specify the number of elements along each dimension of `A1` and `A2`:
const int N = 3;

// Perform the symmetric rank 1 operation `A = α*x*x^T + A`:
c_ssyr( CblasColMajor, CblasUpper, N, 1.0f, x, 1, A, N );
c_ssyr( CblasColMajor, CblasUpper, N, 1.0f, x, 1, A1, N );

// Print the result:
for ( int i = 0; i < N*N; i++ ) {
printf( "A[ %i ] = %f\n", i, A[ i ] );
printf( "A1[ %i ] = %f\n", i, A1[ i ] );
}

// Perform the symmetric rank 1 operation `A = α*x*x^T + A`:
c_ssyr_ndarray( CblasUpper, N, 1.0f, x, 1, 0, A, N, 1, 0 );
// Perform the symmetric rank 1 operation `A = α*x*x^T + A` using alternative indexing semantics:
c_ssyr_ndarray( CblasUpper, N, 1.0, x, 1, 0, A2, N, 1, 0 );

// Print the result:
for ( int i = 0; i < N*N; i++ ) {
printf( "A[ %i ] = %f\n", i, A[ i ] );
printf( "A2[ %i ] = %f\n", i, A[ i ] );
}
}
```
Expand All @@ -315,7 +334,7 @@ int main( void ) {

[blas]: http://www.netlib.org/blas

[blas-ssyr]: https://www.netlib.org/lapack/explore-html/dc/d82/group__her_gad7585662770cdd3001ed08c7a864cd21.html#gad7585662770cdd3001ed08c7a864cd21
[blas-ssyr]: https://www.netlib.org/lapack/explore-html-3.6.1/d6/d30/group__single__blas__level2_ga7b8a99048765ed2bf7c1e770bff0b622.html

[mdn-float32array]: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Float32Array

Expand Down
20 changes: 10 additions & 10 deletions lib/node_modules/@stdlib/blas/base/ssyr/benchmark/c/benchmark.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ static double tic( void ) {
* Runs a benchmark.
*
* @param iterations number of iterations
* @param N number of elements along each dimension
* @param N array dimension size
* @return elapsed time in seconds
*/
static double benchmark1( int iterations, int N ) {
Expand All @@ -93,18 +93,18 @@ static double benchmark1( int iterations, int N ) {
double t;
int i;

stdlib_strided_sfill( N, 0.5f, x, 1 );
stdlib_strided_sfill( N, 1.0f, x, 1 );
stdlib_strided_sfill( N*N, 1.0f, A, 1 );
t = tic();
for ( i = 0; i < iterations; i++ ) {
c_ssyr( CblasRowMajor, CblasUpper, N, 1.0, x, 1, A, N );
if ( A[ 0 ] != A[ 0 ] ) {
c_ssyr( CblasRowMajor, CblasUpper, N, 1.0f, x, 1, A, N );
if ( A[ i%(N*2) ] != A[ i%(N*2) ] ) {
printf( "should not return NaN\n" );
break;
}
}
elapsed = tic() - t;
if ( A[ 0 ] != A[ 0 ] ) {
if ( A[ i%(N*2) ] != A[ i%(N*2) ] ) {
printf( "should not return NaN\n" );
}
return elapsed;
Expand All @@ -114,7 +114,7 @@ static double benchmark1( int iterations, int N ) {
* Runs a benchmark.
*
* @param iterations number of iterations
* @param N number of elements along each dimension
* @param N array dimension size
* @return elapsed time in seconds
*/
static double benchmark2( int iterations, int N ) {
Expand All @@ -124,18 +124,18 @@ static double benchmark2( int iterations, int N ) {
double t;
int i;

stdlib_strided_sfill( N, 0.5f, x, 1 );
stdlib_strided_sfill( N, 1.0f, x, 1 );
stdlib_strided_sfill( N*N, 1.0f, A, 1 );
t = tic();
for ( i = 0; i < iterations; i++ ) {
c_ssyr_ndarray( CblasUpper, N, 1.0, x, 1, 0, A, N, 1, 0 );
if ( A[ 0 ] != A[ 0 ] ) {
c_ssyr_ndarray( CblasUpper, N, 1.0f, x, 1, 0, A, N, 1, 0 );
if ( A[ i%(N*2) ] != A[ i%(N*2) ] ) {
printf( "should not return NaN\n" );
break;
}
}
elapsed = tic() - t;
if ( A[ 0 ] != A[ 0 ] ) {
if ( A[ i%(N*2) ] != A[ i%(N*2) ] ) {
printf( "should not return NaN\n" );
}
return elapsed;
Expand Down
22 changes: 18 additions & 4 deletions lib/node_modules/@stdlib/blas/base/ssyr/docs/repl.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,10 +46,24 @@

Examples
--------
// Standard usage:
> var x = new {{alias:@stdlib/array/float32}}( [ 1.0, 1.0 ] );
> var A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 0.0, 2.0 ] );
> var A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 2.0, 1.0 ] );
> {{alias}}( 'row-major', 'upper', 2, 1.0, x, 1, A, 2 )
<Float32Array>[ 2.0, 3.0, 0.0, 3.0 ]
<Float32Array>[ 2.0, 3.0, 2.0, 2.0 ]

// Advanced indexing:
> var x = new {{alias:@stdlib/array/float32}}( [ 1.0, 1.0 ] );
> var A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 2.0, 1.0 ] );
> {{alias}}( 'row-major', 'upper', 2, 1.0, x, -1, A, 2 )
<Float32Array>[ 2.0, 3.0, 2.0, 2.0 ]

// Using typed array views:
> var x0 = new {{alias:@stdlib/array/float32}}( [ 0.0, 1.0, 1.0 ] );
> A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 2.0, 1.0 ] );
> var x1 = new {{alias:@stdlib/array/float32}}( x0.buffer, x0.BYTES_PER_ELEMENT*1 );
> {{alias}}( 'row-major', 'upper', 2, 1.0, x, 1, A, 2 )
<Float32Array>[ 2.0, 3.0, 2.0, 2.0 ]


{{alias}}.ndarray( uplo, N, α, x, sx, ox, A, sa1, sa2, oa )
Expand Down Expand Up @@ -102,9 +116,9 @@
Examples
--------
> var x = new {{alias:@stdlib/array/float32}}( [ 1.0, 1.0 ] );
> var A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 0.0, 2.0 ] );
> var A = new {{alias:@stdlib/array/float32}}( [ 1.0, 2.0, 2.0, 1.0 ] );
> {{alias}}.ndarray( 'upper', 2, 1.0, x, 1, 0, A, 2, 1, 0 )
<Float32Array>[ 2.0, 3.0, 0.0, 3.0 ]
<Float32Array>[ 2.0, 3.0, 2.0, 2.0 ]

See Also
--------
Loading