Function xAx

jni/stasm/classicdesc.cpp:211–229

// mat is symmetric to roughly halve the number of operations.

static double xAx(
    const VEC& x,   // in
    const MAT& mat) //in: must be symmetric
{
    const int n = NSIZE(x);
    CV_Assert(mat.rows == n && mat.cols == n && x.isContinuous());
    const double* px = Buf(x);
    double diagsum = 0, sum = 0;
    int i = n;
    while (i--)
    {
        const double xi = px[i];
        const double* const rowbuf = mat.ptr<double>(i);
        diagsum += rowbuf[i] * SQ(xi);    // sum diag elements
        for (int j = i+1; j < n; j++)     // sum upper right triangle elements
            sum += rowbuf[j] * xi * px[j];
    }
    return diagsum + 2 * sum; // "2 *" to include lower left triangle elements
}

// If OpenMP is enabled, multiple instances of this function will be called
// concurrently (each call will have a different value of x and y). Thus this