"Fossies" - the Fresh Open Source Software Archive

Member "armadillo-9.800.3/include/armadillo_bits/op_trimat_meat.hpp" (16 Jun 2016, 9932 Bytes) of package /linux/misc/armadillo-9.800.3.tar.xz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "op_trimat_meat.hpp" see the Fossies "Dox" file reference documentation.

    1 // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
    2 // Copyright 2008-2016 National ICT Australia (NICTA)
    3 // 
    4 // Licensed under the Apache License, Version 2.0 (the "License");
    5 // you may not use this file except in compliance with the License.
    6 // You may obtain a copy of the License at
    7 // http://www.apache.org/licenses/LICENSE-2.0
    8 // 
    9 // Unless required by applicable law or agreed to in writing, software
   10 // distributed under the License is distributed on an "AS IS" BASIS,
   11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   12 // See the License for the specific language governing permissions and
   13 // limitations under the License.
   14 // ------------------------------------------------------------------------
   15 
   16 
   17 //! \addtogroup op_trimat
   18 //! @{
   19 
   20 
   21 
   22 template<typename eT>
   23 inline
   24 void
   25 op_trimat::fill_zeros(Mat<eT>& out, const bool upper)
   26   {
   27   arma_extra_debug_sigprint();
   28   
   29   const uword N = out.n_rows;
   30   
   31   if(upper)
   32     {
   33     // upper triangular: set all elements below the diagonal to zero
   34     
   35     for(uword i=0; i<N; ++i)
   36       {
   37       eT* data = out.colptr(i);
   38       
   39       arrayops::inplace_set( &data[i+1], eT(0), (N-(i+1)) );
   40       }
   41     }
   42   else
   43     {
   44     // lower triangular: set all elements above the diagonal to zero
   45     
   46     for(uword i=1; i<N; ++i)
   47       {
   48       eT* data = out.colptr(i);
   49       
   50       arrayops::inplace_set( data, eT(0), i );
   51       }
   52     }
   53   }
   54 
   55 
   56 
   57 template<typename T1>
   58 inline
   59 void
   60 op_trimat::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trimat>& in)
   61   {
   62   arma_extra_debug_sigprint();
   63   
   64   typedef typename T1::elem_type eT;
   65   
   66   const unwrap<T1>   tmp(in.m);
   67   const Mat<eT>& A = tmp.M;
   68   
   69   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" );
   70   
   71   const uword N     = A.n_rows;
   72   const bool  upper = (in.aux_uword_a == 0);
   73   
   74   if(&out != &A)
   75     {
   76     out.copy_size(A);
   77     
   78     if(upper)
   79       {
   80       // upper triangular: copy the diagonal and the elements above the diagonal
   81       for(uword i=0; i<N; ++i)
   82         {
   83         const eT* A_data   = A.colptr(i);
   84               eT* out_data = out.colptr(i);
   85         
   86         arrayops::copy( out_data, A_data, i+1 );
   87         }
   88       }
   89     else
   90       {
   91       // lower triangular: copy the diagonal and the elements below the diagonal
   92       for(uword i=0; i<N; ++i)
   93         {
   94         const eT* A_data   = A.colptr(i);
   95               eT* out_data = out.colptr(i);
   96         
   97         arrayops::copy( &out_data[i], &A_data[i], N-i );
   98         }
   99       }
  100     }
  101   
  102   op_trimat::fill_zeros(out, upper);
  103   }
  104 
  105 
  106 
  107 template<typename T1>
  108 inline
  109 void
  110 op_trimat::apply(Mat<typename T1::elem_type>& out, const Op<Op<T1, op_htrans>, op_trimat>& in)
  111   {
  112   arma_extra_debug_sigprint();
  113   
  114   typedef typename T1::elem_type eT;
  115   
  116   const unwrap<T1>   tmp(in.m.m);
  117   const Mat<eT>& A = tmp.M;
  118   
  119   const bool upper = (in.aux_uword_a == 0);
  120   
  121   op_trimat::apply_htrans(out, A, upper);
  122   }
  123 
  124 
  125 
  126 template<typename eT>
  127 inline
  128 void
  129 op_trimat::apply_htrans
  130   (
  131         Mat<eT>& out,
  132   const Mat<eT>& A,
  133   const bool     upper,
  134   const typename arma_not_cx<eT>::result* junk
  135   )
  136   {
  137   arma_extra_debug_sigprint();
  138   arma_ignore(junk);
  139   
  140   // This specialisation is for trimatl(trans(X)) = trans(trimatu(X)) and also
  141   // trimatu(trans(X)) = trans(trimatl(X)).  We want to avoid the creation of an
  142   // extra temporary.
  143   
  144   // It doesn't matter if the input and output matrices are the same; we will
  145   // pull data from the upper or lower triangular to the lower or upper
  146   // triangular (respectively) and then set the rest to 0, so overwriting issues
  147   // aren't present.
  148   
  149   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" );
  150   
  151   const uword N = A.n_rows;
  152   
  153   if(&out != &A)
  154     {
  155     out.copy_size(A);
  156     }
  157   
  158   // We can't really get away with any array copy operations here,
  159   // unfortunately...
  160   
  161   if(upper)
  162     {
  163     // Upper triangular: but since we're transposing, we're taking the lower
  164     // triangular and putting it in the upper half.
  165     for(uword row = 0; row < N; ++row)
  166       {
  167       eT* out_colptr = out.colptr(row);
  168       
  169       for(uword col = 0; col <= row; ++col)
  170         {
  171         //out.at(col, row) = A.at(row, col);
  172         out_colptr[col] = A.at(row, col);
  173         }
  174       }
  175     }
  176   else
  177     {
  178     // Lower triangular: but since we're transposing, we're taking the upper
  179     // triangular and putting it in the lower half.
  180     for(uword row = 0; row < N; ++row)
  181       {
  182       for(uword col = row; col < N; ++col)
  183         {
  184         out.at(col, row) = A.at(row, col);
  185         }
  186       }
  187     }
  188   
  189   op_trimat::fill_zeros(out, upper);
  190   }
  191 
  192 
  193 
  194 template<typename eT>
  195 inline
  196 void
  197 op_trimat::apply_htrans
  198   (
  199         Mat<eT>& out,
  200   const Mat<eT>& A,
  201   const bool     upper,
  202   const typename arma_cx_only<eT>::result* junk
  203   )
  204   {
  205   arma_extra_debug_sigprint();
  206   arma_ignore(junk);
  207   
  208   arma_debug_check( (A.is_square() == false), "trimatu()/trimatl(): given matrix must be square sized" );
  209   
  210   const uword N = A.n_rows;
  211   
  212   if(&out != &A)
  213     {
  214     out.copy_size(A);
  215     }
  216   
  217   if(upper)
  218     {
  219     // Upper triangular: but since we're transposing, we're taking the lower
  220     // triangular and putting it in the upper half.
  221     for(uword row = 0; row < N; ++row)
  222       {
  223       eT* out_colptr = out.colptr(row);
  224       
  225       for(uword col = 0; col <= row; ++col)
  226         {
  227         //out.at(col, row) = std::conj( A.at(row, col) );
  228         out_colptr[col] = std::conj( A.at(row, col) );
  229         }
  230       }
  231     }
  232   else
  233     {
  234     // Lower triangular: but since we're transposing, we're taking the upper
  235     // triangular and putting it in the lower half.
  236     for(uword row = 0; row < N; ++row)
  237       {
  238       for(uword col = row; col < N; ++col)
  239         {
  240         out.at(col, row) = std::conj( A.at(row, col) );
  241         }
  242       }
  243     }
  244   
  245   op_trimat::fill_zeros(out, upper);
  246   }
  247 
  248 
  249 
  250 //
  251 
  252 
  253 
  254 template<typename T1>
  255 inline
  256 void
  257 op_trimatu_ext::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trimatu_ext>& in)
  258   {
  259   arma_extra_debug_sigprint();
  260   
  261   typedef typename T1::elem_type eT;
  262   
  263   const unwrap<T1>   tmp(in.m);
  264   const Mat<eT>& A = tmp.M;
  265   
  266   arma_debug_check( (A.is_square() == false), "trimatu(): given matrix must be square sized" );
  267   
  268   const uword row_offset = in.aux_uword_a;
  269   const uword col_offset = in.aux_uword_b;
  270   
  271   const uword n_rows = A.n_rows;
  272   const uword n_cols = A.n_cols;
  273   
  274   arma_debug_check( ((row_offset > 0) && (row_offset >= n_rows)) || ((col_offset > 0) && (col_offset >= n_cols)), "trimatu(): requested diagonal is out of bounds" );
  275   
  276   if(&out != &A)
  277     {
  278     out.copy_size(A);
  279     
  280     const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
  281     
  282     for(uword i=0; i < n_cols; ++i)
  283       {
  284       const uword col = i + col_offset;
  285         
  286       if(i < N)
  287         {
  288         const uword end_row = i + row_offset;
  289         
  290         for(uword row=0; row <= end_row; ++row)
  291           {
  292           out.at(row,col) = A.at(row,col);
  293           }
  294         }
  295       else
  296         {
  297         if(col < n_cols)
  298           {
  299           arrayops::copy(out.colptr(col), A.colptr(col), n_rows);
  300           }
  301         }
  302       }
  303     }
  304   
  305   op_trimatu_ext::fill_zeros(out, row_offset, col_offset);
  306   }
  307 
  308 
  309 
  310 template<typename eT>
  311 inline
  312 void
  313 op_trimatu_ext::fill_zeros(Mat<eT>& out, const uword row_offset, const uword col_offset)
  314   {
  315   arma_extra_debug_sigprint();
  316   
  317   const uword n_rows = out.n_rows;
  318   const uword n_cols = out.n_cols;
  319   
  320   const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
  321   
  322   for(uword col=0; col < col_offset; ++col)
  323     {
  324     arrayops::fill_zeros(out.colptr(col), n_rows);
  325     }
  326   
  327   for(uword i=0; i < N; ++i)
  328     {
  329     const uword start_row = i + row_offset + 1;
  330     const uword col       = i + col_offset;
  331     
  332     for(uword row=start_row; row < n_rows; ++row)
  333       {
  334       out.at(row,col) = eT(0);
  335       }
  336     }
  337   }
  338 
  339 
  340 
  341 //
  342 
  343 
  344 
  345 template<typename T1>
  346 inline
  347 void
  348 op_trimatl_ext::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_trimatl_ext>& in)
  349   {
  350   arma_extra_debug_sigprint();
  351   
  352   typedef typename T1::elem_type eT;
  353   
  354   const unwrap<T1>   tmp(in.m);
  355   const Mat<eT>& A = tmp.M;
  356   
  357   arma_debug_check( (A.is_square() == false), "trimatl(): given matrix must be square sized" );
  358   
  359   const uword row_offset = in.aux_uword_a;
  360   const uword col_offset = in.aux_uword_b;
  361   
  362   const uword n_rows = A.n_rows;
  363   const uword n_cols = A.n_cols;
  364   
  365   arma_debug_check( ((row_offset > 0) && (row_offset >= n_rows)) || ((col_offset > 0) && (col_offset >= n_cols)), "trimatl(): requested diagonal is out of bounds" );
  366   
  367   if(&out != &A)
  368     {
  369     out.copy_size(A);
  370     
  371     const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
  372     
  373     for(uword col=0; col < col_offset; ++col)
  374       {
  375       arrayops::copy( out.colptr(col), A.colptr(col), n_rows );
  376       }
  377     
  378     for(uword i=0; i<N; ++i)
  379       {
  380       const uword start_row = i + row_offset;
  381       const uword       col = i + col_offset;
  382       
  383       for(uword row=start_row; row < n_rows; ++row)
  384         {
  385         out.at(row,col) = A.at(row,col);
  386         }
  387       }
  388     }
  389   
  390   op_trimatl_ext::fill_zeros(out, row_offset, col_offset);
  391   }
  392 
  393 
  394 
  395 template<typename eT>
  396 inline
  397 void
  398 op_trimatl_ext::fill_zeros(Mat<eT>& out, const uword row_offset, const uword col_offset)
  399   {
  400   arma_extra_debug_sigprint();
  401   
  402   const uword n_rows = out.n_rows;
  403   const uword n_cols = out.n_cols;
  404   
  405   const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
  406   
  407   for(uword i=0; i < n_cols; ++i)
  408     {
  409     const uword col = i + col_offset;
  410       
  411     if(i < N)
  412       {
  413       const uword end_row = i + row_offset;
  414       
  415       for(uword row=0; row < end_row; ++row)
  416         {
  417         out.at(row,col) = eT(0);
  418         }
  419       }
  420     else
  421       {
  422       if(col < n_cols)
  423         {
  424         arrayops::fill_zeros(out.colptr(col), n_rows);
  425         }
  426       }
  427     }
  428   }
  429 
  430 
  431 
  432 //! @}