Method mv_32

src/ialglib.cpp:44–90

This subroutine calculates fast 32x32 real matrix-vector product: y := beta*y + alpha*A*x using either generic C code or native optimizations (if available) IMPORTANT: * A must be stored in row-major order, stride is alglib_r_block, aligned on alglib_simd_alignment boundary * X must be aligned on alglib_simd_alignment boundary * Y may be non-aligned *************************************

* Y may be non-aligned
********************************************************************/
void ialglib::mv_32(const double *a, const double *x, double *y, int stride, double alpha, double beta)
{
    int i, k;
    const double *pa0, *pa1, *pb;

    pa0 = a;
    pa1 = a+alglib_r_block;
    pb = x;
    for(i=0; i<16; i++)
    {
        double v0 = 0, v1 = 0;
        for(k=0; k<4; k++)
        {
            v0 += pa0[0]*pb[0];
            v1 += pa1[0]*pb[0];
            v0 += pa0[1]*pb[1];
            v1 += pa1[1]*pb[1];
            v0 += pa0[2]*pb[2];
            v1 += pa1[2]*pb[2];
            v0 += pa0[3]*pb[3];
            v1 += pa1[3]*pb[3];
            v0 += pa0[4]*pb[4];
            v1 += pa1[4]*pb[4];
            v0 += pa0[5]*pb[5];
            v1 += pa1[5]*pb[5];
            v0 += pa0[6]*pb[6];
            v1 += pa1[6]*pb[6];
            v0 += pa0[7]*pb[7];
            v1 += pa1[7]*pb[7];
            pa0 += 8;
            pa1 += 8;
            pb  += 8;
        }
        y[0] = beta*y[0]+alpha*v0;
        y[stride] = beta*y[stride]+alpha*v1;

        //
        // now we've processed rows I and I+1,
        // pa0 and pa1 are pointing to rows I+1 and I+2.
        // move to I+2 and I+3.
        //
        pa0 += alglib_r_block;
        pa1 += alglib_r_block;
        pb = x;
        y+=2*stride;
    }
}


/********************************************************************