@@ -39,20 +39,20 @@ var abs = require( '@stdlib/math/base/special/abs' );
3939* @private
4040* @param {NonNegativeInteger } N - order of matrix A.
4141* @param {Float64Array } DL - sub diagonal elements of A.
42- * @param {integer } sdl - stride length for `DL`
43- * @param {NonNegativeInteger } odl - starting index of `DL`
42+ * @param {integer } strideDL - stride length for `DL`
43+ * @param {NonNegativeInteger } offsetDL - starting index of `DL`
4444* @param {Float64Array } D - diagonal elements of A.
45- * @param {integer } sd - stride length for `D`
46- * @param {NonNegativeInteger } od - starting index of `D`
45+ * @param {integer } strideD - stride length for `D`
46+ * @param {NonNegativeInteger } offsetD - starting index of `D`
4747* @param {Float64Array } DU - super diagonal elements of A.
48- * @param {integer } sdu - stride length for `DU`
49- * @param {NonNegativeInteger } odu - starting index of `DU`
48+ * @param {integer } strideDU - stride length for `DU`
49+ * @param {NonNegativeInteger } offsetDU - starting index of `DU`
5050* @param {Float64Array } DU2 - vector to store the second super diagonal of `U`
51- * @param {integer } sdu2 - stride length for `DU2`
52- * @param {NonNegativeInteger } odu2 - starting index of `DU2`
51+ * @param {integer } strideDU2 - stride length for `DU2`
52+ * @param {NonNegativeInteger } offsetDU2 - starting index of `DU2`
5353* @param {Int32Array } IPIV - vector of pivot indices
54- * @param {integer } si - stride length for `IPIV`
55- * @param {NonNegativeInteger } oi - starting index of `IPIV`
54+ * @param {integer } strideIPIV - stride length for `IPIV`
55+ * @param {NonNegativeInteger } offsetIPIV - starting index of `IPIV`
5656* @returns {integer } status code
5757*
5858* @example
@@ -72,7 +72,7 @@ var abs = require( '@stdlib/math/base/special/abs' );
7272* // DU2 => <Float64Array>[ 0 ]
7373* // IPIV => <Int32Array>[ 0, 1, 2 ]
7474*/
75- function dgttrf ( N , DL , sdl , odl , D , sd , od , DU , sdu , odu , DU2 , sdu2 , odu2 , IPIV , si , oi ) { // eslint-disable-line max-len, max-params
75+ function dgttrf ( N , DL , strideDL , offsetDL , D , strideD , offsetD , DU , strideDU , offsetDU , DU2 , strideDU2 , offsetDU2 , IPIV , strideIPIV , offsetIPIV ) { // eslint-disable-line max-len, max-params
7676 var fact ;
7777 var temp ;
7878 var idu2 ;
@@ -91,76 +91,86 @@ function dgttrf( N, DL, sdl, odl, D, sd, od, DU, sdu, odu, DU2, sdu2, odu2, IPIV
9191 return 0 ;
9292 }
9393
94- idl = odl ;
95- id = od ;
96- idu = odu ;
97- idu2 = odu2 ;
98- ip = oi ;
94+ idl = offsetDL ;
95+ id = offsetD ;
96+ idu = offsetDU ;
97+ idu2 = offsetDU2 ;
98+ ip = offsetIPIV ;
9999
100100 // Initialise ith element of IPIV as i
101101 for ( i = 0 ; i < N ; i ++ ) {
102- IPIV [ ip + ( si * i ) ] = i ;
103- }
102+ IPIV [ ip ] = i ;
104103
105- // Initialise ith element of DU2 as 0
106- for ( i = 0 ; i < N - 2 ; i ++ ) {
107- DU2 [ idu + ( sdu * i ) ] = 0 ;
104+ if ( i < N - 2 ) {
105+ // Initialise ith element of DU2 as 0
106+ DU2 [ idu2 ] = 0 ;
107+ }
108+
109+ ip += strideIPIV ;
110+ idu2 += strideDU2 ;
108111 }
109112
113+ // Set the pointers to the starting indices
114+ idu2 = offsetDU2 ;
115+ ip = offsetIPIV ;
116+
110117 for ( i = 0 ; i < N - 2 ; i ++ ) {
111- di = sd * i ;
112- dli = sdl * i ;
113- ii = si * i ;
114- dui = sdu * i ;
118+ di = strideD * i ;
119+ dli = strideDL * i ;
120+ ii = strideIPIV * i ;
121+ dui = strideDU * i ;
115122
116123 if ( abs ( D [ id + di ] ) >= abs ( DL [ idl + dli ] ) ) { // No row interchange required, eleminate ith element of DL
117124 if ( D [ id ] !== 0.0 ) {
118125 fact = DL [ idl + dli ] / D [ id + di ] ;
119126 DL [ idl + dli ] = fact ;
120- D [ id + di + sd ] = D [ id + di + sd ] - ( fact * DU [ idu + dui ] ) ; // eslint-disable-line max-len
127+ D [ id + di + strideD ] = D [ id + di + strideD ] - ( fact * DU [ idu + dui ] ) ; // eslint-disable-line max-len
121128 }
122129 } else { // Interchange the ith and (i+1)th rows and eliminate ith element of DL
123130 fact = D [ id + di ] / DL [ idl + dli ] ;
124131 D [ id + di ] = DL [ idl + dli ] ;
125132 DL [ idl + dli ] = fact ;
126133 temp = DU [ idu + dui ] ;
127- DU [ idu + dui ] = D [ id + di + sd ] ;
128- D [ id + di + sd ] = temp - ( fact * D [ id + di + sd ] ) ;
129- DU2 [ idu2 + ( sdu2 * i ) ] = DU [ idu + dui + sdu ] ;
130- DU [ idu + dui + sdu ] = - fact * DU [ idu + dui + sdu ] ;
134+ DU [ idu + dui ] = D [ id + di + strideD ] ;
135+ D [ id + di + strideD ] = temp - ( fact * D [ id + di + strideD ] ) ;
136+ DU2 [ idu2 + ( strideDU2 * i ) ] = DU [ idu + dui + strideDU ] ;
137+ DU [ idu + dui + strideDU ] = - fact * DU [ idu + dui + strideDU ] ;
131138 IPIV [ ip + ii ] = i + 1 ;
132139 }
133140 }
134141
135142 if ( N > 1 ) {
136143 i = N - 2 ;
137- di = sd * i ;
138- dli = sdl * i ;
139- ii = si * i ;
140- dui = sdu * i ;
144+ di = strideD * i ;
145+ dli = strideDL * i ;
146+ ii = strideIPIV * i ;
147+ dui = strideDU * i ;
141148
142149 if ( abs ( D [ id + di ] ) >= abs ( DL [ idl + dli ] ) ) {
143150 if ( D [ id + di ] !== 0 ) {
144151 fact = DL [ idl + dli ] / D [ id + di ] ;
145152 DL [ idl + dli ] = fact ;
146- D [ id + di + sd ] = D [ id + di + sd ] - ( fact * DU [ idu + dui ] ) ;
153+ D [ id + di + strideD ] = D [ id + di + strideD ] - ( fact * DU [ idu + dui ] ) ; // eslint-disable-line max-len
147154 }
148155 } else {
149156 fact = D [ id + di ] / DL [ idl + dli ] ;
150157 D [ id + di ] = DL [ idl + dli ] ;
151158 DL [ idl + dli ] = fact ;
152159 temp = DU [ idu + dui ] ;
153- DU [ idu + dui ] = D [ id + di + sd ] ;
154- D [ id + di + sd ] = temp - ( fact * D [ id + di + sd ] ) ;
160+ DU [ idu + dui ] = D [ id + di + strideD ] ;
161+ D [ id + di + strideD ] = temp - ( fact * D [ id + di + strideD ] ) ;
155162 IPIV [ ip + ii ] = i + 1 ;
156163 }
157164 }
158165
166+ id = offsetD ;
167+
159168 // Check for a 0 on the diagonal of U
160169 for ( i = 0 ; i < N ; i ++ ) {
161- if ( D [ sd * i ] === 0.0 ) {
170+ if ( D [ id ] === 0.0 ) {
162171 return i ;
163172 }
173+ id += strideD ;
164174 }
165175
166176 return 0 ;
0 commit comments