Skip to content

Commit f3b6595

Browse files
DivitJain26DivitJain26
committed
feat: add base implementation
--- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown_pkg_readmes status: na - task: lint_markdown_docs status: na - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 32ebb87 commit f3b6595

1 file changed

Lines changed: 163 additions & 1 deletion

File tree

  • lib/node_modules/@stdlib/blas/base/cgerc/lib

‎lib/node_modules/@stdlib/blas/base/cgerc/lib/base.js‎

Lines changed: 163 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,166 @@
1414
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
1515
* See the License for the specific language governing permissions and
1616
* limitations under the License.
17-
*/
17+
*/
18+
19+
'use strict';
20+
21+
// MODULES //
22+
23+
var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
24+
var realf = require( '@stdlib/complex/float32/real' );
25+
var imagf = require( '@stdlib/complex/float32/imag' );
26+
var f32 = require( '@stdlib/number/float64/base/to-float32' );
27+
var reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
28+
var muladd = require( '@stdlib/complex/float32/base/mul-add' ).assign;
29+
30+
31+
// MAIN //
32+
33+
/**
34+
* Performs the rank 1 operation `A = α*x*y^T + A`, where `α` is a scalar, `x` is an `M` element vector, `y` is an `N` element vector, and `A` is an `M` by `N` matrix.
35+
*
36+
* @private
37+
* @param {NonNegativeInteger} M - number of rows in the matrix `A`
38+
* @param {NonNegativeInteger} N - number of columns in the matrix `A`
39+
* @param {Complex64} alpha - scalar constant
40+
* @param {Complex64Array} x - first input vector
41+
* @param {integer} strideX - `x` stride length
42+
* @param {NonNegativeInteger} offsetX - starting index for `x`
43+
* @param {Complex64Array} y - second input vector
44+
* @param {integer} strideY - `y` stride length
45+
* @param {NonNegativeInteger} offsetY - starting index for `y`
46+
* @param {Complex64Array} A - input matrix
47+
* @param {integer} strideA1 - stride of the first dimension of `A`
48+
* @param {integer} strideA2 - stride of the second dimension of `A`
49+
* @param {NonNegativeInteger} offsetA - starting index for `A`
50+
* @returns {Complex64Array} `A`
51+
*
52+
* @example
53+
* var Complex64Array = require( '@stdlib/array/complex64' );
54+
* var Complex64 = require( '@stdlib/complex/float32/ctor' );
55+
*
56+
* var A = new Complex64Array( [ 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 5.0, 5.0, 6.0, 6.0 ] );
57+
* var x = new Complex64Array( [ 1.0, 1.0, 2.0, 2.0 ] );
58+
* var y = new Complex64Array( [ 3.0, 3.0, 2.0, 2.0, 1.0, 1.0 ] );
59+
* var alpha = new Complex64( 0.5, 0.5 );
60+
*
61+
* cgerc( 2, 3, alpha, x, 1, 0, y, 1, 0, A, 3, 1, 0 );
62+
* // A => <Complex64Array>[ 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 10.0, 10.0, 9.0, 9.0, 8.0, 8.0 ]
63+
*/
64+
function cgerc( M, N, alpha, x, strideX, offsetX, y, strideY, offsetY, A, strideA1, strideA2, offsetA ) { // eslint-disable-line max-params, max-len
65+
var realpha;
66+
var imalpha;
67+
var viewA;
68+
var viewX;
69+
var viewY;
70+
var retmp;
71+
var imtmp;
72+
var sign;
73+
var cimy;
74+
var tmp;
75+
var rex;
76+
var imx;
77+
var rey;
78+
var imy;
79+
var sa1;
80+
var sa2;
81+
var da0;
82+
var da1;
83+
var S0;
84+
var S1;
85+
var oa;
86+
var ox;
87+
var oy;
88+
var sx;
89+
var sy;
90+
var ia;
91+
var ix;
92+
var iy;
93+
var i0;
94+
var i1;
95+
96+
// Note on variable naming convention: S#, da#, ia#, i# where # corresponds to the loop number, with `0` being the innermost loop...
97+
98+
// Decompose scalars into real and imaginary components:
99+
realpha = realf( alpha );
100+
imalpha = imagf( alpha );
101+
if ( realpha === 0.0 && imalpha === 0.0 ) {
102+
return A;
103+
}
104+
oa = offsetA * 2;
105+
sa1 = strideA1 * 2;
106+
sa2 = strideA2 * 2;
107+
ox = offsetX * 2;
108+
oy = offsetY * 2;
109+
sx = strideX * 2;
110+
sy = strideY * 2;
111+
112+
// Extract loop variables for purposes of loop interchange: dimensions and loop offset (pointer) increments...
113+
if ( isRowMajor( [ strideA1, strideA2 ] ) ) {
114+
// For row-major matrices, the last dimension has the fastest changing index...
115+
S0 = N;
116+
S1 = M;
117+
da0 = sa2; // offset increment for innermost loop
118+
da1 = sa1 - (S0*sa2); // offset increment for outermost loop
119+
120+
// Swap the vectors...
121+
tmp = x;
122+
x = y;
123+
y = tmp;
124+
125+
tmp = sx;
126+
sx = sy;
127+
sy = tmp;
128+
129+
tmp = ox;
130+
ox = oy;
131+
oy = tmp;
132+
133+
sign = -1;
134+
} else { // order === 'column-major'
135+
// For column-major matrices, the first dimension has the fastest changing index...
136+
S0 = M;
137+
S1 = N;
138+
da0 = sa1; // offset increment for innermost loop
139+
da1 = sa2 - (S0*sa1); // offset increment for outermost loop
140+
141+
sign = 1;
142+
}
143+
// Reinterpret arrays as real-valued views of interleaved real and imaginary components:
144+
viewA = reinterpret( A, 0 );
145+
viewX = reinterpret( x, 0 );
146+
viewY = reinterpret( y, 0 );
147+
ix = ox;
148+
iy = oy;
149+
ia = oa;
150+
for ( i1 = 0; i1 < S1; i1++ ) {
151+
rey = viewY[ iy ];
152+
imy = viewY[ iy+1 ];
153+
154+
// Check whether we can avoid the inner loop entirely...
155+
if ( rey === 0.0 && imy === 0.0 ) {
156+
ia += da0 * S0;
157+
} else {
158+
cimy = sign * imy;
159+
retmp = f32(f32(realpha*rey) - f32(imalpha*cimy));
160+
imtmp = f32(f32(realpha*cimy) + f32(imalpha*rey));
161+
ix = ox;
162+
for ( i0 = 0; i0 < S0; i0++ ) {
163+
rex = viewX[ ix ];
164+
imx = viewX[ ix+1 ];
165+
muladd( rex, imx, retmp, imtmp, viewA[ ia ], viewA[ ia+1 ], viewA, 1, ia ); // eslint-disable-line max-len
166+
ix += sx;
167+
ia += da0;
168+
}
169+
}
170+
iy += sy;
171+
ia += da1;
172+
}
173+
return A;
174+
}
175+
176+
177+
// EXPORTS //
178+
179+
module.exports = cgerc;

0 commit comments

Comments
 (0)