1 // SPDX-License-Identifier: Apache-2.0
2 //
3 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
4 // Copyright 2008-2016 National ICT Australia (NICTA)
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 // You may obtain a copy of the License at
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 // ------------------------------------------------------------------------
17
18
19 //! \addtogroup op_norm
20 //! @{
21
22
23
24 template<typename T1>
25 arma_hot
26 inline
27 typename T1::pod_type
vec_norm_1(const Proxy<T1> & P,const typename arma_not_cx<typename T1::elem_type>::result * junk)28 op_norm::vec_norm_1(const Proxy<T1>& P, const typename arma_not_cx<typename T1::elem_type>::result* junk)
29 {
30 arma_extra_debug_sigprint();
31 arma_ignore(junk);
32
33 const bool use_direct_mem = (is_Mat<typename Proxy<T1>::stored_type>::value) || (is_subview_col<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp);
34
35 if(use_direct_mem)
36 {
37 const quasi_unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
38
39 return op_norm::vec_norm_1_direct_std(tmp.M);
40 }
41
42 typedef typename T1::pod_type T;
43
44 T acc = T(0);
45
46 if(Proxy<T1>::use_at == false)
47 {
48 typename Proxy<T1>::ea_type A = P.get_ea();
49
50 const uword N = P.get_n_elem();
51
52 T acc1 = T(0);
53 T acc2 = T(0);
54
55 uword i,j;
56 for(i=0, j=1; j<N; i+=2, j+=2)
57 {
58 acc1 += std::abs(A[i]);
59 acc2 += std::abs(A[j]);
60 }
61
62 if(i < N)
63 {
64 acc1 += std::abs(A[i]);
65 }
66
67 acc = acc1 + acc2;
68 }
69 else
70 {
71 const uword n_rows = P.get_n_rows();
72 const uword n_cols = P.get_n_cols();
73
74 if(n_rows == 1)
75 {
76 for(uword col=0; col<n_cols; ++col)
77 {
78 acc += std::abs(P.at(0,col));
79 }
80 }
81 else
82 {
83 T acc1 = T(0);
84 T acc2 = T(0);
85
86 for(uword col=0; col<n_cols; ++col)
87 {
88 uword i,j;
89
90 for(i=0, j=1; j<n_rows; i+=2, j+=2)
91 {
92 acc1 += std::abs(P.at(i,col));
93 acc2 += std::abs(P.at(j,col));
94 }
95
96 if(i < n_rows)
97 {
98 acc1 += std::abs(P.at(i,col));
99 }
100 }
101
102 acc = acc1 + acc2;
103 }
104 }
105
106 return acc;
107 }
108
109
110
111 template<typename T1>
112 arma_hot
113 inline
114 typename T1::pod_type
vec_norm_1(const Proxy<T1> & P,const typename arma_cx_only<typename T1::elem_type>::result * junk)115 op_norm::vec_norm_1(const Proxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::result* junk)
116 {
117 arma_extra_debug_sigprint();
118 arma_ignore(junk);
119
120 typedef typename T1::elem_type eT;
121 typedef typename T1::pod_type T;
122
123 T acc = T(0);
124
125 if(Proxy<T1>::use_at == false)
126 {
127 typename Proxy<T1>::ea_type A = P.get_ea();
128
129 const uword N = P.get_n_elem();
130
131 for(uword i=0; i<N; ++i)
132 {
133 const std::complex<T>& X = A[i];
134
135 const T a = X.real();
136 const T b = X.imag();
137
138 acc += std::sqrt( (a*a) + (b*b) );
139 }
140 }
141 else
142 {
143 const uword n_rows = P.get_n_rows();
144 const uword n_cols = P.get_n_cols();
145
146 if(n_rows == 1)
147 {
148 for(uword col=0; col<n_cols; ++col)
149 {
150 const std::complex<T>& X = P.at(0,col);
151
152 const T a = X.real();
153 const T b = X.imag();
154
155 acc += std::sqrt( (a*a) + (b*b) );
156 }
157 }
158 else
159 {
160 for(uword col=0; col<n_cols; ++col)
161 for(uword row=0; row<n_rows; ++row)
162 {
163 const std::complex<T>& X = P.at(row,col);
164
165 const T a = X.real();
166 const T b = X.imag();
167
168 acc += std::sqrt( (a*a) + (b*b) );
169 }
170 }
171 }
172
173 if( (acc != T(0)) && arma_isfinite(acc) )
174 {
175 return acc;
176 }
177 else
178 {
179 arma_extra_debug_print("op_norm::vec_norm_1(): detected possible underflow or overflow");
180
181 const quasi_unwrap<typename Proxy<T1>::stored_type> R(P.Q);
182
183 const uword N = R.M.n_elem;
184 const eT* R_mem = R.M.memptr();
185
186 T max_val = priv::most_neg<T>();
187
188 for(uword i=0; i<N; ++i)
189 {
190 const std::complex<T>& X = R_mem[i];
191
192 const T a = std::abs(X.real());
193 const T b = std::abs(X.imag());
194
195 if(a > max_val) { max_val = a; }
196 if(b > max_val) { max_val = b; }
197 }
198
199 if(max_val == T(0)) { return T(0); }
200
201 T alt_acc = T(0);
202
203 for(uword i=0; i<N; ++i)
204 {
205 const std::complex<T>& X = R_mem[i];
206
207 const T a = X.real() / max_val;
208 const T b = X.imag() / max_val;
209
210 alt_acc += std::sqrt( (a*a) + (b*b) );
211 }
212
213 return ( alt_acc * max_val );
214 }
215 }
216
217
218
219 template<typename eT>
220 arma_hot
221 inline
222 eT
vec_norm_1_direct_std(const Mat<eT> & X)223 op_norm::vec_norm_1_direct_std(const Mat<eT>& X)
224 {
225 arma_extra_debug_sigprint();
226
227 const uword N = X.n_elem;
228 const eT* A = X.memptr();
229
230 if(N < uword(32))
231 {
232 return op_norm::vec_norm_1_direct_mem(N,A);
233 }
234 else
235 {
236 #if defined(ARMA_USE_ATLAS)
237 {
238 return atlas::cblas_asum(N,A);
239 }
240 #elif defined(ARMA_USE_BLAS)
241 {
242 return blas::asum(N,A);
243 }
244 #else
245 {
246 return op_norm::vec_norm_1_direct_mem(N,A);
247 }
248 #endif
249 }
250 }
251
252
253
254 template<typename eT>
255 arma_hot
256 inline
257 eT
vec_norm_1_direct_mem(const uword N,const eT * A)258 op_norm::vec_norm_1_direct_mem(const uword N, const eT* A)
259 {
260 arma_extra_debug_sigprint();
261
262 #if defined(ARMA_SIMPLE_LOOPS) || (defined(__FINITE_MATH_ONLY__) && (__FINITE_MATH_ONLY__ > 0))
263 {
264 eT acc1 = eT(0);
265
266 if(memory::is_aligned(A))
267 {
268 memory::mark_as_aligned(A);
269
270 for(uword i=0; i<N; ++i) { acc1 += std::abs(A[i]); }
271 }
272 else
273 {
274 for(uword i=0; i<N; ++i) { acc1 += std::abs(A[i]); }
275 }
276
277 return acc1;
278 }
279 #else
280 {
281 eT acc1 = eT(0);
282 eT acc2 = eT(0);
283
284 uword j;
285
286 for(j=1; j<N; j+=2)
287 {
288 const eT tmp_i = (*A); A++;
289 const eT tmp_j = (*A); A++;
290
291 acc1 += std::abs(tmp_i);
292 acc2 += std::abs(tmp_j);
293 }
294
295 if((j-1) < N)
296 {
297 acc1 += std::abs(*A);
298 }
299
300 return (acc1 + acc2);
301 }
302 #endif
303 }
304
305
306
307 template<typename T1>
308 arma_hot
309 inline
310 typename T1::pod_type
vec_norm_2(const Proxy<T1> & P,const typename arma_not_cx<typename T1::elem_type>::result * junk)311 op_norm::vec_norm_2(const Proxy<T1>& P, const typename arma_not_cx<typename T1::elem_type>::result* junk)
312 {
313 arma_extra_debug_sigprint();
314 arma_ignore(junk);
315
316 const bool use_direct_mem = (is_Mat<typename Proxy<T1>::stored_type>::value) || (is_subview_col<typename Proxy<T1>::stored_type>::value) || (arma_config::openmp && Proxy<T1>::use_mp);
317
318 if(use_direct_mem)
319 {
320 const quasi_unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
321
322 return op_norm::vec_norm_2_direct_std(tmp.M);
323 }
324
325 typedef typename T1::pod_type T;
326
327 T acc = T(0);
328
329 if(Proxy<T1>::use_at == false)
330 {
331 typename Proxy<T1>::ea_type A = P.get_ea();
332
333 const uword N = P.get_n_elem();
334
335 T acc1 = T(0);
336 T acc2 = T(0);
337
338 uword i,j;
339
340 for(i=0, j=1; j<N; i+=2, j+=2)
341 {
342 const T tmp_i = A[i];
343 const T tmp_j = A[j];
344
345 acc1 += tmp_i * tmp_i;
346 acc2 += tmp_j * tmp_j;
347 }
348
349 if(i < N)
350 {
351 const T tmp_i = A[i];
352
353 acc1 += tmp_i * tmp_i;
354 }
355
356 acc = acc1 + acc2;
357 }
358 else
359 {
360 const uword n_rows = P.get_n_rows();
361 const uword n_cols = P.get_n_cols();
362
363 if(n_rows == 1)
364 {
365 for(uword col=0; col<n_cols; ++col)
366 {
367 const T tmp = P.at(0,col);
368
369 acc += tmp * tmp;
370 }
371 }
372 else
373 {
374 for(uword col=0; col<n_cols; ++col)
375 {
376 uword i,j;
377 for(i=0, j=1; j<n_rows; i+=2, j+=2)
378 {
379 const T tmp_i = P.at(i,col);
380 const T tmp_j = P.at(j,col);
381
382 acc += tmp_i * tmp_i;
383 acc += tmp_j * tmp_j;
384 }
385
386 if(i < n_rows)
387 {
388 const T tmp_i = P.at(i,col);
389
390 acc += tmp_i * tmp_i;
391 }
392 }
393 }
394 }
395
396
397 const T sqrt_acc = std::sqrt(acc);
398
399 if( (sqrt_acc != T(0)) && arma_isfinite(sqrt_acc) )
400 {
401 return sqrt_acc;
402 }
403 else
404 {
405 arma_extra_debug_print("op_norm::vec_norm_2(): detected possible underflow or overflow");
406
407 const quasi_unwrap<typename Proxy<T1>::stored_type> tmp(P.Q);
408
409 return op_norm::vec_norm_2_direct_robust(tmp.M);
410 }
411 }
412
413
414
415 template<typename T1>
416 arma_hot
417 inline
418 typename T1::pod_type
vec_norm_2(const Proxy<T1> & P,const typename arma_cx_only<typename T1::elem_type>::result * junk)419 op_norm::vec_norm_2(const Proxy<T1>& P, const typename arma_cx_only<typename T1::elem_type>::result* junk)
420 {
421 arma_extra_debug_sigprint();
422 arma_ignore(junk);
423
424 typedef typename T1::elem_type eT;
425 typedef typename T1::pod_type T;
426
427 T acc = T(0);
428
429 if(Proxy<T1>::use_at == false)
430 {
431 typename Proxy<T1>::ea_type A = P.get_ea();
432
433 const uword N = P.get_n_elem();
434
435 for(uword i=0; i<N; ++i)
436 {
437 const std::complex<T>& X = A[i];
438
439 const T a = X.real();
440 const T b = X.imag();
441
442 acc += (a*a) + (b*b);
443 }
444 }
445 else
446 {
447 const uword n_rows = P.get_n_rows();
448 const uword n_cols = P.get_n_cols();
449
450 if(n_rows == 1)
451 {
452 for(uword col=0; col<n_cols; ++col)
453 {
454 const std::complex<T>& X = P.at(0,col);
455
456 const T a = X.real();
457 const T b = X.imag();
458
459 acc += (a*a) + (b*b);
460 }
461 }
462 else
463 {
464 for(uword col=0; col<n_cols; ++col)
465 for(uword row=0; row<n_rows; ++row)
466 {
467 const std::complex<T>& X = P.at(row,col);
468
469 const T a = X.real();
470 const T b = X.imag();
471
472 acc += (a*a) + (b*b);
473 }
474 }
475 }
476
477 const T sqrt_acc = std::sqrt(acc);
478
479 if( (sqrt_acc != T(0)) && arma_isfinite(sqrt_acc) )
480 {
481 return sqrt_acc;
482 }
483 else
484 {
485 arma_extra_debug_print("op_norm::vec_norm_2(): detected possible underflow or overflow");
486
487 const quasi_unwrap<typename Proxy<T1>::stored_type> R(P.Q);
488
489 const uword N = R.M.n_elem;
490 const eT* R_mem = R.M.memptr();
491
492 T max_val = priv::most_neg<T>();
493
494 for(uword i=0; i<N; ++i)
495 {
496 const T val_i = std::abs(R_mem[i]);
497
498 if(val_i > max_val) { max_val = val_i; }
499 }
500
501 if(max_val == T(0)) { return T(0); }
502
503 T alt_acc = T(0);
504
505 for(uword i=0; i<N; ++i)
506 {
507 const T val_i = std::abs(R_mem[i]) / max_val;
508
509 alt_acc += val_i * val_i;
510 }
511
512 return ( std::sqrt(alt_acc) * max_val );
513 }
514 }
515
516
517
518 template<typename eT>
519 arma_hot
520 inline
521 eT
vec_norm_2_direct_std(const Mat<eT> & X)522 op_norm::vec_norm_2_direct_std(const Mat<eT>& X)
523 {
524 arma_extra_debug_sigprint();
525
526 const uword N = X.n_elem;
527 const eT* A = X.memptr();
528
529 eT result;
530
531 if(N < uword(32))
532 {
533 result = op_norm::vec_norm_2_direct_mem(N,A);
534 }
535 else
536 {
537 #if defined(ARMA_USE_ATLAS)
538 {
539 result = atlas::cblas_nrm2(N,A);
540 }
541 #elif defined(ARMA_USE_BLAS)
542 {
543 result = blas::nrm2(N,A);
544 }
545 #else
546 {
547 result = op_norm::vec_norm_2_direct_mem(N,A);
548 }
549 #endif
550 }
551
552 if( (result != eT(0)) && arma_isfinite(result) )
553 {
554 return result;
555 }
556 else
557 {
558 arma_extra_debug_print("op_norm::vec_norm_2_direct_std(): detected possible underflow or overflow");
559
560 return op_norm::vec_norm_2_direct_robust(X);
561 }
562 }
563
564
565
566 template<typename eT>
567 arma_hot
568 inline
569 eT
vec_norm_2_direct_mem(const uword N,const eT * A)570 op_norm::vec_norm_2_direct_mem(const uword N, const eT* A)
571 {
572 arma_extra_debug_sigprint();
573
574 eT acc;
575
576 #if defined(ARMA_SIMPLE_LOOPS) || (defined(__FINITE_MATH_ONLY__) && (__FINITE_MATH_ONLY__ > 0))
577 {
578 eT acc1 = eT(0);
579
580 if(memory::is_aligned(A))
581 {
582 memory::mark_as_aligned(A);
583
584 for(uword i=0; i<N; ++i) { const eT tmp_i = A[i]; acc1 += tmp_i * tmp_i; }
585 }
586 else
587 {
588 for(uword i=0; i<N; ++i) { const eT tmp_i = A[i]; acc1 += tmp_i * tmp_i; }
589 }
590
591 acc = acc1;
592 }
593 #else
594 {
595 eT acc1 = eT(0);
596 eT acc2 = eT(0);
597
598 uword j;
599
600 for(j=1; j<N; j+=2)
601 {
602 const eT tmp_i = (*A); A++;
603 const eT tmp_j = (*A); A++;
604
605 acc1 += tmp_i * tmp_i;
606 acc2 += tmp_j * tmp_j;
607 }
608
609 if((j-1) < N)
610 {
611 const eT tmp_i = (*A);
612
613 acc1 += tmp_i * tmp_i;
614 }
615
616 acc = acc1 + acc2;
617 }
618 #endif
619
620 return std::sqrt(acc);
621 }
622
623
624
625 template<typename eT>
626 arma_hot
627 inline
628 eT
vec_norm_2_direct_robust(const Mat<eT> & X)629 op_norm::vec_norm_2_direct_robust(const Mat<eT>& X)
630 {
631 arma_extra_debug_sigprint();
632
633 const uword N = X.n_elem;
634 const eT* A = X.memptr();
635
636 eT max_val = priv::most_neg<eT>();
637
638 uword j;
639
640 for(j=1; j<N; j+=2)
641 {
642 eT val_i = (*A); A++;
643 eT val_j = (*A); A++;
644
645 val_i = std::abs(val_i);
646 val_j = std::abs(val_j);
647
648 if(val_i > max_val) { max_val = val_i; }
649 if(val_j > max_val) { max_val = val_j; }
650 }
651
652 if((j-1) < N)
653 {
654 const eT val_i = std::abs(*A);
655
656 if(val_i > max_val) { max_val = val_i; }
657 }
658
659 if(max_val == eT(0)) { return eT(0); }
660
661 const eT* B = X.memptr();
662
663 eT acc1 = eT(0);
664 eT acc2 = eT(0);
665
666 for(j=1; j<N; j+=2)
667 {
668 eT val_i = (*B); B++;
669 eT val_j = (*B); B++;
670
671 val_i /= max_val;
672 val_j /= max_val;
673
674 acc1 += val_i * val_i;
675 acc2 += val_j * val_j;
676 }
677
678 if((j-1) < N)
679 {
680 const eT val_i = (*B) / max_val;
681
682 acc1 += val_i * val_i;
683 }
684
685 return ( std::sqrt(acc1 + acc2) * max_val );
686 }
687
688
689
690 template<typename T1>
691 arma_hot
692 inline
693 typename T1::pod_type
vec_norm_k(const Proxy<T1> & P,const int k)694 op_norm::vec_norm_k(const Proxy<T1>& P, const int k)
695 {
696 arma_extra_debug_sigprint();
697
698 typedef typename T1::pod_type T;
699
700 T acc = T(0);
701
702 if(Proxy<T1>::use_at == false)
703 {
704 typename Proxy<T1>::ea_type A = P.get_ea();
705
706 const uword N = P.get_n_elem();
707
708 for(uword i=0; i<N; ++i)
709 {
710 acc += std::pow(std::abs(A[i]), k);
711 }
712 }
713 else
714 {
715 const uword n_rows = P.get_n_rows();
716 const uword n_cols = P.get_n_cols();
717
718 if(n_rows != 1)
719 {
720 for(uword col=0; col < n_cols; ++col)
721 for(uword row=0; row < n_rows; ++row)
722 {
723 acc += std::pow(std::abs(P.at(row,col)), k);
724 }
725 }
726 else
727 {
728 for(uword col=0; col < n_cols; ++col)
729 {
730 acc += std::pow(std::abs(P.at(0,col)), k);
731 }
732 }
733 }
734
735 return std::pow(acc, T(1)/T(k));
736 }
737
738
739
740 template<typename T1>
741 arma_hot
742 inline
743 typename T1::pod_type
vec_norm_max(const Proxy<T1> & P)744 op_norm::vec_norm_max(const Proxy<T1>& P)
745 {
746 arma_extra_debug_sigprint();
747
748 typedef typename T1::pod_type T;
749
750 const uword N = P.get_n_elem();
751
752 T max_val = (N != 1) ? priv::most_neg<T>() : std::abs(P[0]);
753
754 if(Proxy<T1>::use_at == false)
755 {
756 typename Proxy<T1>::ea_type A = P.get_ea();
757
758 uword i,j;
759 for(i=0, j=1; j<N; i+=2, j+=2)
760 {
761 const T tmp_i = std::abs(A[i]);
762 const T tmp_j = std::abs(A[j]);
763
764 if(max_val < tmp_i) { max_val = tmp_i; }
765 if(max_val < tmp_j) { max_val = tmp_j; }
766 }
767
768 if(i < N)
769 {
770 const T tmp_i = std::abs(A[i]);
771
772 if(max_val < tmp_i) { max_val = tmp_i; }
773 }
774 }
775 else
776 {
777 const uword n_rows = P.get_n_rows();
778 const uword n_cols = P.get_n_cols();
779
780 if(n_rows != 1)
781 {
782 for(uword col=0; col < n_cols; ++col)
783 for(uword row=0; row < n_rows; ++row)
784 {
785 const T tmp = std::abs(P.at(row,col));
786
787 if(max_val < tmp) { max_val = tmp; }
788 }
789 }
790 else
791 {
792 for(uword col=0; col < n_cols; ++col)
793 {
794 const T tmp = std::abs(P.at(0,col));
795
796 if(max_val < tmp) { max_val = tmp; }
797 }
798 }
799 }
800
801 return max_val;
802 }
803
804
805
806 template<typename T1>
807 arma_hot
808 inline
809 typename T1::pod_type
vec_norm_min(const Proxy<T1> & P)810 op_norm::vec_norm_min(const Proxy<T1>& P)
811 {
812 arma_extra_debug_sigprint();
813
814 typedef typename T1::pod_type T;
815
816 const uword N = P.get_n_elem();
817
818 T min_val = (N != 1) ? priv::most_pos<T>() : std::abs(P[0]);
819
820 if(Proxy<T1>::use_at == false)
821 {
822 typename Proxy<T1>::ea_type A = P.get_ea();
823
824 uword i,j;
825 for(i=0, j=1; j<N; i+=2, j+=2)
826 {
827 const T tmp_i = std::abs(A[i]);
828 const T tmp_j = std::abs(A[j]);
829
830 if(min_val > tmp_i) { min_val = tmp_i; }
831 if(min_val > tmp_j) { min_val = tmp_j; }
832 }
833
834 if(i < N)
835 {
836 const T tmp_i = std::abs(A[i]);
837
838 if(min_val > tmp_i) { min_val = tmp_i; }
839 }
840 }
841 else
842 {
843 const uword n_rows = P.get_n_rows();
844 const uword n_cols = P.get_n_cols();
845
846 if(n_rows != 1)
847 {
848 for(uword col=0; col < n_cols; ++col)
849 for(uword row=0; row < n_rows; ++row)
850 {
851 const T tmp = std::abs(P.at(row,col));
852
853 if(min_val > tmp) { min_val = tmp; }
854 }
855 }
856 else
857 {
858 for(uword col=0; col < n_cols; ++col)
859 {
860 const T tmp = std::abs(P.at(0,col));
861
862 if(min_val > tmp) { min_val = tmp; }
863 }
864 }
865 }
866
867 return min_val;
868 }
869
870
871
872 template<typename eT>
873 inline
874 typename get_pod_type<eT>::result
mat_norm_1(const Mat<eT> & X)875 op_norm::mat_norm_1(const Mat<eT>& X)
876 {
877 arma_extra_debug_sigprint();
878
879 // TODO: this can be sped up with a dedicated implementation
880 return as_scalar( max( sum(abs(X), 0), 1) );
881 }
882
883
884
885 template<typename eT>
886 inline
887 typename get_pod_type<eT>::result
mat_norm_2(const Mat<eT> & X)888 op_norm::mat_norm_2(const Mat<eT>& X)
889 {
890 arma_extra_debug_sigprint();
891
892 typedef typename get_pod_type<eT>::result T;
893
894 if(X.is_finite() == false) { arma_debug_warn_level(1, "norm(): given matrix has non-finite elements"); }
895
896 Col<T> S;
897 svd(S, X);
898
899 return (S.n_elem > 0) ? S[0] : T(0);
900 }
901
902
903
904 template<typename eT>
905 inline
906 typename get_pod_type<eT>::result
mat_norm_inf(const Mat<eT> & X)907 op_norm::mat_norm_inf(const Mat<eT>& X)
908 {
909 arma_extra_debug_sigprint();
910
911 // TODO: this can be sped up with a dedicated implementation
912 return as_scalar( max( sum(abs(X), 1), 0) );
913 }
914
915
916
917 //! @}
918