@@ -121,7 +121,7 @@ auxlib::inv_rcond(Mat<eT>& A, typename get_pod_type<eT>::result& out_rcond)
121121 podarray<blas_int> ipiv (A.n_rows );
122122
123123 arma_extra_debug_print (" lapack::lange()"
124- norm_val = lapack::lange<eT>(&norm_id, &n, &n, A.memptr (), &lda, junk.memptr ());
124+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen (A) : lapack::lange<eT>(&norm_id, &n, &n, A.memptr (), &lda, junk.memptr ());
125125
126126 arma_extra_debug_print (" lapack::getrf()"
127127 lapack::getrf (&n, &n, A.memptr (), &lda, ipiv.memptr (), &info);
@@ -333,7 +333,7 @@ auxlib::inv_sympd_rcond(Mat<eT>& A, bool& out_sympd_state, eT& out_rcond)
333333 podarray<T> work (A.n_rows );
334334
335335 arma_extra_debug_print (" lapack::lansy()"
336- norm_val = lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
336+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym (A) : lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
337337
338338 arma_extra_debug_print (" lapack::potrf()"
339339 lapack::potrf (&uplo, &n, A.memptr (), &n, &info);
@@ -399,7 +399,7 @@ auxlib::inv_sympd_rcond(Mat< std::complex<T> >& A, bool& out_sympd_state, T& out
399399 podarray<T> work (A.n_rows );
400400
401401 arma_extra_debug_print (" lapack::lanhe()"
402- norm_val = lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
402+ norm_val = (has_blas_float_bug<T>::value) ? auxlib::norm1_sym (A) : lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
403403
404404 arma_extra_debug_print (" lapack::potrf()"
405405 lapack::potrf (&uplo, &n, A.memptr (), &n, &info);
@@ -4039,7 +4039,7 @@ auxlib::solve_square_rcond(Mat<typename T1::elem_type>& out, typename T1::pod_ty
40394039 podarray<blas_int> ipiv (A.n_rows + 2 ); // +2 for paranoia
40404040
40414041 arma_extra_debug_print (" lapack::lange()"
4042- norm_val = lapack::lange<eT>(&norm_id, &n, &n, A.memptr (), &lda, junk.memptr ());
4042+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen (A) : lapack::lange<eT>(&norm_id, &n, &n, A.memptr (), &lda, junk.memptr ());
40434043
40444044 arma_extra_debug_print (" lapack::getrf()"
40454045 lapack::getrf<eT>(&n, &n, A.memptr (), &n, ipiv.memptr (), &info);
@@ -4368,7 +4368,7 @@ auxlib::solve_sympd_rcond(Mat<typename T1::pod_type>& out, bool& out_sympd_state
43684368 podarray<T> work (A.n_rows );
43694369
43704370 arma_extra_debug_print (" lapack::lansy()"
4371- norm_val = lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
4371+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym (A) : lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
43724372
43734373 arma_extra_debug_print (" lapack::potrf()"
43744374 lapack::potrf<eT>(&uplo, &n, A.memptr (), &n, &info);
@@ -4446,7 +4446,7 @@ auxlib::solve_sympd_rcond(Mat< std::complex<typename T1::pod_type> >& out, bool&
44464446 podarray<T> work (A.n_rows );
44474447
44484448 arma_extra_debug_print (" lapack::lanhe()"
4449- norm_val = lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
4449+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym (A) : lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &n, work.memptr ());
44504450
44514451 arma_extra_debug_print (" lapack::potrf()"
44524452 lapack::potrf<eT>(&uplo, &n, A.memptr (), &n, &info);
@@ -5387,7 +5387,7 @@ auxlib::solve_band_rcond_common(Mat<typename T1::elem_type>& out, typename T1::p
53875387
53885388 arma_debug_assert_blas_size (AB,out);
53895389
5390- char norm_id = ' 1'
5390+ // char norm_id = '1';
53915391 char trans = ' N'
53925392 blas_int n = blas_int (N); // assuming square matrix
53935393 blas_int kl = blas_int (KL);
@@ -5398,11 +5398,14 @@ auxlib::solve_band_rcond_common(Mat<typename T1::elem_type>& out, typename T1::p
53985398 blas_int info = blas_int (0 );
53995399 T norm_val = T (0 );
54005400
5401- podarray<T> junk (1 );
5401+ // podarray<T> junk(1);
54025402 podarray<blas_int> ipiv (N + 2 ); // +2 for paranoia
54035403
5404- arma_extra_debug_print (" lapack::langb()"
5405- norm_val = lapack::langb<eT>(&norm_id, &n, &kl, &ku, AB.memptr (), &ldab, junk.memptr ());
5404+ // // NOTE: lapack::langb() and lapack::gbtrf() use incompatible storage formats for the band matrix
5405+ // arma_extra_debug_print("lapack::langb()");
5406+ // norm_val = lapack::langb<eT>(&norm_id, &n, &kl, &ku, AB.memptr(), &ldab, junk.memptr());
5407+
5408+ norm_val = auxlib::norm1_band (A,KL,KU);
54065409
54075410 arma_extra_debug_print (" lapack::gbtrf()"
54085411 lapack::gbtrf<eT>(&n, &n, &kl, &ku, AB.memptr (), &ldab, ipiv.memptr (), &info);
@@ -6098,7 +6101,7 @@ auxlib::rcond(Mat<eT>& A)
60986101 podarray<blas_int> ipiv ( (std::min)(A.n_rows , A.n_cols ) );
60996102
61006103 arma_extra_debug_print (" lapack::lange()"
6101- norm_val = lapack::lange (&norm_id, &m, &n, A.memptr (), &lda, work.memptr ());
6104+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen (A) : lapack::lange (&norm_id, &m, &n, A.memptr (), &lda, work.memptr ());
61026105
61036106 arma_extra_debug_print (" lapack::getrf()"
61046107 lapack::getrf (&m, &n, A.memptr (), &lda, ipiv.memptr (), &info);
@@ -6148,7 +6151,7 @@ auxlib::rcond(Mat< std::complex<T> >& A)
61486151 podarray<blas_int> ipiv ( (std::min)(A.n_rows , A.n_cols ) );
61496152
61506153 arma_extra_debug_print (" lapack::lange()"
6151- norm_val = lapack::lange (&norm_id, &m, &n, A.memptr (), &lda, junk.memptr ());
6154+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_gen (A) : lapack::lange (&norm_id, &m, &n, A.memptr (), &lda, junk.memptr ());
61526155
61536156 arma_extra_debug_print (" lapack::getrf()"
61546157 lapack::getrf (&m, &n, A.memptr (), &lda, ipiv.memptr (), &info);
@@ -6196,7 +6199,7 @@ auxlib::rcond_sympd(Mat<eT>& A, bool& calc_ok)
61966199 podarray<blas_int> iwork ( A.n_rows );
61976200
61986201 arma_extra_debug_print (" lapack::lansy()"
6199- norm_val = lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &lda, work.memptr ());
6202+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym (A) : lapack::lansy (&norm_id, &uplo, &n, A.memptr (), &lda, work.memptr ());
62006203
62016204 arma_extra_debug_print (" lapack::potrf()"
62026205 lapack::potrf (&uplo, &n, A.memptr (), &lda, &info);
@@ -6257,7 +6260,7 @@ auxlib::rcond_sympd(Mat< std::complex<T> >& A, bool& calc_ok)
62576260 podarray< T> rwork ( A.n_rows );
62586261
62596262 arma_extra_debug_print (" lapack::lanhe()"
6260- norm_val = lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &lda, rwork.memptr ());
6263+ norm_val = (has_blas_float_bug<eT>::value) ? auxlib::norm1_sym (A) : lapack::lanhe (&norm_id, &uplo, &n, A.memptr (), &lda, rwork.memptr ());
62616264
62626265 arma_extra_debug_print (" lapack::potrf()"
62636266 lapack::potrf (&uplo, &n, A.memptr (), &lda, &info);
@@ -6701,6 +6704,105 @@ auxlib::rudimentary_sym_check(const Mat< std::complex<T> >& X)
67016704
67026705
67036706
6707+ template <typename eT>
6708+ inline
6709+ typename get_pod_type<eT>::result
6710+ auxlib::norm1_gen (const Mat<eT>& A)
6711+ {
6712+ arma_extra_debug_sigprint ();
6713+
6714+ typedef typename get_pod_type<eT>::result T;
6715+
6716+ if (A.n_elem == 0 ) { return T (0 ); }
6717+
6718+ const uword n_rows = A.n_rows ;
6719+ const uword n_cols = A.n_cols ;
6720+
6721+ T max_val = T (0 );
6722+
6723+ for (uword c=0 ; c < n_cols; ++c)
6724+ {
6725+ const eT* colmem = A.colptr (c);
6726+ T acc_val = T (0 );
6727+
6728+ for (uword r=0 ; r < n_rows; ++r) { acc_val += std::abs (colmem[r]); }
6729+
6730+ max_val = (acc_val > max_val) ? acc_val : max_val;
6731+ }
6732+
6733+ return max_val;
6734+ }
6735+
6736+
6737+
6738+ template <typename eT>
6739+ inline
6740+ typename get_pod_type<eT>::result
6741+ auxlib::norm1_sym (const Mat<eT>& A)
6742+ {
6743+ arma_extra_debug_sigprint ();
6744+
6745+ typedef typename get_pod_type<eT>::result T;
6746+
6747+ if (A.n_elem == 0 ) { return T (0 ); }
6748+
6749+ const uword N = (std::min)(A.n_rows , A.n_cols );
6750+
6751+ T max_val = T (0 );
6752+
6753+ for (uword col=0 ; col < N; ++col)
6754+ {
6755+ const eT* colmem = A.colptr (col);
6756+ T acc_val = T (0 );
6757+
6758+ for (uword c=0 ; c < col; ++c) { acc_val += std::abs (A.at (col,c)); }
6759+
6760+ for (uword r=col; r < N; ++r) { acc_val += std::abs (colmem[r]); }
6761+
6762+ max_val = (acc_val > max_val) ? acc_val : max_val;
6763+ }
6764+
6765+ return max_val;
6766+ }
6767+
6768+
6769+
6770+ template <typename eT>
6771+ inline
6772+ typename get_pod_type<eT>::result
6773+ auxlib::norm1_band (const Mat<eT>& A, const uword KL, const uword KU)
6774+ {
6775+ arma_extra_debug_sigprint ();
6776+
6777+ typedef typename get_pod_type<eT>::result T;
6778+
6779+ if (A.n_elem == 0 ) { return T (0 ); }
6780+
6781+ const uword n_rows = A.n_rows ;
6782+ const uword n_cols = A.n_cols ;
6783+
6784+ T max_val = T (0 );
6785+
6786+ for (uword c=0 ; c < n_cols; ++c)
6787+ {
6788+ const eT* colmem = A.colptr (c);
6789+ T acc_val = T (0 );
6790+
6791+ // use values only from main diagonal + KU upper diagonals + KL lower diagonals
6792+
6793+ const uword start = ( c > KU ) ? (c - KU) : 0 ;
6794+ const uword end = ((c + KL) < n_rows) ? (c + KL) : (n_rows-1 );
6795+
6796+ for (uword r=start; r <= end; ++r) { acc_val += std::abs (colmem[r]); }
6797+
6798+ max_val = (acc_val > max_val) ? acc_val : max_val;
6799+ }
6800+
6801+ return max_val;
6802+ }
6803+
6804+
6805+
67046806//
67056807
67066808
0 commit comments