|
| 1 | +/** |
| 2 | +* @license Apache-2.0 |
| 3 | +* |
| 4 | +* Copyright (c) 2025 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 | +/* eslint-disable max-len */ |
| 22 | + |
| 23 | +// MODULES // |
| 24 | + |
| 25 | +var iladlc = require( '@stdlib/lapack/base/iladlc' ).ndarray; |
| 26 | +var iladlr = require( '@stdlib/lapack/base/iladlr' ).ndarray; |
| 27 | +var dgemv = require( '@stdlib/blas/base/dgemv' ).ndarray; |
| 28 | +var dger = require( '@stdlib/blas/base/dger' ).ndarray; |
| 29 | +var daxpy = require( '@stdlib/blas/base/daxpy' ).ndarray; |
| 30 | +var dscal = require( '@stdlib/blas/base/dscal' ).ndarray; |
| 31 | + |
| 32 | + |
| 33 | +// MAIN // |
| 34 | + |
| 35 | +/** |
| 36 | +* Applies a real elementary reflector `H = I - tau * v * v ^ T` to a real M by N matrix `C`. |
| 37 | +* |
| 38 | +* ## Notes |
| 39 | +* |
| 40 | +* - `work` should have `N` indexed elements if side = `left` and `M` indexed elements if side = `right`. |
| 41 | +* - `V` should have `1 + (M-1) * abs(strideV)` indexed elements if side = `left` and `1 + (N-1) * abs(strideV)` indexed elements if side = `right`. |
| 42 | +* - `C` is overwritten by `H * C` if side = `left` and `C * H` if side = `right`. |
| 43 | +* |
| 44 | +* @private |
| 45 | +* @param {string} side - specifies the side of multiplication with `C`. Use `left` to form `H * C` and `right` to from `C * H`. |
| 46 | +* @param {NonNegativeInteger} M - number of rows in `C` |
| 47 | +* @param {NonNegativeInteger} N - number of columns in `C` |
| 48 | +* @param {Float64Array} V - the vector `v` in the representation of `H` |
| 49 | +* @param {integer} strideV - stride length for `V` |
| 50 | +* @param {NonNegativeInteger} offsetV - starting index for `V` |
| 51 | +* @param {number} tau - the value of `tau` in representation of `H` |
| 52 | +* @param {Float64Array} C - input matrix |
| 53 | +* @param {integer} strideC1 - stride of the first dimension of `C` |
| 54 | +* @param {integer} strideC2 - stride of the second dimension of `C` |
| 55 | +* @param {NonNegativeInteger} offsetC - starting index for `C` |
| 56 | +* @param {Float64Array} work - workspace array |
| 57 | +* @param {integer} strideWork - stride length for `work` |
| 58 | +* @param {NonNegativeInteger} offsetWork - starting index for `work` |
| 59 | +* @returns {Float64Array} `C * H` or `H * C` |
| 60 | +* |
| 61 | +* @example |
| 62 | +* var Float64Array = require( '@stdlib/array/float64' ); |
| 63 | +* |
| 64 | +* var C = new Float64Array( [ 1.0, 5.0, 9.0, 2.0, 6.0, 10.0, 3.0, 7.0, 11.0, 4.0, 8.0, 12.0 ] ); |
| 65 | +* var V = new Float64Array( [ 0.5, 0.5, 0.5, 0.5 ] ); |
| 66 | +* var work = new Float64Array( 3 ); |
| 67 | +* |
| 68 | +* var out = dlarf1f( 'left', 4, 3, V, 1, 0, 1.0, C, 3, 1, 0, work, 1, 0 ); |
| 69 | +* // returns <Float64Array>[ -4.5, -10.5, -16.5, -0.75, -1.75, -2.75, 0.25, -0.75, -1.75, 1.25, 0.25, -0.75 ] |
| 70 | +*/ |
| 71 | +function dlarf1f( side, M, N, V, strideV, offsetV, tau, C, strideC1, strideC2, offsetC, work, strideWork, offsetWork ) { // eslint-disable-line max-params |
| 72 | + var lastv; |
| 73 | + var lastc; |
| 74 | + var i; |
| 75 | + |
| 76 | + lastv = 1; |
| 77 | + lastc = 0; |
| 78 | + if ( tau !== 0.0 ) { |
| 79 | + if ( side === 'left' ) { |
| 80 | + lastv = M; |
| 81 | + } else { |
| 82 | + lastv = N; |
| 83 | + } |
| 84 | + |
| 85 | + // i points to the last element in V |
| 86 | + i = offsetV + ( ( lastv - 1 ) * strideV ); |
| 87 | + |
| 88 | + // Move i to the last non-zero element in V |
| 89 | + while ( lastv > 0 && V[ i ] === 0.0 ) { |
| 90 | + lastv -= 1; |
| 91 | + i -= strideV; |
| 92 | + } |
| 93 | + if ( side === 'left' ) { |
| 94 | + lastc = iladlc( lastv + 1, N, C, strideC1, strideC2, offsetC ) + 1; // to account for the difference between zero-based and one-based indexing |
| 95 | + } else { |
| 96 | + lastc = iladlr( M, lastv + 1, C, strideC1, strideC2, offsetC ) + 1; // to account for the difference between zero-based and one-based indexing |
| 97 | + } |
| 98 | + // lastc is zero if a matrix is filled with zeros |
| 99 | + } |
| 100 | + if ( lastc === 0 ) { |
| 101 | + // Returns C unchanged if tau is zero or all elements in C are zero |
| 102 | + return C; |
| 103 | + } |
| 104 | + if ( side === 'left' ) { |
| 105 | + if ( lastv === 0 ) { |
| 106 | + dscal( lastc, 1.0 - tau, C, strideC2, offsetC ); // scale the first row |
| 107 | + } else { |
| 108 | + dgemv( 'transpose', lastv-1, lastc, 1.0, C, strideC1, strideC2, offsetC + strideC1, V, strideV, offsetV + strideV, 0.0, work, strideWork, offsetWork ); // C( 1, 0 ) is accessed here |
| 109 | + daxpy( lastc, 1.0, C, strideC2, offsetC, work, strideWork, offsetWork ); // operates on the first row of C |
| 110 | + daxpy( lastc, -tau, work, strideWork, offsetWork, C, strideC2, offsetC ); // operates on the first row of C |
| 111 | + dger( lastv-1, lastc, -tau, V, strideV, offsetV + strideV, work, strideWork, offsetWork, C, strideC1, strideC2, offsetC + strideC1 ); // C( 1, 0 ) is accessed here |
| 112 | + } |
| 113 | + } else if ( lastv === 0 ) { |
| 114 | + dscal( lastc, 1.0 - tau, C, strideC1, offsetC ); // scale the first column |
| 115 | + } else { |
| 116 | + dgemv( 'no-transpose', lastc, lastv-1, 1.0, C, strideC1, strideC2, offsetC + strideC2, V, strideV, offsetV + strideV, 0.0, work, strideWork, offsetWork ); // C( 0, 1 ) is accessed here |
| 117 | + daxpy( lastc, 1.0, C, strideC1, offsetC, work, strideWork, offsetWork ); // operates on the first column of C |
| 118 | + daxpy( lastc, -tau, work, strideWork, offsetWork, C, strideC1, offsetC ); // operates on the first column of C |
| 119 | + dger( lastc, lastv-1, -tau, work, strideWork, offsetWork, V, strideV, offsetV + strideV, C, strideC1, strideC2, offsetC + strideC2 ); // C( 0, 1 ) is accessed here |
| 120 | + } |
| 121 | + |
| 122 | + return C; |
| 123 | +} |
| 124 | + |
| 125 | + |
| 126 | +// EXPORTS // |
| 127 | + |
| 128 | +module.exports = dlarf1f; |
0 commit comments