sequential.cpp

#include <ginac/ginac.h>

#include "profiling.h"
#include "tensormatrix.h"
#include "products.h"

namespace gi = GiNaC;

/*******************************************************************************
 *                              Big multiplication                             *
 *******************************************************************************/

gi::ex multiply_seq( tensor3D_t& T, matrix_int_t& J, int size ) {  // simpler: same dimension everywhere

    gi::ex Tens = 0;
	int a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3;
    gi::ex TAB, TABB, TABC, TABCC, TABCD, TABCDD; 
    gi::ex A;

 	int i, j;
    i = 0;
    j = 0;
    
    const char timerA[] = "A";
    const char timerB[] = "B";
    const char timeradd[] = "add";
    double timeA, timeB;
    uint64_t t_start, t_end;

    for( a1 = 0 ; a1 < size; a1++ ){
		i=i+1; 
        //  std::cout << "Tens: " << Tens << std::endl;
        //	printf("i = %d\n", i); 
		for( a2 = 0; a2 < size ; a2++ ){
			j=j+1; 
            //			printf("j = %d\n", j);
			for( a3 = 0 ; a3 < size ; a3++ ){
			    TAU_START( timerA );
			    A = T[a1][a2][a3];
			    /* Beyond this point, a2 and a3 are only used in the simplectic matrix */
				for( b1 = 0 ; b1 < size ; b1++ ){
					TAB = J[a1][b1]; 
					for( b2 = 0 ; b2 < size ; b2++ ){
						for( b3 = 0 ; b3 < size ; b3++ ){
						  TAU_START( timerB );
						  /* Beyond this point, b1 is not used anymore */
							TABB =  TAB * A*T[b1][b2][b3];
							for( c1 = 0 ; c1 < size ; c1++ ){
								for( c2 = 0 ; c2 < size ; c2++ ){
									TABC = TABB * J[a2][c2];
									for( c3 = 0 ; c3 < size ; c3++ ){
										TABCC = TABC * T[c1][c2][c3] * J[b3][c3] ; 
										for( d1 = 0 ; d1 < size ; d1++ ){
											TABCD = TABCC * J[c1][d1];
											for( d2 = 0 ; d2 < size ; d2++ ){
												TABCDD = TABCD * J[b2][d2];
												for( d3 = 0 ; d3 < size ; d3++ ){
												  TAU_START( timeradd );
												  t_start = rdtsc();
												  Tens = Tens + TABCDD * T[d1][d2][d3]*J[a3][d3];
												  t_end = rdtsc();
												  
												  TAU_STOP( timeradd );
#ifdef TAUPROF
						    //                                                    std::cout << "add " << getTimeSpent( timeradd ) << " len " << Tens.nops() << std::endl;
												  printf( "add %lf %lu len %d\n", getTimeSpent( timeradd ), t_end - t_start, Tens.nops() );

						    //                                                    std::cout << Tens << std::endl;
#endif // TAUPROF
   												}
											}
										}
									}
								}
							}
							TAU_STOP( timerB );
#ifdef TAUPROF
							std::cout << "B " << getTimeSpent( timeradd ) << " len " << Tens.nops() << std::endl;
#endif // TAUPROF
						}
					}
				}
				TAU_STOP( timerA );
#ifdef TAUPROF
				std::cout << "A " << getTimeSpent( timeradd ) << " len " << Tens.nops() << std::endl;
#endif // TAUPROF
			}
		}
    }
    return Tens;
}


gi::ex multiply_seq2( tensor3D_t& T, matrix_int_t& J, int size ) { 

    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;
}

/*******************************************************************************
 *                             1-level decomposition                           *
 *******************************************************************************/

gi::ex multiply_1level( tensor3D_t& T, matrix_int_t& J, int size ) {  // simpler: same dimension everywhere
    
    gi::ex Tens = 0;
    gi::ex Tn;
	int a1, a2, a3, b1;
    gi::ex A;
    
 	int i, j;
    i = 0;
    j = 0;

    double time;

#ifdef TAUPROF
    void* ptr;
    TAU_PROFILER_CREATE( ptr, "b","", TAU_USER );
#endif

    for( a1 = 0 ; a1 < size; a1++ ){
        i=i+1; 
        //   std::cout << "Tens: " << Tens << std::endl;
        //		printf("i = %d\n", i); 
		for( a2 = 0; a2 < size ; a2++ ){
			j=j+1; 
            //	printf("j = %d\n", j);
			for( a3 = 0 ; a3 < size ; a3++ ){
#ifdef TAUPROF
                TAU_START( "b" );
                TAU_PROFILER_START( ptr );
#endif
				A = T[a1][a2][a3];
                /* Beyond this point, a2 and a3 are only used in the simplectic matrix */
				for( b1 = 0 ; b1 < size ; b1++ ){
                    Tn = one_level1_product( &T, &J, A, size, a1, a2, a3, b1 );
                    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" );
#endif // TAUPROF                
			}
		}
	}
    
    return Tens;
}

gi::ex one_level1_product( tensor3D_t* T, matrix_int_t *J, gi::ex A, int size, int a1, int a2, int a3, int b1){
    gi::ex TAB, TABB, TABC, TABCC, TABCD, TABCDD;
    int b2, b3, c1, c2, c3, d1, d2, d3;

    gi::ex Tens = 0;
    gi::ex T0, T1, T2, T3, T4, T5;
    const char timerB[] = "B";

    TAB = (*J)[a1][b1]; 
    for( b2 = 0 ; b2 < size ; b2++ ){
        for( b3 = 0 ; b3 < size ; b3++ ){
            TABB =  TAB * A*(*T)[b1][b2][b3];
            T5 = 0;
            /* Beyond this point, b1 is not used anymore */
	    //            TAU_START( timerB );
            for( c1 = 0 ; c1 < size ; c1++ ){
                T4 = 0;
                for( c2 = 0 ; c2 < size ; c2++ ){
                    TABC = TABB * (*J)[a2][c2];
                    T3 = 0;
                    for( c3 = 0 ; c3 < size ; c3++ ){
                        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
                        T2 = 0;
                        for( d1 = 0 ; d1 < size ; d1++ ){
                            TABCD = TABCC * (*J)[c1][d1];
                            T1 = 0;
                            for( d2 = 0 ; d2 < size ; d2++ ){
                                TABCDD = TABCD * (*J)[b2][d2];
                                T0 = 0;
                                for( d3 = 0 ; d3 < size ; d3++ ){
                                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                                }
                                T1 += T0;
                            }
                            T2 += T1;
                        }
                        T3 += T2;
                    }
                    T4 += T3;
                }
                T5 += T4;
            }
            Tens += T5;
	    //            TAU_STOP( timerB );
        }
    }

#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, matrix_int_t *J, gi::ex A, int size, int a1, int a2, int a3, int b1, int b2, int b3, int c1, int c2 ){
    gi::ex TAB, TABB, TABC, TABCC, TABCD, TABCDD;
    int c3, d1, d2, d3;

    gi::ex Tens = 0;
    gi::ex T0, T1, T2;
    
    TAB = (*J)[a1][b1]; 
    TABB =  TAB * A*(*T)[b1][b2][b3];
    TABC = TABB * (*J)[a2][c2];
    for( c3 = 0 ; c3 < size ; c3++ ){
        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
        T2 = 0;
        for( d1 = 0 ; d1 < size ; d1++ ){
            TABCD = TABCC * (*J)[c1][d1];
            T1 = 0;
            for( d2 = 0 ; d2 < size ; d2++ ){
                TABCDD = TABCD * (*J)[b2][d2];
                T0 = 0;
                for( d3 = 0 ; d3 < size ; d3++ ){
                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                }
                T1 += T0;
            }
            T2 += T1;
        }
        Tens += T2;
    }
    
    return Tens;
}

/*******************************************************************************
 *                             2-level decomposition                           *
 *******************************************************************************/

gi::ex multiply_2levels( tensor3D_t& T, matrix_int_t& J, int size ) {  // simpler: same dimension everywhere
    
    gi::ex Tens = 0;
	int a1, a2, a3, b1;
    gi::ex A;
    
 	int i, j;
    i = 0;
    j = 0;
    
    for( a1 = 0 ; a1 < size; a1++ ){
        i=i+1; 
		for( a2 = 0; a2 < size ; a2++ ){
			j=j+1; 
			for( a3 = 0 ; a3 < size ; a3++ ){
				A = T[a1][a2][a3];
                /* Beyond this point, a2 and a3 are only used in the simplectic matrix */
				for( b1 = 0 ; b1 < size ; b1++ ){
                    Tens += two_level1_product( &T, &J, A, size, a1, a2, a3, b1 );
                }
			}
		}
	}
    return Tens;
}

gi::ex two_level1_product( tensor3D_t* T, matrix_int_t *J, gi::ex A, int size, int a1, int a2, int a3, int b1){
    gi::ex TAB, TABB, TABC;
    int b2, b3, c1, c2;

    gi::ex Tens = 0;
    gi::ex T0, T1, T2;

    TAB = (*J)[a1][b1]; 
    for( b2 = 0 ; b2 < size ; b2++ ){
        T2 = 0;
        for( b3 = 0 ; b3 < size ; b3++ ){
            TABB =  TAB * A*(*T)[b1][b2][b3];
            T1 = 0;
            /* Beyond this point, b1 is not used anymore */
            for( c1 = 0 ; c1 < size ; c1++ ){
                T0 = 0;
                for( c2 = 0 ; c2 < size ; c2++ ){
                    TABC = TABB * (*J)[a2][c2];
                    T0 += two_level2_product( T, J, TABC, size, a3, b2, b3, c1, c2 );
                }
                T1 += T0;
            }
            T2 += T1;
        }
        Tens += T2;
    }
    return Tens;
}

gi::ex two_level2_product( tensor3D_t* T, matrix_int_t *J, gi::ex TABC, int size, int a3, int b2, int b3, int c1, int c2 ){

    int c3, d1, d2, d3;
    gi::ex TABCC, TABCD, TABCDD;
    
    gi::ex Tens = 0;
    gi::ex T0, T1, T2;

    for( c3 = 0 ; c3 < size ; c3++ ){
        TABCC = TABC * (*T)[c1][c2][c3] * (*J)[b3][c3] ;
        T2 = 0;
        for( d1 = 0 ; d1 < size ; d1++ ){
            TABCD = TABCC * (*J)[c1][d1];
            T1 = 0;
            for( d2 = 0 ; d2 < size ; d2++ ){
                TABCDD = TABCD * (*J)[b2][d2];
                T0 = 0;
                for( d3 = 0 ; d3 < size ; d3++ ){
                    T0 += TABCDD * (*T)[d1][d2][d3]*(*J)[a3][d3];
                }
                T1 += T0;
            }
            T2 += T1;
        }
        Tens += T2;
    } 
    return Tens;
}