|  | 
|  | 1 | +<!-- | 
|  | 2 | +
 | 
|  | 3 | +@license Apache-2.0 | 
|  | 4 | +
 | 
|  | 5 | +Copyright (c) 2025 The Stdlib Authors. | 
|  | 6 | +
 | 
|  | 7 | +Licensed under the Apache License, Version 2.0 (the "License"); | 
|  | 8 | +you may not use this file except in compliance with the License. | 
|  | 9 | +You may obtain a copy of the License at | 
|  | 10 | +
 | 
|  | 11 | +   http://www.apache.org/licenses/LICENSE-2.0 | 
|  | 12 | +
 | 
|  | 13 | +Unless required by applicable law or agreed to in writing, software | 
|  | 14 | +distributed under the License is distributed on an "AS IS" BASIS, | 
|  | 15 | +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | 
|  | 16 | +See the License for the specific language governing permissions and | 
|  | 17 | +limitations under the License. | 
|  | 18 | +
 | 
|  | 19 | +--> | 
|  | 20 | + | 
|  | 21 | +# dlatrs | 
|  | 22 | + | 
|  | 23 | +> LAPACK routine to solve a triangular system of equations with the scale factor set to prevent overflow. | 
|  | 24 | +
 | 
|  | 25 | +<section class="intro"> | 
|  | 26 | + | 
|  | 27 | +The `dlatrs` routine solves a triangular system of linear equations with a scaling factor to prevent overflow. It can handle either the original system or its transpose: | 
|  | 28 | + | 
|  | 29 | +<!-- <equation class="equation" label="eq:triangular_system_linear_equations" align="center" raw="A X = s B \quad \text{or} \quad A^{T} X = s B" alt="Triangular system of linear equations solved by `dlatrs`."> --> | 
|  | 30 | + | 
|  | 31 | +```math | 
|  | 32 | +A X = s B \quad \text{or} \quad A^{T} X = s B | 
|  | 33 | +``` | 
|  | 34 | + | 
|  | 35 | +<!-- </equation> --> | 
|  | 36 | + | 
|  | 37 | +where: | 
|  | 38 | + | 
|  | 39 | +-   `A` is an upper or lower triangular matrix, | 
|  | 40 | +-   `X` is the solution vector, | 
|  | 41 | +-   `B` is the right-hand side vector, | 
|  | 42 | +-   `s` is a scaling factor ( 0 < `s` < 1 ) chosen to avoid numerical overflow. | 
|  | 43 | + | 
|  | 44 | +The routine aims to compute a solution `x` such that the magnitude of all elements in `X` remains safely below the machine overflow threshold. | 
|  | 45 | + | 
|  | 46 | +If it determines that solving the unscaled system would not lead to overflow, it calls the Level 2 BLAS routine `dtrsv` for a standard triangular solve. Otherwise, `dlatrs` applies additional scaling and iterative techniques to avoid overflow and maintain numerical stability. | 
|  | 47 | + | 
|  | 48 | +In the presence of a singular matrix `A` (i.e., `A( j, j )` = 0 for some `j`), the scaling factor `s` is set to zero, and a non trivial solution to the homogeneous system: | 
|  | 49 | + | 
|  | 50 | +<!-- <equation class="equation" label="eq:homogeneous_system" align="center" raw="A X = 0" alt="Non trivial solution to homogeneous system."> --> | 
|  | 51 | + | 
|  | 52 | +```math | 
|  | 53 | +A X = 0 | 
|  | 54 | +``` | 
|  | 55 | + | 
|  | 56 | +<!-- </equation> --> | 
|  | 57 | + | 
|  | 58 | +is returned instead. | 
|  | 59 | + | 
|  | 60 | +To ensure stability during the solve, `dlatrs` relies on several internal estimates, including: | 
|  | 61 | + | 
|  | 62 | +-   the off diagonal column norms, | 
|  | 63 | +-   the scaling factor `s`, | 
|  | 64 | +-   and an estimate of the maximum element in the intermediate vector `X`. | 
|  | 65 | + | 
|  | 66 | +The algorithm may scale `X` at various stages using a factor `scale`, such that the final computed result satisfies: | 
|  | 67 | + | 
|  | 68 | +<!-- <equation class="equation" label="eq:scale_factor" align="center" raw="x \leftarrow scale \cdot x" alt="Vector `X` being scaled by a scale factor `s`."> --> | 
|  | 69 | + | 
|  | 70 | +```math | 
|  | 71 | +x \leftarrow scale \cdot x | 
|  | 72 | +``` | 
|  | 73 | + | 
|  | 74 | +<!-- </equation> --> | 
|  | 75 | + | 
|  | 76 | +This ensures both correctness and numerical safety when solving the triangular system under limited-precision arithmetic. | 
|  | 77 | + | 
|  | 78 | +</section> | 
|  | 79 | + | 
|  | 80 | +<!-- /.intro --> | 
|  | 81 | + | 
|  | 82 | +<section class="usage"> | 
|  | 83 | + | 
|  | 84 | +## Usage | 
|  | 85 | + | 
|  | 86 | +```javascript | 
|  | 87 | +var dlatrs = require( '@stdlib/lapack/base/dlatrs' ); | 
|  | 88 | +``` | 
|  | 89 | + | 
|  | 90 | +#### dlatrs( order, uplo, trans, diag, normin, N, A, LDA, X, CNORM ) | 
|  | 91 | + | 
|  | 92 | +Solves a triangular system of equations with the scale factor set to prevent overflow. | 
|  | 93 | + | 
|  | 94 | +```javascript | 
|  | 95 | +var Float64Array = require( '@stdlib/array/float64' ); | 
|  | 96 | + | 
|  | 97 | +var A = new Float64Array( [ 2.0, 1.0, -1.0, 0.0, 3.0, 2.0, 0.0, 0.0, 4.0 ] ); // => [ [ 2.0, 1.0, -1.0 ], [ 0.0, 3.0, 2.0 ], [ 0.0, 0.0, 4.0 ] ] | 
|  | 98 | +var X = new Float64Array( [ 5.0, 10.0, 20.0 ] ); | 
|  | 99 | +var CNORM = new Float64Array( 3 ); | 
|  | 100 | + | 
|  | 101 | +var scale = dlatrs( 'row-major', 'upper', 'no-transpose', 'non-unit', 'no', 3, A, 3, X, CNORM ); | 
|  | 102 | +// returns 1.0 | 
|  | 103 | +// X => <Float64Array>[ 5.0, 0.0, 5.0 ] | 
|  | 104 | +// CNORM => <Float64Array>[ 0.0, 1.0, 3.0 ] | 
|  | 105 | +``` | 
|  | 106 | + | 
|  | 107 | +The function has the following parameters: | 
|  | 108 | + | 
|  | 109 | +-   **order**: storage layout. | 
|  | 110 | +-   **uplo**: specifies whether `A` is an upper or lower triangular matrix. | 
|  | 111 | +-   **trans**: specifies whether `A` should be transposed, conjugate-transposed, or not transposed. | 
|  | 112 | +-   **diag**: specifies whether `A` has a unit diagonal. | 
|  | 113 | +-   **normin**: specifies whether `CNORM` has the column norms on entry or not, should be either `yes` or `no`. | 
|  | 114 | +-   **N**: number of elements along each dimension of `A`. | 
|  | 115 | +-   **A**: input matrix stored in linear memory as a [`Float64Array`][mdn-float64array]. | 
|  | 116 | +-   **LDA**: stride of the first dimension of `A` (a.k.a., leading dimension of the matrix `A`). | 
|  | 117 | +-   **X**: input vector [`Float64Array`][mdn-float64array] describing the right hand side `B`, should have `N` elements. | 
|  | 118 | +-   **CNORM**: [`Float64Array`][mdn-float64array] used to store the column norms of `A`, should have `N` elements. | 
|  | 119 | + | 
|  | 120 | +`X` is overwritten by the solution vector on the left hand side. If `normin` is `yes` then `CNORM` is an input parameter and it should contain the off diagonal column norms of `A`, else `CNORM` is overwritten by the off diagonal column norms computed by the routine | 
|  | 121 | + | 
|  | 122 | +Note that indexing is relative to the first index. To introduce an offset, use [`typed array`][mdn-typed-array] views. | 
|  | 123 | + | 
|  | 124 | +<!-- eslint-disable stdlib/capitalized-comments, max-len --> | 
|  | 125 | + | 
|  | 126 | +```javascript | 
|  | 127 | +var Float64Array = require( '@stdlib/array/float64' ); | 
|  | 128 | +var Int32Array = require( '@stdlib/array/int32' ); | 
|  | 129 | +var dlatrs = require( '@stdlib/lapack/base/dlatrs' ); | 
|  | 130 | + | 
|  | 131 | +// Initial arrays... | 
|  | 132 | +var A0 = new Float64Array( [ 9999.0, 2.0, 1.0, -1.0, 0.0, 3.0, 2.0, 0.0, 0.0, 4.0 ] ); // => [ [ 2.0, 1.0, -1.0 ], [ 0.0, 3.0, 2.0 ], [ 0.0, 0.0, 4.0 ] ] | 
|  | 133 | +var X0 = new Float64Array( [ 9999.0, 5.0, 10.0, 20.0 ] ); | 
|  | 134 | +var CNORM0 = new Float64Array( 4 ); | 
|  | 135 | + | 
|  | 136 | +// Create offset views... | 
|  | 137 | +var A = new Float64Array( A0.buffer, A0.BYTES_PER_ELEMENT*1 ); // start at 2nd element | 
|  | 138 | +var X = new Float64Array( X0.buffer, X0.BYTES_PER_ELEMENT*1 ); // start at 2nd element | 
|  | 139 | +var CNORM = new Float64Array( CNORM0.buffer, CNORM0.BYTES_PER_ELEMENT*1 ); // start at 2nd element | 
|  | 140 | + | 
|  | 141 | +var scale = dlatrs( 'row-major', 'upper', 'no-transpose', 'non-unit', 'no', 3, A, 3, X, CNORM ); | 
|  | 142 | +// returns 1.0 | 
|  | 143 | +// X0 => <Float64Array>[ 9999.0, 5.0, 0.0, 5.0 ] | 
|  | 144 | +// CNORM0 => <Float64Array>[ 0.0, 0.0, 1.0, 3.0 ] | 
|  | 145 | +``` | 
|  | 146 | + | 
|  | 147 | +<!-- lint disable maximum-heading-length --> | 
|  | 148 | + | 
|  | 149 | +#### dlatrs.ndarray( uplo, trans, diag, normin, N, A, strideA1, strideA2, offsetA, X, strideX, offsetX, CNORM, strideCNORM, offsetCNORM ) | 
|  | 150 | + | 
|  | 151 | +Solves a triangular system of equations with the scale factor set to prevent overflow using alternative indexing semantics. | 
|  | 152 | + | 
|  | 153 | +```javascript | 
|  | 154 | +var Float64Array = require( '@stdlib/array/float64' ); | 
|  | 155 | + | 
|  | 156 | +var A = new Float64Array( [ 2.0, 1.0, -1.0, 0.0, 3.0, 2.0, 0.0, 0.0, 4.0 ] ); // => [ [ 2.0, 1.0, -1.0 ], [ 0.0, 3.0, 2.0 ], [ 0.0, 0.0, 4.0 ] ] | 
|  | 157 | +var X = new Float64Array( [ 5.0, 10.0, 20,0 ] ); | 
|  | 158 | +var CNORM = new Float64Array( 3 ); | 
|  | 159 | + | 
|  | 160 | +var scale = dlatrs( 'upper', 'no-transpose', 'non-unit', 'no', 3, A, 3, 1, 0, X, 1, 0, CNORM, 1, 0 ); | 
|  | 161 | +// returns 1.0 | 
|  | 162 | +// X => <Float64Array>[ 5.0, 0.0, 5.0 ] | 
|  | 163 | +// CNORM => <Float64Array>[ 0.0, 1.0, 3.0 ] | 
|  | 164 | +``` | 
|  | 165 | + | 
|  | 166 | +The function has the following additional parameters: | 
|  | 167 | + | 
|  | 168 | +-   **strideA1**: stride of the first dimensiosn of `A`. | 
|  | 169 | +-   **strideA2**: stride of the second dimension of `A`. | 
|  | 170 | +-   **offsetA**: starting index for `A`. | 
|  | 171 | +-   **strideX**: stride length for `X`. | 
|  | 172 | +-   **offsetX**:  starting index for `X`. | 
|  | 173 | +-   **strideCNORM**: stride length for `CNORM`. | 
|  | 174 | +-   **offsetCNORM**:  starting index for `CNORM`. | 
|  | 175 | + | 
|  | 176 | +While [`typed array`][mdn-typed-array] views mandate a view offset based on the underlying buffer, the offset parameters support indexing semantics based on starting indices. For example, | 
|  | 177 | + | 
|  | 178 | +<!-- eslint-disable max-len --> | 
|  | 179 | + | 
|  | 180 | +```javascript | 
|  | 181 | +var Float64Array = require( '@stdlib/array/float64' ); | 
|  | 182 | + | 
|  | 183 | +var A = new Float64Array( [ 9999.0, 2.0, 1.0, -1.0, 0.0, 3.0, 2.0, 0.0, 0.0, 4.0 ] ); // => [ [ 2.0, 1.0, -1.0 ], [ 0.0, 3.0, 2.0 ], [ 0.0, 0.0, 4.0 ] ] | 
|  | 184 | +var X = new Float64Array( [ 9999.0, 5.0, 10.0, 20,0 ] ); | 
|  | 185 | +var CNORM = new Float64Array( 4 ); | 
|  | 186 | + | 
|  | 187 | +var scale = dlatrs( 'upper', 'no-transpose', 'non-unit', 'no', 3, A, 3, 1, 1, X, 1, 1, CNORM, 1, 1 ); | 
|  | 188 | +// returns 1.0 | 
|  | 189 | +// X => <Float64Array>[ 9999.0, 5.0, 0.0, 5.0 ] | 
|  | 190 | +// CNORM => <Float64Array>[ 0.0, 0.0, 1.0, 3.0 ] | 
|  | 191 | +``` | 
|  | 192 | + | 
|  | 193 | +</section> | 
|  | 194 | + | 
|  | 195 | +<!-- /.usage --> | 
|  | 196 | + | 
|  | 197 | +<section class="notes"> | 
|  | 198 | + | 
|  | 199 | +## Notes | 
|  | 200 | + | 
|  | 201 | +-   Both functions mutate the input arrays `X` and `CNORM` (if `normin` = `no`). | 
|  | 202 | + | 
|  | 203 | +-   `dlatrs()` corresponds to the [LAPACK][LAPACK] routine [`dlatrs`][lapack-dlatrs]. | 
|  | 204 | + | 
|  | 205 | +</section> | 
|  | 206 | + | 
|  | 207 | +<!-- /.notes --> | 
|  | 208 | + | 
|  | 209 | +<section class="examples"> | 
|  | 210 | + | 
|  | 211 | +## Examples | 
|  | 212 | + | 
|  | 213 | +<!-- eslint no-undef: "error" --> | 
|  | 214 | + | 
|  | 215 | +```javascript | 
|  | 216 | +var Float64Array = require( '@stdlib/array/float64' ); | 
|  | 217 | +var numel = require( '@stdlib/ndarray/base/numel' ); | 
|  | 218 | +var uniform = require( '@stdlib/random/array/uniform' ); | 
|  | 219 | +var ndarray2array = require( '@stdlib/ndarray/base/to-array' ); | 
|  | 220 | +var dlatrs = require( '@stdlib/lapack/base/dlatrs' ); | 
|  | 221 | + | 
|  | 222 | +// Specify matrix meta data: | 
|  | 223 | +var shape = [ 3, 3 ]; | 
|  | 224 | +var strides = [ 3, 1 ]; | 
|  | 225 | +var offset = 0; | 
|  | 226 | +var N = numel( shape ); | 
|  | 227 | +var order = 'row-major'; | 
|  | 228 | + | 
|  | 229 | +// Create a matrix stored in linear memory: | 
|  | 230 | +var A = uniform( N, 2.0, 10.0, { | 
|  | 231 | +    'dtype': 'float64' | 
|  | 232 | +}); | 
|  | 233 | + | 
|  | 234 | +console.log( ndarray2array( A, shape, strides, offset, order ) ); | 
|  | 235 | + | 
|  | 236 | +var X = uniform( shape[ 0 ], 0.0, 10.0, { | 
|  | 237 | +    'dtype': 'float64' | 
|  | 238 | +}); | 
|  | 239 | + | 
|  | 240 | +var CNORM = new Float64Array( shape[ 0 ] ); | 
|  | 241 | + | 
|  | 242 | +var scale = dlatrs( order, 'upper', 'no-transpose', 'non-unit', 'no', shape[ 0 ], A, strides[ 0 ], X, CNORM ); | 
|  | 243 | +console.log( ndarray2array( A, shape, strides, offset, order ) ); | 
|  | 244 | +console.log( 'scaling factor: ', scale ); | 
|  | 245 | +``` | 
|  | 246 | + | 
|  | 247 | +</section> | 
|  | 248 | + | 
|  | 249 | +<!-- /.examples --> | 
|  | 250 | + | 
|  | 251 | +<!-- C interface documentation. --> | 
|  | 252 | + | 
|  | 253 | +* * * | 
|  | 254 | + | 
|  | 255 | +<section class="c"> | 
|  | 256 | + | 
|  | 257 | +## C APIs | 
|  | 258 | + | 
|  | 259 | +<!-- Section to include introductory text. Make sure to keep an empty line after the intro `section` element and another before the `/section` close. --> | 
|  | 260 | + | 
|  | 261 | +<section class="intro"> | 
|  | 262 | + | 
|  | 263 | +</section> | 
|  | 264 | + | 
|  | 265 | +<!-- /.intro --> | 
|  | 266 | + | 
|  | 267 | +<!-- C usage documentation. --> | 
|  | 268 | + | 
|  | 269 | +<section class="usage"> | 
|  | 270 | + | 
|  | 271 | +### Usage | 
|  | 272 | + | 
|  | 273 | +```c | 
|  | 274 | +TODO | 
|  | 275 | +``` | 
|  | 276 | + | 
|  | 277 | +#### TODO | 
|  | 278 | + | 
|  | 279 | +TODO. | 
|  | 280 | + | 
|  | 281 | +```c | 
|  | 282 | +TODO | 
|  | 283 | +``` | 
|  | 284 | + | 
|  | 285 | +TODO | 
|  | 286 | + | 
|  | 287 | +```c | 
|  | 288 | +TODO | 
|  | 289 | +``` | 
|  | 290 | + | 
|  | 291 | +</section> | 
|  | 292 | + | 
|  | 293 | +<!-- /.usage --> | 
|  | 294 | + | 
|  | 295 | +<!-- C API usage notes. Make sure to keep an empty line after the `section` element and another before the `/section` close. --> | 
|  | 296 | + | 
|  | 297 | +<section class="notes"> | 
|  | 298 | + | 
|  | 299 | +</section> | 
|  | 300 | + | 
|  | 301 | +<!-- /.notes --> | 
|  | 302 | + | 
|  | 303 | +<!-- C API usage examples. --> | 
|  | 304 | + | 
|  | 305 | +<section class="examples"> | 
|  | 306 | + | 
|  | 307 | +### Examples | 
|  | 308 | + | 
|  | 309 | +```c | 
|  | 310 | +TODO | 
|  | 311 | +``` | 
|  | 312 | + | 
|  | 313 | +</section> | 
|  | 314 | + | 
|  | 315 | +<!-- /.examples --> | 
|  | 316 | + | 
|  | 317 | +</section> | 
|  | 318 | + | 
|  | 319 | +<!-- /.c --> | 
|  | 320 | + | 
|  | 321 | +<!-- Section for related `stdlib` packages. Do not manually edit this section, as it is automatically populated. --> | 
|  | 322 | + | 
|  | 323 | +<section class="related"> | 
|  | 324 | + | 
|  | 325 | +</section> | 
|  | 326 | + | 
|  | 327 | +<!-- /.related --> | 
|  | 328 | + | 
|  | 329 | +<!-- Section for all links. Make sure to keep an empty line after the `section` element and another before the `/section` close. --> | 
|  | 330 | + | 
|  | 331 | +<section class="links"> | 
|  | 332 | + | 
|  | 333 | +[lapack]: https://www.netlib.org/lapack/explore-html/ | 
|  | 334 | + | 
|  | 335 | +[lapack-dlatrs]: https://www.netlib.org/lapack/explore-html-3.6.1/d8/dc3/dlatrs_8f_aab36439db8d657622a1d296b89d4acc0.html | 
|  | 336 | + | 
|  | 337 | +[mdn-float64array]: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Float64Array | 
|  | 338 | + | 
|  | 339 | +[mdn-typed-array]: https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/TypedArray | 
|  | 340 | + | 
|  | 341 | +</section> | 
|  | 342 | + | 
|  | 343 | +<!-- /.links --> | 
0 commit comments