Skip to content

Commit c305bbc

Browse files
committed
feat: add base implementation
--- 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: na - 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: na - 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 22f1408 commit c305bbc

File tree

2 files changed

+183
-0
lines changed

2 files changed

+183
-0
lines changed
Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,116 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2024 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 dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray;
24+
var sign = require( '@stdlib/math/base/special/copysign' );
25+
var dlamch = require( '@stdlib/lapack/base/dlamch' );
26+
var abs = require( '@stdlib/math/base/special/abs' );
27+
var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
28+
var dlapy2 = require( './dlapy2.js' );
29+
30+
31+
// MAIN //
32+
33+
/**
34+
* Generates a real elementary reflector `H` of order `n` such that applying `H` to a vector `[alpha; x]` zeros out `x`.
35+
*
36+
* ## Notes
37+
*
38+
* - `H` is a Householder matrix with the form `H = I - tau * [1; v] * [1, v^T]`, where `tau` is a scalar and `v` is a vector.
39+
* - the input vector is `[alpha; x]`, where `alpha` is a scalar and `x` is a real `(n-1)`-element vector.
40+
* - the result of applying `H` to `[alpha; x]` is `[beta; 0]`, with `beta` being a scalar and the rest of the vector zeroed.
41+
* - if all elements of `x` are zero, then `tau = 0` and `H` is the identity matrix.
42+
* - otherwise, `1 <= tau <= 2` and `H` is orthogonal, i.e., `H^T * H = I`.
43+
*
44+
* @private
45+
* @param {NonNegativeInteger} N - order of matrix `A`
46+
* @param {Float64Array} X - overwritten by the vector `V` on exit
47+
* @param {integer} strideX - stride length for `X`
48+
* @param {NonNegativeInteger} offsetX - starting index of `X`
49+
* @param {Float64Array} out - array to store `alpha` and `tau`, first indexed element stores `alpha` and the second indexed element stores `tau`
50+
* @param {integer} strideOut - stride length for `out`
51+
* @param {NonNegativeInteger} offsetOut - starting index of `out`
52+
* @returns {void} overwrites the array `X` and `out` in place
53+
*
54+
* @example
55+
* var Float64Array = require( '@stdlib/array/float64' );
56+
*
57+
* var X = new Float64Array( [ 2.0, 3.0, 4.0 ] );
58+
* var out = new Float64Array( [ 4.0, 0.0 ] );
59+
*
60+
* dlarfg( 4, X, 1, 0, out, 1, 0 );
61+
* // X => <Float64Array>[ ~0.19, ~0.28, ~0.37 ]
62+
* // out => <Float64Array>[ ~6.7, ~1.6 ]
63+
*/
64+
function dlarfg( N, X, strideX, offsetX, out, strideOut, offsetOut ) {
65+
var safemin;
66+
var rsafmin;
67+
var xnorm;
68+
var alpha;
69+
var beta;
70+
var tau;
71+
var knt;
72+
var i;
73+
74+
if ( N <= 1 ) {
75+
tau = 0.0; // tau = 0.0
76+
out[ offsetOut + strideOut ] = tau;
77+
return;
78+
}
79+
80+
xnorm = dnrm2( N - 1, X, strideX, offsetX );
81+
alpha = out[ offsetOut ];
82+
83+
if ( xnorm === 0.0 ) {
84+
tau = 0.0; // tau = 0.0
85+
out[ strideOut + offsetOut ] = tau;
86+
} else {
87+
beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
88+
safemin = dlamch( 'S' ) / dlamch( 'E' );
89+
knt = 0;
90+
if ( abs( beta ) < safemin ) {
91+
rsafmin = 1.0 / safemin;
92+
while ( abs( beta ) < safemin && knt < 20 ) {
93+
knt += 1;
94+
dscal( N-1, rsafmin, X, strideX, offsetX );
95+
beta *= rsafmin;
96+
alpha *= rsafmin; // alpha *= rsafmin
97+
}
98+
xnorm = dnrm2( N - 1, X, strideX, offsetX );
99+
beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
100+
}
101+
tau = ( beta - alpha ) / beta; // tau = (beta - alpha) / beta
102+
dscal( N-1, 1.0 / ( alpha - beta ), X, strideX, offsetX );
103+
for ( i = 0; i < knt; i++ ) {
104+
beta *= safemin;
105+
}
106+
alpha = beta;
107+
108+
out[ offsetOut ] = alpha;
109+
out[ strideOut + offsetOut ] = tau;
110+
}
111+
}
112+
113+
114+
// EXPORTS //
115+
116+
module.exports = dlarfg;
Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2024 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 abs = require( '@stdlib/math/base/special/abs' );
24+
var min = require( '@stdlib/math/base/special/min' );
25+
var max = require( '@stdlib/math/base/special/max' );
26+
var sqrt = require( '@stdlib/math/base/special/sqrt' );
27+
var pow = require( '@stdlib/math/base/special/pow' );
28+
29+
30+
// MAIN //
31+
32+
/**
33+
* Returns sqrt( x^2 + y^2 ) in a manner which doesn't cause unnecessary overflow.
34+
*
35+
* @private
36+
* @param {number} x - order of matrix `A`
37+
* @param {number} y - order of matrix `A`
38+
* @returns {integer} status code
39+
*
40+
* @example
41+
* var out = dlapy2( 3.0, 4.0, );
42+
* // returns 5
43+
*/
44+
function dlapy2( x, y ) {
45+
// TODO - make a separate package for this
46+
var xabs;
47+
var yabs;
48+
var w;
49+
var z;
50+
51+
xabs = abs( x );
52+
yabs = abs( y );
53+
54+
w = max( xabs, yabs );
55+
z = min( xabs, yabs );
56+
57+
if ( z === 0.0 ) {
58+
return w;
59+
}
60+
61+
return w * ( sqrt( 1.0 + pow( z / w, 2.0 ) ) );
62+
}
63+
64+
65+
// EXPORTS //
66+
67+
module.exports = dlapy2;

0 commit comments

Comments
 (0)