From 94d1264e11945a33f4fda2d4b2a4fe78680405bb Mon Sep 17 00:00:00 2001 From: Pranavchiku Date: Thu, 1 Aug 2024 20:52:14 +0530 Subject: [PATCH] refactor: breakdown base implementation --- .../@stdlib/lapack/base/dlange/lib/base.js | 304 +++++++++++------- 1 file changed, 189 insertions(+), 115 deletions(-) diff --git a/lib/node_modules/@stdlib/lapack/base/dlange/lib/base.js b/lib/node_modules/@stdlib/lapack/base/dlange/lib/base.js index ffccad57e5a..c4eac7a8c01 100644 --- a/lib/node_modules/@stdlib/lapack/base/dlange/lib/base.js +++ b/lib/node_modules/@stdlib/lapack/base/dlange/lib/base.js @@ -32,6 +32,191 @@ var min = require( '@stdlib/math/base/special/min' ); var isnan = require( '@stdlib/assert/is-nan' ); +// FUNCTIONS // + +/** +* Computes maximum absolute norm for given matrix `A`. +* +* @private +* @param {string} order - storage layout +* @param {NonNegativeInteger} M - number of rows of `A` +* @param {NonNegativeInteger} N - number of columns of `A` +* @param {Float64Array} A - input array +* @param {integer} strideA1 - stride of the first dimension of `A` +* @param {integer} strideA2 - stride of the second dimension of `A` +* @param {PositiveInteger} offsetA - starting index for `A` +* @returns {number} maximum absolute norm +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* +* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* +* var out = maxAbsNorm( 'row-major', 2, 2, A, 2, 1, 0 ); +* // returns 4.0 +*/ +function maxAbsNorm( order, M, N, A, strideA1, strideA2, offsetA ) { + var value; + var tmp; + var sa; + var i; + var j; + + // Find max( abs( A[ i, j ] ) ) + value = 0.0; + if ( order === 'column-major' ) { + for ( j = 0; j < N; j++ ) { + sa = offsetA + ( j * strideA1 ); + for ( i = 0; i < M; i++ ) { + tmp = abs( A[ sa + ( i * strideA2 ) ] ); + if ( value <= tmp || isnan( tmp ) ) { + value = tmp; + } + } + } + return value; + } + // order === 'row-major' + for ( j = 0; j < N; j++ ) { + sa = offsetA + ( j * strideA2 ); + for ( i = 0; i < M; i++ ) { + tmp = abs( A[ sa + ( i * strideA1 ) ] ); + if ( value <= tmp || isnan( tmp ) ) { + value = tmp; + } + } + } + return value; +} + +/** +* Computes one norm for a given matrix `A`. +* +* @private +* @param {NonNegativeInteger} M - number of rows of `A` +* @param {NonNegativeInteger} N - number of columns of `A` +* @param {Float64Array} A - input array +* @param {integer} strideA1 - stride of the first dimension of `A` +* @param {integer} strideA2 - stride of the second dimension of `A` +* @param {PositiveInteger} offsetA - starting index for `A` +* @returns {number} one norm +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* +* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* +* var out = oneNorm( 2, 2, A, 2, 1, 0 ); +* // returns 6.0 +*/ +function oneNorm( M, N, A, strideA1, strideA2, offsetA ) { + var value; + var sum; + var sa; + var i; + var j; + + // Find norm1( A ) + value = 0.0; + for ( j = 0; j < N; j++ ) { + sum = 0.0; + sa = offsetA + ( j * strideA2 ); + for ( i = 0; i < M; i++ ) { + sum += abs( A[ sa + ( i * strideA1 ) ] ); + } + if ( value < sum || isnan( sum ) ) { + value = sum; + } + } + return value; +} + +/** +* Computes infinity norm for a given matrix `A`. +* +* @private +* @param {NonNegativeInteger} M - number of rows of `A` +* @param {NonNegativeInteger} N - number of columns of `A` +* @param {Float64Array} A - input array +* @param {integer} strideA1 - stride of the first dimension of `A` +* @param {integer} strideA2 - stride of the second dimension of `A` +* @param {PositiveInteger} offsetA - starting index for `A` +* @param {Float64Array} work - workspace array +* @param {integer} strideW - stride length for `work` +* @param {PositiveInteger} offsetW - starting index for `work` +* @returns {number} one norm +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* +* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* var work = new Float64Array( 2 ); +* +* var out = infNorm( 2, 2, A, 2, 1, 0, work, 1, 0 ); +* // returns 7.0 +*/ +function infNorm( M, N, A, strideA1, strideA2, offsetA, work, strideW, offsetW ) { + var value; + var tmp; + var sa; + var i; + var j; + + // Find normI( A ) + value = 0.0; + for ( i = 0; i < M; i++ ) { + work[ offsetW + ( i * strideW ) ] = 0.0; + } + for ( j = 0; j < N; j++ ) { + sa = offsetA + ( j * strideA2 ); + for ( i = 0; i < M; i++ ) { + work[ offsetW + ( i * strideW ) ] += abs( A[ sa + ( i * strideA1 ) ] ); + } + } + for ( i = 0; i < M; i++ ) { + tmp = work[ offsetW + ( i * strideW ) ]; + if ( value < tmp || isnan( tmp ) ) { + value = tmp; + } + } + return value; +} + +/** +* Computes Frobenius norm for a given matrix `A`. +* +* @private +* @param {NonNegativeInteger} M - number of rows of `A` +* @param {NonNegativeInteger} N - number of columns of `A` +* @param {Float64Array} A - input array +* @param {integer} strideA1 - stride of the first dimension of `A` +* @param {integer} strideA2 - stride of the second dimension of `A` +* @param {PositiveInteger} offsetA - starting index for `A` +* @returns {number} Frobenius norm +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* +* var A = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* +* var out = frobeniusNorm( 2, 2, A, 2, 1, 0 ); +* // returns 5.477225575051661 +*/ +function frobeniusNorm( M, N, A, strideA1, strideA2, offsetA ) { + var value; + var out; + var j; + + // Find normF( A ) + out = new Float64Array( [ 0.0, 1.0 ] ); + for ( j = 0; j < N; j++ ) { + out = dlassq( M, A, strideA1, offsetA + ( j * strideA2 ), out[ 0 ], out[ 1 ], out, 1, 0 ); + } + value = out[ 0 ] * sqrt( out[ 1 ] ); + return value; +} + + // MAIN // /** @@ -62,138 +247,27 @@ var isnan = require( '@stdlib/assert/is-nan' ); function dlange( norm, M, N, A, strideA1, strideA2, offsetA, work, strideW, offsetW ) { var order; var value; - var sum; - var tmp; - var out; - var sa0; - var sa1; - var sa; - var i; - var j; norm = lowercase( norm ); - if ( isRowMajor( [ strideA1, strideA2 ] ) ) { - // For row-major matrices, the last dimension has the fastest changing index... order = 'row-major'; - sa0 = strideA2; // stride for innermost loop - sa1 = strideA1; // stride for outermost loop } else { // order == 'column-major' - // For column-major matrices, the first dimension has the fastest changing index... order = 'column-major'; - sa0 = strideA1; - sa1 = strideA2; } if ( min( M, N ) === 0.0 ) { value = 0.0; } else if ( norm === 'm' ) { - // Find max( abs( A[ i, j ] ) ) - value = 0.0; - if ( order === 'column-major' ) { - for ( j = 0; j < N; j++ ) { - sa = offsetA + ( j * sa0 ); - for ( i = 0; i < M; i++ ) { - tmp = abs( A[ sa + ( i * sa1 ) ] ); - if ( value <= tmp || isnan( tmp ) ) { - value = tmp; - } - } - } - return value; - } - // order === 'row-major' - for ( j = 0; j < M; j++ ) { - sa = offsetA + ( j * sa1 ); - for ( i = 0; i < N; i++ ) { - tmp = abs( A[ sa + ( i * sa0 ) ] ); - if ( value <= tmp || isnan( tmp ) ) { - value = tmp; - } - } - } - return value; + return maxAbsNorm( order, M, N, A, strideA1, strideA2, offsetA ); } if ( norm === 'o' || norm === '1' ) { - // Find norm1( A ) - value = 0.0; - if ( order === 'column-major' ) { - for ( j = 0; j < N; j++ ) { - sum = 0.0; - sa = offsetA + ( j * sa1 ); - for ( i = 0; i < M; i++ ) { - sum += abs( A[ sa + ( i * sa0 ) ] ); - } - if ( value < sum || isnan( sum ) ) { - value = sum; - } - } - return value; - } - // order === 'row-major' - // For computing one norm, first dimension has the fastest changing index... - sa0 = strideA1; - sa1 = strideA2; - for ( j = 0; j < N; j++ ) { - sum = 0.0; - sa = offsetA + ( j * sa1 ); - for ( i = 0; i < M; i++ ) { - sum += abs( A[ sa + ( i * sa0 ) ] ); - } - if ( value < sum || isnan( sum ) ) { - value = sum; - } - } - return value; + return oneNorm( M, N, A, strideA1, strideA2, offsetA ); } if ( norm === 'i' ) { - // Find normI( A ) - for ( i = 0; i < M; i++ ) { - work[ offsetW + ( i * strideW ) ] = 0.0; - } - if ( order === 'column-major' ) { - for ( j = 0; j < N; j++ ) { - sa = offsetA + ( j * sa1 ); - for ( i = 0; i < M; i++ ) { - work[ offsetW + ( i * strideW ) ] += abs( A[ sa + ( i * sa0 ) ] ); - } - } - } else { // order === 'row-major' - // For computing infinity norm, first dimension has the fastest changing index... - sa0 = strideA1; - sa1 = strideA2; - for ( j = 0; j < N; j++ ) { - sa = offsetA + ( j * sa1 ); - for ( i = 0; i < M; i++ ) { - work[ offsetW + ( i * strideW ) ] += abs( A[ sa + ( i * sa0 ) ] ); - } - } - } - value = 0.0; - for ( i = 0; i < M; i++ ) { - tmp = work[ offsetW + ( i * strideW ) ]; - if ( value < tmp || isnan( tmp ) ) { - value = tmp; - } - } - return value; + return infNorm( M, N, A, strideA1, strideA2, offsetA, work, strideW, offsetW ); } if ( norm === 'f' ) { - // Find normF( A ) - out = new Float64Array( [ 0.0, 1.0 ] ); - if ( order === 'column-major' ) { - for ( j = 0; j < N; j++ ) { - out = dlassq( M, A, sa0, offsetA + ( j * sa1 ), out[ 0 ], out[ 1 ], out, 1, 0 ); - } - value = out[ 0 ] * sqrt( out[ 1 ] ); - return value; - } - // order === 'row-major' - for ( j = 0; j < N; j++ ) { - out = dlassq( M, A, sa1, offsetA + ( j * sa0 ), out[ 0 ], out[ 1 ], out, 1, 0 ); - } - value = out[ 0 ] * sqrt( out[ 1 ] ); - return value; + return frobeniusNorm( M, N, A, strideA1, strideA2, offsetA ); } return value; }