Skip to content

Commit 8104aed

Browse files
committed
feat: add base algorithm
1 parent 5768926 commit 8104aed

File tree

1 file changed

+126
-0
lines changed
  • lib/node_modules/@stdlib/lapack/base/dppequ/lib

1 file changed

+126
-0
lines changed
Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2024 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+
// MODULES //
22+
23+
var sqrt = require( '@stdlib/math/base/special/sqrt' );
24+
var max = require( '@stdlib/math/base/special/max' );
25+
var min = require( '@stdlib/math/base/special/min' );
26+
27+
28+
// MAIN //
29+
30+
/**
31+
* Computes the row and column scalings intended to equilibrate a symmetric positive definite matrix `A` in packed storage and reduce it's condition number (with respect to the two-norm).
32+
*
33+
* @param {string} order - specifies whether `AP` is packed in row-major or column-major order
34+
* @param {string} uplo - 'Upper' or 'Lower' triangle of `A` is stored
35+
* @param {NonNegativeInteger} N - order of the matrix `A`
36+
* @param {Float64Array} AP - array containing the upper or lower triangle of `A` in packed form
37+
* @param {integer} strideAP - stride length for `AP`
38+
* @param {integer} offsetAP - starting index for `AP`
39+
* @param {Float64Array} S - array to store the scale factors of `A`
40+
* @param {integer} strideS - stride length for `S`
41+
* @param {integer} offsetS - starting index for `S`
42+
* @param {Float64Array} out - array to store the output
43+
* @param {integer} strideOut - stride length for `out`
44+
* @param {integer} offsetOut - starting index for `out`
45+
* @returns {integer} status code
46+
*/
47+
function dppequ( order, uplo, N, AP, strideAP, offsetAP, S, strideS, offsetS, out, strideOut, offsetOut ) { // eslint-disable-line max-len, max-params
48+
var info;
49+
var smin;
50+
var amax;
51+
var jj;
52+
var i;
53+
54+
if ( N === 0 ) {
55+
out[ offsetOut ] = 1.0; // scond
56+
out[ offsetOut + strideOut ] = 0.0; // amax
57+
return 0; // info
58+
}
59+
60+
S[ offsetAP ] = AP[ offsetAP ];
61+
smin = S[ offsetAP ];
62+
amax = S[ offsetAP ];
63+
64+
if ( order === 'row-major' ) {
65+
if ( uplo === 'U' ) { // uplo === 'U'
66+
jj = 0;
67+
for ( i = 1; i < N; i++ ) {
68+
jj += N - i + 1;
69+
S[ offsetS + (i * strideS) ] = AP[ offsetAP + (jj * strideAP) ];
70+
smin = min( smin, S[ offsetS + (i * strideS) ] );
71+
amax = max( amax, S[ offsetS + (i * strideS) ] );
72+
}
73+
} else { // uplo === 'L'
74+
jj = 0;
75+
for ( i = 1; i < N; i++ ) {
76+
jj += i + 1;
77+
S[ offsetS + (i * strideS) ] = AP[ offsetAP + (jj * strideAP) ];
78+
smin = min( smin, S[ offsetS + (i * strideS) ] );
79+
amax = max( amax, S[ offsetS + (i * strideS) ] );
80+
}
81+
}
82+
} else { // order === 'col-major'
83+
if ( uplo === 'U' ) { // uplo === 'U'
84+
jj = 0;
85+
for ( i = 1; i < N; i++ ) {
86+
jj += i + 1;
87+
S[ offsetS + (i * strideS) ] = AP[ offsetAP + (jj * strideAP) ];
88+
smin = min( smin, S[ offsetS + (i * strideS) ] );
89+
amax = max( amax, S[ offsetS + (i * strideS) ] );
90+
}
91+
} else { // uplo === 'L'
92+
jj = 0;
93+
for ( i = 1; i < N; i++ ) {
94+
jj += N - i + 1;
95+
S[ offsetS + (i * strideS) ] = AP[ offsetAP + (jj * strideAP) ];
96+
smin = min( smin, S[ offsetS + (i * strideS) ] );
97+
amax = max( amax, S[ offsetS + (i * strideS) ] );
98+
}
99+
}
100+
}
101+
102+
if ( smin <= 0.0 ) {
103+
for ( i = 0; i < N; i++ ) {
104+
if ( S[ offsetS + (i * strideS) ] <= 0.0 ) {
105+
out[ offsetOut ] = 0.0; // scond
106+
out[ offsetOut + strideOut ] = amax; // amax
107+
info = i;
108+
return info;
109+
}
110+
}
111+
} else {
112+
for ( i = 0; i < N; i++ ) {
113+
S[ offsetS + (i * strideS) ] = 1.0 / sqrt( S[ offsetS + (i * strideS) ] ); // eslint-disable-line max-len
114+
}
115+
}
116+
117+
out[ offsetOut ] = sqrt( smin ) / sqrt( amax ); // scond
118+
out[ offsetOut + strideOut ] = amax; // amax
119+
info = 0;
120+
return info;
121+
}
122+
123+
124+
// EXPORTS //
125+
126+
module.exports = dppequ;

0 commit comments

Comments
 (0)