Cube_meat.hpp   Cube_meat.hpp 
skipping to change at line 164 skipping to change at line 164
if(new_n_elem <= Cube_prealloc::mem_n_elem) if(new_n_elem <= Cube_prealloc::mem_n_elem)
{ {
access::rw(mem) = mem_local; access::rw(mem) = mem_local;
} }
else else
{ {
arma_extra_debug_print("Cube::init(): allocating memory"); arma_extra_debug_print("Cube::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), "Cube::init(): out of memory" ); arma_check_bad_alloc( (mem == 0), "Cube::init(): out of memory" );
} }
if(new_n_elem > 0) if(new_n_elem > 0)
{ {
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(n_elem_slice) = in_n_rows*in_n_cols; access::rw(n_elem_slice) = in_n_rows*in_n_cols;
access::rw(n_slices) = in_n_slices; access::rw(n_slices) = in_n_slices;
access::rw(n_elem) = new_n_elem; access::rw(n_elem) = new_n_elem;
skipping to change at line 322 skipping to change at line 322
if(mem_state <= 2) if(mem_state <= 2)
{ {
if(n_slices <= Cube_prealloc::mat_ptrs_size) if(n_slices <= Cube_prealloc::mat_ptrs_size)
{ {
access::rw(mat_ptrs) = const_cast< const Mat<eT>** >(mat_ptrs_local); access::rw(mat_ptrs) = const_cast< const Mat<eT>** >(mat_ptrs_local);
} }
else else
{ {
access::rw(mat_ptrs) = new(std::nothrow) const Mat<eT>*[n_slices]; access::rw(mat_ptrs) = new(std::nothrow) const Mat<eT>*[n_slices];
arma_check( (mat_ptrs == 0), "Cube::create_mat(): out of memory" ); arma_check_bad_alloc( (mat_ptrs == 0), "Cube::create_mat(): out of me mory" );
} }
} }
for(u32 slice = 0; slice < n_slices; ++slice) for(u32 slice = 0; slice < n_slices; ++slice)
{ {
mat_ptrs[slice] = new Mat<eT>('j', slice_memptr(slice), n_rows, n_cols) ; mat_ptrs[slice] = new Mat<eT>('j', slice_memptr(slice), n_rows, n_cols) ;
} }
} }
//! Set the cube to be equal to the specified scalar. //! Set the cube to be equal to the specified scalar.
skipping to change at line 2405 skipping to change at line 2405
default: default:
arma_warn(print_status, "Cube::load(): unsupported file type"); arma_warn(print_status, "Cube::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("Cube::load(): ", err_msg, name); arma_warn(true, "Cube::load(): ", err_msg, name);
} }
else else
{ {
arma_print("Cube::load(): couldn't read ", name); arma_warn(true, "Cube::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 2467 skipping to change at line 2467
default: default:
arma_warn(print_status, "Cube::load(): unsupported file type"); arma_warn(print_status, "Cube::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("Cube::load(): ", err_msg, "the given stream"); arma_warn(true, "Cube::load(): ", err_msg, "the given stream");
} }
else else
{ {
arma_print("Cube::load(): couldn't load from the given stream"); arma_warn(true, "Cube::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. 6 change blocks. 
6 lines changed or deleted 6 lines changed or added


 Mat_bones.hpp   Mat_bones.hpp 
skipping to change at line 244 skipping to change at line 244
arma_inline arma_warn_unused eT operator() (const u32 in_row, const u32 in_col) const; arma_inline arma_warn_unused eT operator() (const u32 in_row, const u32 in_col) const;
arma_inline const Mat& operator++(); arma_inline const Mat& operator++();
arma_inline void operator++(int); arma_inline void operator++(int);
arma_inline const Mat& operator--(); arma_inline const Mat& operator--();
arma_inline void operator--(int); arma_inline void operator--(int);
arma_inline arma_warn_unused bool is_empty() const; arma_inline arma_warn_unused bool is_empty() const;
arma_inline arma_warn_unused bool is_vec() const; arma_inline arma_warn_unused bool is_vec() const;
arma_inline arma_warn_unused bool is_rowvec() const;
arma_inline arma_warn_unused bool is_colvec() const;
arma_inline arma_warn_unused bool is_square() const; arma_inline arma_warn_unused bool is_square() const;
inline arma_warn_unused bool is_finite() const; inline arma_warn_unused bool is_finite() const;
arma_inline arma_warn_unused bool in_range(const u32 i) const; arma_inline arma_warn_unused bool in_range(const u32 i) const;
arma_inline arma_warn_unused bool in_range(const span& x) const; arma_inline arma_warn_unused bool in_range(const span& x) const;
arma_inline arma_warn_unused bool in_range(const u32 in_row, const u3 2 in_col ) const; arma_inline arma_warn_unused bool in_range(const u32 in_row, const u3 2 in_col ) const;
arma_inline arma_warn_unused bool in_range(const span& row_span, const u3 2 in_col ) const; arma_inline arma_warn_unused bool in_range(const span& row_span, const u3 2 in_col ) const;
arma_inline arma_warn_unused bool in_range(const u32 in_row, const sp an& col_span) const; arma_inline arma_warn_unused bool in_range(const u32 in_row, const sp an& col_span) const;
arma_inline arma_warn_unused bool in_range(const span& row_span, const sp an& col_span) const; arma_inline arma_warn_unused bool in_range(const span& row_span, const sp an& col_span) const;
 End of changes. 1 change blocks. 
0 lines changed or deleted 2 lines changed or added


 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


 arma_version.hpp   arma_version.hpp 
skipping to change at line 18 skipping to change at line 18
// 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 arma_version //! \addtogroup arma_version
//! @{ //! @{
#define ARMA_VERSION_MAJOR 1 #define ARMA_VERSION_MAJOR 1
#define ARMA_VERSION_MINOR 99 #define ARMA_VERSION_MINOR 99
#define ARMA_VERSION_PATCH 4 #define ARMA_VERSION_PATCH 5
#define ARMA_VERSION_NAME "v2.0 beta 4" #define ARMA_VERSION_NAME "v2.0 beta 5"
struct arma_version struct arma_version
{ {
static const unsigned int major = ARMA_VERSION_MAJOR; static const unsigned int major = ARMA_VERSION_MAJOR;
static const unsigned int minor = ARMA_VERSION_MINOR; static const unsigned int minor = ARMA_VERSION_MINOR;
static const unsigned int patch = ARMA_VERSION_PATCH; static const unsigned int patch = ARMA_VERSION_PATCH;
static static
inline inline
std::string std::string
 End of changes. 1 change blocks. 
2 lines changed or deleted 2 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_bones.hpp   auxlib_bones.hpp 
skipping to change at line 62 skipping to change at line 62
template<typename eT, typename T1> template<typename eT, typename T1>
inline static bool inv_tr(Mat<eT>& out, const Base<eT,T1>& X, const u32 l ayout); inline static bool inv_tr(Mat<eT>& out, const Base<eT,T1>& X, const u32 l ayout);
// //
// inv_sym // inv_sym
template<typename eT, typename T1> template<typename eT, typename T1>
inline static bool inv_sym(Mat<eT>& out, const Base<eT,T1>& X, const u32 layout); inline static bool inv_sym(Mat<eT>& out, const Base<eT,T1>& X, const u32 layout);
// //
// inv_sympd
template<typename eT, typename T1>
inline static bool inv_sympd(Mat<eT>& out, const Base<eT,T1>& X, const u3
2 layout);
//
// det // det
template<typename eT, typename T1> template<typename eT, typename T1>
inline static eT det(const Base<eT,T1>& X, const bool slow = false); inline static eT det(const Base<eT,T1>& X, const bool slow = false);
template<typename eT> template<typename eT>
inline static eT det_tinymat(const Mat<eT>& X, const u32 N); inline static eT det_tinymat(const Mat<eT>& X, const u32 N);
template<typename eT> template<typename eT>
inline static eT det_lapack(const Mat<eT>& X, const bool make_copy); inline static eT det_lapack(const Mat<eT>& X, const bool make_copy);
// //
// log_det // log_det
template<typename eT, typename T1> template<typename eT, typename T1>
inline static void log_det(eT& out_val, typename get_pod_type<eT>::result & out_sign, const Base<eT,T1>& X); inline static bool log_det(eT& out_val, typename get_pod_type<eT>::result & out_sign, const Base<eT,T1>& X);
// //
// lu // lu
template<typename eT, typename T1> template<typename eT, typename T1>
inline static void lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, c onst Base<eT,T1>& X); inline static bool lu(Mat<eT>& L, Mat<eT>& U, podarray<blas_int>& ipiv, c onst Base<eT,T1>& X);
template<typename eT, typename T1> template<typename eT, typename T1>
inline static void lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T 1>& X); inline static bool lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Base<eT,T 1>& X);
template<typename eT, typename T1> template<typename eT, typename T1>
inline static void lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X); inline static bool lu(Mat<eT>& L, Mat<eT>& U, const Base<eT,T1>& X);
// //
// eig // eig
template<typename eT, typename T1> template<typename eT, typename T1>
inline static bool eig_sym(Col<eT>& eigval, const Base<eT,T1>& X); inline static bool eig_sym(Col<eT>& eigval, const Base<eT,T1>& X);
template<typename T, typename T1> template<typename T, typename T1>
inline static bool eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1> & X); inline static bool eig_sym(Col<T>& eigval, const Base<std::complex<T>,T1> & X);
 End of changes. 5 change blocks. 
4 lines changed or deleted 11 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


 config.hpp   config.hpp 
skipping to change at line 58 skipping to change at line 58
#define ARMA_USE_ATLAS #define ARMA_USE_ATLAS
#define ARMA_ATLAS_INCLUDE_DIR /usr/include/ #define ARMA_ATLAS_INCLUDE_DIR /usr/include/
//// If you're using ATLAS and the compiler can't find cblas.h and/or clapa ck.h //// If you're using ATLAS and the compiler can't find cblas.h and/or clapa ck.h
//// uncomment the above define and specify the appropriate include directo ry. //// uncomment the above define and specify the appropriate include directo ry.
//// Make sure the directory has a trailing / //// Make sure the directory has a trailing /
#define ARMA_USE_BOOST #define ARMA_USE_BOOST
#define ARMA_USE_BOOST_DATE #define ARMA_USE_BOOST_DATE
#define ARMA_PRINT_LOGIC_ERRORS
#define ARMA_PRINT_RUNTIME_ERRORS
#define ARMA_HAVE_STD_ISFINITE #define ARMA_HAVE_STD_ISFINITE
#define ARMA_HAVE_STD_ISINF #define ARMA_HAVE_STD_ISINF
#define ARMA_HAVE_STD_ISNAN #define ARMA_HAVE_STD_ISNAN
#define ARMA_HAVE_STD_SNPRINTF #define ARMA_HAVE_STD_SNPRINTF
#define ARMA_HAVE_LOG1P #define ARMA_HAVE_LOG1P
#define ARMA_HAVE_GETTIMEOFDAY #define ARMA_HAVE_GETTIMEOFDAY
// #define ARMA_EXTRA_DEBUG // #define ARMA_EXTRA_DEBUG
// #define ARMA_NO_DEBUG // #define ARMA_NO_DEBUG
#if defined(ARMA_DONT_USE_ATLAS)
#undef ARMA_USE_ATLAS
#undef ARMA_ATLAS_INCLUDE_DIR
#endif
 End of changes. 2 change blocks. 
0 lines changed or deleted 3 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


 field_meat.hpp   field_meat.hpp 
skipping to change at line 754 skipping to change at line 754
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
std::string err_msg; std::string err_msg;
const bool save_okay = field_aux::save(*this, name, type, err_msg); const bool save_okay = field_aux::save(*this, name, type, err_msg);
if( (print_status == true) && (save_okay == false) ) if( (print_status == true) && (save_okay == false) )
{ {
if(err_msg.length() > 0) if(err_msg.length() > 0)
{ {
arma_print("field::save(): ", err_msg, name); arma_warn(true, "field::save(): ", err_msg, name);
} }
else else
{ {
arma_print("field::save(): couldn't write to ", name); arma_warn(true, "field::save(): couldn't write to ", name);
} }
} }
return save_okay; return save_okay;
} }
template<typename oT> template<typename oT>
inline inline
bool bool
field<oT>::save(std::ostream& os, const file_type type, const bool print_st atus) const field<oT>::save(std::ostream& os, const file_type type, const bool print_st atus) const
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
std::string err_msg; std::string err_msg;
const bool save_okay = field_aux::save(*this, os, type, err_msg); const bool save_okay = field_aux::save(*this, os, type, err_msg);
if( (print_status == true) && (save_okay == false) ) if( (print_status == true) && (save_okay == false) )
{ {
if(err_msg.length() > 0) if(err_msg.length() > 0)
{ {
arma_print("field::save(): ", err_msg, "[ostream]"); arma_warn(true, "field::save(): ", err_msg, "[ostream]");
} }
else else
{ {
arma_print("field::save(): couldn't write to [ostream]"); arma_warn(true, "field::save(): couldn't write to [ostream]");
} }
} }
return save_okay; return save_okay;
} }
template<typename oT> template<typename oT>
inline inline
bool bool
field<oT>::load(const std::string name, const file_type type, const bool pr int_status) field<oT>::load(const std::string name, const file_type type, const bool pr int_status)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
std::string err_msg; std::string err_msg;
const bool load_okay = field_aux::load(*this, name, type, err_msg); const bool load_okay = field_aux::load(*this, name, type, err_msg);
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("field::load(): ", err_msg, name); arma_warn(true, "field::load(): ", err_msg, name);
} }
else else
{ {
arma_print("field::load(): couldn't read from ", name); arma_warn(true, "field::load(): couldn't read from ", name);
} }
} }
if(load_okay == false) if(load_okay == false)
{ {
(*this).reset(); (*this).reset();
} }
return load_okay; return load_okay;
} }
skipping to change at line 834 skipping to change at line 834
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
std::string err_msg; std::string err_msg;
const bool load_okay = field_aux::load(*this, is, type, err_msg); const bool load_okay = field_aux::load(*this, is, type, err_msg);
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("field::load(): ", err_msg, "[istream]"); arma_warn(true, "field::load(): ", err_msg, "[istream]");
} }
else else
{ {
arma_print("field::load(): couldn't read from [istream]"); arma_warn(true, "field::load(): couldn't read from [istream]");
} }
} }
if(load_okay == false) if(load_okay == false)
{ {
(*this).reset(); (*this).reset();
} }
return load_okay; return load_okay;
} }
skipping to change at line 946 skipping to change at line 946
delete [] mem; delete [] mem;
} }
if(n_elem_new <= sizeof(mem_local)/sizeof(oT*) ) if(n_elem_new <= sizeof(mem_local)/sizeof(oT*) )
{ {
mem = mem_local; mem = mem_local;
} }
else else
{ {
mem = new(std::nothrow) oT* [n_elem_new]; mem = new(std::nothrow) oT* [n_elem_new];
arma_check( (mem == 0), "field::init(): out of memory" ); arma_check_bad_alloc( (mem == 0), "field::init(): out of memory" );
} }
access::rw(n_elem) = n_elem_new; access::rw(n_elem) = n_elem_new;
if(n_elem_new == 0) if(n_elem_new == 0)
{ {
access::rw(n_rows) = 0; access::rw(n_rows) = 0;
access::rw(n_cols) = 0; access::rw(n_cols) = 0;
} }
else else
 End of changes. 9 change blocks. 
9 lines changed or deleted 9 lines changed or added


 fn_chol.hpp   fn_chol.hpp 
// Copyright (C) 2009-2010 NICTA (www.nicta.com.au) // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
// Copyright (C) 2009-2010 Conrad Sanderson // Copyright (C) 2009-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 fn_chol //! \addtogroup fn_chol
//! @{ //! @{
template<typename T1> template<typename T1>
inline inline
const Op<T1, op_chol> const Op<T1, op_chol>
chol(const Base<typename T1::elem_type,T1>& X, const typename arma_blas_typ chol
e_only<typename T1::elem_type>::result* junk = 0) (
const Base<typename T1::elem_type,T1>& X,
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);
return Op<T1, op_chol>(X.get_ref()); return Op<T1, op_chol>(X.get_ref());
} }
template<typename T1> template<typename T1>
inline inline
bool bool
chol(Mat<typename T1::elem_type>& out, const Base<typename T1::elem_type,T1 chol
>& X, const typename arma_blas_type_only<typename T1::elem_type>::result* j (
unk = 0) Mat<typename T1::elem_type>& out,
const Base<typename T1::elem_type,T1>& X,
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);
out = chol(X); try
{
out = chol(X);
}
catch(std::runtime_error&)
{
return false;
}
return (out.n_elem == 0) ? false : true; return true;
} }
//! @} //! @}
 End of changes. 7 change blocks. 
11 lines changed or deleted 24 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_inv.hpp   fn_inv.hpp 
skipping to change at line 63 skipping to change at line 63
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(slow); arma_ignore(slow);
arma_ignore(junk); arma_ignore(junk);
return Op<T1, op_inv_tr>(X.m, X.aux_u32_a, 0); return Op<T1, op_inv_tr>(X.m, X.aux_u32_a, 0);
} }
//! delayed matrix inverse (symmetric matrices) //! delayed matrix inverse (symmetric positive definite matrices)
template<typename T1> template<typename T1>
arma_inline arma_inline
const Op<T1, op_inv_sym> const Op<T1, op_inv_sympd>
inv inv
( (
const Op<T1, op_symmat>& X, const Op<T1, op_sympd>& X,
const bool slow = false, const bool slow = false,
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(slow); arma_ignore(slow);
arma_ignore(junk); arma_ignore(junk);
return Op<T1, op_inv_sym>(X.m, X.aux_u32_a, 0); return Op<T1, op_inv_sympd>(X.m, 0, 0);
}
template<typename T1>
inline
bool
inv
(
Mat<typename T1::elem_type>& out,
const Base<typename T1::elem_type,T1>& X,
const bool slow = false,
const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0
)
{
arma_extra_debug_sigprint();
arma_ignore(junk);
try
{
out = inv(X,slow);
}
catch(std::runtime_error&)
{
return false;
}
return true;
} }
//! @} //! @}
 End of changes. 4 change blocks. 
4 lines changed or deleted 31 lines changed or added


 fn_log_det.hpp   fn_log_det.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
// //
// 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 fn_log_det //! \addtogroup fn_log_det
//! @{ //! @{
//! log determinant of mat //! log determinant of mat
template<typename T1> template<typename T1>
inline inline
void bool
log_det log_det
( (
typename T1::elem_type& out_val, typename T1::elem_type& out_val,
typename T1::pod_type& out_sign, typename T1::pod_type& out_sign,
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);
auxlib::log_det(out_val, out_sign, X); return auxlib::log_det(out_val, out_sign, X);
} }
template<typename T1> template<typename T1>
inline inline
void void
log_det log_det
( (
typename T1::elem_type& out_val, typename T1::elem_type& out_val,
typename T1::pod_type& out_sign, typename T1::pod_type& out_sign,
const Op<T1,op_diagmat>& X, const Op<T1,op_diagmat>& 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);
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
typedef typename T1::pod_type T; typedef typename T1::pod_type T;
const diagmat_proxy<T1> A(X.m); const diagmat_proxy<T1> A(X.m);
const u32 N = A.n_elem; const u32 N = A.n_elem;
if(N == 0) if(N == 0)
{ {
 End of changes. 5 change blocks. 
4 lines changed or deleted 6 lines changed or added


 fn_lu.hpp   fn_lu.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 fn_lu //! \addtogroup fn_lu
//! @{ //! @{
//! immediate lower upper decomposition, permutation info is embedded into L (similar to Matlab/Octave) //! immediate lower upper decomposition, permutation info is embedded into L (similar to Matlab/Octave)
template<typename T1> template<typename T1>
inline inline
void bool
lu lu
( (
Mat<typename T1::elem_type>& L, Mat<typename T1::elem_type>& L,
Mat<typename T1::elem_type>& U, Mat<typename T1::elem_type>& U,
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( (&L == &U), "lu(): L and U are the same object"); arma_debug_check( (&L == &U), "lu(): L and U are the same object");
auxlib::lu(L, U, X); const bool status = auxlib::lu(L, U, X);
if(status == false)
{
L.reset();
U.reset();
arma_bad("lu(): failed to converge", false);
}
return status;
} }
//! immediate lower upper decomposition, also providing the permutation mat rix //! immediate lower upper decomposition, also providing the permutation mat rix
template<typename T1> template<typename T1>
inline inline
void bool
lu lu
( (
Mat<typename T1::elem_type>& L, Mat<typename T1::elem_type>& L,
Mat<typename T1::elem_type>& U, Mat<typename T1::elem_type>& U,
Mat<typename T1::elem_type>& P, Mat<typename T1::elem_type>& P,
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( ( (&L == &U) || (&L == &P) || (&U == &P) ), "lu(): two or more output objects are the same object"); arma_debug_check( ( (&L == &U) || (&L == &P) || (&U == &P) ), "lu(): two or more output objects are the same object");
auxlib::lu(L, U, P, X); const bool status = auxlib::lu(L, U, P, X);
if(status == false)
{
L.reset();
U.reset();
P.reset();
arma_bad("lu(): failed to converge", false);
}
return status;
} }
//! @} //! @}
 End of changes. 7 change blocks. 
8 lines changed or deleted 25 lines changed or added


 fn_misc.hpp   fn_misc.hpp 
skipping to change at line 154 skipping to change at line 154
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
const unwrap_cube<T1> tmp(X.get_ref()); const unwrap_cube<T1> tmp(X.get_ref());
const Cube<eT>& A = tmp.M; const Cube<eT>& A = tmp.M;
return A.is_finite(); return A.is_finite();
} }
template<typename T1>
arma_inline
Op<T1, op_sympd>
sympd(const Base<typename T1::elem_type,T1>& X)
{
arma_extra_debug_sigprint();
return Op<T1, op_sympd>(X.get_ref());
}
//! @} //! @}
 End of changes. 1 change blocks. 
0 lines changed or deleted 10 lines changed or added


 fn_pinv.hpp   fn_pinv.hpp 
// Copyright (C) 2009-2010 NICTA (www.nicta.com.au) // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
// Copyright (C) 2009-2010 Conrad Sanderson // Copyright (C) 2009-2011 Conrad Sanderson
// Copyright (C) 2009-2010 Dimitrios Bouzas // Copyright (C) 2009-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 28 skipping to change at line 28
inline inline
const Op<T1, op_pinv> const Op<T1, op_pinv>
pinv pinv
( (
const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T1>& X,
const typename T1::elem_type tol = 0.0, const typename T1::elem_type tol = 0.0,
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);
return Op<T1, op_pinv>(X.get_ref(), tol); return Op<T1, op_pinv>(X.get_ref(), tol);
} }
template<typename T1>
inline
bool
pinv
(
Mat<typename T1::elem_type>& out,
const Base<typename T1::elem_type,T1>& X,
const typename T1::elem_type tol = 0.0,
const typename arma_blas_type_only<typename T1::elem_type>::result* junk
= 0
)
{
arma_extra_debug_sigprint();
arma_ignore(junk);
try
{
out = pinv(X,tol);
}
catch(std::runtime_error&)
{
return false;
}
return true;
}
//! @} //! @}
 End of changes. 3 change blocks. 
3 lines changed or deleted 29 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


 fn_qr.hpp   fn_qr.hpp 
// Copyright (C) 2009-2010 NICTA (www.nicta.com.au) // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
// Copyright (C) 2009-2010 Conrad Sanderson // Copyright (C) 2009-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 fn_qr //! \addtogroup fn_qr
//! @{ //! @{
//! QR decomposition //! QR decomposition
template<typename T1> template<typename T1>
inline inline
void bool
qr qr
( (
Mat<typename T1::elem_type>& Q, Mat<typename T1::elem_type>& Q,
Mat<typename T1::elem_type>& R, Mat<typename T1::elem_type>& R,
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( (&Q == &R), "qr(): Q and R are the same object"); arma_debug_check( (&Q == &R), "qr(): Q and R are the same object");
const bool status = auxlib::qr(Q, R, X); const bool status = auxlib::qr(Q, R, X);
if(status == false) if(status == false)
{ {
arma_print("qr(): failed to converge");
Q.reset(); Q.reset();
R.reset(); R.reset();
arma_bad("qr(): failed to converge", false);
} }
return false;
} }
//! @} //! @}
 End of changes. 6 change blocks. 
5 lines changed or deleted 6 lines changed or added


 fn_rank.hpp   fn_rank.hpp 
// Copyright (C) 2009-2010 NICTA (www.nicta.com.au) // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
// Copyright (C) 2009-2010 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 29 skipping to change at line 30
arma_warn_unused arma_warn_unused
u32 u32
rank rank
( (
const Base<typename T1::elem_type,T1>& X, const Base<typename T1::elem_type,T1>& X,
typename T1::pod_type tol = 0.0, typename T1::pod_type tol = 0.0,
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);
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
typedef typename T1::pod_type T; typedef typename T1::pod_type T;
u32 X_n_rows; u32 X_n_rows;
u32 X_n_cols; u32 X_n_cols;
Col<T> s; Col<T> s;
const bool status = auxlib::svd(s, X, X_n_rows, X_n_cols); const bool status = auxlib::svd(s, X, X_n_rows, X_n_cols);
const u32 n_elem = s.n_elem;
if(status == true) if(status == true)
{ {
if(tol == T(0)) if( (tol == T(0)) && (n_elem > 0) )
{ {
tol = (std::max)(X_n_rows, X_n_cols) * eop_aux::direct_eps(max(s)); tol = (std::max)(X_n_rows, X_n_cols) * eop_aux::direct_eps(max(s));
} }
// count non zero valued elements in s // count non zero valued elements in s
const T* s_mem = s.memptr(); const T* s_mem = s.memptr();
const u32 n_elem = s.n_elem;
u32 count = 0; u32 count = 0;
for(u32 i=0; i<n_elem; ++i) for(u32 i=0; i<n_elem; ++i)
{ {
if(s_mem[i] > tol) if(s_mem[i] > tol)
{ {
++count; ++count;
} }
} }
return count; return count;
} }
else else
{ {
arma_print("rank(): failed to converge"); arma_bad("rank(): failed to converge");
return u32(0); return u32(0);
} }
} }
//! @} //! @}
 End of changes. 7 change blocks. 
6 lines changed or deleted 7 lines changed or added


 fn_solve.hpp   fn_solve.hpp 
skipping to change at line 64 skipping to change at line 64
arma_ignore(junk); arma_ignore(junk);
return Glue<T1, T2, glue_solve_tr>(A.m, B.get_ref(), A.aux_u32_a); return Glue<T1, T2, glue_solve_tr>(A.m, B.get_ref(), A.aux_u32_a);
} }
template<typename T1, typename T2> template<typename T1, typename T2>
inline inline
bool bool
solve solve
( (
Mat<typename T1::elem_type>& out, Mat<typename T1::elem_type>& out,
const Base<typename T1::elem_type,T1>& A, const Base<typename T1::elem_type,T1>& A,
const Base<typename T1::elem_type,T2>& B, const Base<typename T1::elem_type,T2>& B,
const bool slow = false, const bool slow = false,
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);
out = solve( A.get_ref(), B.get_ref(), slow ); try
{
out = solve( A.get_ref(), B.get_ref(), slow );
}
catch(std::runtime_error&)
{
return false;
}
return (out.n_elem == 0) ? false : true; return true;
} }
//! @} //! @}
 End of changes. 3 change blocks. 
3 lines changed or deleted 10 lines changed or added


 fn_svd.hpp   fn_svd.hpp 
// Copyright (C) 2009-2010 NICTA (www.nicta.com.au) // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
// Copyright (C) 2009-2010 Conrad Sanderson // Copyright (C) 2009-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)
skipping to change at line 27 skipping to change at line 27
inline inline
bool bool
svd svd
( (
Col<typename T1::pod_type>& S, Col<typename T1::pod_type>& S,
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);
// it doesn't matter if X is an alias of S, as auxlib::svd() makes a copy of X // it doesn't matter if X is an alias of S, as auxlib::svd() makes a copy of X
const bool status = auxlib::svd(S, X); const bool status = auxlib::svd(S, X);
if(status == false) if(status == false)
{ {
arma_print("svd(): failed to converge");
S.reset(); S.reset();
arma_bad("svd(): failed to converge", false);
} }
return status; return status;
} }
template<typename T1> template<typename T1>
inline inline
Col<typename T1::pod_type> Col<typename T1::pod_type>
svd svd
( (
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;
svd(out, X); const bool status = auxlib::svd(out, X);
if(status == false)
{
out.reset();
arma_bad("svd(): failed to converge");
}
return out; return out;
} }
template<typename T1> template<typename T1>
inline inline
bool bool
svd svd
( (
Mat<typename T1::elem_type>& U, Mat<typename T1::elem_type>& U,
Col<typename T1::pod_type >& S, Col<typename T1::pod_type >& S,
Mat<typename T1::elem_type>& V, Mat<typename T1::elem_type>& V,
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);
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
arma_debug_check( ( ((void*)(&U) == (void*)(&S)) || (&U == &V) || ((void* arma_debug_check
)(&S) == (void*)(&V)) ), "svd(): two or more output objects are the same ob (
ject" ); ( ((void*)(&U) == (void*)(&S)) || (&U == &V) || ((void*)(&S) == (void*)
(&V)) ),
"svd(): two or more output objects are the same object"
);
// auxlib::svd() makes an internal copy of X // auxlib::svd() makes an internal copy of X
const bool status = auxlib::svd(U, S, V, X); const bool status = auxlib::svd(U, S, V, X);
if(status == false) if(status == false)
{ {
arma_print("svd(): failed to converge");
U.reset(); U.reset();
S.reset(); S.reset();
V.reset(); V.reset();
arma_bad("svd(): failed to converge", false);
} }
return status; return status;
} }
//! @} //! @}
 End of changes. 10 change blocks. 
11 lines changed or deleted 17 lines changed or added


 fn_syl_lyap.hpp   fn_syl_lyap.hpp 
skipping to change at line 42 skipping to change at line 42
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
const unwrap_check<T1> tmp_A(in_A.get_ref(), out); const unwrap_check<T1> tmp_A(in_A.get_ref(), out);
const unwrap_check<T2> tmp_B(in_B.get_ref(), out); const unwrap_check<T2> tmp_B(in_B.get_ref(), out);
const unwrap_check<T3> tmp_C(in_C.get_ref(), out); const unwrap_check<T3> tmp_C(in_C.get_ref(), out);
const Mat<eT>& A = tmp_A.M; const Mat<eT>& A = tmp_A.M;
const Mat<eT>& B = tmp_B.M; const Mat<eT>& B = tmp_B.M;
const Mat<eT>& C = tmp_C.M; const Mat<eT>& C = tmp_C.M;
return auxlib::syl(out, A, B, C); const bool status = auxlib::syl(out, A, B, C);
if(status == false)
{
out.reset();
arma_bad("syl(): equation appears to be singular", false);
}
return false;
} }
template<typename T1, typename T2, typename T3> template<typename T1, typename T2, typename T3>
inline inline
Mat<typename T1::elem_type> Mat<typename T1::elem_type>
syl syl
( (
const Base<typename T1::elem_type,T1>& in_A, const Base<typename T1::elem_type,T1>& in_A,
const Base<typename T1::elem_type,T2>& in_B, const Base<typename T1::elem_type,T2>& in_B,
const Base<typename T1::elem_type,T3>& in_C, const Base<typename T1::elem_type,T3>& in_C,
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);
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
Mat<eT> out; Mat<eT> out;
syl(out, in_A, in_B, in_C); const bool status = syl(out, in_A, in_B, in_C);
if(status == false)
{
out.reset();
arma_bad("syl(): equation appears to be singular");
}
return out; return out;
} }
//! @} //! @}
 End of changes. 2 change blocks. 
2 lines changed or deleted 16 lines changed or added


 forward_bones.hpp   forward_bones.hpp 
skipping to change at line 109 skipping to change at line 109
static const injector_end_of_row endr = injector_end_of_row(); static const injector_end_of_row endr = injector_end_of_row();
//!< endr indicates "end of row" when using the << operator; //!< endr indicates "end of row" when using the << operator;
//!< similar conceptual meaning to std::endl //!< similar conceptual meaning to std::endl
//! @} //! @}
//! \addtogroup diskio //! \addtogroup diskio
//! @{ //! @{
//! file types supported by Armadillo
enum file_type enum file_type
{ {
file_type_unknown,
auto_detect, //!< Automatically detect the file type (file must be one o f the following types) auto_detect, //!< Automatically detect the file type (file must be one o f the following types)
raw_ascii, //!< ASCII format (text), without any other information. raw_ascii, //!< ASCII format (text), without any other information.
arma_ascii, //!< Armadillo ASCII format (text), with information about matrix type and size arma_ascii, //!< Armadillo ASCII format (text), with information about matrix type and size
csv_ascii, //!< comma separated values (CSV), without any other inform ation
raw_binary, //!< raw binary format, without any other information. raw_binary, //!< raw binary format, without any other information.
arma_binary, //!< Armadillo binary format, with information about matrix type and size arma_binary, //!< Armadillo binary format, with information about matrix type and size
pgm_binary, //!< Portable Grey Map (greyscale image) pgm_binary, //!< Portable Grey Map (greyscale image)
ppm_binary //!< Portable Pixel Map (colour image), used by the field a nd cube classes ppm_binary //!< Portable Pixel Map (colour image), used by the field a nd cube classes
}; };
//! @} //! @}
 End of changes. 3 change blocks. 
1 lines changed or deleted 2 lines changed or added


 glue_join_meat.hpp   glue_join_meat.hpp 
skipping to change at line 31 skipping to change at line 31
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
typedef typename T1::elem_type eT; typedef typename T1::elem_type eT;
const unwrap<T1> A_tmp(X.A); const unwrap<T1> A_tmp(X.A);
const unwrap<T2> B_tmp(X.B); const unwrap<T2> B_tmp(X.B);
const Mat<eT>& A = A_tmp.M; const Mat<eT>& A = A_tmp.M;
const Mat<eT>& B = B_tmp.M; const Mat<eT>& B = B_tmp.M;
const u32 A_n_rows = A.n_rows;
const u32 A_n_cols = A.n_cols;
const u32 B_n_rows = B.n_rows;
const u32 B_n_cols = B.n_cols;
const u32 join_type = X.aux_u32; const u32 join_type = X.aux_u32;
if(join_type == 0) if(join_type == 0)
{ {
arma_debug_check( (A.n_cols != B.n_cols), "join_cols(): number of colum arma_debug_check
ns must be the same" ); (
( (A_n_cols != B_n_cols) && ( (A_n_rows > 0) || (A_n_cols > 0) ) && (
(B_n_rows > 0) || (B_n_cols > 0) ) ),
"join_cols(): number of columns must be the same"
);
} }
else else
{ {
arma_debug_check( (A.n_rows != B.n_rows), "join_rows(): number of rows arma_debug_check
must be the same" ); (
( (A_n_rows != B.n_rows) && ( (A_n_rows > 0) || (A_n_cols > 0) ) && (
(B_n_rows > 0) || (B_n_cols > 0) ) ),
"join_rows(): number of rows must be the same"
);
} }
if( (&out != &A) && (&out != &B) ) if( (&out != &A) && (&out != &B) )
{ {
if(join_type == 0) // join columns (i.e. result matrix has more rows) if(join_type == 0) // join columns (i.e. result matrix has more rows)
{ {
out.set_size(A.n_rows + B.n_rows, A.n_cols); out.set_size( A_n_rows + B_n_rows, (std::max)(A_n_cols, B_n_cols) );
if( out.n_elem > 0 ) if( out.n_elem > 0 )
{ {
out.submat(0, 0, A.n_rows-1, out.n_cols-1) = A; if(A.is_empty() == false)
out.submat(A.n_rows, 0, out.n_rows-1, out.n_cols-1) = B; {
out.submat(0, 0, A_n_rows-1, out.n_cols-1) = A;
}
if(B.is_empty() == false)
{
out.submat(A_n_rows, 0, out.n_rows-1, out.n_cols-1) = B;
}
} }
} }
else // join rows (i.e. result matrix has more columns) else // join rows (i.e. result matrix has more columns)
{ {
out.set_size(A.n_rows, A.n_cols + B.n_cols); out.set_size( (std::max)(A_n_rows, B_n_rows), A_n_cols + B_n_cols );
if( out.n_elem > 0 ) if( out.n_elem > 0 )
{ {
out.submat(0, 0, out.n_rows-1, A.n_cols-1) = A; if(A.is_empty() == false)
out.submat(0, A.n_cols, out.n_rows-1, out.n_cols-1) = B; {
out.submat(0, 0, out.n_rows-1, A.n_cols-1) = A;
}
if(B.is_empty() == false)
{
out.submat(0, A_n_cols, out.n_rows-1, out.n_cols-1) = B;
}
} }
} }
} }
else // we have aliasing else // we have aliasing
{ {
Mat<eT> C; Mat<eT> C;
if(join_type == 0) if(join_type == 0)
{ {
C.set_size(A.n_rows + B.n_rows, A.n_cols); C.set_size( A_n_rows + B_n_rows, (std::max)(A_n_cols, B_n_cols) );
if( C.n_elem > 0 ) if( C.n_elem > 0 )
{ {
C.submat(0, 0, A.n_rows-1, C.n_cols-1) = A; if(A.is_empty() == false)
C.submat(A.n_rows, 0, C.n_rows-1, C.n_cols-1) = B; {
C.submat(0, 0, A_n_rows-1, C.n_cols-1) = A;
}
if(B.is_empty() == false)
{
C.submat(A_n_rows, 0, C.n_rows-1, C.n_cols-1) = B;
}
} }
} }
else else
{ {
C.set_size(A.n_rows, A.n_cols + B.n_cols); C.set_size( (std::max)(A_n_rows, B_n_rows), A_n_cols + B_n_cols );
if( C.n_elem > 0 ) if( C.n_elem > 0 )
{ {
C.submat(0, 0, C.n_rows-1, A.n_cols-1) = A; if(A.is_empty() == false)
C.submat(0, A.n_cols, C.n_rows-1, C.n_cols-1) = B; {
C.submat(0, 0, C.n_rows-1, A_n_cols-1) = A;
}
if(B.is_empty() == false)
{
C.submat(0, A_n_cols, C.n_rows-1, C.n_cols-1) = B;
}
} }
} }
out.steal_mem(C); out.steal_mem(C);
} }
} }
template<typename T1, typename T2> template<typename T1, typename T2>
inline inline
 End of changes. 11 change blocks. 
16 lines changed or deleted 58 lines changed or added


 glue_solve_meat.hpp   glue_solve_meat.hpp 
skipping to change at line 38 skipping to change at line 38
const Mat<eT>& B = B_tmp.M; const Mat<eT>& B = B_tmp.M;
arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in A a nd B must be the same" ); arma_debug_check( (A.n_rows != B.n_rows), "solve(): number of rows in A a nd B must be the same" );
bool status; bool status;
if(A.n_rows == A.n_cols) if(A.n_rows == A.n_cols)
{ {
const u32 mode = X.aux_u32; const u32 mode = X.aux_u32;
if(mode == 0) status = (mode == 0) ? auxlib::solve(out, A, B) : auxlib::solve(out, A,
{ B, true);
status = auxlib::solve(out, A, B);
}
else
{
status = auxlib::solve(out, A, B, true);
}
} }
else else
if(A.n_rows > A.n_cols) if(A.n_rows > A.n_cols)
{ {
arma_extra_debug_print("solve(): detected over-determined system"); arma_extra_debug_print("solve(): detected over-determined system");
status = auxlib::solve_od(out, A, B); status = auxlib::solve_od(out, A, B);
} }
else else
{ {
arma_extra_debug_print("solve(): detected under-determined system"); arma_extra_debug_print("solve(): detected under-determined system");
status = auxlib::solve_ud(out, A, B); status = auxlib::solve_ud(out, A, B);
} }
if(status == false) if(status == false)
{ {
out.reset(); out.reset();
arma_print("solve(): solution not found"); arma_bad("solve(): solution not found");
} }
} }
template<typename T1, typename T2> template<typename T1, typename T2>
inline inline
void void
glue_solve_tr::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glu e_solve_tr>& X) glue_solve_tr::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glu e_solve_tr>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 95 skipping to change at line 88
arma_debug_set_error( err_state, err_msg, (A.n_rows != B.n_rows), "solv e(): number of rows in A and B must be the same" ); arma_debug_set_error( err_state, err_msg, (A.n_rows != B.n_rows), "solv e(): number of rows in A and B must be the same" );
arma_debug_set_error( err_state, err_msg, (A.is_square() == false), "solv e(): A is not a square matrix" ); arma_debug_set_error( err_state, err_msg, (A.is_square() == false), "solv e(): A is not a square matrix" );
arma_debug_check(err_state, err_msg); arma_debug_check(err_state, err_msg);
const bool status = auxlib::solve_tr(out, A, B, X.aux_u32); const bool status = auxlib::solve_tr(out, A, B, X.aux_u32);
if(status == false) if(status == false)
{ {
out.reset(); out.reset();
arma_print("solve(): solution not found"); arma_bad("solve(): solution not found");
} }
} }
//! @} //! @}
 End of changes. 3 change blocks. 
10 lines changed or deleted 4 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_chol_meat.hpp   op_chol_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)
skipping to change at line 27 skipping to change at line 27
inline inline
void void
op_chol::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_chol>& X) op_chol::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_chol>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
const bool status = auxlib::chol(out, X.m); const bool status = auxlib::chol(out, X.m);
if(status == false) if(status == false)
{ {
arma_print("chol(): failed to converge");
out.reset(); out.reset();
arma_bad("chol(): failed to converge");
} }
} }
//! @} //! @}
 End of changes. 3 change blocks. 
3 lines changed or deleted 3 lines changed or added


 op_inv_bones.hpp   op_inv_bones.hpp 
skipping to change at line 40 skipping to change at line 40
//! 'invert matrix' operation (triangular matrices) //! 'invert matrix' operation (triangular matrices)
class op_inv_tr class op_inv_tr
{ {
public: public:
template<typename T1> template<typename T1>
inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op _inv_tr>& in); inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op _inv_tr>& in);
}; };
//! 'invert matrix' operation (positive definite symmetric matrices) //! 'invert matrix' operation (symmetric positive definite matrices)
class op_inv_sym class op_inv_sympd
{ {
public: public:
template<typename T1> template<typename T1>
inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op _inv_sym>& in); inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op _inv_sympd>& in);
}; };
//! @} //! @}
 End of changes. 2 change blocks. 
3 lines changed or deleted 3 lines changed or added


 op_inv_meat.hpp   op_inv_meat.hpp 
skipping to change at line 28 skipping to change at line 28
inline inline
void void
op_inv::apply(Mat<eT>& out, const Mat<eT>& A, const bool slow) op_inv::apply(Mat<eT>& out, const Mat<eT>& A, const bool slow)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
// no need to check for aliasing, due to: // no need to check for aliasing, due to:
// - auxlib::inv() copies A to out before inversion // - auxlib::inv() copies A to out before inversion
// - for 2x2 and 3x3 matrices the code is alias safe // - for 2x2 and 3x3 matrices the code is alias safe
auxlib::inv(out, A, slow); bool status = auxlib::inv(out, A, slow);
if(status == false)
{
out.reset();
arma_bad("inv(): matrix appears to be singular");
}
} }
//! immediate inverse of T1, storing the result in a dense matrix //! immediate inverse of T1, storing the result in a dense matrix
template<typename T1> template<typename T1>
inline inline
void void
op_inv::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv>& X) op_inv::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 51 skipping to change at line 57
const strip_diagmat<T1> strip(X.m); const strip_diagmat<T1> strip(X.m);
if(strip.do_diagmat == true) if(strip.do_diagmat == true)
{ {
op_inv::apply_diag(out, strip.M); op_inv::apply_diag(out, strip.M);
} }
else else
{ {
const u32 mode = X.aux_u32_a; const u32 mode = X.aux_u32_a;
if(mode == 0) const bool status = (mode == 0) ? auxlib::inv(out, X.m) : auxlib::inv(o
{ ut, X.m, true);
auxlib::inv(out, X.m);
} if(status == false)
else
{ {
auxlib::inv(out, X.m, true); out.reset();
arma_bad("inv(): matrix appears to be singular");
} }
} }
} }
template<typename T1> template<typename T1>
inline inline
void void
op_inv::apply_diag(Mat<typename T1::elem_type>& out, const Base<typename T1 ::elem_type, T1>& X) op_inv::apply_diag(Mat<typename T1::elem_type>& out, const Base<typename T1 ::elem_type, T1>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
skipping to change at line 96 skipping to change at line 101
} }
//! inverse of T1 (triangular matrices) //! inverse of T1 (triangular matrices)
template<typename T1> template<typename T1>
inline inline
void void
op_inv_tr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_tr>& X) op_inv_tr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_tr>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
auxlib::inv_tr(out, X.m, X.aux_u32_a); const bool status = auxlib::inv_tr(out, X.m, X.aux_u32_a);
if(status == false)
{
out.reset();
arma_bad("inv(): matrix appears to be singular");
}
} }
//! inverse of T1 (symmetric matrices) //! inverse of T1 (symmetric positive definite matrices)
template<typename T1> template<typename T1>
inline inline
void void
op_inv_sym::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_sym> & X) op_inv_sympd::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_inv_sy mpd>& X)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
auxlib::inv_sym(out, X.m, X.aux_u32_a); const bool status = auxlib::inv_sympd(out, X.m, X.aux_u32_a);
if(status == false)
{
out.reset();
arma_bad("inv(): matrix appears to be singular");
}
} }
//! @} //! @}
 End of changes. 7 change blocks. 
11 lines changed or deleted 29 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_bones.hpp   op_misc_bones.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)
skipping to change at line 49 skipping to change at line 49
{ {
public: public:
template<typename T1> template<typename T1>
inline static void apply( Mat<typename T1::pod_type>& out, const mtOp<typ ename T1::pod_type, T1, op_abs>& X); inline static void apply( Mat<typename T1::pod_type>& out, const mtOp<typ ename T1::pod_type, T1, op_abs>& X);
template<typename T1> template<typename T1>
inline static void apply( Cube<typename T1::pod_type>& out, const mtOpCub e<typename T1::pod_type, T1, op_abs>& X); inline static void apply( Cube<typename T1::pod_type>& out, const mtOpCub e<typename T1::pod_type, T1, op_abs>& X);
}; };
class op_sympd
{
public:
template<typename T1>
inline static void apply( Mat<typename T1::elem_type>& out, const Op<T1,
op_sympd>& X);
};
//! @} //! @}
 End of changes. 2 change blocks. 
2 lines changed or deleted 11 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_bones.hpp   op_princomp_bones.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 24 skipping to change at line 24
//! \addtogroup op_princomp //! \addtogroup op_princomp
//! @{ //! @{
class op_princomp class op_princomp
{ {
public: public:
// real element versions // real element versions
template<typename eT> template<typename eT>
inline static void inline static bool
direct_princomp direct_princomp
( (
Mat<eT>& coeff_out, Mat<eT>& coeff_out,
const Mat<eT>& in const Mat<eT>& in
); );
template<typename eT> template<typename eT>
inline static void inline static bool
direct_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
); );
template<typename eT> template<typename eT>
inline static void inline static bool
direct_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
); );
template<typename eT> template<typename eT>
inline static void inline static bool
direct_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
); );
// complex element versions // complex element versions
template<typename T> template<typename T>
inline static void inline static bool
direct_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
); );
template<typename T> template<typename T>
inline static void inline static bool
direct_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
); );
template<typename T> template<typename T>
inline static void inline static bool
direct_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
); );
template<typename T> template<typename T>
inline static void inline static bool
direct_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
); );
template<typename T1> template<typename T1>
 End of changes. 9 change blocks. 
10 lines changed or deleted 10 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


 op_prod_meat.hpp   op_prod_meat.hpp 
skipping to change at line 40 skipping to change at line 40
arma_debug_check( (dim > 1), "prod(): incorrect usage. dim must be 0 or 1 "); arma_debug_check( (dim > 1), "prod(): incorrect usage. dim must be 0 or 1 ");
const unwrap_check<T1> tmp(in.m, out); const unwrap_check<T1> tmp(in.m, out);
const Mat<eT>& X = tmp.M; const Mat<eT>& X = tmp.M;
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) // traverse across rows (i.e. find the product in each colu mn) if(dim == 0) // traverse across rows (i.e. find the product in each colu mn)
{ {
out.set_size((X_n_rows > 0) ? 1 : 0, X_n_cols); out.set_size(1, X_n_cols);
if(X_n_rows > 0) 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] = arrayops::product(X.colptr(col), X_n_rows); out_mem[col] = arrayops::product(X.colptr(col), X_n_rows);
}
} }
} }
else // traverse across columns (i.e. find the product in each row) else // traverse across columns (i.e. find the product in each row)
{ {
out.set_size(X_n_rows, (X_n_cols > 0) ? 1 : 0); out.set_size(X_n_rows, 1);
if(X_n_cols > 0) eT* out_mem = out.memptr();
for(u32 row=0; row<X_n_rows; ++row)
{ {
eT* out_mem = out.memptr(); eT val = eT(1);
for(u32 row=0; row < X_n_rows; ++row) for(u32 col=0; col<X_n_cols; ++col)
{ {
eT val = X.at(row,0); val *= X.at(row,col);
for(u32 col=1; col < X_n_cols; ++col)
{
val *= X.at(row,col);
}
out_mem[row] = val;
} }
out_mem[row] = val;
} }
} }
} }
//! @} //! @}
 End of changes. 9 change blocks. 
20 lines changed or deleted 14 lines changed or added


 op_stddev_meat.hpp   op_stddev_meat.hpp 
skipping to change at line 46 skipping to change at line 46
arma_debug_check( (norm_type > 1), "stddev(): incorrect usage. norm_type must be 0 or 1"); arma_debug_check( (norm_type > 1), "stddev(): incorrect usage. norm_type must be 0 or 1");
arma_debug_check( (dim > 1), "stddev(): incorrect usage. dim must b e 0 or 1" ); arma_debug_check( (dim > 1), "stddev(): incorrect usage. dim must b e 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_stddev::apply(), dim = 0"); arma_extra_debug_print("op_stddev::apply(), dim = 0");
out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); arma_debug_check( (X_n_rows == 0), "stddev(): given object has zero row s" );
if(X_n_rows > 0) out.set_size(1, X_n_cols);
{
out_eT* out_mem = out.memptr(); out_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] = std::sqrt( op_var::direct_var( X.colptr(col), X_n_ro out_mem[col] = std::sqrt( op_var::direct_var( X.colptr(col), X_n_rows
ws, norm_type ) ); , norm_type ) );
}
} }
} }
else else
if(dim == 1) if(dim == 1)
{ {
arma_extra_debug_print("op_stddev::apply(), dim = 1"); arma_extra_debug_print("op_stddev::apply(), dim = 1");
out.set_size( X_n_rows, (X_n_cols > 0) ? 1 : 0 ); arma_debug_check( (X_n_cols == 0), "stddev(): given object has zero col umns" );
if(X_n_cols > 0) out.set_size(X_n_rows, 1);
{
podarray<in_eT> tmp(X_n_cols); podarray<in_eT> tmp(X_n_cols);
in_eT* tmp_mem = tmp.memptr(); in_eT* tmp_mem = tmp.memptr();
out_eT* out_mem = out.memptr(); out_eT* out_mem = out.memptr();
for(u32 row=0; row<X_n_rows; ++row) for(u32 row=0; row<X_n_rows; ++row)
{ {
tmp.copy_row(X, row); tmp.copy_row(X, row);
out_mem[row] = std::sqrt( op_var::direct_var( tmp_mem, X_n_cols, no out_mem[row] = std::sqrt( op_var::direct_var( tmp_mem, X_n_cols, norm
rm_type) ); _type) );
}
} }
} }
} }
//! @} //! @}
 End of changes. 8 change blocks. 
21 lines changed or deleted 19 lines changed or added


 op_sum_meat.hpp   op_sum_meat.hpp 
skipping to change at line 40 skipping to change at line 40
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>& X = tmp.M; const Mat<eT>& X = tmp.M;
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) // traverse across rows (i.e. find the sum in each column) if(dim == 0) // traverse across rows (i.e. find the sum in each column)
{ {
out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); out.set_size(1, X_n_cols);
if(X_n_rows > 0) 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] = arrayops::accumulate( X.colptr(col), X_n_rows ); out_mem[col] = arrayops::accumulate(X.colptr(col), X_n_rows);
}
} }
} }
else // traverse across columns (i.e. find the sum in each row) else // traverse across columns (i.e. find the sum in each row)
{ {
out.set_size( X_n_rows, (X_n_cols > 0) ? 1 : 0 ); out.set_size(X_n_rows, 1);
if(X_n_cols > 0) eT* out_mem = out.memptr();
for(u32 row=0; row<X_n_rows; ++row)
{ {
eT* out_mem = out.memptr(); eT val = eT(0);
for(u32 row=0; row<X_n_rows; ++row) for(u32 col=0; col<X_n_cols; ++col)
{ {
eT val = eT(0); val += X.at(row,col);
for(u32 col=0; col<X_n_cols; ++col)
{
val += X.at(row,col);
}
out_mem[row] = val;
} }
out_mem[row] = val;
} }
} }
} }
//! @} //! @}
 End of changes. 9 change blocks. 
20 lines changed or deleted 14 lines changed or added


 op_var_meat.hpp   op_var_meat.hpp 
skipping to change at line 187 skipping to change at line 187
arma_debug_check( (norm_type > 1), "var(): incorrect usage. norm_type mus t be 0 or 1"); arma_debug_check( (norm_type > 1), "var(): incorrect usage. norm_type mus t be 0 or 1");
arma_debug_check( (dim > 1), "var(): incorrect usage. dim must be 0 or 1" ); arma_debug_check( (dim > 1), "var(): 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_var::apply(), dim = 0"); arma_extra_debug_print("op_var::apply(), dim = 0");
out.set_size( (X_n_rows > 0) ? 1 : 0, X_n_cols ); arma_debug_check( (X_n_rows == 0), "var(): given object has zero rows" );
if(X_n_rows > 0) out.set_size(1, X_n_cols);
{
out_eT* out_mem = out.memptr(); out_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_var::direct_var( X.colptr(col), X_n_rows, norm_ty out_mem[col] = op_var::direct_var( X.colptr(col), X_n_rows, norm_type
pe ); );
}
} }
} }
else else
if(dim == 1) if(dim == 1)
{ {
arma_extra_debug_print("op_var::apply(), dim = 1"); arma_extra_debug_print("op_var::apply(), dim = 1");
out.set_size( X_n_rows, (X_n_cols > 0) ? 1 : 0 ); arma_debug_check( (X_n_cols == 0), "var(): given object has zero column s" );
if(X_n_cols > 0) out.set_size(X_n_rows, 1);
{
podarray<in_eT> tmp(X_n_cols); podarray<in_eT> tmp(X_n_cols);
in_eT* tmp_mem = tmp.memptr(); in_eT* tmp_mem = tmp.memptr();
out_eT* out_mem = out.memptr(); out_eT* out_mem = out.memptr();
for(u32 row=0; row<X_n_rows; ++row) for(u32 row=0; row<X_n_rows; ++row)
{ {
tmp.copy_row(X, row); tmp.copy_row(X, row);
out_mem[row] = op_var::direct_var( tmp_mem, X_n_cols, norm_type ); out_mem[row] = op_var::direct_var( tmp_mem, X_n_cols, norm_type );
}
} }
} }
} }
//! find the variance of an array (robust but slow) //! find the variance of an array (robust but slow)
template<typename eT> template<typename eT>
inline inline
eT eT
op_var::direct_var_robust(const eT* const X, const u32 n_elem, const u32 no rm_type) op_var::direct_var_robust(const eT* const X, const u32 n_elem, const u32 no rm_type)
{ {
 End of changes. 8 change blocks. 
20 lines changed or deleted 18 lines changed or added


 running_stat_meat.hpp   running_stat_meat.hpp 
skipping to change at line 139 skipping to change at line 139
//! update statistics to reflect new sample //! update statistics to reflect new sample
template<typename eT> template<typename eT>
inline inline
void void
running_stat<eT>::operator() (const typename running_stat<eT>::T sample) running_stat<eT>::operator() (const typename running_stat<eT>::T sample)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
if( arma_isfinite(sample) == false ) if( arma_isfinite(sample) == false )
{ {
arma_print("running_stat: sample ignored as it is non-finite" ); arma_warn(true, "running_stat: sample ignored as it is non-finite" );
return; return;
} }
running_stat_aux::update_stats(*this, sample); running_stat_aux::update_stats(*this, sample);
} }
//! update statistics to reflect new sample (version for complex numbers) //! update statistics to reflect new sample (version for complex numbers)
template<typename eT> template<typename eT>
inline inline
void void
running_stat<eT>::operator() (const std::complex< typename running_stat<eT> ::T >& sample) running_stat<eT>::operator() (const std::complex< typename running_stat<eT> ::T >& sample)
{ {
arma_extra_debug_sigprint(); arma_extra_debug_sigprint();
arma_type_check< is_same_type<eT, std::complex< typename running_stat<eT> ::T > >::value == false >::apply(); arma_type_check< is_same_type<eT, std::complex< typename running_stat<eT> ::T > >::value == false >::apply();
if( arma_isfinite(sample) == false ) if( arma_isfinite(sample) == false )
{ {
arma_print("running_stat: sample ignored as it is non-finite" ); arma_warn(true, "running_stat: sample ignored as it is non-finite" );
return; return;
} }
running_stat_aux::update_stats(*this, sample); running_stat_aux::update_stats(*this, sample);
} }
//! set all statistics to zero //! set all statistics to zero
template<typename eT> template<typename eT>
inline inline
void void
 End of changes. 2 change blocks. 
2 lines changed or deleted 2 lines changed or added


 running_stat_vec_meat.hpp   running_stat_vec_meat.hpp 
skipping to change at line 86 skipping to change at line 86
const unwrap<T1> tmp(X.get_ref()); const unwrap<T1> tmp(X.get_ref());
const Mat<eT>& sample = tmp.M; const Mat<eT>& sample = tmp.M;
if( sample.is_empty() ) if( sample.is_empty() )
{ {
return; return;
} }
if( sample.is_finite() == false ) if( sample.is_finite() == false )
{ {
arma_print("running_stat_vec: sample ignored as it has non-finite eleme nts"); arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements");
return; return;
} }
running_stat_vec_aux::update_stats(*this, sample); running_stat_vec_aux::update_stats(*this, sample);
} }
//! update statistics to reflect new sample (version for complex numbers) //! update statistics to reflect new sample (version for complex numbers)
template<typename eT> template<typename eT>
template<typename T1> template<typename T1>
arma_hot arma_hot
skipping to change at line 115 skipping to change at line 115
const unwrap<T1> tmp(X.get_ref()); const unwrap<T1> tmp(X.get_ref());
const Mat<eT>& sample = tmp.M; const Mat<eT>& sample = tmp.M;
if( sample.is_empty() ) if( sample.is_empty() )
{ {
return; return;
} }
if( sample.is_finite() == false ) if( sample.is_finite() == false )
{ {
arma_print("running_stat_vec: sample ignored as it has non-finite eleme nts"); arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements");
return; return;
} }
running_stat_vec_aux::update_stats(*this, sample); running_stat_vec_aux::update_stats(*this, sample);
} }
//! set all statistics to zero //! set all statistics to zero
template<typename eT> template<typename eT>
inline inline
void void
 End of changes. 2 change blocks. 
2 lines changed or deleted 2 lines changed or added

This html diff was produced by rfcdiff 1.41. The latest version is available from http://tools.ietf.org/tools/rfcdiff/