Skip to content
Snippets Groups Projects
sequential.cpp 21.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • Camille Coti's avatar
    Camille Coti committed
    #include <ginac/ginac.h>
    
    
    Camille Coti's avatar
    Camille Coti committed
    #include "tensormatrix.h"
    #include "products.h"
    
    namespace gi = GiNaC;
    
    /*******************************************************************************
     *                              Big multiplication                             *
     *******************************************************************************/
    
    
    Camille Coti's avatar
    Camille Coti committed
    gi::ex multiply_seq( tensor3D_t& T, int size ) { 
    
    Camille Coti's avatar
    Camille Coti committed
    
        gi::ex Tens = 0;
    	int a1, a2, a3, a4, a5, a6;
        int A1, A2, A3, A4, A5, A6;
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
    
        
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
    
        int N = size/2;
        
        for( a4 = 0 ; a4 < N ; a4++ ) {
            A4 = a4 + N;
            for( a2 = 0 ; a2 < N ; a2++ ) {
                A2 = a2 + N;
                for( a6 = 0 ; a6 < N ; a6++ ) {
                    A6 = a6 + N;
                    W1 = T[a4][a2][a6];
                    W2 = T[a4][A2][a6];
                    W3 = T[a4][a2][A6];
                    W4 = T[A4][A2][a6];
                    W5 = T[a4][A2][A6];
                    W6 = T[A4][a2][A6];
                    W7 = T[A4][A2][A6];
                    for( a1 = 0 ; a1 < N ; a1++ ) {
                        A1 = a1 + N;
                        for( a5 = 0 ; a5 < N ; a5++ ) {
                            A5 = a5 + N;
                            Z1 = T[a1][a5][a6];
                            Z2 = T[A1][a5][a6];
                            Z6 = T[A1][a5][A6];
                            t5 = W3*T[a1][A5][a6];
                            tE = W4*T[A1][A5][A6];
                            t1 = W3*Z2;
                            t13 = t1;
                            t2 = W5*Z1;
                            t23 = t2;
                            t3 = W3*Z1;
                            t4 = W6*Z1;
                            t12 = W5*Z2;
                            t14 = W6*Z2;
                            t134 = t14 ;
                            t16 = W1*Z6;
                            t24 = W7*Z1;
                            t26 = W2*T[a1][a5][A6];
                            t123 = W5*Z2;
                            t126 = W2*Z6;
                            
                            for( a3 = 0 ; a3 < N ; a3++ ) {
                                A3 = a3 + N;
                                TE = TE + tE*T[a1][a2][a3]*T[a4][a5][A3];
                                T5 = T5 + t5*T[A1][A2][A3]*T[A4][a5][a3];
                                X7Y5 = T[a1][A2][A3]*T[A4][A5][a3];
                                T1 = T1 + t1*X7Y5;
                                T16 = T16 + t16*X7Y5;
                                T2 = T2 + t2*T[A1][a2][A3]*T[A4][A5][a3];
                                T3 = T3 + t3*T[A1][A2][a3]*T[A4][A5][A3];
                                T4 = T4 + t4*T[A1][A2][A3]*T[a4][A5][a3];
                                T12 = T12 + t12*T[a1][a2][A3]*T[A4][A5][a3];
                                T13 = T13 + t13*T[a1][A2][a3]*T[A4][A5][A3];
                                T14 = T14 + t14*T[a1][A2][A3]*T[a4][A5][a3];
                                T23 = T23 + t23*T[A1][a2][a3]*T[A4][A5][A3];
                                T24 = T24 + t24*T[A1][a2][A3]*T[a4][A5][a3];
                                T26 = T26 + t26*T[A1][a2][A3]*T[A4][A5][a3];
                                T123 = T123 + t123*T[a1][a2][a3]*T[A4][A5][A3];
                                T126 = T126 + t126*T[a1][a2][A3]*T[A4][A5][a3];
                                T134 = T134 + t134*T[a1][A2][a3]*T[a4][A5][A3];
                            }
                        }
                    }
                }
            }
        }
        Tens = 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134));
    
        return Tens;
    }
    
    
    Camille Coti's avatar
    Camille Coti committed
    /*******************************************************************************
     *                             1-level decomposition                           *
     *******************************************************************************/
    
    
    Camille Coti's avatar
    Camille Coti committed
    gi::ex multiply_1level( tensor3D_t& T, int size ) {  // simpler: same dimension everywhere
    
    Camille Coti's avatar
    Camille Coti committed
        
        gi::ex Tens = 0;
        gi::ex Tn;
    
    
    #ifdef TAUPROF
        void* ptr;
        TAU_PROFILER_CREATE( ptr, "b","", TAU_USER );
    #endif
    
        for( a4 = 0 ; a4 < N ; a4++ ) {
            for( a2 = 0 ; a2 < N ; a2++ ) {
    
    #ifdef TAUPROF
                    TAU_START( "b" );
                    TAU_PROFILER_START( ptr );
    #endif
    
    Camille Coti's avatar
    Camille Coti committed
                    /* Beyond this point, a2 and a3 are only used in the simplectic matrix */
    
    				for( a1 = 0 ; a1 < N ; a1++ ){
                        Tn = one_level1_product( &T, size, a4, a2, a1 );
    
    Camille Coti's avatar
    Camille Coti committed
                        Tens += Tn;
                    }
    
    #ifdef TAUPROF
                    TAU_STOP( "b" );
                    //                time = getTimeSpent( ptr );
                    
                      long calls, childcalls;
      double incl[TAU_MAX_COUNTERS], excl[TAU_MAX_COUNTERS];
      const char **counters;
      int numcounters, i, j; 
    
      TAU_PROFILER_GET_CALLS(ptr, &calls);
      TAU_PROFILER_GET_CHILD_CALLS(ptr, &childcalls);
      TAU_PROFILER_GET_INCLUSIVE_VALUES(ptr, &incl);
      TAU_PROFILER_GET_EXCLUSIVE_VALUES(ptr, &excl);
    
      TAU_PROFILER_GET_COUNTER_INFO(&counters, &numcounters);
      printf("Calls = %ld, child = %ld\n", calls, childcalls);
      printf("numcounters = %d\n", numcounters);
      for (j = 0; j < numcounters ; j++)  {
        printf(">>>");
        printf("counter [%d] = %s\n", j, counters[j]);
        printf(" excl [%d] = %g, incl [%d] = %g\n", j, excl[j], j, incl[j]);
      }
    
    //                time = getTimeSpent( "b" );
    
    Camille Coti's avatar
    Camille Coti committed
        return Tens;
    }
    
    
    gi::ex one_level1_product( tensor3D_t* T, int size, int a4, int a2, int a1 ){
    
    Camille Coti's avatar
    Camille Coti committed
    
        gi::ex Tens = 0;
    
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
    
        
    	int a3, a5, a6;
        int A1, A2, A3, A4, A5, A6;
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
    
        int N = size/2;
        
        A1 = a1 + N;
        A4 = a4 + N;
        A2 = a2 + N;
        for( a6 = 0 ; a6 < N ; a6++ ) {
            A6 = a6 + N;
            
            W1 = (*T)[a4][a2][a6];
            W2 = (*T)[a4][A2][a6];
            W3 = (*T)[a4][a2][A6];
            W4 = (*T)[A4][A2][a6];
            W5 = (*T)[a4][A2][A6];
            W6 = (*T)[A4][a2][A6];
            W7 = (*T)[A4][A2][A6];
            
            for( a5 = 0 ; a5 < N ; a5++ ) {
                A5 = a5 + N;
                Z1 = (*T)[a1][a5][a6];
                Z2 = (*T)[A1][a5][a6];
                Z6 = (*T)[A1][a5][A6];
                t5 = W3*(*T)[a1][A5][a6];
                tE = W4*(*T)[A1][A5][A6];
                t1 = W3*Z2;
                t13 = t1;
                t2 = W5*Z1;
                t23 = t2;
                t3 = W3*Z1;
                t4 = W6*Z1;
                t12 = W5*Z2;
                t14 = W6*Z2;
                t134 = t14 ;
                t16 = W1*Z6;
                t24 = W7*Z1;
                t26 = W2*(*T)[a1][a5][A6];
                t123 = W5*Z2;
                t126 = W2*Z6;
                
                for( a3 = 0 ; a3 < N ; a3++ ) {
                    A3 = a3 + N;
                    TE = TE + tE*(*T)[a1][a2][a3]*(*T)[a4][a5][A3];
                    T5 = T5 + t5*(*T)[A1][A2][A3]*(*T)[A4][a5][a3];
                    X7Y5 = (*T)[a1][A2][A3]*(*T)[A4][A5][a3];
                    T1 = T1 + t1*X7Y5;
                    T16 = T16 + t16*X7Y5;
                    T2 = T2 + t2*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                    T3 = T3 + t3*(*T)[A1][A2][a3]*(*T)[A4][A5][A3];
                    T4 = T4 + t4*(*T)[A1][A2][A3]*(*T)[a4][A5][a3];
                    T12 = T12 + t12*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                    T13 = T13 + t13*(*T)[a1][A2][a3]*(*T)[A4][A5][A3];
                    T14 = T14 + t14*(*T)[a1][A2][A3]*(*T)[a4][A5][a3];
                    T23 = T23 + t23*(*T)[A1][a2][a3]*(*T)[A4][A5][A3];
                    T24 = T24 + t24*(*T)[A1][a2][A3]*(*T)[a4][A5][a3];
                    T26 = T26 + t26*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                    T123 = T123 + t123*(*T)[a1][a2][a3]*(*T)[A4][A5][A3];
                    T126 = T126 + t126*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                    T134 = T134 + t134*(*T)[a1][A2][a3]*(*T)[a4][A5][A3];
    
        Tens = 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134));
    
    
    #if 0
        std::ostringstream oss;
        oss << "output_" << getpid();
        std::ofstream fd;
        fd.open( oss.str(), std::ios::app	);
        fd << "T" << A << "=" << Tens << ";";
        //    fd << " with " << A << " " << a1 <<  " " << a2 <<  " " << a3 <<  " " << b1;
        fd << std::endl;
        fd.close();
    #endif
        
        return Tens;
    }
    
    gi::ex one_level1_product( tensor3D_t* T, int size, int a4 ){
    
        gi::ex Tens = 0;
        gi::ex Ti0, Ti1, Ti2;
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
        const char timerB[] = "B";
        
    	int a1, a2, a3, a5, a6;
        int A1, A2, A3, A4, A5, A6;
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
        Ti0 = Ti1 = Ti2 = 0;
        
        int N = size/2;
        
        A4 = a4 + N;
        for( a2 = 0 ; a2 < N ; a2++ ) {
            A2 = a2 + N;
            Ti2 = 0;
            for( a6 = 0 ; a6 < N ; a6++ ) {
                A6 = a6 + N;
                
                W1 = (*T)[a4][a2][a6];
                W2 = (*T)[a4][A2][a6];
                W3 = (*T)[a4][a2][A6];
                W4 = (*T)[A4][A2][a6];
                W5 = (*T)[a4][A2][A6];
                W6 = (*T)[A4][a2][A6];
                W7 = (*T)[A4][A2][A6];
                
                Ti1 = 0;
                for( a1 = 0 ; a1 < N ; a1++ ) {
                    A1 = a1 + N;
                    Ti0 = TE = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T1 = T2 = T3 = T4 = T5 = T123 = T126 = T134 = 0;
                    for( a5 = 0 ; a5 < N ; a5++ ) {
                        A5 = a5 + N;
                        Z1 = (*T)[a1][a5][a6];
                        Z2 = (*T)[A1][a5][a6];
                        Z6 = (*T)[A1][a5][A6];
                        t5 = W3*(*T)[a1][A5][a6];
                        tE = W4*(*T)[A1][A5][A6];
                        t1 = W3*Z2;
                        t13 = t1;
                        t2 = W5*Z1;
                        t23 = t2;
                        t3 = W3*Z1;
                        t4 = W6*Z1;
                        t12 = W5*Z2;
                        t14 = W6*Z2;
                        t134 = t14 ;
                        t16 = W1*Z6;
                        t24 = W7*Z1;
                        t26 = W2*(*T)[a1][a5][A6];
                        t123 = W5*Z2;
                        t126 = W2*Z6;
                        
                        for( a3 = 0 ; a3 < N ; a3++ ) {
                            A3 = a3 + N;
                            TE = TE + tE*(*T)[a1][a2][a3]*(*T)[a4][a5][A3];
                            T5 = T5 + t5*(*T)[A1][A2][A3]*(*T)[A4][a5][a3];
                            X7Y5 = (*T)[a1][A2][A3]*(*T)[A4][A5][a3];
                            T1 = T1 + t1*X7Y5;
                            T16 = T16 + t16*X7Y5;
                            T2 = T2 + t2*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                            T3 = T3 + t3*(*T)[A1][A2][a3]*(*T)[A4][A5][A3];
                            T4 = T4 + t4*(*T)[A1][A2][A3]*(*T)[a4][A5][a3];
                            T12 = T12 + t12*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                            T13 = T13 + t13*(*T)[a1][A2][a3]*(*T)[A4][A5][A3];
                            T14 = T14 + t14*(*T)[a1][A2][A3]*(*T)[a4][A5][a3];
                            T23 = T23 + t23*(*T)[A1][a2][a3]*(*T)[A4][A5][A3];
                            T24 = T24 + t24*(*T)[A1][a2][A3]*(*T)[a4][A5][a3];
                            T26 = T26 + t26*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                            T123 = T123 + t123*(*T)[a1][a2][a3]*(*T)[A4][A5][A3];
                            T126 = T126 + t126*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                            T134 = T134 + t134*(*T)[a1][A2][a3]*(*T)[a4][A5][A3];
                        }
                        Ti0 += ( 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134)) );
                    }
                    Ti1 += Ti0;
                }
                Ti2 += Ti1;
            }
            Tens += Ti2;
        }
    
    
    Camille Coti's avatar
    Camille Coti committed
    #if 0
        std::ostringstream oss;
        oss << "output_" << getpid();
        std::ofstream fd;
        fd.open( oss.str(), std::ios::app	);
        fd << "T" << A << "=" << Tens << ";";
        //    fd << " with " << A << " " << a1 <<  " " << a2 <<  " " << a3 <<  " " << b1;
        fd << std::endl;
        fd.close();
    #endif
        
        return Tens;
    }
    
    
    gi::ex one_level2_product( tensor3D_t* T, int size, int a4, int a2, int a1, int a6 ){
    
        int a3, a5;
        int A1, A2, A3, A4, A5, A6;
    
    Camille Coti's avatar
    Camille Coti committed
    
        gi::ex Tens = 0;
    
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
    
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
    
        int N = size/2;
    
        A1 = a1 + N;
        A4 = a4 + N;
        A2 = a2 + N;
        A6 = a6 + N;
    
        W1 = (*T)[a4][a2][a6];
        W2 = (*T)[a4][A2][a6];
        W3 = (*T)[a4][a2][A6];
        W4 = (*T)[A4][A2][a6];
        W5 = (*T)[a4][A2][A6];
        W6 = (*T)[A4][a2][A6];
        W7 = (*T)[A4][A2][A6];
        
        
        for( a5 = 0 ; a5 < N ; a5++ ) {
            A5 = a5 + N;
            Z1 = (*T)[a1][a5][a6];
            Z2 = (*T)[A1][a5][a6];
            Z6 = (*T)[A1][a5][A6];
            t5 = W3*(*T)[a1][A5][a6];
            tE = W4*(*T)[A1][A5][A6];
            t1 = W3*Z2;
            t13 = t1;
            t2 = W5*Z1;
            t23 = t2;
            t3 = W3*Z1;
            t4 = W6*Z1;
            t12 = W5*Z2;
            t14 = W6*Z2;
            t134 = t14 ;
            t16 = W1*Z6;
            t24 = W7*Z1;
            t26 = W2*(*T)[a1][a5][A6];
            t123 = W5*Z2;
            t126 = W2*Z6;
            
            for( a3 = 0 ; a3 < N ; a3++ ) {
                A3 = a3 + N;
                TE = TE + tE*(*T)[a1][a2][a3]*(*T)[a4][a5][A3];
                T5 = T5 + t5*(*T)[A1][A2][A3]*(*T)[A4][a5][a3];
                X7Y5 = (*T)[a1][A2][A3]*(*T)[A4][A5][a3];
                T1 = T1 + t1*X7Y5;
                T16 = T16 + t16*X7Y5;
                T2 = T2 + t2*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                T3 = T3 + t3*(*T)[A1][A2][a3]*(*T)[A4][A5][A3];
                T4 = T4 + t4*(*T)[A1][A2][A3]*(*T)[a4][A5][a3];
                T12 = T12 + t12*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                T13 = T13 + t13*(*T)[a1][A2][a3]*(*T)[A4][A5][A3];
                T14 = T14 + t14*(*T)[a1][A2][A3]*(*T)[a4][A5][a3];
                T23 = T23 + t23*(*T)[A1][a2][a3]*(*T)[A4][A5][A3];
                T24 = T24 + t24*(*T)[A1][a2][A3]*(*T)[a4][A5][a3];
                T26 = T26 + t26*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                T123 = T123 + t123*(*T)[a1][a2][a3]*(*T)[A4][A5][A3];
                T126 = T126 + t126*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                T134 = T134 + t134*(*T)[a1][A2][a3]*(*T)[a4][A5][A3];
    
        Tens = 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134));
        
    
    Camille Coti's avatar
    Camille Coti committed
        return Tens;
    }
    
    /*******************************************************************************
     *                             2-level decomposition                           *
     *******************************************************************************/
    
    
    Camille Coti's avatar
    Camille Coti committed
    gi::ex multiply_2levels( tensor3D_t& T, int size ) {  // simpler: same dimension everywhere
    
    
    Camille Coti's avatar
    Camille Coti committed
        gi::ex Tens = 0;
    
    Camille Coti's avatar
    Camille Coti committed
    	int a2, a4;
    
    Camille Coti's avatar
    Camille Coti committed
        int N = size/2;
    
    Camille Coti's avatar
    Camille Coti committed
        for( a4 = 0 ; a4 < N ; a4++ ) {
            for( a2 = 0 ; a2 < N ; a2++ ) {
                Tens += two_level1_product( &T, size, a2, a4 );
    
    Camille Coti's avatar
    Camille Coti committed
    		}
    	}
        return Tens;
    }
    
    
    gi::ex two_level1_product( tensor3D_t* T, int size, int a4, int a2 ){
    
    Camille Coti's avatar
    Camille Coti committed
    
        gi::ex Tens = 0;
    
    Camille Coti's avatar
    Camille Coti committed
    	int a1, a6;
    
        int N = size/2;
        for( a6 = 0 ; a6 < N ; a6++ ) {
            for( a1 = 0 ; a1 < N ; a1++ ) {
                
                Tens += two_level2_product( T, size, a2, a4, a6, a1 );
    
    Camille Coti's avatar
    Camille Coti committed
    
    
    Camille Coti's avatar
    Camille Coti committed
        return Tens;
    }
    
    
    gi::ex two_level2_product( tensor3D_t* T, int size, int a4, int a2, int a1, int a6 ) {
    
    Camille Coti's avatar
    Camille Coti committed
        gi::ex Tens = 0;
    	int a3, a5;
        int A1, A2, A3, A4, A5, A6;
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
        
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
    
    Camille Coti's avatar
    Camille Coti committed
        int N = size/2;
    
        A4 = a4 + N;
        A2 = a2 + N;
        A1 = a1 + N;
        A6 = a6 + N;
    
    Camille Coti's avatar
    Camille Coti committed
        W1 = (*T)[a4][a2][a6];
        W2 = (*T)[a4][A2][a6];
        W3 = (*T)[a4][a2][A6];
        W4 = (*T)[A4][A2][a6];
        W5 = (*T)[a4][A2][A6];
        W6 = (*T)[A4][a2][A6];
        W7 = (*T)[A4][A2][A6];
        
        for( a5 = 0 ; a5 < N ; a5++ ) {
            A5 = a5 + N;
            Z1 = (*T)[a1][a5][a6];
            Z2 = (*T)[A1][a5][a6];
            Z6 = (*T)[A1][a5][A6];
            t5 = W3*(*T)[a1][A5][a6];
            tE = W4*(*T)[A1][A5][A6];
            t1 = W3*Z2;
            t13 = t1;
            t2 = W5*Z1;
            t23 = t2;
            t3 = W3*Z1;
            t4 = W6*Z1;
            t12 = W5*Z2;
            t14 = W6*Z2;
            t134 = t14 ;
            t16 = W1*Z6;
            t24 = W7*Z1;
            t26 = W2*(*T)[a1][a5][A6];
            t123 = W5*Z2;
            t126 = W2*Z6;
                            
            for( a3 = 0 ; a3 < N ; a3++ ) {
                A3 = a3 + N;
                TE = TE + tE*(*T)[a1][a2][a3]*(*T)[a4][a5][A3];
                T5 = T5 + t5*(*T)[A1][A2][A3]*(*T)[A4][a5][a3];
                X7Y5 = (*T)[a1][A2][A3]*(*T)[A4][A5][a3];
                T1 = T1 + t1*X7Y5;
                T16 = T16 + t16*X7Y5;
                T2 = T2 + t2*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                T3 = T3 + t3*(*T)[A1][A2][a3]*(*T)[A4][A5][A3];
                T4 = T4 + t4*(*T)[A1][A2][A3]*(*T)[a4][A5][a3];
                T12 = T12 + t12*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                T13 = T13 + t13*(*T)[a1][A2][a3]*(*T)[A4][A5][A3];
                T14 = T14 + t14*(*T)[a1][A2][A3]*(*T)[a4][A5][a3];
                T23 = T23 + t23*(*T)[A1][a2][a3]*(*T)[A4][A5][A3];
                T24 = T24 + t24*(*T)[A1][a2][A3]*(*T)[a4][A5][a3];
                T26 = T26 + t26*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                T123 = T123 + t123*(*T)[a1][a2][a3]*(*T)[A4][A5][A3];
                T126 = T126 + t126*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                T134 = T134 + t134*(*T)[a1][A2][a3]*(*T)[a4][A5][A3];
    
    Camille Coti's avatar
    Camille Coti committed
            }
    
    Camille Coti's avatar
    Camille Coti committed
        }
    
        Tens = 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134));
    
    Camille Coti's avatar
    Camille Coti committed
        return Tens;
    }
    
    
    gi::ex two_level2_product( tensor3D_t* T, int size, int a4, int a2 ) {
        gi::ex Tens = 0;
    	int a3, a5, a1, a6;
        int A1, A2, A3, A4, A5, A6;
        gi::ex W1, W2, W3, W4, W5, W6, W7;
        gi::ex Z1, Z2, Z6, t5, tE, t1, t12, t123, t126, t13, t134, t14, t16, t2, t23, t24, t26, t3, t4, X7Y5;
        
        gi::ex TE, T1, T2, T3, T4, T5, T12, T13, T14, T16, T23, T24, T26, T123, T126, T134;
        TE = T1 = T2 = T3 = T4 = T5 = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T123 = T126 = T134 = 0;
        gi::ex Ti0, Ti1;
    
        int N = size/2;
    
        A4 = a4 + N;
        A2 = a2 + N;
        
        for( a6 = 0 ; a6 < N ; a6++ ) {
            A6 = a6 + N;
            
            W1 = (*T)[a4][a2][a6];
            W2 = (*T)[a4][A2][a6];
            W3 = (*T)[a4][a2][A6];
            W4 = (*T)[A4][A2][a6];
            W5 = (*T)[a4][A2][A6];
            W6 = (*T)[A4][a2][A6];
            W7 = (*T)[A4][A2][A6];
            
            Ti1 = 0;
            for( a1 = 0 ; a1 < N ; a1++ ) {
                A1 = a1 + N;
                Ti0 = TE = T12 = T13 = T14 = T16 = T23 = T24 = T26 = T1 = T2 = T3 = T4 = T5 = T123 = T126 = T134 = 0;
                for( a5 = 0 ; a5 < N ; a5++ ) {
                    A5 = a5 + N;
                    Z1 = (*T)[a1][a5][a6];
                    Z2 = (*T)[A1][a5][a6];
                    Z6 = (*T)[A1][a5][A6];
                    t5 = W3*(*T)[a1][A5][a6];
                    tE = W4*(*T)[A1][A5][A6];
                    t1 = W3*Z2;
                    t13 = t1;
                    t2 = W5*Z1;
                    t23 = t2;
                    t3 = W3*Z1;
                    t4 = W6*Z1;
                    t12 = W5*Z2;
                    t14 = W6*Z2;
                    t134 = t14 ;
                    t16 = W1*Z6;
                    t24 = W7*Z1;
                    t26 = W2*(*T)[a1][a5][A6];
                    t123 = W5*Z2;
                    t126 = W2*Z6;
                    
                    for( a3 = 0 ; a3 < N ; a3++ ) {
                        A3 = a3 + N;
                        TE = TE + tE*(*T)[a1][a2][a3]*(*T)[a4][a5][A3];
                        T5 = T5 + t5*(*T)[A1][A2][A3]*(*T)[A4][a5][a3];
                        X7Y5 = (*T)[a1][A2][A3]*(*T)[A4][A5][a3];
                        T1 = T1 + t1*X7Y5;
                        T16 = T16 + t16*X7Y5;
                        T2 = T2 + t2*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                        T3 = T3 + t3*(*T)[A1][A2][a3]*(*T)[A4][A5][A3];
                        T4 = T4 + t4*(*T)[A1][A2][A3]*(*T)[a4][A5][a3];
                        T12 = T12 + t12*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                        T13 = T13 + t13*(*T)[a1][A2][a3]*(*T)[A4][A5][A3];
                        T14 = T14 + t14*(*T)[a1][A2][A3]*(*T)[a4][A5][a3];
                        T23 = T23 + t23*(*T)[A1][a2][a3]*(*T)[A4][A5][A3];
                        T24 = T24 + t24*(*T)[A1][a2][A3]*(*T)[a4][A5][a3];
                        T26 = T26 + t26*(*T)[A1][a2][A3]*(*T)[A4][A5][a3];
                        T123 = T123 + t123*(*T)[a1][a2][a3]*(*T)[A4][A5][A3];
                        T126 = T126 + t126*(*T)[a1][a2][A3]*(*T)[A4][A5][a3];
                        T134 = T134 + t134*(*T)[a1][A2][a3]*(*T)[a4][A5][A3];
                    }
                    Ti0 += ( 4*(TE+T12+T13+T14+T16+T23+T24+T26 - (T1 + T2 + T3 + T4 + T5 +T123 + T126 + T134)) );
                }
                Ti1 += Ti0;
            }
            Tens += Ti1;
        }
    
        return Tens;
    }