C++ Mat::colptr方法代码示例

inline
void
op_trimat::apply_htrans
  (
        Mat<eT>& out,
  const Mat<eT>& A,
  const bool     upper,
  const typename arma_cx_only<eT>::result* junk
  )
  {
  arma_extra_debug_sigprint();
  arma_ignore(junk);
  
  arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" );
  
  const uword N = A.n_rows;
  
  if(&out != &A)
    {
    out.copy_size(A);
    }
  
  if(upper)
    {
    // Upper triangular: but since we're transposing, we're taking the lower
    // triangular and putting it in the upper half.
    for(uword row = 0; row < N; ++row)
      {
      eT* out_colptr = out.colptr(row);
      
      for(uword col = 0; col <= row; ++col)
        {
        //out.at(col, row) = std::conj( A.at(row, col) );
        out_colptr[col] = std::conj( A.at(row, col) );
        }
      }
    }
  else
    {
    // Lower triangular: but since we're transposing, we're taking the upper
    // triangular and putting it in the lower half.
    for(uword row = 0; row < N; ++row)
      {
      for(uword col = row; col < N; ++col)
        {
        out.at(col, row) = std::conj( A.at(row, col) );
        }
      }
    }
  
  op_trimat::fill_zeros(out, upper);
  }

inline
void
op_var::apply(Mat< typename get_pod_type<eT>::result >& out, const Mat<eT>& X, const u32 norm_type, const u32 dim)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_check( (X.n_elem == 0), "var(): given matrix has no elements" );
  
  arma_debug_check( (norm_type > 1), "var(): incorrect usage. norm_type must be 0 or 1");
  arma_debug_check( (dim > 1),       "var(): incorrect usage. dim must be 0 or 1"      );
  
  
  if(dim == 0)
    {
    arma_extra_debug_print("op_var::apply(), dim = 0");
    
    out.set_size(1, X.n_cols);
    
    for(u32 col=0; col<X.n_cols; ++col)
      {
      out[col] = op_var::direct_var( X.colptr(col), X.n_rows, norm_type );
      }
    }
  else
  if(dim == 1)
    {
    arma_extra_debug_print("op_var::apply(), dim = 1");
    
    const u32 n_rows = X.n_rows;
    const u32 n_cols = X.n_cols;
    
    out.set_size(n_rows, 1);
    
    podarray<eT> tmp(n_cols);
    
    eT* tmp_mem = tmp.memptr();
    
    for(u32 row=0; row<n_rows; ++row)
      {
      for(u32 col=0; col<n_cols; ++col)
        {
        tmp_mem[col] = X.at(row,col);
        }
      
      out[row] = op_var::direct_var(tmp_mem, n_cols, norm_type);
      }
    
    }
  
  }

inline
void
op_prod::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim)
  {
  arma_extra_debug_sigprint();
  
  const uword X_n_rows = X.n_rows;
  const uword X_n_cols = X.n_cols;
    
  if(dim == 0)  // traverse across rows (i.e. find the product in each column)
    {
    out.set_size(1, X_n_cols);
    
    eT* out_mem = out.memptr();
    
    for(uword col=0; col < X_n_cols; ++col)
      {
      out_mem[col] = arrayops::product(X.colptr(col), X_n_rows);
      }
    }
  else  // traverse across columns (i.e. find the product in each row)
    {
    out.ones(X_n_rows, 1);
    
    eT* out_mem = out.memptr();
    
    for(uword col=0; col < X_n_cols; ++col)
      {
      const eT* X_col_mem = X.colptr(col);
      
      for(uword row=0; row < X_n_rows; ++row)
        {
        out_mem[row] *= X_col_mem[row];
        }
      }
    }
  }

arma_hot
inline
void
op_strans::apply(Mat<eT>& out, const TA& A)
  {
  arma_extra_debug_sigprint();
  
  if(&out != &A)
    {
    op_strans::apply_noalias(out, A);
    }
  else
    {
    const uword n_rows = A.n_rows;
    const uword n_cols = A.n_cols;
    
    if(n_rows == n_cols)
      {
      arma_extra_debug_print("op_strans::apply(): doing in-place transpose of a square matrix");
      
      const uword N = n_rows;
      
      for(uword k=0; k < N; ++k)
        {
        eT* colptr = out.colptr(k);
        
        uword i,j;
        
        for(i=(k+1), j=(k+2); j < N; i+=2, j+=2)
          {
          std::swap(out.at(k,i), colptr[i]);
          std::swap(out.at(k,j), colptr[j]);
          }
        
        if(i < N)
          {
          std::swap(out.at(k,i), colptr[i]);
          }
        }
      }
    else
      {
      Mat<eT> tmp;
      op_strans::apply_noalias(tmp, A);
      
      out.steal_mem(tmp);
      }
    }
  }

 arma_hot
 inline
 static
 void
 apply
   (
         Mat<eT>& C,
   const TA&      A,
   const TB&      B,
   const eT       alpha = eT(1),
   const eT       beta  = eT(0)
   )
   {
   arma_extra_debug_sigprint();
   
   switch(A.n_rows)
     {
     case  4:  gemv_emul_tinysq<do_trans_A, use_alpha, use_beta>::apply( C.colptr(3), A, B.colptr(3), alpha, beta );
     case  3:  gemv_emul_tinysq<do_trans_A, use_alpha, use_beta>::apply( C.colptr(2), A, B.colptr(2), alpha, beta );
     case  2:  gemv_emul_tinysq<do_trans_A, use_alpha, use_beta>::apply( C.colptr(1), A, B.colptr(1), alpha, beta );
     case  1:  gemv_emul_tinysq<do_trans_A, use_alpha, use_beta>::apply( C.colptr(0), A, B.colptr(0), alpha, beta );
     default:  ;
     }
   }

inline
void
subview_cube<eT>::schur_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_assert_same_size(out, in, "matrix schur product");
  
  const u32 in_n_rows     = in.n_rows;
  const u32 in_n_cols     = in.n_cols;
  const u32 in_aux_slice1 = in.aux_slice1;
  
  for(u32 col = 0; col < in_n_cols; ++col)
    {
    arrayops::inplace_mul( out.colptr(col), in.slice_colptr(in_aux_slice1, col), in_n_rows );
    }
  }

inline
void
op_htrans::apply(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A)
  {
  arma_extra_debug_sigprint();
  
  typedef typename std::complex<T> eT;

  if(&out != &A)
    {
    op_htrans::apply_noalias(out, A);
    }
  else
    {
    if(out.n_rows == out.n_cols)
      {
      arma_extra_debug_print("doing in-place hermitian transpose of a square matrix");
      
      const u32 n_rows = out.n_rows;
      const u32 n_cols = out.n_cols;
      
      for(u32 col=0; col<n_cols; ++col)
        {
        eT* coldata = out.colptr(col);
        
        out.at(col,col) = std::conj( out.at(col,col) );
        
        for(u32 row=(col+1); row<n_rows; ++row)
          {
          const eT val1 = std::conj(coldata[row]);
          const eT val2 = std::conj(out.at(col,row));
          
          out.at(col,row) = val1;
          coldata[row]    = val2;
          }
        }
      }
    else
      {
      const Mat<eT> A_copy = A;
      op_htrans::apply_noalias(out, A_copy);
      }
    }
  
  }

static inline void pwhichmax_core(Mat<double> posteriors, Vec<int> clusters, int nthreads){
	int ncol = posteriors.ncol;
	int nrow = posteriors.nrow;
	#pragma omp parallel for num_threads(nthreads)
	for (int col = 0; col < ncol; ++col){
		double* postCol = posteriors.colptr(col);
		double max = *postCol;
		int whichmax = 1;
		for (int row = 1; row < nrow; ++row){
			++postCol;
			if (*postCol > max){
				max = *postCol;
				whichmax = row + 1;
			}
		}
		clusters[col] = whichmax;
	}
}

inline
void
op_min::apply_noalias(Mat<eT>& out, const Mat<eT>& X, const uword dim, const typename arma_cx_only<eT>::result* junk)
  {
  arma_extra_debug_sigprint();
  arma_ignore(junk);
  
  const uword X_n_rows = X.n_rows;
  const uword X_n_cols = X.n_cols;
  
  if(dim == 0)
    {
    arma_extra_debug_print("op_min::apply(): dim = 0");
    
    out.set_size((X_n_rows > 0) ? 1 : 0, X_n_cols);
    
    if(X_n_rows == 0)  { return; }
    
    eT* out_mem = out.memptr();
    
    for(uword col=0; col<X_n_cols; ++col)
      {
      out_mem[col] = op_min::direct_min( X.colptr(col), X_n_rows );
      }
    }
  else
  if(dim == 1)
    {
    arma_extra_debug_print("op_min::apply(): dim = 1");
    
    out.set_size(X_n_rows, (X_n_cols > 0) ? 1 : 0);
    
    if(X_n_cols == 0)  { return; }
    
    eT* out_mem = out.memptr();
    
    for(uword row=0; row<X_n_rows; ++row)
      {
      out_mem[row] = op_min::direct_min( X, row );
      }
    }
  }

inline
void
subview_cube<eT>::div_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_assert_same_size(out, in, "matrix element-wise division");
  
  for(u32 col = 0; col < in.n_cols; ++col)
    {
    const eT* in_coldata  = in.slice_colptr(in.aux_slice1, col);
          eT* out_coldata = out.colptr(col);
    
    for(u32 row = 0; row < in.n_rows; ++row)
      {
      out_coldata[row] /= in_coldata[row];
      }
    }
    
  }

inline
void
op_fliplr::apply_proxy_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const uword P_n_rows = P.get_n_rows();
  const uword P_n_cols = P.get_n_cols();
  
  const uword P_n_cols_m1 = P_n_cols - 1;
  
  out.set_size(P_n_rows, P_n_cols);
  
  if( ((T1::is_row) || (P_n_rows == 1)) && (Proxy<T1>::use_at == false) )
    {
    eT* out_mem = out.memptr();
    
    const typename Proxy<T1>::ea_type P_ea = P.get_ea();
    
    for(uword col=0; col < P_n_cols; ++col)
      {
      out_mem[P_n_cols_m1 - col] = P_ea[col];
      }
    }
  else
    {
    for(uword col=0; col < P_n_cols; ++col)
      {
      eT* out_colmem = out.colptr(P_n_cols_m1 - col);
      
      for(uword row=0; row < P_n_rows; ++row)
        {
        out_colmem[row] = P.at(row,col);
        }
      }
    }
  }

inline
void
glue_toeplitz::apply(Mat<typename T1::elem_type>& out, const Glue<T1, T2, glue_toeplitz>& in)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  if( ((void*)(&in.A)) == ((void*)(&in.B)) )
    {
    arma_extra_debug_print("glue_toeplitz::apply(): one argument version");
    
    const unwrap_check<T1>  tmp(in.A, out);
    const Mat<eT>& A      = tmp.M;
    
    arma_debug_check( (A.is_vec() == false), "toeplitz(): input argument must be a vector" );
    
    const u32 N     = A.n_elem;
    const eT* A_mem = A.memptr();
    
    out.set_size(N,N);
    
    for(u32 col=0; col<N; ++col)
      {
      eT* col_mem = out.colptr(col);
      
      u32 i;
      
      i = col;
      for(u32 row=0; row<col; ++row, --i)
        {
        col_mem[row] = A_mem[i];
        }
      
      i = 0;
      for(u32 row=col; row<N; ++row, ++i)
        {
        col_mem[row] = A_mem[i];
        }      
      }
    }
  else
    {
    arma_extra_debug_print("glue_toeplitz::apply(): two argument version");
    
    const unwrap_check<T1> tmp1(in.A, out);
    const unwrap_check<T2> tmp2(in.B, out);
    
    const Mat<eT>& A = tmp1.M;
    const Mat<eT>& B = tmp2.M;
    
    arma_debug_check( ( (A.is_vec() == false) || (B.is_vec() == false) ), "toeplitz(): input arguments must be vectors" );
    
    const u32 A_N = A.n_elem;
    const u32 B_N = B.n_elem;
    
    const eT* A_mem = A.memptr();
    const eT* B_mem = B.memptr();
    
    out.set_size(A_N, B_N);
    
    for(u32 col=0; col<B_N; ++col)
      {
      eT* col_mem = out.colptr(col);
      
      u32 i = 0;
      for(u32 row=col; row<A_N; ++row, ++i)
        {
        col_mem[row] = A_mem[i];
        }
      }
    
    for(u32 row=0; row<A_N; ++row)
      {
      u32 i = 1;
      for(u32 col=(row+1); col<B_N; ++col, ++i)
        {
        out.at(row,col) = B_mem[i];
        }
      }
    
    }
  
  
  }

inline
void
subview_cube<eT>::div_inplace(Mat<eT>& out, const subview_cube<eT>& in)
  {
  arma_extra_debug_sigprint();
  
  arma_debug_assert_cube_as_mat(out, in, "element-wise division", true);
  
  const uword in_n_rows   = in.n_rows;
  const uword in_n_cols   = in.n_cols;
  const uword in_n_slices = in.n_slices;
  
  const uword out_n_rows    = out.n_rows;
  const uword out_n_cols    = out.n_cols;
  const uword out_vec_state = out.vec_state;
  
  if(in_n_slices == 1)
    {
    for(uword col=0; col < in_n_cols; ++col)
      {
      arrayops::inplace_div( out.colptr(col), in.slice_colptr(0, col), in_n_rows );
      }
    }
  else
    {
    if(out_vec_state == 0)
      {
      if( (in_n_rows == out_n_rows) && (in_n_cols == 1) && (in_n_slices == out_n_cols) )
        {
        for(uword i=0; i < in_n_slices; ++i)
          {
          arrayops::inplace_div( out.colptr(i), in.slice_colptr(i, 0), in_n_rows );
          }
        }
      else
      if( (in_n_rows == 1) && (in_n_cols == out_n_rows) && (in_n_slices == out_n_cols) )
        {
        const Cube<eT>& Q = in.m;
        
        const uword in_aux_row1   = in.aux_row1;
        const uword in_aux_col1   = in.aux_col1;
        const uword in_aux_slice1 = in.aux_slice1;
        
        for(uword slice=0; slice < in_n_slices; ++slice)
          {
          const uword mod_slice = in_aux_slice1 + slice;
          
          eT* out_colptr = out.colptr(slice);
          
          uword i,j;
          for(i=0, j=1; j < in_n_cols; i+=2, j+=2)
            {
            const eT tmp_i = Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
            const eT tmp_j = Q.at(in_aux_row1, in_aux_col1 + j, mod_slice);
            
            out_colptr[i] /= tmp_i;
            out_colptr[j] /= tmp_j;
            }
          
          if(i < in_n_cols)
            {
            out_colptr[i] /= Q.at(in_aux_row1, in_aux_col1 + i, mod_slice);
            }
          }
        }
      }
    else
      {
      eT* out_mem = out.memptr();
      
      const Cube<eT>& Q = in.m;
      
      const uword in_aux_row1   = in.aux_row1;
      const uword in_aux_col1   = in.aux_col1;
      const uword in_aux_slice1 = in.aux_slice1;
      
      for(uword i=0; i<in_n_slices; ++i)
        {
        out_mem[i] /= Q.at(in_aux_row1, in_aux_col1, in_aux_slice1 + i);
        }
      }
    }
  }

inline
void
op_fft_cx::apply_noalias(Mat<typename T1::elem_type>& out, const Proxy<T1>& P, const uword a, const uword b)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const uword n_rows = P.get_n_rows();
  const uword n_cols = P.get_n_cols();
  const uword n_elem = P.get_n_elem();
  
  const bool is_vec = ( (n_rows == 1) || (n_cols == 1) );
  
  const uword N_orig = (is_vec) ? n_elem : n_rows;
  const uword N_user = (b == 0) ? a      : N_orig;
  
  fft_engine<eT,inverse> worker(N_user);
  
  if(is_vec)
    {
    (n_cols == 1) ? out.set_size(N_user, 1) : out.set_size(1, N_user);
    
    if( (out.n_elem == 0) || (N_orig == 0) )
      {
      out.zeros();
      return;
      }
    
    if( (N_user == 1) && (N_orig >= 1) )
      {
      out[0] = P[0];
      return;
      }
    
    if( (N_user > N_orig) || (is_Mat<typename Proxy<T1>::stored_type>::value == false) )
      {
      podarray<eT> data(N_user);
      
      eT* data_mem = data.memptr();
      
      if(N_user > N_orig)  { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
      
      op_fft_cx::copy_vec( data_mem, P, (std::min)(N_user, N_orig) );
      
      worker.run( out.memptr(), data_mem );
      }
    else
      {
      const unwrap< typename Proxy<T1>::stored_type > tmp(P.Q);
      
      worker.run( out.memptr(), tmp.M.memptr() );
      }
    }
  else
    {
    // process each column seperately
    
    out.set_size(N_user, n_cols);
    
    if( (out.n_elem == 0) || (N_orig == 0) )
      {
      out.zeros();
      return;
      }
    
    if( (N_user == 1) && (N_orig >= 1) )
      {
      for(uword col=0; col < n_cols; ++col)  { out.at(0,col) = P.at(0,col); }
      
      return;
      }
    
    if( (N_user > N_orig) || (is_Mat<typename Proxy<T1>::stored_type>::value == false) )
      {
      podarray<eT> data(N_user);
      
      eT* data_mem = data.memptr();
      
      if(N_user > N_orig)  { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
      
      const uword N = (std::min)(N_user, N_orig);
      
      for(uword col=0; col < n_cols; ++col)
        {
        for(uword i=0; i < N; ++i)  { data_mem[i] = P.at(i, col); }
        
        worker.run( out.colptr(col), data_mem );
        }
      }
    else
      {
      const unwrap< typename Proxy<T1>::stored_type > tmp(P.Q);
      
      for(uword col=0; col < n_cols; ++col)
        {
        worker.run( out.colptr(col), tmp.M.colptr(col) );
        }
      }
    }
//.........这里部分代码省略.........

inline
void
op_fft_real::apply( Mat< std::complex<typename T1::pod_type> >& out, const mtOp<std::complex<typename T1::pod_type>,T1,op_fft_real>& in )
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::pod_type         in_eT;
  typedef typename std::complex<in_eT> out_eT;
  
  const Proxy<T1> P(in.m);
  
  const uword n_rows = P.get_n_rows();
  const uword n_cols = P.get_n_cols();
  const uword n_elem = P.get_n_elem();
  
  const bool is_vec = ( (n_rows == 1) || (n_cols == 1) );
  
  const uword N_orig = (is_vec)              ? n_elem         : n_rows;
  const uword N_user = (in.aux_uword_b == 0) ? in.aux_uword_a : N_orig;
  
  fft_engine<out_eT,false> worker(N_user);
  
  // no need to worry about aliasing, as we're going from a real object to complex complex, which by definition cannot alias
  
  if(is_vec)
    {
    (n_cols == 1) ? out.set_size(N_user, 1) : out.set_size(1, N_user);
    
    if( (out.n_elem == 0) || (N_orig == 0) )
      {
      out.zeros();
      return;
      }
    
    if( (N_user == 1) && (N_orig >= 1) )
      {
      out[0] = out_eT( P[0] );
      return;
      }
    
    podarray<out_eT> data(N_user);
    
    out_eT* data_mem = data.memptr();
    
    if(N_user > N_orig)  { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
    
    const uword N = (std::min)(N_user, N_orig);
    
    if(Proxy<T1>::use_at == false)
      {
      typename Proxy<T1>::ea_type X = P.get_ea();
      
      for(uword i=0; i < N; ++i)  { data_mem[i] = out_eT( X[i], in_eT(0) ); }
      }
    else
      {
      if(n_cols == 1)
        {
        for(uword i=0; i < N; ++i)  { data_mem[i] = out_eT( P.at(i,0), in_eT(0) ); }
        }
      else
        {
        for(uword i=0; i < N; ++i)  { data_mem[i] = out_eT( P.at(0,i), in_eT(0) ); }
        }
      }
    
    worker.run( out.memptr(), data_mem );
    }
  else
    {
    // process each column seperately
    
    out.set_size(N_user, n_cols);
    
    if( (out.n_elem == 0) || (N_orig == 0) )
      {
      out.zeros();
      return;
      }
    
    if( (N_user == 1) && (N_orig >= 1) )
      {
      for(uword col=0; col < n_cols; ++col)  { out.at(0,col) = out_eT( P.at(0,col) ); }
      
      return;
      }
    
    podarray<out_eT> data(N_user);
    
    out_eT* data_mem = data.memptr();
    
    if(N_user > N_orig)  { arrayops::fill_zeros( &data_mem[N_orig], (N_user - N_orig) ); }
    
    const uword N = (std::min)(N_user, N_orig);
    
    for(uword col=0; col < n_cols; ++col)
      {
      for(uword i=0; i < N; ++i)  { data_mem[i] = P.at(i, col); }
      
      worker.run( out.colptr(col), data_mem );
//.........这里部分代码省略.........