Move into version control

1 files changed, 456 insertions, 0 deletions

diff --git a/src/armadillo/include/armadillo_bits/spop_diagmat_meat.hpp b/src/armadillo/include/armadillo_bits/spop_diagmat_meat.hpp
new file mode 100644
index 0000000..a6f9faf
--- /dev/null
+++ b/src/armadillo/include/armadillo_bits/spop_diagmat_meat.hpp
@@ -0,0 +1,456 @@
+// SPDX-License-Identifier: Apache-2.0
+// 
+// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
+// Copyright 2008-2016 National ICT Australia (NICTA)
+// 
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+// http://www.apache.org/licenses/LICENSE-2.0
+// 
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+// ------------------------------------------------------------------------
+
+
+//! \addtogroup spop_diagmat
+//! @{
+
+
+
+template<typename T1>
+inline
+void
+spop_diagmat::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_diagmat>& in)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  if(in.is_alias(out) == false)
+    {
+    spop_diagmat::apply_noalias(out, in.m);
+    }
+  else
+    {
+    SpMat<eT> tmp;
+    
+    spop_diagmat::apply_noalias(tmp, in.m);
+    
+    out.steal_mem(tmp);
+    }
+  }
+
+
+
+template<typename T1>
+inline
+void
+spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpBase<typename T1::elem_type, T1>& expr)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const SpProxy<T1> P(expr.get_ref());
+  
+  const uword P_n_rows = P.get_n_rows();
+  const uword P_n_cols = P.get_n_cols();
+  const uword P_n_nz   = P.get_n_nonzero();
+  
+  const bool P_is_vec = (P_n_rows == 1) || (P_n_cols == 1);
+  
+  if(P_is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const uword N = (P_n_rows == 1) ? P_n_cols : P_n_rows;
+    
+    out.zeros(N, N);
+    
+    if(P_n_nz == 0)  { return; }
+    
+    typename SpProxy<T1>::const_iterator_type it = P.begin();
+    
+    if(P_n_cols == 1)
+      {
+      for(uword i=0; i < P_n_nz; ++i)
+        {
+        const uword row = it.row();
+        
+        out.at(row,row) = (*it);
+        
+        ++it;
+        }
+      }
+    else
+    if(P_n_rows == 1)
+      {
+      for(uword i=0; i < P_n_nz; ++i)
+        {
+        const uword col = it.col();
+        
+        out.at(col,col) = (*it);
+        
+        ++it;
+        }
+      }
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    out.zeros(P_n_rows, P_n_cols);
+    
+    const uword N = (std::min)(P_n_rows, P_n_cols);
+    
+    if( (is_SpMat<typename SpProxy<T1>::stored_type>::value) && (P_n_nz >= 5*N) )
+      {
+      const unwrap_spmat<typename SpProxy<T1>::stored_type> U(P.Q);
+      
+      const SpMat<eT>& X = U.M;
+      
+      for(uword i=0; i < N; ++i)
+        {
+        const eT val = X.at(i,i);  // use binary search
+        
+        if(val != eT(0))  { out.at(i,i) = val; }
+        }
+      }
+    else
+      {
+      if(P_n_nz == 0)  { return; }
+      
+      typename SpProxy<T1>::const_iterator_type it = P.begin();
+      
+      for(uword i=0; i < P_n_nz; ++i)
+        {
+        const uword row = it.row();
+        const uword col = it.col();
+        
+        if(row == col)  { out.at(row,row) = (*it); }
+        
+        ++it;
+        }
+      }
+    }
+  }
+
+
+
+template<typename T1, typename T2>
+inline
+void
+spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_plus>& expr)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const unwrap_spmat<T1> UA(expr.A);
+  const unwrap_spmat<T2> UB(expr.B);
+  
+  const SpMat<eT>& A = UA.M;
+  const SpMat<eT>& B = UB.M;
+  
+  arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "addition");
+  
+  const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
+  
+  if(is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
+    
+    out.zeros(N,N);
+    
+    if(A.n_rows == 1)
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) + B.at(0,i); }
+      }
+    else
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) + B.at(i,0); }
+      }
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    SpMat<eT> AA;  spop_diagmat::apply_noalias(AA, A);
+    SpMat<eT> BB;  spop_diagmat::apply_noalias(BB, B);
+    
+    out = AA + BB;
+    }
+  }
+
+
+
+template<typename T1, typename T2>
+inline
+void
+spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_minus>& expr)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const unwrap_spmat<T1> UA(expr.A);
+  const unwrap_spmat<T2> UB(expr.B);
+  
+  const SpMat<eT>& A = UA.M;
+  const SpMat<eT>& B = UB.M;
+  
+  arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "subtraction");
+  
+  const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
+  
+  if(is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
+    
+    out.zeros(N,N);
+    
+    if(A.n_rows == 1)
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) - B.at(0,i); }
+      }
+    else
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) - B.at(i,0); }
+      }
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    SpMat<eT> AA;  spop_diagmat::apply_noalias(AA, A);
+    SpMat<eT> BB;  spop_diagmat::apply_noalias(BB, B);
+    
+    out = AA - BB;
+    }
+  }
+
+
+
+template<typename T1, typename T2>
+inline
+void
+spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_schur>& expr)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const unwrap_spmat<T1> UA(expr.A);
+  const unwrap_spmat<T2> UB(expr.B);
+  
+  const SpMat<eT>& A = UA.M;
+  const SpMat<eT>& B = UB.M;
+  
+  arma_debug_assert_same_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "element-wise multiplication");
+  
+  const bool is_vec = (A.n_rows == 1) || (A.n_cols == 1);
+  
+  if(is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const uword N = (A.n_rows == 1) ? A.n_cols : A.n_rows;
+    
+    out.zeros(N,N);
+    
+    if(A.n_rows == 1)
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(0,i) * B.at(0,i); }
+      }
+    else
+      {
+      for(uword i=0; i < N; ++i) { out.at(i,i) = A.at(i,0) * B.at(i,0); }
+      }
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    SpMat<eT> AA;  spop_diagmat::apply_noalias(AA, A);
+    SpMat<eT> BB;  spop_diagmat::apply_noalias(BB, B);
+    
+    out = AA % BB;
+    }
+  }
+
+
+
+template<typename T1, typename T2>
+inline
+void
+spop_diagmat::apply_noalias(SpMat<typename T1::elem_type>& out, const SpGlue<T1,T2,spglue_times>& expr)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const unwrap_spmat<T1> UA(expr.A);
+  const unwrap_spmat<T2> UB(expr.B);
+  
+  const SpMat<eT>& A = UA.M;
+  const SpMat<eT>& B = UB.M;
+  
+  arma_debug_assert_mul_size(A.n_rows, A.n_cols, B.n_rows, B.n_cols, "matrix multiplication");
+  
+  const uword C_n_rows = A.n_rows;
+  const uword C_n_cols = B.n_cols;
+  
+  const bool is_vec = (C_n_rows == 1) || (C_n_cols == 1);
+  
+  if(is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const SpMat<eT> C = A*B;
+    
+    spop_diagmat::apply_noalias(out, C);
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    const uword N = (std::min)(C_n_rows, C_n_cols);
+    
+    if( (A.n_nonzero >= 5*N) || (B.n_nonzero >= 5*N) )
+      {
+      out.zeros(C_n_rows, C_n_cols);
+      
+      for(uword k=0; k < N; ++k)
+        {
+        typename SpMat<eT>::const_col_iterator B_it     = B.begin_col_no_sync(k);
+        typename SpMat<eT>::const_col_iterator B_it_end = B.end_col_no_sync(k);
+        
+        eT acc = eT(0);
+        
+        while(B_it != B_it_end)
+          {
+          const eT    B_val = (*B_it);
+          const uword i     = B_it.row();
+          
+          acc += A.at(k,i) * B_val;
+          
+          ++B_it;
+          }
+        
+        out(k,k) = acc;
+        }
+      }
+    else
+      {
+      const SpMat<eT> C = A*B;
+      
+      spop_diagmat::apply_noalias(out, C);
+      }
+    }
+  }
+
+
+
+//
+//
+
+
+
+template<typename T1>
+inline
+void
+spop_diagmat2::apply(SpMat<typename T1::elem_type>& out, const SpOp<T1, spop_diagmat2>& in)
+  {
+  arma_extra_debug_sigprint();
+  
+  typedef typename T1::elem_type eT;
+  
+  const uword row_offset = in.aux_uword_a;
+  const uword col_offset = in.aux_uword_b;
+  
+  const unwrap_spmat<T1> U(in.m);
+  
+  if(U.is_alias(out))
+    {
+    SpMat<eT> tmp;
+    
+    spop_diagmat2::apply_noalias(tmp, U.M, row_offset, col_offset);
+    
+    out.steal_mem(tmp);
+    }
+  else
+    {
+    spop_diagmat2::apply_noalias(out, U.M, row_offset, col_offset);
+    }
+  }
+
+
+
+template<typename eT>
+inline
+void
+spop_diagmat2::apply_noalias(SpMat<eT>& out, const SpMat<eT>& X, const uword row_offset, const uword col_offset)
+  {
+  arma_extra_debug_sigprint();
+  
+  const uword n_rows = X.n_rows;
+  const uword n_cols = X.n_cols;
+  const uword n_elem = X.n_elem;
+  
+  if(n_elem == 0)  { out.reset(); return; }
+  
+  const bool X_is_vec = (n_rows == 1) || (n_cols == 1);
+  
+  if(X_is_vec)    // generate a diagonal matrix out of a vector
+    {
+    const uword n_pad = (std::max)(row_offset, col_offset);
+    
+    out.zeros(n_elem + n_pad, n_elem + n_pad);
+    
+    const uword X_n_nz = X.n_nonzero;
+    
+    if(X_n_nz == 0)  { return; }
+    
+    typename SpMat<eT>::const_iterator it = X.begin();
+      
+    if(n_cols == 1)
+      {
+      for(uword i=0; i < X_n_nz; ++i)
+        {
+        const uword row = it.row();
+        
+        out.at(row_offset + row, col_offset + row) = (*it);
+        
+        ++it;
+        }
+      }
+    else
+    if(n_rows == 1)
+      {
+      for(uword i=0; i < X_n_nz; ++i)
+        {
+        const uword col = it.col();
+        
+        out.at(row_offset + col, col_offset + col) = (*it);
+        
+        ++it;
+        }
+      }
+    }
+  else   // generate a diagonal matrix out of a matrix
+    {
+    arma_debug_check_bounds
+      (
+      ((row_offset > 0) && (row_offset >= n_rows)) || ((col_offset > 0) && (col_offset >= n_cols)),
+      "diagmat(): requested diagonal out of bounds"
+      );
+    
+    out.zeros(n_rows, n_cols);
+    
+    if(X.n_nonzero == 0)  { return; }
+    
+    const uword N = (std::min)(n_rows - row_offset, n_cols - col_offset);
+    
+    for(uword i=0; i<N; ++i)
+      {
+      const uword row = i + row_offset;
+      const uword col = i + col_offset;
+      
+      const eT val = X.at(row,col);
+      
+      if(val != eT(0))  { out.at(row,col) = val; }
+      }
+    }
+  }
+
+
+
+//! @}