All files / lapack/base/dorm2r/lib base.js

55.24% Statements 79/143
100% Branches 1/1
0% Functions 0/1
55.24% Lines 79/143

Press n or j to go to the next uncovered block, b, p or k for the previous block.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 1443x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x 3x                                                                                                                                 3x 3x 3x 3x 3x  
/**
* @license Apache-2.0
*
* Copyright (c) 2025 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 dlarf1f = require( './dlarf1f.js' );
 
 
// MAIN //
 
/**
* Multiplies a matrix `C` by an orthogonal matrix `Q` from a QR factorization.
*
* ## Notes
*
* `dorm2r` overwrites the general real M by N matrix C with
*
* -   `Q * C`  if side = 'left' and trans = 'no-transpose'
* -   `Q^T * C`  if side = 'left' and trans = 'transpose'
* -   `C * Q`  if side = 'right' and trans = 'no-transpose'
* -   `C * Q^T` if side = 'right' and trans = 'transpose'
*
* where Q is a real orthogonal matrix defined as the product of K elementary reflectors `Q = H(1) H(2) ... H(K)` as returned by `dgeqrf`. Q is of order M if side = 'left' and of order N if side = 'right'.
*
* @private
* @param {string} side - specifies the side of multiplication with `C`
* @param {string} trans - specifies the operation to be performed
* @param {NonNegativeInteger} M - number of rows in matrix `C`
* @param {NonNegativeInteger} N - number of columns in matrix `C`
* @param {NonNegativeInteger} K - number of elementary reflectors whose product defines the matrix `Q`
* @param {Float64Array} A - input matrix containing the elementary reflectors
* @param {integer} strideA1 - stride length for the first dimension of `A`
* @param {integer} strideA2 - stride length for the second dimension of `A`
* @param {NonNegativeInteger} offsetA - starting index for `A`
* @param {Float64Array} tau - array containing the scalar factors of the elementary reflectors
* @param {integer} strideTau - stride length for `tau`
* @param {NonNegativeInteger} offsetTau - starting index for `tau`
* @param {Float64Array} C - input/output matrix to be multiplied
* @param {integer} strideC1 - stride length for the first dimension of `C`
* @param {integer} strideC2 - stride length for the second dimension of `C`
* @param {NonNegativeInteger} offsetC - starting index for `C`
* @param {Float64Array} work - workspace array for intermediate calculations
* @param {integer} strideWork - stride length for `work`
* @param {NonNegativeInteger} offsetWork - starting index for `work`
* @returns {Float64Array} the modified matrix `C`
*
* @example
* var Float64Array = require( '@stdlib/array/float64' );
*
* var A = new Float64Array( [ 1.0, 0.0, 0.0, 2.0, 4.0, 0.0, 3.0, 5.0, 6.0 ] );
* var tau = new Float64Array( [ 7.0, 8.0, 9.0 ] );
* var C = new Float64Array( [ 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0 ] );
* var work = new Float64Array( 3 );
*
* dorm2r( 'left', 'no-transpose', 3, 3, 3, A, 3, 1, 0, tau, 1, 0, C, 3, 1, 0, work, 1, 0 );
* // C => <Float64Array>[ -261638.0, -298618.0, -335598.0, -521066.0, -594715.0, -668364.0, -773933.0, -883324.0, -992715.0 ]
*/
function dorm2r( side, trans, M, N, K, A, strideA1, strideA2, offsetA, tau, strideTau, offsetTau, C, strideC1, strideC2, offsetC, work, strideWork, offsetWork ) { // eslint-disable-line max-params, max-len
	var aOffset;
	var cOffset;
	var tauVal;
	var notran;
	var left;
	var i1;
	var i2;
	var i3;
	var mi;
	var ni;
	var ic;
	var jc;
	var i;

	// Quick return if possible
	if ( M === 0 || N === 0 || K === 0 ) {
		return C;
	}

	left = ( side === 'left' );
	notran = ( trans === 'no-transpose' );

	// Determine loop indices based on side and trans
	if ( ( left && !notran ) || ( !left && notran ) ) {
		i1 = 1;
		i2 = K;
		i3 = 1;
	} else {
		i1 = K;
		i2 = 1;
		i3 = -1;
	}

	// Initialize dimensions
	if ( left ) {
		ni = N;
		jc = 0; // 0-based indexing
	} else {
		mi = M;
		ic = 0; // 0-based indexing
	}

	// Main loop
	for ( i = i1; ( i3 > 0 ) ? ( i <= i2 ) : ( i >= i2 ); i += i3 ) {
		if ( left ) {
			// H(i) is applied to C(i:m,1:n)
			mi = M - i + 1;
			ic = i - 1; // Convert to 0-based indexing
		} else {
			// H(i) is applied to C(1:m,i:n)
			ni = N - i + 1;
			jc = i - 1; // Convert to 0-based indexing
		}

		// Apply H(i) using `dlarf1f`
		tauVal = tau[ offsetTau + ( ( i - 1 ) * strideTau ) ];
		aOffset = offsetA + ( ( i - 1 ) * strideA1 ) + ( ( i - 1 ) * strideA2 );
		cOffset = offsetC + ( ic * strideC1 ) + ( jc * strideC2 );
		dlarf1f( side, mi, ni, A, strideA1, aOffset, tauVal, C, strideC1, strideC2, cOffset, work, strideWork, offsetWork ); // eslint-disable-line max-len
	}

	return C;
}
 
 
// EXPORTS //
 
module.exports = dorm2r;