Skip to content

Commit dca43a3

Browse files
committed
chore: add implementation
1 parent 63756c5 commit dca43a3

File tree

6 files changed

+474
-0
lines changed

6 files changed

+474
-0
lines changed
Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
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+
var discreteUniform = require( '@stdlib/random/base/discrete-uniform' );
22+
var filledarrayBy = require( '@stdlib/array/filled-by' );
23+
var Complex64 = require( '@stdlib/complex/float32/ctor' );
24+
var logEach = require( '@stdlib/console/log-each' );
25+
var cher = require( './../lib' );
26+
27+
function rand() {
28+
return new Complex64( discreteUniform( 0, 10 ), discreteUniform( -5, 5 ) );
29+
}
30+
31+
var x = filledarrayBy( 3, 'complex64', rand );
32+
console.log( x.get( 0 ).toString() );
33+
34+
var A = filledarrayBy( 9, 'complex64', rand );
35+
console.log( A.get( 0 ).toString() );
36+
37+
cher( 'row-major', 'upper', 3, 1.0, x, 1, A, 3 );
38+
39+
// Print the results:
40+
logEach( '(%s)', A );
41+
42+
cher.ndarray( 'upper', 3, 1.0, x, 1, 0, A, 3, 1, 0 );
43+
44+
// Print the results:
45+
logEach( '(%s)', A );
Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2020 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 reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
24+
var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
25+
26+
27+
// MAIN //
28+
29+
/**
30+
* Performs the hermitian rank 1 operation `A = alpha*x*x**H + A`, where `alpha` is a real scalar, `X` is an `N` element vector and `A` is an `N` by `N` hermitian matrix.
31+
*
32+
* @private
33+
* @param {string} uplo - specifies whether the upper or lower triangular part of the symmetric matrix `A` should be referenced
34+
* @param {NonNegativeInteger} N - number of elements along each dimension of `A`
35+
* @param {number} alpha - scalar
36+
* @param {Complex64Array} x - input array
37+
* @param {integer} strideX - `x` stride length
38+
* @param {NonNegativeInteger} offsetX - starting `x` index
39+
* @param {Complex64Array} A - input matrix
40+
* @param {integer} strideA1 - stride of the first dimension of `A`
41+
* @param {integer} strideA2 - stride of the second dimension of `A`
42+
* @param {NonNegativeInteger} offsetA - starting index for `A`
43+
* @returns {Float64Array} `A`
44+
*
45+
* @example
46+
* var Complex64Array = require( '@stdlib/array/complex64' );
47+
*
48+
* var x = new Complex64Array( [ 1.0, 2.0, 3.0, 4.0 ] );
49+
* var y = new Complex64Array( [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ] );
50+
*
51+
* cher( 'lower', x.length, 2.0, x, 1, 0, A, 3, 1, 0 );
52+
* // y => <Complex64Array>[ 11.0, 0.0, 22.0, -4.0, 0.0, 0.0, 51.0, 0.0 ]
53+
*/
54+
function cher( uplo, N, alpha, x, strideX, offsetX, A, strideA1, strideA2, offsetA ) {
55+
var viewX;
56+
var viewA;
57+
var isrm;
58+
var tmp1;
59+
var tmp2;
60+
var ix1;
61+
var ix0;
62+
var re0;
63+
var re1;
64+
var im0;
65+
var im1;
66+
var idx;
67+
var sa0;
68+
var sa1;
69+
var i1;
70+
var i0;
71+
var ix;
72+
var ia;
73+
var re;
74+
var im;
75+
var sx;
76+
77+
viewX = reinterpret( x, 0 );
78+
viewA = reinterpret( A, 0 );
79+
80+
isrm = isRowMajor( [ strideA1, strideA2 ] );
81+
82+
if ( N <= 0 || alpha === 0.0 ){
83+
return A;
84+
}
85+
86+
if ( isrm ) {
87+
// For row-major matrices, the last dimension has the fastest changing index...
88+
sa0 = strideA2 * 2; // stride for innermost loop
89+
sa1 = strideA1 * 2; // stride for outermost loop
90+
} else { // isColMajor
91+
// For column-major matrices, the first dimension has the fastest changing index...
92+
sa0 = strideA1 * 2; // stride for innermost loop
93+
sa1 = strideA2 * 2; // stride for outermost loop
94+
}
95+
96+
ix = offsetX * 2;
97+
ia = offsetA * 2;
98+
sx = strideX * 2;
99+
100+
if (
101+
( isrm && uplo === 'upper' ) ||
102+
( !isrm && uplo === 'lower' )
103+
) {
104+
for ( i1 = 0; i1 < N; i1++ ) {
105+
ix1 = ix + ( i1 * sx );
106+
tmp1 = viewX[ ix1 ];
107+
tmp2 = -viewX[ ix1 + 1 ];
108+
re1 = alpha * tmp1;
109+
im1 = alpha * tmp2;
110+
111+
for ( i0 = 0; i0 <= i1; i0++ ) {
112+
ix0 = ix + ( i0 * sx );
113+
re0 = viewX[ ix0 ];
114+
im0 = viewX[ ix0 + 1 ];
115+
116+
re = ( re0 * re1 ) - ( im0 * im1 );
117+
im = ( im0 * re1 ) + ( re0 * im1 );
118+
119+
idx = ia + ( i0 * sa0 ) + ( i1 * sa1 );
120+
viewA[ idx ] += re;
121+
viewA[ idx + 1 ] += im;
122+
}
123+
}
124+
}
125+
// ( isrm && uplo === 'lower' ) || ( !isrm && uplo === 'upper' )
126+
for ( i1 = 0; i1 < N; i1++ ) {
127+
ix1 = ix + ( i1 * sx );
128+
tmp1 = viewX[ ix1 ];
129+
tmp2 = -viewX[ ix1 + 1 ];
130+
re1 = alpha * tmp1;
131+
im1 = alpha * tmp2;
132+
for ( i0 = i1; i0 < N; i0++) {
133+
ix0 = ix + ( i0 * sx );
134+
re0 = viewX[ ix0 ];
135+
im0 = viewX[ ix0 + 1 ];
136+
137+
re = ( re0 * re1 ) - ( im0 * im1 );
138+
im = ( im0 * re1 ) + ( re0 * im1 );
139+
140+
idx = ia + ( i0 * sa0 ) + ( i1 * sa1 );
141+
viewA[ idx ] += re;
142+
viewA[ idx + 1 ] += im;
143+
}
144+
}
145+
return A;
146+
}
147+
148+
149+
// EXPORTS //
150+
151+
module.exports = cher;
Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
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 max = require( '@stdlib/math/base/special/fast/max' );
24+
var isLayout = require( '@stdlib/blas/base/assert/is-layout' );
25+
var isMatrixTriangle = require( '@stdlib/blas/base/assert/is-matrix-triangle' );
26+
var isColumnMajor = require( '@stdlib/ndarray/base/assert/is-column-major-string' );
27+
var stride2offset = require( '@stdlib/strided/base/stride2offset' );
28+
var format = require( '@stdlib/string/format' );
29+
var base = require( './base.js' );
30+
31+
32+
// MAIN //
33+
34+
/**
35+
* Performs the hermitian rank 1 operation `A = alpha*x*x**H + A`, where `alpha` is a real scalar, `X` is an `N` element vector and `A` is an `N` by `N` hermitian matrix.
36+
*
37+
* @param {string} order - storage layout
38+
* @param {string} uplo - specifies whether the upper or lower triangular part of the symmetric matrix `A` should be referenced
39+
* @param {NonNegativeInteger} N - number of elements along each dimension of `A`
40+
* @param {number} alpha - scalar
41+
* @param {Complex64Array} x - input array
42+
* @param {integer} strideX - `x` stride length
43+
* @param {NonNegativeInteger} offsetX - starting `x` index
44+
* @param {Complex64Array} A - input matrix
45+
* @param {integer} strideA1 - stride of the first dimension of `A`
46+
* @param {integer} strideA2 - stride of the second dimension of `A`
47+
* @param {NonNegativeInteger} offsetA - starting index for `A`
48+
* @returns {Float64Array} `A`
49+
*
50+
* @example
51+
* var Complex64Array = require( '@stdlib/array/complex64' );
52+
*
53+
* var x = new Complex64Array( [ 1.0, 2.0, 3.0, 4.0 ] );
54+
* var y = new Complex64Array( [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ] );
55+
*
56+
* cher( 'lower', x.length, 2.0, x, 1, 0, A, 3, 1, 0 );
57+
* // y => <Complex64Array>[ 11.0, 0.0, 22.0, -4.0, 0.0, 0.0, 51.0, 0.0 ]
58+
*/
59+
function cher( order, uplo, N, alpha, x, strideX, A, LDA ) {
60+
var sa1;
61+
var sa2;
62+
var ox;
63+
64+
if ( !isLayout( order ) ) {
65+
throw new TypeError( format( 'invalid argument. First argument must be a valid order. Value: `%s`.', order ) );
66+
}
67+
if ( !isMatrixTriangle( uplo ) ) {
68+
throw new TypeError( format( 'invalid argument. Second argument must specify whether to reference the lower or upper triangular matrix. Value: `%s`.', uplo ) );
69+
}
70+
if ( N < 0 ) {
71+
throw new RangeError( format( 'invalid argument. Third argument must be a nonnegative integer. Value: `%d`.', N ) );
72+
}
73+
if ( strideX === 0 ) {
74+
throw new RangeError( format( 'invalid argument. Sixth argument must be non-zero. Value: `%d`.', strideX ) );
75+
}
76+
if ( LDA < max( 1, N ) ) {
77+
throw new RangeError( format( 'invalid argument. Eighth argument must be greater than or equal to max(1,%d). Value: `%d`.', N, LDA ) );
78+
}
79+
if ( N === 0 || alpha === 0.0 ) {
80+
return A;
81+
}
82+
if ( isColumnMajor( order ) ) {
83+
sa1 = 1;
84+
sa2 = LDA;
85+
} else { // order === 'row-major'
86+
sa1 = LDA;
87+
sa2 = 1;
88+
}
89+
ox = stride2offset( N, strideX );
90+
return base( uplo, N, alpha, x, strideX, ox, A, sa1, sa2, 0 );
91+
}
92+
93+
94+
// EXPORTS //
95+
96+
module.exports = cher;
Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
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+
/**
22+
* Performs the hermitian rank 1 operation `A = alpha*x*x**H + A`, where `alpha` is a real scalar, `X` is an `N` element vector and `A` is an `N` by `N` hermitian matrix.
23+
*
24+
* @module @stdlib/blas/base/cher
25+
*
26+
* @example
27+
* var Complex64Array = require( '@stdlib/array/complex64' );
28+
*
29+
* var x = new Complex64Array( [ 1.0, 2.0, 3.0, 4.0 ] );
30+
* var y = new Complex64Array( [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ] );
31+
*
32+
* cher( 'lower', x.length, 2.0, x, 1, 0, A, 3, 1, 0 );
33+
* // y => <Complex64Array>[ 11.0, 0.0, 22.0, -4.0, 0.0, 0.0, 51.0, 0.0 ]
34+
*
35+
* @example
36+
* var Complex64Array = require( '@stdlib/array/complex64' );
37+
*
38+
* var x = new Complex64Array( [ 1.0, 2.0, 3.0, 4.0 ] );
39+
* var y = new Complex64Array( [ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0 ] );
40+
*
41+
* cher( 'lower', x.length, 2.0, x, 1, 0, A, 3, 1, 0 );
42+
* // y => <Complex64Array>[ 11.0, 0.0, 22.0, -4.0, 0.0, 0.0, 51.0, 0.0 ]
43+
*/
44+
45+
// MODULES //
46+
47+
var join = require( 'path' ).join;
48+
var tryRequire = require( '@stdlib/utils/try-require' );
49+
var isError = require( '@stdlib/assert/is-error' );
50+
var main = require( './main.js' );
51+
52+
53+
// MAIN //
54+
55+
var cher;
56+
var tmp = tryRequire( join( __dirname, './native.js' ) );
57+
if ( isError( tmp ) ) {
58+
cher = main;
59+
} else {
60+
cher = tmp;
61+
}
62+
63+
64+
// EXPORTS //
65+
66+
module.exports = cher;
67+
68+
// exports: { "ndarray": "cher.ndarray" }
Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
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 setReadOnly = require( '@stdlib/utils/define-nonenumerable-read-only-property' );
24+
var cher = require( './cher.js' );
25+
var ndarray = require( './ndarray.js' );
26+
27+
28+
// MAIN //
29+
30+
setReadOnly( cher, 'ndarray', ndarray );
31+
32+
33+
// EXPORTS //
34+
35+
module.exports = cher;

0 commit comments

Comments
 (0)