You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

235 lines
6.2KB

  1. ;******************************************************************************
  2. ;* linear least squares model
  3. ;*
  4. ;* Copyright (c) 2013 Loren Merritt
  5. ;*
  6. ;* This file is part of Libav.
  7. ;*
  8. ;* Libav is free software; you can redistribute it and/or
  9. ;* modify it under the terms of the GNU Lesser General Public
  10. ;* License as published by the Free Software Foundation; either
  11. ;* version 2.1 of the License, or (at your option) any later version.
  12. ;*
  13. ;* Libav is distributed in the hope that it will be useful,
  14. ;* but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. ;* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  16. ;* Lesser General Public License for more details.
  17. ;*
  18. ;* You should have received a copy of the GNU Lesser General Public
  19. ;* License along with Libav; if not, write to the Free Software
  20. ;* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  21. ;******************************************************************************
  22. %include "x86util.asm"
  23. SECTION .text
  24. %define MAX_VARS 32
  25. %define MAX_VARS_ALIGN (MAX_VARS+4)
  26. %define COVAR_STRIDE MAX_VARS_ALIGN*8
  27. %define COVAR(x,y) [covarq + (x)*8 + (y)*COVAR_STRIDE]
  28. struc LLSModel
  29. .covariance: resq MAX_VARS_ALIGN*MAX_VARS_ALIGN
  30. .coeff: resq MAX_VARS*MAX_VARS
  31. .variance: resq MAX_VARS
  32. .indep_count: resd 1
  33. endstruc
  34. %macro ADDPD_MEM 2
  35. %if cpuflag(avx)
  36. vaddpd %2, %2, %1
  37. %else
  38. addpd %2, %1
  39. %endif
  40. mova %1, %2
  41. %endmacro
  42. INIT_XMM sse2
  43. %define movdqa movaps
  44. cglobal update_lls, 2,5,8, ctx, var, i, j, covar2
  45. %define covarq ctxq
  46. mov id, [ctxq + LLSModel.indep_count]
  47. lea varq, [varq + iq*8]
  48. neg iq
  49. mov covar2q, covarq
  50. .loopi:
  51. ; Compute all 3 pairwise products of a 2x2 block that lies on the diagonal
  52. mova m1, [varq + iq*8]
  53. mova m3, [varq + iq*8 + 16]
  54. pshufd m4, m1, q1010
  55. pshufd m5, m1, q3232
  56. pshufd m6, m3, q1010
  57. pshufd m7, m3, q3232
  58. mulpd m0, m1, m4
  59. mulpd m1, m1, m5
  60. lea covarq, [covar2q + 16]
  61. ADDPD_MEM COVAR(-2,0), m0
  62. ADDPD_MEM COVAR(-2,1), m1
  63. lea jq, [iq + 2]
  64. cmp jd, -2
  65. jg .skip4x4
  66. .loop4x4:
  67. ; Compute all 16 pairwise products of a 4x4 block
  68. mulpd m0, m4, m3
  69. mulpd m1, m5, m3
  70. mulpd m2, m6, m3
  71. mulpd m3, m3, m7
  72. ADDPD_MEM COVAR(0,0), m0
  73. ADDPD_MEM COVAR(0,1), m1
  74. ADDPD_MEM COVAR(0,2), m2
  75. ADDPD_MEM COVAR(0,3), m3
  76. mova m3, [varq + jq*8 + 16]
  77. mulpd m0, m4, m3
  78. mulpd m1, m5, m3
  79. mulpd m2, m6, m3
  80. mulpd m3, m3, m7
  81. ADDPD_MEM COVAR(2,0), m0
  82. ADDPD_MEM COVAR(2,1), m1
  83. ADDPD_MEM COVAR(2,2), m2
  84. ADDPD_MEM COVAR(2,3), m3
  85. mova m3, [varq + jq*8 + 32]
  86. add covarq, 32
  87. add jq, 4
  88. cmp jd, -2
  89. jle .loop4x4
  90. .skip4x4:
  91. test jd, jd
  92. jg .skip2x4
  93. mulpd m4, m3
  94. mulpd m5, m3
  95. mulpd m6, m3
  96. mulpd m7, m3
  97. ADDPD_MEM COVAR(0,0), m4
  98. ADDPD_MEM COVAR(0,1), m5
  99. ADDPD_MEM COVAR(0,2), m6
  100. ADDPD_MEM COVAR(0,3), m7
  101. .skip2x4:
  102. add iq, 4
  103. add covar2q, 4*COVAR_STRIDE+32
  104. cmp id, -2
  105. jle .loopi
  106. test id, id
  107. jg .ret
  108. mov jq, iq
  109. %define covarq covar2q
  110. .loop2x1:
  111. movsd m0, [varq + iq*8]
  112. movlhps m0, m0
  113. mulpd m0, [varq + jq*8]
  114. ADDPD_MEM COVAR(0,0), m0
  115. inc iq
  116. add covarq, COVAR_STRIDE
  117. test id, id
  118. jle .loop2x1
  119. .ret:
  120. REP_RET
  121. INIT_YMM avx
  122. cglobal update_lls, 3,6,8, ctx, var, count, i, j, count2
  123. %define covarq ctxq
  124. mov countd, [ctxq + LLSModel.indep_count]
  125. lea count2d, [countq-2]
  126. xor id, id
  127. .loopi:
  128. ; Compute all 10 pairwise products of a 4x4 block that lies on the diagonal
  129. mova ymm1, [varq + iq*8]
  130. vbroadcastsd ymm4, [varq + iq*8]
  131. vbroadcastsd ymm5, [varq + iq*8 + 8]
  132. vbroadcastsd ymm6, [varq + iq*8 + 16]
  133. vbroadcastsd ymm7, [varq + iq*8 + 24]
  134. vextractf128 xmm3, ymm1, 1
  135. vmulpd ymm0, ymm1, ymm4
  136. vmulpd ymm1, ymm1, ymm5
  137. vmulpd xmm2, xmm3, xmm6
  138. vmulpd xmm3, xmm3, xmm7
  139. ADDPD_MEM COVAR(iq ,0), ymm0
  140. ADDPD_MEM COVAR(iq ,1), ymm1
  141. ADDPD_MEM COVAR(iq+2,2), xmm2
  142. ADDPD_MEM COVAR(iq+2,3), xmm3
  143. lea jd, [iq + 4]
  144. cmp jd, count2d
  145. jg .skip4x4
  146. .loop4x4:
  147. ; Compute all 16 pairwise products of a 4x4 block
  148. mova ymm3, [varq + jq*8]
  149. vmulpd ymm0, ymm3, ymm4
  150. vmulpd ymm1, ymm3, ymm5
  151. vmulpd ymm2, ymm3, ymm6
  152. vmulpd ymm3, ymm3, ymm7
  153. ADDPD_MEM COVAR(jq,0), ymm0
  154. ADDPD_MEM COVAR(jq,1), ymm1
  155. ADDPD_MEM COVAR(jq,2), ymm2
  156. ADDPD_MEM COVAR(jq,3), ymm3
  157. add jd, 4
  158. cmp jd, count2d
  159. jle .loop4x4
  160. .skip4x4:
  161. cmp jd, countd
  162. jg .skip2x4
  163. mova xmm3, [varq + jq*8]
  164. vmulpd xmm0, xmm3, xmm4
  165. vmulpd xmm1, xmm3, xmm5
  166. vmulpd xmm2, xmm3, xmm6
  167. vmulpd xmm3, xmm3, xmm7
  168. ADDPD_MEM COVAR(jq,0), xmm0
  169. ADDPD_MEM COVAR(jq,1), xmm1
  170. ADDPD_MEM COVAR(jq,2), xmm2
  171. ADDPD_MEM COVAR(jq,3), xmm3
  172. .skip2x4:
  173. add id, 4
  174. add covarq, 4*COVAR_STRIDE
  175. cmp id, count2d
  176. jle .loopi
  177. cmp id, countd
  178. jg .ret
  179. mov jd, id
  180. .loop2x1:
  181. vmovddup xmm0, [varq + iq*8]
  182. vmulpd xmm0, [varq + jq*8]
  183. ADDPD_MEM COVAR(jq,0), xmm0
  184. inc id
  185. add covarq, COVAR_STRIDE
  186. cmp id, countd
  187. jle .loop2x1
  188. .ret:
  189. REP_RET
  190. INIT_XMM sse2
  191. cglobal evaluate_lls, 3,4,2, ctx, var, order, i
  192. ; This function is often called on the same buffer as update_lls, but with
  193. ; an offset. They can't both be aligned.
  194. ; Load halves rather than movu to avoid store-forwarding stalls, since the
  195. ; input was initialized immediately prior to this function using scalar math.
  196. %define coefsq ctxq
  197. mov id, orderd
  198. imul orderd, MAX_VARS
  199. lea coefsq, [ctxq + LLSModel.coeff + orderq*8]
  200. movsd m0, [varq]
  201. movhpd m0, [varq + 8]
  202. mulpd m0, [coefsq]
  203. lea coefsq, [coefsq + iq*8]
  204. lea varq, [varq + iq*8]
  205. neg iq
  206. add iq, 2
  207. .loop:
  208. movsd m1, [varq + iq*8]
  209. movhpd m1, [varq + iq*8 + 8]
  210. mulpd m1, [coefsq + iq*8]
  211. addpd m0, m1
  212. add iq, 2
  213. jl .loop
  214. jg .skip1
  215. movsd m1, [varq + iq*8]
  216. mulsd m1, [coefsq + iq*8]
  217. addpd m0, m1
  218. .skip1:
  219. movhlps m1, m0
  220. addsd m0, m1
  221. %if ARCH_X86_32
  222. movsd r0m, m0
  223. fld qword r0m
  224. %endif
  225. RET