How to Guides -- How to optimize GEMM on CPU

最新推荐文章于 2024-07-29 17:17:10 发布
最新推荐文章于 2024-07-29 17:17:10 发布
阅读量297 收藏
点赞数
分类专栏： TVM文档翻译 文章标签： TVM 深度学习编译器
原文链接：https://tvm.apache.org/docs/how_to/optimize_operators/opt_gemm.html#sphx-glr-how-to-optimize-operators-opt-gemm-py
版权

Author: Jian Weng, Ruofei Yu

TVM 提供抽象接口，允许用户分别描述算法和算法的实现组织（所谓的调度）。 通常，以高性能调度编写算法会破坏算法的可读性和模块化。 此外，尝试各种看似有希望的时间表非常耗时。 在 TVM 的帮助下，我们可以有效地尝试这些时间表以提高性能。

在本教程中，我们将演示如何使用 TVM 优化方阵乘法，并通过简单地添加 18 行额外的代码实现比基线快 200 倍。

在 CPU 上执行的密集计算应用程序有两个重要的优化：

  • 提高内存访问的缓存命中率。 高缓存命中率可以加速复杂的数值计算和热点内存访问。 这需要我们将原始内存访问模式转换为适合缓存策略的模式。
  • SIMD（单指令多数据），或者我们称之为向量处理单元。 每次都会处理一小批数据，而不是单个网格。 这需要我们将循环体中的数据访问模式转换为统一模式，以便 LLVM 后端可以将其降低到 SIMD。

实际上，本教程中使用的所有方法都是这个 repo 中提到的技巧的一个子集。 其中一些已被 TVM 抽象自动应用，但由于 TVM 的限制，其中一些不能简单地应用。
下面提到的所有实验结果，都是在配备 Intel i7-4770HQ CPU 的 2015 年 15 英寸 MacBook 上执行的。 所有 x86 CPU 的高速缓存行大小应为 64 字节。

Preparation and Baseline

在本教程中，我们将演示如何使用 TVM 优化矩阵乘法。 在实际演示之前，我们首先定义这些变量。 然后我们编写一个基线实现，这是在 TVM 中编写矩阵乘法的最简单方法。

import tvm
import tvm.testing
from tvm import te
import numpy
import timeit

# The size of the matrix
# (M, K) x (K, N)
# You are free to try out different shapes, sometimes TVM optimization outperforms numpy with MKL.
M = 1024
K = 1024
N = 1024

# The default tensor type in tvm
dtype = "float32"

# using Intel AVX2(Advanced Vector Extensions) ISA for SIMD
# To get the best performance, please change the following line
# to llvm -mcpu=core-avx2, or specific type of CPU you use
target = "llvm"
dev = tvm.device(target, 0)

# Random generated tensor for testing
a = tvm.nd.array(numpy.random.rand(M, K).astype(dtype), dev)
b = tvm.nd.array(numpy.random.rand(K, N).astype(dtype), dev)

np_repeat = 100
np_runing_time = timeit.timeit(
    setup="import numpy\n"
    "M = " + str(M) + "\n"
    "K = " + str(K) + "\n"
    "N = " + str(N) + "\n"
    'dtype = "float32"\n'
    "a = numpy.random.rand(M, K).astype(dtype)\n"
    "b = numpy.random.rand(K, N).astype(dtype)\n",
    stmt="answer = numpy.dot(a, b)",
    number=np_repeat,
)
print("Numpy running time: %f" % (np_runing_time / np_repeat))

answer = numpy.dot(a.numpy(), b.numpy())

# Algorithm
k = te.reduce_axis((0, K), "k")
A = te.placeholder((M, K), name="A")
B = te.placeholder((K, N), name="B")
C = te.compute((M, N), lambda m, n: te.sum(A[m, k] * B[k, n], axis=k), name="C")

# Default schedule
s = te.create_schedule(C.op)
func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=1)
print("Baseline: %f" % evaluator(a, b, c).mean)

Out:
Numpy running time: 0.014836
Baseline: 3.256232

在 TVM 中，我们始终可以检查较低级别的 IR 以调试或优化我们的调度。 这是使用我们的基线计划生成的 IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  for (m: int32, 0, 1024) {
    for (n: int32, 0, 1024) {
      C_2[((m*1024) + n)] = 0f32
      for (k: int32, 0, 1024) {
        let cse_var_2: int32 = (m*1024)
        let cse_var_1: int32 = (cse_var_2 + n)
        C_2[cse_var_1] = ((float32*)C_2[cse_var_1] + ((float32*)A_2[(cse_var_2 + k)]*(float32*)B_2[((k*1024) + n)]))
      }
    }
  }
}

Blocking

提高缓存命中率的一个重要技巧是阻塞——数据块将逐块计算。 块内的内存访问是一个具有高内存局部性的小邻域。 在本教程中，我选择了 32 作为阻塞因子。 因此该块将填充 32 * 32 * sizeof(float) ，即总大小为 32KB 的缓存中的 4KB（L1 数据缓存）

bn = 32
kfactor = 4
s = te.create_schedule(C.op)

# Blocking by loop tiling
mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)
(kaxis,) = s[C].op.reduce_axis
ko, ki = s[C].split(kaxis, factor=kfactor)

# Hoist reduction domain outside the blocking loop
s[C].reorder(mo, no, ko, ki, mi, ni)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

# By simply tiling the loop 32x32, and hoisting ko, ki outside the blocking loops,
# we can see big speedup compared with the baseline.
evaluator = func.time_evaluator(func.entry_name, dev, number=10)
print("Opt1: %f" % evaluator(a, b, c).mean)

Opt1: 0.289397

这是阻塞后生成的IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  for (m.outer: int32, 0, 32) {
    for (n.outer: int32, 0, 32) {
      for (m.inner.init: int32, 0, 32) {
        for (n.inner.init: int32, 0, 32) {
          C_2[((((m.outer*32768) + (m.inner.init*1024)) + (n.outer*32)) + n.inner.init)] = 0f32
        }
      }
      for (k.outer: int32, 0, 256) {
        for (k.inner: int32, 0, 4) {
          for (m.inner: int32, 0, 32) {
            for (n.inner: int32, 0, 32) {
              let cse_var_3: int32 = (n.outer*32)
              let cse_var_2: int32 = ((m.outer*32768) + (m.inner*1024))
              let cse_var_1: int32 = ((cse_var_2 + cse_var_3) + n.inner)
              C_2[cse_var_1] = ((float32*)C_2[cse_var_1] + ((float32*)A_2[((cse_var_2 + (k.outer*4)) + k.inner)]*(float32*)B_2[((((k.outer*4096) + (k.inner*1024)) + cse_var_3) + n.inner)]))
            }
          }
        }
      }
    }
  }
}

Vectorization

另一个重要技巧是矢量化。 当内存访问模式一致时，编译器可以检测到这种模式并将连续内存传递给向量处理器。 在 TVM 中，我们可以使用 vectorize 接口来提示编译器这种模式，这样我们就可以大大加速它。
在本教程中，我们选择向量化内部循环行数据，因为它是缓存友好的。

s = te.create_schedule(C.op)
mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)
(kaxis,) = s[C].op.reduce_axis
ko, ki = s[C].split(kaxis, factor=kfactor)

s[C].reorder(mo, no, ko, ki, mi, ni)

# Vectorization
s[C].vectorize(ni)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=10)
print("Opt2: %f" % evaluator(a, b, c).mean)

Opt2: 0.328705

这是矢量化后生成的 IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  for (m.outer: int32, 0, 32) {
    for (n.outer: int32, 0, 32) {
      for (m.inner.init: int32, 0, 32) {
        C_2[ramp((((m.outer*32768) + (m.inner.init*1024)) + (n.outer*32)), 1, 32)] = broadcast(0f32, 32)
      }
      for (k.outer: int32, 0, 256) {
        for (k.inner: int32, 0, 4) {
          for (m.inner: int32, 0, 32) {
            let cse_var_3: int32 = (n.outer*32)
            let cse_var_2: int32 = ((m.outer*32768) + (m.inner*1024))
            let cse_var_1: int32 = (cse_var_2 + cse_var_3)
            C_2[ramp(cse_var_1, 1, 32)] = ((float32x32*)C_2[ramp(cse_var_1, 1, 32)] + (broadcast((float32*)A_2[((cse_var_2 + (k.outer*4)) + k.inner)], 32)*(float32x32*)B_2[ramp((((k.outer*4096) + (k.inner*1024)) + cse_var_3), 1, 32)]))
          }
        }
      }
    }
  }
}

Loop Permutation

如果我们查看上面的 IR，我们可以看到内部循环行数据针对 B 和 C 进行了矢量化。接下来我们将查看 A 的访问模式。在当前调度中，A 被逐列访问，即 不缓存友好。 如果我们改变 ki 和内轴 mi 的嵌套循环顺序，A 矩阵的访问模式对缓存更友好。

s = te.create_schedule(C.op)
mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)
(kaxis,) = s[C].op.reduce_axis
ko, ki = s[C].split(kaxis, factor=kfactor)

# re-ordering
s[C].reorder(mo, no, ko, mi, ki, ni)
s[C].vectorize(ni)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=10)
print("Opt3: %f" % evaluator(a, b, c).mean)

Opt3: 0.113126

这是循环置换后生成的 IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  for (m.outer: int32, 0, 32) {
    for (n.outer: int32, 0, 32) {
      for (m.inner.init: int32, 0, 32) {
        C_2[ramp((((m.outer*32768) + (m.inner.init*1024)) + (n.outer*32)), 1, 32)] = broadcast(0f32, 32)
      }
      for (k.outer: int32, 0, 256) {
        for (m.inner: int32, 0, 32) {
          for (k.inner: int32, 0, 4) {
            let cse_var_3: int32 = (n.outer*32)
            let cse_var_2: int32 = ((m.outer*32768) + (m.inner*1024))
            let cse_var_1: int32 = (cse_var_2 + cse_var_3)
            C_2[ramp(cse_var_1, 1, 32)] = ((float32x32*)C_2[ramp(cse_var_1, 1, 32)] + (broadcast((float32*)A_2[((cse_var_2 + (k.outer*4)) + k.inner)], 32)*(float32x32*)B_2[ramp((((k.outer*4096) + (k.inner*1024)) + cse_var_3), 1, 32)]))
          }
        }
      }
    }
  }
}

Array Packing

另一个重要的技巧是数组打包。 诀窍是对多维数组的存储进行重新排序，以便在展平并存储在一维内存中之后顺序访问它。
在这里插入图片描述

注意：此图是阵列打包如何工作的一般说明。

我们可以使用数组打包来解决 B 的访问模式。观察展平后 B 的数组访问模式，当我们迭代 K 维时，它不是顺序的。 我们可以用维度 [K][N] 对 B 重新排序，使其具有 [N/bn][K][bn] 维度，其中 bn 是阻塞因子，也是内循环中 B 的向量大小。 这种重新排序将 N 拆分为两个维度 - bigN (N/bn) 和 littleN (bn) - 新维度 [N/bn][K][bn] 匹配 B 从外部到内部循环的索引（no, ko, ki, ni) 在展平后导致 B 的顺序访问模式。

# We have to re-write the algorithm slightly.
packedB = te.compute(
    (N / bn, K, bn), lambda bigN, k, littleN: B[k, bigN * bn + littleN], name="packedB"
)
C = te.compute(
    (M, N),
    lambda m, n: te.sum(A[m, k] * packedB[n // bn, k, tvm.tir.indexmod(n, bn)], axis=k),
    name="C",
)

s = te.create_schedule(C.op)

mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)
(kaxis,) = s[C].op.reduce_axis
ko, ki = s[C].split(kaxis, factor=kfactor)

s[C].reorder(mo, no, ko, mi, ki, ni)
s[C].vectorize(ni)

bigN, _, littleN = s[packedB].op.axis
s[packedB].vectorize(littleN)
s[packedB].parallel(bigN)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=10)
print("Opt4: %f" % evaluator(a, b, c).mean)

Opt4: 0.106034

这是数组打包后生成的 IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  allocate(packedB: Pointer(global float32x32), float32x32, [32768]), storage_scope = global {
    for (bigN: int32, 0, 32) "parallel" {
      for (k: int32, 0, 1024) {
        packedB[ramp(((bigN*32768) + (k*32)), 1, 32)] = (float32x32*)B_2[ramp(((k*1024) + (bigN*32)), 1, 32)]
      }
    }
    for (m.outer: int32, 0, 32) {
      for (n.outer: int32, 0, 32) {
        for (m.inner.init: int32, 0, 32) {
          C_2[ramp((((m.outer*32768) + (m.inner.init*1024)) + (n.outer*32)), 1, 32)] = broadcast(0f32, 32)
        }
        for (k.outer: int32, 0, 256) {
          for (m.inner: int32, 0, 32) {
            for (k.inner: int32, 0, 4) {
              let cse_var_2: int32 = ((m.outer*32768) + (m.inner*1024))
              let cse_var_1: int32 = (cse_var_2 + (n.outer*32))
              C_2[ramp(cse_var_1, 1, 32)] = ((float32x32*)C_2[ramp(cse_var_1, 1, 32)] + (broadcast((float32*)A_2[((cse_var_2 + (k.outer*4)) + k.inner)], 32)*(float32x32*)packedB[ramp((((n.outer*32768) + (k.outer*128)) + (k.inner*32)), 1, 32)]))
            }
          }
        }
      }
    }
  }
}

Write cache for block

阻塞后，程序会逐块将结果写入C，访问模式不是顺序的。 因此，我们可以使用顺序缓存数组来保存块结果，并在所有块结果准备好时写入 C。

s = te.create_schedule(C.op)

# Allocate write cache
CC = s.cache_write(C, "global")

mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)

# Write cache is computed at no
s[CC].compute_at(s[C], no)

# New inner axes
mc, nc = s[CC].op.axis

(kaxis,) = s[CC].op.reduce_axis
ko, ki = s[CC].split(kaxis, factor=kfactor)
s[CC].reorder(ko, mc, ki, nc)
s[CC].vectorize(nc)

# TODO: Add separate optimization step to discuss loop unrolloing
# unrolling is a loop optimization strategy which can reduce branch
# prediction failures and increases the chance of concurrent execution
# unroll kfactor loops
s[CC].unroll(ki)

bigN, _, littleN = s[packedB].op.axis
s[packedB].vectorize(littleN)
s[packedB].parallel(bigN)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=10)
print("Opt5: %f" % evaluator(a, b, c).mean)

Opt5: 0.096499

这是阻塞后生成的IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  allocate(packedB: Pointer(global float32x32), float32x32, [32768]), storage_scope = global;
  allocate(C.global: Pointer(global float32), float32, [1024]), storage_scope = global {
    for (bigN: int32, 0, 32) "parallel" {
      for (k: int32, 0, 1024) {
        packedB[ramp(((bigN*32768) + (k*32)), 1, 32)] = (float32x32*)B_2[ramp(((k*1024) + (bigN*32)), 1, 32)]
      }
    }
    for (m.outer: int32, 0, 32) {
      for (n.outer: int32, 0, 32) {
        for (m.c.init: int32, 0, 32) {
          C.global[ramp((m.c.init*32), 1, 32)] = broadcast(0f32, 32)
        }
        for (k.outer: int32, 0, 256) {
          for (m.c: int32, 0, 32) {
            let cse_var_3: int32 = (m.c*32)
            let cse_var_2: int32 = ((n.outer*32768) + (k.outer*128))
            let cse_var_1: int32 = (((m.outer*32768) + (m.c*1024)) + (k.outer*4))
             {
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[cse_var_1], 32)*(float32x32*)packedB[ramp(cse_var_2, 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 1)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 32), 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 2)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 64), 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 3)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 96), 1, 32)]))
            }
          }
        }
        for (m.inner: int32, 0, 32) {
          for (n.inner: int32, 0, 32) {
            C_2[((((m.outer*32768) + (m.inner*1024)) + (n.outer*32)) + n.inner)] = (float32*)C.global[((m.inner*32) + n.inner)]
          }
        }
      }
    }
  }
}

Parallel

此外，我们还可以利用多核处理器进行线程级并行化。

s = te.create_schedule(C.op)

CC = s.cache_write(C, "global")

mo, no, mi, ni = s[C].tile(C.op.axis[0], C.op.axis[1], bn, bn)

s[CC].compute_at(s[C], no)

mc, nc = s[CC].op.axis

(kaxis,) = s[CC].op.reduce_axis
ko, ki = s[CC].split(kaxis, factor=kfactor)
s[CC].reorder(ko, mc, ki, nc)
s[CC].vectorize(nc)
s[CC].unroll(ki)

# parallel
s[C].parallel(mo)

bigN, _, littleN = s[packedB].op.axis
s[packedB].vectorize(littleN)
s[packedB].parallel(bigN)

func = tvm.build(s, [A, B, C], target=target, name="mmult")
assert func

c = tvm.nd.array(numpy.zeros((M, N), dtype=dtype), dev)
func(a, b, c)
tvm.testing.assert_allclose(c.numpy(), answer, rtol=1e-5)

evaluator = func.time_evaluator(func.entry_name, dev, number=50)
opt6_time = evaluator(a, b, c).mean
print("Opt6: %f" % opt6_time)

Opt6: 0.125081

这是并行化后生成的 IR。

print(tvm.lower(s, [A, B, C], simple_mode=True))

@main = primfn(A_1: handle, B_1: handle, C_1: handle) -> ()
  attr = {"from_legacy_te_schedule": True, "global_symbol": "main", "tir.noalias": True}
  buffers = {C: Buffer(C_2: Pointer(float32), float32, [1024, 1024], []),
             A: Buffer(A_2: Pointer(float32), float32, [1024, 1024], []),
             B: Buffer(B_2: Pointer(float32), float32, [1024, 1024], [])}
  buffer_map = {A_1: A, B_1: B, C_1: C} {
  allocate(packedB: Pointer(global float32x32), float32x32, [32768]), storage_scope = global {
    for (bigN: int32, 0, 32) "parallel" {
      for (k: int32, 0, 1024) {
        packedB[ramp(((bigN*32768) + (k*32)), 1, 32)] = (float32x32*)B_2[ramp(((k*1024) + (bigN*32)), 1, 32)]
      }
    }
    for (m.outer: int32, 0, 32) "parallel" {
      allocate(C.global: Pointer(global float32), float32, [1024]), storage_scope = global;
      for (n.outer: int32, 0, 32) {
        for (m.c.init: int32, 0, 32) {
          C.global[ramp((m.c.init*32), 1, 32)] = broadcast(0f32, 32)
        }
        for (k.outer: int32, 0, 256) {
          for (m.c: int32, 0, 32) {
            let cse_var_3: int32 = (m.c*32)
            let cse_var_2: int32 = ((n.outer*32768) + (k.outer*128))
            let cse_var_1: int32 = (((m.outer*32768) + (m.c*1024)) + (k.outer*4))
             {
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[cse_var_1], 32)*(float32x32*)packedB[ramp(cse_var_2, 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 1)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 32), 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 2)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 64), 1, 32)]))
              C.global[ramp(cse_var_3, 1, 32)] = ((float32x32*)C.global[ramp(cse_var_3, 1, 32)] + (broadcast((float32*)A_2[(cse_var_1 + 3)], 32)*(float32x32*)packedB[ramp((cse_var_2 + 96), 1, 32)]))
            }
          }
        }
        for (m.inner: int32, 0, 32) {
          for (n.inner: int32, 0, 32) {
            C_2[((((m.outer*32768) + (m.inner*1024)) + (n.outer*32)) + n.inner)] = (float32*)C.global[((m.inner*32) + n.inner)]
          }
        }
      }
    }
  }
}

Summary

应用上述简单优化后，仅用 18 行代码，我们生成的代码可以达到使用 MKL 的 60% 的 numpy 性能。 请注意，网页上的输出反映了非专有 Docker 容器上的运行时间，因此它们是不可靠的。 强烈建议您自己运行本教程，以观察 TVM 所获得的性能提升。

姆克儿
关注
专栏目录
如何在CPU上优化GEMM（上）
吴建明wujianming_110117
12-09 658
如何在CPU上优化GEMM（上） （TL；DR）TVM提供了抽象接口，用户分别描述算法和算法的实现组织（所谓的调度）。通常，在高性能调度中编写算法会破坏算法的可读性和模块性。尝试各种看似有希望的时间表是很耗时的。在TVM的帮助下，可以有效地尝试这些调度来提高性能。 本文将演示如何使用TVM优化平方矩阵乘法，并通过简单地添加18行额外的代码来实现比baseline基线快200倍的速度。 在CPU上执行的高强度计算应用程序有两个重要的优化： 提高内存访问的缓存命中率。高速缓存命中率可以加速复杂的数值计算和热点内
【TVM 巡礼】How to optimize cpu(x86) gemm串讲
I good vegetable a!
05-02 958
0x0. 前言 本文主要梳理一下在21年接触到优化gemm的知识，做一个学习总结。行文的顺序大概为： 介绍本文依赖的硬件环境和本文要完成的任务。 回顾gflops的计算方法然后计算本地环境的gflops。 How to optimize gemm on cpu ...
【darknet源码解析-06】gemm.h和gemm.c解析
最新发布
生活需要深度
07-29 756
@param K A的列数（不做转置）或者A^T的列数（做转置）,B的行数（不做转置）或者B^T行数（做转置），此处A未转置，B转置，故为A^T的列数，B的行数。* @param K A的列数（不做转置）或者A^T的列数（做转置）,B的行数（不做转置）或者B^T（做转置），此处A转置，B转置，故为A^T的列数，B^T的行数。* @param K A的列数（不做转置）或者A^T的列数（做转置）,B的行数（不做转置）或者B^T（做转置），此处A未转置，B转置，故为A的列数，B^T的行数。
基于how-to-optimize-gemm初探矩阵乘法优化
I good vegetable a!
11-05 1218
1. 前言 这次，我们来聊一个轻松一点的话题，那就是给你一个矩阵A和一个矩阵B，使用矩阵乘法获得目标矩阵C，相信大家都不难写出下面的代码： #define A( i, j ) a[ (i)*lda + (j) ] #define B( i, j ) b[ (i)*ldb + (j) ] #define C( i, j ) c[ (i)*ldc + (j) ] // gemm C = A * B + C void MatrixMultiply(int m, int n, int k, float *a, in
探索极致性能：how-to-optimize-gemm — GEMM优化教程与实战
gitblog_00091的博客
05-21 661
探索极致性能：how-to-optimize-gemm — GEMM优化教程与实战 how-to-optimize-gemmrow-major matmul optimization项目地址:https://gitcode.com/gh_mirrors/how/how-to-optimize-gemm 项目简介 how-to-optimize-gemm 是一个全面的矩阵乘法（Generalize...
如何加速矩阵乘法——优化GEMM (CPU单线程篇)
sinat_22510827的博客
07-22 6049
矩阵乘法GEMM(Generalmatrixmultiply)是一个被广泛使用的基础算法，各种领域都需要应用，例如神经网络的核心计算任务就是矩阵乘法，交易中的各种信号计算也可能用到矩阵乘法。因此矩阵乘法的效率是极其关键的。关于如何优化矩阵乘法，我准备写一个较短的系列博文，包括CPU单线程篇、CPU多线程篇、GPU篇。原本计划还有一个稀疏矩阵乘法篇，由于这学期毕业前也没有时间把GPU篇做到满意，因此，稀疏矩阵篇就没了，GPU篇也很不完整。如果日后有时间且有那冲动可能会补上（大概率没有）。...
how-to-optimize-gemm
05-29
How To Optimize GEMM 行主序的GEMM优化方案 x86 cpufp 克隆高叔叔的X86浮点峰值测试工具，使用方法见README。 本机（Intel(R) Xeon(R) CPU E5-2678 v3 @ 2.50GHz 12 Cores）测试结果如下： Thread(s): 1 fma fp32 ...
bp网络实现matlab代码-Imporved-weed-algorithm-to-optimize-BP-network:采用改进的杂草算法
05-28
Imporved-weed-algorithm-to-optimize-BP-network In the matlab environment, by combining the differential evolution algorithm and the weed optimization algorithm, an improved differential evolution weed...
How To Optimize Content For Featured Snippets-crx插件
03-19
如何优化精选片段的内容 清单指南，可帮助您优化精选片段的内容并增加内容访问量。 支持语言:English
genetic-to-optimize-lifetime-of-wsn.rar_GA wsn_The Network_WSN
07-14
there is a new approach Using in Wireless sensor network to Clustering, in this approach we use GA algorithm to optimize the lifetime of network.
如何针对奔腾微处理器进行优化How to Optimize for the Pentium Microprocessors
11-15
本手册介绍了有关为英特尔微处理器优化代码所需的一切，包括奔腾，奔腾MMX，奔腾Pro，奔腾II，奔腾III，奔腾4，赛扬和至强。
如何优化CPU GEMM
四哥的专栏
05-24 201
截止到目前为止，已经具有和两份实现，虽然在性能上有一定的优化，但距离真正的高性能计算库的要求还相差甚远。即使抛开性能问题不谈，目前 optimize_row_major_sgemm 也很难视为一个合格的库函数，因为该函数在接口定义上太过随意，别人很难直接复用。众所周知，矩阵乘优化已经是非常成熟的课题了，其中自然衍生了许多标准，以方便不同开发者或者研究人员之间工作的交流和复用，其中最基础的便是 BLAS接口规范。是一系列基本线性代数运算函数1的接口（interface）标准。
GEMM算法及优化流程详解
热门推荐
独行侠
02-16 2万+
前言 神经网络前向耗时主要由卷积的耗时决定，参考賈杨青毕业论文，那么如何对卷积加速便成了重要的一个点，主流的加速方法有 以下几种： im2col+GEMM：目前几乎所有的主流计算框架包括 Caffe, MXNet 等都实现了该方法. 该方法把整个卷积过程转化成了GEMM过程，而GEMM在各种 BLAS 库中都是被极致优化的，一般来说，速度较快。 Winograd： Winograd 是存在...
如何在 CPU 上优化 GEMM
吴建明wujianming_110117
10-05 399
如何在 CPU 上优化 GEMM (TL;DR) TVM 提供抽象接口，允许用户分别描述算法和算法的实施组织（所谓的调度）。通常，在高性能调度中编写算法，会破坏算法的可读性和模块化。尝试各种看似有前途的调度也很耗时。在 TVM 的帮助下，可以有效地尝试这些调度，提高性能。 将演示如何使用 TVM 优化矩阵乘法，通过简单地添加 18 行额外代码，实现比基线快 200 倍。 在 CPU 上执行的密集计算应用程序，有两个重要的优化： 提高内存访问的缓存命中率。复杂的数值计算和热点内存访问，都可以通过高缓存命中率
tvm tutorial (1.2)
@fei
02-09 221
紧接着上一篇实现了array packing对矩阵乘法进行优化的策略，接下来记录cache写回优化（block writing throung caching) cache的写回说头比较多，有两种情况，一种叫write back，另一种是write through，详细内容可以参考下面链接： https://blog.csdn.net/xingzhe22222/article/details/81988101 矩阵C的写回操作也是跳跃式的，不利于 ...
darknet中gemm实现
永远飞翔的鸟
10-24 707
darknet对卷积计算的处理实际上是： 先使用im2col将input_channel*(height*width)的输入特征图（实际存储是按照行存储的，即是1*(input_channel*height*width)的一维数组）转化成(input_channel*kernel_size*kernel_size)*(out_height*out_width)的特征矩阵，这里同样是按行存储的。 之后通过gemm函数实现通用的矩阵乘法实现卷积计算，即让卷积核矩阵*im2col后的输入特征矩阵。其中卷积核的
caffe中的矩阵运算函数caffe_cpu_gemm,cblas_sgemm等解析
小小码农的博客
05-31 5835
caffe中最典型且常用的卷积运算，是通过将卷积操作转化成矩阵乘法来实现的，因此，卷积层的一系列程序实际上就是在为矩阵的卷积式展开和矩阵乘法函数做准备，caffe_cpu_gemm也就是在调用矩阵乘法函数cblas_sgemm。1.矩阵间乘法：这里我们先来看一下这些矩阵乘法函数有输入需求，以单精度矩阵乘法函数cblas_sgemm为例：cblas_sgemm(CblasRowMajor, Cbla...
how to learn linux
05-16
There are many ways to learn Linux, but here are some steps to get started: 1. Choose a Linux distribution: There are many distributions to choose from, such as Ubuntu, Debian, Fedora, or CentOS. Choose one that best suits your needs and interests. 2. Install Linux: You can install Linux on your computer as a dual-boot system, on a virtual machine, or on a dedicated machine. 3. Familiarize yourself with the Linux command line: Linux is often used through the command line interface (CLI), so it's important to learn basic commands such as ls, cd, mkdir, and rm. 4. Learn basic system administration tasks: Learn how to manage users, install packages, configure network settings, and monitor system performance. 5. Join a Linux community: Join forums, attend meetups, and contribute to open source projects to learn from others and get help when needed. 6. Practice, practice, practice: The best way to learn Linux is to use it regularly. Install software, create scripts, and optimize your system to gain experience and confidence. 7. Read books and online resources: There are many books and online resources available to help you learn Linux, such as "The Linux Command Line" by William Shotts and Linux documentation sites like Linux.org and Linux.com. Remember that learning Linux takes time and practice, so be patient and keep learning!
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值