Skip to content

Commit 5afed2b

Browse files
committed
feat: add branch 1 logic in base.js
1 parent 550593f commit 5afed2b

File tree

1 file changed

+129
-1
lines changed
  • lib/node_modules/@stdlib/blas/base/ctpmv/lib

1 file changed

+129
-1
lines changed

lib/node_modules/@stdlib/blas/base/ctpmv/lib/base.js

Lines changed: 129 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,132 @@
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-string' );
24+
var f32 = require( '@stdlib/number/float64/base/to-float32' );
25+
var reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
26+
var muladd = require( '@stdlib/complex/float32/base/mul-add' ).assign;
27+
28+
29+
// MAIN //
30+
31+
/**
32+
* Performs one of the matrix-vector operations `x = A*x`, or `x = A**T*x`, or `x = A**H*x`, where `x` is an `N` element vector and `A` is an `N` by `N` unit, or non-unit, upper or lower triangular matrix, supplied in packed form.
33+
*
34+
* @private
35+
* @param {string} order - storage layout
36+
* @param {string} uplo - specifies whether `A` is an upper or lower triangular matrix
37+
* @param {string} trans - specifies whether `A` should be transposed, conjugate-transposed, or not transposed
38+
* @param {string} diag - specifies whether `A` has a unit diagonal
39+
* @param {NonNegativeInteger} N - number of elements along each dimension of `A`
40+
* @param {Complex64Array} AP - packed form of a symmetric matrix `A`
41+
* @param {integer} strideAP - `AP` stride length
42+
* @param {NonNegativeInteger} offsetAP - starting index for `AP`
43+
* @param {Complex64Array} x - input vector
44+
* @param {integer} strideX - `x` stride length
45+
* @param {NonNegativeInteger} offsetX - starting index for `x`
46+
* @returns {Complex64Array} `x`
47+
*
48+
* @example
49+
* var Complex64Array = require( '@stdlib/array/float64' );
50+
*
51+
* var AP = new Complex64Array( [ ] );
52+
* var x = new Complex64Array( [ ] );
53+
*
54+
* ctpmv( 'row-major', 'upper', 'no-transpose', 'unit', 3, AP, 1, 0, x, 1, 0 );
55+
* // x => <Complex64Array>[ ]
56+
*/
57+
function ctpmv( order, uplo, trans, diag, N, AP, strideAP, offsetAP, x, strideX, offsetX ) { // eslint-disable-line max-params, max-len
58+
var nonunit;
59+
var viewAP;
60+
var viewX;
61+
var retmp;
62+
var imtmp;
63+
var isrm;
64+
var sign;
65+
var rea;
66+
var ima;
67+
var rex;
68+
var imx;
69+
var iap;
70+
var ix0;
71+
var ix1;
72+
var sap;
73+
var kk;
74+
var ox;
75+
var sx;
76+
var re;
77+
var im;
78+
var i0;
79+
var i1;
80+
81+
// Layout
82+
isrm = isRowMajor( order );
83+
84+
// Diagonal
85+
nonunit = ( diag === 'non-unit' );
86+
87+
// Reinterpret arrays to raw numeric views
88+
viewAP = reinterpret( AP, 0 );
89+
viewX = reinterpret( x, 0 );
90+
91+
// Adjust sign to account for conjugation
92+
if ( trans === 'conjugate-transpose' ) {
93+
sign = -1;
94+
} else {
95+
sign = 1;
96+
}
97+
98+
// Vector indexing base
99+
ox = offsetX * 2;
100+
kk = offsetAP * 2;
101+
102+
// Vector strides
103+
sx = strideX * 2;
104+
sap = strideAP * 2;
105+
106+
if (
107+
( !isrm && trans === 'no-transpose' && uplo === 'upper' ) ||
108+
( isrm && trans !== 'no-transpose' && uplo === 'lower' )
109+
) {
110+
ix1 = ox;
111+
for ( i1 = 0; i1 < N; i1++ ) {
112+
retmp = viewX[ ix1 ];
113+
imtmp = viewX[ ix1 + 1 ];
114+
if (retmp !== 0.0 || imtmp !== 0.0) {
115+
iap = kk;
116+
ix0 = ox;
117+
for ( i0 = 0; i0 < i1; i0++ ) {
118+
rea = viewAP[ iap ];
119+
ima = sign * viewAP[ iap + 1 ];
120+
rex = viewX[ ix0 ];
121+
imx = viewX[ ix0 + 1 ];
122+
muladd( rea, ima, retmp, imtmp, rex, imx, viewX, 1, ix0 );
123+
iap += sap;
124+
ix0 += sx;
125+
}
126+
if ( nonunit ) {
127+
rea = viewAP[ iap ];
128+
ima = sign * viewAP[ iap + 1 ];
129+
rex = viewX[ ix0 ];
130+
imx = viewX[ ix0 + 1 ];
131+
viewX[ ix0 ] = f32( ( rea * rex ) - ( ima * imx ) );
132+
viewX[ ix0 + 1 ] = f32( ( rea * imx ) + ( ima * rex ) );
133+
}
134+
}
135+
ix1 += sx;
136+
kk += ( i1 + 1 ) * sap;
137+
}
138+
return x;
139+
}
140+
}
141+
142+
143+
// EXPORTS //
144+
145+
module.exports = ctpmv;

0 commit comments

Comments
 (0)