@@ -75,6 +75,10 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
7575 CBLAS_INT child ;
7676 CBLAS_INT ix ;
7777 CBLAS_INT iy ;
78+ CBLAS_INT sx ;
79+ CBLAS_INT sy ;
80+ CBLAS_INT ox ;
81+ CBLAS_INT oy ;
7882 CBLAS_INT n ;
7983 CBLAS_INT j ;
8084 CBLAS_INT k ;
@@ -88,10 +92,15 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
8892 }
8993 // For a positive stride, sorting in decreasing order is equivalent to providing a negative stride and sorting in increasing order, and, for a negative stride, sorting in decreasing order is equivalent to providing a positive stride and sorting in increasing order...
9094 if ( order < 0.0 ) {
91- strideX *= -1 ;
92- strideY *= -1 ;
93- offsetX -= (N - 1 ) * strideX ;
94- offsetY -= (N - 1 ) * strideY ;
95+ sx = - strideX ;
96+ sy = - strideY ;
97+ ox = offsetX - (N - 1 ) * sx ;
98+ oy = offsetY - (N - 1 ) * sy ;
99+ } else {
100+ sx = strideX ;
101+ sy = strideY ;
102+ ox = offsetX ;
103+ oy = offsetY ;
95104 }
96105 // Set the initial heap size:
97106 n = N ;
@@ -104,8 +113,8 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
104113 if ( parent > 0 ) {
105114 // We need to build the heap...
106115 parent -= 1 ;
107- tx = X [ offsetX + (parent * strideX ) ];
108- ty = Y [ offsetY + (parent * strideY ) ];
116+ tx = X [ ox + (parent * sx ) ];
117+ ty = Y [ oy + (parent * sy ) ];
109118 } else {
110119 // Reduce the heap size:
111120 n -= 1 ;
@@ -115,14 +124,14 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
115124 return ;
116125 }
117126 // Store the last heap value in a temporary variable in order to make room for the heap root being placed into its sorted position:
118- ix = offsetX + (n * strideX );
127+ ix = ox + (n * sx );
119128 tx = X [ ix ];
120- iy = offsetY + (n * strideY );
129+ iy = oy + (n * sy );
121130 ty = Y [ iy ];
122131
123132 // Move the heap root to its sorted position:
124- X [ ix ] = X [ offsetX ];
125- Y [ iy ] = Y [ offsetY ];
133+ X [ ix ] = X [ ox ];
134+ Y [ iy ] = Y [ oy ];
126135 }
127136 // We need to "sift down", pushing `t` down the heap to in order to replace the parent and satisfy the heap property...
128137
@@ -136,20 +145,20 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
136145 // Find the largest child...
137146 k = child + 1 ;
138147 if ( k < n ) {
139- v1 = X [ offsetX + (k * strideX ) ];
140- v2 = X [ offsetX + (child * strideX ) ];
148+ v1 = X [ ox + (k * sx ) ];
149+ v2 = X [ ox + (child * sx ) ];
141150
142151 // Check if a "right" child exists and is "bigger"...
143152 if ( v1 > v2 || stdlib_base_is_nan ( v1 ) || (v1 == v2 && stdlib_base_is_positive_zero ( v1 ) ) ) {
144153 child += 1 ;
145154 }
146155 }
147156 // Check if the largest child is bigger than `t`...
148- v1 = X [ offsetX + (child * strideX ) ];
157+ v1 = X [ ox + (child * sx ) ];
149158 if ( v1 > tx || stdlib_base_is_nan ( v1 ) || ( v1 == tx && stdlib_base_is_positive_zero ( v1 ) ) ) {
150159 // Insert the larger child value:
151- X [ offsetX + (j * strideX ) ] = v1 ;
152- Y [ offsetY + (j * strideY ) ] = Y [ offsetY + (child * strideY ) ];
160+ X [ ox + (j * sx ) ] = v1 ;
161+ Y [ oy + (j * sy ) ] = Y [ oy + (child * sy ) ];
153162
154163 // Update `j` to point to the child index:
155164 j = child ;
@@ -162,7 +171,7 @@ void API_SUFFIX(stdlib_strided_dsort2hp_ndarray)( const CBLAS_INT N, const doubl
162171 }
163172 }
164173 // Insert `t` into the heap:
165- X [ offsetX + (j * strideX ) ] = tx ;
166- Y [ offsetY + (j * strideY ) ] = ty ;
174+ X [ ox + (j * sx ) ] = tx ;
175+ Y [ oy + (j * sy ) ] = ty ;
167176 }
168177}
0 commit comments