中文第一计算机图形学社区OpenGPU 版权所有2007-2018

 找回密码
 注册

扫一扫,访问微社区

搜索
查看: 350|回复: 3

OpenCl : clBlas 基于gtest框架测试fail

[复制链接]
发表于 2018-5-24 15:56:21 | 显示全部楼层 |阅读模式
大家好,请教个问题:

CUDA下clblas已经编译成功了,在做googletest的时候有一小部分fail了,主要fail在了ztrmm,“复数三角矩阵相乘”。

从代码看来是blas::trmm和clblas::trmm做的对比,前者调用了blas库中的函数没法看到,后者可以看到kernel。

由于test的时候矩阵比较大(63x63,128x128),所以我用小维度的矩阵试了下,发现不知道最后的结果是怎么算出来的,(和我想的“矩阵乘按矩阵乘的规则,具体计算时按照复数乘的规则”不一样)。

大家遇到过这个问题吗?具体fail的原因是因为clblas code错误吗?需要根据具体的平台修改clblas code吗?

一个fail实例,矩阵128x128,kernel有点复杂,不知道是怎么计算的?附上:
#ifdef cl_khr_fp64
#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#else
#pragma OPENCL EXTENSION cl_amd_fp64 : enable
#endif

typedef union GPtr {
    __global float *f;
    __global double *d;
    __global float2 *f2v;
    __global double2 *d2v;
    __global float4 *f4v;
    __global double4 *d4v;
    __global float8 *f8v;
    __global double8 *d8v;
    __global float16 *f16v;
    __global double16 *d16v;
} GPtr;

typedef union LPtr {
    __local float *f;
    __local double *d;
    __local float2 *f2v;
    __local double2 *d2v;
    __local float4 *f4v;
    __local double4 *d4v;
    __local float8 *f8v;
    __local double8 *d8v;
    __local float16 *f16v;
    __local double16 *d16v;
} LPtr;

typedef union PPtr {
    float *f;
    double *d;
    float2 *f2v;
    double2 *d2v;
    float4 *f4v;
    double4 *d4v;
    float8 *f8v;
    double8 *d8v;
    float16 *f16v;
    double16 *d16v;
} PPtr;

__attribute__((reqd_work_group_size(64, 1, 1)))
void __kernel
ztrmmBlock(
    uint M,
    uint N,
    double2 alpha,
    const __global double2 *restrict A,
    uint lda,
    const __global double2 *restrict B,
    __global double2 *C,
    uint ldb)
{
    double16 a0, a1;
    double4 b0, b1, b2, b3, b4, b5, b6, b7;
    double16 c0, c1, c2, c3, c4, c5, c6, c7;
    uint currM, currN;
    uint4 coord = 0; /* contains coordB, coordA, k */
    const int lid = get_local_id(0);
    const int gid = get_global_id(0) / 64;

    currN = gid * 32;
    currM = 0;

    coord.x = currN + lid % 4 * 8;
    for (uint m0 = 0; m0 < M; m0 += 128) {
        uint kBegin = currM;
        coord.z = kBegin;
        coord.y = currM + lid / 4 * 8;
        c0 = 0;
        c1 = 0;
        c2 = 0;
        c3 = 0;
        c4 = 0;
        c5 = 0;
        c6 = 0;
        c7 = 0;

        uint k0;
        __global double4 *pA = (__global double4 *)&A[mad24(coord.z, lda, coord.y)];
        __global double4 *pB = (__global double4 *)&B[mad24(coord.x, ldb, coord.z)];
        for (k0 = kBegin; (k0 <= (kBegin + 128u))&&(k0 < M); k0 += 2) {
            coord.z = k0;
            /* -- Tiles multiplier -- */
            b0 = pB[0];
            b1 = pB[(ldb >> 1)];
            b2 = pB[ldb];
            b3 = pB[mad24(3u, (ldb >> 1), 0u)];
            b4 = pB[(ldb << 1)];
            b5 = pB[mad24(5u, (ldb >> 1), 0u)];
            b6 = pB[mad24(6u, (ldb >> 1), 0u)];
            b7 = pB[mad24(7u, (ldb >> 1), 0u)];


            a0.s0123 = pA[0];
            a0.s4567 = pA[1];
            a0.s89ab = pA[2];
            a0.scdef = pA[3];
            a1.s0123 = pA[(lda >> 1)];
            a1.s4567 = pA[(lda >> 1) + 1];
            a1.s89ab = pA[(lda >> 1) + 2];
            a1.scdef = pA[(lda >> 1) + 3];

            // post fetch A
            {
                uint zy = coord.z;
                a0.s01 = zy < coord.y ? 0 : a0.s01;
                a0.s01 = zy == coord.y ? (double2)(1, 0) : a0.s01;
                a0.s23 = zy < coord.y + 1 ? 0 : a0.s23;
                a0.s23 = zy == coord.y + 1 ? (double2)(1, 0) : a0.s23;
                a0.s45 = zy < coord.y + 2 ? 0 : a0.s45;
                a0.s45 = zy == coord.y + 2 ? (double2)(1, 0) : a0.s45;
                a0.s67 = zy < coord.y + 3 ? 0 : a0.s67;
                a0.s67 = zy == coord.y + 3 ? (double2)(1, 0) : a0.s67;
                a0.s89 = zy < coord.y + 4 ? 0 : a0.s89;
                a0.s89 = zy == coord.y + 4 ? (double2)(1, 0) : a0.s89;
                a0.sab = zy < coord.y + 5 ? 0 : a0.sab;
                a0.sab = zy == coord.y + 5 ? (double2)(1, 0) : a0.sab;
                a0.scd = zy < coord.y + 6 ? 0 : a0.scd;
                a0.scd = zy == coord.y + 6 ? (double2)(1, 0) : a0.scd;
                a0.sef = zy < coord.y + 7 ? 0 : a0.sef;
                a0.sef = zy == coord.y + 7 ? (double2)(1, 0) : a0.sef;
                zy++;
                a1.s01 = zy < coord.y ? 0 : a1.s01;
                a1.s01 = zy == coord.y ? (double2)(1, 0) : a1.s01;
                a1.s23 = zy < coord.y + 1 ? 0 : a1.s23;
                a1.s23 = zy == coord.y + 1 ? (double2)(1, 0) : a1.s23;
                a1.s45 = zy < coord.y + 2 ? 0 : a1.s45;
                a1.s45 = zy == coord.y + 2 ? (double2)(1, 0) : a1.s45;
                a1.s67 = zy < coord.y + 3 ? 0 : a1.s67;
                a1.s67 = zy == coord.y + 3 ? (double2)(1, 0) : a1.s67;
                a1.s89 = zy < coord.y + 4 ? 0 : a1.s89;
                a1.s89 = zy == coord.y + 4 ? (double2)(1, 0) : a1.s89;
                a1.sab = zy < coord.y + 5 ? 0 : a1.sab;
                a1.sab = zy == coord.y + 5 ? (double2)(1, 0) : a1.sab;
                a1.scd = zy < coord.y + 6 ? 0 : a1.scd;
                a1.scd = zy == coord.y + 6 ? (double2)(1, 0) : a1.scd;
                a1.sef = zy < coord.y + 7 ? 0 : a1.sef;
                a1.sef = zy == coord.y + 7 ? (double2)(1, 0) : a1.sef;
            }

            c0.s01 += b0.s01 * a0.s0 + b0.s10 * (double2)(-a0.s1, a0.s1);
            c1.s01 += b1.s01 * a0.s0 + b1.s10 * (double2)(-a0.s1, a0.s1);
            c2.s01 += b2.s01 * a0.s0 + b2.s10 * (double2)(-a0.s1, a0.s1);
            c3.s01 += b3.s01 * a0.s0 + b3.s10 * (double2)(-a0.s1, a0.s1);
            c4.s01 += b4.s01 * a0.s0 + b4.s10 * (double2)(-a0.s1, a0.s1);
            c5.s01 += b5.s01 * a0.s0 + b5.s10 * (double2)(-a0.s1, a0.s1);
            c6.s01 += b6.s01 * a0.s0 + b6.s10 * (double2)(-a0.s1, a0.s1);
            c7.s01 += b7.s01 * a0.s0 + b7.s10 * (double2)(-a0.s1, a0.s1);
            c0.s23 += b0.s01 * a0.s2 + b0.s10 * (double2)(-a0.s3, a0.s3);
            c1.s23 += b1.s01 * a0.s2 + b1.s10 * (double2)(-a0.s3, a0.s3);
            c2.s23 += b2.s01 * a0.s2 + b2.s10 * (double2)(-a0.s3, a0.s3);
            c3.s23 += b3.s01 * a0.s2 + b3.s10 * (double2)(-a0.s3, a0.s3);
            c4.s23 += b4.s01 * a0.s2 + b4.s10 * (double2)(-a0.s3, a0.s3);
            c5.s23 += b5.s01 * a0.s2 + b5.s10 * (double2)(-a0.s3, a0.s3);
            c6.s23 += b6.s01 * a0.s2 + b6.s10 * (double2)(-a0.s3, a0.s3);
            c7.s23 += b7.s01 * a0.s2 + b7.s10 * (double2)(-a0.s3, a0.s3);
            c0.s45 += b0.s01 * a0.s4 + b0.s10 * (double2)(-a0.s5, a0.s5);
            c1.s45 += b1.s01 * a0.s4 + b1.s10 * (double2)(-a0.s5, a0.s5);
            c2.s45 += b2.s01 * a0.s4 + b2.s10 * (double2)(-a0.s5, a0.s5);
            c3.s45 += b3.s01 * a0.s4 + b3.s10 * (double2)(-a0.s5, a0.s5);
            c4.s45 += b4.s01 * a0.s4 + b4.s10 * (double2)(-a0.s5, a0.s5);
            c5.s45 += b5.s01 * a0.s4 + b5.s10 * (double2)(-a0.s5, a0.s5);
            c6.s45 += b6.s01 * a0.s4 + b6.s10 * (double2)(-a0.s5, a0.s5);
            c7.s45 += b7.s01 * a0.s4 + b7.s10 * (double2)(-a0.s5, a0.s5);
            c0.s67 += b0.s01 * a0.s6 + b0.s10 * (double2)(-a0.s7, a0.s7);
            c1.s67 += b1.s01 * a0.s6 + b1.s10 * (double2)(-a0.s7, a0.s7);
            c2.s67 += b2.s01 * a0.s6 + b2.s10 * (double2)(-a0.s7, a0.s7);
            c3.s67 += b3.s01 * a0.s6 + b3.s10 * (double2)(-a0.s7, a0.s7);
            c4.s67 += b4.s01 * a0.s6 + b4.s10 * (double2)(-a0.s7, a0.s7);
            c5.s67 += b5.s01 * a0.s6 + b5.s10 * (double2)(-a0.s7, a0.s7);
            c6.s67 += b6.s01 * a0.s6 + b6.s10 * (double2)(-a0.s7, a0.s7);
            c7.s67 += b7.s01 * a0.s6 + b7.s10 * (double2)(-a0.s7, a0.s7);
            c0.s89 += b0.s01 * a0.s8 + b0.s10 * (double2)(-a0.s9, a0.s9);
            c1.s89 += b1.s01 * a0.s8 + b1.s10 * (double2)(-a0.s9, a0.s9);
            c2.s89 += b2.s01 * a0.s8 + b2.s10 * (double2)(-a0.s9, a0.s9);
            c3.s89 += b3.s01 * a0.s8 + b3.s10 * (double2)(-a0.s9, a0.s9);
            c4.s89 += b4.s01 * a0.s8 + b4.s10 * (double2)(-a0.s9, a0.s9);
            c5.s89 += b5.s01 * a0.s8 + b5.s10 * (double2)(-a0.s9, a0.s9);
            c6.s89 += b6.s01 * a0.s8 + b6.s10 * (double2)(-a0.s9, a0.s9);
            c7.s89 += b7.s01 * a0.s8 + b7.s10 * (double2)(-a0.s9, a0.s9);
            c0.sab += b0.s01 * a0.sa + b0.s10 * (double2)(-a0.sb, a0.sb);
            c1.sab += b1.s01 * a0.sa + b1.s10 * (double2)(-a0.sb, a0.sb);
            c2.sab += b2.s01 * a0.sa + b2.s10 * (double2)(-a0.sb, a0.sb);
            c3.sab += b3.s01 * a0.sa + b3.s10 * (double2)(-a0.sb, a0.sb);
            c4.sab += b4.s01 * a0.sa + b4.s10 * (double2)(-a0.sb, a0.sb);
            c5.sab += b5.s01 * a0.sa + b5.s10 * (double2)(-a0.sb, a0.sb);
            c6.sab += b6.s01 * a0.sa + b6.s10 * (double2)(-a0.sb, a0.sb);
            c7.sab += b7.s01 * a0.sa + b7.s10 * (double2)(-a0.sb, a0.sb);
            c0.scd += b0.s01 * a0.sc + b0.s10 * (double2)(-a0.sd, a0.sd);
            c1.scd += b1.s01 * a0.sc + b1.s10 * (double2)(-a0.sd, a0.sd);
            c2.scd += b2.s01 * a0.sc + b2.s10 * (double2)(-a0.sd, a0.sd);
            c3.scd += b3.s01 * a0.sc + b3.s10 * (double2)(-a0.sd, a0.sd);
            c4.scd += b4.s01 * a0.sc + b4.s10 * (double2)(-a0.sd, a0.sd);
            c5.scd += b5.s01 * a0.sc + b5.s10 * (double2)(-a0.sd, a0.sd);
            c6.scd += b6.s01 * a0.sc + b6.s10 * (double2)(-a0.sd, a0.sd);
            c7.scd += b7.s01 * a0.sc + b7.s10 * (double2)(-a0.sd, a0.sd);
            c0.sef += b0.s01 * a0.se + b0.s10 * (double2)(-a0.sf, a0.sf);
            c1.sef += b1.s01 * a0.se + b1.s10 * (double2)(-a0.sf, a0.sf);
            c2.sef += b2.s01 * a0.se + b2.s10 * (double2)(-a0.sf, a0.sf);
            c3.sef += b3.s01 * a0.se + b3.s10 * (double2)(-a0.sf, a0.sf);
            c4.sef += b4.s01 * a0.se + b4.s10 * (double2)(-a0.sf, a0.sf);
            c5.sef += b5.s01 * a0.se + b5.s10 * (double2)(-a0.sf, a0.sf);
            c6.sef += b6.s01 * a0.se + b6.s10 * (double2)(-a0.sf, a0.sf);
            c7.sef += b7.s01 * a0.se + b7.s10 * (double2)(-a0.sf, a0.sf);

            c0.s01 += b0.s23 * a1.s0 + b0.s32 * (double2)(-a1.s1, a1.s1);
            c1.s01 += b1.s23 * a1.s0 + b1.s32 * (double2)(-a1.s1, a1.s1);
            c2.s01 += b2.s23 * a1.s0 + b2.s32 * (double2)(-a1.s1, a1.s1);
            c3.s01 += b3.s23 * a1.s0 + b3.s32 * (double2)(-a1.s1, a1.s1);
            c4.s01 += b4.s23 * a1.s0 + b4.s32 * (double2)(-a1.s1, a1.s1);
            c5.s01 += b5.s23 * a1.s0 + b5.s32 * (double2)(-a1.s1, a1.s1);
            c6.s01 += b6.s23 * a1.s0 + b6.s32 * (double2)(-a1.s1, a1.s1);
            c7.s01 += b7.s23 * a1.s0 + b7.s32 * (double2)(-a1.s1, a1.s1);
            c0.s23 += b0.s23 * a1.s2 + b0.s32 * (double2)(-a1.s3, a1.s3);
            c1.s23 += b1.s23 * a1.s2 + b1.s32 * (double2)(-a1.s3, a1.s3);
            c2.s23 += b2.s23 * a1.s2 + b2.s32 * (double2)(-a1.s3, a1.s3);
            c3.s23 += b3.s23 * a1.s2 + b3.s32 * (double2)(-a1.s3, a1.s3);
            c4.s23 += b4.s23 * a1.s2 + b4.s32 * (double2)(-a1.s3, a1.s3);
            c5.s23 += b5.s23 * a1.s2 + b5.s32 * (double2)(-a1.s3, a1.s3);
            c6.s23 += b6.s23 * a1.s2 + b6.s32 * (double2)(-a1.s3, a1.s3);
            c7.s23 += b7.s23 * a1.s2 + b7.s32 * (double2)(-a1.s3, a1.s3);
            c0.s45 += b0.s23 * a1.s4 + b0.s32 * (double2)(-a1.s5, a1.s5);
            c1.s45 += b1.s23 * a1.s4 + b1.s32 * (double2)(-a1.s5, a1.s5);
            c2.s45 += b2.s23 * a1.s4 + b2.s32 * (double2)(-a1.s5, a1.s5);
            c3.s45 += b3.s23 * a1.s4 + b3.s32 * (double2)(-a1.s5, a1.s5);
            c4.s45 += b4.s23 * a1.s4 + b4.s32 * (double2)(-a1.s5, a1.s5);
            c5.s45 += b5.s23 * a1.s4 + b5.s32 * (double2)(-a1.s5, a1.s5);
            c6.s45 += b6.s23 * a1.s4 + b6.s32 * (double2)(-a1.s5, a1.s5);
            c7.s45 += b7.s23 * a1.s4 + b7.s32 * (double2)(-a1.s5, a1.s5);
            c0.s67 += b0.s23 * a1.s6 + b0.s32 * (double2)(-a1.s7, a1.s7);
            c1.s67 += b1.s23 * a1.s6 + b1.s32 * (double2)(-a1.s7, a1.s7);
            c2.s67 += b2.s23 * a1.s6 + b2.s32 * (double2)(-a1.s7, a1.s7);
            c3.s67 += b3.s23 * a1.s6 + b3.s32 * (double2)(-a1.s7, a1.s7);
            c4.s67 += b4.s23 * a1.s6 + b4.s32 * (double2)(-a1.s7, a1.s7);
            c5.s67 += b5.s23 * a1.s6 + b5.s32 * (double2)(-a1.s7, a1.s7);
            c6.s67 += b6.s23 * a1.s6 + b6.s32 * (double2)(-a1.s7, a1.s7);
            c7.s67 += b7.s23 * a1.s6 + b7.s32 * (double2)(-a1.s7, a1.s7);
            c0.s89 += b0.s23 * a1.s8 + b0.s32 * (double2)(-a1.s9, a1.s9);
            c1.s89 += b1.s23 * a1.s8 + b1.s32 * (double2)(-a1.s9, a1.s9);
            c2.s89 += b2.s23 * a1.s8 + b2.s32 * (double2)(-a1.s9, a1.s9);
            c3.s89 += b3.s23 * a1.s8 + b3.s32 * (double2)(-a1.s9, a1.s9);
            c4.s89 += b4.s23 * a1.s8 + b4.s32 * (double2)(-a1.s9, a1.s9);
            c5.s89 += b5.s23 * a1.s8 + b5.s32 * (double2)(-a1.s9, a1.s9);
            c6.s89 += b6.s23 * a1.s8 + b6.s32 * (double2)(-a1.s9, a1.s9);
            c7.s89 += b7.s23 * a1.s8 + b7.s32 * (double2)(-a1.s9, a1.s9);
            c0.sab += b0.s23 * a1.sa + b0.s32 * (double2)(-a1.sb, a1.sb);
            c1.sab += b1.s23 * a1.sa + b1.s32 * (double2)(-a1.sb, a1.sb);
            c2.sab += b2.s23 * a1.sa + b2.s32 * (double2)(-a1.sb, a1.sb);
            c3.sab += b3.s23 * a1.sa + b3.s32 * (double2)(-a1.sb, a1.sb);
            c4.sab += b4.s23 * a1.sa + b4.s32 * (double2)(-a1.sb, a1.sb);
            c5.sab += b5.s23 * a1.sa + b5.s32 * (double2)(-a1.sb, a1.sb);
            c6.sab += b6.s23 * a1.sa + b6.s32 * (double2)(-a1.sb, a1.sb);
            c7.sab += b7.s23 * a1.sa + b7.s32 * (double2)(-a1.sb, a1.sb);
            c0.scd += b0.s23 * a1.sc + b0.s32 * (double2)(-a1.sd, a1.sd);
            c1.scd += b1.s23 * a1.sc + b1.s32 * (double2)(-a1.sd, a1.sd);
            c2.scd += b2.s23 * a1.sc + b2.s32 * (double2)(-a1.sd, a1.sd);
            c3.scd += b3.s23 * a1.sc + b3.s32 * (double2)(-a1.sd, a1.sd);
            c4.scd += b4.s23 * a1.sc + b4.s32 * (double2)(-a1.sd, a1.sd);
            c5.scd += b5.s23 * a1.sc + b5.s32 * (double2)(-a1.sd, a1.sd);
            c6.scd += b6.s23 * a1.sc + b6.s32 * (double2)(-a1.sd, a1.sd);
            c7.scd += b7.s23 * a1.sc + b7.s32 * (double2)(-a1.sd, a1.sd);
            c0.sef += b0.s23 * a1.se + b0.s32 * (double2)(-a1.sf, a1.sf);
            c1.sef += b1.s23 * a1.se + b1.s32 * (double2)(-a1.sf, a1.sf);
            c2.sef += b2.s23 * a1.se + b2.s32 * (double2)(-a1.sf, a1.sf);
            c3.sef += b3.s23 * a1.se + b3.s32 * (double2)(-a1.sf, a1.sf);
            c4.sef += b4.s23 * a1.se + b4.s32 * (double2)(-a1.sf, a1.sf);
            c5.sef += b5.s23 * a1.se + b5.s32 * (double2)(-a1.sf, a1.sf);
            c6.sef += b6.s23 * a1.se + b6.s32 * (double2)(-a1.sf, a1.sf);
            c7.sef += b7.s23 * a1.se + b7.s32 * (double2)(-a1.sf, a1.sf);

            pA += lda;
            pB += 1;
            /* ---------------------- */
        }
        for (; k0 <= max(0, (int)M - 128); k0 += 2) {
            /* -- Tiles multiplier -- */
            b0 = pB[0];
            b1 = pB[(ldb >> 1)];
            b2 = pB[ldb];
            b3 = pB[mad24(3u, (ldb >> 1), 0u)];
            b4 = pB[(ldb << 1)];
            b5 = pB[mad24(5u, (ldb >> 1), 0u)];
            b6 = pB[mad24(6u, (ldb >> 1), 0u)];
            b7 = pB[mad24(7u, (ldb >> 1), 0u)];

            a0.s0123 = pA[0];
            a0.s4567 = pA[1];
            a0.s89ab = pA[2];
            a0.scdef = pA[3];
            a1.s0123 = pA[(lda >> 1)];
            a1.s4567 = pA[(lda >> 1) + 1];
            a1.s89ab = pA[(lda >> 1) + 2];
            a1.scdef = pA[(lda >> 1) + 3];

            c0.s01 += b0.s01 * a0.s0 + b0.s10 * (double2)(-a0.s1, a0.s1);
            c1.s01 += b1.s01 * a0.s0 + b1.s10 * (double2)(-a0.s1, a0.s1);
            c2.s01 += b2.s01 * a0.s0 + b2.s10 * (double2)(-a0.s1, a0.s1);
            c3.s01 += b3.s01 * a0.s0 + b3.s10 * (double2)(-a0.s1, a0.s1);
            c4.s01 += b4.s01 * a0.s0 + b4.s10 * (double2)(-a0.s1, a0.s1);
            c5.s01 += b5.s01 * a0.s0 + b5.s10 * (double2)(-a0.s1, a0.s1);
            c6.s01 += b6.s01 * a0.s0 + b6.s10 * (double2)(-a0.s1, a0.s1);
            c7.s01 += b7.s01 * a0.s0 + b7.s10 * (double2)(-a0.s1, a0.s1);
            c0.s23 += b0.s01 * a0.s2 + b0.s10 * (double2)(-a0.s3, a0.s3);
            c1.s23 += b1.s01 * a0.s2 + b1.s10 * (double2)(-a0.s3, a0.s3);
            c2.s23 += b2.s01 * a0.s2 + b2.s10 * (double2)(-a0.s3, a0.s3);
            c3.s23 += b3.s01 * a0.s2 + b3.s10 * (double2)(-a0.s3, a0.s3);
            c4.s23 += b4.s01 * a0.s2 + b4.s10 * (double2)(-a0.s3, a0.s3);
            c5.s23 += b5.s01 * a0.s2 + b5.s10 * (double2)(-a0.s3, a0.s3);
            c6.s23 += b6.s01 * a0.s2 + b6.s10 * (double2)(-a0.s3, a0.s3);
            c7.s23 += b7.s01 * a0.s2 + b7.s10 * (double2)(-a0.s3, a0.s3);
            c0.s45 += b0.s01 * a0.s4 + b0.s10 * (double2)(-a0.s5, a0.s5);
            c1.s45 += b1.s01 * a0.s4 + b1.s10 * (double2)(-a0.s5, a0.s5);
            c2.s45 += b2.s01 * a0.s4 + b2.s10 * (double2)(-a0.s5, a0.s5);
            c3.s45 += b3.s01 * a0.s4 + b3.s10 * (double2)(-a0.s5, a0.s5);
            c4.s45 += b4.s01 * a0.s4 + b4.s10 * (double2)(-a0.s5, a0.s5);
            c5.s45 += b5.s01 * a0.s4 + b5.s10 * (double2)(-a0.s5, a0.s5);
            c6.s45 += b6.s01 * a0.s4 + b6.s10 * (double2)(-a0.s5, a0.s5);
            c7.s45 += b7.s01 * a0.s4 + b7.s10 * (double2)(-a0.s5, a0.s5);
            c0.s67 += b0.s01 * a0.s6 + b0.s10 * (double2)(-a0.s7, a0.s7);
            c1.s67 += b1.s01 * a0.s6 + b1.s10 * (double2)(-a0.s7, a0.s7);
            c2.s67 += b2.s01 * a0.s6 + b2.s10 * (double2)(-a0.s7, a0.s7);
            c3.s67 += b3.s01 * a0.s6 + b3.s10 * (double2)(-a0.s7, a0.s7);
            c4.s67 += b4.s01 * a0.s6 + b4.s10 * (double2)(-a0.s7, a0.s7);
            c5.s67 += b5.s01 * a0.s6 + b5.s10 * (double2)(-a0.s7, a0.s7);
            c6.s67 += b6.s01 * a0.s6 + b6.s10 * (double2)(-a0.s7, a0.s7);
            c7.s67 += b7.s01 * a0.s6 + b7.s10 * (double2)(-a0.s7, a0.s7);
            c0.s89 += b0.s01 * a0.s8 + b0.s10 * (double2)(-a0.s9, a0.s9);
            c1.s89 += b1.s01 * a0.s8 + b1.s10 * (double2)(-a0.s9, a0.s9);
            c2.s89 += b2.s01 * a0.s8 + b2.s10 * (double2)(-a0.s9, a0.s9);
            c3.s89 += b3.s01 * a0.s8 + b3.s10 * (double2)(-a0.s9, a0.s9);
            c4.s89 += b4.s01 * a0.s8 + b4.s10 * (double2)(-a0.s9, a0.s9);
            c5.s89 += b5.s01 * a0.s8 + b5.s10 * (double2)(-a0.s9, a0.s9);
            c6.s89 += b6.s01 * a0.s8 + b6.s10 * (double2)(-a0.s9, a0.s9);
            c7.s89 += b7.s01 * a0.s8 + b7.s10 * (double2)(-a0.s9, a0.s9);
            c0.sab += b0.s01 * a0.sa + b0.s10 * (double2)(-a0.sb, a0.sb);
            c1.sab += b1.s01 * a0.sa + b1.s10 * (double2)(-a0.sb, a0.sb);
            c2.sab += b2.s01 * a0.sa + b2.s10 * (double2)(-a0.sb, a0.sb);
            c3.sab += b3.s01 * a0.sa + b3.s10 * (double2)(-a0.sb, a0.sb);
            c4.sab += b4.s01 * a0.sa + b4.s10 * (double2)(-a0.sb, a0.sb);
            c5.sab += b5.s01 * a0.sa + b5.s10 * (double2)(-a0.sb, a0.sb);
            c6.sab += b6.s01 * a0.sa + b6.s10 * (double2)(-a0.sb, a0.sb);
            c7.sab += b7.s01 * a0.sa + b7.s10 * (double2)(-a0.sb, a0.sb);
            c0.scd += b0.s01 * a0.sc + b0.s10 * (double2)(-a0.sd, a0.sd);
            c1.scd += b1.s01 * a0.sc + b1.s10 * (double2)(-a0.sd, a0.sd);
            c2.scd += b2.s01 * a0.sc + b2.s10 * (double2)(-a0.sd, a0.sd);
            c3.scd += b3.s01 * a0.sc + b3.s10 * (double2)(-a0.sd, a0.sd);
            c4.scd += b4.s01 * a0.sc + b4.s10 * (double2)(-a0.sd, a0.sd);
            c5.scd += b5.s01 * a0.sc + b5.s10 * (double2)(-a0.sd, a0.sd);
            c6.scd += b6.s01 * a0.sc + b6.s10 * (double2)(-a0.sd, a0.sd);
            c7.scd += b7.s01 * a0.sc + b7.s10 * (double2)(-a0.sd, a0.sd);
            c0.sef += b0.s01 * a0.se + b0.s10 * (double2)(-a0.sf, a0.sf);
            c1.sef += b1.s01 * a0.se + b1.s10 * (double2)(-a0.sf, a0.sf);
            c2.sef += b2.s01 * a0.se + b2.s10 * (double2)(-a0.sf, a0.sf);
            c3.sef += b3.s01 * a0.se + b3.s10 * (double2)(-a0.sf, a0.sf);
            c4.sef += b4.s01 * a0.se + b4.s10 * (double2)(-a0.sf, a0.sf);
            c5.sef += b5.s01 * a0.se + b5.s10 * (double2)(-a0.sf, a0.sf);
            c6.sef += b6.s01 * a0.se + b6.s10 * (double2)(-a0.sf, a0.sf);
            c7.sef += b7.s01 * a0.se + b7.s10 * (double2)(-a0.sf, a0.sf);

            c0.s01 += b0.s23 * a1.s0 + b0.s32 * (double2)(-a1.s1, a1.s1);
            c1.s01 += b1.s23 * a1.s0 + b1.s32 * (double2)(-a1.s1, a1.s1);
            c2.s01 += b2.s23 * a1.s0 + b2.s32 * (double2)(-a1.s1, a1.s1);
            c3.s01 += b3.s23 * a1.s0 + b3.s32 * (double2)(-a1.s1, a1.s1);
            c4.s01 += b4.s23 * a1.s0 + b4.s32 * (double2)(-a1.s1, a1.s1);
            c5.s01 += b5.s23 * a1.s0 + b5.s32 * (double2)(-a1.s1, a1.s1);
            c6.s01 += b6.s23 * a1.s0 + b6.s32 * (double2)(-a1.s1, a1.s1);
            c7.s01 += b7.s23 * a1.s0 + b7.s32 * (double2)(-a1.s1, a1.s1);
            c0.s23 += b0.s23 * a1.s2 + b0.s32 * (double2)(-a1.s3, a1.s3);
            c1.s23 += b1.s23 * a1.s2 + b1.s32 * (double2)(-a1.s3, a1.s3);
            c2.s23 += b2.s23 * a1.s2 + b2.s32 * (double2)(-a1.s3, a1.s3);
            c3.s23 += b3.s23 * a1.s2 + b3.s32 * (double2)(-a1.s3, a1.s3);
            c4.s23 += b4.s23 * a1.s2 + b4.s32 * (double2)(-a1.s3, a1.s3);
            c5.s23 += b5.s23 * a1.s2 + b5.s32 * (double2)(-a1.s3, a1.s3);
            c6.s23 += b6.s23 * a1.s2 + b6.s32 * (double2)(-a1.s3, a1.s3);
            c7.s23 += b7.s23 * a1.s2 + b7.s32 * (double2)(-a1.s3, a1.s3);
            c0.s45 += b0.s23 * a1.s4 + b0.s32 * (double2)(-a1.s5, a1.s5);
            c1.s45 += b1.s23 * a1.s4 + b1.s32 * (double2)(-a1.s5, a1.s5);
            c2.s45 += b2.s23 * a1.s4 + b2.s32 * (double2)(-a1.s5, a1.s5);
            c3.s45 += b3.s23 * a1.s4 + b3.s32 * (double2)(-a1.s5, a1.s5);
            c4.s45 += b4.s23 * a1.s4 + b4.s32 * (double2)(-a1.s5, a1.s5);
            c5.s45 += b5.s23 * a1.s4 + b5.s32 * (double2)(-a1.s5, a1.s5);
            c6.s45 += b6.s23 * a1.s4 + b6.s32 * (double2)(-a1.s5, a1.s5);
            c7.s45 += b7.s23 * a1.s4 + b7.s32 * (double2)(-a1.s5, a1.s5);
            c0.s67 += b0.s23 * a1.s6 + b0.s32 * (double2)(-a1.s7, a1.s7);
            c1.s67 += b1.s23 * a1.s6 + b1.s32 * (double2)(-a1.s7, a1.s7);
            c2.s67 += b2.s23 * a1.s6 + b2.s32 * (double2)(-a1.s7, a1.s7);
            c3.s67 += b3.s23 * a1.s6 + b3.s32 * (double2)(-a1.s7, a1.s7);
            c4.s67 += b4.s23 * a1.s6 + b4.s32 * (double2)(-a1.s7, a1.s7);
            c5.s67 += b5.s23 * a1.s6 + b5.s32 * (double2)(-a1.s7, a1.s7);
            c6.s67 += b6.s23 * a1.s6 + b6.s32 * (double2)(-a1.s7, a1.s7);
            c7.s67 += b7.s23 * a1.s6 + b7.s32 * (double2)(-a1.s7, a1.s7);
            c0.s89 += b0.s23 * a1.s8 + b0.s32 * (double2)(-a1.s9, a1.s9);
            c1.s89 += b1.s23 * a1.s8 + b1.s32 * (double2)(-a1.s9, a1.s9);
            c2.s89 += b2.s23 * a1.s8 + b2.s32 * (double2)(-a1.s9, a1.s9);
            c3.s89 += b3.s23 * a1.s8 + b3.s32 * (double2)(-a1.s9, a1.s9);
            c4.s89 += b4.s23 * a1.s8 + b4.s32 * (double2)(-a1.s9, a1.s9);
            c5.s89 += b5.s23 * a1.s8 + b5.s32 * (double2)(-a1.s9, a1.s9);
            c6.s89 += b6.s23 * a1.s8 + b6.s32 * (double2)(-a1.s9, a1.s9);
            c7.s89 += b7.s23 * a1.s8 + b7.s32 * (double2)(-a1.s9, a1.s9);
            c0.sab += b0.s23 * a1.sa + b0.s32 * (double2)(-a1.sb, a1.sb);
            c1.sab += b1.s23 * a1.sa + b1.s32 * (double2)(-a1.sb, a1.sb);
            c2.sab += b2.s23 * a1.sa + b2.s32 * (double2)(-a1.sb, a1.sb);
            c3.sab += b3.s23 * a1.sa + b3.s32 * (double2)(-a1.sb, a1.sb);
            c4.sab += b4.s23 * a1.sa + b4.s32 * (double2)(-a1.sb, a1.sb);
            c5.sab += b5.s23 * a1.sa + b5.s32 * (double2)(-a1.sb, a1.sb);
            c6.sab += b6.s23 * a1.sa + b6.s32 * (double2)(-a1.sb, a1.sb);
            c7.sab += b7.s23 * a1.sa + b7.s32 * (double2)(-a1.sb, a1.sb);
            c0.scd += b0.s23 * a1.sc + b0.s32 * (double2)(-a1.sd, a1.sd);
            c1.scd += b1.s23 * a1.sc + b1.s32 * (double2)(-a1.sd, a1.sd);
            c2.scd += b2.s23 * a1.sc + b2.s32 * (double2)(-a1.sd, a1.sd);
            c3.scd += b3.s23 * a1.sc + b3.s32 * (double2)(-a1.sd, a1.sd);
            c4.scd += b4.s23 * a1.sc + b4.s32 * (double2)(-a1.sd, a1.sd);
            c5.scd += b5.s23 * a1.sc + b5.s32 * (double2)(-a1.sd, a1.sd);
            c6.scd += b6.s23 * a1.sc + b6.s32 * (double2)(-a1.sd, a1.sd);
            c7.scd += b7.s23 * a1.sc + b7.s32 * (double2)(-a1.sd, a1.sd);
            c0.sef += b0.s23 * a1.se + b0.s32 * (double2)(-a1.sf, a1.sf);
            c1.sef += b1.s23 * a1.se + b1.s32 * (double2)(-a1.sf, a1.sf);
            c2.sef += b2.s23 * a1.se + b2.s32 * (double2)(-a1.sf, a1.sf);
            c3.sef += b3.s23 * a1.se + b3.s32 * (double2)(-a1.sf, a1.sf);
            c4.sef += b4.s23 * a1.se + b4.s32 * (double2)(-a1.sf, a1.sf);
            c5.sef += b5.s23 * a1.se + b5.s32 * (double2)(-a1.sf, a1.sf);
            c6.sef += b6.s23 * a1.se + b6.s32 * (double2)(-a1.sf, a1.sf);
            c7.sef += b7.s23 * a1.se + b7.s32 * (double2)(-a1.sf, a1.sf);

            pA += lda;
            pB += 1;
            /* ---------------------- */
        }
        for (; k0 < M; k0 += 2) {
            coord.z = k0;
            /* -- Tiles multiplier -- */
            b0 = pB[0];
            b1 = pB[(ldb >> 1)];
            b2 = pB[ldb];
            b3 = pB[mad24(3u, (ldb >> 1), 0u)];
            b4 = pB[(ldb << 1)];
            b5 = pB[mad24(5u, (ldb >> 1), 0u)];
            b6 = pB[mad24(6u, (ldb >> 1), 0u)];
            b7 = pB[mad24(7u, (ldb >> 1), 0u)];


            a0.s0123 = pA[0];
            a0.s4567 = pA[1];
            a0.s89ab = pA[2];
            a0.scdef = pA[3];
            a1.s0123 = pA[(lda >> 1)];
            a1.s4567 = pA[(lda >> 1) + 1];
            a1.s89ab = pA[(lda >> 1) + 2];
            a1.scdef = pA[(lda >> 1) + 3];

            // post fetch A
            {
                uint zy = coord.z;
                a0.s01 = zy < coord.y ? 0 : a0.s01;
                a0.s01 = zy == coord.y ? (double2)(1, 0) : a0.s01;
                a0.s23 = zy < coord.y + 1 ? 0 : a0.s23;
                a0.s23 = zy == coord.y + 1 ? (double2)(1, 0) : a0.s23;
                a0.s45 = zy < coord.y + 2 ? 0 : a0.s45;
                a0.s45 = zy == coord.y + 2 ? (double2)(1, 0) : a0.s45;
                a0.s67 = zy < coord.y + 3 ? 0 : a0.s67;
                a0.s67 = zy == coord.y + 3 ? (double2)(1, 0) : a0.s67;
                a0.s89 = zy < coord.y + 4 ? 0 : a0.s89;
                a0.s89 = zy == coord.y + 4 ? (double2)(1, 0) : a0.s89;
                a0.sab = zy < coord.y + 5 ? 0 : a0.sab;
                a0.sab = zy == coord.y + 5 ? (double2)(1, 0) : a0.sab;
                a0.scd = zy < coord.y + 6 ? 0 : a0.scd;
                a0.scd = zy == coord.y + 6 ? (double2)(1, 0) : a0.scd;
                a0.sef = zy < coord.y + 7 ? 0 : a0.sef;
                a0.sef = zy == coord.y + 7 ? (double2)(1, 0) : a0.sef;
                zy++;
                a1.s01 = zy < coord.y ? 0 : a1.s01;
                a1.s01 = zy == coord.y ? (double2)(1, 0) : a1.s01;
                a1.s23 = zy < coord.y + 1 ? 0 : a1.s23;
                a1.s23 = zy == coord.y + 1 ? (double2)(1, 0) : a1.s23;
                a1.s45 = zy < coord.y + 2 ? 0 : a1.s45;
                a1.s45 = zy == coord.y + 2 ? (double2)(1, 0) : a1.s45;
                a1.s67 = zy < coord.y + 3 ? 0 : a1.s67;
                a1.s67 = zy == coord.y + 3 ? (double2)(1, 0) : a1.s67;
                a1.s89 = zy < coord.y + 4 ? 0 : a1.s89;
                a1.s89 = zy == coord.y + 4 ? (double2)(1, 0) : a1.s89;
                a1.sab = zy < coord.y + 5 ? 0 : a1.sab;
                a1.sab = zy == coord.y + 5 ? (double2)(1, 0) : a1.sab;
                a1.scd = zy < coord.y + 6 ? 0 : a1.scd;
                a1.scd = zy == coord.y + 6 ? (double2)(1, 0) : a1.scd;
                a1.sef = zy < coord.y + 7 ? 0 : a1.sef;
                a1.sef = zy == coord.y + 7 ? (double2)(1, 0) : a1.sef;
            }

            c0.s01 += b0.s01 * a0.s0 + b0.s10 * (double2)(-a0.s1, a0.s1);
            c1.s01 += b1.s01 * a0.s0 + b1.s10 * (double2)(-a0.s1, a0.s1);
            c2.s01 += b2.s01 * a0.s0 + b2.s10 * (double2)(-a0.s1, a0.s1);
            c3.s01 += b3.s01 * a0.s0 + b3.s10 * (double2)(-a0.s1, a0.s1);
            c4.s01 += b4.s01 * a0.s0 + b4.s10 * (double2)(-a0.s1, a0.s1);
            c5.s01 += b5.s01 * a0.s0 + b5.s10 * (double2)(-a0.s1, a0.s1);
            c6.s01 += b6.s01 * a0.s0 + b6.s10 * (double2)(-a0.s1, a0.s1);
            c7.s01 += b7.s01 * a0.s0 + b7.s10 * (double2)(-a0.s1, a0.s1);
            c0.s23 += b0.s01 * a0.s2 + b0.s10 * (double2)(-a0.s3, a0.s3);
            c1.s23 += b1.s01 * a0.s2 + b1.s10 * (double2)(-a0.s3, a0.s3);
            c2.s23 += b2.s01 * a0.s2 + b2.s10 * (double2)(-a0.s3, a0.s3);
            c3.s23 += b3.s01 * a0.s2 + b3.s10 * (double2)(-a0.s3, a0.s3);
            c4.s23 += b4.s01 * a0.s2 + b4.s10 * (double2)(-a0.s3, a0.s3);
            c5.s23 += b5.s01 * a0.s2 + b5.s10 * (double2)(-a0.s3, a0.s3);
            c6.s23 += b6.s01 * a0.s2 + b6.s10 * (double2)(-a0.s3, a0.s3);
            c7.s23 += b7.s01 * a0.s2 + b7.s10 * (double2)(-a0.s3, a0.s3);
            c0.s45 += b0.s01 * a0.s4 + b0.s10 * (double2)(-a0.s5, a0.s5);
            c1.s45 += b1.s01 * a0.s4 + b1.s10 * (double2)(-a0.s5, a0.s5);
            c2.s45 += b2.s01 * a0.s4 + b2.s10 * (double2)(-a0.s5, a0.s5);
            c3.s45 += b3.s01 * a0.s4 + b3.s10 * (double2)(-a0.s5, a0.s5);
            c4.s45 += b4.s01 * a0.s4 + b4.s10 * (double2)(-a0.s5, a0.s5);
            c5.s45 += b5.s01 * a0.s4 + b5.s10 * (double2)(-a0.s5, a0.s5);
            c6.s45 += b6.s01 * a0.s4 + b6.s10 * (double2)(-a0.s5, a0.s5);
            c7.s45 += b7.s01 * a0.s4 + b7.s10 * (double2)(-a0.s5, a0.s5);
            c0.s67 += b0.s01 * a0.s6 + b0.s10 * (double2)(-a0.s7, a0.s7);
            c1.s67 += b1.s01 * a0.s6 + b1.s10 * (double2)(-a0.s7, a0.s7);
            c2.s67 += b2.s01 * a0.s6 + b2.s10 * (double2)(-a0.s7, a0.s7);
            c3.s67 += b3.s01 * a0.s6 + b3.s10 * (double2)(-a0.s7, a0.s7);
            c4.s67 += b4.s01 * a0.s6 + b4.s10 * (double2)(-a0.s7, a0.s7);
            c5.s67 += b5.s01 * a0.s6 + b5.s10 * (double2)(-a0.s7, a0.s7);
            c6.s67 += b6.s01 * a0.s6 + b6.s10 * (double2)(-a0.s7, a0.s7);
            c7.s67 += b7.s01 * a0.s6 + b7.s10 * (double2)(-a0.s7, a0.s7);
            c0.s89 += b0.s01 * a0.s8 + b0.s10 * (double2)(-a0.s9, a0.s9);
            c1.s89 += b1.s01 * a0.s8 + b1.s10 * (double2)(-a0.s9, a0.s9);
            c2.s89 += b2.s01 * a0.s8 + b2.s10 * (double2)(-a0.s9, a0.s9);
            c3.s89 += b3.s01 * a0.s8 + b3.s10 * (double2)(-a0.s9, a0.s9);
            c4.s89 += b4.s01 * a0.s8 + b4.s10 * (double2)(-a0.s9, a0.s9);
            c5.s89 += b5.s01 * a0.s8 + b5.s10 * (double2)(-a0.s9, a0.s9);
            c6.s89 += b6.s01 * a0.s8 + b6.s10 * (double2)(-a0.s9, a0.s9);
            c7.s89 += b7.s01 * a0.s8 + b7.s10 * (double2)(-a0.s9, a0.s9);
            c0.sab += b0.s01 * a0.sa + b0.s10 * (double2)(-a0.sb, a0.sb);
            c1.sab += b1.s01 * a0.sa + b1.s10 * (double2)(-a0.sb, a0.sb);
            c2.sab += b2.s01 * a0.sa + b2.s10 * (double2)(-a0.sb, a0.sb);
            c3.sab += b3.s01 * a0.sa + b3.s10 * (double2)(-a0.sb, a0.sb);
            c4.sab += b4.s01 * a0.sa + b4.s10 * (double2)(-a0.sb, a0.sb);
            c5.sab += b5.s01 * a0.sa + b5.s10 * (double2)(-a0.sb, a0.sb);
            c6.sab += b6.s01 * a0.sa + b6.s10 * (double2)(-a0.sb, a0.sb);
            c7.sab += b7.s01 * a0.sa + b7.s10 * (double2)(-a0.sb, a0.sb);
            c0.scd += b0.s01 * a0.sc + b0.s10 * (double2)(-a0.sd, a0.sd);
            c1.scd += b1.s01 * a0.sc + b1.s10 * (double2)(-a0.sd, a0.sd);
            c2.scd += b2.s01 * a0.sc + b2.s10 * (double2)(-a0.sd, a0.sd);
            c3.scd += b3.s01 * a0.sc + b3.s10 * (double2)(-a0.sd, a0.sd);
            c4.scd += b4.s01 * a0.sc + b4.s10 * (double2)(-a0.sd, a0.sd);
            c5.scd += b5.s01 * a0.sc + b5.s10 * (double2)(-a0.sd, a0.sd);
            c6.scd += b6.s01 * a0.sc + b6.s10 * (double2)(-a0.sd, a0.sd);
            c7.scd += b7.s01 * a0.sc + b7.s10 * (double2)(-a0.sd, a0.sd);
            c0.sef += b0.s01 * a0.se + b0.s10 * (double2)(-a0.sf, a0.sf);
            c1.sef += b1.s01 * a0.se + b1.s10 * (double2)(-a0.sf, a0.sf);
            c2.sef += b2.s01 * a0.se + b2.s10 * (double2)(-a0.sf, a0.sf);
            c3.sef += b3.s01 * a0.se + b3.s10 * (double2)(-a0.sf, a0.sf);
            c4.sef += b4.s01 * a0.se + b4.s10 * (double2)(-a0.sf, a0.sf);
            c5.sef += b5.s01 * a0.se + b5.s10 * (double2)(-a0.sf, a0.sf);
            c6.sef += b6.s01 * a0.se + b6.s10 * (double2)(-a0.sf, a0.sf);
            c7.sef += b7.s01 * a0.se + b7.s10 * (double2)(-a0.sf, a0.sf);

            c0.s01 += b0.s23 * a1.s0 + b0.s32 * (double2)(-a1.s1, a1.s1);
            c1.s01 += b1.s23 * a1.s0 + b1.s32 * (double2)(-a1.s1, a1.s1);
            c2.s01 += b2.s23 * a1.s0 + b2.s32 * (double2)(-a1.s1, a1.s1);
            c3.s01 += b3.s23 * a1.s0 + b3.s32 * (double2)(-a1.s1, a1.s1);
            c4.s01 += b4.s23 * a1.s0 + b4.s32 * (double2)(-a1.s1, a1.s1);
            c5.s01 += b5.s23 * a1.s0 + b5.s32 * (double2)(-a1.s1, a1.s1);
            c6.s01 += b6.s23 * a1.s0 + b6.s32 * (double2)(-a1.s1, a1.s1);
            c7.s01 += b7.s23 * a1.s0 + b7.s32 * (double2)(-a1.s1, a1.s1);
            c0.s23 += b0.s23 * a1.s2 + b0.s32 * (double2)(-a1.s3, a1.s3);
            c1.s23 += b1.s23 * a1.s2 + b1.s32 * (double2)(-a1.s3, a1.s3);
            c2.s23 += b2.s23 * a1.s2 + b2.s32 * (double2)(-a1.s3, a1.s3);
            c3.s23 += b3.s23 * a1.s2 + b3.s32 * (double2)(-a1.s3, a1.s3);
            c4.s23 += b4.s23 * a1.s2 + b4.s32 * (double2)(-a1.s3, a1.s3);
            c5.s23 += b5.s23 * a1.s2 + b5.s32 * (double2)(-a1.s3, a1.s3);
            c6.s23 += b6.s23 * a1.s2 + b6.s32 * (double2)(-a1.s3, a1.s3);
            c7.s23 += b7.s23 * a1.s2 + b7.s32 * (double2)(-a1.s3, a1.s3);
            c0.s45 += b0.s23 * a1.s4 + b0.s32 * (double2)(-a1.s5, a1.s5);
            c1.s45 += b1.s23 * a1.s4 + b1.s32 * (double2)(-a1.s5, a1.s5);
            c2.s45 += b2.s23 * a1.s4 + b2.s32 * (double2)(-a1.s5, a1.s5);
            c3.s45 += b3.s23 * a1.s4 + b3.s32 * (double2)(-a1.s5, a1.s5);
            c4.s45 += b4.s23 * a1.s4 + b4.s32 * (double2)(-a1.s5, a1.s5);
            c5.s45 += b5.s23 * a1.s4 + b5.s32 * (double2)(-a1.s5, a1.s5);
            c6.s45 += b6.s23 * a1.s4 + b6.s32 * (double2)(-a1.s5, a1.s5);
            c7.s45 += b7.s23 * a1.s4 + b7.s32 * (double2)(-a1.s5, a1.s5);
            c0.s67 += b0.s23 * a1.s6 + b0.s32 * (double2)(-a1.s7, a1.s7);
            c1.s67 += b1.s23 * a1.s6 + b1.s32 * (double2)(-a1.s7, a1.s7);
            c2.s67 += b2.s23 * a1.s6 + b2.s32 * (double2)(-a1.s7, a1.s7);
            c3.s67 += b3.s23 * a1.s6 + b3.s32 * (double2)(-a1.s7, a1.s7);
            c4.s67 += b4.s23 * a1.s6 + b4.s32 * (double2)(-a1.s7, a1.s7);
            c5.s67 += b5.s23 * a1.s6 + b5.s32 * (double2)(-a1.s7, a1.s7);
            c6.s67 += b6.s23 * a1.s6 + b6.s32 * (double2)(-a1.s7, a1.s7);
            c7.s67 += b7.s23 * a1.s6 + b7.s32 * (double2)(-a1.s7, a1.s7);
            c0.s89 += b0.s23 * a1.s8 + b0.s32 * (double2)(-a1.s9, a1.s9);
            c1.s89 += b1.s23 * a1.s8 + b1.s32 * (double2)(-a1.s9, a1.s9);
            c2.s89 += b2.s23 * a1.s8 + b2.s32 * (double2)(-a1.s9, a1.s9);
            c3.s89 += b3.s23 * a1.s8 + b3.s32 * (double2)(-a1.s9, a1.s9);
            c4.s89 += b4.s23 * a1.s8 + b4.s32 * (double2)(-a1.s9, a1.s9);
            c5.s89 += b5.s23 * a1.s8 + b5.s32 * (double2)(-a1.s9, a1.s9);
            c6.s89 += b6.s23 * a1.s8 + b6.s32 * (double2)(-a1.s9, a1.s9);
            c7.s89 += b7.s23 * a1.s8 + b7.s32 * (double2)(-a1.s9, a1.s9);
            c0.sab += b0.s23 * a1.sa + b0.s32 * (double2)(-a1.sb, a1.sb);
            c1.sab += b1.s23 * a1.sa + b1.s32 * (double2)(-a1.sb, a1.sb);
            c2.sab += b2.s23 * a1.sa + b2.s32 * (double2)(-a1.sb, a1.sb);
            c3.sab += b3.s23 * a1.sa + b3.s32 * (double2)(-a1.sb, a1.sb);
            c4.sab += b4.s23 * a1.sa + b4.s32 * (double2)(-a1.sb, a1.sb);
            c5.sab += b5.s23 * a1.sa + b5.s32 * (double2)(-a1.sb, a1.sb);
            c6.sab += b6.s23 * a1.sa + b6.s32 * (double2)(-a1.sb, a1.sb);
            c7.sab += b7.s23 * a1.sa + b7.s32 * (double2)(-a1.sb, a1.sb);
            c0.scd += b0.s23 * a1.sc + b0.s32 * (double2)(-a1.sd, a1.sd);
            c1.scd += b1.s23 * a1.sc + b1.s32 * (double2)(-a1.sd, a1.sd);
            c2.scd += b2.s23 * a1.sc + b2.s32 * (double2)(-a1.sd, a1.sd);
            c3.scd += b3.s23 * a1.sc + b3.s32 * (double2)(-a1.sd, a1.sd);
            c4.scd += b4.s23 * a1.sc + b4.s32 * (double2)(-a1.sd, a1.sd);
            c5.scd += b5.s23 * a1.sc + b5.s32 * (double2)(-a1.sd, a1.sd);
            c6.scd += b6.s23 * a1.sc + b6.s32 * (double2)(-a1.sd, a1.sd);
            c7.scd += b7.s23 * a1.sc + b7.s32 * (double2)(-a1.sd, a1.sd);
            c0.sef += b0.s23 * a1.se + b0.s32 * (double2)(-a1.sf, a1.sf);
            c1.sef += b1.s23 * a1.se + b1.s32 * (double2)(-a1.sf, a1.sf);
            c2.sef += b2.s23 * a1.se + b2.s32 * (double2)(-a1.sf, a1.sf);
            c3.sef += b3.s23 * a1.se + b3.s32 * (double2)(-a1.sf, a1.sf);
            c4.sef += b4.s23 * a1.se + b4.s32 * (double2)(-a1.sf, a1.sf);
            c5.sef += b5.s23 * a1.se + b5.s32 * (double2)(-a1.sf, a1.sf);
            c6.sef += b6.s23 * a1.se + b6.s32 * (double2)(-a1.sf, a1.sf);
            c7.sef += b7.s23 * a1.se + b7.s32 * (double2)(-a1.sf, a1.sf);

            pA += lda;
            pB += 1;
            /* ---------------------- */
        }
        barrier(CLK_GLOBAL_MEM_FENCE);

        GPtr uC;
        double2 alphaR = (double2)(alpha.x);
        double2 alphaI = (double2)(-alpha.y, alpha.y);

        uC.d2v = C + coord.x * ldb + coord.y;

        __global double4 *pC = uC.d2v;

        double16 tempC0, tempC1, tempC2, tempC3;

        tempC0.s01 = alpha * c0.s0 + alpha.s10 * (double2)(-c0.s1, c0.s1);
        tempC0.s23 = alpha * c0.s2 + alpha.s10 * (double2)(-c0.s3, c0.s3);
        tempC0.s45 = alpha * c0.s4 + alpha.s10 * (double2)(-c0.s5, c0.s5);
        tempC0.s67 = alpha * c0.s6 + alpha.s10 * (double2)(-c0.s7, c0.s7);
        tempC0.s89 = alpha * c0.s8 + alpha.s10 * (double2)(-c0.s9, c0.s9);
        tempC0.sab = alpha * c0.sa + alpha.s10 * (double2)(-c0.sb, c0.sb);
        tempC0.scd = alpha * c0.sc + alpha.s10 * (double2)(-c0.sd, c0.sd);
        tempC0.sef = alpha * c0.se + alpha.s10 * (double2)(-c0.sf, c0.sf);
        tempC1.s01 = alpha * c1.s0 + alpha.s10 * (double2)(-c1.s1, c1.s1);
        tempC1.s23 = alpha * c1.s2 + alpha.s10 * (double2)(-c1.s3, c1.s3);
        tempC1.s45 = alpha * c1.s4 + alpha.s10 * (double2)(-c1.s5, c1.s5);
        tempC1.s67 = alpha * c1.s6 + alpha.s10 * (double2)(-c1.s7, c1.s7);
        tempC1.s89 = alpha * c1.s8 + alpha.s10 * (double2)(-c1.s9, c1.s9);
        tempC1.sab = alpha * c1.sa + alpha.s10 * (double2)(-c1.sb, c1.sb);
        tempC1.scd = alpha * c1.sc + alpha.s10 * (double2)(-c1.sd, c1.sd);
        tempC1.sef = alpha * c1.se + alpha.s10 * (double2)(-c1.sf, c1.sf);
        tempC2.s01 = alpha * c2.s0 + alpha.s10 * (double2)(-c2.s1, c2.s1);
        tempC2.s23 = alpha * c2.s2 + alpha.s10 * (double2)(-c2.s3, c2.s3);
        tempC2.s45 = alpha * c2.s4 + alpha.s10 * (double2)(-c2.s5, c2.s5);
        tempC2.s67 = alpha * c2.s6 + alpha.s10 * (double2)(-c2.s7, c2.s7);
        tempC2.s89 = alpha * c2.s8 + alpha.s10 * (double2)(-c2.s9, c2.s9);
        tempC2.sab = alpha * c2.sa + alpha.s10 * (double2)(-c2.sb, c2.sb);
        tempC2.scd = alpha * c2.sc + alpha.s10 * (double2)(-c2.sd, c2.sd);
        tempC2.sef = alpha * c2.se + alpha.s10 * (double2)(-c2.sf, c2.sf);
        tempC3.s01 = alpha * c3.s0 + alpha.s10 * (double2)(-c3.s1, c3.s1);
        tempC3.s23 = alpha * c3.s2 + alpha.s10 * (double2)(-c3.s3, c3.s3);
        tempC3.s45 = alpha * c3.s4 + alpha.s10 * (double2)(-c3.s5, c3.s5);
        tempC3.s67 = alpha * c3.s6 + alpha.s10 * (double2)(-c3.s7, c3.s7);
        tempC3.s89 = alpha * c3.s8 + alpha.s10 * (double2)(-c3.s9, c3.s9);
        tempC3.sab = alpha * c3.sa + alpha.s10 * (double2)(-c3.sb, c3.sb);
        tempC3.scd = alpha * c3.sc + alpha.s10 * (double2)(-c3.sd, c3.sd);
        tempC3.sef = alpha * c3.se + alpha.s10 * (double2)(-c3.sf, c3.sf);
        pC[0] = tempC0.s0123;
        pC[1] = tempC0.s4567;
        pC[2] = tempC0.s89ab;
        pC[3] = tempC0.scdef;
        pC[(ldb >> 1)] = tempC1.s0123;
        pC[(ldb >> 1) + 1] = tempC1.s4567;
        pC[(ldb >> 1) + 2] = tempC1.s89ab;
        pC[(ldb >> 1) + 3] = tempC1.scdef;
        pC[ldb] = tempC2.s0123;
        pC[ldb + 1] = tempC2.s4567;
        pC[ldb + 2] = tempC2.s89ab;
        pC[ldb + 3] = tempC2.scdef;
        pC[mad24(3u, (ldb >> 1), 0u)] = tempC3.s0123;
        pC[mad24(3u, (ldb >> 1), 1u)] = tempC3.s4567;
        pC[mad24(3u, (ldb >> 1), 2u)] = tempC3.s89ab;
        pC[mad24(3u, (ldb >> 1), 3u)] = tempC3.scdef;

        tempC0.s01 = alpha * c4.s0 + alpha.s10 * (double2)(-c4.s1, c4.s1);
        tempC0.s23 = alpha * c4.s2 + alpha.s10 * (double2)(-c4.s3, c4.s3);
        tempC0.s45 = alpha * c4.s4 + alpha.s10 * (double2)(-c4.s5, c4.s5);
        tempC0.s67 = alpha * c4.s6 + alpha.s10 * (double2)(-c4.s7, c4.s7);
        tempC0.s89 = alpha * c4.s8 + alpha.s10 * (double2)(-c4.s9, c4.s9);
        tempC0.sab = alpha * c4.sa + alpha.s10 * (double2)(-c4.sb, c4.sb);
        tempC0.scd = alpha * c4.sc + alpha.s10 * (double2)(-c4.sd, c4.sd);
        tempC0.sef = alpha * c4.se + alpha.s10 * (double2)(-c4.sf, c4.sf);
        tempC1.s01 = alpha * c5.s0 + alpha.s10 * (double2)(-c5.s1, c5.s1);
        tempC1.s23 = alpha * c5.s2 + alpha.s10 * (double2)(-c5.s3, c5.s3);
        tempC1.s45 = alpha * c5.s4 + alpha.s10 * (double2)(-c5.s5, c5.s5);
        tempC1.s67 = alpha * c5.s6 + alpha.s10 * (double2)(-c5.s7, c5.s7);
        tempC1.s89 = alpha * c5.s8 + alpha.s10 * (double2)(-c5.s9, c5.s9);
        tempC1.sab = alpha * c5.sa + alpha.s10 * (double2)(-c5.sb, c5.sb);
        tempC1.scd = alpha * c5.sc + alpha.s10 * (double2)(-c5.sd, c5.sd);
        tempC1.sef = alpha * c5.se + alpha.s10 * (double2)(-c5.sf, c5.sf);
        tempC2.s01 = alpha * c6.s0 + alpha.s10 * (double2)(-c6.s1, c6.s1);
        tempC2.s23 = alpha * c6.s2 + alpha.s10 * (double2)(-c6.s3, c6.s3);
        tempC2.s45 = alpha * c6.s4 + alpha.s10 * (double2)(-c6.s5, c6.s5);
        tempC2.s67 = alpha * c6.s6 + alpha.s10 * (double2)(-c6.s7, c6.s7);
        tempC2.s89 = alpha * c6.s8 + alpha.s10 * (double2)(-c6.s9, c6.s9);
        tempC2.sab = alpha * c6.sa + alpha.s10 * (double2)(-c6.sb, c6.sb);
        tempC2.scd = alpha * c6.sc + alpha.s10 * (double2)(-c6.sd, c6.sd);
        tempC2.sef = alpha * c6.se + alpha.s10 * (double2)(-c6.sf, c6.sf);
        tempC3.s01 = alpha * c7.s0 + alpha.s10 * (double2)(-c7.s1, c7.s1);
        tempC3.s23 = alpha * c7.s2 + alpha.s10 * (double2)(-c7.s3, c7.s3);
        tempC3.s45 = alpha * c7.s4 + alpha.s10 * (double2)(-c7.s5, c7.s5);
        tempC3.s67 = alpha * c7.s6 + alpha.s10 * (double2)(-c7.s7, c7.s7);
        tempC3.s89 = alpha * c7.s8 + alpha.s10 * (double2)(-c7.s9, c7.s9);
        tempC3.sab = alpha * c7.sa + alpha.s10 * (double2)(-c7.sb, c7.sb);
        tempC3.scd = alpha * c7.sc + alpha.s10 * (double2)(-c7.sd, c7.sd);
        tempC3.sef = alpha * c7.se + alpha.s10 * (double2)(-c7.sf, c7.sf);
        pC[(ldb << 1)] = tempC0.s0123;
        pC[mad24(2u, ldb, 1u)] = tempC0.s4567;
        pC[mad24(2u, ldb, 2u)] = tempC0.s89ab;
        pC[mad24(2u, ldb, 3u)] = tempC0.scdef;
        pC[mad24(5u, (ldb >> 1), 0u)] = tempC1.s0123;
        pC[mad24(5u, (ldb >> 1), 1u)] = tempC1.s4567;
        pC[mad24(5u, (ldb >> 1), 2u)] = tempC1.s89ab;
        pC[mad24(5u, (ldb >> 1), 3u)] = tempC1.scdef;
        pC[mad24(6u, (ldb >> 1), 0u)] = tempC2.s0123;
        pC[mad24(6u, (ldb >> 1), 1u)] = tempC2.s4567;
        pC[mad24(6u, (ldb >> 1), 2u)] = tempC2.s89ab;
        pC[mad24(6u, (ldb >> 1), 3u)] = tempC2.scdef;
        pC[mad24(7u, (ldb >> 1), 0u)] = tempC3.s0123;
        pC[mad24(7u, (ldb >> 1), 1u)] = tempC3.s4567;
        pC[mad24(7u, (ldb >> 1), 2u)] = tempC3.s89ab;
        pC[mad24(7u, (ldb >> 1), 3u)] = tempC3.scdef;
        currM += 128;
    }
}


发表于 2018-5-25 13:45:35 | 显示全部楼层
有源码?想看下感兴趣的几个数值算法。
 楼主| 发表于 2018-5-25 15:17:03 | 显示全部楼层
同志您好

CLBLAS是开源的,里面有许多的矩阵、向量相关的算法,按照网上的步骤用cmake在windows上可以编译出需要的库来。

里面自带的测试过不了啊,ztrmm算出来的结果不知道怎么算的,https://www.math.utah.edu/software/lapack/lapack-blas/ztrmm.html
这里有介绍,但是我带进数去还是得到的跟想的不一样,帮忙看下啊!
发表于 2018-6-1 14:01:33 | 显示全部楼层
本坛几乎没有对数值类的感兴趣的人,我也只是面对几个工程问题才对此关注。话再说开来,国内几乎没有论坛讨论数学问题。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

QQ|关于我们|小黑屋|Archiver|手机版|中文第一计算机图形学社区OpenGPU

GMT+8, 2018-10-24 01:17 , Processed in 0.080148 second(s), 19 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表