Commit c81edcf3 authored by Fabio Roberto Vitello's avatar Fabio Roberto Vitello
Browse files

Update aprod.c

parent 8e801e93
Loading
Loading
Loading
Loading
+15 −86
Original line number Diff line number Diff line
@@ -13,15 +13,11 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
           double *systemMatrix, long int *matrixIndex, int *instrCol, int *instrConstrIlung, struct comData comlsqr, time_t *ompSec)
{
   

    
    // Parallel definitions
    int myid, nproc;
    long int *mapNoss, *mapNcoeff;
    int nthreads, tid, ntasks;
    long **mapForThread;
    ///    struct comData *comlsqr;
    //

    FILE *fk, *fk0;

@@ -77,18 +73,8 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
    long offsetGlobParam = comlsqr.offsetGlobParam;


   


    
    /* for(int n=0;n<ntasks;n++)
        comlsqr->mapForThread[n]=(long *) calloc(3,sizeof(long));
*/
    //nthreads = 1;
    tid = 0;
    FILE *fp1, *fp2;
    //        fp1=fopen("test1_aprod","w");
    //        fp2=fopen("test2_aprod","w");


    if (mode != 1 && mode != 2)
    {
@@ -103,7 +89,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
        
        
        time_t startTime = time(NULL);
////        #pragma omp parallel private(myid, sum, k, l1, l2, l, j, tid, nthreads, i2, na) shared(mapNoss, instrCol, comlsqr, vVect, systemMatrix, matrixIndex, knownTerms, j2)
        {
            myid = comlsqr.myid;
            
@@ -128,8 +113,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
            int nGlobVal = nAstroPSolved + nAttP + nInstrPSolved;
            jstartAstro = miValAstro - offLocalAstro;
          
            //FV_ EDIT ompSs
          
            
            for(int nt=0; nt < ntasks; nt++ )
            {
@@ -221,8 +204,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                {
                    sum = 0.0;
                    offExtConstr = mapNcoeff[myid] + iexc * nOfElextObs;
                    //FV_ EDIT ompSs
                    //#pragma omp for
                    for (int j3 = 0; j3 < numOfExtStar * nAstroPSolved; j3++)
                        sum += systemMatrix[offExtConstr + j3] * vVect[j3];

@@ -230,16 +211,13 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                    {
                        offExtAtt = offExtConstr + numOfExtStar * nAstroPSolved + nax * numOfExtAttCol;
                        vVIx = offExtAttConstr + nax * nDegFreedomAtt;
                        //FV_ EDIT ompSs
                        //#pragma omp for
                        for (int j3 = 0; j3 < numOfExtAttCol; j3++)
                            sum += systemMatrix[offExtAtt + j3] * vVect[vVIx + j3];
                    }
                    //FV_ EDIT ompSs
                    //#pragma omp atomic
                    knownTerms[ktIx + iexc] += sum;
                } //for iexc
            }

            //////////////////////////////////////////////////////
            /// Mode 1 BarConstr
            if (nEqBarConstr)
@@ -251,12 +229,8 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                {
                    sum = 0.0;
                    offBarConstrIx = offBarConstr + iexc * nOfElBarObs;
                    //FV_ EDIT ompSs
                    //#pragma omp for
                    for (int j3 = 0; j3 < numOfBarStar * nAstroPSolved; j3++)
                        sum += systemMatrix[offBarConstrIx + j3] * vVect[j3];
                    //FV_ EDIT ompSs
                    //#pragma omp atomic
                    knownTerms[ktIx + iexc] += sum;
                } //for iexc
            }
@@ -283,15 +257,11 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                        offSetInstr += instrConstrIlung[m];
                    }
                    offvV = mapNoss[myid] * nInstrPSolved + offSetInstr;
                    //FV_ EDIT ompSs
                    //#pragma omp for
                    for (int j3 = 0; j3 < instrConstrIlung[i1]; j3++)
                    {
                        vVix = instrCol[offvV + j3];
                        sum += systemMatrix[offSetInstrInc + j3] * vVect[offSetInstrConstr1 + vVix];
                    }
                    //FV_ EDIT ompSs
                    //#pragma omp atomic
                    knownTerms[ktIx + i1] += sum;
                }
            }
@@ -308,7 +278,7 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
    }
    
    else //mode==2
    {    //if(mode
    {    
        
        time_t startTime = time(NULL);
      
@@ -316,19 +286,12 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
        myid = comlsqr.myid;
        printf("vVect[1]: %f:\n", vVect[1] );

            //qui lavorare
            /*
            #ifdef OMP
            tid = omp_get_thread_num();
            nthreads = omp_get_num_threads();
            #endif
            */
         
            /////////////////////////////////////////////////////
            /// Mode 2 Astrometric Sect
            if (nAstroPSolved)
            {
                long offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved;
                //                for (long ix = mapForThread[tid][0]; ix < mapForThread[tid][2]; ix++)
                long lset;
                long jstartAstro;
                               
@@ -373,10 +336,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                            aprodM2AstroP(ix, nt, comlsqr, vVect, systemMatrix, matrixIndex, knownTerms);
                        } //for ix 
                    }
                  //  #pragma omp taskwait

                    

                }
        
              
@@ -395,12 +354,9 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
            long lset;


          //  printf("ntask: %d\n",ntasks);
            
            for(int nt=0; nt < ntasks; nt++ )
            {
                // #pragma omp task shared(nparam,nInstrPSolved,matrixIndex) in (nt) commutative(vVect)
//                #pragma omp task commutative(vVect)
                #pragma omp task label(nAttP) commutative(vVect)//out(vVect[mix : mix + (nAttAxes-1) * nDegFreedomAtt])
                {
                    for (long ix = mapForThread[nt][0]; ix < mapForThread[nt][2]; ix++)
@@ -436,11 +392,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,

            }
       
           // long lset;
            ///#pragma omp for
            
            
       
        }
        /////////////////////////////////////////////////////
        /// Mode 2 Instrument Sect
@@ -457,8 +408,6 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,

            for(int nt=0; nt < ntasks; nt++ )
            {
                // #pragma omp task shared(vVect,nparam,offLset,nInstrPSolved) in (nt)
               // #pragma omp task commutative(vVect)
              #pragma omp task label(nOfInstrConstr) commutative(vVect) //out(vVect[mix : mix + (nAttAxes-1) * nDegFreedomAtt])
                {
                    long lset;
@@ -487,15 +436,13 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
                        aprodM2InstrP(ix, nt, comlsqr, vVect, systemMatrix, instrCol, knownTerms);
                    }     //for ix
                }
              //  #pragma omp taskwait
            }
           
        }
        /////////////////////////////////////////////////////
        ///} //pragma
        *ompSec += time(NULL) - startTime;
        if (comlsqr.itn <= 2 && (myid == 0 || myid == nproc - 1 || debugMode == 1))
            printf("PE=%d AprodTiming: mode=2.1  OmpSec=%ld\n", myid, time(NULL) - startTime);
        
        /////////////////////////////////////////////////////
        /// Mode 2 Global Sect
        if (nGlobP )
@@ -533,29 +480,15 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
        }
        ////////////////////////////////////////////////////
        startTime = time(NULL);       
        //FV_ EDIT ompSs
        //  #pragma omp parallel private(myid, yi, localSum, tid, nthreads, i2, j2, na) shared(mapNoss, vVect, systemMatrix, knownTerms, k, j3)
       
        
        {
        myid = comlsqr.myid;
            //FV_ EDIT ompSs
            /*
            //#ifdef OMP
            tid = omp_get_thread_num();
            nthreads = omp_get_num_threads();
            #endif
*/
            localSum = 0.0;
            
           
        localSum = 0.0;
        printf("\t\t\t\nEqExtConstr: %d nEqBarConstr: %d nOfInstrConstr: %d", nEqExtConstr,nEqBarConstr,nOfInstrConstr);


        #pragma omp taskwait
            
            //////////////////////////////////////////////////////
            /// Mode 2 ExtConstr
        if (nEqExtConstr && true)
            {
                time_t startTime_nEqExtConstr = time(NULL);
@@ -679,13 +612,9 @@ void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms,
            }
            ////////////

        } //pragma
        *ompSec += time(NULL) - startTime;
        if (comlsqr.itn <= 2 && (myid == 0 || myid == nproc - 1 || debugMode == 1))
            printf("PE=%d AprodTiming: mode=2.2  OmpSec=%ld\n", myid, time(NULL) - startTime);
    } // else if(mode==2)
    
    
    
  
}