Chris@49: // Copyright (C) 2010-2012 NICTA (www.nicta.com.au)
Chris@49: // Copyright (C) 2010-2012 Conrad Sanderson
Chris@49: // Copyright (C) 2011      Ryan Curtin
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_trimat
Chris@49: //! @{
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: op_trimat::fill_zeros(Mat<eT>& out, const bool upper)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   const uword N = out.n_rows;
Chris@49:   
Chris@49:   if(upper)
Chris@49:     {
Chris@49:     // upper triangular: set all elements below the diagonal to zero
Chris@49:     
Chris@49:     for(uword i=0; i<N; ++i)
Chris@49:       {
Chris@49:       eT* data = out.colptr(i);
Chris@49:       
Chris@49:       arrayops::inplace_set( &data[i+1], eT(0), (N-(i+1)) );
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // lower triangular: set all elements above the diagonal to zero
Chris@49:     
Chris@49:     for(uword i=1; i<N; ++i)
Chris@49:       {
Chris@49:       eT* data = out.colptr(i);
Chris@49:       
Chris@49:       arrayops::inplace_set( data, eT(0), i );
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_trimat::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trimat>& in)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
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), "trimatu()/trimatl(): 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:       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:       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:   op_trimat::fill_zeros(out, upper);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename T1>
Chris@49: inline
Chris@49: void
Chris@49: op_trimat::apply(Mat<typename T1::elem_type>& out, const Op<Op<T1, op_htrans>, op_trimat>& in)
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   
Chris@49:   typedef typename T1::elem_type eT;
Chris@49:   
Chris@49:   const unwrap<T1>   tmp(in.m.m);
Chris@49:   const Mat<eT>& A = tmp.M;
Chris@49:   
Chris@49:   const bool upper = (in.aux_uword_a == 0);
Chris@49:   
Chris@49:   op_trimat::apply_htrans(out, A, upper);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: op_trimat::apply_htrans
Chris@49:   (
Chris@49:         Mat<eT>& out,
Chris@49:   const Mat<eT>& A,
Chris@49:   const bool     upper,
Chris@49:   const typename arma_not_cx<eT>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   // This specialisation is for trimatl(trans(X)) = trans(trimatu(X)) and also
Chris@49:   // trimatu(trans(X)) = trans(trimatl(X)).  We want to avoid the creation of an
Chris@49:   // extra temporary.
Chris@49:   
Chris@49:   // It doesn't matter if the input and output matrices are the same; we will
Chris@49:   // pull data from the upper or lower triangular to the lower or upper
Chris@49:   // triangular (respectively) and then set the rest to 0, so overwriting issues
Chris@49:   // aren't present.
Chris@49:   
Chris@49:   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square" );
Chris@49:   
Chris@49:   const uword N = A.n_rows;
Chris@49:   
Chris@49:   if(&out != &A)
Chris@49:     {
Chris@49:     out.copy_size(A);
Chris@49:     }
Chris@49:   
Chris@49:   // We can't really get away with any array copy operations here,
Chris@49:   // unfortunately...
Chris@49:   
Chris@49:   if(upper)
Chris@49:     {
Chris@49:     // Upper triangular: but since we're transposing, we're taking the lower
Chris@49:     // triangular and putting it in the upper half.
Chris@49:     for(uword row = 0; row < N; ++row)
Chris@49:       {
Chris@49:       eT* out_colptr = out.colptr(row);
Chris@49:       
Chris@49:       for(uword col = 0; col <= row; ++col)
Chris@49:         {
Chris@49:         //out.at(col, row) = A.at(row, col);
Chris@49:         out_colptr[col] = A.at(row, col);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // Lower triangular: but since we're transposing, we're taking the upper
Chris@49:     // triangular and putting it in the lower half.
Chris@49:     for(uword row = 0; row < N; ++row)
Chris@49:       {
Chris@49:       for(uword col = row; col < N; ++col)
Chris@49:         {
Chris@49:         out.at(col, row) = A.at(row, col);
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   op_trimat::fill_zeros(out, upper);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: template<typename eT>
Chris@49: inline
Chris@49: void
Chris@49: op_trimat::apply_htrans
Chris@49:   (
Chris@49:         Mat<eT>& out,
Chris@49:   const Mat<eT>& A,
Chris@49:   const bool     upper,
Chris@49:   const typename arma_cx_only<eT>::result* junk
Chris@49:   )
Chris@49:   {
Chris@49:   arma_extra_debug_sigprint();
Chris@49:   arma_ignore(junk);
Chris@49:   
Chris@49:   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square" );
Chris@49:   
Chris@49:   const uword N = A.n_rows;
Chris@49:   
Chris@49:   if(&out != &A)
Chris@49:     {
Chris@49:     out.copy_size(A);
Chris@49:     }
Chris@49:   
Chris@49:   if(upper)
Chris@49:     {
Chris@49:     // Upper triangular: but since we're transposing, we're taking the lower
Chris@49:     // triangular and putting it in the upper half.
Chris@49:     for(uword row = 0; row < N; ++row)
Chris@49:       {
Chris@49:       eT* out_colptr = out.colptr(row);
Chris@49:       
Chris@49:       for(uword col = 0; col <= row; ++col)
Chris@49:         {
Chris@49:         //out.at(col, row) = std::conj( A.at(row, col) );
Chris@49:         out_colptr[col] = std::conj( A.at(row, col) );
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   else
Chris@49:     {
Chris@49:     // Lower triangular: but since we're transposing, we're taking the upper
Chris@49:     // triangular and putting it in the lower half.
Chris@49:     for(uword row = 0; row < N; ++row)
Chris@49:       {
Chris@49:       for(uword col = row; col < N; ++col)
Chris@49:         {
Chris@49:         out.at(col, row) = std::conj( A.at(row, col) );
Chris@49:         }
Chris@49:       }
Chris@49:     }
Chris@49:   
Chris@49:   op_trimat::fill_zeros(out, upper);
Chris@49:   }
Chris@49: 
Chris@49: 
Chris@49: 
Chris@49: //! @}