Skip to content

Commit 6c90e8f

Browse files
authored
feat: add accessor arrays support to blas/ext/base/gcusumpw
PR-URL: #5012 Reviewed-by: Athan Reines <[email protected]>
1 parent 1dd8adf commit 6c90e8f

File tree

7 files changed

+596
-41
lines changed

7 files changed

+596
-41
lines changed

Diff for: lib/node_modules/@stdlib/blas/ext/base/gcusumpw/README.md

+3
Original file line numberDiff line numberDiff line change
@@ -129,6 +129,7 @@ gcusumpw.ndarray( 4, 0.0, x, 2, 1, y, -1, y.length-1 );
129129

130130
- If `N <= 0`, both functions return `y` unchanged.
131131
- In general, pairwise summation is more numerically stable than ordinary recursive summation (i.e., "simple" summation), with slightly worse performance. While not the most numerically stable summation technique (e.g., compensated summation techniques such as the Kahan–Babuška-Neumaier algorithm are generally more numerically stable), pairwise summation strikes a reasonable balance between numerical stability and performance. If either numerical stability or performance is more desirable for your use case, consider alternative summation techniques.
132+
- Both functions support array-like objects having getter and setter accessors for array element access (e.g., [`@stdlib/array/base/accessor`][@stdlib/array/base/accessor])
132133
- Depending on the environment, the typed versions ([`dcusumpw`][@stdlib/blas/ext/base/dcusumpw], [`scusumpw`][@stdlib/blas/ext/base/scusumpw], etc.) are likely to be significantly more performant.
133134

134135
</section>
@@ -197,6 +198,8 @@ console.log( y );
197198

198199
[mdn-typed-array]: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/TypedArray
199200

201+
[@stdlib/array/base/accessor]: https://github.com/stdlib-js/stdlib/tree/develop/lib/node_modules/%40stdlib/array/base/accessor
202+
200203
[@higham:1993a]: https://doi.org/10.1137/0914050
201204

202205
<!-- <related-links> -->

Diff for: lib/node_modules/@stdlib/blas/ext/base/gcusumpw/docs/types/index.d.ts

+13-3
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,17 @@
2020

2121
/// <reference types="@stdlib/types"/>
2222

23-
import { NumericArray } from '@stdlib/types/array';
23+
import { NumericArray, Collection, AccessorArrayLike } from '@stdlib/types/array';
24+
25+
/**
26+
* Input array.
27+
*/
28+
type InputArray = NumericArray | Collection<number> | AccessorArrayLike<number>;
29+
30+
/**
31+
* Output array.
32+
*/
33+
type OutputArray = NumericArray | Collection<number> | AccessorArrayLike<number>;
2434

2535
/**
2636
* Interface describing `gcusumpw`.
@@ -44,7 +54,7 @@ interface Routine {
4454
* gcusumpw( x.length, 0.0, x, 1, y, 1 );
4555
* // y => [ 1.0, -1.0, 1.0 ]
4656
*/
47-
( N: number, sum: number, x: NumericArray, strideX: number, y: NumericArray, strideY: number ): NumericArray;
57+
<T extends OutputArray>( N: number, sum: number, x: InputArray, strideX: number, y: T, strideY: number ): T;
4858

4959
/**
5060
* Computes the cumulative sum of strided array elements using pairwise summation and alternative indexing semantics.
@@ -66,7 +76,7 @@ interface Routine {
6676
* gcusumpw.ndarray( x.length, 0.0, x, 1, 0, y, 1, 0 );
6777
* // y => [ 1.0, -1.0, 1.0 ]
6878
*/
69-
ndarray( N: number, sum: number, x: NumericArray, strideX: number, offsetX: number, y: NumericArray, strideY: number, offsetY: number ): NumericArray;
79+
ndarray<T extends OutputArray>( N: number, sum: number, x: InputArray, strideX: number, offsetX: number, y: T, strideY: number, offsetY: number ): T;
7080
}
7181

7282
/**

Diff for: lib/node_modules/@stdlib/blas/ext/base/gcusumpw/docs/types/test.ts

+5-2
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616
* limitations under the License.
1717
*/
1818

19+
import AccessorArray = require( '@stdlib/array/base/accessor' );
1920
import gcusumpw = require( './index' );
2021

2122

@@ -26,7 +27,8 @@ import gcusumpw = require( './index' );
2627
const x = new Float64Array( 10 );
2728
const y = new Float64Array( 10 );
2829

29-
gcusumpw( x.length, 0.0, x, 1, y, 1 ); // $ExpectType NumericArray
30+
gcusumpw( x.length, 0.0, x, 1, y, 1 ); // $ExpectType Float64Array
31+
gcusumpw( x.length, 0.0, new AccessorArray( x ), 1, new AccessorArray( y ), 1 ); // $ExpectType AccessorArray<number>
3032
}
3133

3234
// The compiler throws an error if the function is provided a first argument which is not a number...
@@ -139,7 +141,8 @@ import gcusumpw = require( './index' );
139141
const x = new Float64Array( 10 );
140142
const y = new Float64Array( 10 );
141143

142-
gcusumpw.ndarray( x.length, 0.0, x, 1, 0, y, 1, 0 ); // $ExpectType NumericArray
144+
gcusumpw.ndarray( x.length, 0.0, x, 1, 0, y, 1, 0 ); // $ExpectType Float64Array
145+
gcusumpw.ndarray( x.length, 0.0, new AccessorArray( x ), 1, 0, new AccessorArray( y ), 1, 0 ); // $ExpectType AccessorArray<number>
143146
}
144147

145148
// The compiler throws an error if the `ndarray` method is provided a first argument which is not a number...
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2025 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 floor = require( '@stdlib/math/base/special/floor' );
24+
25+
26+
// VARIABLES //
27+
28+
// Blocksize for pairwise summation:
29+
var BLOCKSIZE = 128;
30+
31+
32+
// MAIN //
33+
34+
/**
35+
* Computes the cumulative sum of strided array elements using pairwise summation.
36+
*
37+
* ## Method
38+
*
39+
* - This implementation uses pairwise summation, which accrues rounding error `O(log2 N)` instead of `O(N)`. The recursion depth is also `O(log2 N)`.
40+
*
41+
* ## References
42+
*
43+
* - Higham, Nicholas J. 1993. "The Accuracy of Floating Point Summation." _SIAM Journal on Scientific Computing_ 14 (4): 783–99. doi:[10.1137/0914050](https://doi.org/10.1137/0914050).
44+
*
45+
* @private
46+
* @param {PositiveInteger} N - number of indexed elements
47+
* @param {number} sum - initial sum
48+
* @param {Object} x - input array object
49+
* @param {Collection} x.data - input array data
50+
* @param {Array<Function>} x.accessors - array element accessors
51+
* @param {integer} strideX - stride length for `x`
52+
* @param {NonNegativeInteger} offsetX - starting index for `x`
53+
* @param {Object} y - output array object
54+
* @param {Collection} y.data - output array data
55+
* @param {Array<Function>} y.accessors - array element accessors
56+
* @param {integer} strideY - stride length for `y`
57+
* @param {NonNegativeInteger} offsetY - starting index for `y`
58+
* @returns {Object} output array object
59+
*
60+
* @example
61+
* var toAccessorArray = require( '@stdlib/array/base/to-accessor-array' );
62+
* var arraylike2object = require( '@stdlib/array/base/arraylike2object' );
63+
*
64+
* var x = [ 2.0, 1.0, 2.0, -2.0, -2.0, 2.0, 3.0, 4.0 ];
65+
* var y = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ];
66+
*
67+
* gcusumpw( 4, 0.0, arraylike2object( toAccessorArray( x ) ), 2, 1, arraylike2object( toAccessorArray( y ) ), 1, 0 );
68+
* // y => [ 1.0, -1.0, 1.0, 5.0, 0.0, 0.0, 0.0, 0.0 ]
69+
*/
70+
function gcusumpw( N, sum, x, strideX, offsetX, y, strideY, offsetY ) {
71+
var xbuf;
72+
var ybuf;
73+
var xget;
74+
var yget;
75+
var yset;
76+
var ix;
77+
var iy;
78+
var s;
79+
var n;
80+
var i;
81+
82+
// Cache reference to array data:
83+
xbuf = x.data;
84+
ybuf = y.data;
85+
86+
// Cache reference to the element accessors:
87+
xget = x.accessors[ 0 ];
88+
yget = y.accessors[ 0 ];
89+
yset = y.accessors[ 1 ];
90+
91+
ix = offsetX;
92+
iy = offsetY;
93+
if ( N <= BLOCKSIZE ) {
94+
s = 0.0;
95+
for ( i = 0; i < N; i++ ) {
96+
s += xget( xbuf, ix );
97+
yset( ybuf, iy, sum + s);
98+
ix += strideX;
99+
iy += strideY;
100+
}
101+
return y;
102+
}
103+
n = floor( N/2 );
104+
gcusumpw( n, sum, x, strideX, ix, y, strideY, iy );
105+
iy += (n-1) * strideY;
106+
gcusumpw( N-n, yget( ybuf, iy ), x, strideX, ix+(n*strideX), y, strideY, iy+strideY ); // eslint-disable-line max-len
107+
return y;
108+
}
109+
110+
111+
// EXPORTS //
112+
113+
module.exports = gcusumpw;

Diff for: lib/node_modules/@stdlib/blas/ext/base/gcusumpw/lib/ndarray.js

+10
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,9 @@
2020

2121
// MODULES //
2222

23+
var arraylike2object = require( '@stdlib/array/base/arraylike2object' );
2324
var floor = require( '@stdlib/math/base/special/floor' );
25+
var accessors = require( './accessors.js' );
2426

2527

2628
// VARIABLES //
@@ -62,13 +64,21 @@ var BLOCKSIZE = 128;
6264
function gcusumpw( N, sum, x, strideX, offsetX, y, strideY, offsetY ) {
6365
var ix;
6466
var iy;
67+
var ox;
68+
var oy;
6569
var s;
6670
var n;
6771
var i;
6872

6973
if ( N <= 0 ) {
7074
return y;
7175
}
76+
ox = arraylike2object( x );
77+
oy = arraylike2object( y );
78+
if ( ox.accessorProtocol || oy.accessorProtocol ) {
79+
accessors( N, sum, ox, strideX, offsetX, oy, strideY, offsetY );
80+
return y;
81+
}
7282
ix = offsetX;
7383
iy = offsetY;
7484
if ( N <= BLOCKSIZE ) {

0 commit comments

Comments
 (0)