Chris@49: // Copyright (C) 2011-2012 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2011-2012 Conrad Sanderson
Chris@49: // 
Chris@49: // This Source Code Form is subject to the terms of the Mozilla Public
Chris@49: // License, v. 2.0. If a copy of the MPL was not distributed with this
Chris@49: // file, You can obtain one at http://mozilla.org/MPL/2.0/.
Chris@49: 
Chris@49: 
Chris@49: //! \addtogroup op_symmat
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: void
Chris@49: op_symmat::apply
Chris@49:   (
Chris@49:         Mat<typename T1::elem_type>& out,
Chris@49:   const Op<T1,op_symmat>&            in,
Chris@49:   const typename arma_not_cx<typename T1::elem_type>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const unwrap<T1>   tmp(in.m);
Chris@49:   const Mat<eT>& A = tmp.M;
Chris@49:   
Chris@49:   arma_debug_check( (A.is_square() == false), "symmatu()/symmatl(): given matrix must be square" );
Chris@49:   
Chris@49:   const uword N     = A.n_rows;
Chris@49:   const bool  upper = (in.aux_uword_a == 0);
Chris@49:   
Chris@49:   if(&out != &A)
Chris@49:     {
Chris@49:     out.copy_size(A);
Chris@49:     
Chris@49:     if(upper)
Chris@49:       {
Chris@49:       // upper triangular: copy the diagonal and the elements above the diagonal
Chris@49:       
Chris@49:       for(uword i=0; i<N; ++i)
Chris@49:         {
Chris@49:         const eT* A_data   = A.colptr(i);
Chris@49:               eT* out_data = out.colptr(i);
Chris@49:         
Chris@49:         arrayops::copy( out_data, A_data, i+1 );
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       // lower triangular: copy the diagonal and the elements below the diagonal
Chris@49:       
Chris@49:       for(uword i=0; i<N; ++i)
Chris@49:         {
Chris@49:         const eT* A_data   = A.colptr(i);
Chris@49:               eT* out_data = out.colptr(i);
Chris@49:         
Chris@49:         arrayops::copy( &out_data[i], &A_data[i], N-i );
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   if(upper)
Chris@49:     {
Chris@49:     // reflect elements across the diagonal from upper triangle to lower triangle
Chris@49:     
Chris@49:     for(uword col=1; col < N; ++col)
Chris@49:       {
Chris@49:       const eT* coldata = out.colptr(col);
Chris@49:       
Chris@49:       for(uword row=0; row < col; ++row)
Chris@49:         {
Chris@49:         out.at(col,row) = coldata[row];
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // reflect elements across the diagonal from lower triangle to upper triangle
Chris@49:     
Chris@49:     for(uword col=0; col < N; ++col)
Chris@49:       {
Chris@49:       const eT* coldata = out.colptr(col);
Chris@49:       
Chris@49:       for(uword row=(col+1); row < N; ++row)
Chris@49:         {
Chris@49:         out.at(col,row) = coldata[row];
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: void
Chris@49: op_symmat::apply
Chris@49:   (
Chris@49:         Mat<typename T1::elem_type>& out,
Chris@49:   const Op<T1,op_symmat>&            in,
Chris@49:   const typename arma_cx_only<typename T1::elem_type>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const unwrap<T1>   tmp(in.m);
Chris@49:   const Mat<eT>& A = tmp.M;
Chris@49:   
Chris@49:   arma_debug_check( (A.is_square() == false), "symmatu()/symmatl(): given matrix must be square" );
Chris@49:   
Chris@49:   const uword N     = A.n_rows;
Chris@49:   const bool  upper = (in.aux_uword_a == 0);
Chris@49:   
Chris@49:   if(&out != &A)
Chris@49:     {
Chris@49:     out.copy_size(A);
Chris@49:     
Chris@49:     if(upper)
Chris@49:       {
Chris@49:       // upper triangular: copy the diagonal and the elements above the diagonal
Chris@49:       
Chris@49:       for(uword i=0; i<N; ++i)
Chris@49:         {
Chris@49:         const eT* A_data   = A.colptr(i);
Chris@49:               eT* out_data = out.colptr(i);
Chris@49:         
Chris@49:         arrayops::copy( out_data, A_data, i+1 );
Chris@49:         }
Chris@49:       }
Chris@49:     else
Chris@49:       {
Chris@49:       // lower triangular: copy the diagonal and the elements below the diagonal
Chris@49:       
Chris@49:       for(uword i=0; i<N; ++i)
Chris@49:         {
Chris@49:         const eT* A_data   = A.colptr(i);
Chris@49:               eT* out_data = out.colptr(i);
Chris@49:         
Chris@49:         arrayops::copy( &out_data[i], &A_data[i], N-i );
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   
Chris@49:   if(upper)
Chris@49:     {
Chris@49:     // reflect elements across the diagonal from upper triangle to lower triangle
Chris@49:     
Chris@49:     for(uword col=1; col < N; ++col)
Chris@49:       {
Chris@49:       const eT* coldata = out.colptr(col);
Chris@49:       
Chris@49:       for(uword row=0; row < col; ++row)
Chris@49:         {
Chris@49:         out.at(col,row) = std::conj(coldata[row]);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // reflect elements across the diagonal from lower triangle to upper triangle
Chris@49:     
Chris@49:     for(uword col=0; col < N; ++col)
Chris@49:       {
Chris@49:       const eT* coldata = out.colptr(col);
Chris@49:       
Chris@49:       for(uword row=(col+1); row < N; ++row)
Chris@49:         {
Chris@49:         out.at(col,row) = std::conj(coldata[row]);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}