| Mat_meat.hpp | | Mat_meat.hpp | |
| | | | |
| skipping to change at line 175 | | skipping to change at line 175 | |
| if(new_n_elem <= arma_config::mat_prealloc) | | if(new_n_elem <= arma_config::mat_prealloc) | |
| { | | { | |
| access::rw(mem) = mem_local; | | access::rw(mem) = mem_local; | |
| } | | } | |
| else | | else | |
| { | | { | |
| arma_extra_debug_print("Mat::init(): allocating memory"); | | arma_extra_debug_print("Mat::init(): allocating memory"); | |
| | | | |
| access::rw(mem) = new(std::nothrow) eT[new_n_elem]; | | access::rw(mem) = new(std::nothrow) eT[new_n_elem]; | |
| | | | |
|
| arma_check( (mem == 0), "Mat::init(): out of memory" ); | | arma_check_bad_alloc( (mem == 0), "Mat::init(): out of memory" ); | |
| } | | } | |
| | | | |
| access::rw(n_elem) = new_n_elem; | | access::rw(n_elem) = new_n_elem; | |
| access::rw(n_rows) = in_n_rows; | | access::rw(n_rows) = in_n_rows; | |
| access::rw(n_cols) = in_n_cols; | | access::rw(n_cols) = in_n_cols; | |
| access::rw(mem_state) = 0; | | access::rw(mem_state) = 0; | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| | | | |
| skipping to change at line 2210 | | skipping to change at line 2210 | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| Mat<eT>::insert_rows(const u32 row_num, const Base<eT,T1>& X) | | Mat<eT>::insert_rows(const u32 row_num, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| const Mat<eT>& C = tmp.M; | | const Mat<eT>& C = tmp.M; | |
| | | | |
|
| const u32 N = C.n_rows; | | const u32 C_n_rows = C.n_rows; | |
| | | const u32 C_n_cols = C.n_cols; | |
| | | | |
| const u32 t_n_rows = n_rows; | | const u32 t_n_rows = n_rows; | |
| const u32 t_n_cols = n_cols; | | const u32 t_n_cols = n_cols; | |
| | | | |
| const u32 A_n_rows = row_num; | | const u32 A_n_rows = row_num; | |
| const u32 B_n_rows = t_n_rows - row_num; | | const u32 B_n_rows = t_n_rows - row_num; | |
| | | | |
|
| | | bool err_state = false; | |
| | | char* err_msg = 0; | |
| | | | |
| // insertion at row_num == n_rows is in effect an append operation | | // insertion at row_num == n_rows is in effect an append operation | |
|
| arma_debug_check( (row_num > t_n_rows), "Mat::insert_rows(): out of bou | | | |
| nds"); | | | |
| arma_debug_check( (C.n_cols != t_n_cols), "Mat::insert_rows(): given obje | | | |
| ct has an incompatible number of columns"); | | | |
| | | | |
|
| if(N > 0) | | arma_debug_set_error | |
| | | ( | |
| | | err_state, | |
| | | err_msg, | |
| | | (row_num > t_n_rows), | |
| | | "Mat::insert_rows(): out of bounds" | |
| | | ); | |
| | | | |
| | | arma_debug_set_error | |
| | | ( | |
| | | err_state, | |
| | | err_msg, | |
| | | ( (C_n_cols != t_n_cols) && ( (t_n_rows > 0) || (t_n_cols > 0) ) && ( ( | |
| | | C_n_rows > 0) || (C_n_cols > 0) ) ), | |
| | | "Mat::insert_rows(): given object has an incompatible number of columns | |
| | | " | |
| | | ); | |
| | | | |
| | | arma_debug_check(err_state, err_msg); | |
| | | | |
| | | if(C_n_rows > 0) | |
| { | | { | |
|
| Mat<eT> out(t_n_rows + N, t_n_cols); | | Mat<eT> out( t_n_rows + C_n_rows, (std::max)(t_n_cols, C_n_cols) ); | |
| | | | |
|
| if(A_n_rows > 0) | | if(t_n_cols > 0) | |
| { | | { | |
|
| out.rows(0, A_n_rows-1) = rows(0, A_n_rows-1); | | if(A_n_rows > 0) | |
| | | { | |
| | | out.rows(0, A_n_rows-1) = rows(0, A_n_rows-1); | |
| | | } | |
| | | | |
| | | if( (t_n_cols > 0) && (B_n_rows > 0) ) | |
| | | { | |
| | | out.rows(row_num + C_n_rows, t_n_rows + C_n_rows - 1) = rows(row_nu | |
| | | m, t_n_rows - 1); | |
| | | } | |
| } | | } | |
| | | | |
|
| if(B_n_rows > 0) | | if(C_n_cols > 0) | |
| { | | { | |
|
| out.rows(row_num + N, t_n_rows + N - 1) = rows(row_num, t_n_rows - 1)
; | | out.rows(row_num, row_num + C_n_rows - 1) = C; | |
| } | | } | |
| | | | |
|
| out.rows(row_num, row_num + N - 1) = C; | | | |
| | | | |
| steal_mem(out); | | steal_mem(out); | |
| } | | } | |
| } | | } | |
| | | | |
| //! insert the given object at the specified column position; | | //! insert the given object at the specified column position; | |
| //! the given object must have the same number of rows as the matrix | | //! the given object must have the same number of rows as the matrix | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| Mat<eT>::insert_cols(const u32 col_num, const Base<eT,T1>& X) | | Mat<eT>::insert_cols(const u32 col_num, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| const Mat<eT>& C = tmp.M; | | const Mat<eT>& C = tmp.M; | |
| | | | |
|
| const u32 N = C.n_cols; | | const u32 C_n_rows = C.n_rows; | |
| | | const u32 C_n_cols = C.n_cols; | |
| | | | |
| const u32 t_n_rows = n_rows; | | const u32 t_n_rows = n_rows; | |
| const u32 t_n_cols = n_cols; | | const u32 t_n_cols = n_cols; | |
| | | | |
| const u32 A_n_cols = col_num; | | const u32 A_n_cols = col_num; | |
| const u32 B_n_cols = t_n_cols - col_num; | | const u32 B_n_cols = t_n_cols - col_num; | |
| | | | |
|
| | | bool err_state = false; | |
| | | char* err_msg = 0; | |
| | | | |
| // insertion at col_num == n_cols is in effect an append operation | | // insertion at col_num == n_cols is in effect an append operation | |
|
| arma_debug_check( (col_num > t_n_cols), "Mat::insert_cols(): out of bou | | | |
| nds"); | | | |
| arma_debug_check( (C.n_rows != t_n_rows), "Mat::insert_cols(): given obje | | | |
| ct has an incompatible number of rows"); | | | |
| | | | |
|
| if(N > 0) | | arma_debug_set_error | |
| | | ( | |
| | | err_state, | |
| | | err_msg, | |
| | | (col_num > t_n_cols), | |
| | | "Mat::insert_cols(): out of bounds" | |
| | | ); | |
| | | | |
| | | arma_debug_set_error | |
| | | ( | |
| | | err_state, | |
| | | err_msg, | |
| | | ( (C_n_rows != t_n_rows) && ( (t_n_rows > 0) || (t_n_cols > 0) ) && ( ( | |
| | | C_n_rows > 0) || (C_n_cols > 0) ) ), | |
| | | "Mat::insert_cols(): given object has an incompatible number of rows" | |
| | | ); | |
| | | | |
| | | arma_debug_check(err_state, err_msg); | |
| | | | |
| | | if(C_n_cols > 0) | |
| { | | { | |
|
| Mat<eT> out(t_n_rows, t_n_cols + N); | | Mat<eT> out( (std::max)(t_n_rows, C_n_rows), t_n_cols + C_n_cols ); | |
| | | | |
|
| if(A_n_cols > 0) | | if(t_n_rows > 0) | |
| { | | { | |
|
| out.cols(0, A_n_cols-1) = cols(0, A_n_cols-1); | | if(A_n_cols > 0) | |
| | | { | |
| | | out.cols(0, A_n_cols-1) = cols(0, A_n_cols-1); | |
| | | } | |
| | | | |
| | | if(B_n_cols > 0) | |
| | | { | |
| | | out.cols(col_num + C_n_cols, t_n_cols + C_n_cols - 1) = cols(col_nu | |
| | | m, t_n_cols - 1); | |
| | | } | |
| } | | } | |
| | | | |
|
| if(B_n_cols > 0) | | if(C_n_rows > 0) | |
| { | | { | |
|
| out.cols(col_num + N, t_n_cols + N - 1) = cols(col_num, t_n_cols - 1)
; | | out.cols(col_num, col_num + C_n_cols - 1) = C; | |
| } | | } | |
| | | | |
|
| out.cols(col_num, col_num + N - 1) = C; | | | |
| | | | |
| steal_mem(out); | | steal_mem(out); | |
| } | | } | |
| } | | } | |
| | | | |
| //! create a matrix from Op, i.e. run the previously delayed unary operatio
ns | | //! create a matrix from Op, i.e. run the previously delayed unary operatio
ns | |
| template<typename eT> | | template<typename eT> | |
| template<typename T1, typename op_type> | | template<typename T1, typename op_type> | |
| inline | | inline | |
| Mat<eT>::Mat(const Op<T1, op_type>& X) | | Mat<eT>::Mat(const Op<T1, op_type>& X) | |
| : n_rows(0) | | : n_rows(0) | |
| | | | |
| skipping to change at line 3148 | | skipping to change at line 3200 | |
| //! returns true if the object can be interpreted as a column or row vector | | //! returns true if the object can be interpreted as a column or row vector | |
| template<typename eT> | | template<typename eT> | |
| arma_inline | | arma_inline | |
| arma_warn_unused | | arma_warn_unused | |
| bool | | bool | |
| Mat<eT>::is_vec() const | | Mat<eT>::is_vec() const | |
| { | | { | |
| return ( (n_rows == 1) || (n_cols == 1) ); | | return ( (n_rows == 1) || (n_cols == 1) ); | |
| } | | } | |
| | | | |
|
| | | //! returns true if the object can be interpreted as a row vector | |
| | | template<typename eT> | |
| | | arma_inline | |
| | | arma_warn_unused | |
| | | bool | |
| | | Mat<eT>::is_rowvec() const | |
| | | { | |
| | | return (n_rows == 1); | |
| | | } | |
| | | | |
| | | //! returns true if the object can be interpreted as a column vector | |
| | | template<typename eT> | |
| | | arma_inline | |
| | | arma_warn_unused | |
| | | bool | |
| | | Mat<eT>::is_colvec() const | |
| | | { | |
| | | return (n_cols == 1); | |
| | | } | |
| | | | |
| //! returns true if the object has the same number of non-zero rows and col
umnns | | //! returns true if the object has the same number of non-zero rows and col
umnns | |
| template<typename eT> | | template<typename eT> | |
| arma_inline | | arma_inline | |
| arma_warn_unused | | arma_warn_unused | |
| bool | | bool | |
| Mat<eT>::is_square() const | | Mat<eT>::is_square() const | |
| { | | { | |
| return (n_rows == n_cols); | | return (n_rows == n_cols); | |
| } | | } | |
| | | | |
| | | | |
| skipping to change at line 3888 | | skipping to change at line 3960 | |
| switch(type) | | switch(type) | |
| { | | { | |
| case raw_ascii: | | case raw_ascii: | |
| save_okay = diskio::save_raw_ascii(*this, name); | | save_okay = diskio::save_raw_ascii(*this, name); | |
| break; | | break; | |
| | | | |
| case arma_ascii: | | case arma_ascii: | |
| save_okay = diskio::save_arma_ascii(*this, name); | | save_okay = diskio::save_arma_ascii(*this, name); | |
| break; | | break; | |
| | | | |
|
| | | case csv_ascii: | |
| | | save_okay = diskio::save_csv_ascii(*this, name); | |
| | | break; | |
| | | | |
| case raw_binary: | | case raw_binary: | |
| save_okay = diskio::save_raw_binary(*this, name); | | save_okay = diskio::save_raw_binary(*this, name); | |
| break; | | break; | |
| | | | |
| case arma_binary: | | case arma_binary: | |
| save_okay = diskio::save_arma_binary(*this, name); | | save_okay = diskio::save_arma_binary(*this, name); | |
| break; | | break; | |
| | | | |
| case pgm_binary: | | case pgm_binary: | |
| save_okay = diskio::save_pgm_binary(*this, name); | | save_okay = diskio::save_pgm_binary(*this, name); | |
| | | | |
| skipping to change at line 3930 | | skipping to change at line 4006 | |
| switch(type) | | switch(type) | |
| { | | { | |
| case raw_ascii: | | case raw_ascii: | |
| save_okay = diskio::save_raw_ascii(*this, os); | | save_okay = diskio::save_raw_ascii(*this, os); | |
| break; | | break; | |
| | | | |
| case arma_ascii: | | case arma_ascii: | |
| save_okay = diskio::save_arma_ascii(*this, os); | | save_okay = diskio::save_arma_ascii(*this, os); | |
| break; | | break; | |
| | | | |
|
| | | case csv_ascii: | |
| | | save_okay = diskio::save_csv_ascii(*this, os); | |
| | | break; | |
| | | | |
| case raw_binary: | | case raw_binary: | |
| save_okay = diskio::save_raw_binary(*this, os); | | save_okay = diskio::save_raw_binary(*this, os); | |
| break; | | break; | |
| | | | |
| case arma_binary: | | case arma_binary: | |
| save_okay = diskio::save_arma_binary(*this, os); | | save_okay = diskio::save_arma_binary(*this, os); | |
| break; | | break; | |
| | | | |
| case pgm_binary: | | case pgm_binary: | |
| save_okay = diskio::save_pgm_binary(*this, os); | | save_okay = diskio::save_pgm_binary(*this, os); | |
| | | | |
| skipping to change at line 3977 | | skipping to change at line 4057 | |
| break; | | break; | |
| | | | |
| case raw_ascii: | | case raw_ascii: | |
| load_okay = diskio::load_raw_ascii(*this, name, err_msg); | | load_okay = diskio::load_raw_ascii(*this, name, err_msg); | |
| break; | | break; | |
| | | | |
| case arma_ascii: | | case arma_ascii: | |
| load_okay = diskio::load_arma_ascii(*this, name, err_msg); | | load_okay = diskio::load_arma_ascii(*this, name, err_msg); | |
| break; | | break; | |
| | | | |
|
| | | case csv_ascii: | |
| | | load_okay = diskio::load_csv_ascii(*this, name, err_msg); | |
| | | break; | |
| | | | |
| case raw_binary: | | case raw_binary: | |
| load_okay = diskio::load_raw_binary(*this, name, err_msg); | | load_okay = diskio::load_raw_binary(*this, name, err_msg); | |
| break; | | break; | |
| | | | |
| case arma_binary: | | case arma_binary: | |
| load_okay = diskio::load_arma_binary(*this, name, err_msg); | | load_okay = diskio::load_arma_binary(*this, name, err_msg); | |
| break; | | break; | |
| | | | |
| case pgm_binary: | | case pgm_binary: | |
| load_okay = diskio::load_pgm_binary(*this, name, err_msg); | | load_okay = diskio::load_pgm_binary(*this, name, err_msg); | |
| | | | |
| skipping to change at line 3998 | | skipping to change at line 4082 | |
| | | | |
| default: | | default: | |
| arma_warn(print_status, "Mat::load(): unsupported file type"); | | arma_warn(print_status, "Mat::load(): unsupported file type"); | |
| load_okay = false; | | load_okay = false; | |
| } | | } | |
| | | | |
| if( (print_status == true) && (load_okay == false) ) | | if( (print_status == true) && (load_okay == false) ) | |
| { | | { | |
| if(err_msg.length() > 0) | | if(err_msg.length() > 0) | |
| { | | { | |
|
| arma_print("Mat::load(): ", err_msg, name); | | arma_warn(true, "Mat::load(): ", err_msg, name); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| arma_print("Mat::load(): couldn't read ", name); | | arma_warn(true, "Mat::load(): couldn't read ", name); | |
| } | | } | |
| } | | } | |
| | | | |
| if(load_okay == false) | | if(load_okay == false) | |
| { | | { | |
| (*this).reset(); | | (*this).reset(); | |
| } | | } | |
| | | | |
| return load_okay; | | return load_okay; | |
| } | | } | |
| | | | |
| skipping to change at line 4039 | | skipping to change at line 4123 | |
| break; | | break; | |
| | | | |
| case raw_ascii: | | case raw_ascii: | |
| load_okay = diskio::load_raw_ascii(*this, is, err_msg); | | load_okay = diskio::load_raw_ascii(*this, is, err_msg); | |
| break; | | break; | |
| | | | |
| case arma_ascii: | | case arma_ascii: | |
| load_okay = diskio::load_arma_ascii(*this, is, err_msg); | | load_okay = diskio::load_arma_ascii(*this, is, err_msg); | |
| break; | | break; | |
| | | | |
|
| | | case csv_ascii: | |
| | | load_okay = diskio::load_csv_ascii(*this, is, err_msg); | |
| | | break; | |
| | | | |
| case raw_binary: | | case raw_binary: | |
| load_okay = diskio::load_raw_binary(*this, is, err_msg); | | load_okay = diskio::load_raw_binary(*this, is, err_msg); | |
| break; | | break; | |
| | | | |
| case arma_binary: | | case arma_binary: | |
| load_okay = diskio::load_arma_binary(*this, is, err_msg); | | load_okay = diskio::load_arma_binary(*this, is, err_msg); | |
| break; | | break; | |
| | | | |
| case pgm_binary: | | case pgm_binary: | |
| load_okay = diskio::load_pgm_binary(*this, is, err_msg); | | load_okay = diskio::load_pgm_binary(*this, is, err_msg); | |
| | | | |
| skipping to change at line 4060 | | skipping to change at line 4148 | |
| | | | |
| default: | | default: | |
| arma_warn(print_status, "Mat::load(): unsupported file type"); | | arma_warn(print_status, "Mat::load(): unsupported file type"); | |
| load_okay = false; | | load_okay = false; | |
| } | | } | |
| | | | |
| if( (print_status == true) && (load_okay == false) ) | | if( (print_status == true) && (load_okay == false) ) | |
| { | | { | |
| if(err_msg.length() > 0) | | if(err_msg.length() > 0) | |
| { | | { | |
|
| arma_print("Mat::load(): ", err_msg, "the given stream"); | | arma_warn(true, "Mat::load(): ", err_msg, "the given stream"); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| arma_print("Mat::load(): couldn't load from the given stream"); | | arma_warn(true, "Mat::load(): couldn't load from the given stream"); | |
| } | | } | |
| } | | } | |
| | | | |
| if(load_okay == false) | | if(load_okay == false) | |
| { | | { | |
| (*this).reset(); | | (*this).reset(); | |
| } | | } | |
| | | | |
| return load_okay; | | return load_okay; | |
| } | | } | |
| | | | |
End of changes. 30 change blocks. |
| 31 lines changed or deleted | | 120 lines changed or added | |
|
| armadillo | | armadillo | |
| | | | |
| skipping to change at line 192 | | skipping to change at line 192 | |
| #include "armadillo_bits/op_repmat_bones.hpp" | | #include "armadillo_bits/op_repmat_bones.hpp" | |
| #include "armadillo_bits/op_reshape_bones.hpp" | | #include "armadillo_bits/op_reshape_bones.hpp" | |
| #include "armadillo_bits/op_cov_bones.hpp" | | #include "armadillo_bits/op_cov_bones.hpp" | |
| #include "armadillo_bits/op_cor_bones.hpp" | | #include "armadillo_bits/op_cor_bones.hpp" | |
| #include "armadillo_bits/op_shuffle_bones.hpp" | | #include "armadillo_bits/op_shuffle_bones.hpp" | |
| #include "armadillo_bits/op_prod_bones.hpp" | | #include "armadillo_bits/op_prod_bones.hpp" | |
| #include "armadillo_bits/op_pinv_bones.hpp" | | #include "armadillo_bits/op_pinv_bones.hpp" | |
| #include "armadillo_bits/op_dotext_bones.hpp" | | #include "armadillo_bits/op_dotext_bones.hpp" | |
| #include "armadillo_bits/op_flip_bones.hpp" | | #include "armadillo_bits/op_flip_bones.hpp" | |
| #include "armadillo_bits/op_princomp_bones.hpp" | | #include "armadillo_bits/op_princomp_bones.hpp" | |
|
| #include "armadillo_bits/op_princomp_cov_bones.hpp" | | | |
| #include "armadillo_bits/op_misc_bones.hpp" | | #include "armadillo_bits/op_misc_bones.hpp" | |
| #include "armadillo_bits/op_relational_bones.hpp" | | #include "armadillo_bits/op_relational_bones.hpp" | |
| #include "armadillo_bits/op_find_bones.hpp" | | #include "armadillo_bits/op_find_bones.hpp" | |
| #include "armadillo_bits/op_chol_bones.hpp" | | #include "armadillo_bits/op_chol_bones.hpp" | |
| #include "armadillo_bits/op_cx_scalar_bones.hpp" | | #include "armadillo_bits/op_cx_scalar_bones.hpp" | |
| #include "armadillo_bits/op_trimat_bones.hpp" | | #include "armadillo_bits/op_trimat_bones.hpp" | |
| #include "armadillo_bits/op_cumsum_bones.hpp" | | #include "armadillo_bits/op_cumsum_bones.hpp" | |
| #include "armadillo_bits/op_symmat_bones.hpp" | | #include "armadillo_bits/op_symmat_bones.hpp" | |
| | | | |
| #include "armadillo_bits/glue_times_bones.hpp" | | #include "armadillo_bits/glue_times_bones.hpp" | |
| | | | |
| skipping to change at line 336 | | skipping to change at line 335 | |
| #include "armadillo_bits/fn_cor.hpp" | | #include "armadillo_bits/fn_cor.hpp" | |
| #include "armadillo_bits/fn_shuffle.hpp" | | #include "armadillo_bits/fn_shuffle.hpp" | |
| #include "armadillo_bits/fn_prod.hpp" | | #include "armadillo_bits/fn_prod.hpp" | |
| #include "armadillo_bits/fn_eps.hpp" | | #include "armadillo_bits/fn_eps.hpp" | |
| #include "armadillo_bits/fn_pinv.hpp" | | #include "armadillo_bits/fn_pinv.hpp" | |
| #include "armadillo_bits/fn_rank.hpp" | | #include "armadillo_bits/fn_rank.hpp" | |
| #include "armadillo_bits/fn_kron.hpp" | | #include "armadillo_bits/fn_kron.hpp" | |
| #include "armadillo_bits/fn_flip.hpp" | | #include "armadillo_bits/fn_flip.hpp" | |
| #include "armadillo_bits/fn_as_scalar.hpp" | | #include "armadillo_bits/fn_as_scalar.hpp" | |
| #include "armadillo_bits/fn_princomp.hpp" | | #include "armadillo_bits/fn_princomp.hpp" | |
|
| #include "armadillo_bits/fn_princomp_cov.hpp" | | | |
| #include "armadillo_bits/fn_cross.hpp" | | #include "armadillo_bits/fn_cross.hpp" | |
| #include "armadillo_bits/fn_join.hpp" | | #include "armadillo_bits/fn_join.hpp" | |
| #include "armadillo_bits/fn_conv.hpp" | | #include "armadillo_bits/fn_conv.hpp" | |
| #include "armadillo_bits/fn_trunc_exp.hpp" | | #include "armadillo_bits/fn_trunc_exp.hpp" | |
| #include "armadillo_bits/fn_trunc_log.hpp" | | #include "armadillo_bits/fn_trunc_log.hpp" | |
| #include "armadillo_bits/fn_toeplitz.hpp" | | #include "armadillo_bits/fn_toeplitz.hpp" | |
| #include "armadillo_bits/fn_trimat.hpp" | | #include "armadillo_bits/fn_trimat.hpp" | |
| #include "armadillo_bits/fn_cumsum.hpp" | | #include "armadillo_bits/fn_cumsum.hpp" | |
| #include "armadillo_bits/fn_symmat.hpp" | | #include "armadillo_bits/fn_symmat.hpp" | |
| #include "armadillo_bits/fn_syl_lyap.hpp" | | #include "armadillo_bits/fn_syl_lyap.hpp" | |
| | | | |
| skipping to change at line 404 | | skipping to change at line 402 | |
| #include "armadillo_bits/op_repmat_meat.hpp" | | #include "armadillo_bits/op_repmat_meat.hpp" | |
| #include "armadillo_bits/op_reshape_meat.hpp" | | #include "armadillo_bits/op_reshape_meat.hpp" | |
| #include "armadillo_bits/op_cov_meat.hpp" | | #include "armadillo_bits/op_cov_meat.hpp" | |
| #include "armadillo_bits/op_cor_meat.hpp" | | #include "armadillo_bits/op_cor_meat.hpp" | |
| #include "armadillo_bits/op_shuffle_meat.hpp" | | #include "armadillo_bits/op_shuffle_meat.hpp" | |
| #include "armadillo_bits/op_prod_meat.hpp" | | #include "armadillo_bits/op_prod_meat.hpp" | |
| #include "armadillo_bits/op_pinv_meat.hpp" | | #include "armadillo_bits/op_pinv_meat.hpp" | |
| #include "armadillo_bits/op_dotext_meat.hpp" | | #include "armadillo_bits/op_dotext_meat.hpp" | |
| #include "armadillo_bits/op_flip_meat.hpp" | | #include "armadillo_bits/op_flip_meat.hpp" | |
| #include "armadillo_bits/op_princomp_meat.hpp" | | #include "armadillo_bits/op_princomp_meat.hpp" | |
|
| #include "armadillo_bits/op_princomp_cov_meat.hpp" | | | |
| #include "armadillo_bits/op_misc_meat.hpp" | | #include "armadillo_bits/op_misc_meat.hpp" | |
| #include "armadillo_bits/op_relational_meat.hpp" | | #include "armadillo_bits/op_relational_meat.hpp" | |
| #include "armadillo_bits/op_find_meat.hpp" | | #include "armadillo_bits/op_find_meat.hpp" | |
| #include "armadillo_bits/op_chol_meat.hpp" | | #include "armadillo_bits/op_chol_meat.hpp" | |
| #include "armadillo_bits/op_cx_scalar_meat.hpp" | | #include "armadillo_bits/op_cx_scalar_meat.hpp" | |
| #include "armadillo_bits/op_trimat_meat.hpp" | | #include "armadillo_bits/op_trimat_meat.hpp" | |
| #include "armadillo_bits/op_cumsum_meat.hpp" | | #include "armadillo_bits/op_cumsum_meat.hpp" | |
| #include "armadillo_bits/op_symmat_meat.hpp" | | #include "armadillo_bits/op_symmat_meat.hpp" | |
| | | | |
| #include "armadillo_bits/glue_times_meat.hpp" | | #include "armadillo_bits/glue_times_meat.hpp" | |
| | | | |
End of changes. 3 change blocks. |
| 3 lines changed or deleted | | 0 lines changed or added | |
|
| auxlib_meat.hpp | | auxlib_meat.hpp | |
| // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2008-2011 Conrad Sanderson | | // Copyright (C) 2008-2011 Conrad Sanderson | |
|
| // Copyright (C) 2009 Edmund Highcock | | // Copyright (C) 2009 Edmund Highcock | |
| // Copyright (C) 2011 James Sanders | | // Copyright (C) 2011 James Sanders | |
| | | // Copyright (C) 2011 Stanislav Funiak | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| | | | |
| skipping to change at line 44 | | skipping to change at line 45 | |
| if( (N <= 4) && (slow == false) ) | | if( (N <= 4) && (slow == false) ) | |
| { | | { | |
| status = auxlib::inv_inplace_tinymat(out, N); | | status = auxlib::inv_inplace_tinymat(out, N); | |
| } | | } | |
| | | | |
| if( (N > 4) || (status == false) ) | | if( (N > 4) || (status == false) ) | |
| { | | { | |
| status = auxlib::inv_inplace_lapack(out); | | status = auxlib::inv_inplace_lapack(out); | |
| } | | } | |
| | | | |
|
| if(status == false) | | | |
| { | | | |
| arma_print("inv(): matrix appears to be singular" ); | | | |
| out.reset(); | | | |
| } | | | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::inv(Mat<eT>& out, const Mat<eT>& X, const bool slow) | | auxlib::inv(Mat<eT>& out, const Mat<eT>& X, const bool slow) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 77 | | skipping to change at line 72 | |
| { | | { | |
| status = (&out != &X) ? auxlib::inv_noalias_tinymat(out, X, N) : auxlib
::inv_inplace_tinymat(out, N); | | status = (&out != &X) ? auxlib::inv_noalias_tinymat(out, X, N) : auxlib
::inv_inplace_tinymat(out, N); | |
| } | | } | |
| | | | |
| if( (N > 4) || (status == false) ) | | if( (N > 4) || (status == false) ) | |
| { | | { | |
| out = X; | | out = X; | |
| status = auxlib::inv_inplace_lapack(out); | | status = auxlib::inv_inplace_lapack(out); | |
| } | | } | |
| | | | |
|
| if(status == false) | | | |
| { | | | |
| arma_print("inv(): matrix appears to be singular" ); | | | |
| out.reset(); | | | |
| } | | | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::inv_noalias_tinymat(Mat<eT>& out, const Mat<eT>& X, const u32 N) | | auxlib::inv_noalias_tinymat(Mat<eT>& out, const Mat<eT>& X, const u32 N) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 357 | | skipping to change at line 346 | |
| return det_ok; | | return det_ok; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::inv_inplace_lapack(Mat<eT>& out) | | auxlib::inv_inplace_lapack(Mat<eT>& out) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| | | if(out.is_empty()) | |
| | | { | |
| | | return true; | |
| | | } | |
| | | | |
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_ATLAS) | |
| { | | { | |
| podarray<int> ipiv(out.n_rows); | | podarray<int> ipiv(out.n_rows); | |
| | | | |
| int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n
_cols, out.memptr(), out.n_rows, ipiv.memptr()); | | int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n
_cols, out.memptr(), out.n_rows, ipiv.memptr()); | |
| | | | |
| if(info == 0) | | if(info == 0) | |
| { | | { | |
| info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.mem
ptr(), out.n_rows, ipiv.memptr()); | | info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.mem
ptr(), out.n_rows, ipiv.memptr()); | |
| } | | } | |
| | | | |
| skipping to change at line 459 | | skipping to change at line 453 | |
| status = (info == 0); | | status = (info == 0); | |
| } | | } | |
| #else | | #else | |
| { | | { | |
| arma_ignore(layout); | | arma_ignore(layout); | |
| arma_stop("inv(): use of LAPACK needs to be enabled"); | | arma_stop("inv(): use of LAPACK needs to be enabled"); | |
| status = false; | | status = false; | |
| } | | } | |
| #endif | | #endif | |
| | | | |
|
| if(status == false) | | if(status == true) | |
| { | | | |
| arma_print("inv(): matrix appears to be singular" ); | | | |
| out.reset(); | | | |
| } | | | |
| else | | | |
| { | | { | |
| if(layout == 0) | | if(layout == 0) | |
| { | | { | |
| // upper triangular | | // upper triangular | |
| out = trimatu(out); | | out = trimatu(out); | |
| } | | } | |
| else | | else | |
| { | | { | |
| // lower triangular | | // lower triangular | |
| out = trimatl(out); | | out = trimatl(out); | |
| | | | |
| skipping to change at line 501 | | skipping to change at line 490 | |
| | | | |
| if(out.is_empty()) | | if(out.is_empty()) | |
| { | | { | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| bool status; | | bool status; | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
|
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; | |
| blas_int n = blas_int(out.n_rows); | | blas_int n = blas_int(out.n_rows); | |
| blas_int info = 0; | | blas_int lwork = n*n; // TODO: use lwork = -1 to determine optimal size | |
| | | blas_int info = 0; | |
| | | | |
|
| lapack::potrf(&uplo, &n, out.memptr(), &n, &info); | | podarray<blas_int> ipiv; | |
| lapack::potri(&uplo, &n, out.memptr(), &n, &info); | | ipiv.set_size(out.n_rows); | |
| | | | |
|
| out = (layout == 0) ? symmatu(out) : symmatl(out); | | podarray<eT> work; | |
| | | work.set_size( u32(lwork) ); | |
| | | | |
| | | lapack::sytrf(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr() | |
| | | , &lwork, &info); | |
| | | | |
| status = (info == 0); | | status = (info == 0); | |
|
| | | | |
| | | if(status == true) | |
| | | { | |
| | | lapack::sytri(&uplo, &n, out.memptr(), &n, ipiv.memptr(), work.memptr | |
| | | (), &info); | |
| | | | |
| | | out = (layout == 0) ? symmatu(out) : symmatl(out); | |
| | | | |
| | | status = (info == 0); | |
| | | } | |
| } | | } | |
| #else | | #else | |
| { | | { | |
| arma_ignore(layout); | | arma_ignore(layout); | |
| arma_stop("inv(): use of LAPACK needs to be enabled"); | | arma_stop("inv(): use of LAPACK needs to be enabled"); | |
| status = false; | | status = false; | |
| } | | } | |
| #endif | | #endif | |
| | | | |
|
| if(status == false) | | return status; | |
| | | } | |
| | | | |
| | | template<typename eT, typename T1> | |
| | | inline | |
| | | bool | |
| | | auxlib::inv_sympd(Mat<eT>& out, const Base<eT,T1>& X, const u32 layout) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | out = X.get_ref(); | |
| | | | |
| | | arma_debug_check( (out.is_square() == false), "inv(): given matrix is not | |
| | | square" ); | |
| | | | |
| | | if(out.is_empty()) | |
| | | { | |
| | | return true; | |
| | | } | |
| | | | |
| | | bool status; | |
| | | | |
| | | #if defined(ARMA_USE_LAPACK) | |
| | | { | |
| | | char uplo = (layout == 0) ? 'U' : 'L'; | |
| | | blas_int n = blas_int(out.n_rows); | |
| | | blas_int info = 0; | |
| | | | |
| | | lapack::potrf(&uplo, &n, out.memptr(), &n, &info); | |
| | | | |
| | | status = (info == 0); | |
| | | | |
| | | if(status == true) | |
| | | { | |
| | | lapack::potri(&uplo, &n, out.memptr(), &n, &info); | |
| | | | |
| | | out = (layout == 0) ? symmatu(out) : symmatl(out); | |
| | | | |
| | | status = (info == 0); | |
| | | } | |
| | | } | |
| | | #else | |
| { | | { | |
|
| arma_print("inv(): matrix appears to be singular" ); | | arma_ignore(layout); | |
| out.reset(); | | arma_stop("inv(): use of LAPACK needs to be enabled"); | |
| | | status = false; | |
| } | | } | |
|
| | | #endif | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
| eT | | eT | |
| auxlib::det(const Base<eT,T1>& X, const bool slow) | | auxlib::det(const Base<eT,T1>& X, const bool slow) | |
| { | | { | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| | | | |
| skipping to change at line 692 | | skipping to change at line 736 | |
| | | | |
| if(tmp.is_empty()) | | if(tmp.is_empty()) | |
| { | | { | |
| return eT(1); | | return eT(1); | |
| } | | } | |
| | | | |
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_ATLAS) | |
| { | | { | |
| podarray<int> ipiv(tmp.n_rows); | | podarray<int> ipiv(tmp.n_rows); | |
| | | | |
|
| | | //const int info = | |
| atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.
memptr(), tmp.n_rows, ipiv.memptr()); | | atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.
memptr(), tmp.n_rows, ipiv.memptr()); | |
| | | | |
| // on output tmp appears to be L+U_alt, where U_alt is U with the main
diagonal set to zero | | // on output tmp appears to be L+U_alt, where U_alt is U with the main
diagonal set to zero | |
| eT val = tmp.at(0,0); | | eT val = tmp.at(0,0); | |
| for(u32 i=1; i < tmp.n_rows; ++i) | | for(u32 i=1; i < tmp.n_rows; ++i) | |
| { | | { | |
| val *= tmp.at(i,i); | | val *= tmp.at(i,i); | |
| } | | } | |
| | | | |
| int sign = +1; | | int sign = +1; | |
| | | | |
| skipping to change at line 754 | | skipping to change at line 799 | |
| arma_ignore(tmp); | | arma_ignore(tmp); | |
| arma_stop("det(): use of ATLAS or LAPACK needs to be enabled"); | | arma_stop("det(): use of ATLAS or LAPACK needs to be enabled"); | |
| return eT(0); | | return eT(0); | |
| } | | } | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| //! immediate log determinant of a matrix using ATLAS or LAPACK | | //! immediate log determinant of a matrix using ATLAS or LAPACK | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, c
onst Base<eT,T1>& X) | | auxlib::log_det(eT& out_val, typename get_pod_type<eT>::result& out_sign, c
onst Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename get_pod_type<eT>::result T; | | typedef typename get_pod_type<eT>::result T; | |
| | | | |
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_ATLAS) | |
| { | | { | |
| Mat<eT> tmp(X.get_ref()); | | Mat<eT> tmp(X.get_ref()); | |
| arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix
is not square" ); | | arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix
is not square" ); | |
| | | | |
| if(tmp.is_empty()) | | if(tmp.is_empty()) | |
| { | | { | |
| out_val = eT(0); | | out_val = eT(0); | |
| out_sign = T(1); | | out_sign = T(1); | |
| return; | | return; | |
| } | | } | |
| | | | |
| podarray<int> ipiv(tmp.n_rows); | | podarray<int> ipiv(tmp.n_rows); | |
| | | | |
|
| atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.
memptr(), tmp.n_rows, ipiv.memptr()); | | const int info = atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows,
tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr()); | |
| | | | |
| // on output tmp appears to be L+U_alt, where U_alt is U with the main
diagonal set to zero | | // on output tmp appears to be L+U_alt, where U_alt is U with the main
diagonal set to zero | |
| | | | |
| s32 sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.
at(0,0) ) < T(0)) ? -1 : +1 ) : +1; | | s32 sign = (is_complex<eT>::value == false) ? ( (access::tmp_real( tmp.
at(0,0) ) < T(0)) ? -1 : +1 ) : +1; | |
| eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_re
al( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( t
mp.at(0,0) ); | | eT val = (is_complex<eT>::value == false) ? std::log( (access::tmp_re
al( tmp.at(0,0) ) < T(0)) ? tmp.at(0,0)*T(-1) : tmp.at(0,0) ) : std::log( t
mp.at(0,0) ); | |
| | | | |
| for(u32 i=1; i < tmp.n_rows; ++i) | | for(u32 i=1; i < tmp.n_rows; ++i) | |
| { | | { | |
| const eT x = tmp.at(i,i); | | const eT x = tmp.at(i,i); | |
| | | | |
| | | | |
| skipping to change at line 800 | | skipping to change at line 845 | |
| for(u32 i=0; i < tmp.n_rows; ++i) | | for(u32 i=0; i < tmp.n_rows; ++i) | |
| { | | { | |
| if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the
clapack version of getrf() assumes counting from 0 | | if( int(i) != ipiv.mem[i] ) // NOTE: no adjustment required, as the
clapack version of getrf() assumes counting from 0 | |
| { | | { | |
| sign *= -1; | | sign *= -1; | |
| } | | } | |
| } | | } | |
| | | | |
| out_val = val; | | out_val = val; | |
| out_sign = T(sign); | | out_sign = T(sign); | |
|
| | | | |
| | | return (info == 0); | |
| } | | } | |
| #elif defined(ARMA_USE_LAPACK) | | #elif defined(ARMA_USE_LAPACK) | |
| { | | { | |
| Mat<eT> tmp(X.get_ref()); | | Mat<eT> tmp(X.get_ref()); | |
| arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix
is not square" ); | | arma_debug_check( (tmp.is_square() == false), "log_det(): given matrix
is not square" ); | |
| | | | |
| if(tmp.is_empty()) | | if(tmp.is_empty()) | |
| { | | { | |
| out_val = eT(0); | | out_val = eT(0); | |
| out_sign = T(1); | | out_sign = T(1); | |
| | | | |
| skipping to change at line 844 | | skipping to change at line 891 | |
| for(u32 i=0; i < tmp.n_rows; ++i) | | for(u32 i=0; i < tmp.n_rows; ++i) | |
| { | | { | |
| if( blas_int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is
required as Fortran counts from 1 | | if( blas_int(i) != (ipiv.mem[i] - 1) ) // NOTE: adjustment of -1 is
required as Fortran counts from 1 | |
| { | | { | |
| sign *= -1; | | sign *= -1; | |
| } | | } | |
| } | | } | |
| | | | |
| out_val = val; | | out_val = val; | |
| out_sign = T(sign); | | out_sign = T(sign); | |
|
| | | | |
| | | return (info == 0); | |
| } | | } | |
| #else | | #else | |
| { | | { | |
|
| arma_stop("log_det(): use of ATLAS or LAPACK needs to be enabled"); | | | |
| | | | |
| out_val = eT(0); | | out_val = eT(0); | |
| out_sign = T(0); | | out_sign = T(0); | |
|
| | | | |
| | | arma_stop("log_det(): use of ATLAS or LAPACK needs to be enabled"); | |
| | | | |
| | | return false; | |
| } | | } | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| //! immediate LU decomposition of a matrix using ATLAS or LAPACK | | //! immediate LU decomposition of a matrix using ATLAS or LAPACK | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,
T1>& X) | | auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, const Base<eT,
T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| U = X.get_ref(); | | U = X.get_ref(); | |
| | | | |
|
| | | const u32 U_n_rows = U.n_rows; | |
| | | const u32 U_n_cols = U.n_cols; | |
| | | | |
| if(U.is_empty()) | | if(U.is_empty()) | |
| { | | { | |
|
| L.reset(); | | L.set_size(U_n_rows, 0); | |
| U.reset(); | | U.set_size(0, U_n_cols); | |
| ipiv.reset(); | | ipiv.reset(); | |
|
| return; | | return true; | |
| } | | } | |
| | | | |
|
| const u32 U_n_rows = U.n_rows; | | | |
| const u32 U_n_cols = U.n_cols; | | | |
| | | | |
| #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK) | |
| { | | { | |
|
| | | bool status; | |
| | | | |
| #if defined(ARMA_USE_ATLAS) | | #if defined(ARMA_USE_ATLAS) | |
| { | | { | |
| ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | | ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | |
| | | | |
|
| //int info = | | int info = atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_c | |
| atlas::clapack_getrf(atlas::CblasColMajor, U_n_rows, U_n_cols, U.memp | | ols, U.memptr(), U_n_rows, ipiv.memptr()); | |
| tr(), U_n_rows, ipiv.memptr()); | | | |
| | | status = (info == 0); | |
| } | | } | |
| #elif defined(ARMA_USE_LAPACK) | | #elif defined(ARMA_USE_LAPACK) | |
| { | | { | |
| ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | | ipiv.set_size( (std::min)(U_n_rows, U_n_cols) ); | |
| | | | |
| blas_int info = 0; | | blas_int info = 0; | |
| | | | |
| blas_int n_rows = U_n_rows; | | blas_int n_rows = U_n_rows; | |
| blas_int n_cols = U_n_cols; | | blas_int n_cols = U_n_cols; | |
| | | | |
| lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &
info); | | lapack::getrf(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &
info); | |
| | | | |
| // take into account that Fortran counts from 1 | | // take into account that Fortran counts from 1 | |
| arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem); | | arrayops::inplace_minus(ipiv.memptr(), blas_int(1), ipiv.n_elem); | |
|
| | | | |
| | | status = (info == 0); | |
| } | | } | |
| #endif | | #endif | |
| | | | |
| L.copy_size(U); | | L.copy_size(U); | |
| | | | |
| for(u32 col=0; col < U_n_cols; ++col) | | for(u32 col=0; col < U_n_cols; ++col) | |
| { | | { | |
| for(u32 row=0; (row < col) && (row < U_n_rows); ++row) | | for(u32 row=0; (row < col) && (row < U_n_rows); ++row) | |
| { | | { | |
| L.at(row,col) = eT(0); | | L.at(row,col) = eT(0); | |
| | | | |
| skipping to change at line 921 | | skipping to change at line 977 | |
| { | | { | |
| L.at(col,col) = eT(1); | | L.at(col,col) = eT(1); | |
| } | | } | |
| | | | |
| for(u32 row = (col+1); row < U_n_rows; ++row) | | for(u32 row = (col+1); row < U_n_rows; ++row) | |
| { | | { | |
| L.at(row,col) = U.at(row,col); | | L.at(row,col) = U.at(row,col); | |
| U.at(row,col) = eT(0); | | U.at(row,col) = eT(0); | |
| } | | } | |
| } | | } | |
|
| | | | |
| | | return status; | |
| } | | } | |
| #else | | #else | |
| { | | { | |
|
| arma_ignore(U_n_rows); | | | |
| arma_ignore(U_n_cols); | | | |
| arma_stop("lu(): use of ATLAS or LAPACK needs to be enabled"); | | arma_stop("lu(): use of ATLAS or LAPACK needs to be enabled"); | |
|
| | | | |
| | | return false; | |
| } | | } | |
| #endif | | #endif | |
| } | | } | |
| | | | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X) | | auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| podarray<blas_int> ipiv1; | | podarray<blas_int> ipiv1; | |
|
| auxlib::lu(L, U, ipiv1, X); | | const bool status = auxlib::lu(L, U, ipiv1, X); | |
| | | | |
|
| if(U.is_empty()) | | if(status == true) | |
| { | | { | |
|
| L.reset(); | | if(U.is_empty()) | |
| U.reset(); | | { | |
| P.reset(); | | // L and U have been already set to the correct empty matrices | |
| return; | | P.eye(L.n_rows, L.n_rows); | |
| } | | return true; | |
| | | } | |
| const u32 n = ipiv1.n_elem; | | | |
| const u32 P_rows = U.n_rows; | | | |
| | | | |
|
| podarray<blas_int> ipiv2(P_rows); | | const u32 n = ipiv1.n_elem; | |
| | | const u32 P_rows = U.n_rows; | |
| | | | |
|
| const blas_int* ipiv1_mem = ipiv1.memptr(); | | podarray<blas_int> ipiv2(P_rows); | |
| blas_int* ipiv2_mem = ipiv2.memptr(); | | | |
| | | | |
|
| for(u32 i=0; i<P_rows; ++i) | | const blas_int* ipiv1_mem = ipiv1.memptr(); | |
| { | | blas_int* ipiv2_mem = ipiv2.memptr(); | |
| ipiv2_mem[i] = blas_int(i); | | | |
| } | | | |
| | | | |
|
| for(u32 i=0; i<n; ++i) | | for(u32 i=0; i<P_rows; ++i) | |
| { | | { | |
| const u32 k = static_cast<u32>(ipiv1_mem[i]); | | ipiv2_mem[i] = blas_int(i); | |
| | | } | |
| | | | |
|
| if( ipiv2_mem[i] != ipiv2_mem[k] ) | | for(u32 i=0; i<n; ++i) | |
| { | | { | |
|
| std::swap( ipiv2_mem[i], ipiv2_mem[k] ); | | const u32 k = static_cast<u32>(ipiv1_mem[i]); | |
| | | | |
| | | if( ipiv2_mem[i] != ipiv2_mem[k] ) | |
| | | { | |
| | | std::swap( ipiv2_mem[i], ipiv2_mem[k] ); | |
| | | } | |
| } | | } | |
|
| } | | | |
| | | | |
|
| P.zeros(P_rows, P_rows); | | P.zeros(P_rows, P_rows); | |
| | | | |
|
| for(u32 row=0; row<P_rows; ++row) | | for(u32 row=0; row<P_rows; ++row) | |
| { | | { | |
| P.at(row, static_cast<u32>(ipiv2_mem[row])) = eT(1); | | P.at(row, static_cast<u32>(ipiv2_mem[row])) = eT(1); | |
| } | | } | |
| | | | |
|
| if(L.n_cols > U.n_rows) | | if(L.n_cols > U.n_rows) | |
| { | | { | |
| L.shed_cols(U.n_rows, L.n_cols-1); | | L.shed_cols(U.n_rows, L.n_cols-1); | |
| } | | } | |
| | | | |
|
| if(U.n_rows > L.n_cols) | | if(U.n_rows > L.n_cols) | |
| { | | { | |
| U.shed_rows(L.n_cols, U.n_rows-1); | | U.shed_rows(L.n_cols, U.n_rows-1); | |
| | | } | |
| } | | } | |
|
| | | | |
| | | return status; | |
| } | | } | |
| | | | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X) | | auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| podarray<blas_int> ipiv1; | | podarray<blas_int> ipiv1; | |
|
| auxlib::lu(L, U, ipiv1, X); | | const bool status = auxlib::lu(L, U, ipiv1, X); | |
| | | | |
|
| if(U.is_empty()) | | if(status == true) | |
| { | | { | |
|
| L.reset(); | | if(U.is_empty()) | |
| U.reset(); | | { | |
| return; | | // L and U have been already set to the correct empty matrices | |
| } | | return true; | |
| | | } | |
| | | | |
|
| const u32 n = ipiv1.n_elem; | | const u32 n = ipiv1.n_elem; | |
| const u32 P_rows = U.n_rows; | | const u32 P_rows = U.n_rows; | |
| | | | |
|
| podarray<blas_int> ipiv2(P_rows); | | podarray<blas_int> ipiv2(P_rows); | |
| | | | |
|
| const blas_int* ipiv1_mem = ipiv1.memptr(); | | const blas_int* ipiv1_mem = ipiv1.memptr(); | |
| blas_int* ipiv2_mem = ipiv2.memptr(); | | blas_int* ipiv2_mem = ipiv2.memptr(); | |
| | | | |
|
| for(u32 i=0; i<P_rows; ++i) | | for(u32 i=0; i<P_rows; ++i) | |
| { | | { | |
| ipiv2_mem[i] = blas_int(i); | | ipiv2_mem[i] = blas_int(i); | |
| } | | } | |
| | | | |
|
| for(u32 i=0; i<n; ++i) | | for(u32 i=0; i<n; ++i) | |
| { | | { | |
| const u32 k = static_cast<u32>(ipiv1_mem[i]); | | const u32 k = static_cast<u32>(ipiv1_mem[i]); | |
| | | | |
| | | if( ipiv2_mem[i] != ipiv2_mem[k] ) | |
| | | { | |
| | | std::swap( ipiv2_mem[i], ipiv2_mem[k] ); | |
| | | L.swap_rows( static_cast<u32>(ipiv2_mem[i]), static_cast<u32>(ipiv2 | |
| | | _mem[k]) ); | |
| | | } | |
| | | } | |
| | | | |
|
| if( ipiv2_mem[i] != ipiv2_mem[k] ) | | if(L.n_cols > U.n_rows) | |
| { | | { | |
|
| std::swap( ipiv2_mem[i], ipiv2_mem[k] ); | | L.shed_cols(U.n_rows, L.n_cols-1); | |
| L.swap_rows( static_cast<u32>(ipiv2_mem[i]), static_cast<u32>(ipiv2_m | | | |
| em[k]) ); | | | |
| } | | } | |
|
| } | | | |
| | | | |
|
| if(L.n_cols > U.n_rows) | | if(U.n_rows > L.n_cols) | |
| { | | { | |
| L.shed_cols(U.n_rows, L.n_cols-1); | | U.shed_rows(L.n_cols, U.n_rows-1); | |
| | | } | |
| } | | } | |
| | | | |
|
| if(U.n_rows > L.n_cols) | | return status; | |
| { | | | |
| U.shed_rows(L.n_cols, U.n_rows-1); | | | |
| } | | | |
| } | | } | |
| | | | |
| //! immediate eigenvalues of a symmetric real matrix using LAPACK | | //! immediate eigenvalues of a symmetric real matrix using LAPACK | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::eig_sym(Col<eT>& eigval, const Base<eT,T1>& X) | | auxlib::eig_sym(Col<eT>& eigval, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 1501 | | skipping to change at line 1567 | |
| inline | | inline | |
| bool | | bool | |
| auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X) | | auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| R = X.get_ref(); | | R = X.get_ref(); | |
| | | | |
|
| | | const u32 R_n_rows = R.n_rows; | |
| | | const u32 R_n_cols = R.n_cols; | |
| | | | |
| if(R.is_empty()) | | if(R.is_empty()) | |
| { | | { | |
|
| Q.reset(); | | Q.eye(R_n_rows, R_n_rows); | |
| R.reset(); | | | |
| return true; | | return true; | |
| } | | } | |
| | | | |
|
| const u32 R_n_rows = R.n_rows; | | | |
| const u32 R_n_cols = R.n_cols; | | | |
| | | | |
| blas_int m = static_cast<blas_int>(R_n_rows); | | blas_int m = static_cast<blas_int>(R_n_rows); | |
| blas_int n = static_cast<blas_int>(R_n_cols); | | blas_int n = static_cast<blas_int>(R_n_cols); | |
| blas_int work_len = (std::max)(blas_int(1),n); | | blas_int work_len = (std::max)(blas_int(1),n); | |
| blas_int work_len_tmp; | | blas_int work_len_tmp; | |
| blas_int k = (std::min)(m,n); | | blas_int k = (std::min)(m,n); | |
| blas_int info; | | blas_int info; | |
| | | | |
| podarray<eT> tau( static_cast<u32>(k) ); | | podarray<eT> tau( static_cast<u32>(k) ); | |
| podarray<eT> work( static_cast<u32>(work_len) ); | | podarray<eT> work( static_cast<u32>(work_len) ); | |
| | | | |
| | | | |
| skipping to change at line 1608 | | skipping to change at line 1673 | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| Mat<eT> A(X.get_ref()); | | Mat<eT> A(X.get_ref()); | |
| | | | |
| X_n_rows = A.n_rows; | | X_n_rows = A.n_rows; | |
| X_n_cols = A.n_cols; | | X_n_cols = A.n_cols; | |
| | | | |
| if(A.is_empty()) | | if(A.is_empty()) | |
| { | | { | |
| S.reset(); | | S.reset(); | |
|
| X_n_rows = 0; | | | |
| X_n_cols = 0; | | | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| Mat<eT> U(1, 1); | | Mat<eT> U(1, 1); | |
| Mat<eT> V(1, A.n_cols); | | Mat<eT> V(1, A.n_cols); | |
| | | | |
| char jobu = 'N'; | | char jobu = 'N'; | |
| char jobvt = 'N'; | | char jobvt = 'N'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| | | | |
| skipping to change at line 1702 | | skipping to change at line 1765 | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| Mat<eT> A(X.get_ref()); | | Mat<eT> A(X.get_ref()); | |
| | | | |
| X_n_rows = A.n_rows; | | X_n_rows = A.n_rows; | |
| X_n_cols = A.n_cols; | | X_n_cols = A.n_cols; | |
| | | | |
| if(A.is_empty()) | | if(A.is_empty()) | |
| { | | { | |
| S.reset(); | | S.reset(); | |
|
| X_n_rows = 0; | | | |
| X_n_cols = 0; | | | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| Mat<eT> U(1, 1); | | Mat<eT> U(1, 1); | |
| Mat<eT> V(1, A.n_cols); | | Mat<eT> V(1, A.n_cols); | |
| | | | |
| char jobu = 'N'; | | char jobu = 'N'; | |
| char jobvt = 'N'; | | char jobvt = 'N'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| | | | |
| skipping to change at line 1815 | | skipping to change at line 1876 | |
| auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X) | | auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Base<eT,T1>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| Mat<eT> A(X.get_ref()); | | Mat<eT> A(X.get_ref()); | |
| | | | |
| if(A.is_empty()) | | if(A.is_empty()) | |
| { | | { | |
|
| U.reset(); | | U.eye(A.n_rows, A.n_rows); | |
| S.reset(); | | S.reset(); | |
|
| V.reset(); | | V.eye(A.n_cols, A.n_cols); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); | |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); | |
| | | | |
| char jobu = 'A'; | | char jobu = 'A'; | |
| char jobvt = 'A'; | | char jobvt = 'A'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| | | | |
| skipping to change at line 1906 | | skipping to change at line 1967 | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| Mat<eT> A(X.get_ref()); | | Mat<eT> A(X.get_ref()); | |
| | | | |
| if(A.is_empty()) | | if(A.is_empty()) | |
| { | | { | |
|
| U.reset(); | | U.eye(A.n_rows, A.n_rows); | |
| S.reset(); | | S.reset(); | |
|
| V.reset(); | | V.eye(A.n_cols, A.n_cols); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| U.set_size(A.n_rows, A.n_rows); | | U.set_size(A.n_rows, A.n_rows); | |
| V.set_size(A.n_cols, A.n_cols); | | V.set_size(A.n_cols, A.n_cols); | |
| | | | |
| char jobu = 'A'; | | char jobu = 'A'; | |
| char jobvt = 'A'; | | char jobvt = 'A'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| | | | |
| skipping to change at line 1997 | | skipping to change at line 2058 | |
| //! Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows | | //! Assumes that A.n_rows = A.n_cols and B.n_rows = A.n_rows | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::solve(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B, const bool slow) | | auxlib::solve(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B, const bool slow) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) | |
| { | | { | |
|
| out.reset(); | | out.zeros(A.n_cols, B.n_cols); | |
| A.reset(); | | | |
| return true; | | return true; | |
| } | | } | |
| else | | else | |
| { | | { | |
| const u32 A_n_rows = A.n_rows; | | const u32 A_n_rows = A.n_rows; | |
| | | | |
| bool status = false; | | bool status = false; | |
| | | | |
| if( (A_n_rows <= 4) && (slow == false) ) | | if( (A_n_rows <= 4) && (slow == false) ) | |
| { | | { | |
| | | | |
| skipping to change at line 2066 | | skipping to change at line 2126 | |
| inline | | inline | |
| bool | | bool | |
| auxlib::solve_od(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B) | | auxlib::solve_od(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) | |
| { | | { | |
|
| out.reset(); | | out.zeros(A.n_cols, B.n_cols); | |
| A.reset(); | | | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| char trans = 'N'; | | char trans = 'N'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| blas_int n = A.n_cols; | | blas_int n = A.n_cols; | |
| blas_int lda = A.n_rows; | | blas_int lda = A.n_rows; | |
| blas_int ldb = A.n_rows; | | blas_int ldb = A.n_rows; | |
| blas_int nrhs = B.n_cols; | | blas_int nrhs = B.n_cols; | |
| | | | |
| skipping to change at line 2133 | | skipping to change at line 2192 | |
| inline | | inline | |
| bool | | bool | |
| auxlib::solve_ud(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B) | | auxlib::solve_ud(Mat<eT>& out, Mat<eT>& A, const Mat<eT>& B) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) | |
| { | | { | |
|
| out.reset(); | | out.zeros(A.n_cols, B.n_cols); | |
| A.reset(); | | | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| char trans = 'N'; | | char trans = 'N'; | |
| | | | |
| blas_int m = A.n_rows; | | blas_int m = A.n_rows; | |
| blas_int n = A.n_cols; | | blas_int n = A.n_cols; | |
| blas_int lda = A.n_rows; | | blas_int lda = A.n_rows; | |
| blas_int ldb = A.n_cols; | | blas_int ldb = A.n_cols; | |
| blas_int nrhs = B.n_cols; | | blas_int nrhs = B.n_cols; | |
| | | | |
| skipping to change at line 2214 | | skipping to change at line 2272 | |
| inline | | inline | |
| bool | | bool | |
| auxlib::solve_tr(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const u3
2 layout) | | auxlib::solve_tr(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const u3
2 layout) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| #if defined(ARMA_USE_LAPACK) | | #if defined(ARMA_USE_LAPACK) | |
| { | | { | |
| if(A.is_empty() || B.is_empty()) | | if(A.is_empty() || B.is_empty()) | |
| { | | { | |
|
| out.reset(); | | out.zeros(A.n_cols, B.n_cols); | |
| return true; | | return true; | |
| } | | } | |
| | | | |
| out = B; | | out = B; | |
| | | | |
| char uplo = (layout == 0) ? 'U' : 'L'; | | char uplo = (layout == 0) ? 'U' : 'L'; | |
| char trans = 'N'; | | char trans = 'N'; | |
| char diag = 'N'; | | char diag = 'N'; | |
| blas_int n = blas_int(A.n_rows); | | blas_int n = blas_int(A.n_rows); | |
| blas_int nrhs = blas_int(B.n_cols); | | blas_int nrhs = blas_int(B.n_cols); | |
| | | | |
| skipping to change at line 2420 | | skipping to change at line 2478 | |
| | | | |
| status = (info == 0); | | status = (info == 0); | |
| } | | } | |
| #else | | #else | |
| { | | { | |
| arma_stop("syl(): use of LAPACK needs to be enabled"); | | arma_stop("syl(): use of LAPACK needs to be enabled"); | |
| return false; | | return false; | |
| } | | } | |
| #endif | | #endif | |
| | | | |
|
| if(status == false) | | | |
| { | | | |
| arma_print("syl(): equation appears to be singular" ); | | | |
| X.reset(); | | | |
| } | | | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| // | | // | |
| // lyap (solution of the continuous Lyapunov equation AX + XA^H + Q = 0) | | // lyap (solution of the continuous Lyapunov equation AX + XA^H + Q = 0) | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| auxlib::lyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q) | | auxlib::lyap(Mat<eT>& X, const Mat<eT>& A, const Mat<eT>& Q) | |
| | | | |
End of changes. 74 change blocks. |
| 138 lines changed or deleted | | 193 lines changed or added | |
|
| debug.hpp | | debug.hpp | |
| // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2008-2011 Conrad Sanderson | | // Copyright (C) 2008-2011 Conrad Sanderson | |
|
| | | // Copyright (C) 2011 Stanislav Funiak | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| //! \addtogroup debug | | //! \addtogroup debug | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::ostream& | | std::ostream& | |
|
| arma_log_stream(std::ostream* user_stream) | | arma_stream_err1(std::ostream* user_stream) | |
| { | | { | |
|
| static std::ostream* log_stream = &(std::cout); | | static std::ostream* stream_err1 = &(std::cout); | |
| | | | |
| if(user_stream != NULL) | | if(user_stream != NULL) | |
| { | | { | |
|
| log_stream = user_stream; | | stream_err1 = user_stream; | |
| } | | } | |
| | | | |
|
| return *log_stream; | | return *stream_err1; | |
| | | } | |
| | | | |
| | | template<typename T> | |
| | | inline | |
| | | std::ostream& | |
| | | arma_stream_err2(std::ostream* user_stream) | |
| | | { | |
| | | static std::ostream* stream_err2 = &(std::cout); | |
| | | | |
| | | if(user_stream != NULL) | |
| | | { | |
| | | stream_err2 = user_stream; | |
| | | } | |
| | | | |
| | | return *stream_err2; | |
| | | } | |
| | | | |
| | | inline | |
| | | void | |
| | | set_stream_err1(std::ostream& user_stream) | |
| | | { | |
| | | arma_stream_err1<char>(&user_stream); | |
| } | | } | |
| | | | |
| inline | | inline | |
| void | | void | |
|
| set_log_stream(std::ostream& user_stream) | | set_stream_err2(std::ostream& user_stream) | |
| { | | { | |
|
| arma_log_stream<char>(&user_stream); | | arma_stream_err2<char>(&user_stream); | |
| } | | } | |
| | | | |
| inline | | inline | |
| std::ostream& | | std::ostream& | |
|
| get_log_stream() | | get_stream_err1() | |
| { | | { | |
|
| return arma_log_stream<char>(NULL); | | return arma_stream_err1<char>(NULL); | |
| | | } | |
| | | | |
| | | inline | |
| | | std::ostream& | |
| | | get_stream_err2() | |
| | | { | |
| | | return arma_stream_err2<char>(NULL); | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_stop | | // arma_stop | |
| | | | |
|
| //! print a message to get_log_stream() and/or throw a run-time error excep
tion | | //! print a message to get_stream_err1() and/or throw a logic_error excepti
on | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| arma_cold | | arma_cold | |
| arma_stop(const T1& x) | | arma_stop(const T1& x) | |
| { | | { | |
|
| std::ostream& log_stream = get_log_stream(); | | #if defined(ARMA_PRINT_LOGIC_ERRORS) | |
| | | { | |
| | | std::ostream& out = get_stream_err1(); | |
| | | | |
| | | out.flush(); | |
| | | | |
| | | out << '\n'; | |
| | | out << "error: " << x << '\n'; | |
| | | out << '\n'; | |
| | | out.flush(); | |
| | | } | |
| | | #else | |
| | | { | |
| | | arma_ignore(x); | |
| | | } | |
| | | #endif | |
| | | | |
| | | throw std::logic_error(""); | |
| | | } | |
| | | | |
| | | template<typename T1> | |
| | | inline | |
| | | void | |
| | | arma_cold | |
| | | arma_stop_bad_alloc(const T1& x) | |
| | | { | |
| | | std::ostream& out = get_stream_err1(); | |
| | | | |
| | | out.flush(); | |
| | | | |
| | | out << '\n'; | |
| | | out << "error: " << x << '\n'; | |
| | | out << '\n'; | |
| | | out.flush(); | |
| | | | |
| | | throw std::bad_alloc(); | |
| | | } | |
| | | | |
| | | // | |
| | | // arma_bad | |
| | | | |
| | | //! print a message to get_stream_err2() and/or throw a run-time error exce | |
| | | ption | |
| | | template<typename T1> | |
| | | inline | |
| | | void | |
| | | arma_cold | |
| | | arma_bad(const T1& x, const bool hurl = true) | |
| | | { | |
| | | #if defined(ARMA_PRINT_RUNTIME_ERRORS) | |
| | | { | |
| | | std::ostream& out = get_stream_err2(); | |
| | | | |
|
| log_stream.flush(); | | out.flush(); | |
| | | | |
|
| log_stream << '\n'; | | out << '\n'; | |
| log_stream << "run-time error: " << x << '\n'; | | out << "error: " << x << '\n'; | |
| log_stream << '\n'; | | out << '\n'; | |
| log_stream.flush(); | | out.flush(); | |
| | | } | |
| | | #else | |
| | | { | |
| | | arma_ignore(x); | |
| | | } | |
| | | #endif | |
| | | | |
|
| throw std::runtime_error(""); | | if(hurl == true) | |
| | | { | |
| | | throw std::runtime_error(""); | |
| | | } | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_print | | // arma_print | |
| | | | |
| inline | | inline | |
| void | | void | |
| arma_cold | | arma_cold | |
| arma_print() | | arma_print() | |
| { | | { | |
|
| get_log_stream() << std::endl; | | get_stream_err1() << std::endl; | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| arma_cold | | arma_cold | |
| arma_print(const T1& x) | | arma_print(const T1& x) | |
| { | | { | |
|
| get_log_stream() << x << std::endl; | | get_stream_err1() << x << std::endl; | |
| } | | } | |
| | | | |
| template<typename T1, typename T2> | | template<typename T1, typename T2> | |
| inline | | inline | |
| void | | void | |
| arma_cold | | arma_cold | |
| arma_print(const T1& x, const T2& y) | | arma_print(const T1& x, const T2& y) | |
| { | | { | |
|
| get_log_stream() << x << y << std::endl; | | get_stream_err1() << x << y << std::endl; | |
| } | | } | |
| | | | |
| template<typename T1, typename T2, typename T3> | | template<typename T1, typename T2, typename T3> | |
| inline | | inline | |
| void | | void | |
| arma_cold | | arma_cold | |
| arma_print(const T1& x, const T2& y, const T3& z) | | arma_print(const T1& x, const T2& y, const T3& z) | |
| { | | { | |
|
| get_log_stream() << x << y << z << std::endl; | | get_stream_err1() << x << y << z << std::endl; | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_sigprint | | // arma_sigprint | |
| | | | |
| //! print a message the the log stream with a preceding @ character. | | //! print a message the the log stream with a preceding @ character. | |
| //! by default the log stream is cout. | | //! by default the log stream is cout. | |
| //! used for printing the signature of a function | | //! used for printing the signature of a function | |
| //! (see the arma_extra_debug_sigprint macro) | | //! (see the arma_extra_debug_sigprint macro) | |
| inline | | inline | |
| void | | void | |
| arma_sigprint(const char* x) | | arma_sigprint(const char* x) | |
| { | | { | |
|
| get_log_stream() << "@ " << x; | | get_stream_err1() << "@ " << x; | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_bktprint | | // arma_bktprint | |
| | | | |
| inline | | inline | |
| void | | void | |
| arma_bktprint() | | arma_bktprint() | |
| { | | { | |
|
| get_log_stream() << std::endl; | | get_stream_err1() << std::endl; | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| arma_bktprint(const T1& x) | | arma_bktprint(const T1& x) | |
| { | | { | |
|
| get_log_stream() << " [" << x << ']' << std::endl; | | get_stream_err1() << " [" << x << ']' << std::endl; | |
| } | | } | |
| | | | |
| template<typename T1, typename T2> | | template<typename T1, typename T2> | |
| inline | | inline | |
| void | | void | |
| arma_bktprint(const T1& x, const T2& y) | | arma_bktprint(const T1& x, const T2& y) | |
| { | | { | |
|
| get_log_stream() << " [" << x << y << ']' << std::endl; | | get_stream_err1() << " [" << x << y << ']' << std::endl; | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_thisprint | | // arma_thisprint | |
| | | | |
| inline | | inline | |
| void | | void | |
|
| arma_thisprint(void* this_ptr) | | arma_thisprint(const void* this_ptr) | |
| { | | { | |
|
| get_log_stream() << " [this = " << this_ptr << ']' << std::endl; | | get_stream_err1() << " [this = " << this_ptr << ']' << std::endl; | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_warn | | // arma_warn | |
| | | | |
|
| //! if state is true, print a message on cout | | //! print a message to the warn stream | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
|
| arma_hot | | arma_cold | |
| arma_warn(const bool state, const T1& x) | | arma_warn(const bool state, const T1& x) | |
| { | | { | |
| if(state==true) | | if(state==true) | |
| { | | { | |
|
| arma_print(x); | | get_stream_err2() << x << std::endl; | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1, typename T2> | | template<typename T1, typename T2> | |
| inline | | inline | |
| void | | void | |
|
| arma_hot | | arma_cold | |
| arma_warn(const bool state, const T1& x, const T2& y) | | arma_warn(const bool state, const T1& x, const T2& y) | |
| { | | { | |
| if(state==true) | | if(state==true) | |
| { | | { | |
|
| arma_print(x,y); | | get_stream_err2() << x << y << std::endl; | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename T1, typename T2, typename T3> | |
| | | inline | |
| | | void | |
| | | arma_cold | |
| | | arma_warn(const bool state, const T1& x, const T2& y, const T3& z) | |
| | | { | |
| | | if(state==true) | |
| | | { | |
| | | get_stream_err2() << x << y << z << std::endl; | |
| } | | } | |
| } | | } | |
| | | | |
| // | | // | |
| // arma_check | | // arma_check | |
| | | | |
| //! if state is true, abort program | | //! if state is true, abort program | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| | | | |
| skipping to change at line 211 | | skipping to change at line 312 | |
| void | | void | |
| arma_hot | | arma_hot | |
| arma_check(const bool state, const T1& x, const T2& y) | | arma_check(const bool state, const T1& x, const T2& y) | |
| { | | { | |
| if(state==true) | | if(state==true) | |
| { | | { | |
| arma_stop( std::string(x) + std::string(y) ); | | arma_stop( std::string(x) + std::string(y) ); | |
| } | | } | |
| } | | } | |
| | | | |
|
| | | template<typename T1> | |
| | | inline | |
| | | void | |
| | | arma_hot | |
| | | arma_check_bad_alloc(const bool state, const T1& x) | |
| | | { | |
| | | if(state==true) | |
| | | { | |
| | | arma_stop_bad_alloc(x); | |
| | | } | |
| | | } | |
| | | | |
| // | | // | |
| // arma_set_error | | // arma_set_error | |
| | | | |
| arma_inline | | arma_inline | |
| void | | void | |
| arma_hot | | arma_hot | |
| arma_set_error(bool& err_state, char*& err_msg, const bool expression, cons
t char* message) | | arma_set_error(bool& err_state, char*& err_msg, const bool expression, cons
t char* message) | |
| { | | { | |
| if(expression == true) | | if(expression == true) | |
| { | | { | |
| | | | |
| skipping to change at line 901 | | skipping to change at line 1014 | |
| { | | { | |
| unsigned short a; | | unsigned short a; | |
| unsigned char b[sizeof(unsigned short)]; | | unsigned char b[sizeof(unsigned short)]; | |
| } endian_test; | | } endian_test; | |
| | | | |
| endian_test.a = 1; | | endian_test.a = 1; | |
| | | | |
| const bool little_endian = (endian_test.b[0] == 1); | | const bool little_endian = (endian_test.b[0] == 1); | |
| const char* nickname = ARMA_VERSION_NAME; | | const char* nickname = ARMA_VERSION_NAME; | |
| | | | |
|
| get_log_stream() << "@ ---" << '\n'; | | std::ostream& out = get_stream_err1(); | |
| get_log_stream() << "@ Armadillo " | | | |
| << arma_version::major << '.' << arma_version::minor << ' | | out << "@ ---" << '\n'; | |
| .' << arma_version::patch | | out << "@ Armadillo " | |
| << " (" << nickname << ")\n"; | | << arma_version::major << '.' << arma_version::minor << '.' << | |
| | | arma_version::patch | |
| get_log_stream() << "@ arma_config::mat_prealloc = " << arma_config | | << " (" << nickname << ")\n"; | |
| ::mat_prealloc << " element(s)\n"; | | | |
| get_log_stream() << "@ arma_config::atlas = " << arma_config | | out << "@ arma_config::mat_prealloc = " << arma_config::mat_preallo | |
| ::atlas << '\n'; | | c << " element(s)\n"; | |
| get_log_stream() << "@ arma_config::lapack = " << arma_config | | out << "@ arma_config::atlas = " << arma_config::atlas | |
| ::lapack << '\n'; | | << '\n'; | |
| get_log_stream() << "@ arma_config::blas = " << arma_config | | out << "@ arma_config::lapack = " << arma_config::lapack | |
| ::blas << '\n'; | | << '\n'; | |
| get_log_stream() << "@ arma_config::boost = " << arma_config | | out << "@ arma_config::blas = " << arma_config::blas | |
| ::boost << '\n'; | | << '\n'; | |
| get_log_stream() << "@ arma_config::boost_date = " << arma_config | | out << "@ arma_config::boost = " << arma_config::boost | |
| ::boost_date << '\n'; | | << '\n'; | |
| get_log_stream() << "@ arma_config::good_comp = " << arma_config | | out << "@ arma_config::boost_date = " << arma_config::boost_date | |
| ::good_comp << '\n'; | | << '\n'; | |
| get_log_stream() << "@ arma_config::extra_code = " << arma_config | | out << "@ arma_config::good_comp = " << arma_config::good_comp | |
| ::extra_code << '\n'; | | << '\n'; | |
| get_log_stream() << "@ sizeof(void*) = " << sizeof(void*) << | | out << "@ arma_config::extra_code = " << arma_config::extra_code | |
| '\n'; | | << '\n'; | |
| get_log_stream() << "@ sizeof(int) = " << sizeof(int) << | | out << "@ sizeof(void*) = " << sizeof(void*) << '\n'; | |
| '\n'; | | out << "@ sizeof(int) = " << sizeof(int) << '\n'; | |
| get_log_stream() << "@ sizeof(long) = " << sizeof(long) << | | out << "@ sizeof(long) = " << sizeof(long) << '\n'; | |
| '\n'; | | out << "@ sizeof(blas_int) = " << sizeof(blas_int) << '\n'; | |
| get_log_stream() << "@ sizeof(blas_int) = " << sizeof(blas_int) << | | out << "@ little_endian = " << little_endian << '\n'; | |
| '\n'; | | out << "@ ---" << std::endl; | |
| get_log_stream() << "@ little_endian = " << little_endian << | | | |
| '\n'; | | | |
| get_log_stream() << "@ ---" << std::endl; | | | |
| } | | } | |
| | | | |
| }; | | }; | |
| | | | |
| static arma_first_extra_debug_message arma_first_extra_debug_message_ru
n; | | static arma_first_extra_debug_message arma_first_extra_debug_message_ru
n; | |
| } | | } | |
| | | | |
| #endif | | #endif | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 31 change blocks. |
| 64 lines changed or deleted | | 175 lines changed or added | |
|
| diskio_bones.hpp | | diskio_bones.hpp | |
| | | | |
| skipping to change at line 28 | | skipping to change at line 28 | |
| class diskio | | class diskio | |
| { | | { | |
| public: | | public: | |
| | | | |
| template<typename eT> inline static std::string gen_txt_header(const Mat<
eT>& x); | | template<typename eT> inline static std::string gen_txt_header(const Mat<
eT>& x); | |
| template<typename eT> inline static std::string gen_bin_header(const Mat<
eT>& x); | | template<typename eT> inline static std::string gen_bin_header(const Mat<
eT>& x); | |
| | | | |
| template<typename eT> inline static std::string gen_txt_header(const Cube
<eT>& x); | | template<typename eT> inline static std::string gen_txt_header(const Cube
<eT>& x); | |
| template<typename eT> inline static std::string gen_bin_header(const Cube
<eT>& x); | | template<typename eT> inline static std::string gen_bin_header(const Cube
<eT>& x); | |
| | | | |
|
| inline static bool is_raw_binary(std::istream& f); | | inline static file_type guess_file_type(std::istream& f); | |
| | | | |
| inline static char conv_to_hex_char(const u8 x); | | inline static char conv_to_hex_char(const u8 x); | |
| inline static void conv_to_hex(char* out, const u8 x); | | inline static void conv_to_hex(char* out, const u8 x); | |
| | | | |
| inline static std::string gen_tmp_name(const std::string& x); | | inline static std::string gen_tmp_name(const std::string& x); | |
| | | | |
| inline static bool safe_rename(const std::string& old_name, const std::st
ring& new_name); | | inline static bool safe_rename(const std::string& old_name, const std::st
ring& new_name); | |
| | | | |
| // | | // | |
| // matrix saving | | // matrix saving | |
| | | | |
| template<typename eT> inline static bool save_raw_ascii (const Mat<eT>&
x, const std::string& final_name); | | template<typename eT> inline static bool save_raw_ascii (const Mat<eT>&
x, const std::string& final_name); | |
| template<typename eT> inline static bool save_raw_binary (const Mat<eT>&
x, const std::string& final_name); | | template<typename eT> inline static bool save_raw_binary (const Mat<eT>&
x, const std::string& final_name); | |
| template<typename eT> inline static bool save_arma_ascii (const Mat<eT>&
x, const std::string& final_name); | | template<typename eT> inline static bool save_arma_ascii (const Mat<eT>&
x, const std::string& final_name); | |
|
| | | template<typename eT> inline static bool save_csv_ascii (const Mat<eT>&
x, const std::string& final_name); | |
| template<typename eT> inline static bool save_arma_binary(const Mat<eT>&
x, const std::string& final_name); | | template<typename eT> inline static bool save_arma_binary(const Mat<eT>&
x, const std::string& final_name); | |
| template<typename eT> inline static bool save_pgm_binary (const Mat<eT>&
x, const std::string& final_name); | | template<typename eT> inline static bool save_pgm_binary (const Mat<eT>&
x, const std::string& final_name); | |
| template<typename T> inline static bool save_pgm_binary (const Mat< std:
:complex<T> >& x, const std::string& final_name); | | template<typename T> inline static bool save_pgm_binary (const Mat< std:
:complex<T> >& x, const std::string& final_name); | |
| | | | |
| template<typename eT> inline static bool save_raw_ascii (const Mat<eT>&
x, std::ostream& f); | | template<typename eT> inline static bool save_raw_ascii (const Mat<eT>&
x, std::ostream& f); | |
| template<typename eT> inline static bool save_raw_binary (const Mat<eT>&
x, std::ostream& f); | | template<typename eT> inline static bool save_raw_binary (const Mat<eT>&
x, std::ostream& f); | |
| template<typename eT> inline static bool save_arma_ascii (const Mat<eT>&
x, std::ostream& f); | | template<typename eT> inline static bool save_arma_ascii (const Mat<eT>&
x, std::ostream& f); | |
|
| | | template<typename eT> inline static bool save_csv_ascii (const Mat<eT>&
x, std::ostream& f); | |
| template<typename eT> inline static bool save_arma_binary(const Mat<eT>&
x, std::ostream& f); | | template<typename eT> inline static bool save_arma_binary(const Mat<eT>&
x, std::ostream& f); | |
| template<typename eT> inline static bool save_pgm_binary (const Mat<eT>&
x, std::ostream& f); | | template<typename eT> inline static bool save_pgm_binary (const Mat<eT>&
x, std::ostream& f); | |
| template<typename T> inline static bool save_pgm_binary (const Mat< std:
:complex<T> >& x, std::ostream& f); | | template<typename T> inline static bool save_pgm_binary (const Mat< std:
:complex<T> >& x, std::ostream& f); | |
| | | | |
| // | | // | |
| // matrix loading | | // matrix loading | |
| | | | |
| template<typename eT> inline static bool load_raw_ascii (Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_raw_ascii (Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| template<typename eT> inline static bool load_raw_binary (Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_raw_binary (Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| template<typename eT> inline static bool load_arma_ascii (Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_arma_ascii (Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
|
| | | template<typename eT> inline static bool load_csv_ascii (Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| template<typename eT> inline static bool load_arma_binary(Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_arma_binary(Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| template<typename eT> inline static bool load_pgm_binary (Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_pgm_binary (Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| template<typename T> inline static bool load_pgm_binary (Mat< std::compl
ex<T> >& x, const std::string& name, std::string& err_msg); | | template<typename T> inline static bool load_pgm_binary (Mat< std::compl
ex<T> >& x, const std::string& name, std::string& err_msg); | |
| template<typename eT> inline static bool load_auto_detect(Mat<eT>&
x, const std::string& name, std::string& err_msg); | | template<typename eT> inline static bool load_auto_detect(Mat<eT>&
x, const std::string& name, std::string& err_msg); | |
| | | | |
| template<typename eT> inline static bool load_raw_ascii (Mat<eT>&
x, std::istream& f, std::string& err_msg); | | template<typename eT> inline static bool load_raw_ascii (Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
| template<typename eT> inline static bool load_raw_binary (Mat<eT>&
x, std::istream& f, std::string& err_msg); | | template<typename eT> inline static bool load_raw_binary (Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
| template<typename eT> inline static bool load_arma_ascii (Mat<eT>&
x, std::istream& f, std::string& err_msg); | | template<typename eT> inline static bool load_arma_ascii (Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
|
| | | template<typename eT> inline static bool load_csv_ascii (Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
| template<typename eT> inline static bool load_arma_binary(Mat<eT>&
x, std::istream& f, std::string& err_msg); | | template<typename eT> inline static bool load_arma_binary(Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
| template<typename eT> inline static bool load_pgm_binary (Mat<eT>&
x, std::istream& is, std::string& err_msg); | | template<typename eT> inline static bool load_pgm_binary (Mat<eT>&
x, std::istream& is, std::string& err_msg); | |
| template<typename T> inline static bool load_pgm_binary (Mat< std::compl
ex<T> >& x, std::istream& is, std::string& err_msg); | | template<typename T> inline static bool load_pgm_binary (Mat< std::compl
ex<T> >& x, std::istream& is, std::string& err_msg); | |
| template<typename eT> inline static bool load_auto_detect(Mat<eT>&
x, std::istream& f, std::string& err_msg); | | template<typename eT> inline static bool load_auto_detect(Mat<eT>&
x, std::istream& f, std::string& err_msg); | |
| | | | |
| inline static void pnm_skip_comments(std::istream& f); | | inline static void pnm_skip_comments(std::istream& f); | |
| | | | |
| // | | // | |
| // cube saving | | // cube saving | |
| | | | |
| | | | |
End of changes. 5 change blocks. |
| 1 lines changed or deleted | | 5 lines changed or added | |
|
| diskio_meat.hpp | | diskio_meat.hpp | |
| | | | |
| skipping to change at line 318 | | skipping to change at line 318 | |
| return std::string("ARMA_CUB_BIN_FC016"); | | return std::string("ARMA_CUB_BIN_FC016"); | |
| } | | } | |
| else | | else | |
| { | | { | |
| return std::string(); | | return std::string(); | |
| } | | } | |
| | | | |
| } | | } | |
| | | | |
| inline | | inline | |
|
| bool | | file_type | |
| diskio::is_raw_binary(std::istream& f) | | diskio::guess_file_type(std::istream& f) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| f.clear(); | | f.clear(); | |
| const std::fstream::pos_type pos1 = f.tellg(); | | const std::fstream::pos_type pos1 = f.tellg(); | |
| | | | |
| f.clear(); | | f.clear(); | |
| f.seekg(0, ios::end); | | f.seekg(0, ios::end); | |
| | | | |
| f.clear(); | | f.clear(); | |
| | | | |
| skipping to change at line 344 | | skipping to change at line 344 | |
| f.clear(); | | f.clear(); | |
| f.seekg(pos1); | | f.seekg(pos1); | |
| | | | |
| podarray<unsigned char> data(N); | | podarray<unsigned char> data(N); | |
| | | | |
| unsigned char* ptr = data.memptr(); | | unsigned char* ptr = data.memptr(); | |
| | | | |
| f.clear(); | | f.clear(); | |
| f.read( reinterpret_cast<char*>(ptr), std::streamsize(N) ); | | f.read( reinterpret_cast<char*>(ptr), std::streamsize(N) ); | |
| | | | |
|
| bool has_non_text_val = false; | | const bool load_okay = f.good(); | |
| | | | |
| | | f.clear(); | |
| | | f.seekg(pos1); | |
| | | | |
|
| if(f.good() == true) | | bool has_binary = false; | |
| | | bool has_comma = false; | |
| | | | |
| | | if(load_okay == true) | |
| { | | { | |
|
| for(u32 i=0; i<N; ++i) | | u32 i = 0; | |
| | | u32 j = (N >= 2) ? 1 : 0; | |
| | | | |
| | | for(; j<N; i+=2, j+=2) | |
| { | | { | |
|
| const unsigned char val = ptr[i]; | | const unsigned char val_i = ptr[i]; | |
| | | const unsigned char val_j = ptr[j]; | |
| | | | |
| // the range checking can be made more elaborate | | // the range checking can be made more elaborate | |
|
| if( (val <= 8) || (val >= 123) ) | | if( ((val_i <= 8) || (val_i >= 123)) || ((val_j <= 8) || (val_j >= 12
3)) ) | |
| { | | { | |
|
| has_non_text_val = true; | | has_binary = true; | |
| | | break; | |
| | | } | |
| | | | |
| | | if( (val_i == ',') || (val_j == ',') ) | |
| | | { | |
| | | has_comma = true; | |
| break; | | break; | |
| } | | } | |
| } | | } | |
| } | | } | |
|
| | | else | |
| | | { | |
| | | return file_type_unknown; | |
| | | } | |
| | | | |
|
| f.clear(); | | if(has_binary) | |
| f.seekg(pos1); | | { | |
| | | return raw_binary; | |
| | | } | |
| | | | |
| | | if(has_comma) | |
| | | { | |
| | | return csv_ascii; | |
| | | } | |
| | | | |
|
| return has_non_text_val; | | return raw_ascii; | |
| } | | } | |
| | | | |
| inline | | inline | |
| char | | char | |
| diskio::conv_to_hex_char(const u8 x) | | diskio::conv_to_hex_char(const u8 x) | |
| { | | { | |
| char out; | | char out; | |
| | | | |
| switch(x) | | switch(x) | |
| { | | { | |
| | | | |
| skipping to change at line 690 | | skipping to change at line 717 | |
| f.put('\n'); | | f.put('\n'); | |
| } | | } | |
| | | | |
| const bool save_okay = f.good(); | | const bool save_okay = f.good(); | |
| | | | |
| f.flags(orig_flags); | | f.flags(orig_flags); | |
| | | | |
| return save_okay; | | return save_okay; | |
| } | | } | |
| | | | |
|
| | | //! Save a matrix in CSV text format (human readable) | |
| | | template<typename eT> | |
| | | inline | |
| | | bool | |
| | | diskio::save_csv_ascii(const Mat<eT>& x, const std::string& final_name) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const std::string tmp_name = diskio::gen_tmp_name(final_name); | |
| | | | |
| | | std::ofstream f(tmp_name.c_str()); | |
| | | | |
| | | bool save_okay = f.is_open(); | |
| | | | |
| | | if(save_okay == true) | |
| | | { | |
| | | save_okay = diskio::save_csv_ascii(x, f); | |
| | | | |
| | | f.flush(); | |
| | | f.close(); | |
| | | | |
| | | if(save_okay == true) | |
| | | { | |
| | | save_okay = diskio::safe_rename(tmp_name, final_name); | |
| | | } | |
| | | } | |
| | | | |
| | | return save_okay; | |
| | | } | |
| | | | |
| | | //! Save a matrix in CSV text format (human readable) | |
| | | template<typename eT> | |
| | | inline | |
| | | bool | |
| | | diskio::save_csv_ascii(const Mat<eT>& x, std::ostream& f) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | const ios::fmtflags orig_flags = f.flags(); | |
| | | | |
| | | // TODO: need sane values for complex numbers | |
| | | | |
| | | if( (is_float<eT>::value == true) || (is_double<eT>::value == true) ) | |
| | | { | |
| | | f.setf(ios::scientific); | |
| | | f.precision(10); | |
| | | } | |
| | | | |
| | | u32 x_n_rows = x.n_rows; | |
| | | u32 x_n_cols = x.n_cols; | |
| | | | |
| | | for(u32 row=0; row < x_n_rows; ++row) | |
| | | { | |
| | | for(u32 col=0; col < x_n_cols; ++col) | |
| | | { | |
| | | f << x.at(row,col); | |
| | | | |
| | | if( col < (x_n_cols-1) ) | |
| | | { | |
| | | f.put(','); | |
| | | } | |
| | | } | |
| | | | |
| | | f.put('\n'); | |
| | | } | |
| | | | |
| | | const bool save_okay = f.good(); | |
| | | | |
| | | f.flags(orig_flags); | |
| | | | |
| | | return save_okay; | |
| | | } | |
| | | | |
| //! Save a matrix in binary format, | | //! Save a matrix in binary format, | |
| //! with a header that stores the matrix type as well as its dimensions | | //! with a header that stores the matrix type as well as its dimensions | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| diskio::save_arma_binary(const Mat<eT>& x, const std::string& final_name) | | diskio::save_arma_binary(const Mat<eT>& x, const std::string& final_name) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const std::string tmp_name = diskio::gen_tmp_name(final_name); | | const std::string tmp_name = diskio::gen_tmp_name(final_name); | |
| | | | |
| skipping to change at line 864 | | skipping to change at line 964 | |
| //! Load a matrix as raw text (no header, human readable). | | //! Load a matrix as raw text (no header, human readable). | |
| //! Can read matrices saved as text in Matlab and Octave. | | //! Can read matrices saved as text in Matlab and Octave. | |
| //! NOTE: this is much slower than reading a file with a header. | | //! NOTE: this is much slower than reading a file with a header. | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| diskio::load_raw_ascii(Mat<eT>& x, std::istream& f, std::string& err_msg) | | diskio::load_raw_ascii(Mat<eT>& x, std::istream& f, std::string& err_msg) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| bool load_okay = true; | | bool load_okay = f.good(); | |
| | | | |
|
| //std::fstream::pos_type start = f.tellg(); | | f.clear(); | |
| | | const std::fstream::pos_type pos1 = f.tellg(); | |
| | | | |
| // | | // | |
| // work out the size | | // work out the size | |
| | | | |
| u32 f_n_rows = 0; | | u32 f_n_rows = 0; | |
| u32 f_n_cols = 0; | | u32 f_n_cols = 0; | |
| | | | |
| bool f_n_cols_found = false; | | bool f_n_cols_found = false; | |
| | | | |
| std::string line_string; | | std::string line_string; | |
| | | | |
| skipping to change at line 917 | | skipping to change at line 1018 | |
| load_okay = false; | | load_okay = false; | |
| } | | } | |
| } | | } | |
| | | | |
| ++f_n_rows; | | ++f_n_rows; | |
| } | | } | |
| | | | |
| if(load_okay == true) | | if(load_okay == true) | |
| { | | { | |
| f.clear(); | | f.clear(); | |
|
| f.seekg(0, ios::beg); | | f.seekg(pos1); | |
| //f.seekg(start); | | | |
| | | | |
| x.set_size(f_n_rows, f_n_cols); | | x.set_size(f_n_rows, f_n_cols); | |
| | | | |
| eT val; | | eT val; | |
| | | | |
| for(u32 row=0; (row < x.n_rows) && (load_okay == true); ++row) | | for(u32 row=0; (row < x.n_rows) && (load_okay == true); ++row) | |
| { | | { | |
| for(u32 col=0; (col < x.n_cols) && (load_okay == true); ++col) | | for(u32 col=0; (col < x.n_cols) && (load_okay == true); ++col) | |
| { | | { | |
| f >> val; | | f >> val; | |
| | | | |
| skipping to change at line 1071 | | skipping to change at line 1171 | |
| } | | } | |
| else | | else | |
| { | | { | |
| load_okay = false; | | load_okay = false; | |
| err_msg = "incorrect header in "; | | err_msg = "incorrect header in "; | |
| } | | } | |
| | | | |
| return load_okay; | | return load_okay; | |
| } | | } | |
| | | | |
|
| | | //! Load a matrix in CSV text format (human readable) | |
| | | template<typename eT> | |
| | | inline | |
| | | bool | |
| | | diskio::load_csv_ascii(Mat<eT>& x, const std::string& name, std::string& er | |
| | | r_msg) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | std::fstream f; | |
| | | f.open(name.c_str(), std::fstream::in); | |
| | | | |
| | | bool load_okay = f.is_open(); | |
| | | | |
| | | if(load_okay == true) | |
| | | { | |
| | | load_okay = diskio::load_csv_ascii(x, f, err_msg); | |
| | | f.close(); | |
| | | } | |
| | | | |
| | | return load_okay; | |
| | | } | |
| | | | |
| | | //! Load a matrix in CSV text format (human readable) | |
| | | template<typename eT> | |
| | | inline | |
| | | bool | |
| | | diskio::load_csv_ascii(Mat<eT>& x, std::istream& f, std::string& err_msg) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | bool load_okay = f.good(); | |
| | | | |
| | | f.clear(); | |
| | | const std::fstream::pos_type pos1 = f.tellg(); | |
| | | | |
| | | // | |
| | | // work out the size | |
| | | | |
| | | u32 f_n_rows = 0; | |
| | | u32 f_n_cols = 0; | |
| | | | |
| | | std::string line_string; | |
| | | std::string token; | |
| | | | |
| | | while( (f.good() == true) && (load_okay == true) ) | |
| | | { | |
| | | std::getline(f, line_string); | |
| | | | |
| | | if(line_string.size() == 0) | |
| | | { | |
| | | break; | |
| | | } | |
| | | | |
| | | std::stringstream line_stream(line_string); | |
| | | | |
| | | u32 line_n_cols = 0; | |
| | | | |
| | | while(line_stream.good() == true) | |
| | | { | |
| | | getline(line_stream, token, ','); | |
| | | ++line_n_cols; | |
| | | } | |
| | | | |
| | | if(f_n_cols < line_n_cols) | |
| | | { | |
| | | f_n_cols = line_n_cols; | |
| | | } | |
| | | | |
| | | ++f_n_rows; | |
| | | } | |
| | | | |
| | | f.clear(); | |
| | | f.seekg(pos1); | |
| | | | |
| | | x.zeros(f_n_rows, f_n_cols); | |
| | | | |
| | | u32 row = 0; | |
| | | | |
| | | while(f.good() == true) | |
| | | { | |
| | | std::getline(f, line_string); | |
| | | | |
| | | if(line_string.size() == 0) | |
| | | { | |
| | | break; | |
| | | } | |
| | | | |
| | | std::stringstream line_stream(line_string); | |
| | | | |
| | | u32 col = 0; | |
| | | | |
| | | while(line_stream.good() == true) | |
| | | { | |
| | | getline(line_stream, token, ','); | |
| | | | |
| | | eT val; | |
| | | | |
| | | std::stringstream ss(token); | |
| | | | |
| | | ss >> val; | |
| | | | |
| | | if(ss.fail() == false) | |
| | | { | |
| | | x.at(row,col) = val; | |
| | | } | |
| | | | |
| | | ++col; | |
| | | } | |
| | | | |
| | | ++row; | |
| | | } | |
| | | | |
| | | return load_okay; | |
| | | } | |
| | | | |
| //! Load a matrix in binary format, | | //! Load a matrix in binary format, | |
| //! with a header that indicates the matrix type as well as its dimensions | | //! with a header that indicates the matrix type as well as its dimensions | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| diskio::load_arma_binary(Mat<eT>& x, const std::string& name, std::string&
err_msg) | | diskio::load_arma_binary(Mat<eT>& x, const std::string& name, std::string&
err_msg) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| std::ifstream f; | | std::ifstream f; | |
| | | | |
| skipping to change at line 1362 | | skipping to change at line 1577 | |
| { | | { | |
| return load_arma_binary(x, f, err_msg); | | return load_arma_binary(x, f, err_msg); | |
| } | | } | |
| else | | else | |
| if(P5 == header.substr(0,P5.length())) | | if(P5 == header.substr(0,P5.length())) | |
| { | | { | |
| return load_pgm_binary(x, f, err_msg); | | return load_pgm_binary(x, f, err_msg); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| if(is_raw_binary(f) == true) | | const file_type ft = guess_file_type(f); | |
| { | | | |
| return load_raw_binary(x, f, err_msg); | | switch(ft) | |
| } | | | |
| else | | | |
| { | | { | |
|
| return load_raw_ascii(x, f, err_msg); | | case csv_ascii: | |
| | | return load_csv_ascii(x, f, err_msg); | |
| | | break; | |
| | | | |
| | | case raw_binary: | |
| | | return load_raw_binary(x, f, err_msg); | |
| | | break; | |
| | | | |
| | | case raw_ascii: | |
| | | return load_raw_ascii(x, f, err_msg); | |
| | | break; | |
| | | | |
| | | default: | |
| | | err_msg = "unknown data in "; | |
| | | return false; | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| // cubes | | // cubes | |
| | | | |
| //! Save a cube as raw text (no header, human readable). | | //! Save a cube as raw text (no header, human readable). | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| bool | | bool | |
| | | | |
| skipping to change at line 1924 | | skipping to change at line 2151 | |
| { | | { | |
| return load_arma_binary(x, f, err_msg); | | return load_arma_binary(x, f, err_msg); | |
| } | | } | |
| else | | else | |
| if(P6 == header.substr(0, P6.length())) | | if(P6 == header.substr(0, P6.length())) | |
| { | | { | |
| return load_ppm_binary(x, f, err_msg); | | return load_ppm_binary(x, f, err_msg); | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| if(is_raw_binary(f) == true) | | const file_type ft = guess_file_type(f); | |
| { | | | |
| return load_raw_binary(x, f, err_msg); | | switch(ft) | |
| } | | | |
| else | | | |
| { | | { | |
|
| return load_raw_ascii(x, f, err_msg); | | // case csv_ascii: | |
| | | // return load_csv_ascii(x, f, err_msg); | |
| | | // break; | |
| | | | |
| | | case raw_binary: | |
| | | return load_raw_binary(x, f, err_msg); | |
| | | break; | |
| | | | |
| | | case raw_ascii: | |
| | | return load_raw_ascii(x, f, err_msg); | |
| | | break; | |
| | | | |
| | | default: | |
| | | err_msg = "unknown data in "; | |
| | | return false; | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| // fields | | // fields | |
| | | | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| bool | | bool | |
| diskio::save_arma_binary(const field<T1>& x, const std::string& final_name) | | diskio::save_arma_binary(const field<T1>& x, const std::string& final_name) | |
| | | | |
End of changes. 19 change blocks. |
| 27 lines changed or deleted | | 267 lines changed or added | |
|
| fn_eig.hpp | | fn_eig.hpp | |
|
| // Copyright (C) 2008-2010 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2008-2010 Conrad Sanderson | | // Copyright (C) 2008-2011 Conrad Sanderson | |
| // Copyright (C) 2009 Edmund Highcock | | // Copyright (C) 2009 Edmund Highcock | |
| | | // Copyright (C) 2011 Stanislav Funiak | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| | | | |
| skipping to change at line 33 | | skipping to change at line 34 | |
| inline | | inline | |
| bool | | bool | |
| eig_sym | | eig_sym | |
| ( | | ( | |
| Col<typename T1::pod_type>& eigval, | | Col<typename T1::pod_type>& eigval, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| // unwrap_check not used as T1::elem_type and T1::pod_type may not be the
same. | | // unwrap_check not used as T1::elem_type and T1::pod_type may not be the
same. | |
| // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::
eig_sym() makes a copy of X | | // furthermore, it doesn't matter if X is an alias of eigval, as auxlib::
eig_sym() makes a copy of X | |
| | | | |
| const bool status = auxlib::eig_sym(eigval, X); | | const bool status = auxlib::eig_sym(eigval, X); | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
|
| arma_print("eig_sym(): failed to converge"); | | | |
| eigval.reset(); | | eigval.reset(); | |
|
| | | arma_bad("eig_sym(): failed to converge", false); | |
| } | | } | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| //! Eigenvalues of real/complex symmetric/hermitian matrix X | | //! Eigenvalues of real/complex symmetric/hermitian matrix X | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| Col<typename T1::pod_type> | | Col<typename T1::pod_type> | |
| eig_sym | | eig_sym | |
| ( | | ( | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| Col<typename T1::pod_type> out; | | Col<typename T1::pod_type> out; | |
|
| eig_sym(out, X); | | const bool status = auxlib::eig_sym(out, X); | |
| | | | |
| | | if(status == false) | |
| | | { | |
| | | out.reset(); | |
| | | arma_bad("eig_sym(): failed to converge"); | |
| | | } | |
| | | | |
| return out; | | return out; | |
| } | | } | |
| | | | |
| //! Eigenvalues and eigenvectors of real/complex symmetric/hermitian matrix
X | | //! Eigenvalues and eigenvectors of real/complex symmetric/hermitian matrix
X | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| bool | | bool | |
| eig_sym | | eig_sym | |
| ( | | ( | |
| Col<typename T1::pod_type>& eigval, | | Col<typename T1::pod_type>& eigval, | |
| Mat<typename T1::elem_type>& eigvec, | | Mat<typename T1::elem_type>& eigvec, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_sym(
): eigval is an alias of eigvec" ); | | arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_sym(
): eigval is an alias of eigvec" ); | |
| | | | |
| const bool status = auxlib::eig_sym(eigval, eigvec, X); | | const bool status = auxlib::eig_sym(eigval, eigvec, X); | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
|
| arma_print("eig_sym(): failed to converge"); | | | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
|
| | | arma_bad("eig_sym(): failed to converge", false); | |
| } | | } | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| // | | // | |
| // general matrices | | // general matrices | |
| // | | // | |
| | | | |
| //! Eigenvalues and eigenvectors (both left and right) of general real/comp
lex square matrix X | | //! Eigenvalues and eigenvectors (both left and right) of general real/comp
lex square matrix X | |
| | | | |
| skipping to change at line 118 | | skipping to change at line 122 | |
| eig_gen | | eig_gen | |
| ( | | ( | |
| Col< std::complex<typename T1::pod_type> >& eigval, | | Col< std::complex<typename T1::pod_type> >& eigval, | |
| Mat<typename T1::elem_type>& l_eigvec, | | Mat<typename T1::elem_type>& l_eigvec, | |
| Mat<typename T1::elem_type>& r_eigvec, | | Mat<typename T1::elem_type>& r_eigvec, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| ((&l_eigvec) == (&r_eigvec)), | | ((&l_eigvec) == (&r_eigvec)), | |
| "eig_gen(): l_eigvec is an alias of r_eigvec" | | "eig_gen(): l_eigvec is an alias of r_eigvec" | |
| ); | | ); | |
| | | | |
| arma_debug_check | | arma_debug_check | |
| ( | | ( | |
| | | | |
| skipping to change at line 141 | | skipping to change at line 144 | |
| || | | || | |
| (((void*)(&eigval)) == ((void*)(&r_eigvec))) | | (((void*)(&eigval)) == ((void*)(&r_eigvec))) | |
| ), | | ), | |
| "eig_gen(): eigval is an alias of l_eigvec or r_eigvec" | | "eig_gen(): eigval is an alias of l_eigvec or r_eigvec" | |
| ); | | ); | |
| | | | |
| const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b'); | | const bool status = auxlib::eig_gen(eigval, l_eigvec, r_eigvec, X, 'b'); | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
|
| arma_print("eig_gen(): failed to converge"); | | | |
| eigval.reset(); | | eigval.reset(); | |
| l_eigvec.reset(); | | l_eigvec.reset(); | |
| r_eigvec.reset(); | | r_eigvec.reset(); | |
|
| | | arma_bad("eig_gen(): failed to converge", false); | |
| } | | } | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| //! Eigenvalues and eigenvectors of general real square matrix X. | | //! Eigenvalues and eigenvectors of general real square matrix X. | |
| //! Optional argument 'side' specifies which eigenvectors should be compute
d: | | //! Optional argument 'side' specifies which eigenvectors should be compute
d: | |
| //! 'r' for right (default) and 'l' for left. | | //! 'r' for right (default) and 'l' for left. | |
| template<typename eT, typename T1> | | template<typename eT, typename T1> | |
| inline | | inline | |
| | | | |
| skipping to change at line 166 | | skipping to change at line 169 | |
| eig_gen | | eig_gen | |
| ( | | ( | |
| Col< std::complex<eT> >& eigval, | | Col< std::complex<eT> >& eigval, | |
| Mat< std::complex<eT> >& eigvec, | | Mat< std::complex<eT> >& eigvec, | |
| const Base<eT, T1>& X, | | const Base<eT, T1>& X, | |
| const char side = 'r', | | const char side = 'r', | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| //std::cout << "real" << std::endl; | | //std::cout << "real" << std::endl; | |
| | | | |
| arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(
): eigval is an alias of eigvec" ); | | arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(
): eigval is an alias of eigvec" ); | |
| | | | |
| Mat<eT> dummy_eigvec; | | Mat<eT> dummy_eigvec; | |
| Mat<eT> tmp_eigvec; | | Mat<eT> tmp_eigvec; | |
| | | | |
| bool status; | | bool status; | |
| | | | |
| skipping to change at line 195 | | skipping to change at line 197 | |
| status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side); | | status = auxlib::eig_gen(eigval, tmp_eigvec, dummy_eigvec, X, side); | |
| break; | | break; | |
| | | | |
| default: | | default: | |
| arma_stop("eig_gen(): parameter 'side' is invalid"); | | arma_stop("eig_gen(): parameter 'side' is invalid"); | |
| status = false; | | status = false; | |
| } | | } | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
|
| arma_print("eig_gen(): failed to converge"); | | | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
|
| return false; | | arma_bad("eig_gen(): failed to converge", false); | |
| } | | } | |
|
| | | else | |
| const u32 n = eigval.n_elem; | | | |
| | | | |
| if(n > 0) | | | |
| { | | { | |
|
| eigvec.set_size(n,n); | | const u32 n = eigval.n_elem; | |
| | | | |
|
| for(u32 j=0; j<n; ++j) | | if(n > 0) | |
| { | | { | |
|
| if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) ) | | eigvec.set_size(n,n); | |
| { | | | |
| // eigvec.col(j) = Mat< std::complex<eT> >( tmp_eigvec.col(j), t | | | |
| mp_eigvec.col(j+1) ); | | | |
| // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), -t | | | |
| mp_eigvec.col(j+1) ); | | | |
| | | | |
|
| for(u32 i=0; i<n; ++i) | | for(u32 j=0; j<n; ++j) | |
| | | { | |
| | | if( (j < n-1) && (eigval[j] == std::conj(eigval[j+1])) ) | |
| { | | { | |
|
| eigvec.at(i,j) = std::complex<eT>( tmp_eigvec.at(i,j), tmp_eig | | // eigvec.col(j) = Mat< std::complex<eT> >( tmp_eigvec.col(j), | |
| vec.at(i,j+1) ); | | tmp_eigvec.col(j+1) ); | |
| eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_eig | | // eigvec.col(j+1) = Mat< std::complex<eT> >( tmp_eigvec.col(j), | |
| vec.at(i,j+1) ); | | -tmp_eigvec.col(j+1) ); | |
| } | | | |
| | | | |
|
| ++j; | | for(u32 i=0; i<n; ++i) | |
| } | | { | |
| else | | eigvec.at(i,j) = std::complex<eT>( tmp_eigvec.at(i,j), tmp_e | |
| { | | igvec.at(i,j+1) ); | |
| // eigvec.col(i) = tmp_eigvec.col(i); | | eigvec.at(i,j+1) = std::complex<eT>( tmp_eigvec.at(i,j), -tmp_e | |
| | | igvec.at(i,j+1) ); | |
| | | } | |
| | | | |
|
| for(u32 i=0; i<n; ++i) | | ++j; | |
| { | | | |
| eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0)); | | | |
| } | | } | |
|
| | | else | |
| | | { | |
| | | // eigvec.col(i) = tmp_eigvec.col(i); | |
| | | | |
| | | for(u32 i=0; i<n; ++i) | |
| | | { | |
| | | eigvec.at(i,j) = std::complex<eT>(tmp_eigvec.at(i,j), eT(0)); | |
| | | } | |
| | | | |
|
| | | } | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
|
| return true; | | return status; | |
| } | | } | |
| | | | |
| //! Eigenvalues and eigenvectors of general complex square matrix X | | //! Eigenvalues and eigenvectors of general complex square matrix X | |
| //! Optional argument 'side' specifies which eigenvectors should be compute
d: | | //! Optional argument 'side' specifies which eigenvectors should be compute
d: | |
| //! 'r' for right (default) and 'l' for left. | | //! 'r' for right (default) and 'l' for left. | |
| template<typename T, typename T1> | | template<typename T, typename T1> | |
| inline | | inline | |
| bool | | bool | |
| eig_gen | | eig_gen | |
| ( | | ( | |
| Col<std::complex<T> >& eigval, | | Col<std::complex<T> >& eigval, | |
| Mat<std::complex<T> >& eigvec, | | Mat<std::complex<T> >& eigvec, | |
| const Base<std::complex<T>, T1>& X, | | const Base<std::complex<T>, T1>& X, | |
| const char side = 'r', | | const char side = 'r', | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
|
| | | | |
| arma_ignore(junk); | | arma_ignore(junk); | |
| | | | |
| //std::cout << "complex" << std::endl; | | //std::cout << "complex" << std::endl; | |
| | | | |
| arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(
): eigval is an alias of eigvec" ); | | arma_debug_check( ( ((void*)(&eigval)) == ((void*)(&eigvec)) ), "eig_gen(
): eigval is an alias of eigvec" ); | |
| | | | |
| Mat< std::complex<T> > dummy_eigvec; | | Mat< std::complex<T> > dummy_eigvec; | |
| | | | |
| bool status; | | bool status; | |
| | | | |
| | | | |
| skipping to change at line 282 | | skipping to change at line 284 | |
| status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side); | | status = auxlib::eig_gen(eigval, eigvec, dummy_eigvec, X, side); | |
| break; | | break; | |
| | | | |
| default: | | default: | |
| arma_stop("eig_gen(): parameter 'side' is invalid"); | | arma_stop("eig_gen(): parameter 'side' is invalid"); | |
| status = false; | | status = false; | |
| } | | } | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
|
| arma_print("eig_gen(): failed to converge"); | | | |
| eigval.reset(); | | eigval.reset(); | |
| eigvec.reset(); | | eigvec.reset(); | |
|
| | | arma_bad("eig_gen(): failed to converge", false); | |
| } | | } | |
| | | | |
| return status; | | return status; | |
| } | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 29 change blocks. |
| 43 lines changed or deleted | | 45 lines changed or added | |
|
| fn_princomp.hpp | | fn_princomp.hpp | |
|
| // Copyright (C) 2010 NICTA (www.nicta.com.au) | | // Copyright (C) 2010-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2010 Conrad Sanderson | | // Copyright (C) 2010-2011 Conrad Sanderson | |
| // Copyright (C) 2010 Dimitrios Bouzas | | // Copyright (C) 2010 Dimitrios Bouzas | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| skipping to change at line 25 | | skipping to change at line 25 | |
| //! @{ | | //! @{ | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 4 arguments version | | //! principal component analysis -- 4 arguments version | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| //! tsquared_out -> Hotelling's T^2 statistic | | //! tsquared_out -> Hotelling's T^2 statistic | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| princomp | | princomp | |
| ( | | ( | |
| Mat<typename T1::elem_type>& coeff_out, | | Mat<typename T1::elem_type>& coeff_out, | |
| Mat<typename T1::elem_type>& score_out, | | Mat<typename T1::elem_type>& score_out, | |
| Col<typename T1::pod_type>& latent_out, | | Col<typename T1::pod_type>& latent_out, | |
| Col<typename T1::elem_type>& tsquared_out, | | Col<typename T1::elem_type>& tsquared_out, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; | |
| | | | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
|
| op_princomp::direct_princomp(coeff_out, score_out, latent_out, tsquared_o | | const bool status = op_princomp::direct_princomp(coeff_out, score_out, la | |
| ut, A); | | tent_out, tsquared_out, A); | |
| | | | |
| | | if(status == false) | |
| | | { | |
| | | coeff_out.reset(); | |
| | | score_out.reset(); | |
| | | latent_out.reset(); | |
| | | tsquared_out.reset(); | |
| | | | |
| | | arma_bad("princomp(): failed to converge", false); | |
| | | } | |
| | | | |
| | | return status; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 3 arguments version | | //! principal component analysis -- 3 arguments version | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| princomp | | princomp | |
| ( | | ( | |
| Mat<typename T1::elem_type>& coeff_out, | | Mat<typename T1::elem_type>& coeff_out, | |
| Mat<typename T1::elem_type>& score_out, | | Mat<typename T1::elem_type>& score_out, | |
| Col<typename T1::pod_type>& latent_out, | | Col<typename T1::pod_type>& latent_out, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; | |
| | | | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
|
| op_princomp::direct_princomp(coeff_out, score_out, latent_out, A); | | const bool status = op_princomp::direct_princomp(coeff_out, score_out, la | |
| | | tent_out, A); | |
| | | | |
| | | if(status == false) | |
| | | { | |
| | | coeff_out.reset(); | |
| | | score_out.reset(); | |
| | | latent_out.reset(); | |
| | | | |
| | | arma_bad("princomp(): failed to converge", false); | |
| | | } | |
| | | | |
| | | return status; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 2 arguments version | | //! principal component analysis -- 2 arguments version | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
|
| void | | bool | |
| princomp | | princomp | |
| ( | | ( | |
| Mat<typename T1::elem_type>& coeff_out, | | Mat<typename T1::elem_type>& coeff_out, | |
| Mat<typename T1::elem_type>& score_out, | | Mat<typename T1::elem_type>& score_out, | |
| const Base<typename T1::elem_type,T1>& X, | | const Base<typename T1::elem_type,T1>& X, | |
| const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | | const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0 | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; | |
| | | | |
| const unwrap<T1> tmp(X.get_ref()); | | const unwrap<T1> tmp(X.get_ref()); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
|
| op_princomp::direct_princomp(coeff_out, score_out, A); | | const bool status = op_princomp::direct_princomp(coeff_out, score_out, A) | |
| | | ; | |
| | | | |
| | | if(status == false) | |
| | | { | |
| | | coeff_out.reset(); | |
| | | score_out.reset(); | |
| | | | |
| | | arma_bad("princomp(): failed to converge", false); | |
| | | } | |
| | | | |
| | | return status; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 1 argument version | | //! principal component analysis -- 1 argument version | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| const Op<T1, op_princomp> | | const Op<T1, op_princomp> | |
| princomp | | princomp | |
| ( | | ( | |
| | | | |
End of changes. 7 change blocks. |
| 9 lines changed or deleted | | 44 lines changed or added | |
|
| lapack_bones.hpp | | lapack_bones.hpp | |
| | | | |
| skipping to change at line 98 | | skipping to change at line 98 | |
| #define arma_sgees sgees | | #define arma_sgees sgees | |
| #define arma_dgees dgees | | #define arma_dgees dgees | |
| #define arma_cgees cgees | | #define arma_cgees cgees | |
| #define arma_zgees zgees | | #define arma_zgees zgees | |
| | | | |
| #define arma_strsyl strsyl | | #define arma_strsyl strsyl | |
| #define arma_dtrsyl dtrsyl | | #define arma_dtrsyl dtrsyl | |
| #define arma_ctrsyl ctrsyl | | #define arma_ctrsyl ctrsyl | |
| #define arma_ztrsyl ztrsyl | | #define arma_ztrsyl ztrsyl | |
| | | | |
|
| | | #define arma_ssytrf ssytrf | |
| | | #define arma_dsytrf dsytrf | |
| | | #define arma_csytrf csytrf | |
| | | #define arma_zsytrf zsytrf | |
| | | | |
| | | #define arma_ssytri ssytri | |
| | | #define arma_dsytri dsytri | |
| | | #define arma_csytri csytri | |
| | | #define arma_zsytri zsytri | |
| | | | |
| #else | | #else | |
| | | | |
| #define arma_sgetrf SGETRF | | #define arma_sgetrf SGETRF | |
| #define arma_dgetrf DGETRF | | #define arma_dgetrf DGETRF | |
| #define arma_cgetrf CGETRF | | #define arma_cgetrf CGETRF | |
| #define arma_zgetrf ZGETRF | | #define arma_zgetrf ZGETRF | |
| | | | |
| #define arma_sgetri SGETRI | | #define arma_sgetri SGETRI | |
| #define arma_dgetri DGETRI | | #define arma_dgetri DGETRI | |
| #define arma_cgetri CGETRI | | #define arma_cgetri CGETRI | |
| | | | |
| skipping to change at line 179 | | skipping to change at line 189 | |
| #define arma_sgees SGEES | | #define arma_sgees SGEES | |
| #define arma_dgees DGEES | | #define arma_dgees DGEES | |
| #define arma_cgees CGEES | | #define arma_cgees CGEES | |
| #define arma_zgees ZGEES | | #define arma_zgees ZGEES | |
| | | | |
| #define arma_strsyl STRSYL | | #define arma_strsyl STRSYL | |
| #define arma_dtrsyl DTRSYL | | #define arma_dtrsyl DTRSYL | |
| #define arma_ctrsyl CTRSYL | | #define arma_ctrsyl CTRSYL | |
| #define arma_ztrsyl ZTRSYL | | #define arma_ztrsyl ZTRSYL | |
| | | | |
|
| | | #define arma_ssytrf SSYTRF | |
| | | #define arma_dsytrf DSYTRF | |
| | | #define arma_csytrf CSYTRF | |
| | | #define arma_zsytrf ZSYTRF | |
| | | | |
| | | #define arma_ssytri SSYTRI | |
| | | #define arma_dsytri DSYTRI | |
| | | #define arma_csytri CSYTRI | |
| | | #define arma_zsytri ZSYTRI | |
| | | | |
| #endif | | #endif | |
| | | | |
| //! \namespace lapack namespace for LAPACK functions | | //! \namespace lapack namespace for LAPACK functions | |
| namespace lapack | | namespace lapack | |
| { | | { | |
| //! \addtogroup LAPACK | | //! \addtogroup LAPACK | |
| //! @{ | | //! @{ | |
| | | | |
| extern "C" | | extern "C" | |
| { | | { | |
| | | | |
| skipping to change at line 289 | | skipping to change at line 309 | |
| // Schur decomposition (complex matrices) | | // Schur decomposition (complex matrices) | |
| void arma_fortran(arma_cgees)(char* jobvs, char* sort, blas_int* select
, blas_int* n, void* a, blas_int* lda, blas_int* sdim, void* w, void* vs, b
las_int* ldvs, void* work, blas_int* lwork, float* rwork, blas_int* bwork,
blas_int* info); | | void arma_fortran(arma_cgees)(char* jobvs, char* sort, blas_int* select
, blas_int* n, void* a, blas_int* lda, blas_int* sdim, void* w, void* vs, b
las_int* ldvs, void* work, blas_int* lwork, float* rwork, blas_int* bwork,
blas_int* info); | |
| void arma_fortran(arma_zgees)(char* jobvs, char* sort, blas_int* select
, blas_int* n, void* a, blas_int* lda, blas_int* sdim, void* w, void* vs, b
las_int* ldvs, void* work, blas_int* lwork, double* rwork, blas_int* bwork,
blas_int* info); | | void arma_fortran(arma_zgees)(char* jobvs, char* sort, blas_int* select
, blas_int* n, void* a, blas_int* lda, blas_int* sdim, void* w, void* vs, b
las_int* ldvs, void* work, blas_int* lwork, double* rwork, blas_int* bwork,
blas_int* info); | |
| | | | |
| // solve a Sylvester equation ax + xb = c, with a and b assumed to be i
n Schur form | | // solve a Sylvester equation ax + xb = c, with a and b assumed to be i
n Schur form | |
| void arma_fortran(arma_strsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const float* a, blas_int* lda, const float*
b, blas_int* ldb, float* c, blas_int* ldc, float* scale, blas_int* info)
; | | void arma_fortran(arma_strsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const float* a, blas_int* lda, const float*
b, blas_int* ldb, float* c, blas_int* ldc, float* scale, blas_int* info)
; | |
| void arma_fortran(arma_dtrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const double* a, blas_int* lda, const double*
b, blas_int* ldb, double* c, blas_int* ldc, double* scale, blas_int* info)
; | | void arma_fortran(arma_dtrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const double* a, blas_int* lda, const double*
b, blas_int* ldb, double* c, blas_int* ldc, double* scale, blas_int* info)
; | |
| void arma_fortran(arma_ctrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const void* a, blas_int* lda, const void*
b, blas_int* ldb, void* c, blas_int* ldc, float* scale, blas_int* info)
; | | void arma_fortran(arma_ctrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const void* a, blas_int* lda, const void*
b, blas_int* ldb, void* c, blas_int* ldc, float* scale, blas_int* info)
; | |
| void arma_fortran(arma_ztrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const void* a, blas_int* lda, const void*
b, blas_int* ldb, void* c, blas_int* ldc, double* scale, blas_int* info)
; | | void arma_fortran(arma_ztrsyl)(char* transa, char* transb, blas_int* is
gn, blas_int* m, blas_int* n, const void* a, blas_int* lda, const void*
b, blas_int* ldb, void* c, blas_int* ldc, double* scale, blas_int* info)
; | |
| | | | |
|
| | | void arma_fortran(arma_ssytrf)(char* uplo, blas_int* n, float* a, blas | |
| | | _int* lda, blas_int* ipiv, float* work, blas_int* lwork, blas_int* info); | |
| | | void arma_fortran(arma_dsytrf)(char* uplo, blas_int* n, double* a, blas | |
| | | _int* lda, blas_int* ipiv, double* work, blas_int* lwork, blas_int* info); | |
| | | void arma_fortran(arma_csytrf)(char* uplo, blas_int* n, void* a, blas | |
| | | _int* lda, blas_int* ipiv, void* work, blas_int* lwork, blas_int* info); | |
| | | void arma_fortran(arma_zsytrf)(char* uplo, blas_int* n, void* a, blas | |
| | | _int* lda, blas_int* ipiv, void* work, blas_int* lwork, blas_int* info); | |
| | | | |
| | | void arma_fortran(arma_ssytri)(char* uplo, blas_int* n, float* a, blas | |
| | | _int* lda, blas_int* ipiv, float* work, blas_int* info); | |
| | | void arma_fortran(arma_dsytri)(char* uplo, blas_int* n, double* a, blas | |
| | | _int* lda, blas_int* ipiv, double* work, blas_int* info); | |
| | | void arma_fortran(arma_csytri)(char* uplo, blas_int* n, void* a, blas | |
| | | _int* lda, blas_int* ipiv, void* work, blas_int* info); | |
| | | void arma_fortran(arma_zsytri)(char* uplo, blas_int* n, void* a, blas | |
| | | _int* lda, blas_int* ipiv, void* work, blas_int* info); | |
| | | | |
| // void arma_fortran(arma_dgeqp3)(blas_int* m, blas_int* n, double* a,
blas_int* lda, blas_int* jpvt, double* tau, double* work, blas_int* lwork,
blas_int* info); | | // void arma_fortran(arma_dgeqp3)(blas_int* m, blas_int* n, double* a,
blas_int* lda, blas_int* jpvt, double* tau, double* work, blas_int* lwork,
blas_int* info); | |
| // void arma_fortran(arma_dormqr)(char* side, char* trans, blas_int* m,
blas_int* n, blas_int* k, double* a, blas_int* lda, double* tau, double* c
, blas_int* ldc, double* work, blas_int* lwork, blas_int* info); | | // void arma_fortran(arma_dormqr)(char* side, char* trans, blas_int* m,
blas_int* n, blas_int* k, double* a, blas_int* lda, double* tau, double* c
, blas_int* ldc, double* work, blas_int* lwork, blas_int* info); | |
| // void arma_fortran(arma_dposv)(char* uplo, blas_int* n, blas_int* nr
hs, double* a, blas_int* lda, double* b, blas_int* ldb, blas_int* info); | | // void arma_fortran(arma_dposv)(char* uplo, blas_int* n, blas_int* nr
hs, double* a, blas_int* lda, double* b, blas_int* ldb, blas_int* info); | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| void | | void | |
| getrf(blas_int* m, blas_int* n, eT* a, blas_int* lda, blas_int* ipiv, bla
s_int* info) | | getrf(blas_int* m, blas_int* n, eT* a, blas_int* lda, blas_int* ipiv, bla
s_int* info) | |
| { | | { | |
| | | | |
| skipping to change at line 866 | | skipping to change at line 896 | |
| arma_fortran(ctrsyl)(transa, transb, isgn, m, n, (T*)a, lda, (T*)b, l
db, (T*)c, ldc, (float*)scale, info); | | arma_fortran(ctrsyl)(transa, transb, isgn, m, n, (T*)a, lda, (T*)b, l
db, (T*)c, ldc, (float*)scale, info); | |
| } | | } | |
| else | | else | |
| if(is_supported_complex_double<eT>::value == true) | | if(is_supported_complex_double<eT>::value == true) | |
| { | | { | |
| typedef std::complex<double> T; | | typedef std::complex<double> T; | |
| arma_fortran(ztrsyl)(transa, transb, isgn, m, n, (T*)a, lda, (T*)b, l
db, (T*)c, ldc, (double*)scale, info); | | arma_fortran(ztrsyl)(transa, transb, isgn, m, n, (T*)a, lda, (T*)b, l
db, (T*)c, ldc, (double*)scale, info); | |
| } | | } | |
| } | | } | |
| | | | |
|
| | | template<typename eT> | |
| | | inline | |
| | | void | |
| | | sytrf(char* uplo, blas_int* n, eT* a, blas_int* lda, blas_int* ipiv, eT* | |
| | | work, blas_int* lwork, blas_int* info) | |
| | | { | |
| | | arma_type_check<is_supported_blas_type<eT>::value == false>::apply(); | |
| | | | |
| | | if(is_float<eT>::value == true) | |
| | | { | |
| | | typedef float T; | |
| | | arma_fortran(arma_ssytrf)(uplo, n, (T*)a, lda, ipiv, (T*)work, lwork, | |
| | | info); | |
| | | } | |
| | | else | |
| | | if(is_double<eT>::value == true) | |
| | | { | |
| | | typedef double T; | |
| | | arma_fortran(arma_dsytrf)(uplo, n, (T*)a, lda, ipiv, (T*)work, lwork, | |
| | | info); | |
| | | } | |
| | | else | |
| | | if(is_supported_complex_float<eT>::value == true) | |
| | | { | |
| | | typedef std::complex<float> T; | |
| | | arma_fortran(arma_csytrf)(uplo, n, (T*)a, lda, ipiv, (T*)work, lwork, | |
| | | info); | |
| | | } | |
| | | else | |
| | | if(is_supported_complex_double<eT>::value == true) | |
| | | { | |
| | | typedef std::complex<double> T; | |
| | | arma_fortran(arma_zsytrf)(uplo, n, (T*)a, lda, ipiv, (T*)work, lwork, | |
| | | info); | |
| | | } | |
| | | } | |
| | | | |
| | | template<typename eT> | |
| | | inline | |
| | | void | |
| | | sytri(char* uplo, blas_int* n, eT* a, blas_int* lda, blas_int* ipiv, eT* | |
| | | work, blas_int* info) | |
| | | { | |
| | | arma_type_check<is_supported_blas_type<eT>::value == false>::apply(); | |
| | | | |
| | | if(is_float<eT>::value == true) | |
| | | { | |
| | | typedef float T; | |
| | | arma_fortran(arma_ssytri)(uplo, n, (T*)a, lda, ipiv, (T*)work, info); | |
| | | } | |
| | | else | |
| | | if(is_double<eT>::value == true) | |
| | | { | |
| | | typedef double T; | |
| | | arma_fortran(arma_dsytri)(uplo, n, (T*)a, lda, ipiv, (T*)work, info); | |
| | | } | |
| | | else | |
| | | if(is_supported_complex_float<eT>::value == true) | |
| | | { | |
| | | typedef std::complex<float> T; | |
| | | arma_fortran(arma_csytri)(uplo, n, (T*)a, lda, ipiv, (T*)work, info); | |
| | | } | |
| | | else | |
| | | if(is_supported_complex_double<eT>::value == true) | |
| | | { | |
| | | typedef std::complex<double> T; | |
| | | arma_fortran(arma_zsytri)(uplo, n, (T*)a, lda, ipiv, (T*)work, info); | |
| | | } | |
| | | } | |
| | | | |
| //! @} | | //! @} | |
| } | | } | |
| | | | |
| #endif | | #endif | |
| | | | |
End of changes. 4 change blocks. |
| 0 lines changed or deleted | | 108 lines changed or added | |
|
| op_max_meat.hpp | | op_max_meat.hpp | |
| | | | |
| skipping to change at line 24 | | skipping to change at line 24 | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_pure | | arma_pure | |
| inline | | inline | |
| eT | | eT | |
| op_max::direct_max(const eT* const X, const u32 n_elem) | | op_max::direct_max(const eT* const X, const u32 n_elem) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| eT max_val = (n_elem != 1) ? priv::most_neg<eT>() : X[0]; | | eT max_val = priv::most_neg<eT>(); | |
| | | | |
| u32 i,j; | | u32 i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| const eT X_i = X[i]; | | const eT X_i = X[i]; | |
| const eT X_j = X[j]; | | const eT X_j = X[j]; | |
| | | | |
| if(X_i > max_val) | | if(X_i > max_val) | |
| { | | { | |
| | | | |
| skipping to change at line 64 | | skipping to change at line 64 | |
| return max_val; | | return max_val; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_max::direct_max(const eT* const X, const u32 n_elem, u32& index_of_max_v
al) | | op_max::direct_max(const eT* const X, const u32 n_elem, u32& index_of_max_v
al) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| eT max_val = (n_elem != 1) ? priv::most_neg<eT>() : X[0]; | | eT max_val = priv::most_neg<eT>(); | |
| | | | |
| u32 best_index = 0; | | u32 best_index = 0; | |
| | | | |
| u32 i,j; | | u32 i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| const eT X_i = X[i]; | | const eT X_i = X[i]; | |
| const eT X_j = X[j]; | | const eT X_j = X[j]; | |
| | | | |
| | | | |
| skipping to change at line 113 | | skipping to change at line 113 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_max::direct_max(const Mat<eT>& X, const u32 row) | | op_max::direct_max(const Mat<eT>& X, const u32 row) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
|
| eT max_val = (X_n_cols != 1) ? priv::most_neg<eT>() : X.at(row,0); | | eT max_val = priv::most_neg<eT>(); | |
| | | | |
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| const eT tmp_val = X.at(row,col); | | const eT tmp_val = X.at(row,col); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 137 | | skipping to change at line 137 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_max::direct_max(const subview<eT>& X) | | op_max::direct_max(const subview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
|
| eT max_val = (X_n_elem != 1) ? priv::most_neg<eT>() : X[0]; | | eT max_val = priv::most_neg<eT>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| eT tmp_val = X[i]; | | eT tmp_val = X[i]; | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 161 | | skipping to change at line 161 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_max::direct_max(const diagview<eT>& X) | | op_max::direct_max(const diagview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
|
| eT max_val = (X_n_elem != 1) ? priv::most_neg<eT>() : X[0]; | | eT max_val = priv::most_neg<eT>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| eT tmp_val = X[i]; | | eT tmp_val = X[i]; | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 202 | | skipping to change at line 202 | |
| const u32 dim = in.aux_u32_a; | | const u32 dim = in.aux_u32_a; | |
| arma_debug_check( (dim > 1), "max(): incorrect usage. dim must be 0 or 1"
); | | arma_debug_check( (dim > 1), "max(): incorrect usage. dim must be 0 or 1"
); | |
| | | | |
| const u32 X_n_rows = X.n_rows; | | const u32 X_n_rows = X.n_rows; | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| if(dim == 0) | | if(dim == 0) | |
| { | | { | |
| arma_extra_debug_print("op_max::apply(), dim = 0"); | | arma_extra_debug_print("op_max::apply(), dim = 0"); | |
| | | | |
|
| out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); | | arma_debug_check( (X_n_rows == 0), "max(): given object has zero rows"
); | |
| | | | |
|
| if(X_n_rows > 0) | | out.set_size(1, X_n_cols); | |
| { | | | |
| eT* out_mem = out.memptr(); | | eT* out_mem = out.memptr(); | |
| | | | |
|
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| out_mem[col] = op_max::direct_max( X.colptr(col), X_n_rows ); | | out_mem[col] = op_max::direct_max( X.colptr(col), X_n_rows ); | |
| } | | | |
| } | | } | |
| } | | } | |
| else | | else | |
| if(dim == 1) | | if(dim == 1) | |
| { | | { | |
| arma_extra_debug_print("op_max::apply(), dim = 1"); | | arma_extra_debug_print("op_max::apply(), dim = 1"); | |
| | | | |
|
| out.set_size( X_n_rows, (X_n_cols > 0) ? 1 : 0 ); | | arma_debug_check( (X_n_cols == 0), "max(): given object has zero column
s" ); | |
| | | | |
|
| if(X_n_cols > 0) | | out.set_size(X_n_rows, 1); | |
| { | | | |
| eT* out_mem = out.memptr(); | | eT* out_mem = out.memptr(); | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | for(u32 row=0; row<X_n_rows; ++row) | |
| { | | { | |
| out_mem[row] = op_max::direct_max( X, row ); | | out_mem[row] = op_max::direct_max( X, row ); | |
| } | | | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_max::direct_max(const std::complex<T>* const X, const u32 n_elem) | | op_max::direct_max(const std::complex<T>* const X, const u32 n_elem) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T max_val = (n_elem != 1) ? priv::most_neg<T>() : std::abs(X[0]); | | T max_val = priv::most_neg<T>(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 265 | | skipping to change at line 263 | |
| } | | } | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_max::direct_max(const std::complex<T>* const X, const u32 n_elem, u32& i
ndex_of_max_val) | | op_max::direct_max(const std::complex<T>* const X, const u32 n_elem, u32& i
ndex_of_max_val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T max_val = (n_elem != 1) ? priv::most_neg<T>() : std::abs(X[0]); | | T max_val = priv::most_neg<T>(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 293 | | skipping to change at line 291 | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_max::direct_max(const Mat< std::complex<T> >& X, const u32 row) | | op_max::direct_max(const Mat< std::complex<T> >& X, const u32 row) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T max_val = (X_n_cols != 1) ? priv::most_neg<T>() : std::abs(X.at(row,0
)); | | T max_val = priv::most_neg<T>(); | |
| | | | |
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| const T tmp_val = std::abs(X.at(row,col)); | | const T tmp_val = std::abs(X.at(row,col)); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| index = col; | | index = col; | |
| } | | } | |
| | | | |
| skipping to change at line 319 | | skipping to change at line 317 | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_max::direct_max(const subview< std::complex<T> >& X) | | op_max::direct_max(const subview< std::complex<T> >& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T max_val = (X_n_elem != 1) ? priv::most_neg<T>() : std::abs(X[0]); | | T max_val = priv::most_neg<T>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 345 | | skipping to change at line 343 | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_max::direct_max(const diagview< std::complex<T> >& X) | | op_max::direct_max(const diagview< std::complex<T> >& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T max_val = (X_n_elem != 1) ? priv::most_neg<T>() : std::abs(X[0]); | | T max_val = priv::most_neg<T>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val > max_val) | | if(tmp_val > max_val) | |
| { | | { | |
| max_val = tmp_val; | | max_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
End of changes. 16 change blocks. |
| 26 lines changed or deleted | | 24 lines changed or added | |
|
| op_median_meat.hpp | | op_median_meat.hpp | |
| | | | |
| skipping to change at line 123 | | skipping to change at line 123 | |
| const u32 X_n_rows = X.n_rows; | | const u32 X_n_rows = X.n_rows; | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| const u32 dim = in.aux_u32_a; | | const u32 dim = in.aux_u32_a; | |
| arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or
1"); | | arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or
1"); | |
| | | | |
| if(dim == 0) // in each column | | if(dim == 0) // in each column | |
| { | | { | |
| arma_extra_debug_print("op_median::apply(), dim = 0"); | | arma_extra_debug_print("op_median::apply(), dim = 0"); | |
| | | | |
|
| out.set_size((X_n_rows > 0) ? 1 : 0, X_n_cols); | | arma_debug_check( (X_n_rows == 0), "median(): given object has zero row
s" ); | |
| | | | |
|
| if(X_n_rows > 0) | | out.set_size(1, X_n_cols); | |
| { | | | |
| std::vector<eT> tmp_vec(X_n_rows); | | | |
| | | | |
|
| for(u32 col=0; col<X_n_cols; ++col) | | std::vector<eT> tmp_vec(X_n_rows); | |
| { | | | |
| const eT* colmem = X.colptr(col); | | | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| tmp_vec[row] = colmem[row]; | | const eT* colmem = X.colptr(col); | |
| } | | | |
| | | | |
|
| out[col] = op_median::direct_median(tmp_vec); | | for(u32 row=0; row<X_n_rows; ++row) | |
| | | { | |
| | | tmp_vec[row] = colmem[row]; | |
| } | | } | |
|
| | | | |
| | | out[col] = op_median::direct_median(tmp_vec); | |
| } | | } | |
| } | | } | |
| else | | else | |
| if(dim == 1) // in each row | | if(dim == 1) // in each row | |
| { | | { | |
| arma_extra_debug_print("op_median::apply(), dim = 1"); | | arma_extra_debug_print("op_median::apply(), dim = 1"); | |
| | | | |
|
| out.set_size(X_n_rows, (X_n_cols > 0) ? 1 : 0); | | arma_debug_check( (X_n_cols == 0), "median(): given object has zero col
umns" ); | |
| | | | |
|
| if(X_n_cols > 0) | | out.set_size(X_n_rows, 1); | |
| { | | | |
| std::vector<eT> tmp_vec(X_n_cols); | | | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | std::vector<eT> tmp_vec(X_n_cols); | |
| { | | | |
| for(u32 col=0; col<X_n_cols; ++col) | | | |
| { | | | |
| tmp_vec[col] = X.at(row,col); | | | |
| } | | | |
| | | | |
|
| out[row] = op_median::direct_median(tmp_vec); | | for(u32 row=0; row<X_n_rows; ++row) | |
| | | { | |
| | | for(u32 col=0; col<X_n_cols; ++col) | |
| | | { | |
| | | tmp_vec[col] = X.at(row,col); | |
| } | | } | |
|
| | | | |
| | | out[row] = op_median::direct_median(tmp_vec); | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T> | | template<typename T> | |
| arma_inline | | arma_inline | |
| std::complex<T> | | std::complex<T> | |
| op_median::robust_mean(const std::complex<T>& A, const std::complex<T>& B) | | op_median::robust_mean(const std::complex<T>& A, const std::complex<T>& B) | |
| { | | { | |
| return A + (B - A)/T(2); | | return A + (B - A)/T(2); | |
| | | | |
| skipping to change at line 302 | | skipping to change at line 300 | |
| const u32 X_n_rows = X.n_rows; | | const u32 X_n_rows = X.n_rows; | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| const u32 dim = in.aux_u32_a; | | const u32 dim = in.aux_u32_a; | |
| arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or
1"); | | arma_debug_check( (dim > 1), "median(): incorrect usage. dim must be 0 or
1"); | |
| | | | |
| if(dim == 0) // in each column | | if(dim == 0) // in each column | |
| { | | { | |
| arma_extra_debug_print("op_median::apply(), dim = 0"); | | arma_extra_debug_print("op_median::apply(), dim = 0"); | |
| | | | |
|
| out.set_size((X_n_rows > 0) ? 1 : 0, X_n_cols); | | arma_debug_check( (X_n_rows == 0), "median(): given object has zero row
s" ); | |
| | | | |
|
| if(X_n_rows > 0) | | out.set_size(1, X_n_cols); | |
| | | | |
| | | std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_rows); | |
| | | | |
| | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
|
| std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_rows); | | const eT* colmem = X.colptr(col); | |
| | | | |
|
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 row=0; row<X_n_rows; ++row) | |
| { | | { | |
|
| const eT* colmem = X.colptr(col); | | tmp_vec[row].val = std::abs(colmem[row]); | |
| | | tmp_vec[row].index = row; | |
| | | } | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | u32 index1; | |
| { | | u32 index2; | |
| tmp_vec[row].val = std::abs(colmem[row]); | | op_median::direct_cx_median_index(index1, index2, tmp_vec); | |
| tmp_vec[row].index = row; | | | |
| } | | | |
| | | | |
| u32 index1; | | | |
| u32 index2; | | | |
| op_median::direct_cx_median_index(index1, index2, tmp_vec); | | | |
| | | | |
|
| out[col] = op_median::robust_mean(colmem[index1], colmem[index2]); | | out[col] = op_median::robust_mean(colmem[index1], colmem[index2]); | |
| } | | | |
| } | | } | |
| } | | } | |
| else | | else | |
| if(dim == 1) // in each row | | if(dim == 1) // in each row | |
| { | | { | |
| arma_extra_debug_print("op_median::apply(), dim = 1"); | | arma_extra_debug_print("op_median::apply(), dim = 1"); | |
| | | | |
|
| out.set_size(X_n_rows, (X_n_cols > 0) ? 1 : 0); | | arma_debug_check( (X_n_cols == 0), "median(): given object has zero col
umns" ); | |
| | | | |
|
| if(X_n_cols > 0) | | out.set_size(X_n_rows, 1); | |
| { | | | |
| std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_cols); | | | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | std::vector< arma_cx_median_packet<T> > tmp_vec(X_n_cols); | |
| | | | |
| | | for(u32 row=0; row<X_n_rows; ++row) | |
| | | { | |
| | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
|
| for(u32 col=0; col<X_n_cols; ++col) | | tmp_vec[col].val = std::abs(X.at(row,col)); | |
| { | | tmp_vec[row].index = col; | |
| tmp_vec[col].val = std::abs(X.at(row,col)); | | | |
| tmp_vec[row].index = col; | | | |
| } | | | |
| | | | |
| if(X_n_cols > 0) | | | |
| { | | | |
| u32 index1; | | | |
| u32 index2; | | | |
| op_median::direct_cx_median_index(index1, index2, tmp_vec); | | | |
| | | | |
| out[row] = op_median::robust_mean( X.at(row,index1), X.at(row,ind | | | |
| ex2) ); | | | |
| } | | | |
| else | | | |
| { | | | |
| out[row] = eT(0); | | | |
| } | | | |
| } | | } | |
|
| | | | |
| | | u32 index1; | |
| | | u32 index2; | |
| | | op_median::direct_cx_median_index(index1, index2, tmp_vec); | |
| | | | |
| | | out[row] = op_median::robust_mean( X.at(row,index1), X.at(row,index2) | |
| | | ); | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 23 change blocks. |
| 63 lines changed or deleted | | 52 lines changed or added | |
|
| op_min_meat.hpp | | op_min_meat.hpp | |
| | | | |
| skipping to change at line 24 | | skipping to change at line 24 | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| arma_pure | | arma_pure | |
| inline | | inline | |
| eT | | eT | |
| op_min::direct_min(const eT* const X, const u32 n_elem) | | op_min::direct_min(const eT* const X, const u32 n_elem) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| eT min_val = (n_elem != 1) ? priv::most_pos<eT>() : X[0]; | | eT min_val = priv::most_pos<eT>(); | |
| | | | |
| u32 i,j; | | u32 i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| const eT X_i = X[i]; | | const eT X_i = X[i]; | |
| const eT X_j = X[j]; | | const eT X_j = X[j]; | |
| | | | |
| if(X_i < min_val) | | if(X_i < min_val) | |
| { | | { | |
| | | | |
| skipping to change at line 64 | | skipping to change at line 64 | |
| return min_val; | | return min_val; | |
| } | | } | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_min::direct_min(const eT* const X, const u32 n_elem, u32& index_of_min_v
al) | | op_min::direct_min(const eT* const X, const u32 n_elem, u32& index_of_min_v
al) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
|
| eT min_val = (n_elem != 1) ? priv::most_pos<eT>() : X[0]; | | eT min_val = priv::most_pos<eT>(); | |
| | | | |
| u32 best_index = 0; | | u32 best_index = 0; | |
| | | | |
| u32 i,j; | | u32 i,j; | |
| | | | |
| for(i=0, j=1; j<n_elem; i+=2, j+=2) | | for(i=0, j=1; j<n_elem; i+=2, j+=2) | |
| { | | { | |
| const eT X_i = X[i]; | | const eT X_i = X[i]; | |
| const eT X_j = X[j]; | | const eT X_j = X[j]; | |
| | | | |
| | | | |
| skipping to change at line 113 | | skipping to change at line 113 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_min::direct_min(const Mat<eT>& X, const u32 row) | | op_min::direct_min(const Mat<eT>& X, const u32 row) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
|
| eT min_val = (X_n_cols != 1) ? priv::most_pos<eT>() : X.at(row,0); | | eT min_val = priv::most_pos<eT>(); | |
| | | | |
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| const eT tmp_val = X.at(row,col); | | const eT tmp_val = X.at(row,col); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 137 | | skipping to change at line 137 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_min::direct_min(const subview<eT>& X) | | op_min::direct_min(const subview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
|
| eT min_val = (X_n_elem != 1) ? priv::most_pos<eT>() : X[0]; | | eT min_val = priv::most_pos<eT>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| eT tmp_val = X[i]; | | eT tmp_val = X[i]; | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 161 | | skipping to change at line 161 | |
| | | | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
| eT | | eT | |
| op_min::direct_min(const diagview<eT>& X) | | op_min::direct_min(const diagview<eT>& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| | | | |
|
| eT min_val = (X_n_elem != 1) ? priv::most_pos<eT>() : X[0]; | | eT min_val = priv::most_pos<eT>();; | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| eT tmp_val = X[i]; | | eT tmp_val = X[i]; | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| } | | } | |
| } | | } | |
| | | | |
| skipping to change at line 200 | | skipping to change at line 200 | |
| const u32 dim = in.aux_u32_a; | | const u32 dim = in.aux_u32_a; | |
| arma_debug_check( (dim > 1), "min(): incorrect usage. dim must be 0 or 1"
); | | arma_debug_check( (dim > 1), "min(): incorrect usage. dim must be 0 or 1"
); | |
| | | | |
| const u32 X_n_rows = X.n_rows; | | const u32 X_n_rows = X.n_rows; | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| if(dim == 0) // min in each column | | if(dim == 0) // min in each column | |
| { | | { | |
| arma_extra_debug_print("op_min::apply(), dim = 0"); | | arma_extra_debug_print("op_min::apply(), dim = 0"); | |
| | | | |
|
| out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); | | arma_debug_check( (X_n_rows == 0), "min(): given object has zero rows"
); | |
| | | | |
|
| if(X_n_rows > 0) | | out.set_size(1, X_n_cols); | |
| { | | | |
| eT* out_mem = out.memptr(); | | eT* out_mem = out.memptr(); | |
| | | | |
|
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| out_mem[col] = op_min::direct_min( X.colptr(col), X_n_rows ); | | out_mem[col] = op_min::direct_min( X.colptr(col), X_n_rows ); | |
| } | | | |
| } | | } | |
| } | | } | |
| else | | else | |
| if(dim == 1) // min in each row | | if(dim == 1) // min in each row | |
| { | | { | |
| arma_extra_debug_print("op_min::apply(), dim = 1"); | | arma_extra_debug_print("op_min::apply(), dim = 1"); | |
| | | | |
|
| out.set_size( X_n_rows, (X_n_cols > 0) ? 1 : 0 ); | | arma_debug_check( (X_n_cols == 0), "min(): given object has zero column
s" ); | |
| | | | |
|
| if(X_n_cols > 0) | | out.set_size(X_n_rows, 1); | |
| { | | | |
| eT* out_mem = out.memptr(); | | eT* out_mem = out.memptr(); | |
| | | | |
|
| for(u32 row=0; row<X_n_rows; ++row) | | for(u32 row=0; row<X_n_rows; ++row) | |
| { | | { | |
| out_mem[row] = op_min::direct_min( X, row ); | | out_mem[row] = op_min::direct_min( X, row ); | |
| } | | | |
| } | | } | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_min::direct_min(const std::complex<T>* const X, const u32 n_elem) | | op_min::direct_min(const std::complex<T>* const X, const u32 n_elem) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T min_val = (n_elem != 1) ? priv::most_pos<T>() : std::abs(X[0]); | | T min_val = priv::most_pos<T>(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 263 | | skipping to change at line 261 | |
| } | | } | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_min::direct_min(const std::complex<T>* const X, const u32 n_elem, u32& i
ndex_of_min_val) | | op_min::direct_min(const std::complex<T>* const X, const u32 n_elem, u32& i
ndex_of_min_val) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T min_val = (n_elem != 1) ? priv::most_pos<T>() : std::abs(X[0]); | | T min_val = priv::most_pos<T>(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 291 | | skipping to change at line 289 | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_min::direct_min(const Mat< std::complex<T> >& X, const u32 row) | | op_min::direct_min(const Mat< std::complex<T> >& X, const u32 row) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_cols = X.n_cols; | | const u32 X_n_cols = X.n_cols; | |
| | | | |
| u32 index = 0; | | u32 index = 0; | |
|
| T min_val = (X_n_cols != 1) ? priv::most_pos<T>() : std::abs(X.at(row,0
)); | | T min_val = priv::most_pos<T>(); | |
| | | | |
| for(u32 col=0; col<X_n_cols; ++col) | | for(u32 col=0; col<X_n_cols; ++col) | |
| { | | { | |
| const T tmp_val = std::abs(X.at(row,col)); | | const T tmp_val = std::abs(X.at(row,col)); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| index = col; | | index = col; | |
| } | | } | |
| | | | |
| skipping to change at line 316 | | skipping to change at line 314 | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_min::direct_min(const subview< std::complex<T> >& X) | | op_min::direct_min(const subview< std::complex<T> >& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| u32 index = 0; | | u32 index = 0; | |
|
| T min_val = (X_n_elem != 1) ? priv::most_pos<T>() : std::abs(X[0
]); | | T min_val = priv::most_pos<T>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
| skipping to change at line 341 | | skipping to change at line 339 | |
| | | | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
| std::complex<T> | | std::complex<T> | |
| op_min::direct_min(const diagview< std::complex<T> >& X) | | op_min::direct_min(const diagview< std::complex<T> >& X) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 X_n_elem = X.n_elem; | | const u32 X_n_elem = X.n_elem; | |
| u32 index = 0; | | u32 index = 0; | |
|
| T min_val = (X_n_elem != 1) ? priv::most_pos<T>() : std::abs(X[0
]); | | T min_val = priv::most_pos<T>(); | |
| | | | |
| for(u32 i=0; i<X_n_elem; ++i) | | for(u32 i=0; i<X_n_elem; ++i) | |
| { | | { | |
| const T tmp_val = std::abs(X[i]); | | const T tmp_val = std::abs(X[i]); | |
| | | | |
| if(tmp_val < min_val) | | if(tmp_val < min_val) | |
| { | | { | |
| min_val = tmp_val; | | min_val = tmp_val; | |
| index = i; | | index = i; | |
| } | | } | |
| | | | |
End of changes. 16 change blocks. |
| 26 lines changed or deleted | | 24 lines changed or added | |
|
| op_misc_meat.hpp | | op_misc_meat.hpp | |
|
| // Copyright (C) 2008-2010 NICTA (www.nicta.com.au) | | // Copyright (C) 2008-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2008-2010 Conrad Sanderson | | // Copyright (C) 2008-2011 Conrad Sanderson | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| //! \addtogroup op_misc | | //! \addtogroup op_misc | |
| //! @{ | | //! @{ | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_real::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::po | | void | |
| d_type, T1, op_real>& X) | | op_real::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::po | |
| | | d_type, T1, op_real>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const Proxy<T1> A(X.m); | | const Proxy<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols()); | | out.set_size(A.get_n_rows(), A.get_n_cols()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::real(A[i]); | | out_mem[i] = std::real(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_real::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T | | void | |
| 1::pod_type, T1, op_real>& X) | | op_real::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T | |
| | | 1::pod_type, T1, op_real>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const ProxyCube<T1> A(X.m); | | const ProxyCube<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | | out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::real(A[i]); | | out_mem[i] = std::real(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_imag::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::po | | void | |
| d_type, T1, op_imag>& X) | | op_imag::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::po | |
| | | d_type, T1, op_imag>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const Proxy<T1> A(X.m); | | const Proxy<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols()); | | out.set_size(A.get_n_rows(), A.get_n_cols()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::imag(A[i]); | | out_mem[i] = std::imag(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_imag::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T | | void | |
| 1::pod_type, T1, op_imag>& X) | | op_imag::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T | |
| | | 1::pod_type, T1, op_imag>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const ProxyCube<T1> A(X.m); | | const ProxyCube<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | | out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::imag(A[i]); | | out_mem[i] = std::imag(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_abs::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::pod | | void | |
| _type, T1, op_abs>& X) | | op_abs::apply( Mat<typename T1::pod_type>& out, const mtOp<typename T1::pod | |
| | | _type, T1, op_abs>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const Proxy<T1> A(X.m); | | const Proxy<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols()); | | out.set_size(A.get_n_rows(), A.get_n_cols()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::abs(A[i]); | | out_mem[i] = std::abs(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
|
| inline void | | inline | |
| op_abs::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1 | | void | |
| ::pod_type, T1, op_abs>& X) | | op_abs::apply( Cube<typename T1::pod_type>& out, const mtOpCube<typename T1 | |
| | | ::pod_type, T1, op_abs>& X ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| | | | |
| const ProxyCube<T1> A(X.m); | | const ProxyCube<T1> A(X.m); | |
| | | | |
| out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | | out.set_size(A.get_n_rows(), A.get_n_cols(), A.get_n_slices()); | |
| | | | |
| const u32 n_elem = out.n_elem; | | const u32 n_elem = out.n_elem; | |
| T* out_mem = out.memptr(); | | T* out_mem = out.memptr(); | |
| | | | |
| for(u32 i=0; i<n_elem; ++i) | | for(u32 i=0; i<n_elem; ++i) | |
| { | | { | |
| out_mem[i] = std::abs(A[i]); | | out_mem[i] = std::abs(A[i]); | |
| } | | } | |
| } | | } | |
| | | | |
|
| | | template<typename T1> | |
| | | inline | |
| | | void | |
| | | op_sympd::apply( Mat<typename T1::elem_type>& out, const Op<T1, op_sympd>& | |
| | | X ) | |
| | | { | |
| | | arma_extra_debug_sigprint(); | |
| | | | |
| | | out = X.m; | |
| | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 8 change blocks. |
| 20 lines changed or deleted | | 37 lines changed or added | |
|
| op_pinv_meat.hpp | | op_pinv_meat.hpp | |
| // Copyright (C) 2009-2011 NICTA (www.nicta.com.au) | | // Copyright (C) 2009-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2009-2011 Conrad Sanderson | | // Copyright (C) 2009-2011 Conrad Sanderson | |
| // Copyright (C) 2009-2010 Dimitrios Bouzas | | // Copyright (C) 2009-2010 Dimitrios Bouzas | |
|
| | | // Copyright (C) 2011 Stanislav Funiak | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| | | | |
| skipping to change at line 32 | | skipping to change at line 33 | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 n_rows = A.n_rows; | | const u32 n_rows = A.n_rows; | |
| const u32 n_cols = A.n_cols; | | const u32 n_cols = A.n_cols; | |
| | | | |
| // SVD decomposition | | // SVD decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<eT> s; | | Col<eT> s; | |
| Mat<eT> V; | | Mat<eT> V; | |
| | | | |
|
| const bool status = (n_cols > n_rows) ? svd(U,s,V,trans(A)) : svd(U,s,V,A | | const bool status = (n_cols > n_rows) ? auxlib::svd(U,s,V,trans(A)) : a | |
| ); | | uxlib::svd(U,s,V,A); | |
| | | const u32 s_n_elem = s.n_elem; | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
| out.reset(); | | out.reset(); | |
|
| | | arma_bad("pinv(): svd failed"); | |
| return; | | return; | |
| } | | } | |
| | | | |
| // set tolerance to default if it hasn't been specified as an argument | | // set tolerance to default if it hasn't been specified as an argument | |
|
| if(tol == eT(0)) | | if( (tol == eT(0)) && (s_n_elem > 0) ) | |
| { | | { | |
| tol = (std::max)(n_rows,n_cols) * eop_aux::direct_eps(max(s)); | | tol = (std::max)(n_rows,n_cols) * eop_aux::direct_eps(max(s)); | |
| } | | } | |
| | | | |
| // count non zero valued elements in s | | // count non zero valued elements in s | |
| | | | |
|
| const u32 s_n_rows = s.n_rows; | | u32 count = 0; | |
| u32 count = 0; | | | |
| | | | |
|
| for(u32 i=0; i < s_n_rows; ++i) | | for(u32 i=0; i < s_n_elem; ++i) | |
| { | | { | |
| if(s[i] > tol) | | if(s[i] > tol) | |
| { | | { | |
| ++count; | | ++count; | |
| } | | } | |
| } | | } | |
| | | | |
| if(count != 0) | | if(count != 0) | |
| { | | { | |
| // reduce the length of s in order to contain only the values above tol
erance | | // reduce the length of s in order to contain only the values above tol
erance | |
|
| if(count < s_n_rows) | | if(count < s_n_elem) | |
| { | | { | |
| //s = s.rows(0,count-1); | | //s = s.rows(0,count-1); | |
|
| s.shed_rows(count, s_n_rows-1); | | s.shed_rows(count, s_n_elem-1); | |
| } | | } | |
| | | | |
| // set the elements of s equal to their reciprocals | | // set the elements of s equal to their reciprocals | |
| s = eT(1) / s; | | s = eT(1) / s; | |
| | | | |
| if(A.n_cols <= A.n_rows) | | if(A.n_cols <= A.n_rows) | |
| { | | { | |
| out = ( V.n_cols > count ? V.cols(0,count-1) : V ) * diagmat(s) * tra
ns( U.n_cols > count ? U.cols(0,count-1) : U ); | | out = ( V.n_cols > count ? V.cols(0,count-1) : V ) * diagmat(s) * tra
ns( U.n_cols > count ? U.cols(0,count-1) : U ); | |
| } | | } | |
| else | | else | |
| | | | |
| skipping to change at line 103 | | skipping to change at line 105 | |
| const u32 n_rows = A.n_rows; | | const u32 n_rows = A.n_rows; | |
| const u32 n_cols = A.n_cols; | | const u32 n_cols = A.n_cols; | |
| | | | |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; | |
| | | | |
| // SVD decomposition | | // SVD decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col< T> s; | | Col< T> s; | |
| Mat<eT> V; | | Mat<eT> V; | |
| | | | |
|
| const bool status = (n_cols > n_rows) ? svd(U,s,V,trans(A)) : svd(U,s,V,A | | const bool status = (n_cols > n_rows) ? auxlib::svd(U,s,V,trans(A)) : au | |
| ); | | xlib::svd(U,s,V,A); | |
| | | const u32 s_n_elem = s.n_elem; | |
| | | | |
| if(status == false) | | if(status == false) | |
| { | | { | |
| out.reset(); | | out.reset(); | |
|
| | | arma_bad("pinv(): svd failed"); | |
| return; | | return; | |
| } | | } | |
| | | | |
| // set tolerance to default if it hasn't been specified as an argument | | // set tolerance to default if it hasn't been specified as an argument | |
|
| if(tol == T(0)) | | if( (tol == T(0)) && (s_n_elem > 0) ) | |
| { | | { | |
| tol = (std::max)(n_rows,n_cols) * eop_aux::direct_eps(max(s)); | | tol = (std::max)(n_rows,n_cols) * eop_aux::direct_eps(max(s)); | |
| } | | } | |
| | | | |
| // count non zero valued elements in s | | // count non zero valued elements in s | |
| | | | |
|
| const u32 s_n_rows = s.n_rows; | | u32 count = 0; | |
| u32 count = 0; | | | |
| | | | |
|
| for(u32 i = 0; i < s_n_rows; ++i) | | for(u32 i = 0; i < s_n_elem; ++i) | |
| { | | { | |
| if(s[i] > tol) | | if(s[i] > tol) | |
| { | | { | |
| ++count; | | ++count; | |
| } | | } | |
| } | | } | |
| | | | |
| if(count != 0) | | if(count != 0) | |
| { | | { | |
| // reduce the length of s in order to contain only the values above tol
erance | | // reduce the length of s in order to contain only the values above tol
erance | |
|
| if(count < s_n_rows) | | if(count < s_n_elem) | |
| { | | { | |
| // s = s.rows(0,count-1); | | // s = s.rows(0,count-1); | |
|
| s.shed_rows(count, s_n_rows-1); | | s.shed_rows(count, s_n_elem-1); | |
| } | | } | |
| | | | |
| // set the elements of s equal to their reciprocals | | // set the elements of s equal to their reciprocals | |
| s = T(1) / s; | | s = T(1) / s; | |
| | | | |
| if(n_rows >= n_cols) | | if(n_rows >= n_cols) | |
| { | | { | |
| out = ( V.n_cols > count ? V.cols(0,count-1) : V ) * diagmat(s) * tra
ns( U.n_cols > count ? U.cols(0,count-1) : U ); | | out = ( V.n_cols > count ? V.cols(0,count-1) : V ) * diagmat(s) * tra
ns( U.n_cols > count ? U.cols(0,count-1) : U ); | |
| } | | } | |
| else | | else | |
| | | | |
| skipping to change at line 168 | | skipping to change at line 171 | |
| inline | | inline | |
| void | | void | |
| op_pinv::apply(Mat<typename T1::pod_type>& out, const Op<T1,op_pinv>& in) | | op_pinv::apply(Mat<typename T1::pod_type>& out, const Op<T1,op_pinv>& in) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type eT; | | typedef typename T1::pod_type eT; | |
| | | | |
| const eT tol = in.aux; | | const eT tol = in.aux; | |
| | | | |
|
| arma_debug_check((tol < eT(0)), "pinv(): tol must be >= 0"); | | arma_debug_check((tol < eT(0)), "pinv(): tolerance must be >= 0"); | |
| | | | |
| const unwrap<T1> tmp(in.m); | | const unwrap<T1> tmp(in.m); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
| op_pinv::direct_pinv(out, A, tol); | | op_pinv::direct_pinv(out, A, tol); | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| op_pinv::apply(Mat< std::complex<typename T1::pod_type> >& out, const Op<T1
,op_pinv>& in) | | op_pinv::apply(Mat< std::complex<typename T1::pod_type> >& out, const Op<T1
,op_pinv>& in) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::pod_type T; | | typedef typename T1::pod_type T; | |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; | |
| | | | |
| const T tol = in.aux.real(); | | const T tol = in.aux.real(); | |
| | | | |
|
| arma_debug_check((tol < T(0)), "pinv(): tol must be >= 0"); | | arma_debug_check((tol < T(0)), "pinv(): tolerance must be >= 0"); | |
| | | | |
| const unwrap<T1> tmp(in.m); | | const unwrap<T1> tmp(in.m); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
| op_pinv::direct_pinv(out, A, tol); | | op_pinv::direct_pinv(out, A, tol); | |
| } | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 17 change blocks. |
| 18 lines changed or deleted | | 21 lines changed or added | |
|
| op_princomp_meat.hpp | | op_princomp_meat.hpp | |
|
| // Copyright (C) 2010 NICTA (www.nicta.com.au) | | // Copyright (C) 2010-2011 NICTA (www.nicta.com.au) | |
| // Copyright (C) 2010 Conrad Sanderson | | // Copyright (C) 2010-2011 Conrad Sanderson | |
| // Copyright (C) 2010 Dimitrios Bouzas | | // Copyright (C) 2010 Dimitrios Bouzas | |
|
| | | // Copyright (C) 2011 Stanislav Funiak | |
| // | | // | |
| // This file is part of the Armadillo C++ library. | | // This file is part of the Armadillo C++ library. | |
| // It is provided without any warranty of fitness | | // It is provided without any warranty of fitness | |
| // for any purpose. You can redistribute this file | | // for any purpose. You can redistribute this file | |
| // and/or modify it under the terms of the GNU | | // and/or modify it under the terms of the GNU | |
| // Lesser General Public License (LGPL) as published | | // Lesser General Public License (LGPL) as published | |
| // by the Free Software Foundation, either version 3 | | // by the Free Software Foundation, either version 3 | |
| // of the License or (at your option) any later version. | | // of the License or (at your option) any later version. | |
| // (see http://www.opensource.org/licenses for more info) | | // (see http://www.opensource.org/licenses for more info) | |
| | | | |
| | | | |
| skipping to change at line 26 | | skipping to change at line 27 | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 4 arguments version | | //! principal component analysis -- 4 arguments version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| //! tsquared_out -> Hotelling's T^2 statistic | | //! tsquared_out -> Hotelling's T^2 statistic | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat<eT>& coeff_out, | | Mat<eT>& coeff_out, | |
| Mat<eT>& score_out, | | Mat<eT>& score_out, | |
| Col<eT>& latent_out, | | Col<eT>& latent_out, | |
| Col<eT>& tsquared_out, | | Col<eT>& tsquared_out, | |
| const Mat<eT>& in | | const Mat<eT>& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| skipping to change at line 54 | | skipping to change at line 55 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<eT> s; | | Col<eT> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| tsquared_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| //U.reset(); // TODO: do we need this ? U will get automatically dele
ted anyway | | //U.reset(); // TODO: do we need this ? U will get automatically dele
ted anyway | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| | | | |
| skipping to change at line 99 | | skipping to change at line 93 | |
| else | | else | |
| { | | { | |
| // compute the Hotelling's T-squared | | // compute the Hotelling's T-squared | |
| const Mat<eT> S = score_out * diagmat(Col<eT>( eT(1) / s)); | | const Mat<eT> S = score_out * diagmat(Col<eT>( eT(1) / s)); | |
| tsquared_out = sum(S%S,1); | | tsquared_out = sum(S%S,1); | |
| } | | } | |
| | | | |
| // compute the eigenvalues of the principal vectors | | // compute the eigenvalues of the principal vectors | |
| latent_out = s%s; | | latent_out = s%s; | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | | |
| | | | |
|
| score_out.copy_size(in); | | score_out.copy_size(in); | |
| score_out.zeros(); | | score_out.zeros(); | |
| | | | |
|
| latent_out.set_size(n_cols); | | latent_out.set_size(n_cols); | |
| latent_out.zeros(); | | latent_out.zeros(); | |
| | | | |
|
| tsquared_out.set_size(1); | | tsquared_out.set_size(n_rows); | |
| tsquared_out.zeros(); | | tsquared_out.zeros(); | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| tsquared_out.reset(); | | | |
| } | | | |
| } | | } | |
| | | | |
|
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 3 arguments version | | //! principal component analysis -- 3 arguments version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat<eT>& coeff_out, | | Mat<eT>& coeff_out, | |
| Mat<eT>& score_out, | | Mat<eT>& score_out, | |
| Col<eT>& latent_out, | | Col<eT>& latent_out, | |
| const Mat<eT>& in | | const Mat<eT>& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 160 | | skipping to change at line 145 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<eT> s; | | Col<eT> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| // U.reset(); | | // U.reset(); | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| | | | |
| skipping to change at line 190 | | skipping to change at line 169 | |
| | | | |
| Col<eT> s_tmp = zeros< Col<eT> >(n_cols); | | Col<eT> s_tmp = zeros< Col<eT> >(n_cols); | |
| s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | | s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | |
| s = s_tmp; | | s = s_tmp; | |
| } | | } | |
| | | | |
| // compute the eigenvalues of the principal vectors | | // compute the eigenvalues of the principal vectors | |
| latent_out = s%s; | | latent_out = s%s; | |
| | | | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | | |
| | | | |
|
| score_out.copy_size(in); | | score_out.copy_size(in); | |
| score_out.zeros(); | | score_out.zeros(); | |
| | | | |
|
| latent_out.set_size(n_cols); | | latent_out.set_size(n_cols); | |
| latent_out.zeros(); | | latent_out.zeros(); | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| } | | | |
| } | | } | |
| | | | |
|
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 2 arguments version | | //! principal component analysis -- 2 arguments version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat<eT>& coeff_out, | | Mat<eT>& coeff_out, | |
| Mat<eT>& score_out, | | Mat<eT>& score_out, | |
| const Mat<eT>& in | | const Mat<eT>& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| const u32 n_rows = in.n_rows; | | const u32 n_rows = in.n_rows; | |
| | | | |
| skipping to change at line 245 | | skipping to change at line 216 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<eT> s; | | Col<eT> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| // U.reset(); | | // U.reset(); | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| if(n_rows <= n_cols) // number of samples is less than their dimensiona
lity | | if(n_rows <= n_cols) // number of samples is less than their dimensiona
lity | |
| { | | { | |
| score_out.cols(n_rows-1,n_cols-1).zeros(); | | score_out.cols(n_rows-1,n_cols-1).zeros(); | |
| | | | |
| Col<eT> s_tmp = zeros< Col<eT> >(n_cols); | | Col<eT> s_tmp = zeros< Col<eT> >(n_cols); | |
| s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | | s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | |
| s = s_tmp; | | s = s_tmp; | |
| } | | } | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | score_out.copy_size(in); | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | score_out.zeros(); | |
| score_out.copy_size(in); | | | |
| score_out.zeros(); | | | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| } | | | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 1 argument version | | //! principal component analysis -- 1 argument version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| template<typename eT> | | template<typename eT> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat<eT>& coeff_out, | | Mat<eT>& coeff_out, | |
| const Mat<eT>& in | | const Mat<eT>& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| if(in.n_elem != 0) | | if(in.n_elem != 0) | |
| { | | { | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<eT> s; | | Col<eT> s; | |
| | | | |
| const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1); | | const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1); | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out, tmp); | | const bool svd_ok = svd(U,s,coeff_out, tmp); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| } | | } | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| coeff_out.reset(); | | coeff_out.eye(in.n_cols, in.n_cols); | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 4 arguments complex version | | //! principal component analysis -- 4 arguments complex version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| //! tsquared_out -> Hotelling's T^2 statistic | | //! tsquared_out -> Hotelling's T^2 statistic | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat< std::complex<T> >& coeff_out, | | Mat< std::complex<T> >& coeff_out, | |
| Mat< std::complex<T> >& score_out, | | Mat< std::complex<T> >& score_out, | |
| Col<T>& latent_out, | | Col<T>& latent_out, | |
| Col< std::complex<T> >& tsquared_out, | | Col< std::complex<T> >& tsquared_out, | |
| const Mat< std::complex<T> >& in | | const Mat< std::complex<T> >& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| skipping to change at line 363 | | skipping to change at line 323 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col<T> s; | | Col<T> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| tsquared_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| //U.reset(); | | //U.reset(); | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| | | | |
| skipping to change at line 405 | | skipping to change at line 358 | |
| { | | { | |
| // compute the Hotelling's T-squared | | // compute the Hotelling's T-squared | |
| const Mat<eT> S = score_out * diagmat(Col<T>(T(1) / s)); | | const Mat<eT> S = score_out * diagmat(Col<T>(T(1) / s)); | |
| tsquared_out = sum(S%S,1); | | tsquared_out = sum(S%S,1); | |
| } | | } | |
| | | | |
| // compute the eigenvalues of the principal vectors | | // compute the eigenvalues of the principal vectors | |
| latent_out = s%s; | | latent_out = s%s; | |
| | | | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | | |
| | | | |
|
| score_out.copy_size(in); | | score_out.copy_size(in); | |
| score_out.zeros(); | | score_out.zeros(); | |
| | | | |
|
| latent_out.set_size(n_cols); | | latent_out.set_size(n_cols); | |
| latent_out.zeros(); | | latent_out.zeros(); | |
| | | | |
|
| tsquared_out.set_size(1); | | tsquared_out.set_size(n_rows); | |
| tsquared_out.zeros(); | | tsquared_out.zeros(); | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| tsquared_out.reset(); | | | |
| } | | | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 3 arguments complex version | | //! principal component analysis -- 3 arguments complex version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| //! latent_out -> eigenvalues of principal vectors | | //! latent_out -> eigenvalues of principal vectors | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat< std::complex<T> >& coeff_out, | | Mat< std::complex<T> >& coeff_out, | |
| Mat< std::complex<T> >& score_out, | | Mat< std::complex<T> >& score_out, | |
| Col<T>& latent_out, | | Col<T>& latent_out, | |
| const Mat< std::complex<T> >& in | | const Mat< std::complex<T> >& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| | | | |
| skipping to change at line 467 | | skipping to change at line 412 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col< T> s; | | Col< T> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| // U.reset(); | | // U.reset(); | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| | | | |
| skipping to change at line 497 | | skipping to change at line 436 | |
| | | | |
| Col<T> s_tmp = zeros< Col<T> >(n_cols); | | Col<T> s_tmp = zeros< Col<T> >(n_cols); | |
| s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | | s_tmp.rows(0,n_rows-2) = s.rows(0,n_rows-2); | |
| s = s_tmp; | | s = s_tmp; | |
| } | | } | |
| | | | |
| // compute the eigenvalues of the principal vectors | | // compute the eigenvalues of the principal vectors | |
| latent_out = s%s; | | latent_out = s%s; | |
| | | | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | score_out.copy_size(in); | |
| score_out.copy_size(in); | | score_out.zeros(); | |
| score_out.zeros(); | | | |
| latent_out.set_size(n_cols); | | latent_out.set_size(n_cols); | |
| latent_out.zeros(); | | latent_out.zeros(); | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| latent_out.reset(); | | | |
| } | | | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 2 arguments complex version | | //! principal component analysis -- 2 arguments complex version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| //! score_out -> projected samples | | //! score_out -> projected samples | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat< std::complex<T> >& coeff_out, | | Mat< std::complex<T> >& coeff_out, | |
| Mat< std::complex<T> >& score_out, | | Mat< std::complex<T> >& score_out, | |
| const Mat< std::complex<T> >& in | | const Mat< std::complex<T> >& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef std::complex<T> eT; | | typedef std::complex<T> eT; | |
| | | | |
| skipping to change at line 551 | | skipping to change at line 485 | |
| score_out = in - repmat(mean(in), n_rows, 1); | | score_out = in - repmat(mean(in), n_rows, 1); | |
| | | | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col< T> s; | | Col< T> s; | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out,score_out); | | const bool svd_ok = svd(U,s,coeff_out,score_out); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| | | | |
| return; | | | |
| } | | } | |
| | | | |
| // U.reset(); | | // U.reset(); | |
| | | | |
| // normalize the eigenvalues | | // normalize the eigenvalues | |
| s /= std::sqrt( double(n_rows - 1) ); | | s /= std::sqrt( double(n_rows - 1) ); | |
| | | | |
| // project the samples to the principals | | // project the samples to the principals | |
| score_out *= coeff_out; | | score_out *= coeff_out; | |
| | | | |
| if(n_rows <= n_cols) // number of samples is less than their dimensiona
lity | | if(n_rows <= n_cols) // number of samples is less than their dimensiona
lity | |
| { | | { | |
| score_out.cols(n_rows-1,n_cols-1).zeros(); | | score_out.cols(n_rows-1,n_cols-1).zeros(); | |
| } | | } | |
| | | | |
| } | | } | |
|
| else // single sample - row | | else // 0 or 1 samples | |
| { | | { | |
|
| if(n_rows == 1) | | coeff_out.eye(n_cols, n_cols); | |
| { | | | |
| coeff_out = eye< Mat<eT> >(n_cols, n_cols); | | | |
| | | | |
|
| score_out.copy_size(in); | | score_out.copy_size(in); | |
| score_out.zeros(); | | score_out.zeros(); | |
| } | | | |
| else | | | |
| { | | | |
| coeff_out.reset(); | | | |
| score_out.reset(); | | | |
| } | | | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| //! \brief | | //! \brief | |
| //! principal component analysis -- 1 argument complex version | | //! principal component analysis -- 1 argument complex version | |
| //! computation is done via singular value decomposition | | //! computation is done via singular value decomposition | |
| //! coeff_out -> principal component coefficients | | //! coeff_out -> principal component coefficients | |
| template<typename T> | | template<typename T> | |
| inline | | inline | |
|
| void | | bool | |
| op_princomp::direct_princomp | | op_princomp::direct_princomp | |
| ( | | ( | |
| Mat< std::complex<T> >& coeff_out, | | Mat< std::complex<T> >& coeff_out, | |
| const Mat< std::complex<T> >& in | | const Mat< std::complex<T> >& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename std::complex<T> eT; | | typedef typename std::complex<T> eT; | |
| | | | |
| | | | |
| skipping to change at line 619 | | skipping to change at line 542 | |
| // singular value decomposition | | // singular value decomposition | |
| Mat<eT> U; | | Mat<eT> U; | |
| Col< T> s; | | Col< T> s; | |
| | | | |
| const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1); | | const Mat<eT> tmp = in - repmat(mean(in), in.n_rows, 1); | |
| | | | |
| const bool svd_ok = svd(U,s,coeff_out, tmp); | | const bool svd_ok = svd(U,s,coeff_out, tmp); | |
| | | | |
| if(svd_ok == false) | | if(svd_ok == false) | |
| { | | { | |
|
| arma_print("princomp(): singular value decomposition failed"); | | return false; | |
| | | | |
| coeff_out.reset(); | | | |
| } | | } | |
| } | | } | |
| else | | else | |
| { | | { | |
|
| coeff_out.reset(); | | coeff_out.eye(in.n_cols, in.n_cols); | |
| } | | } | |
|
| | | | |
| | | return true; | |
| } | | } | |
| | | | |
| template<typename T1> | | template<typename T1> | |
| inline | | inline | |
| void | | void | |
| op_princomp::apply | | op_princomp::apply | |
| ( | | ( | |
| Mat<typename T1::elem_type>& out, | | Mat<typename T1::elem_type>& out, | |
| const Op<T1,op_princomp>& in | | const Op<T1,op_princomp>& in | |
| ) | | ) | |
| { | | { | |
| arma_extra_debug_sigprint(); | | arma_extra_debug_sigprint(); | |
| | | | |
| typedef typename T1::elem_type eT; | | typedef typename T1::elem_type eT; | |
| | | | |
| const unwrap_check<T1> tmp(in.m, out); | | const unwrap_check<T1> tmp(in.m, out); | |
| const Mat<eT>& A = tmp.M; | | const Mat<eT>& A = tmp.M; | |
| | | | |
|
| op_princomp::direct_princomp(out, A); | | const bool status = op_princomp::direct_princomp(out, A); | |
| | | | |
| | | if(status == false) | |
| | | { | |
| | | out.reset(); | |
| | | | |
| | | arma_bad("princomp(): failed to converge"); | |
| | | } | |
| } | | } | |
| | | | |
| //! @} | | //! @} | |
| | | | |
End of changes. 50 change blocks. |
| 151 lines changed or deleted | | 81 lines changed or added | |
|