op_powmat_meat.hpp 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223
  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. //! \addtogroup op_powmat
  16. //! @{
  17. template<typename T1>
  18. inline
  19. void
  20. op_powmat::apply(Mat<typename T1::elem_type>& out, const Op<T1, op_powmat>& expr)
  21. {
  22. arma_extra_debug_sigprint();
  23. typedef typename T1::elem_type eT;
  24. const uword y = expr.aux_uword_a;
  25. const bool y_neg = (expr.aux_uword_b == uword(1));
  26. if(y_neg)
  27. {
  28. if(y == uword(1))
  29. {
  30. const bool inv_status = inv(out, expr.m);
  31. if(inv_status == false)
  32. {
  33. out.soft_reset();
  34. arma_stop_runtime_error("powmat(): matrix inverse failed");
  35. return;
  36. }
  37. }
  38. else
  39. {
  40. Mat<eT> X_inv;
  41. const bool inv_status = inv(X_inv, expr.m);
  42. if(inv_status == false)
  43. {
  44. out.soft_reset();
  45. arma_stop_runtime_error("powmat(): matrix inverse failed");
  46. return;
  47. }
  48. op_powmat::apply(out, X_inv, y);
  49. }
  50. }
  51. else
  52. {
  53. const quasi_unwrap<T1> U(expr.m);
  54. arma_debug_check( (U.M.is_square() == false), "powmat(): given matrix must be square sized" );
  55. op_powmat::apply(out, U.M, y);
  56. }
  57. }
  58. template<typename eT>
  59. inline
  60. void
  61. op_powmat::apply(Mat<eT>& out, const Mat<eT>& X, const uword y)
  62. {
  63. arma_extra_debug_sigprint();
  64. const uword N = X.n_rows;
  65. if(y == uword(0)) { out.eye(N,N); return; }
  66. if(y == uword(1)) { out = X; return; }
  67. if(X.is_diagmat())
  68. {
  69. podarray<eT> tmp(N); // use temporary array in case we have aliasing
  70. for(uword i=0; i<N; ++i) { tmp[i] = eop_aux::pow(X.at(i,i), int(y)); }
  71. out.zeros(N,N);
  72. for(uword i=0; i<N; ++i) { out.at(i,i) = tmp[i]; }
  73. }
  74. else
  75. {
  76. if(y == uword(2)) { out = X*X; }
  77. else if(y == uword(3)) { const Mat<eT> tmp = X*X; out = X*tmp; }
  78. else if(y == uword(4)) { const Mat<eT> tmp = X*X; out = tmp*tmp; }
  79. else if(y == uword(5)) { const Mat<eT> tmp = X*X; out = X*tmp*tmp; }
  80. else
  81. {
  82. Mat<eT> tmp = X;
  83. out = X;
  84. uword z = y-1;
  85. while(z > 0)
  86. {
  87. if(z & 1) { out = tmp * out; }
  88. z /= uword(2);
  89. if(z > 0) { tmp = tmp * tmp; }
  90. }
  91. }
  92. }
  93. }
  94. template<typename T1>
  95. inline
  96. void
  97. op_powmat_cx::apply(Mat< std::complex<typename T1::pod_type> >& out, const mtOp<std::complex<typename T1::pod_type>,T1,op_powmat_cx>& expr)
  98. {
  99. arma_extra_debug_sigprint();
  100. typedef typename T1::elem_type in_eT;
  101. typedef typename T1::pod_type in_T;
  102. typedef std::complex<in_T> out_eT;
  103. const in_T y = std::real(expr.aux_out_eT);
  104. if( y == in_T(int(y)) )
  105. {
  106. arma_extra_debug_print("op_powmat_cx::apply(): integer exponent detected; redirecting to op_powmat");
  107. out = conv_to< Mat<out_eT> >::from( powmat(expr.m,int(y)) );
  108. return;
  109. }
  110. const quasi_unwrap<T1> U(expr.m);
  111. const Mat<in_eT>& A = U.M;
  112. arma_debug_check( (A.is_square() == false), "powmat(): given matrix must be square sized" );
  113. const uword N = A.n_rows;
  114. if(A.is_diagmat())
  115. {
  116. podarray<out_eT> tmp(N); // use temporary array in case we have aliasing
  117. for(uword i=0; i<N; ++i) { tmp[i] = eop_aux::pow( std::complex<in_T>(A.at(i,i)), y) ; }
  118. out.zeros(N,N);
  119. for(uword i=0; i<N; ++i) { out.at(i,i) = tmp[i]; }
  120. return;
  121. }
  122. #if defined(ARMA_OPTIMISE_SYMPD)
  123. const bool try_sympd = sympd_helper::guess_sympd_anysize(A);
  124. #else
  125. const bool try_sympd = false;
  126. #endif
  127. if(try_sympd)
  128. {
  129. Col<in_T> eigval;
  130. Mat<in_eT> eigvec;
  131. const bool eig_status = eig_sym(eigval, eigvec, A);
  132. if(eig_status)
  133. {
  134. eigval = pow(eigval, y);
  135. const Mat<in_eT> tmp = diagmat(eigval) * eigvec.t();
  136. out = conv_to< Mat<out_eT> >::from(eigvec * tmp);
  137. return;
  138. }
  139. // fallthrough
  140. }
  141. bool powmat_status = false;
  142. Col<out_eT> eigval;
  143. Mat<out_eT> eigvec;
  144. const bool eig_status = eig_gen(eigval, eigvec, A);
  145. if(eig_status)
  146. {
  147. eigval = pow(eigval, y);
  148. Mat<out_eT> eigvec_t = trans(eigvec);
  149. Mat<out_eT> tmp = diagmat(conj(eigval)) * eigvec_t;
  150. const bool solve_status = auxlib::solve_square_fast(out, eigvec_t, tmp);
  151. if(solve_status) { out = trans(out); powmat_status = true; }
  152. }
  153. if(powmat_status == false)
  154. {
  155. out.soft_reset();
  156. arma_stop_runtime_error("powmat(): transformation failed");
  157. }
  158. }
  159. //! @}