Skip to content

feat: add lapack/base/dlabrd #12999

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

Draft
prajjwalbajpai wants to merge 1 commit into
base: develop
Choose a base branch
from
+340 −0
Draft

feat: add lapack/base/dlabrd #12999

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
198 changes: 198 additions & 0 deletions lib/node_modules/@stdlib/lapack/base/dlabrd/lib/base.js
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,198 @@
/**
* @license Apache-2.0
*
* Copyright (c) 2026 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.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

/* eslint-disable max-len, max-params */

'use strict';

// MODULES //

var dgemv = require( '@stdlib/blas/base/dgemv' ).ndarray;
var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
var min = require( '@stdlib/math/base/special/min' );
var dlarfg = require( './dlarfg.js' );


// MAIN //

/**
* Reduces the first `NB` rows and columns of a real general matrix `A` to upper or lower bi-diagonal form by an orthogonal transformation `Q**T*A*P`.
*
* ## Notes
*
* - If `M >= N`,
*
* - `A` is reduced to upper bi-diagonal form.
* - Elements on and below the diagonal in the first `NB` columns, with the array `TAUQ`, represent the orthogonal matrix `Q` as a product of elementary reflectors.
* - Elements above the diagonal in the first `NB` rows, with the array `TAUP`, represent the orthogonal matrix `P` as a product of elementary reflectors.
*
* - If `M < N`,
*
* - `A` is reduced to lower bi-diagonal form.
* - Elements below the diagonal in the first `NB` columns, with the array `TAUQ`, represent the orthogonal matrix `Q` as a product of elementary reflectors.
* - Elements on and above the diagonal in the first `NB` rows, with the array `TAUP`, represent the orthogonal matrix `P` as a product of elementary reflectors.
*
* @private
* @param {NonNegativeInteger} M - number of rows in `A`
* @param {NonNegativeInteger} N - number of columns in `A`
* @param {integer} NB - number of leading rows and columns of `A` to reduce
* @param {Float64Array} A - input matrix
* @param {integer} strideA1 - stride of the first dimension of `A`
* @param {integer} strideA2 - stride of the second dimension of `A`
* @param {NonNegativeInteger} offsetA - starting index for `A`
* @param {Float64Array} D - real diagonal elements (length `NB`)
* @param {integer} strideD - stride length for `D`
* @param {NonNegativeInteger} offsetD - starting index for `D`
* @param {Float64Array} E - real off-diagonal elements (length `NB`)
* @param {integer} strideE - stride length for `E`
* @param {NonNegativeInteger} offsetE - starting index for `E`
* @param {Float64Array} TAUQ - scalars factors of the elementary reflectors that represent the orthogonal matrix `Q` (length `NB`)
* @param {integer} strideTAUQ - stride length for `TAUQ`
* @param {NonNegativeInteger} offsetTAUQ - starting index for `TAUQ`
* @param {Float64Array} TAUP - scalars factors of the elementary reflectors that represent the orthogonal matrix `P` (length `NB`)
* @param {integer} strideTAUP - stride length for `TAUP`
* @param {NonNegativeInteger} offsetTAUP - starting index for `TAUP`
* @param {Float64Array} X - output matrix
* @param {integer} strideX1 - stride of the first dimension of `X`
* @param {integer} strideX2 - stride of the second dimension of `X`
* @param {NonNegativeInteger} offsetX - starting index for `X`
* @param {Float64Array} Y - output matrix
* @param {integer} strideY1 - stride of the first dimension of `Y`
* @param {integer} strideY2 - stride of the second dimension of `Y`
* @param {NonNegativeInteger} offsetY - starting index for `Y`
*/
function dlabrd( M, N, NB, A, strideA1, strideA2, offsetA, D, strideD, offsetD, E, strideE, offsetE, TAUQ, strideTAUQ, offsetTAUQ, TAUP, strideTAUP, offsetTAUP, X, strideX1, strideX2, offsetX, Y, strideY1, strideY2, offsetY ) {
var i;

if ( M <= 0 || N <= 0 ) {
return;
}

if ( M >= N ) {
// Reduce to upper bi-diagonal form
for ( i = 0; i < NB; i++ ) {
// Update A(i:M-1, i)
dgemv( 'no-transpose', M - i, i, -1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ), Y, strideY2, offsetY + ( i * strideY1 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) );

dgemv( 'no-transpose', M - i, i, -1.0, X, strideX1, strideX2, offsetX + ( i * strideX1 ), A, strideA1, offsetA + ( i * strideA2 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) );

// Generate elementary reflector H(i) to annihilate A(i+1:M-1, i)
dlarfg( M - i, A, offsetA + ( i * strideA1 ) + ( i * strideA2 ), A, strideA1, offsetA + ( ( i + 1, M - 1 ) * strideA1 ) + ( i * strideA2 ), TAUQ, offsetTAUQ + ( i * strideTAUQ ) );
A[ offsetD + ( i * strideD ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ];

if ( i < N - 1 ) {
A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ] = 1.0;

// Compute Y(i+1:N-1, i)
dgemv( 'transpose', M - i, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ) + ( ( i + 1 ) * strideA2 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dgemv( 'transpose', M - i, i, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );

dgemv( 'no-transpose', N - i - 1, i, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dgemv( 'transpose', M - i, i, 1.0, X, strideX1, strideX2, offsetX + ( i * strideX1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );

dgemv( 'transpose', i, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dscal( N - i - 1, TAUQ[ offsetTAUQ + ( i * strideTAUQ ) ], Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

// Update A(i, i+1:N-1)
dgemv( 'no-transpose', N - i - 1, i + 1, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), A, strideA2, offsetA + ( i * strideA1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 );

dgemv( 'transpose', i, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), X, strideX2, offsetX + ( i * strideX1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 );

// Generate elementary reflector P(i) to annihilate A(i, i+2:N-1)
dlarfg( N - i - 1, A, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, A, strideA2, offsetA + ( i * strideA1 ) + ( min( i + 2, N - 1 ) * strideA2 ), TAUP, offsetTAUP + ( i * strideTAUP ) );
E[ offsetE + ( i * strideE ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 ];
A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2 ] = 1.0;

// Compute X(i+1:M-1, i)

dgemv( 'no-transpose', M - i - 1, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * ( strideA1 + strideA2 ) ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dgemv( 'transpose', N - i - 1, i + 1, 1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * strideX2 ) );

dgemv( 'no-transpose', M - i - 1, i + 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dgemv( 'no-transpose', i, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA2, 0.0, X, strideX1, offsetX + ( i * strideX2 ) );

dgemv( 'no-transpose', M - i - 1, i, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dscal( M - i - 1, TAUP[ offsetTAUP + ( i * strideTAUP ) ], X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );
}
}
} else {
// Reduce to lower bi-diagonal form (M < N)
for ( i = 0; i < NB; i++ ) {
// Update A(i, i:N-1)
dgemv( 'no-transpose', N - i, i, -1.0, Y, strideY1, strideY2, offsetY + ( i * strideY1 ), A, strideA2, offsetA + ( i * strideA1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) );

dgemv( 'transpose', i, N - i, -1.0, A, strideA1, strideA2, offsetA + ( i * strideA2 ), X, strideX2, offsetX + ( i * strideX1 ), 1.0, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) );

// Generate elementary reflector P(i) to annihilate A(i, i+1:N-1)
dlarfg( N - i, A, offsetA + ( i * ( strideA1 + strideA2 ) ), A, strideA2, offsetA + ( i * strideA1 ) + ( min( i + 1, N - 1 ) * strideA2 ), TAUP, offsetTAUP + ( i * strideTAUP ) );
D[ offsetD + ( i * strideD ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ];

if ( i < M - 1 ) {
// Set A(i,i) = 1 for use as reflector vector
A[ offsetA + ( i * ( strideA1 + strideA2 ) ) ] = 1.0;

// Compute X(i+1:M-1, i)
dgemv( 'no-transpose', M - i - 1, N - i, 1.0, A, strideA1, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dgemv( 'transpose', N - i, i, 1.0, Y, strideY1, strideY2, offsetY + ( i * strideY1 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * strideX2 ) );

dgemv( 'no-transpose', M - i - 1, i, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dgemv( 'no-transpose', i, N - i, 1.0, A, strideA1, strideA2, offsetA + ( i * strideA2 ), A, strideA2, offsetA + ( i * ( strideA1 + strideA2 ) ), 0.0, X, strideX1, offsetX + ( i * strideX2 ) );

dgemv( 'no-transpose', M - i - 1, i, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), X, strideX1, offsetX + ( i * strideX2 ), 1.0, X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

dscal( M - i - 1, TAUP[ offsetTAUP + ( i * strideTAUP ) ], X, strideX1, offsetX + ( i * ( strideX1 + strideX2 ) ) + strideX1 );

// Update A(i+1:M-1, i)
dgemv( 'no-transpose', M - i - 1, i, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), Y, strideY2, offsetY + ( i * strideY1 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 );

dgemv( 'no-transpose', M - i - 1, i + 1, -1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), A, strideA1, offsetA + ( i * strideA2 ), 1.0, A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 );

dlarfg( M - i - 1, A, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, A, strideA1, offsetA + ( min( i + 2, M - 1 ) * strideA1 ) + ( i * strideA2 ), TAUQ, offsetTAUQ + ( i * strideTAUQ ) );

E[ offsetE + ( i * strideE ) ] = A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 ];
A[ offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1 ] = 1.0;

// Compute Y(i+1:N-1, i)
dgemv( 'transpose', M - i - 1, N - i - 1, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * ( strideA1 + strideA2 ) ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dgemv( 'transpose', M - i - 1, i, 1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );

dgemv( 'no-transpose', N - i - 1, i, -1.0, Y, strideY1, strideY2, offsetY + ( ( i + 1 ) * strideY1 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dgemv( 'transpose', M - i - 1, i + 1, 1.0, X, strideX1, strideX2, offsetX + ( ( i + 1 ) * strideX1 ), A, strideA1, offsetA + ( i * ( strideA1 + strideA2 ) ) + strideA1, 0.0, Y, strideY1, offsetY + ( i * strideY2 ) );

dgemv( 'transpose', i + 1, N - i - 1, -1.0, A, strideA1, strideA2, offsetA + ( ( i + 1 ) * strideA2 ), Y, strideY1, offsetY + ( i * strideY2 ), 1.0, Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );

dscal( N - i - 1, TAUQ[ offsetTAUQ + ( i * strideTAUQ ) ], Y, strideY1, offsetY + ( i * ( strideY1 + strideY2 ) ) + strideY1 );
}
}
}
}


// EXPORTS //

module.exports = dlabrd;
142 changes: 142 additions & 0 deletions lib/node_modules/@stdlib/lapack/base/dlabrd/lib/dlarfg.js
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
/**
* @license Apache-2.0
*
* Copyright (c) 2026 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.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/

'use strict';

// MODULES //

var dnrm2 = require( '@stdlib/blas/base/dnrm2' ).ndarray;
var sign = require( '@stdlib/math/base/special/copysign' );
var dlamch = require( '@stdlib/lapack/base/dlamch' );
var abs = require( '@stdlib/math/base/special/abs' );
var dscal = require( '@stdlib/blas/base/dscal' ).ndarray;
var dlapy2 = require( '@stdlib/lapack/base/dlapy2' );


// MAIN //

/**
* Generates a real elementary reflector `H` of order `N` such that applying `H` to a vector `[alpha; X]` zeros out `X`.
*
* `H` is a Householder matrix with the form:
*
* ```tex
* H \cdot \begin{bmatrix} \alpha \\ x \end{bmatrix} = \begin{bmatrix} \beta \\ 0 \end{bmatrix}, \quad \text{and} \quad H^T H = I
* ```
*
* where:
*
* - `tau` is a scalar
* - `X` is a vector of length `N-1`
* - `beta` is a scalar value
* - `H` is an orthogonal matrix known as a Householder reflector.
*
* The reflector `H` is constructed in the form:
*
* ```tex
* H = I - \tau \begin{bmatrix}1 \\ v \end{bmatrix} \begin{bmatrix}1 & v^T \end{bmatrix}
* ```
*
* where:
*
* - `tau` is a real scalar
* - `V` is a real vector of length `N-1` that defines the Householder vector
* - The vector `[1; V]` is the Householder direction.
*
* The values of `tau` and `V` are chosen so that applying `H` to the vector `[alpha; X]` results in a new vector `[beta; 0]`, i.e., only the first component remains nonzero. The reflector matrix `H` is symmetric and orthogonal, satisfying `H^T = H` and `H^T H = I`
*
* ## Special cases
*
* - If all elements of `X` are zero, then `tau = 0` and `H = I`, the identity matrix.
* - Otherwise, `tau` satisfies `1 ≤ tau ≤ 2`, ensuring numerical stability in transformations.
*
* ## Notes
*
* - `X` should have `N-1` indexed elements
* - The output array contains the following two elements: `alpha` and `tau`
*
* @private
* @param {NonNegativeInteger} N - number of rows/columns of the elementary reflector `H`
* @param {Float64Array} X - input vector
* @param {integer} strideX - stride length for `X`
* @param {NonNegativeInteger} offsetX - starting index of `X`
* @param {Float64Array} out - output array
* @param {integer} strideOut - stride length for `out`
* @param {NonNegativeInteger} offsetOut - starting index of `out`
* @returns {void}
*
* @example
* var Float64Array = require( '@stdlib/array/float64' );
*
* var X = new Float64Array( [ 2.0, 3.0, 4.0 ] );
* var out = new Float64Array( [ 4.0, 0.0 ] );
*
* dlarfg( 4, X, 1, 0, out, 1, 0 );
* // X => <Float64Array>[ ~0.19, ~0.28, ~0.37 ]
* // out => <Float64Array>[ ~-6.7, ~1.6 ]
*/
function dlarfg( N, X, strideX, offsetX, out, strideOut, offsetOut ) {
var safemin;
var rsafmin;
var xnorm;
var alpha;
var beta;
var tau;
var knt;
var i;

if ( N <= 1 ) {
out[ offsetOut + strideOut ] = 0.0;
return;
}

xnorm = dnrm2( N - 1, X, strideX, offsetX );
alpha = out[ offsetOut ];

if ( xnorm === 0.0 ) {
out[ strideOut + offsetOut ] = 0.0;
} else {
beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
safemin = dlamch( 'safemin' ) / dlamch( 'epsilon' );
knt = 0;
if ( abs( beta ) < safemin ) {
rsafmin = 1.0 / safemin;
while ( abs( beta ) < safemin && knt < 20 ) {
knt += 1;
dscal( N-1, rsafmin, X, strideX, offsetX );
beta *= rsafmin;
alpha *= rsafmin;
}
xnorm = dnrm2( N - 1, X, strideX, offsetX );
beta = -1.0 * sign( dlapy2( alpha, xnorm ), alpha );
}
tau = ( beta - alpha ) / beta;
dscal( N-1, 1.0 / ( alpha - beta ), X, strideX, offsetX );
for ( i = 0; i < knt; i++ ) {
beta *= safemin;
}

out[ offsetOut ] = beta;
out[ strideOut + offsetOut ] = tau;
}
}


// EXPORTS //

module.exports = dlarfg;
Loading