fn_qr.hpp 3.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143
  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 fn_qr
  16. //! @{
  17. //! QR decomposition
  18. template<typename T1>
  19. inline
  20. bool
  21. qr
  22. (
  23. Mat<typename T1::elem_type>& Q,
  24. Mat<typename T1::elem_type>& R,
  25. const Base<typename T1::elem_type,T1>& X,
  26. const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
  27. )
  28. {
  29. arma_extra_debug_sigprint();
  30. arma_ignore(junk);
  31. arma_debug_check( (&Q == &R), "qr(): Q and R are the same object");
  32. const bool status = auxlib::qr(Q, R, X);
  33. if(status == false)
  34. {
  35. Q.soft_reset();
  36. R.soft_reset();
  37. arma_debug_warn("qr(): decomposition failed");
  38. }
  39. return status;
  40. }
  41. //! economical QR decomposition
  42. template<typename T1>
  43. inline
  44. bool
  45. qr_econ
  46. (
  47. Mat<typename T1::elem_type>& Q,
  48. Mat<typename T1::elem_type>& R,
  49. const Base<typename T1::elem_type,T1>& X,
  50. const typename arma_blas_type_only<typename T1::elem_type>::result* junk = 0
  51. )
  52. {
  53. arma_extra_debug_sigprint();
  54. arma_ignore(junk);
  55. arma_debug_check( (&Q == &R), "qr_econ(): Q and R are the same object");
  56. const bool status = auxlib::qr_econ(Q, R, X);
  57. if(status == false)
  58. {
  59. Q.soft_reset();
  60. R.soft_reset();
  61. arma_debug_warn("qr_econ(): decomposition failed");
  62. }
  63. return status;
  64. }
  65. //! QR decomposition with pivoting
  66. template<typename T1>
  67. inline
  68. typename enable_if2< is_supported_blas_type<typename T1::elem_type>::value, bool >::result
  69. qr
  70. (
  71. Mat<typename T1::elem_type>& Q,
  72. Mat<typename T1::elem_type>& R,
  73. Mat<uword>& P,
  74. const Base<typename T1::elem_type,T1>& X,
  75. const char* P_mode = "matrix"
  76. )
  77. {
  78. arma_extra_debug_sigprint();
  79. arma_debug_check( (&Q == &R), "qr(): Q and R are the same object");
  80. const char sig = (P_mode != NULL) ? P_mode[0] : char(0);
  81. arma_debug_check( ((sig != 'm') && (sig != 'v')), "qr(): argument 'P_mode' must be \"vector\" or \"matrix\"" );
  82. bool status = false;
  83. if(sig == 'v')
  84. {
  85. status = auxlib::qr_pivot(Q, R, P, X);
  86. }
  87. else
  88. if(sig == 'm')
  89. {
  90. Mat<uword> P_vec;
  91. status = auxlib::qr_pivot(Q, R, P_vec, X);
  92. if(status)
  93. {
  94. // construct P
  95. const uword N = P_vec.n_rows;
  96. P.zeros(N,N);
  97. for(uword row=0; row < N; ++row) { P.at(P_vec[row], row) = uword(1); }
  98. }
  99. }
  100. if(status == false)
  101. {
  102. Q.soft_reset();
  103. R.soft_reset();
  104. P.soft_reset();
  105. arma_debug_warn("qr(): decomposition failed");
  106. }
  107. return status;
  108. }
  109. //! @}