solvergaiaSim.c

        printf("ssInstrFlag= %7d\n", ssInstrFlag);
        
        if(ssInstrFlag<0 || ssInstrFlag >1){
            printf("Execution aborted  with invalid ssInstrFlag\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"ssInstrFlag= %hd\n", ssInstrFlag);
       if(fscanf(fpSolProps, "nuInstrFlag= %d\n",&nuInstrFlag) != 1)
        {
            printf("Error reading nuInstrFlag  %d \n",nuInstrFlag);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nuInstrFlag= %7d\n", nuInstrFlag);
        
        if(nuInstrFlag<0 || nuInstrFlag >1){
            printf("Execution aborted  with invalid lsInstrFlag\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nuInstrFlag= %hd\n", nuInstrFlag);
        if(fscanf(fpSolProps, "maInstrFlag= %d\n",&maInstrFlag) != 1)
        {
            printf("Error reading maInstrFlag  %d \n",maInstrFlag);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("maInstrFlag= %7d\n", maInstrFlag);
        
        if(maInstrFlag<0 || maInstrFlag >1){
            printf("Execution aborted  with invalid maInstrFlag\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"maInstrFlag= %hd\n", maInstrFlag);

        if(nInstrPSolved != maInstrFlag+nuInstrFlag+ssInstrFlag+3*lsInstrFlag){
            printf("Execution aborted   invalid nInstrPSolved=%d. It should be =%d\n",nInstrPSolved,maInstrFlag+nuInstrFlag+ssInstrFlag+3*lsInstrFlag);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
           
        }
        
        if(fscanf(fpSolProps, "nOfInstrConstr= %d\n",&nOfInstrConstr) != 1)
        {
            printf("Error reading nOfInstrConstr  %d \n",nOfInstrConstr);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nOfInstrConstr= %7d\n", nOfInstrConstr);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nOfInstrConstr= %d\n", nOfInstrConstr);
        
        if(fscanf(fpSolProps, "nElemIC= %d\n",&nElemIC) != 1)
        {
            printf("Error reading nElemIC  %d \n",nElemIC);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nElemIC= %7d\n", nElemIC);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nElemIC= %d\n", nElemIC);
       
        
        if(fscanf(fpSolProps, "nGlobP= %hd\n",&nGlobP) != 1) {
            printf("Error reading nGlobP  %hd \n",nGlobP);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroGlob) nGlobP=0;
        printf("nGlobP= %7hd\n", nGlobP);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nGlobP= %hd\n", nGlobP);
        
        if(fscanf(fpSolProps, "nObs= %ld\n",&nobs) != 1) {
            printf("Error reading nObs  %ld \n",nobs);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nObs= %10ld\n", nobs);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nObs= %ld\n", nobs);
        
        fclose(fpSolProps);
        if(wrFilebin) fclose(fpSolPropsFileBin);
      }//if (autoRun) else
    }
    
    endTime=MPI_Wtime();
    if((nDegFreedomAtt==0 && nAttAxes>0) || (nDegFreedomAtt>0 && nAttAxes==0) ){
        if(myid==0){
            printf("inconsistent values for nDegFreedomAtt=%ld and nAttaxes=%d\n",nDegFreedomAtt,nAttAxes);
        MPI_Abort(MPI_COMM_WORLD, 1);
        exit(EXIT_FAILURE);
        }
    }
 
    if(myid==0) printf("Time to read initial parameters =%f sec.\n",endTime-startTime);
    // Start section for parameter broadcast
    MPI_Bcast( &atol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( &btol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( &conlim, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( &damp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nAstroP, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nAstroPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nInstrP, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nInstrPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &lsInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &ssInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nuInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &maInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nElemIC, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nOfInstrConstr, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nGlobP, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &itnlim, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nStar, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nDegFreedomAtt, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nAttAxes, 1, MPI_SHORT, 0, MPI_COMM_WORLD);
    MPI_Bcast( &instrSetUp, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( instrConst, DEFAULT_NINSTRINDEXES , MPI_INT, 0, MPI_COMM_WORLD); //  errore
    MPI_Bcast( &nobs, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nConstrLat, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nConstrLong, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nConstrMuLat, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( &nConstrMuLong, 1, MPI_LONG, 0, MPI_COMM_WORLD);
    if(myid !=0)
    {
        constrLatId=(long *) calloc(nConstrLat , sizeof(long));
        constrLatW=(double *) calloc(nConstrLat , sizeof(double));
        constrLongId=(long *) calloc(nConstrLong , sizeof(long));
        constrLongW=(double *) calloc(nConstrLong , sizeof(double));
        constrMuLatId=(long *) calloc(nConstrMuLat , sizeof(long));
        constrMuLatW=(double *) calloc(nConstrMuLat , sizeof(double));
        constrMuLongId=(long *) calloc(nConstrMuLong , sizeof(long));
        constrMuLongW=(double *) calloc(nConstrMuLong , sizeof(double));
    }
    MPI_Bcast( constrLatId, nConstrLat, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrLatW, nConstrLat, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrLongId, nConstrLong, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrLongW, nConstrLong, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrMuLatId, nConstrMuLat, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrMuLatW, nConstrMuLat, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrMuLongId, nConstrMuLong, MPI_LONG, 0, MPI_COMM_WORLD);
    MPI_Bcast( constrMuLongW, nConstrMuLong, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    
    if(nInstrPSolved==0)
        zeroInstr=1;
    /////////// Multiplicity of matrixIndex
    
    int multMI;
    
    if(nAstroPSolved)
        multMI=2;
    else{
        multMI=1;
        extConstraint=0;
        barConstraint=0;
    }
    nAstroParam = nStar * nAstroPSolved;
    nAttParam = nDegFreedomAtt * nAttAxes;
    if(nDegFreedomAtt<4) nAttParAxis=nDegFreedomAtt;
    if(nDegFreedomAtt) nAttP = nAttAxes * nAttParAxis;
    long nFoVs=1+instrConst[0];
    long nCCDs=instrConst[1];
    long nPixelColumns=instrConst[2];
    long nTimeIntervals=instrConst[3];
    // tot. instr. param. = nCCDs (Cmag) + nFoVs*nCCDs (Cnu) + nCCDs*nPixelColumns (delta_eta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_eta) + nCCDs*nPixelColumns (delta_zeta) + 3*nFoVs*nCCDs*nTimeIntervals (Delta_zeta)
    // added flags switch on and off the appropriate kind of parameters
    nInstrParam = maInstrFlag*nCCDs+nuInstrFlag*nFoVs*nCCDs+ssInstrFlag*2*nCCDs*nPixelColumns+lsInstrFlag*2*3*nFoVs*nCCDs*nTimeIntervals;
    nInstrParamTot =  nCCDs+ nFoVs*nCCDs+ 2*nCCDs*nPixelColumns+2*3*nFoVs*nCCDs*nTimeIntervals;
    nGlobalParam = nGlobP;
    
    if(nElemIC<0 || nOfInstrConstr <0){
        
        
        nOfInstrConstrLSAL = lsInstrFlag*nTimeIntervals;
        nElemICLSAL = nFoVs*nCCDs;
        nOfInstrConstrLSAC = lsInstrFlag*nFoVs*nTimeIntervals;
        nElemICLSAC = nCCDs;
        nOfInstrConstrSS = lsInstrFlag*ssInstrFlag*2*nCCDs*3; // the factor 2 comes from the AL and AC constraint equations
        nElemICSS = nPixelColumns;
        nOfInstrConstr = nOfInstrConstrLSAL + nOfInstrConstrLSAC + nOfInstrConstrSS;
        nElemIC = nOfInstrConstrLSAL*nElemICLSAL + nOfInstrConstrLSAC*nElemICLSAC + nOfInstrConstrSS*nElemICSS;
        
    }

    
    // Map the solved astrometric parameters, putting their indeces into an array of size nAstroPSolved
    // astroCoeff[0] -> parallax
    // astroCoeff[1] -> alpha
    // astroCoeff[2] -> delta
    // astroCoeff[3] -> mu alpha
    // astroCoeff[4] -> mu delta
    // Therefore, nAstroPSolved=2 => alpha, delta, and the array is mapAstroP=[1,2]
    //            nAstroPSolved=3 => parallax, alpha, delta, and the array is mapAstroP=[0,1,2]
    //            nAstroPSolved=4 => alpha, delta, mu alpha, mu delta and the array is mapAstroP=[1,2,3,4]
    //            nAstroPSolved=5 => parallax, alpha, delta, mu alpha, mu delta and the array is mapAstroP=[0,1,2,3,4]
    
    int LatPos=-1, LongPos=-1,MuLatPos=-1, MuLongPos=-1;
   
    if(nAstroPSolved)
    {
        mapAstroP=(int *) calloc(nAstroPSolved, sizeof(int));
        switch(nAstroPSolved) {
            case 2:
                mapAstroP[0] = 1;
                mapAstroP[1] = 2;
                LongPos=0;
                LatPos=1;
                MuLongPos=-1;
                MuLatPos=-1;
                nConstrMuLat=0;
                nConstrMuLong=0;
                if(extConstraint) nEqExtConstr=3;
                if(barConstraint) nEqBarConstr=3;
                break;
            case 3:
                mapAstroP[0] = 0;
                mapAstroP[1] = 1;
                mapAstroP[2] = 2;
                LongPos=1;
                LatPos=2;
                MuLongPos=-1;
                MuLatPos=-1;
                nConstrMuLat=0;
                nConstrMuLong=0;
                if(extConstraint) nEqExtConstr=3;
                if(barConstraint) nEqBarConstr=3;
               break;
            case 4:
                mapAstroP[0] = 1;
                mapAstroP[1] = 2;
                mapAstroP[2] = 3;
                mapAstroP[3] = 4;
                LongPos=0;
                LatPos=1;
                MuLongPos=2;
                MuLatPos=3;
                break;
            case 5:
                mapAstroP[0] = 0;
                mapAstroP[1] = 1;
                mapAstroP[2] = 2;
                mapAstroP[3] = 3;
                mapAstroP[4] = 4;
                LongPos=1;
                LatPos=2;
                MuLongPos=3;
                MuLatPos=4;
                break;
            default:
                if(myid==0) {
                    printf("nAstroPSolved=%d, invalid value. Aborting.\n", nAstroPSolved);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
        }
    }
    // End map
    

    /////////////////////// end broadcast section
    
    // Initialize the values needed to dimension the vectors
    nConstr = nConstrLong + nConstrLat + nConstrMuLong + nConstrMuLat; // number of constraints to be generated
    
    nObsxStar=nobs/nStar;
    nobsOri=nobs;
    if(noConstr)
    {
        nConstr=0;
        nConstrLong=0;
        nConstrLat=0;
        nConstrMuLong=0;
        nConstrMuLat=0;
    } else
    {
        for(int q=0;q<nConstrLat;q++)
            if(constrLatId[q]>=nStar)
            {
                printf("PE=%d Error invalit Lat Constraint %ld\n",myid,constrLatId[q]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        
        for(int q=0;q<nConstrLong;q++)
            if(constrLongId[q]>=nStar)
            {
                printf("PE=%d Error invalit Long Constraint %ld\n",myid,constrLongId[q]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        for(int q=0;q<nConstrMuLat;q++)
            if(constrMuLatId[q]>=nStar)
            {
                printf("PE=%d Error invalit MuLat Constraint %ld\n",myid,constrMuLatId[q]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        
        for(int q=0;q<nConstrMuLong;q++)
            if(constrMuLongId[q]>=nStar)
            {
                printf("PE=%d Error invalit MuLat Constraint %ld\n",myid,constrMuLongId[q]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        
        nConstr = nConstrLong + nConstrLat + nConstrMuLong + nConstrMuLat; // number of constraints to be generated
    }
    
    
    nobs+=nConstr;
    nunk = ((long) nAstroParam) + nAttParam + nInstrParam + nGlobalParam; // number of unknowns (i.e. columns of the system matrix)
    nparam = nAstroPSolved + nAttP + nInstrPSolved + nGlobP; // number of non-zero coefficients for each observation (i.e. for each system row)
    if(nparam==0){
        printf("Abort. Empty system nparam=0 . nAstroPSolved=%d nAttP=%d nInstrPSolved=%d nGlobP=%d\n",nAstroPSolved,nAttP,nInstrPSolved,nGlobP);
        MPI_Abort(MPI_COMM_WORLD,1);
    }
    ncoeff = nparam * nobs; // total number of non-zero coefficients of the system
    
    if(nobs <=nunk){
        printf("SEVERE ERROR: number of equations=%ld and number of unknown=%ld make solution unsafe\n",nobs,nunk);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);

    }
    
    // Start map distributed array
    mapNoss=(long int *) calloc(nproc,sizeof(long int));
    mapNcoeff=(long int *) calloc(nproc,sizeof(long int));
    mapNossAfter=0;
    mapNcoeffAfter=0;
    mapNossBefore=0;
    mapNcoeffBefore=0;
    for(i=0;i<nproc;i++)
    {
        mapNoss[i]=(nobs)/nproc;
        if(nobs % nproc >=i+1) mapNoss[i]++;
        mapNcoeff[i]=mapNoss[i]*nparam;
        if(i<myid)
        {
            mapNossBefore+=mapNoss[i];
            mapNcoeffBefore+=mapNcoeff[i];
        }
        if(i>myid)
        {
            mapNossAfter+=mapNoss[i];
            mapNcoeffAfter+=mapNcoeff[i];
        }
    }
    ////////////////// Simulating the ... of NObsxStar file
    if(extConstraint){
        int * sumNObsxStar;
        sumNObsxStar=(int *) calloc(nStar, sizeof(int));
        int irest=nobs % nStar;
        for(int i=0;i<nStar;i++){
            sumNObsxStar[i]=nobs/nStar;
            if(i<irest) sumNObsxStar[i]++;
        }
        if(wrFilebin && myid==0){
            chdir(wpath);
            chdir(wrfileDir);
            FILE *fpNObsxStar;
            fpNObsxStar=fopen("NObsStar.bin","wb");
            fwrite(sumNObsxStar,sizeof(int),nStar,fpNObsxStar);
            fclose(fpNObsxStar);
            chdir(wpath);
            
        }
        long counterObsxStar=0;
        for(int i=0;i<nStar;i++){
            counterObsxStar+=sumNObsxStar[i];
            if(counterObsxStar>mapNossBefore && firstStarConstr==-1) firstStarConstr=i;  //first star assigned in  extConstr
            if(counterObsxStar>=mapNossBefore+mapNoss[myid] && lastStarConstr==-1){
                lastStarConstr=i; //last star assigned in  extConstr (it will be eqaul to lastrStar-1 in case of overlap)
                if(counterObsxStar>(mapNossBefore+mapNoss[myid]) && myid!=(nproc-1)){
                    starOverlap=1;
                    lastStarConstr--;
                }
                break;
            }
        }
        numOfExtStar=lastStarConstr-firstStarConstr+1;  //number of stars computed in ext Constr
        
        int counterAttCols=0;
        startingAttColExtConstr=0;  // numero di colonna di assetto escluso l'offset: la prima è 0 per il PE0 x asse
        endingAttColExtConstr=0; // numero di colonna di assetto finale l'offset: la prima è nDegFreedomAtt/nproc-1 (+1 in caso di modulo) per il PE0 x asse
        int attRes=nDegFreedomAtt%nproc;
        startingAttColExtConstr=(nDegFreedomAtt/nproc)*myid;
        if(myid<attRes)startingAttColExtConstr+=myid;
        else startingAttColExtConstr+=attRes;
        endingAttColExtConstr=startingAttColExtConstr+(nDegFreedomAtt/nproc)-1;
        if(myid<attRes)endingAttColExtConstr++;
        
        
        numOfExtAttCol=endingAttColExtConstr-startingAttColExtConstr+1;  //numeroi di colonne x asse
    }
    //////////////////// barConstraint
    if(barConstraint){
        int * sumNObsxStar;
        sumNObsxStar=(int *) calloc(nStar, sizeof(int));
        int irest=nobs % nStar;
        for(int i=0;i<nStar;i++){
            sumNObsxStar[i]=nobs/nStar;
            if(i<irest) sumNObsxStar[i]++;
        }
        if(wrFilebin && myid==0){
            chdir(wpath);
            chdir(wrfileDir);
            FILE *fpNObsxStar;
            fpNObsxStar=fopen("NObsStar.bin","wb");
            fwrite(sumNObsxStar,sizeof(int),nStar,fpNObsxStar);
            fclose(fpNObsxStar);
            chdir(wpath);
            
        }
        int counterObsxStar=0;
        for(int i=0;i<nStar;i++){
            counterObsxStar+=sumNObsxStar[i];
            if(counterObsxStar>mapNossBefore && firstStarConstr==-1) firstStarConstr=i;  //first star assigned in  barConstr
            if(counterObsxStar>=mapNossBefore+mapNoss[myid] && lastStarConstr==-1){
                lastStarConstr=i; //last star assigned in  barConstr (it will be eqaul to lastrStar-1 in case of overlap)
                if(counterObsxStar>(mapNossBefore+mapNoss[myid]) && myid!=(nproc-1)){
                    starOverlap=1;
                    lastStarConstr--;
                }
                break;
            }
        }
        numOfBarStar=lastStarConstr-firstStarConstr+1;  //number of stars computed in bar Constr
        
    }
    ////////////////////// comlsqr
    
    
    
    comlsqr.nStar=nStar;
    comlsqr.nAstroP=nAstroP;
    comlsqr.nAstroPSolved=nAstroPSolved;
    comlsqr.nAttP=nAttP;
    comlsqr.nInstrP=nInstrP;
    comlsqr.nInstrPSolved=nInstrPSolved;
    comlsqr.nGlobP=nGlobP;
    comlsqr.mapNossBefore=mapNossBefore;
    comlsqr.mapNossAfter=mapNossAfter;
    comlsqr.myid=myid;
    comlsqr.nproc=nproc;
    comlsqr.mapNoss=mapNoss;
    comlsqr.mapNcoeff=mapNcoeff;
    comlsqr.multMI=multMI;
    comlsqr.debugMode=debugMode;
    comlsqr.noCPR=noCPR;
    comlsqr.nAttParam=nAttParam;
    comlsqr.extConstraint=extConstraint;
    comlsqr.nEqExtConstr=nEqExtConstr;
    comlsqr.numOfExtStar=numOfExtStar;
    comlsqr.barConstraint=barConstraint;
    comlsqr.nEqBarConstr=nEqBarConstr;
    comlsqr.numOfBarStar=numOfBarStar;
    comlsqr.firstStarConstr=firstStarConstr;
    comlsqr.lastStarConstr=lastStarConstr;
    comlsqr.numOfExtAttCol=numOfExtAttCol;
    comlsqr.startingAttColExtConstr=startingAttColExtConstr;
    comlsqr.setBound[0]=0;
    comlsqr.setBound[1]=nAstroPSolved;
    comlsqr.setBound[2]=nAstroPSolved+nAttP;
    comlsqr.setBound[3]=nAstroPSolved+nAttP+nInstrPSolved;
    comlsqr.nDegFreedomAtt=nDegFreedomAtt;
    comlsqr.nAttParAxis=nAttParAxis;
    comlsqr.nAttAxes=nAttAxes;
    comlsqr.nobs=nobs;
    comlsqr.lsInstrFlag=lsInstrFlag;
    comlsqr.ssInstrFlag=ssInstrFlag;
    comlsqr.nuInstrFlag=nuInstrFlag;
    comlsqr.maInstrFlag=maInstrFlag;
    comlsqr.myid=myid;
    comlsqr.cCDLSAACZP=cCDLSAACZP;
    comlsqr.nOfInstrConstr=nOfInstrConstr;
    comlsqr.nElemIC=nElemIC;
    comlsqr.nElemICLSAL=nElemICLSAL;
    comlsqr.nElemICLSAC=nElemICLSAC;
    comlsqr.nElemICSS=nElemICSS;
    comlsqr.instrConst[0]=instrConst[0];
    comlsqr.instrConst[1]=instrConst[1];
    comlsqr.instrConst[2]=instrConst[2];
    comlsqr.instrConst[3]=instrConst[3];
    comlsqr.nInstrParam=nInstrParam;
    comlsqr.nGlobalParam=nGlobalParam;

    ///////////////////// end buidl map distributed arrays
    
    // Allocate the memory for the vectors and compute the total memory allocated
    totmem = 0;
    nElements = mapNcoeff[myid];
    if (extConstraint){
        addElementextStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved;
        addElementAtt=numOfExtAttCol*nAttAxes;
    }
    if (barConstraint){
        addElementbarStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved;
    }
    nOfElextObs=addElementextStar+addElementAtt;
    nOfElBarObs=addElementbarStar;
    comlsqr.nOfElextObs=nOfElextObs;
    comlsqr.nOfElBarObs=nOfElBarObs;

    nElements+=nOfElextObs*nEqExtConstr+nOfElBarObs*nEqBarConstr+nElemIC;
    
    systemMatrix = (double *) calloc(nElements,sizeof(double));
    if (!systemMatrix)
        exit(err_malloc("systemMatrix",myid));
    totmem += nElements*sizeof(double);
    
    nElements = mapNoss[myid]*multMI;
    matrixIndex = (long int *) calloc(nElements, sizeof(long int));
    if (!matrixIndex)
        exit(err_malloc("matrixIndex",myid));
    totmem += nElements* sizeof(long int);
    
    nElements = mapNoss[myid]*nInstrPSolved+nElemIC;
    instrCol = (int *) calloc(nElements, sizeof(int));
    if (!instrCol)
        exit(err_malloc("instrCol",myid));
    totmem += nElements* sizeof(int);

    nElements = nOfInstrConstr;
    instrConstrIlung = (int *) calloc(nElements, sizeof(int)); // it is the vectorr that for each observation (in this PE) save the INDEX for the values of instr
    if (!instrConstrIlung)
        exit(err_malloc("instrConstrIlung",myid));
    totmem += nElements* sizeof(int);
    
    
    nElements = mapNoss[myid];
    if(extConstraint) nElements+=nEqExtConstr;
    if(barConstraint) nElements+=nEqBarConstr;
    if(nOfInstrConstr>0) nElements+=nOfInstrConstr;
    knownTerms = (double *) calloc(nElements, sizeof(double));
    if (!knownTerms)
        exit(err_malloc("knownTerms",myid));
    totmem += nElements* sizeof(double);
    
    ielem = 0;
    offsetAttParam = nAstroParam;
    offsetInstrParam = offsetAttParam + nAttParam;
    offsetGlobParam = offsetInstrParam + nInstrParam;
    comlsqr.offsetAttParam=offsetAttParam;
    
    
    long rowsxFile=nobsOri/(nproc*nfileProc);
   
    if(withFile)
    {
        char filenameSim_SM[128]="SIM_PE_SM_";
        char filenameSim_MI[128]="SIM_PE_MI_";
        char filenameSim_II[128]="SIM_PE_II_";
        char filenameSim_KT[128]="SIM_PE_KT_";
        char varpe[32]="";
        char varrows[32]="";
        sprintf(varpe,"%d",myid);
        sprintf(varrows,"%ld",rowsxFile);
        strcat(filenameSim_SM,varrows);
        strcat(filenameSim_SM,"_");
        strcat(filenameSim_SM,varpe);
        strcat(filenameSim_SM,".bin");
        strcat(filenameSim_MI,varrows);
        strcat(filenameSim_MI,"_");
        strcat(filenameSim_MI,varpe);
        strcat(filenameSim_MI,".bin");
        strcat(filenameSim_II,varrows);
        strcat(filenameSim_II,"_");
        strcat(filenameSim_II,varpe);
        strcat(filenameSim_II,".bin");
        strcat(filenameSim_KT,varrows);
        strcat(filenameSim_KT,"_");
        strcat(filenameSim_KT,varpe);
        strcat(filenameSim_KT,".bin");
        
        double *smsim;
        smsim=(double *)calloc(rowsxFile*nparam,sizeof(double));
        fpSM=fopen(filenameSim_SM,"wb");
        fwrite(smsim,sizeof(double),rowsxFile*nparam,fpSM);
        fclose(fpSM);
        fpKT=fopen(filenameSim_KT,"wb");
        fwrite(smsim,sizeof(double),rowsxFile,fpKT);
        fclose(fpKT);
        free(smsim);
        int *intsim;
        intsim=(int *)calloc(rowsxFile*nInstrPSolved,sizeof(int));
        fpII=fopen(filenameSim_II,"wb");
        fwrite(intsim,sizeof(int),rowsxFile*nInstrPSolved,fpII);
        fclose(fpII);
        free(intsim);
        long *misim;
        misim=(long *)calloc(rowsxFile*multMI,sizeof(long));
        fpMI=fopen(filenameSim_MI,"wb");
        fwrite(misim,sizeof(long),rowsxFile*multMI,fpMI);
        fclose(fpMI);
        free(misim);
        long nrowToRead;
        if(myid==0) printf(
                           "Found %d SM data files in the directory. Reading the coefficients...\n",
                           (nfileProc)*nproc);
        timeToReadFiles=MPI_Wtime();
        
        for(ii=0;ii<nfileProc+2;ii++) {
            
            if(ii==0 && myid==0) printf("PE=%d Opening  %d  and reading %s  and associated files\n",myid,nfileProc,filenameSim_SM);
            
            fpSM=fopen(filenameSim_SM,"rb");
            if (!fpSM)
            {
                printf("PE=%d Error while open %s\n",myid,filenameSim_SM);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fpMI=fopen(filenameSim_MI,"rb");
            if (!fpMI)
            {
                printf("PE=%d Error while open %s\n",myid,filenameSim_MI);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fpII=fopen(filenameSim_II,"rb");
            if (!fpII)
            {
                printf("PE=%d Error while open %s\n",myid,filenameSim_II);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fpKT=fopen(filenameSim_KT,"rb");
            if (!fpSM)
            {
                printf("PE=%d Error while open %s\n",myid,filenameSim_KT);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            
            
            ///// Reading systemMatrix
            if(ii==0 || ii==nfileProc+1)
            {
                if(ii==0) fseek(fpSM,((rowsxFile*nparam)/2)*sizeof(double),SEEK_SET);
                fread(systemMatrix,sizeof(double),((rowsxFile*nparam)/2)*sizeof(double),fpSM);
                if(ii==0) fseek(fpMI,((rowsxFile*multMI)/2)*sizeof(long),SEEK_SET);
                fread(matrixIndex,sizeof(long),((rowsxFile*multMI)/2)*sizeof(double),fpMI);
                if(ii==0) fseek(fpII,((rowsxFile*nInstrPSolved)/2)*sizeof(long),SEEK_SET);
                fread(instrCol,sizeof(int),((rowsxFile*nInstrPSolved)/2)*sizeof(int),fpII);
                if(ii==0) fseek(fpKT,((rowsxFile)/2)*sizeof(double),SEEK_SET);
                fread(instrCol,sizeof(int),((rowsxFile)/2)*sizeof(double),fpKT);
            } else{
                fread(systemMatrix,sizeof(double),((rowsxFile*nparam))*sizeof(double),fpSM);
                fread(matrixIndex,sizeof(long),((rowsxFile*multMI))*sizeof(double),fpMI);
                fread(instrCol,sizeof(int),((rowsxFile*nInstrPSolved))*sizeof(int),fpII);
                fread(systemMatrix,sizeof(double),((rowsxFile))*sizeof(double),fpKT);
            }
            fclose(fpSM);
            fclose(fpMI);
            fclose(fpII);
            fclose(fpKT);
            
        } //for(ii)
        MPI_Barrier(MPI_COMM_WORLD);
        timeToReadFiles=MPI_Wtime()-timeToReadFiles;
        if(myid==0) printf("PE=%d time seconds=%lf TO READ %d Files\n",myid, timeToReadFiles,nfileProc);
    }// if withFile (no  external and baricentric Contraint are simulated for reading time
    timeToReadFiles=MPI_Wtime();
    
    long startingStar=0;
    long obsTotal=0;
    int residual=0; //represents the number of observation for the staring stars (if equal to zero  the number of observations for the starting  star is nObsxStar)
    /////////
    for(int p=0;p<myid;p++)
    {
        while(obsTotal<mapNoss[p])
        {
            if(residual==0)
            {
                obsTotal+=nObsxStar;
                if(startingStar<nobsOri%nStar) obsTotal++;
                if(nConstr>0)
                {
                    for(int q=0;q<nConstrLat;q++)
                        if(constrLatId[q]==startingStar) obsTotal++;
                    for(int q=0;q<nConstrLong;q++)
                        if(constrLongId[q]==startingStar) obsTotal++;
                    for(int q=0;q<nConstrMuLat;q++)
                        if(constrMuLatId[q]==startingStar) obsTotal++;
                    for(int q=0;q<nConstrMuLong;q++)
                        if(constrMuLongId[q]==startingStar) obsTotal++;
                }
            } else { //if residual
                obsTotal=residual;
                residual=0;
            }
            
            if(obsTotal<=mapNoss[p])startingStar++;
        }//while
        residual=obsTotal-mapNoss[p];
        obsTotal=0;
    } // for in p
    //////////////////////////
    if(debugMode) printf("PE=%d mapNoss[myid]=%ld starting star %ld residual=%d\n",myid,mapNoss[myid],startingStar,residual);
   
    //////////////// filling the system
    long currentStar=startingStar;
    int obsStar=residual;
    int obsStarnow;
    int numOfObslast=0;
    long startFreedom=(nDegFreedomAtt/nproc)*myid;
    long endFreedom=startFreedom+(nDegFreedomAtt/nproc)+1;
    long lastFreedom=startFreedom;
    int freedomReached=0;
    long instrStartFreedom=(nInstrParam/nproc)*myid;
    long instrEndFreedom=instrStartFreedom+(nInstrParam/nproc)+1;
    if(myid==nproc-1)
        instrEndFreedom=nInstrParam-1;
    int instrFreedomReached=0;
    int isConstraint=0;
    int instrLastFreedom=instrStartFreedom;
    srand(myid);
    if(debugMode)
        printf("PE=%d instrStartFreedom=%ld instrEndFreedom=%ld nInstrParam=%d\n",myid,instrStartFreedom,instrEndFreedom,nInstrParam);
    if(obsStar==0)
    {
        obsStar=nObsxStar;
        if(currentStar<nobsOri%nStar) obsStar++;
        if(nConstr>0)
        {
            for(int q=0;q<nConstrLat;q++)
                if(constrLatId[q]==currentStar) obsStar++;
            for(int q=0;q<nConstrLong;q++)
                if(constrLongId[q]==currentStar) obsStar++;
            for(int q=0;q<nConstrMuLat;q++)
                if(constrMuLatId[q]==currentStar) obsStar++;
            for(int q=0;q<nConstrLong;q++)
                if(constrMuLongId[q]==currentStar) obsStar++;
        }
        
    }
    obsStarnow=obsStar;
    if(nConstr>0)
    {
        for(int q=0;q<nConstrLat;q++)
            if(constrLatId[q]==currentStar) isConstraint++;
        for(int q=0;q<nConstrLong;q++)
            if(constrLongId[q]==currentStar) isConstraint++;
        for(int q=0;q<nConstrMuLat;q++)
            if(constrMuLatId[q]==currentStar) isConstraint++;
        for(int q=0;q<nConstrMuLong;q++)
            if(constrMuLongId[q]==currentStar) isConstraint++;
    }
    
    int offsetConstraint=isConstraint-obsStar; // number of constraint alredy computed in the previous PE
    if(offsetConstraint<0)offsetConstraint=0;
    
    
    int counterStarObs=0;
    rowInFile=-1;
    int changedStar=0;
    int counterConstr=0;
    /////////////////////////////////
    ///////////   RUNNING ON ALL OBSERVATIONS
    /////////////////////////////////
    for(ii=0;ii<mapNoss[myid];ii++)
    {
        rowInFile++;
        if(currentStar%1000==0 && changedStar){
            rowInFile=0;
            changedStar=0;
        }
        ///////////// generate MatrixIndex
        if(currentStar==nStar)
        {
            printf("PE=%d Severe Error in currentStar=%ld ii=%ld\n",myid,currentStar,ii);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }

        if(nAstroPSolved) matrixIndex[ii*multMI]=currentStar*nAstroPSolved;
        
        if(!freedomReached && nAstroPSolved)
        {
            if((obsStar-counterStarObs)<=isConstraint)
            {  //constraint
                matrixIndex[ii*multMI+(multMI-1)]=offsetAttParam;
                constraintFound[counterConstr][0]=currentStar/1000;
                constraintFound[counterConstr][1]=rowInFile;
                counterConstr++;
                
                if(counterConstr==MAX_CONSTR){
                    printf("PE=%d Abort increase MAX_CONSTR and recompile =%d \n",myid,counterConstr);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
            }
            else{
                if(lastFreedom>=nDegFreedomAtt-nAttParAxis) lastFreedom=nDegFreedomAtt-nAttParAxis;
                matrixIndex[ii*multMI+(multMI-1)]=offsetAttParam+lastFreedom;
                if(lastFreedom>=endFreedom || lastFreedom>=nDegFreedomAtt-nAttParAxis) freedomReached=1;
                lastFreedom+=nAttParAxis;
            }
        } else {
            lastFreedom=( ( (double)rand() ) / ( ((double)RAND_MAX) ) ) * (nDegFreedomAtt-nAttParAxis+1);
            if(lastFreedom>nDegFreedomAtt-nAttParAxis) lastFreedom=nDegFreedomAtt-nAttParAxis;
            if((obsStar-counterStarObs)<=isConstraint) //constraint
            {
                lastFreedom=0;
                constraintFound[counterConstr][0]=currentStar/1000;
                constraintFound[counterConstr][1]=rowInFile;
                counterConstr++;
            }
            matrixIndex[ii*multMI+(multMI-1)]=offsetAttParam+lastFreedom;
        }
        ///////////// generate InstrIndex
        
        if(!instrFreedomReached  && nInstrPSolved)
        {
            if((obsStar-counterStarObs)<=isConstraint)
            {  //constraint
                for(int kk=0;kk<nInstrPSolved;kk++)
                    instrCol[ii*nInstrPSolved+kk]=0;
            }
            else{
                if(instrLastFreedom>instrEndFreedom) instrLastFreedom=instrEndFreedom;
                instrCol[ii*nInstrPSolved]=instrLastFreedom;
                for(int kk=1;kk<nInstrPSolved;kk++)
                    instrCol[ii*nInstrPSolved+kk]=(((double)rand()) / (((double)RAND_MAX))) * (nInstrParam-1);
                if(instrLastFreedom==instrEndFreedom)
                    instrFreedomReached=1;
                instrLastFreedom++;
            }
        } else {
            if((obsStar-counterStarObs)<=isConstraint)
            { //constraint
                for(int kk=0;kk<nInstrPSolved;kk++)
                    instrCol[ii*nInstrPSolved+kk]=0;
            }else{
                for(int kk=0;kk<nInstrPSolved;kk++)
                    instrCol[ii*nInstrPSolved+kk]=(((double)rand() ) / ( ((double)RAND_MAX) ) ) * (nInstrParam-1);
            }
        }
        ///////////// generate systemMatrix
        if((obsStar-counterStarObs)>isConstraint )
        {
            for(int q=0;q<nAstroPSolved;q++)
                systemMatrix[ii*nparam+q]=(((double)rand())/RAND_MAX)*2 - 1.0;
            for(int q=0;q<nAttP+nInstrPSolved+nGlobP;q++)
                systemMatrix[ii*nparam+nAstroPSolved+q]=(((double)rand())/RAND_MAX)*2 - 1.0;
        } else  // I add a Constraint
        {
            for(int q=0;q<nAstroPSolved+nAttP+nInstrPSolved+nGlobP;q++)
                systemMatrix[ii*nparam+q]=0.;
            if(nAstroPSolved>0)
            {
                if(ii!=0) offsetConstraint=0;
                int foundedConstraint=(obsStar-counterStarObs)+offsetConstraint;
                int itis=0;
                for(int q=0;q<nConstrLong;q++)
                    if(constrLongId[q]==currentStar){
                        itis++;
                        if(itis==foundedConstraint)
                            systemMatrix[ii*nparam+LongPos]=constrLongW[q];
                    }
                for(int q=0;q<nConstrLat;q++)
                    if(constrLatId[q]==currentStar){
                        itis++;
                        if(itis==foundedConstraint)
                            systemMatrix[ii*nparam+LatPos]=constrLatW[q];
                    }
                for(int q=0;q<nConstrMuLong;q++)
                    if(constrMuLongId[q]==currentStar){
                        itis++;
                        if(itis==foundedConstraint)
                            systemMatrix[ii*nparam+MuLongPos]=constrMuLongW[q];
                    }
                for(int q=0;q<nConstrMuLat;q++)
                    if(constrMuLatId[q]==currentStar){
                        itis++;
                        if(itis==foundedConstraint)
                            systemMatrix[ii*nparam+MuLatPos]=constrMuLatW[q];
                    }
            }
        }
        /////////////////////////
        if((obsStar-counterStarObs)<=isConstraint )
        {
            printf("PE=%d isConstraint=%d ii=%ld matrixIndex[ii*2]=%ld matrixIndex[ii*2+1]=%ld\n",myid,isConstraint,ii,matrixIndex[ii*2],matrixIndex[ii*2+1]);
            for(int q=0;q<nparam;q++) printf("PE=%d systemMatrix[%ld]=%lf ",myid,ii*nparam+q,systemMatrix[ii*nparam+q]);
            printf("\n");
        }
        
        /////////////////// Prepare next Obs
        counterStarObs++;
        if(counterStarObs==obsStar)
        {
            if(myid==(nproc-1))
                numOfObslast=counterStarObs;
            counterStarObs=0;
            currentStar++;
            changedStar=1;
            isConstraint=0;
            obsStar=nObsxStar;
            if(currentStar<nobsOri%nStar) obsStar++;
            if(nConstr>0)
            {
                for(int q=0;q<nConstrLat;q++)
                    if(constrLatId[q]==currentStar)
                    {
                        obsStar++;
                        isConstraint++;
                    }
                for(int q=0;q<nConstrLong;q++)
                    if(constrLongId[q]==currentStar)
                    {
                        obsStar++;
                        isConstraint++;
                    }
                for(int q=0;q<nConstrMuLat;q++)
                    if(constrMuLatId[q]==currentStar)
                    {
                        obsStar++;
                        isConstraint++;
                    }
                for(int q=0;q<nConstrMuLong;q++)
                    if(constrMuLongId[q]==currentStar)
                    {
                        obsStar++;
                        isConstraint++;
                    }
                
            }
        }
        ///////////////////////////////// Filling knownTerms  -1.. 1
        if(!idtest) knownTerms[ii]=(((double) rand())/RAND_MAX)*2.0-1.0; //show idtest=1  at the beginning instead of =0
        /////////////////////////////////////////
        
        
        ///////////////////////////////////
    } // for ii=0 mapNoss[myid]
    
    if (!freedomReached && !zeroAtt)
    {
        printf("PE=%d Error ndegFreedomAtt not correctly generated\n",myid);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
    }
    if (!instrFreedomReached && !zeroInstr)
    {
        printf("PE=%d Error instrP not all generated instrLastFreedom=%d\n",myid,instrLastFreedom);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
        
    }
    ////////////////////
    ///////////////////////// generate extConstr on systemMatrix
    if(extConstraint){
        double randVal;
        double *accumulator;
        accumulator=(double *) calloc(nEqExtConstr,sizeof(double));
        attNS=(double *) calloc(nDegFreedomAtt,sizeof(double));
        if (!attNS)
            exit(err_malloc("attNS",myid));
        if(myid==0){
            for(int i=0;i<nDegFreedomAtt;i++)
                attNS[i]=(((double)rand())/RAND_MAX)*2 - 1.0;
        }
        MPI_Bcast( attNS, nDegFreedomAtt, MPI_DOUBLE, 0, MPI_COMM_WORLD);
            
        for(int j=0;j<nEqExtConstr;j++){
            for(int i=0;i<addElementextStar+addElementAtt;i++){
                randVal=(((double)rand())/RAND_MAX)*2 - 1.0;
                if(i<addElementextStar){
                    if(nAstroPSolved==3 && i%nAstroPSolved==0) randVal=0.;
                    if(nAstroPSolved==4 && i%nAstroPSolved>=2) randVal=0.;
                    if(nAstroPSolved==5 && i%nAstroPSolved==0) randVal=0.;
                    if(nAstroPSolved==5 && i%nAstroPSolved>2 && j<3) randVal=0.;
                }
                if(i>=addElementextStar){
                    if(j<3) randVal=1.0;
                    if(j==0 || j==3){
                        if(i>=addElementextStar+addElementAtt/nAttAxes) randVal=0.0;
                    }
                    if(j==1 || j==4){
                        if(i<addElementextStar+addElementAtt/nAttAxes) randVal=0.0;
                        if(i>=addElementextStar+2*addElementAtt/nAttAxes) randVal=0.0;
                    }
                    if(j==2 || j==5){
                        if(i<addElementextStar+2*addElementAtt/nAttAxes) randVal=0.0;
                    }
                    
                    
                }
                systemMatrix[mapNcoeff[myid]+i+j*nOfElextObs]=randVal*extConstrW;
                accumulator[j]+=randVal*extConstrW;

            }
            if(!idtest)
                  knownTerms[mapNoss[myid]+j]=0.;
        }// j=0
        if(idtest)
            MPI_Allreduce(accumulator, &knownTerms[mapNoss[myid]], nEqExtConstr, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        free(accumulator);