Skip to content

Commit 2e1ec09

Browse files
committed
refactor: add comments and optimise pointer arithmetic
--- 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: passed - task: lint_javascript_tests status: passed - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 3d5761d commit 2e1ec09

File tree

6 files changed

+235
-115
lines changed

6 files changed

+235
-115
lines changed

lib/node_modules/@stdlib/lapack/base/dlacn2/README.md

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -242,16 +242,30 @@ var V = new Float64Array( 4 );
242242

243243
var work = new Float64Array( 4 );
244244

245+
// Start the estimation with `KASE` = 0
245246
while ( true ) {
246247
dlacn2( 4, V, X, ISGN, EST, KASE, ISAVE );
247248

248249
if ( KASE[ 0 ] === 0 ) {
250+
// Estimation is complete
249251
break;
250252
}
251253
else if ( KASE[ 0 ] === 1 ) {
254+
/*
255+
* - If `KASE` == 1: `dlacn2` wants us to compute A * X and overwrite X.
256+
* - `dgemv` performs the matrix-vector product: work = A * X.
257+
* - A workspace array is necessary here because of the function signature of `dgemv`.
258+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
259+
*/
252260
dgemv( 'row-major', 'no-transpose', shape[ 0 ], shape[ 1 ], 1.0, A, strides[ 0 ], X, 1, 0, work, 1 );
253261
dcopy( shape[ 0 ], work, 1, X, 1 );
254262
} else if ( KASE[ 0 ] === 2 ) {
263+
/*
264+
* - If `KASE` == 2: `dlacn2` wants us to compute A ^ T * X and overwrite X.
265+
* - `dgemv` performs the matrix-vector product: work = A ^ T * X.
266+
* - A workspace array is necessary here because of the function signature of `dgemv`.
267+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
268+
*/
255269
dgemv( 'row-major', 'transpose', shape[ 0 ], shape[ 1 ], 1.0, A, strides[ 0 ], X, 1, 0, work, 1 );
256270
dcopy( shape[ 0 ], work, 1, X, 1 );
257271
}

lib/node_modules/@stdlib/lapack/base/dlacn2/examples/index.js

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,16 +52,30 @@ var V = new Float64Array( 4 );
5252

5353
var work = new Float64Array( 4 );
5454

55+
// Start the estimation with `KASE` = 0
5556
while ( true ) {
5657
dlacn2( 4, V, X, ISGN, EST, KASE, ISAVE );
5758

5859
if ( KASE[ 0 ] === 0 ) {
60+
// Estimation is complete
5961
break;
6062
}
6163
else if ( KASE[ 0 ] === 1 ) {
64+
/*
65+
* - If `KASE` == 1: `dlacn2` wants us to compute A * X and overwrite X.
66+
* - `dgemv` performs the matrix-vector product: work = A * X.
67+
* - A workspace array is necessary here because of the function signature of `dgemv`.
68+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
69+
*/
6270
dgemv( 'row-major', 'no-transpose', shape[ 0 ], shape[ 1 ], 1.0, A, strides[ 0 ], X, 1, 0, work, 1 );
6371
dcopy( shape[ 0 ], work, 1, X, 1 );
6472
} else if ( KASE[ 0 ] === 2 ) {
73+
/*
74+
* - If `KASE` == 2: `dlacn2` wants us to compute A ^ T * X and overwrite X.
75+
* - `dgemv` performs the matrix-vector product: work = A ^ T * X.
76+
* - A workspace array is necessary here because of the function signature of `dgemv`.
77+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
78+
*/
6579
dgemv( 'row-major', 'transpose', shape[ 0 ], shape[ 1 ], 1.0, A, strides[ 0 ], X, 1, 0, work, 1 );
6680
dcopy( shape[ 0 ], work, 1, X, 1 );
6781
}

lib/node_modules/@stdlib/lapack/base/dlacn2/lib/base.js

Lines changed: 27 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -26,10 +26,11 @@ var abs = require( '@stdlib/math/base/special/abs' );
2626
var idamax = require( '@stdlib/blas/base/idamax' ).ndarray;
2727
var dcopy = require( '@stdlib/blas/base/dcopy' ).ndarray;
2828
var dasum = require( '@stdlib/blas/base/dasum' ).ndarray;
29-
var nint = require( './nint.js' );
29+
var dfill = require( '@stdlib/blas/ext/base/dfill' ).ndarray;
30+
var nint = require( '@stdlib/math/base/special/round-nearest-even' );
3031

3132

32-
// MAIN //
33+
// FUNCTIONS //
3334

3435
/**
3536
* Applies a deterministic fallback vector for final evaluation.
@@ -175,19 +176,17 @@ function isaveIsOne( N, V, offsetV, X, strideX, offsetX, ISGN, strideISGN, offse
175176
* // KASE => <Int32Array>[ 1 ]
176177
*/
177178
function isaveIsTwo( N, X, strideX, offsetX, KASE, offsetKASE, ISAVE, strideISAVE, offsetISAVE ) {
178-
var ix;
179-
var i;
179+
var xmax;
180+
var i1;
180181

181-
ISAVE[ offsetISAVE + strideISAVE ] = idamax( N, X, strideX, offsetX );
182-
ISAVE[ offsetISAVE + ( 2 * strideISAVE ) ] = 2;
182+
i1 = offsetISAVE + strideISAVE;
183+
ISAVE[ i1 ] = idamax( N, X, strideX, offsetX ); // stores the index of the max element in X
184+
ISAVE[ i1 + strideISAVE ] = 2;
185+
xmax = offsetX + ( ISAVE[ i1 ] * strideX ); // pointer to the max element in X
183186

184-
ix = offsetX;
185-
for ( i = 0; i < N; i++ ) {
186-
X[ ix ] = 0.0;
187-
ix += strideX;
188-
}
187+
dfill( N, 0.0, X, strideX, offsetX );
189188

190-
X[ offsetX + ( ISAVE[ offsetISAVE + strideISAVE ] * strideX ) ] = 1.0;
189+
X[ xmax ] = 1.0;
191190
KASE[ offsetKASE ] = 1;
192191
ISAVE[ offsetISAVE ] = 3;
193192
}
@@ -315,21 +314,23 @@ function isaveIsThree( N, X, strideX, offsetX, V, strideV, offsetV, EST, offsetE
315314
* // KASE => <Int32Array>[ 1 ]
316315
*/
317316
function isaveIsFour( N, X, strideX, offsetX, ISAVE, strideISAVE, offsetISAVE, KASE, offsetKASE ) {
317+
var prevxmax;
318318
var jlast;
319-
var ix;
320-
var i;
319+
var xmax;
320+
var i2;
321321

322322
jlast = ISAVE[ offsetISAVE + strideISAVE ];
323+
prevxmax = offsetX + ( jlast * strideX ); // points to X[ isave(1) ]
323324
ISAVE[ offsetISAVE + strideISAVE ] = idamax( N, X, strideX, offsetX );
325+
xmax = offsetX + ( ISAVE[ offsetISAVE + strideISAVE ] * strideX ); // points to the largest value in X
326+
i2 = offsetISAVE + ( 2 * strideISAVE ); // points to isave(2), the number of refinement iterations
324327

325-
if ( X[ offsetX + ( jlast * strideX ) ] !== abs( X[ offsetX + ( ISAVE[ offsetISAVE + strideISAVE ] * strideX ) ] ) && ISAVE[ offsetISAVE + ( 2 * strideISAVE ) ] < 5 ) {
326-
ISAVE[ offsetISAVE + ( 2 * strideISAVE ) ] += 1;
327-
ix = offsetX;
328-
for ( i = 0; i < N; i++ ) {
329-
X[ ix ] = 0.0;
330-
ix += strideX;
331-
}
332-
X[ offsetX + ( ISAVE[ offsetISAVE + strideISAVE ] * strideX ) ] = 1.0;
328+
if ( X[ prevxmax ] !== abs( X[ xmax ] ) && ISAVE[ i2 ] < 5 ) {
329+
ISAVE[ i2 ] += 1;
330+
331+
dfill( N, 0.0, X, strideX, offsetX );
332+
333+
X[ xmax ] = 1.0;
333334
KASE[ offsetKASE ] = 1;
334335
ISAVE[ offsetISAVE ] = 3;
335336
return;
@@ -385,6 +386,9 @@ function isaveIsFive( N, V, strideV, offsetV, X, strideX, offsetX, EST, offsetES
385386
KASE[ offsetKASE ] = 0;
386387
}
387388

389+
390+
// MAIN //
391+
388392
/**
389393
* Estimates the one-norm of a square matrix `A`, using reverse communication for evaluating matrix-vector products.
390394
*
@@ -442,15 +446,8 @@ function isaveIsFive( N, V, strideV, offsetV, X, strideX, offsetX, EST, offsetES
442446
* // KASE => <Int32Array>[ 1 ]
443447
*/
444448
function dlacn2( N, V, strideV, offsetV, X, strideX, offsetX, ISGN, strideISGN, offsetISGN, EST, offsetEST, KASE, offsetKASE, ISAVE, strideISAVE, offsetISAVE ) {
445-
var ix;
446-
var i;
447-
448449
if ( KASE[ offsetKASE ] === 0 ) {
449-
ix = offsetX;
450-
for ( i = 0; i < N; i++ ) {
451-
X[ ix ] = 1 / N;
452-
ix += strideX;
453-
}
450+
dfill( N, 1 / N, X, strideX, offsetX );
454451
KASE[ offsetKASE ] = 1;
455452
ISAVE[ offsetISAVE ] = 1;
456453
return;

lib/node_modules/@stdlib/lapack/base/dlacn2/lib/nint.js

Lines changed: 0 additions & 65 deletions
This file was deleted.

lib/node_modules/@stdlib/lapack/base/dlacn2/test/test.dlacn2.js

Lines changed: 46 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -43,6 +43,7 @@ tape( 'the function has an arity of 7', function test( t ) {
4343
t.end();
4444
});
4545

46+
// Expected values determined via the reference Fortran routines in LAPACK:
4647
tape( 'the function returns expected values when calculating the one-norm of a square matrix for N = 1', function test( t ) {
4748
var expectedISAVE;
4849
var expectedKASE;
@@ -69,16 +70,29 @@ tape( 'the function returns expected values when calculating the one-norm of a s
6970

7071
work = new Float64Array( 1 );
7172

73+
// Start the estimation with `KASE` = 0:
7274
while ( true ) {
7375
dlacn2( 1, V, X, ISGN, EST, KASE, ISAVE );
7476

7577
if ( KASE[ 0 ] === 0 ) {
78+
// Estimation is complete:
7679
break;
77-
}
78-
else if ( KASE[ 0 ] === 1 ) {
80+
} else if ( KASE[ 0 ] === 1 ) {
81+
/*
82+
* - If `KASE` == 1: `dlacn2` wants us to compute A * X and overwrite X.
83+
* - `dgemv` performs the matrix-vector product: work = A * X.
84+
* - A workspace array is necessary here because of the function signature of `dgemv`.
85+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
86+
*/
7987
dgemv( 'row-major', 'no-transpose', 1, 1, 1.0, A, 1, X, 1, 0, work, 1 );
8088
dcopy( 1, work, 1, X, 1 );
8189
} else if ( KASE[ 0 ] === 2 ) {
90+
/*
91+
* - If `KASE` == 2: `dlacn2` wants us to compute A ^ T * X and overwrite X.
92+
* - `dgemv` performs the matrix-vector product: work = A ^ T * X.
93+
* - A workspace array is necessary here because of the function signature of `dgemv`.
94+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
95+
*/
8296
dgemv( 'row-major', 'transpose', 1, 1, 1.0, A, 1, X, 1, 0, work, 1 );
8397
dcopy( 1, work, 1, X, 1 );
8498
}
@@ -130,16 +144,29 @@ tape( 'the function returns expected values when calculating the one-norm of a s
130144

131145
work = new Float64Array( 4 );
132146

147+
// Start the estimation with `KASE` = 0:
133148
while ( true ) {
134149
dlacn2( 4, V, X, ISGN, EST, KASE, ISAVE );
135150

136151
if ( KASE[ 0 ] === 0 ) {
152+
// Estimation is complete:
137153
break;
138-
}
139-
else if ( KASE[ 0 ] === 1 ) {
154+
} else if ( KASE[ 0 ] === 1 ) {
155+
/*
156+
* - If `KASE` == 1: `dlacn2` wants us to compute A * X and overwrite X.
157+
* - `dgemv` performs the matrix-vector product: work = A * X.
158+
* - A workspace array is necessary here because of the function signature of `dgemv`.
159+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
160+
*/
140161
dgemv( 'row-major', 'no-transpose', 4, 4, 1.0, A, 4, X, 1, 0, work, 1 );
141162
dcopy( 4, work, 1, X, 1 );
142163
} else if ( KASE[ 0 ] === 2 ) {
164+
/*
165+
* - If `KASE` == 2: `dlacn2` wants us to compute A ^ T * X and overwrite X.
166+
* - `dgemv` performs the matrix-vector product: work = A ^ T * X.
167+
* - A workspace array is necessary here because of the function signature of `dgemv`.
168+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
169+
*/
143170
dgemv( 'row-major', 'transpose', 4, 4, 1.0, A, 4, X, 1, 0, work, 1 );
144171
dcopy( 4, work, 1, X, 1 );
145172
}
@@ -191,16 +218,29 @@ tape( 'the function returns expected values when calculating the one-norm of a s
191218

192219
work = new Float64Array( 4 );
193220

221+
// Start the estimation with `KASE` = 0:
194222
while ( true ) {
195223
dlacn2( 4, V, X, ISGN, EST, KASE, ISAVE );
196224

197225
if ( KASE[ 0 ] === 0 ) {
226+
// Estimation is complete:
198227
break;
199-
}
200-
else if ( KASE[ 0 ] === 1 ) {
228+
} else if ( KASE[ 0 ] === 1 ) {
229+
/*
230+
* - If `KASE` == 1: `dlacn2` wants us to compute A * X and overwrite X.
231+
* - `dgemv` performs the matrix-vector product: work = A * X.
232+
* - A workspace array is necessary here because of the function signature of `dgemv`.
233+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
234+
*/
201235
dgemv( 'row-major', 'no-transpose', 4, 4, 1.0, A, 4, X, 1, 0, work, 1 );
202236
dcopy( 4, work, 1, X, 1 );
203237
} else if ( KASE[ 0 ] === 2 ) {
238+
/*
239+
* - If `KASE` == 2: `dlacn2` wants us to compute A ^ T * X and overwrite X.
240+
* - `dgemv` performs the matrix-vector product: work = A ^ T * X.
241+
* - A workspace array is necessary here because of the function signature of `dgemv`.
242+
* - `dcopy` copies work into X, so that the updated X is passed back to `dlacn2`.
243+
*/
204244
dgemv( 'row-major', 'transpose', 4, 4, 1.0, A, 4, X, 1, 0, work, 1 );
205245
dcopy( 4, work, 1, X, 1 );
206246
}

0 commit comments

Comments
 (0)