Skip to content

Commit

Permalink
refactor: breakdown base implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@Pranavchiku
Pranavchiku committed Aug 1, 2024
1 parent 983a292 commit 94d1264
Showing 1 changed file with 189 additions and 115 deletions.
304 changes: 189 additions & 115 deletions lib/node_modules/@stdlib/lapack/base/dlange/lib/base.js
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 //

/**
Expand Down Expand Up @@ -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;
}
Expand Down

0 comments on commit 94d1264

Please sign in to comment.