matrix_kernels.h

Go to the documentation of this file.

#ifndef TBCI_MATRIX_KERNELS_H
#define TBCI_MATRIX_KERNELS_H

template <typename T> class   Matrix;
template <typename T> class  TMatrix;
template <typename T> class TSMatrix;
template <typename T> class   Vector;
template <typename T> class  TVector;
template <typename T> class TSVector;

//#include "vec_kern_unr_pref.h"

#ifdef OLD_MAT_MAT_MULT
template <typename T>
/*inline*/ void do_mat_mat_mult (const unsigned start, const unsigned end,
                                 TMatrix<T> *res, const Matrix<T> *a,
                                 const Matrix<T> *b)
{
        const unsigned bc = b->columns(); 
        const unsigned ac = a->columns();
        for (unsigned i=start; i<end; ++i) {
                for (unsigned j=0; j<bc; ++j) {
                        PREFETCH_R(a->getcolptr(i), 2);
                        PREFETCH_W(&res->setval(i,j), 2);
                        ALIGN2 (register T tmp (a->get(i,0) * b->get(0,j)), MIN_ALIGN2);
                        for (register unsigned l=1; l<ac; ++l)
                                tmp += a->get(i,l) * b->get(l,j);
                        res->set(tmp,i,j);
                }
        }
}
#else


/* Following Dowd/Severance, ch. 8: Friendlier mem access pattern */
/* Should be faster in general. However, on iPII, it's faster (+)
 * or slower (-), depending on the data type T:
 * int: -, float: +, double: -, long long double: +, cplx<double>: +
 * AXP: Generally faster
 * (KG, 99/07/28) */
template <typename T>
/*inline*/ void do_mat_mat_mult (const unsigned start, const unsigned end,
                                 TMatrix<T> *res, 
                                 const Matrix<T> *a, const Matrix<T> *b)
{
        const unsigned bc = b->columns(); 
        const unsigned ac = a->columns();
        for (unsigned i=start; i<end; ++i) {
                const T* RESTRICT aptr = a->getcolptr(i);
                PREFETCH_R(aptr,2);
                for (unsigned l=0; l<ac; ++l) {
                        T* RESTRICT rptr = res->getcolptr(i);
                        const T* RESTRICT bptr = b->getcolptr(l);
                        int n = bc - 8;
                        PREFETCH_R(rptr,2); PREFETCH_R(bptr,1);
                        ALIGN3(const register T tmp, *aptr, MIN_ALIGN2);
                        if (LIKELY(n >= 0)) {
                                if (sizeof(T)*2 >= CACHELINE_SZ) {
                                        PREFETCH_R(rptr+2,2); PREFETCH_R(bptr+2,1);
                                }
                                PREFETCH_R(rptr+4,2); PREFETCH_R(bptr+4,1);
                                if (sizeof(T)*2 >= CACHELINE_SZ) {
                                        PREFETCH_R(rptr+6,2); PREFETCH_R(bptr+6,1);
                                }
                                for (; n > 3; n -= 4) {
                                        *rptr += tmp * *bptr;
                                        PREFETCH_R(rptr+ 8,2);
                                        *(rptr+1) += tmp * *(bptr+1);
                                        PREFETCH_R(bptr+ 8,1);
                                        *(rptr+2) += tmp * *(bptr+2);
                                        if (sizeof(T)*2 >= CACHELINE_SZ)
                                                PREFETCH_R(rptr+10,2);
                                        rptr += 4; 
                                        if (sizeof(T)*2 >= CACHELINE_SZ)
                                                PREFETCH_R(bptr+10,1);
                                        *(rptr-1) += tmp * *(bptr+3);
                                        bptr += 4;
                                }
                        }
                        if (LIKELY(l+4 < ac-CACHELINE_SZ/sizeof(T)))
                                PREFETCH_R(aptr+4, 2);
                        for (n += 8; n; --n)
                                *rptr++ += tmp * *bptr++;
                        aptr++;
                }
        }
}
#endif

/* costs of mat vec mult; r = row, c = col */
#define COST_MATVEC(r,c) (r*(COST_UNIT_STORE+COST_LOOP \
 +c*(3*COST_UNIT_LOAD+COST_CALL+COST_ADD+COST_MULT+COST_LOOP)))

template <typename T>
/*inline*/ void do_mat_vec_mult (const unsigned start, const unsigned end,
                                 TVector<T> *res, 
                                 const Matrix<T> *mat, const Vector<T> *vec)
{
#ifdef DEBUG_THREAD
        fprintf (stderr, "do_mat_vec_mult (pid %i): %p %p %p, %i - %i\n", getpid(), res, mat, vec, start, end);
#endif
        PREFETCH_R(vec->vec,3);
        const unsigned mc = mat->col;
#ifdef TBCI_SIMD_SUM
        if (UNLIKELY(sizeof(T) < TBCI_SIMD_ALIGN && mc % (TBCI_SIMD_ALIGN/sizeof(T)))) {
                for (unsigned rw = start; rw < end; ++rw) {
                        res->vec[rw] = T(0.0);
#if 0
                        if (!((unsigned long)mat->mat[rw] % TBCI_SIMD_ALIGN))
#endif                          
                                do_vec_mult<T>(mc, mat->mat[rw], vec->vec, res->vec[rw]);
#if 0
                        else
                                do_vec_mult_unaligned<T>(mc, mat->mat[rw], vec->vec, res->vec[rw]);
#endif
                }
        } else
#endif
        for (unsigned rw = start; rw < end; ++rw) {
#if 1           
                //PREFETCH_W(res->vec+rw,3);
                res->vec[rw] = T(0.0);
                do_vec_mult<T> (mc, mat->mat[rw], vec->vec, res->vec[rw]);
#else
                register const T* matptr = mat->mat[rw];
                PREFETCH_R(matptr, 1);
                register const T* vecptr = vec->vec;
                //PREFETCH_R(vecptr, 3);
                register int i = mc - 9;
                ALIGN2(T el (*matptr++ * *vecptr++), MIN_ALIGN2);
                if (i >= 0) {
                        if (sizeof(T)*2 >= CACHELINE_SZ)
                                PREFETCH_R(matptr+2, 1); 
                        PREFETCH_R(matptr+4, 1); 
                        if (sizeof(T)*2 >= CACHELINE_SZ)
                                PREFETCH_R(matptr+6, 1);
                        //PREFETCH_R(vecptr+2, 3); PREFETCH_R(vecptr+4, 3); PREFETCH_R(vecptr+6, 3);
                        for (; i > 3; i -= 4) {
                                el += *matptr * *vecptr;
                                PREFETCH_R(matptr+ 8,1);
                                el += *(matptr+1) * *(vecptr+1);
                                //PREFETCH_R(vecptr+ 8,3);
                                el += *(matptr+2) * *(vecptr+2);
                                if (sizeof(T)*2 >= CACHELINE_SZ)
                                        PREFETCH_R(vecptr+10,1);
                                el += *(matptr+3) * *(vecptr+3);
                                //PREFETCH_R(vecptr+10,3);
                                matptr += 4; vecptr += 4;
                        }
                }
                if (LIKELY(rw < end-CACHELINE_SZ/sizeof(T))) {
                        PREFETCH_W(res->vec+rw, 2);
                }
                for (i += 8; i; --i)
                        el += *matptr++ * *vecptr++;
                
                res->vec[rw] = el;
#endif
        }
}

template <typename T>
/*inline*/ void do_mat_tsv_mult (const unsigned start, const unsigned end,
                                   TVector<T> *res, 
                                   const Matrix<T> *mat, const TSVector<T> *tsv)
{
        PREFETCH_R(tsv->vec,3);
        const unsigned mc = mat->col;
        const T fact = tsv->getfac();
#ifdef TBCI_SIMD_SUM
        if (UNLIKELY(sizeof(T) < TBCI_SIMD_ALIGN && mc % (TBCI_SIMD_ALIGN/sizeof(T)))) {
                for (unsigned rw = start; rw < end; ++rw) {
                        res->vec[rw] = T(0.0);
#if 0
                        if (!((unsigned long)mat->mat[rw] % TBCI_SIMD_ALIGN))
#endif
                                do_vec_mult<T>(mc, mat->mat[rw], tsv->vec, res->vec[rw]);
#if 0
                        else
                                do_vec_mult_unaligned<T>(mc, mat->mat[rw], tsv->vec, res->vec[rw]);
#endif
                        res->vec[rw] *= fact;
                }
        } else
#endif
        for (unsigned rw = start; rw < end; ++rw) {
                //PREFETCH_W(res->vec+rw,3);
                res->vec[rw] = T(0.0);
                // TODO: May be unaligned!
                do_vec_mult<T> (mc, mat->mat[rw], tsv->vec, res->vec[rw]);
                res->vec[rw] *= fact;
        }
}

template <typename T>
/*inline*/ void do_mat_vec_transmult (const unsigned start, const unsigned end,
                                      TVector<T> *res, 
                                      const Matrix<T> *mat, const Vector<T> *vec)
{
        PREFETCH_R(vec->vec,3);
        const unsigned mr = mat->row;
        for (unsigned cl = start; cl < end; cl++)
        {
                PREFETCH_R(mat->mat[0]+cl,1);
                register T* vecptr = vec->vec;
                register int i = mr - 9;
                int r = 1;
                ALIGN3(T el, mat->mat[0][cl] * *vecptr++, MIN_ALIGN2);
                if (LIKELY(i >= 0)) {
                        PREFETCH_R(mat->mat[1]+cl,1);
                        PREFETCH_R(mat->mat[2]+cl,1);
                        PREFETCH_R(mat->mat[3]+cl,1);
                        PREFETCH_R(mat->mat[4]+cl,1);
                        PREFETCH_R(mat->mat[5]+cl,1);
                        PREFETCH_R(mat->mat[6]+cl,1);
                        PREFETCH_R(mat->mat[7]+cl,1);
                        for (; i > 3; i -= 4) {
                                el += mat->mat[r][cl] * *vecptr;
                                PREFETCH_R(mat->mat[r+ 8]+cl, 1);
                                el += mat->mat[r+1][cl] * *(vecptr+1);
                                PREFETCH_R(mat->mat[r+ 9]+cl, 1);
                                el += mat->mat[r+2][cl] * *(vecptr+2);
                                PREFETCH_R(mat->mat[r+10]+cl, 1);
                                el += mat->mat[r+3][cl] * *(vecptr+3);
                                PREFETCH_R(mat->mat[r+11]+cl, 1);
                                r += 4; vecptr += 4;
                        }
                }
                if (LIKELY(cl < end-CACHELINE_SZ/sizeof(T))) {
                        PREFETCH_W(res->vec + cl, 2);
                }
                for (i += 8; i; --i, r++)
                        el += mat->mat[r][cl] * *vecptr++;

                res->vec[cl] = el;
        }
}

/* TODO:
 * Divide and Conquer algorithm for matrix-matrix multiplication:
 * (Strassen's algorithm)
 *
 * [C11   C12 ]     [A11   A12 ]     [B11   B12 ]
 * |          |  =  |          |  ×  |          |
 * [C21   C22 ]     [A21   A22 ]     [B21   B22 ]
 * 
 *          P1   = (A11 + A22) (B11 + B22)
 *          P2   = (A21 + A22) B11
 *          P3   = A11 (B12 - B22)
 *          P4   = A22( B21 - B11)
 *          P5   = (A11 + A12) B22
 *          P6   = (A21 - A11) (B11 + B12)
 *          P7   = (A12 - A22) (B21 + B22)
 * 
 *         C11   = P1 + P4 - P5 + P7
 *         C12   = P3 + P5
 *         C21   = P2 + P4
 *         C22   = P1 + P3 - P2 + P6
 *
 * O(n^2log7) multiplications
 * 
 *  http://www.cse.ohio-state.edu/~gurari/course/cse693s04/cse693s04su61.html#x84-7200011.4
 *  http://www.brpreiss.com/books/opus5/html/page457.html
 */

#endif  /* TBCI_MATRIX_KERNELS_H */

---