max@0: // Copyright (C) 2010-2011 NICTA (www.nicta.com.au)
max@0: // Copyright (C) 2010-2011 Conrad Sanderson
max@0: // Copyright (C) 2011      Ryan Curtin
max@0: // 
max@0: // This file is part of the Armadillo C++ library.
max@0: // It is provided without any warranty of fitness
max@0: // for any purpose. You can redistribute this file
max@0: // and/or modify it under the terms of the GNU
max@0: // Lesser General Public License (LGPL) as published
max@0: // by the Free Software Foundation, either version 3
max@0: // of the License or (at your option) any later version.
max@0: // (see http://www.opensource.org/licenses for more info)
max@0: 
max@0: 
max@0: //! \addtogroup op_trimat
max@0: //! @{
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: void
max@0: op_trimat::fill_zeros(Mat<eT>& out, const bool upper)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   const uword N = out.n_rows;
max@0:   
max@0:   if(upper)
max@0:     {
max@0:     // upper triangular: set all elements below the diagonal to zero
max@0:     
max@0:     for(uword i=0; i<N; ++i)
max@0:       {
max@0:       eT* data = out.colptr(i);
max@0:       
max@0:       arrayops::inplace_set( &data[i+1], eT(0), (N-(i+1)) );
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     // lower triangular: set all elements above the diagonal to zero
max@0:     
max@0:     for(uword i=1; i<N; ++i)
max@0:       {
max@0:       eT* data = out.colptr(i);
max@0:       
max@0:       arrayops::inplace_set( data, eT(0), i );
max@0:       }
max@0:     }
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: void
max@0: op_trimat::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trimat>& in)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   
max@0:   const unwrap<T1> tmp(in.m);
max@0:   const Mat<eT>& A = tmp.M;
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square" );
max@0:   
max@0:   const uword  N     = A.n_rows;
max@0:   const bool upper = (in.aux_uword_a == 0);
max@0:   
max@0:   if(&out != &A)
max@0:     {
max@0:     out.copy_size(A);
max@0:     
max@0:     if(upper)
max@0:       {
max@0:       // upper triangular: copy the diagonal and the elements above the diagonal
max@0:       for(uword i=0; i<N; ++i)
max@0:         {
max@0:         const eT* A_data   = A.colptr(i);
max@0:               eT* out_data = out.colptr(i);
max@0:         
max@0:         arrayops::copy( out_data, A_data, i+1 );
max@0:         }
max@0:       }
max@0:     else
max@0:       {
max@0:       // lower triangular: copy the diagonal and the elements below the diagonal
max@0:       for(uword i=0; i<N; ++i)
max@0:         {
max@0:         const eT* A_data   = A.colptr(i);
max@0:               eT* out_data = out.colptr(i);
max@0:         
max@0:         arrayops::copy( &out_data[i], &A_data[i], N-i );
max@0:         }
max@0:       }
max@0:     }
max@0:   
max@0:   op_trimat::fill_zeros(out, upper);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename T1>
max@0: inline
max@0: void
max@0: op_trimat::apply(Mat<typename T1::elem_type>& out, const Op<Op<T1, op_htrans>, op_trimat>& in)
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   
max@0:   typedef typename T1::elem_type eT;
max@0:   
max@0:   const unwrap<T1>   tmp(in.m.m);
max@0:   const Mat<eT>& A = tmp.M;
max@0:   
max@0:   const bool upper = (in.aux_uword_a == 0);
max@0:   
max@0:   op_trimat::apply_htrans(out, A, upper);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: void
max@0: op_trimat::apply_htrans
max@0:   (
max@0:         Mat<eT>& out,
max@0:   const Mat<eT>& A,
max@0:   const bool     upper,
max@0:   const typename arma_not_cx<eT>::result* junk
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   // This specialisation is for trimatl(trans(X)) = trans(trimatu(X)) and also
max@0:   // trimatu(trans(X)) = trans(trimatl(X)).  We want to avoid the creation of an
max@0:   // extra temporary.
max@0:   
max@0:   // It doesn't matter if the input and output matrices are the same; we will
max@0:   // pull data from the upper or lower triangular to the lower or upper
max@0:   // triangular (respectively) and then set the rest to 0, so overwriting issues
max@0:   // aren't present.
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square" );
max@0:   
max@0:   const uword N = A.n_rows;
max@0:   
max@0:   if(&out != &A)
max@0:     {
max@0:     out.copy_size(A);
max@0:     }
max@0:   
max@0:   // We can't really get away with any array copy operations here,
max@0:   // unfortunately...
max@0:   
max@0:   if(upper)
max@0:     {
max@0:     // Upper triangular: but since we're transposing, we're taking the lower
max@0:     // triangular and putting it in the upper half.
max@0:     for(uword row = 0; row < N; ++row)
max@0:       {
max@0:       eT* out_colptr = out.colptr(row);
max@0:       
max@0:       for(uword col = 0; col <= row; ++col)
max@0:         {
max@0:         //out.at(col, row) = A.at(row, col);
max@0:         out_colptr[col] = A.at(row, col);
max@0:         }
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     // Lower triangular: but since we're transposing, we're taking the upper
max@0:     // triangular and putting it in the lower half.
max@0:     for(uword row = 0; row < N; ++row)
max@0:       {
max@0:       for(uword col = row; col < N; ++col)
max@0:         {
max@0:         out.at(col, row) = A.at(row, col);
max@0:         }
max@0:       }
max@0:     }
max@0:   
max@0:   op_trimat::fill_zeros(out, upper);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: template<typename eT>
max@0: inline
max@0: void
max@0: op_trimat::apply_htrans
max@0:   (
max@0:         Mat<eT>& out,
max@0:   const Mat<eT>& A,
max@0:   const bool     upper,
max@0:   const typename arma_cx_only<eT>::result* junk
max@0:   )
max@0:   {
max@0:   arma_extra_debug_sigprint();
max@0:   arma_ignore(junk);
max@0:   
max@0:   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square" );
max@0:   
max@0:   const uword N = A.n_rows;
max@0:   
max@0:   if(&out != &A)
max@0:     {
max@0:     out.copy_size(A);
max@0:     }
max@0:   
max@0:   if(upper)
max@0:     {
max@0:     // Upper triangular: but since we're transposing, we're taking the lower
max@0:     // triangular and putting it in the upper half.
max@0:     for(uword row = 0; row < N; ++row)
max@0:       {
max@0:       eT* out_colptr = out.colptr(row);
max@0:       
max@0:       for(uword col = 0; col <= row; ++col)
max@0:         {
max@0:         //out.at(col, row) = std::conj( A.at(row, col) );
max@0:         out_colptr[col] = std::conj( A.at(row, col) );
max@0:         }
max@0:       }
max@0:     }
max@0:   else
max@0:     {
max@0:     // Lower triangular: but since we're transposing, we're taking the upper
max@0:     // triangular and putting it in the lower half.
max@0:     for(uword row = 0; row < N; ++row)
max@0:       {
max@0:       for(uword col = row; col < N; ++col)
max@0:         {
max@0:         out.at(col, row) = std::conj( A.at(row, col) );
max@0:         }
max@0:       }
max@0:     }
max@0:   
max@0:   op_trimat::fill_zeros(out, upper);
max@0:   }
max@0: 
max@0: 
max@0: 
max@0: //! @}