2323var isnan = require ( '@stdlib/math/base/assert/is-nan' ) ;
2424var dlassq = require ( '@stdlib/lapack/base/dlassq' ) . ndarray ;
2525var isRowMajor = require ( '@stdlib/ndarray/base/assert/is-row-major' ) ;
26- var dasumpw = require ( '@stdlib/blas/ext/base/dasumpw' ) . ndarray ;
26+ var loopOrder = require ( '@stdlib/ndarray/base/nullary-loop-interchange-order' ) ;
27+ var dasum = require ( '@stdlib/blas/base/dasum' ) . ndarray ;
2728var Float64Array = require ( '@stdlib/array/float64' ) ;
2829var min = require ( '@stdlib/math/base/special/min' ) ;
29- var abs = require ( '@stdlib/math/base/special/abs' ) ;
3030var sqrt = require ( '@stdlib/math/base/special/sqrt' ) ;
31+ var abs = require ( '@stdlib/math/base/special/abs' ) ;
3132
3233
33- // MAIN //
34+ // FUNCTIONS //
3435
3536/**
36- * Returns the value of the one norm, or the frobenius norm, or the infinity norm, or the element with the largest absolute value of a real matrix `A`.
37+ * Returns the value of the one norm of a real matrix `A`.
3738*
38- * ## Notes
39+ * @private
40+ * @param {NonNegativeInteger } M - number of rows in `A`
41+ * @param {NonNegativeInteger } N - number of columns in `A`
42+ * @param {Float64Array } A - input array
43+ * @param {integer } strideA1 - stride of the first dimension of `A`
44+ * @param {integer } strideA2 - stride of the second dimension of `A`
45+ * @param {NonNegativeInteger } offsetA - starting index of `A`
46+ * @returns {number } required norm value
3947*
40- * - use `norm` = `max`, to calculate the element with the largest absolute value
41- * - use `norm` = `one`, to calculate the one norm
42- * - use `norm` = `infinity`, to calculate the infinity norm
43- * - use `norm` = `frobenius`, to calculate the frobenius norm
48+ * @example
49+ * var Float64Array = require( '@stdlib/array/float64' );
50+ *
51+ * var A = new Float64Array( [ 1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0 ] );
52+ *
53+ * var out = oneNorm( 3, 4, A, 4, 1, 0 );
54+ * // returns 33.0
55+ */
56+ function oneNorm ( M , N , A , strideA1 , strideA2 , offsetA ) {
57+ var value ;
58+ var ia1 ;
59+ var sum ;
60+ var j ;
61+
62+ value = 0.0 ;
63+ ia1 = offsetA ;
64+ for ( j = 0 ; j < N ; j ++ ) {
65+ sum = dasum ( M , A , strideA1 , ia1 ) ;
66+ if ( value < sum || isnan ( sum ) ) {
67+ value = sum ;
68+ }
69+ ia1 += strideA2 ;
70+ }
71+ return value ;
72+ }
73+
74+ /**
75+ * Returns the absolute value of the maximum element of a real matrix `A`.
76+ *
77+ * @private
78+ * @param {NonNegativeInteger } M - number of rows in `A`
79+ * @param {NonNegativeInteger } N - number of columns in `A`
80+ * @param {Float64Array } A - input array
81+ * @param {integer } strideA1 - stride of the first dimension of `A`
82+ * @param {integer } strideA2 - stride of the second dimension of `A`
83+ * @param {NonNegativeInteger } offsetA - starting index of `A`
84+ * @returns {number } required norm value
85+ *
86+ * @example
87+ * var Float64Array = require( '@stdlib/array/float64' );
88+ *
89+ * var A = new Float64Array( [ 1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0 ] );
90+ *
91+ * var out = maxAbs( 3, 4, A, 4, 1, 0 );
92+ * // returns 12.0
93+ */
94+ function maxAbs ( M , N , A , strideA1 , strideA2 , offsetA ) {
95+ var value ;
96+ var temp ;
97+ var da0 ;
98+ var da1 ;
99+ var ia ;
100+ var sa ;
101+ var sh ;
102+ var S0 ;
103+ var S1 ;
104+ var o ;
105+ var i ;
106+ var j ;
107+
108+ value = 0.0 ;
109+
110+ // Resolve the loop interchange order:
111+ o = loopOrder ( [ M , N ] , [ strideA1 , strideA2 ] ) ;
112+ sh = o . sh ;
113+ sa = o . sx ;
114+ S0 = sh [ 0 ] ;
115+ S1 = sh [ 1 ] ;
116+ da0 = sa [ 0 ] ;
117+ da1 = sa [ 1 ] - ( S0 * sa [ 0 ] ) ;
118+ ia = offsetA ;
119+
120+ for ( i = 0 ; i < S1 ; i ++ ) {
121+ for ( j = 0 ; j < S0 ; j ++ ) {
122+ temp = A [ ia ] ;
123+ if ( value < temp || isnan ( temp ) ) {
124+ value = temp ;
125+ }
126+ ia += da0 ;
127+ }
128+ ia += da1 ;
129+ }
130+ return value ;
131+ }
132+
133+ /**
134+ * Returns the value of the infinity norm of a real matrix `A`.
44135*
45136* @private
46- * @param {string } norm - specifies the type of norm to be calculated
47137* @param {NonNegativeInteger } M - number of rows in `A`
48138* @param {NonNegativeInteger } N - number of columns in `A`
49139* @param {Float64Array } A - input array
@@ -61,123 +151,154 @@ var sqrt = require( '@stdlib/math/base/special/sqrt' );
61151* var A = new Float64Array( [ 1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0 ] );
62152* var work = new Float64Array( 3 );
63153*
64- * var out = dlange( 'frobenius', 3, 4, A, 4, 1, 0, work, 1, 0 );
65- * // returns ~25.5
154+ * var out = infinityNorm( 3, 4, A, 4, 1, 0, work, 1, 0 );
155+ * // returns 30.0
66156*/
67- function dlange ( norm , M , N , A , strideA1 , strideA2 , offsetA , work , strideWork , offsetWork ) { // eslint-disable-line max-len
157+ function infinityNorm ( M , N , A , strideA1 , strideA2 , offsetA , work , strideWork , offsetWork ) { // eslint-disable-line max-len
68158 var value ;
69159 var temp ;
70- var sum ;
71- var out ;
72160 var ia1 ;
73161 var ia2 ;
74162 var iw ;
75163 var i ;
76164 var j ;
77165
78- if ( min ( M , N ) === 0 ) {
79- value = 0.0 ;
166+ iw = offsetWork ;
167+ for ( i = 0 ; i < M ; i ++ ) {
168+ work [ iw ] = 0.0 ;
169+ iw += strideWork ;
80170 }
81171
82- else if ( norm === 'max' ) {
83- value = 0.0 ;
84- if ( isRowMajor ( [ strideA1 , strideA2 ] ) ) {
85- ia1 = offsetA ;
86- for ( i = 0 ; i < M ; i ++ ) {
87- ia2 = 0 ;
88- for ( j = 0 ; j < N ; j ++ ) {
89- temp = A [ ia1 + ia2 ] ;
90- if ( value < temp || isnan ( temp ) ) {
91- value = temp ;
92- }
93- ia2 += strideA2 ;
94- }
95- ia1 += strideA1 ;
96- }
97- } else {
98- ia1 = offsetA ;
172+ if ( isRowMajor ( [ strideA1 , strideA2 ] ) ) {
173+ ia1 = offsetA ;
174+ iw = offsetWork ;
175+ for ( i = 0 ; i < M ; i ++ ) {
176+ ia2 = 0 ;
99177 for ( j = 0 ; j < N ; j ++ ) {
100- ia2 = 0 ;
101- for ( i = 0 ; i < M ; i ++ ) {
102- temp = A [ ia1 + ia2 ] ;
103- if ( value < temp || isnan ( temp ) ) {
104- value = temp ;
105- }
106- ia2 += strideA1 ;
107- }
108- ia1 += strideA2 ;
178+ work [ iw ] += abs ( A [ ia1 + ia2 ] ) ;
179+ ia2 += strideA2 ;
109180 }
181+ ia1 += strideA1 ;
182+ iw += strideWork ;
110183 }
111- }
112-
113- else if ( norm === 'one' ) {
114- value = 0.0 ;
184+ } else {
115185 ia1 = offsetA ;
116186 for ( j = 0 ; j < N ; j ++ ) {
117- sum = dasumpw ( M , A , strideA1 , ia1 ) ;
118- if ( value < sum || isnan ( sum ) ) {
119- value = sum ;
187+ ia2 = 0 ;
188+ iw = offsetWork ;
189+ for ( i = 0 ; i < M ; i ++ ) {
190+ work [ iw ] += abs ( A [ ia1 + ia2 ] ) ;
191+ iw += strideWork ;
192+ ia2 += strideA1 ;
120193 }
121194 ia1 += strideA2 ;
122195 }
123196 }
124197
125- else if ( norm === 'infinity' ) {
126- iw = offsetWork ;
127- for ( i = 0 ; i < M ; i ++ ) {
128- work [ iw ] = 0.0 ;
129- iw += strideWork ;
130- }
198+ value = 0.0 ;
131199
132- if ( isRowMajor ( [ strideA1 , strideA2 ] ) ) {
133- ia1 = offsetA ;
134- iw = offsetWork ;
135- for ( i = 0 ; i < M ; i ++ ) {
136- ia2 = 0 ;
137- for ( j = 0 ; j < N ; j ++ ) {
138- work [ iw ] += abs ( A [ ia1 + ia2 ] ) ;
139- ia2 += strideA2 ;
140- }
141- ia1 += strideA1 ;
142- iw += strideWork ;
143- }
144- } else {
145- ia1 = offsetA ;
146- for ( j = 0 ; j < N ; j ++ ) {
147- ia2 = 0 ;
148- iw = offsetWork ;
149- for ( i = 0 ; i < M ; i ++ ) {
150- work [ iw ] += abs ( A [ ia1 + ia2 ] ) ;
151- iw += strideWork ;
152- ia2 += strideA1 ;
153- }
154- ia1 += strideA2 ;
155- }
200+ iw = offsetWork ;
201+ for ( i = 0 ; i < M ; i ++ ) {
202+ temp = work [ iw ] ;
203+ if ( value < temp || isnan ( temp ) ) {
204+ value = temp ;
156205 }
206+ iw += strideWork ;
207+ }
157208
158- value = 0.0 ;
209+ return value ;
210+ }
159211
160- iw = offsetWork ;
161- for ( i = 0 ; i < M ; i ++ ) {
162- temp = work [ iw ] ;
163- if ( value < temp || isnan ( temp ) ) {
164- value = temp ;
165- }
166- iw += strideWork ;
167- }
212+ /**
213+ * Returns the absolute value of the frobenius norm of a real matrix `A`.
214+ *
215+ * @private
216+ * @param {NonNegativeInteger } M - number of rows in `A`
217+ * @param {NonNegativeInteger } N - number of columns in `A`
218+ * @param {Float64Array } A - input array
219+ * @param {integer } strideA1 - stride of the first dimension of `A`
220+ * @param {integer } strideA2 - stride of the second dimension of `A`
221+ * @param {NonNegativeInteger } offsetA - starting index of `A`
222+ * @returns {number } required norm value
223+ *
224+ * @example
225+ * var Float64Array = require( '@stdlib/array/float64' );
226+ *
227+ * var A = new Float64Array( [ 1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0 ] );
228+ *
229+ * var out = frobeniusNorm( 3, 4, A, 4, 1, 0 );
230+ * // returns ~25.5
231+ */
232+ function frobeniusNorm ( M , N , A , strideA1 , strideA2 , offsetA ) {
233+ var ia1 ;
234+ var out ;
235+ var i ;
236+
237+ ia1 = offsetA ;
238+ out = new Float64Array ( [ 0.0 , 1.0 ] ) ;
239+ for ( i = 0 ; i < N ; i ++ ) {
240+ dlassq ( M , A , strideA1 , ia1 , out [ 0 ] , out [ 1 ] , out , 1 , 0 ) ;
241+ ia1 += strideA2 ;
168242 }
169243
170- else if ( norm === 'frobenius' ) {
171- out = new Float64Array ( [ 0.0 , 1.0 ] ) ;
172- ia1 = offsetA ;
173- for ( i = 0 ; i < N ; i ++ ) {
174- dlassq ( M , A , strideA1 , ia1 , out [ 0 ] , out [ 1 ] , out , 1 , 0 ) ;
175- ia1 += strideA2 ;
176- }
177- value = out [ 0 ] * sqrt ( out [ 1 ] ) ;
244+ return out [ 0 ] * sqrt ( out [ 1 ] ) ;
245+ }
246+
247+
248+ // MAIN //
249+
250+ /**
251+ * Returns the value of the one norm, or the frobenius norm, or the infinity norm, or the element with the largest absolute value of a real matrix `A`.
252+ *
253+ * ## Notes
254+ *
255+ * - use `norm` = `max`, to calculate the element with the largest absolute value
256+ * - use `norm` = `one`, to calculate the one norm
257+ * - use `norm` = `infinity`, to calculate the infinity norm
258+ * - use `norm` = `frobenius`, to calculate the frobenius norm
259+ *
260+ * @private
261+ * @param {string } norm - specifies the type of norm to be calculated
262+ * @param {NonNegativeInteger } M - number of rows in `A`
263+ * @param {NonNegativeInteger } N - number of columns in `A`
264+ * @param {Float64Array } A - input array
265+ * @param {integer } strideA1 - stride of the first dimension of `A`
266+ * @param {integer } strideA2 - stride of the second dimension of `A`
267+ * @param {NonNegativeInteger } offsetA - starting index of `A`
268+ * @param {Float64Array } work - only used to compute the infinity norm, expects `M` indexed elements if computing the infinity norm otherwise it's fine to pass a dummy array
269+ * @param {integer } strideWork - stride length of `work`
270+ * @param {NonNegativeInteger } offsetWork - starting index of `work`
271+ * @returns {number } required norm value
272+ *
273+ * @example
274+ * var Float64Array = require( '@stdlib/array/float64' );
275+ *
276+ * var A = new Float64Array( [ 1.0, 4.0, 7.0, 10.0, 2.0, 5.0, 8.0, 11.0, 3.0, 6.0, 9.0, 12.0 ] );
277+ * var work = new Float64Array( 3 );
278+ *
279+ * var out = dlange( 'frobenius', 3, 4, A, 4, 1, 0, work, 1, 0 );
280+ * // returns ~25.5
281+ */
282+ function dlange ( norm , M , N , A , strideA1 , strideA2 , offsetA , work , strideWork , offsetWork ) { // eslint-disable-line max-len
283+ if ( min ( M , N ) === 0 ) {
284+ return 0.0 ;
178285 }
179286
180- return value ;
287+ if ( norm === 'max' ) {
288+ return maxAbs ( M , N , A , strideA1 , strideA2 , offsetA ) ;
289+ }
290+
291+ if ( norm === 'one' ) {
292+ return oneNorm ( M , N , A , strideA1 , strideA2 , offsetA ) ;
293+ }
294+
295+ if ( norm === 'infinity' ) {
296+ return infinityNorm ( M , N , A , strideA1 , strideA2 , offsetA , work , strideWork , offsetWork ) ; // eslint-disable-line max-len
297+ }
298+
299+ if ( norm === 'frobenius' ) {
300+ return frobeniusNorm ( M , N , A , strideA1 , strideA2 , offsetA ) ;
301+ }
181302}
182303
183304
0 commit comments