Skip to content

Commit 3ba563f

Browse files
committed
feat: adds branch 1 logic
1 parent 730fd55 commit 3ba563f

1 file changed

Lines changed: 148 additions & 1 deletion

File tree

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

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

Lines changed: 148 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,4 +14,151 @@
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 reinterpret = require( '@stdlib/strided/base/reinterpret-complex64' );
24+
var isRowMajor = require( '@stdlib/ndarray/base/assert/is-row-major' );
25+
var f32 = require( '@stdlib/number/float64/base/to-float32' );
26+
27+
28+
// MAIN //
29+
30+
/**
31+
* Solves one of the systems of equations `A*x = b` or `A^T*x = b` or `A^H*x = b` where `b` and `x` are `N` element complex vectors and `A` is an `N` by `N` unit, or non-unit, upper or lower triangular complex matrix.
32+
*
33+
* @private
34+
* @param {string} uplo - specifies whether `A` is an upper or lower triangular matrix
35+
* @param {string} trans - specifies whether `A` should be transposed, conjugate-transposed, or not transposed
36+
* @param {string} diag - specifies whether `A` has a unit diagonal
37+
* @param {NonNegativeInteger} N - number of elements along each dimension of `A`
38+
* @param {Complex64Array} A - input matrix
39+
* @param {integer} strideA1 - stride of the first dimension of `A`
40+
* @param {integer} strideA2 - stride of the second dimension of `A`
41+
* @param {NonNegativeInteger} offsetA - starting index for `A`
42+
* @param {Complex64Array} x - input vector
43+
* @param {integer} strideX - `x` stride length
44+
* @param {NonNegativeInteger} offsetX - starting index for `x`
45+
* @returns {Complex64Array} `x`
46+
*
47+
* @example
48+
* var Complex64Array = require( '@stdlib/array/complex64' );
49+
*
50+
* var A = new Complex64Array( [ ] );
51+
* var x = new Complex64Array( [ ] );
52+
*
53+
* ctrsv( 'lower', 'no-transpose', 'non-unit', 3, A, 2, 1, 0, x, 1, 0 );
54+
* // x => <Complex64Array>[ 0.0, 2.0, 0.0, 16.0, 0.0, 46.0 ]
55+
*/
56+
function ctrsv( uplo, trans, diag, N, A, strideA1, strideA2, offsetA, x, strideX, offsetX ) { // eslint-disable-line max-params, max-len
57+
var nonunit;
58+
var viewA;
59+
var viewX;
60+
var retmp;
61+
var imtmp;
62+
var magsq;
63+
var remul;
64+
var immul;
65+
var isrm;
66+
var sign;
67+
var rex;
68+
var imx;
69+
var rea;
70+
var ima;
71+
var ix0;
72+
var ix1;
73+
var sa0;
74+
var sa1;
75+
var oa2;
76+
var ox;
77+
var sx;
78+
var oa;
79+
var i0;
80+
var i1;
81+
var ia;
82+
83+
// Layout
84+
isrm = isRowMajor( [ strideA1, strideA2 ] );
85+
86+
// Diagonal
87+
nonunit = ( diag === 'non-unit' );
88+
89+
// Reinterpret arrays to raw numeric views
90+
viewA = reinterpret( A, 0 );
91+
viewX = reinterpret( x, 0 );
92+
93+
// Set sign to handle conjugation: flip the imaginary part for conjugate-transpose
94+
if ( trans === 'conjugate-transpose' ) {
95+
sign = -1;
96+
} else {
97+
sign = 1;
98+
}
99+
100+
if ( isrm ) {
101+
// For row-major matrices, the last dimension has the fastest changing index...
102+
sa0 = strideA2 * 2; // offset increment for innermost loop
103+
sa1 = strideA1 * 2; // offset increment for outermost loop
104+
} else { // isColMajor
105+
// For column-major matrices, the first dimension has the fastest changing index...
106+
sa0 = strideA1 * 2; // offset increment for innermost loop
107+
sa1 = strideA2 * 2; // offset increment for outermost loop
108+
}
109+
110+
// Vector indexing base
111+
oa = offsetA * 2;
112+
ox = offsetX * 2;
113+
114+
// Vector strides
115+
sx = strideX * 2;
116+
117+
if (
118+
( !isrm && uplo === 'upper' && trans === 'no-transpose' ) ||
119+
( isrm && uplo === 'lower' && trans !== 'no-transpose' )
120+
) {
121+
ix1 = ox + ( ( N - 1 ) * sx );
122+
for ( i1 = N - 1; i1 >= 0; i1-- ) {
123+
rex = viewX[ ix1 ];
124+
imx = viewX[ ix1 + 1 ];
125+
if ( rex !== 0.0 || imx !== 0.0 ) {
126+
oa2 = oa + ( sa1 * i1 );
127+
if ( nonunit ) {
128+
ia = oa2 + ( sa0 * i1 );
129+
rea = viewA[ ia ];
130+
ima = sign * viewA[ ia + 1 ];
131+
magsq = f32( ( rea * rea ) + ( ima * ima ) );
132+
retmp = f32( f32( ( rex * rea ) + ( imx * ima ) ) / magsq );
133+
imtmp = f32( f32( ( imx * rea ) - ( rex * ima ) ) / magsq );
134+
viewX[ ix1 ] = retmp;
135+
viewX[ ix1 + 1 ] = imtmp;
136+
} else {
137+
retmp = rex;
138+
imtmp = imx;
139+
}
140+
ix0 = ix1;
141+
for ( i0 = i1 - 1; i0 >= 0; i0-- ) {
142+
ix0 -= sx;
143+
ia = oa2 + ( sa0 * i0 );
144+
rea = viewA[ ia ];
145+
ima = sign * viewA[ ia + 1 ];
146+
rex = viewX[ ix0 ];
147+
imx = viewX[ ix0 + 1 ];
148+
remul = f32( ( retmp * rea ) - ( imtmp * ima ) );
149+
immul = f32( ( retmp * ima ) + ( imtmp * rea ) );
150+
viewX[ ix0 ] = f32( rex - remul );
151+
viewX[ ix0 + 1 ] = f32( imx - immul );
152+
}
153+
}
154+
ix1 -= sx;
155+
}
156+
return x;
157+
}
158+
159+
}
160+
161+
162+
// EXPORTS //
163+
164+
module.exports = ctrsv;

0 commit comments

Comments
 (0)