qcd.cpp

c++经典教材 Blitz++ v0.8

/*****************************************************************************
 * qcd.cpp     Blitz++ TinyMatrix<T> example.  
 *
 * This example illustrates the TinyMatrix<T,N,M> class.  TinyMatrix can
 * be used for small matrices whose dimensions are known at compile time.
 * The application is from lattice quantum chromodynamics.
 *
 * See: S. Booth, Lattice QCD simulation programs on the Cray T3D,
 * Edinburgh Parallel Computing Centre Technical Report EPCC-TR96-03.  
 * http://www.epcc.ed.ac.uk/t3d/documents/techreports.html
 *****************************************************************************/

#include <blitz/blitz.h>

#ifdef BZ_HAVE_COMPLEX

#include <blitz/tinymat.h>
#include <blitz/vector-et.h>

BZ_USING_NAMESPACE(blitz)

/*
 * This example is drawn from a technical report by Stephen Booth of
 * the Edinburgh Parallel Computing Centre.  In his group's implementation 
 * of lattice QCD, most of the CPU time is taken up by a group of routines 
 * which multiplies 2-spinors by SU(3) gauge elements.  In Fortran, this 
 * code looks like:
 *
 * COMPLEX M(V,3,3) res(V,3,2), src(V,3,2)
 *   DO spin=1,2
 *       DO col=1,3
 *           DO site=1,V
 *               res(site,col,spin)=
 *                        M(site,col,1) * src(site,1,spin)
 *                      + M(site,col,2) * src(site,2,spin)
 *                      + M(site,col,3) * src(site,3,spin)
 *           END DO
 *       END DO
 *   END DO
 *
 * This is a situation in which TinyMatrix can be used to good
 * effect: each element of `src' is a 3x2 complex matrix, and
 * each element of `M' is a 3x3 complex matrix.  
 */

typedef TinyMatrix<complex<double>, 3, 3> gauge;
typedef TinyMatrix<complex<double>, 3, 2> spinor;

/*
 * The following two-line function is equivalent to the Fortran version
 * above.
 */

inline void qcdUpdate(Vector<spinor>& res, Vector<gauge>& M, 
    Vector<spinor>& src)
{
    for (int i=0; i < res.length(); ++i)
        res[i] = product(M[i], src[i]);
}

int main()
{
    const int numElements = 100;   // Some number of elements

    Vector<gauge> M(numElements);
    Vector<spinor> res(numElements), src(numElements);

    qcdUpdate(res, M, src);

    return 0;
}

/*
 * Since the dimensions of gauge and spinor are known at compile
 * time, the Blitz++ library unrolls the matrix multiplication
 * loops.  Here is the assembler code generated for qcdUpdate(..) 
 * using KAI C++ on an RS/6000.
 *
 * fm    floating-point multiply
 * fma   floating-point multiply and add
 * lfd   load double-precision floating point
 * stfd  store double-precision floating point
 *
 * Note the use of pre-fetching and quad operand instructions.
 */

#if 0
        l       r10,4(r28)
        l       r9,0(r28)
        cal     r11,100(r0)
        cmpi    0,r10,0
        ai      r7,r27,-40
        ai      r8,r29,-8
        ai      r9,r9,-72
        mtspr   CTR,r11
__L740:                                 # 0x00000740 (H.10.NO_SYMBOL+0x740)
        lfd     fp2,68(r7)
        lfd     fp25,116(r9)
        lfd     fp26,108(r9)
        fm      fp3,fp2,fp25
        lfd     fp6,60(r7)
        fm      fp0,fp26,fp2
        lfd     fp1,84(r7)
        fms     fp9,fp6,fp26,fp3
        lfd     fp23,148(r9)
        fma     fp10,fp6,fp25,fp0
        lfd     fp24,140(r9)
        fm      fp3,fp1,fp23
        lfd     fp7,76(r7)
        fm      fp0,fp24,fp1
        lfd     fp5,52(r7)
        fms     fp3,fp7,fp24,fp3
        lfd     fp28,84(r9)
        fma     fp0,fp7,fp23,fp0
        lfd     fp27,76(r9)
        fm      fp4,fp5,fp28
        lfd     fp8,44(r7)
        fa      fp3,fp9,fp3
        fms     fp9,fp8,fp27,fp4
        fm      fp4,fp27,fp5
        fa      fp0,fp10,fp0
        fa      fp3,fp9,fp3
        fma     fp4,fp8,fp28,fp4
        stfd    fp3,12(r8)
        fa      fp0,fp0,fp4
        stfd    fp0,20(r8)
        lfd     fp30,124(r9)
        lfd     fp13,132(r9)
        fm      fp0,fp2,fp30
        lfd     fp12,156(r9)
        fm      fp4,fp2,fp13
        lfd     fp29,164(r9)
        fma     fp3,fp6,fp13,fp0
        lfd     fp31,92(r9)
        fm      fp0,fp1,fp12
        lfdu    fp11,100(r9)
        fm      fp2,fp1,fp29
        fma     fp1,fp7,fp29,fp0
        fm      fp0,fp5,fp31
        fms     fp4,fp6,fp30,fp4
        fms     fp2,fp7,fp12,fp2
        fm      fp5,fp5,fp11
        fa      fp1,fp3,fp1
        fma     fp0,fp8,fp11,fp0
        fa      fp2,fp4,fp2
        fms     fp3,fp8,fp31,fp5
        fa      fp0,fp0,fp1
        fa      fp1,fp2,fp3
        stfd    fp0,36(r8)
        stfd    fp1,28(r8)
        lfd     fp9,132(r7)
        lfd     fp22,116(r7)
        fm      fp1,fp23,fp9
        lfd     fp7,100(r7)
        fm      fp3,fp25,fp22
        lfd     fp8,108(r7)
        fm      fp0,fp26,fp22
        lfd     fp10,124(r7)
        fm      fp2,fp24,fp9
        lfd     fp5,92(r7)
        fms     fp3,fp26,fp8,fp3
        fms     fp1,fp24,fp10,fp1
        fma     fp0,fp25,fp8,fp0
        fma     fp2,fp23,fp10,fp2
        fm      fp4,fp28,fp7
        fm      fp6,fp27,fp7
        fa      fp1,fp3,fp1
        fa      fp0,fp0,fp2
        fma     fp3,fp28,fp5,fp6
        fms     fp2,fp27,fp5,fp4
        fa      fp0,fp0,fp3
        fa      fp1,fp1,fp2
        fm      fp2,fp13,fp22
        stfd    fp1,44(r8)
        stfd    fp0,52(r8)
        fms     fp1,fp30,fp8,fp2
        fm      fp0,fp29,fp9
        fm      fp3,fp12,fp9
        fm      fp2,fp30,fp22
        fms     fp0,fp12,fp10,fp0
        fma     fp4,fp29,fp10,fp3
        fma     fp2,fp13,fp8,fp2
        fm      fp3,fp31,fp7
        fa      fp0,fp1,fp0
        fm      fp6,fp11,fp7
        fa      fp1,fp2,fp4
        fma     fp2,fp11,fp5,fp3
        fms     fp3,fp31,fp5,fp6
        fa      fp1,fp1,fp2
        fa      fp0,fp0,fp3
        stfd    fp1,68(r8)
        stfd    fp0,60(r8)
        lfd     fp7,164(r7)
        lfd     fp2,180(r7)
        fm      fp4,fp25,fp7
        lfd     fp0,156(r7)
        fm      fp5,fp23,fp2
        lfd     fp1,172(r7)
        fm      fp3,fp26,fp7
        fms     fp6,fp26,fp0,fp4
        fms     fp5,fp24,fp1,fp5
        fma     fp4,fp25,fp0,fp3
        fa      fp3,fp6,fp5
        fm      fp9,fp24,fp2
        lfd     fp8,140(r7)
        fm      fp6,fp13,fp7
        lfdu    fp5,148(r7)
        fm      fp7,fp30,fp7
        fma     fp25,fp23,fp1,fp9
        fm      fp9,fp29,fp2
        fm      fp26,fp28,fp5
        fm      fp10,fp27,fp5
        fma     fp7,fp13,fp0,fp7
        fm      fp2,fp12,fp2
        fms     fp13,fp27,fp8,fp26
        fms     fp9,fp12,fp1,fp9
        fa      fp4,fp4,fp25
        fma     fp10,fp28,fp8,fp10
        fma     fp1,fp29,fp1,fp2
        fms     fp0,fp30,fp0,fp6
        fa      fp2,fp13,fp3
        fm      fp6,fp11,fp5
        fa      fp3,fp10,fp4
        fm      fp4,fp31,fp5
        fa      fp0,fp0,fp9
        fms     fp5,fp31,fp8,fp6
        fma     fp4,fp11,fp8,fp4
        stfd    fp2,76(r8)
        stfd    fp3,84(r8)
        fa      fp1,fp7,fp1
        fa      fp0,fp5,fp0
        fa      fp1,fp4,fp1
        stfd    fp0,92(r8)
        stfdu   fp1,100(r8)
        bc      BO_dCTR_NZERO,CR0_LT,__L740

#endif

#else // BZ_HAVE_COMPLEX

#include <iostream.h>

int main()
{
    cout << "This example requires <complex> from the ISO/ANSI C++ standard."
         << endl;

    return 0;
}

#endif