aprod.c

#ifdef OMP
#include <omp.h>
#endif
#include <stdlib.h>
#include <stdio.h>
#include "mpi.h"
//#include "pardef.h"
#include "util.h"

// NOTA: vedere come e cosa cambia aumentando i valori nAstroPSolved, nAttP e nInstrConst facendoli non usuali ma molto grossi

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;

    double zero = 0.0;
    double sum, yi;
    long int l1, l2;
    long int i, i1;
    long int l, j, k;
    int i2 = 0, j2 = 0, j3 = 0, na = 0;
    int setBound[4];

    double localSum;
    int nAstroPSolved = comlsqr.nAstroPSolved;

    long localAstro = comlsqr.VrIdAstroPDim * nAstroPSolved;
    long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved;
    //  Initialize.
    myid = comlsqr.myid;
    nproc = comlsqr.nproc;
    mapNcoeff = comlsqr.mapNcoeff;
    mapNoss = comlsqr.mapNoss;
    nthreads = comlsqr.nthreads;
    ntasks= comlsqr.ntasks;
    
    mapForThread = comlsqr.mapForThread;
    int multMI = comlsqr.multMI;
    long nparam = comlsqr.parOss;
    short nAttAxes = comlsqr.nAttAxes;
    int numOfExtStar = comlsqr.numOfExtStar;
    int numOfBarStar = comlsqr.numOfBarStar;
    int numOfExtAttCol = comlsqr.numOfExtAttCol;
    long VrIdAstroPDimMax = comlsqr.VrIdAstroPDimMax;
    int startingAttColExtConstr = comlsqr.startingAttColExtConstr;
    int nOfElextObs = comlsqr.nOfElextObs;
    int nEqExtConstr = comlsqr.nEqExtConstr;
    int nOfElBarObs = comlsqr.nOfElBarObs;
    int nEqBarConstr = comlsqr.nEqBarConstr;
    int debugMode = comlsqr.debugMode;
    int nInstrPSolved = comlsqr.nInstrPSolved;
    int nOfInstrConstr = comlsqr.nOfInstrConstr;
    int nElemIC = comlsqr.nElemIC;
    int nAttP = comlsqr.nAttP;
    int nGlobP = comlsqr.nGlobP;

    setBound[0] = comlsqr.setBound[0];
    setBound[1] = comlsqr.setBound[1];
    setBound[2] = comlsqr.setBound[2];
    setBound[3] = comlsqr.setBound[3];
    long nDegFreedomAtt = comlsqr.nDegFreedomAtt;
    int nAttParAxis = comlsqr.nAttParAxis;
    long offsetAttParam = comlsqr.offsetAttParam;
    long offsetInstrParam = comlsqr.offsetInstrParam;
    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)
    {
        printf("ERROR: Invalid mode=%d in aprod function\n", mode);
        exit(1);
    }
    l1 = 0;
    l2 = 0;
    myid = comlsqr.myid;
    if (mode == 1)
    {
        
        
        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;
            
            
            if (comlsqr.itn == 1 && debugMode)
                printf("PE=%d Aprod1 OpenMP num of threads =%d from thread =%d icycle=%ld  comlsqr.itn=%d\n", myid, nthreads, tid, i, comlsqr.itn);
            long miValAstro = 0;
            long miValAtt = 0;
            long jstartAtt = 0;
            long jstartAstro = 0;
            long lset = 0;
            long offLocalAstro = 0;
            long offLocalAtt = 0;
            long offLocalInstr = 0; //Offset on Instruments
            long ixInstr = 0;
            int nInstrVal = 0;
            offLocalInstr = offsetInstrParam + (localAstroMax - offsetAttParam); //Offset on Instruments
            nInstrVal = nAstroPSolved + nAttP;
            offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved;               //Offset on my mstars
            offLocalAtt = localAstroMax - offsetAttParam;                           //Offset on attitude
            long offLocalGlob = offsetGlobParam + (localAstroMax - offsetAttParam); //Offset on GlobP
            int nGlobVal = nAstroPSolved + nAttP + nInstrPSolved;
            jstartAstro = miValAstro - offLocalAstro;

            //FV_ EDIT ompSs
          
            
            for(int nt=0; nt < ntasks; nt++ )
            {
                #pragma omp task label(mode1)
                {
                    for (long ix = mapForThread[nt][0]; ix < mapForThread[nt][2]; ix++)  //FUNZIONE DA FARE
                    {

/*================= INIZIO
                        sum = 0.;   // FARE UNA FUNZIONE CHIAMATA NEL CICLO FOR
                        /////////////////////////////////////////////////////
                        /// Mode 1 Astrometric Sect
                        if (nAstroPSolved)
                        {

                            lset = ix * nparam;
                            if (matrixIndex[multMI * ix] != miValAstro)
                            {
                                miValAstro = matrixIndex[multMI * ix];
                                jstartAstro = miValAstro - offLocalAstro;
                            }
                            for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++)
                            {
                                sum = sum + systemMatrix[lset] * vVect[jx];
                                lset++;
                            }
                        }
                        //////////////////////////////////////////////////////
                        /// Mode 1 Attitude Sect
                        if (nAttP)
                        {
                            lset = ix * nparam + nAstroPSolved;
                            miValAtt = matrixIndex[multMI * ix + (multMI - 1)];
                            for (int nax = 0; nax < nAttAxes; nax++)
                            {
                                jstartAtt = miValAtt + offLocalAtt + nax * nDegFreedomAtt;
                                for (long inpax = jstartAtt; inpax < jstartAtt + nAttParAxis; inpax++)
                                {
                                    sum = sum + systemMatrix[lset] * vVect[inpax];
                                    lset++;
                                }
                            }
                        }
                        //////////////////////////////////////////////////////
                        /// Mode 1 Instrument Sect
                        if (nInstrPSolved)
                        {

                            lset = ix * nparam + nInstrVal;
                            long iiVal = ix * nInstrPSolved;
                            for (int inInstr = 0; inInstr < nInstrPSolved; inInstr++)
                            {
                                ixInstr = offLocalInstr + instrCol[iiVal + inInstr];
                                sum = sum + systemMatrix[lset] * vVect[ixInstr];
                                lset++;
                            }
                        }
                        //////////////////////////////////////////////////////
                        /// Mode 1 Global sect
                        if (nGlobP)
                        {
                            lset = ix * nparam + nGlobVal;
                            for (long inGlob = offLocalGlob; inGlob < offLocalGlob + nGlobP; inGlob++)
                            {
                                sum = sum + systemMatrix[lset] * vVect[inGlob];
                                lset++;
                            }
                        }