solvergaiaSim.c

    } //if(extConstr
    
    ///////////////////////// generate barConstr on systemMatrix
    if(barConstraint){
        double randVal;
        double *accumulator;
        accumulator=(double *) calloc(nEqBarConstr,sizeof(double));
        
        for(int j=0;j<nEqBarConstr;j++){
            for(int i=0;i<addElementbarStar;i++){
                randVal=(((double)rand())/RAND_MAX)*2 - 1.0;
                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.;
                systemMatrix[mapNcoeff[myid]+nEqExtConstr*nOfElextObs+i+j*nOfElBarObs]=randVal*barConstrW;
                accumulator[j]+=randVal*barConstrW;
                
            }
            if(!idtest)
            knownTerms[mapNoss[myid]+nEqExtConstr+j]=0.;
        }// j=0
        if(idtest)
        MPI_Allreduce(accumulator, &knownTerms[mapNoss[myid]+nEqExtConstr], nEqBarConstr, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        free(accumulator);
    } //if(barConstr
  
    ///////////////////////// generate instrConstr on systemMatrix
    // There are 1 AL + 2 AC constraint equations for each time interval for the large scale parameters (total=3*nTimeIntervals)
    // There are 1 AL + 1 AC constraint equations for each CCD and Legendre polinomial degree (total=2*nCCDs*3)
    // The equations are described by three arrays: instrCoeffConstr, instrColsConstr, instrConstrIlung, which contain the
    // coefficients, the column indexes and the length of the non-zero coefficients of each equation respectively.

    comlsqr.offsetCMag=maInstrFlag*nCCDs; // offest=0 if maInstrFlag=0
    comlsqr.offsetCnu=comlsqr.offsetCMag+nuInstrFlag*nFoVs*nCCDs; // offest=offsetCMag if nuInstrFlag=0
    comlsqr.offsetCdelta_eta=comlsqr.offsetCnu+ssInstrFlag*nCCDs*nPixelColumns; // offest=offsetCnu if ssInstrFlag=0
    comlsqr.offsetCDelta_eta_1=comlsqr.offsetCdelta_eta+lsInstrFlag*nFoVs*nCCDs*nTimeIntervals;
    comlsqr.offsetCDelta_eta_2=comlsqr.offsetCdelta_eta+lsInstrFlag*2*nFoVs*nCCDs*nTimeIntervals;
    comlsqr.offsetCDelta_eta_3=comlsqr.offsetCdelta_eta+lsInstrFlag*3*nFoVs*nCCDs*nTimeIntervals;
    comlsqr.offsetCdelta_zeta=comlsqr.offsetCDelta_eta_3+ssInstrFlag*nCCDs*nPixelColumns;
    comlsqr.offsetCDelta_zeta_1=comlsqr.offsetCdelta_zeta+lsInstrFlag*nFoVs*nCCDs*nTimeIntervals;
    comlsqr.offsetCDelta_zeta_2=comlsqr.offsetCdelta_zeta+lsInstrFlag*2*nFoVs*nCCDs*nTimeIntervals;
    comlsqr.nInstrPSolved=nInstrPSolved;

    
    if(myid==0 && lsInstrFlag && nElemIC!=0 ){
        double * instrCoeffConstr;
        instrCoeffConstr=(double *) calloc(nElemIC, sizeof(double));
        int * instrColsConstr;
        instrColsConstr=(int *) calloc(nElemIC, sizeof(int));
        
        if(!computeInstrConstr(comlsqr, instrCoeffConstr, instrColsConstr, instrConstrIlung))
        {
            printf("SEVERE ERROR PE=0 computeInstrConstr failed\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        //////////////////////////
        for(int k=0;k<nElemIC;k++)
            instrCoeffConstr[k]=instrCoeffConstr[k]*wgInstrCoeff;
        /////////////////////////
        for(int j=0;j<nElemIC;j++){
            systemMatrix[mapNcoeff[myid]+nOfElextObs*nEqExtConstr+nOfElBarObs*nEqBarConstr+j]=instrCoeffConstr[j];
            instrCol[mapNoss[myid]*nInstrPSolved+j]=instrColsConstr[j];
        }
        int counter0=0;
        for(int j=0;j<nOfInstrConstr;j++){
            double sumVal=0.;
            for(int k=0;k<instrConstrIlung[j];k++){
               sumVal+= systemMatrix[mapNcoeff[myid]+nOfElextObs*nEqExtConstr+nOfElBarObs*nEqBarConstr+counter0];
                counter0++;
            }
            if(idtest){
            knownTerms[mapNoss[myid]+nEqExtConstr+nEqBarConstr+j]=sumVal;
            } else{
                knownTerms[mapNoss[myid]+nEqExtConstr+nEqBarConstr+j]=0.;
            }
        }
        if(counter0 != nElemIC){
            printf("SEVERE ERROR PE=0 counter0=%d != nElemIC=%d when computing knownTerms for InstrConstr\n",counter0,nElemIC);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }

        //..... write files
        if(wrFilebin){
            chdir(wpath);
            chdir(wrfileDir);
            FILE *fpICCoeff,*fpICCols,*fpICIlung;
            fpICCoeff=fopen("instrConstrRows_Coeff.bin","wb");
            fpICCols=fopen("instrConstrRows_Cols.bin","wb");
            fpICIlung=fopen("instrConstrRows_Ilung.bin","wb");
        
            fwrite(instrConstrIlung,sizeof(int),nOfInstrConstr,fpICIlung);
            fwrite(instrCoeffConstr,sizeof(double),nElemIC,fpICCoeff);
            fwrite(instrColsConstr,sizeof(int),nElemIC,fpICCols);
        
            fclose(fpICCoeff);
            fclose(fpICCols);
            fclose(fpICIlung);
            chdir(wpath);
        }
        free(instrCoeffConstr);
        free(instrColsConstr);
     } // if(myid==0)
    MPI_Bcast(&systemMatrix[mapNcoeff[myid]+nOfElextObs*nEqExtConstr+nOfElBarObs*nEqBarConstr], nElemIC, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    MPI_Bcast(&instrCol[mapNoss[myid]*nInstrPSolved], nElemIC, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(instrConstrIlung, nOfInstrConstr, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&knownTerms[mapNoss[myid]+nEqExtConstr+nEqBarConstr], nOfInstrConstr, MPI_DOUBLE, 0, MPI_COMM_WORLD);

    /////// Search for map
    if(nAstroPSolved)
    {
        firstStar=matrixIndex[0]/nAstroPSolved;
        lastStar=matrixIndex[mapNoss[myid]*2-2]/nAstroPSolved;
        seqStar=lastStar-firstStar+1;
    }
    else{
        firstStar=0;
        lastStar=0;
    }
    if(extConstraint && (firstStar!=firstStarConstr || lastStar!=lastStarConstr+starOverlap)){
        printf("PE=%d Error extConstraint:  firstStar=%ld firstStarConstr=%ld lastStar=%ld lastStarConstr=%ld\n",myid,firstStar,firstStarConstr,lastStar,lastStarConstr);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
        
    }
    if(barConstraint && (firstStar!=firstStarConstr || lastStar!=lastStarConstr+starOverlap)){
        printf("PE=%d Error barConstraint:  firstStar=%ld firstStarConstr=%ld lastStar=%ld lastStarConstr=%ld\n",myid,firstStar,firstStarConstr,lastStar,lastStarConstr);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
        
    }
   
    ///////////////////////////////////////////////////////
    
    ////////////////////////////////////////////////////////////////////
    chdir(wpath);
    if(outputDirOpt) chdir(outputDir); // go to output dir
   
    //////////// Identity solution test
    // There are two possible IDtest modes: compatible and incompatible.
    // * in "compatible" mode the known term read from the input file is substituted by
    //   the sum of the coefficients of the system row, solving in this way a compatible
    //   system whose solution is exactly a vector of ones apart from numerical approximations
    //   (a further test is done in this case, in order to check that the difference between
    //   the sum and the known term is less than 1e-12)
    // * in "incompatible" mode the known term read from the input file is substituted by
    //   the sum of the coefficients of the system row further perturbed by a quantity x
    //   extracted from a random distribution having 0 mean and a sigma determined by the
    //   variable srIDtest. In this way the problem is reconducted to the solution of an
    //   incompatible system whose solution is close to the previous one. In particular,
    //   srIDtest represents the ratio between the sigma and the unperturbed known term
    //   (i.e. the sum of the coefficients since we are in IDtest mode). Therefore the
    //   perturbed known term is computed as
    //                KT_perturbed = KT_unperturbed*(1+x)
    
    
    if(idtest)   //if Identity test, overwrite the value of the knownterm
    {
        for(ii=0;ii<mapNoss[myid];ii++)
        {
            knownTerms[ii]=0.;
            for(jj=0;jj<nparam;jj++)
                knownTerms[ii] += systemMatrix[ii*nparam+jj];
        }
        if(srIDtest!=0.)
        {
            for(ii=0;ii<mapNoss[myid];ii++)
            {
                pert = gauss(0.0, srIDtest, idum);
                knownTerms[ii] *= (1.0+pert);
            }
        }
        
    }
    //////////////////////////////////////
    endTime=MPI_Wtime();
    timeToReadFiles=MPI_Wtime()-timeToReadFiles;
    printf("PE=%d time seconds=%lf to set system coefficients\n",myid, timeToReadFiles);
    /////  check, Fix map and dim
    if(seqStar<=1 && nAstroPSolved>0)
    {
        printf("ERROR PE=%d Only %d star Run not allowed with this PE numbers .n",myid,seqStar);
        exit(EXIT_FAILURE);
    }
    comlsqr.VrIdAstroPDim=seqStar;
    long tempDimMax=comlsqr.VrIdAstroPDim;
    long tempVrIdAstroPDimMax;
    MPI_Allreduce(&tempDimMax,&tempVrIdAstroPDimMax,1,MPI_LONG,MPI_MAX,MPI_COMM_WORLD);
   
    comlsqr.VrIdAstroPDimMax=tempVrIdAstroPDimMax;
    
    int **tempStarSend, **tempStarRecv;
    tempStarSend=(int **) calloc(nproc,sizeof(int *));
    for(int i=0;i<nproc;i++)
        tempStarSend[i]=(int *) calloc(2,sizeof(int));
    tempStarRecv=(int **) calloc(nproc,sizeof(int *));
    for(int i=0;i<nproc;i++)
        tempStarRecv[i]=(int *) calloc(2,sizeof(int));
    
    int *testVectSend, *testVectRecv;
    testVectSend=(int *) calloc(2*nproc,sizeof(int));
    testVectRecv=(int *) calloc(2*nproc,sizeof(int));
    testVectSend[2*myid]=firstStar;
    testVectSend[2*myid+1]=lastStar;
    
    MPI_Allreduce(testVectSend,testVectRecv,2*nproc,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
  
    
    comlsqr.mapStar=(int **) calloc(nproc,sizeof(int *));
    for(int i=0;i<nproc;i++)
        comlsqr.mapStar[i]=(int *) calloc(2,sizeof(int));
    for(int i=0;i<nproc;i++)
    {
        comlsqr.mapStar[i][0]=testVectRecv[2*i];
        comlsqr.mapStar[i][1]=testVectRecv[2*i+1];
    }
    
    
    for(int i=0;i<nproc;i++) {
        free(tempStarSend[i]);
        free(tempStarRecv[i]);
    }
    
    free(tempStarSend);
    free(tempStarRecv);
    
    if(comlsqr.mapStar[myid][0]==comlsqr.mapStar[myid][1] && nAstroPSolved>0)
    {
        printf("PE=%d ERROR. Only one star in this PE: starting star=%d ending start=%d\n",myid,comlsqr.mapStar[myid][0],comlsqr.mapStar[myid][1]);
        MPI_Abort(MPI_COMM_WORLD,1);
        exit(EXIT_FAILURE);
        
    }

    if(myid==0)
        for(int i=0;i<nproc;i++)printf("mapStar[%d][0]=%d mapStar[%d][1]=%d\n",i,comlsqr.mapStar[i][0],i,comlsqr.mapStar[i][1]);
    
    
    
     ////////  Check Null Space Vector
    if(extConstraint){
        seconds[4] = time(NULL);
        
        nullSpaceCk=cknullSpace(systemMatrix,matrixIndex,attNS,comlsqr);
        if(myid==0){
            printf("NullSpace check\n");
            for(int j=0;j<nEqExtConstr;j++)
                printf("Eq. Constraint %d: Norm=%15.7f Min=%15.7f Max=%15.7f Avg=%15.7f Var=%15.7f\n",j,nullSpaceCk.vectNorm[j],nullSpaceCk.compMin[j],nullSpaceCk.compMax[j],nullSpaceCk.compAvg[j],nullSpaceCk.compVar[j]);
        }
        seconds[5] = time(NULL);
        tot_sec[4] = seconds[5] - seconds[4];
        if(myid==0) printf("Time to check nullspace: %ld\n", tot_sec[4]);
        free(attNS);
    }
    ////////  WRITE BIN FILES and FileConstr_GsrSolProps.dat
    //////////
    if(wrFilebin){
        if(myid==0) printf("Writing bin files...\n");
        writeBinFiles(systemMatrix,matrixIndex,instrCol,knownTerms,wrfileDir,wpath,comlsqr,debugMode);
        MPI_Barrier(MPI_COMM_WORLD);
        if(myid==0)	printf(" finished.");
        
        if(myid==0){
            FILE* fpFilePosConstr;
            MPI_Status statusMpi;
            chdir(wpath);
            chdir(wrfileDir);
            int fileNum=-1;
            int nFiles;
            nFiles = scandir(".", &namelistMI, selMI, alphasort);
            if(debugMode)
                for(ii=0;ii<nFiles ;ii++){
                    printf("%s\n",namelistMI[ii]->d_name);
                }
            if(nFiles <=0)
            {
                printf("error on scandir n=%d\n",nFiles);
                MPI_Abort(MPI_COMM_WORLD,1);
            }
            
            ////////////////// Writing FileConstr_GsrSolProps.dat
            char fileConstr[512]="";
            strcat(fileConstr, "FileConstr_");
            strcat(fileConstr, filenameSolProps);
            fpFilePosConstr=fopen(fileConstr,"w");
            for(int k=0;k<counterConstr;k++){
                fileNum=constraintFound[k][0];
                fprintf(fpFilePosConstr,"%d %s %d\n", fileNum, namelistMI[fileNum]->d_name,constraintFound[k][1]);
            }
            if(debugMode)
                for(int k=0;k<counterConstr;k++){
                    fileNum=constraintFound[k][0];
                    printf("PE=%d %d %s %d\n",myid, fileNum, namelistMI[fileNum]->d_name,constraintFound[k][1]);
                }
            for(int np=1;np<nproc;np++){
                MPI_Recv(&counterConstr, 1, MPI_INT, np, 0, MPI_COMM_WORLD,&statusMpi);
                MPI_Recv(&constraintFound[0][0], counterConstr*2, MPI_INT, np, 1, MPI_COMM_WORLD,&statusMpi);
                
                for(int k=0;k<counterConstr;k++){
                    fileNum=constraintFound[k][0];
                    fprintf(fpFilePosConstr,"%d %s %d\n", fileNum, namelistMI[fileNum]->d_name,constraintFound[k][1]);
                }
                if(debugMode)
                    for(int k=0;k<counterConstr;k++){
                        fileNum=constraintFound[k][0];
                        printf("PE=%d %d %s %d\n",np, fileNum, namelistMI[fileNum]->d_name,constraintFound[k][1]);
                    }
            }
            fclose(fpFilePosConstr);
        } else {
            MPI_Send(&counterConstr, 1, MPI_INT, 0, 0, MPI_COMM_WORLD);
            MPI_Send(&constraintFound[0][0], counterConstr*2, MPI_INT, 0, 1, MPI_COMM_WORLD);
        }
    }
    
    VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax;
    VrIdAstroPDim=comlsqr.VrIdAstroPDim;
    VroffsetAttParam = VrIdAstroPDimMax*nAstroPSolved;
    comlsqr.VroffsetAttParam=VroffsetAttParam;
    
    nunkSplit =   VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam;
    comlsqr.nunkSplit=nunkSplit;
    nElements = nunkSplit ;
    preCondVect = (double *) calloc(nElements,sizeof(double));
    if (!preCondVect)
        exit(err_malloc("preCondVect",myid));
    totmem += nElements*sizeof(double);
    
    nElements = nunkSplit; //  the arrays v, w, x, and se require the same amount of memory (nunk * sizeof(double))
    vVect = (double *) calloc(nElements, sizeof(double));
    if (!vVect)
        exit(err_malloc("vVect",myid));
    totmem += nElements* sizeof(double);
    
    wVect = (double *) calloc(nElements, sizeof(double));
    if (!wVect)
        exit(err_malloc("wVect",myid));
    totmem += nElements* sizeof(double);
    
    xSolution = (double *) calloc(nElements, sizeof(double));
    if (!xSolution)
        exit(err_malloc("xSolution",myid));
    totmem += nElements* sizeof(double);
    
    standardError = (double *) calloc(nElements, sizeof(double));
    if (!standardError)
        exit(err_malloc("standardError",myid));
    totmem += nElements* sizeof(double);
    
    totmem += nElements* sizeof(double);  //dcopy+vAuxVect locally allocated on lsqr.c
    
    totmem=totmem/(1024*1024);  // mem in MB
    
    // Compute and write the total memory allocated
    if(myid==0)
    {
        printf("LOCAL %ld MB of memory allocated on each task.\n", totmem);
        printf("TOTAL MB memory allocated= %ld\n", nproc*totmem );
    }
    
    
    ///////////////////////////////
    
    
    
    
    MPI_Barrier(MPI_COMM_WORLD);
    // Compute the preconditioning vector for the system columns
    
    
    long int aIndex=0;
    long instrLongIndex[6];
    long iIelem=0;
    
    
    
    for (ii = 0; ii < mapNoss[myid]*multMI; ii=ii+multMI)
    {
        long int numOfStarPos=0;
        if(nAstroPSolved>0) numOfStarPos=matrixIndex[ii]/nAstroPSolved; // number of star associated to matrixIndex[ii]
        long int numOfAttPos=matrixIndex[ii+(multMI-1)];
        long int VrIdAstroPValue=-1; //
        
        VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0];
        if(VrIdAstroPValue==-1)
        {
            printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        for(int ns=0;ns<nAstroPSolved;ns++)
        {
            ncolumn=VrIdAstroPValue*nAstroPSolved+ns;
            ////
            // ncolumn = numOfStarPos+ns; // cancellato
            if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii,
                       matrixIndex[ii]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
            if(preCondVect[ncolumn]==0.0)
                printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn);
            aIndex++;
        }
        //
        for(int naxis=0;naxis<nAttAxes;naxis++)
            for(int ns=0;ns<nAttParAxis;ns++)
            {
                ncolumn = numOfAttPos+(VroffsetAttParam-offsetAttParam)+ns+naxis*nDegFreedomAtt;
                if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                    printf("ERROR. PE=%d numOfAttPos=%ld nStar*nAstroPSolved=%ld ncolumn=%ld ns=%d naxis=%d matrixIndex[ii+%d]=%ld\n",myid,numOfAttPos,nStar*nAstroPSolved,ncolumn, ns, naxis,multMI-1,matrixIndex[ii+(multMI-1)]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
                if(preCondVect[ncolumn]==0.0)
                    printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]); //   if aggiunto
                aIndex++;
            }
        ///// End of Attitude preCondVect
        
        
        if(nInstrPSolved>0)
        {
            for(int ns=0;ns<nInstrPSolved;ns++)
            {
                ncolumn = offsetInstrParam+(VroffsetAttParam-offsetAttParam)+instrCol[(ii/multMI)*nInstrPSolved+ns];
                if (ncolumn >= nunkSplit || ncolumn < 0 )
                {
                    printf("ERROR. PE=%d ii=%ld ",myid, ii);
                    for(int ke=0;ke<nInstrPSolved;ke++)
                        printf("instrCol[%d]=%d ",ii/multMI+ke, instrCol[ii/multMI+ke]);
                    printf("ncolumn=%ld   ns=%d\n", ncolumn, ns);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
//                if(ncolumn==55718+8000){
//                    printf("OSS BEFORE aindex=%ld SM[aindex]=%f preCondVect[%ld]=%f\n",aIndex,systemMatrix[aIndex],ncolumn,preCondVect[ncolumn]);
//                }
                preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
//                if(ncolumn==55718+8000){
//                    printf("OSS AFTER aindex=%ld SM[aindex]=%f preCondVect[%ld]=%f\n",aIndex,systemMatrix[aIndex],ncolumn,preCondVect[ncolumn]);
//                }
                if(preCondVect[ncolumn]==0.0)
                    printf("Instrument: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]);//   if aggiunto
                aIndex++;
            }
        }
        ////// End of Instruments preCondVect
        for(int ns=0;ns<nGlobP;ns++)
        {
            ncolumn = offsetGlobParam+(VroffsetAttParam-offsetAttParam)+ns;
            if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii+2]=%ld\n",myid, ncolumn, ii,
                       matrixIndex[ii]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
            if(preCondVect[ncolumn]==0)
                printf("Global: preCondVect[%ld]=0.0\n", ncolumn); //   if aggiunto
            aIndex++;
        }
    }
    
    /////  precondvect for extConstr
    if(extConstraint){
        if(aIndex!=mapNcoeff[myid]){
            printf("PE=%d. Error on aIndex=%ld  different of mapNcoeff[%d]=%ld\n",myid,aIndex,myid,mapNcoeff[myid]);
            MPI_Abort(MPI_COMM_WORLD, 1);
            exit(EXIT_FAILURE);
        }
        for(int i=0;i<nEqExtConstr;i++){
            long int numOfStarPos=0;
            if(nAstroPSolved>0) numOfStarPos=firstStarConstr; // number of star associated to matrixIndex[ii]
            long int numOfAttPos=startingAttColExtConstr;
            long int VrIdAstroPValue=-1; //
            
            VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0];
            if(VrIdAstroPValue==-1)
            {
                printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            for(int ns=0;ns<nAstroPSolved*numOfExtStar;ns++)
            {
                ncolumn=ns;
                ////
                // ncolumn = numOfStarPos+ns; // cancellato
                if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                    printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii,
                           matrixIndex[ii]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
                if(preCondVect[ncolumn]==0.0)
                    printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn);
                aIndex++;
            }
            //
            for(int naxis=0;naxis<nAttAxes;naxis++){
                for(int j=0;j<numOfExtAttCol;j++){
                    ncolumn = VrIdAstroPDimMax*nAstroPSolved+startingAttColExtConstr+j+naxis*nDegFreedomAtt;
                    if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                        printf("ERROR. PE=%d  ncolumn=%ld  naxis=%d j=%d\n",myid,ncolumn, naxis,j);
                        MPI_Abort(MPI_COMM_WORLD,1);
                        exit(EXIT_FAILURE);
                    }
                preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
//                    printf("%d %d %d %d\n",ncolumn, aIndex,j,naxis);
                if(preCondVect[ncolumn]==0.0)
                    printf("Attitude: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]); //   if aggiunto
                aIndex++;
                }
            }
        }
    }
    ///// end precondvect for extConstr
    /////  precondvect for barConstr
    if(barConstraint){
        for(int i=0;i<nEqBarConstr;i++){
            long int numOfStarPos=0;
            if(nAstroPSolved>0) numOfStarPos=firstStarConstr; // number of star associated to matrixIndex[ii]
            long int VrIdAstroPValue=-1; //
            
            VrIdAstroPValue=numOfStarPos-comlsqr.mapStar[myid][0];
            if(VrIdAstroPValue==-1)
            {
                printf("PE=%d ERROR. Can't find gsrId for precondvect.\n",myid);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            for(int ns=0;ns<nAstroPSolved*numOfBarStar;ns++)
            {
                ncolumn=ns;
                ////
                // ncolumn = numOfStarPos+ns; // cancellato
                if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                    printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n",myid, ncolumn, ii,
                           matrixIndex[ii]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
                if(preCondVect[ncolumn]==0.0)
                printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn);
                aIndex++;
            }
            //
        }
    }
    ///// end precondvect for barConstr

    /////  precondvect for instrConstr
    if(nElemIC>0){
        for(int i=0;i<nElemIC;i++){
            ncolumn=offsetInstrParam+(VroffsetAttParam-offsetAttParam)+instrCol[mapNoss[myid]*nInstrPSolved+i];
            if (ncolumn >= nunkSplit || ncolumn < 0 ) {
                printf("ERROR on instrConstr. PE=%d  ncolumn=%ld   i=%d\n",myid,ncolumn, i);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex];
            if(preCondVect[ncolumn]==0.0 && debugMode)
                printf("Instrument: PE=%d preCondVect[%ld]=0.0 aIndex=%ld systemMatrix[aIndex]=%12.9lf\n",myid, ncolumn,aIndex,systemMatrix[aIndex]);
            aIndex++;
            
        }
    }
    ////////////////
    
    
    MPI_Barrier(MPI_COMM_WORLD);
    
    // ACCUMULATE d
    //
    double *dcopy;
    dcopy=(double *)calloc(nAttParam,sizeof(double));
    if(!dcopy) exit(err_malloc("dcopy",myid));
    mpi_allreduce(&preCondVect[ VrIdAstroPDimMax*nAstroPSolved],dcopy,(long int) nAttParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    for(i=0;i<nAttParam;i++)
    {
        preCondVect[VrIdAstroPDimMax*nAstroPSolved+i]=dcopy[i];
        if(preCondVect[VrIdAstroPDimMax*nAstroPSolved+i]==0.0){
            printf("PE=%d preCondVect[%ld]=0!!\n", myid, VrIdAstroPDimMax*nAstroPSolved+i);
        }
    }
    free(dcopy);
    dcopy=(double *)calloc(nInstrParam,sizeof(double));
    if(!dcopy) exit(err_malloc("dcopy",myid));
    mpi_allreduce(&preCondVect[ VrIdAstroPDimMax*nAstroPSolved+nAttParam],dcopy,(long int) nInstrParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    for(i=0;i<nInstrParam;i++)
    {
        preCondVect[VrIdAstroPDimMax*nAstroPSolved+nAttParam+i]=dcopy[i];
        if(preCondVect[ VrIdAstroPDimMax*nAstroPSolved+nAttParam+i]==0.0) printf("PE=%d preCondVect[%d]=0!!\n", myid, i);
    }
    free(dcopy);
    dcopy=(double *)calloc(nGlobalParam,sizeof(double));
    if(!dcopy) exit(err_malloc("dcopy",myid));
    MPI_Allreduce(&preCondVect[VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam],dcopy,(int) nGlobalParam,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    for(i=0;i<nGlobalParam;i++)
    {
        preCondVect[VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam+i]=dcopy[i];
        if(preCondVect[VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam+i]==0.0) printf("PE=%d preCondVect[%d]=0!!\n", myid, i);
    }
    free(dcopy);
    if(nAstroPSolved) SumCirc(preCondVect,comlsqr);
    
    
    ////////// TEST for NO ZERO Column on A matrix
    
    
    if(precond){
        
        if(myid==0) printf("Inverting preCondVect\n");
        for (ii = 0; ii < VrIdAstroPDim*nAstroPSolved; ii++) {
            if(preCondVect[ii]==0.0)
            {
                if(ii-VrIdAstroPDimMax*nAstroPSolved<nAttParam)
                    printf("ERROR Att ZERO column: PE=%d preCondVect[%ld]=0.0 AttParam=%ld \n",myid,ii, ii-VrIdAstroPDimMax*nAstroPSolved);
                if(ii-VrIdAstroPDimMax*nAstroPSolved>nAttParam && ii-VrIdAstroPDimMax*nAstroPSolved<nAttParam+nInstrParam)
                    printf("ERROR Instr ZERO column: PE=%d preCondVect[%ld]=0.0 InstrParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam));
                if(ii-VrIdAstroPDimMax*nAstroPSolved>nAttParam+nInstrParam)
                    printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam));
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ii] = 1.0 / sqrt(preCondVect[ii]);
        }
        for (ii = VrIdAstroPDimMax*nAstroPSolved; ii < VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) {
            if(preCondVect[ii]==0.0)
            {
                printf("ERROR non-Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0\n",myid,ii);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ii] = 1.0 / sqrt(preCondVect[ii]);
        }
    } else{
        if(myid==0) printf("Setting preCondVect to 1.0\n");
        for (ii = 0; ii < VrIdAstroPDim*nAstroPSolved; ii++) {
            if(preCondVect[ii]==0.0)
            {
                printf("ERROR Astrometric ZERO column: PE=%d preCondVect[%ld]=0.0 Star=%ld\n",myid,ii,ii/nAstroPSolved);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ii] = 1.0;
        }
        for (ii = VrIdAstroPDimMax*nAstroPSolved; ii < VrIdAstroPDimMax*nAstroPSolved + nAttParam + nInstrParam + nGlobalParam; ii++) {
            if(preCondVect[ii]==0.0)
            {
                if(ii-VrIdAstroPDimMax*nAstroPSolved<nAttParam)
                    printf("ERROR Att ZERO column: PE=%d preCondVect[%ld]=0.0 AttParam=%ld \n",myid,ii, ii-VrIdAstroPDimMax*nAstroPSolved);
                if(ii-VrIdAstroPDimMax*nAstroPSolved>nAttParam && ii-VrIdAstroPDimMax*nAstroPSolved<nAttParam+nInstrParam)
                    printf("ERROR Instr ZERO column: PE=%d preCondVect[%ld]=0.0 InstrParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam));
                if(ii-VrIdAstroPDimMax*nAstroPSolved>nAttParam+nInstrParam)
                    printf("ERROR Global ZERO column: PE=%d preCondVect[%ld]=0.0 GlobalParam=%ld \n",myid,ii, ii-(VrIdAstroPDimMax*nAstroPSolved+nAttParam+nInstrParam));
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            preCondVect[ii] = 1.0;
        }
        
    }
    
    if(myid==0)
    {
        printf("Computation of the preconditioning vector finished.\n");
    }
    comlsqr.nvinc=0;
    comlsqr.parOss=(long) nparam;
    comlsqr.nunk=nunk;
    comlsqr.offsetAttParam=offsetAttParam;
    comlsqr.offsetInstrParam=offsetInstrParam;
    comlsqr.offsetGlobParam=offsetGlobParam;
    comlsqr.nAttParam=nAttParam;
    comlsqr.instrConst[0]=instrConst[0];
    comlsqr.instrConst[1]=instrConst[1];
    comlsqr.instrConst[2]=instrConst[2];
    comlsqr.instrConst[3]=instrConst[3];
    comlsqr.timeCPR=timeCPR;
    comlsqr.timeLimit=timeLimit;
    comlsqr.itnCPR=itnCPR;
    comlsqr.itnCPRstop=itnCPRstop;
    comlsqr.itnLimit=itnlim;
    comlsqr.Test=0; //it is not a running test but a production run
    
    initThread(matrixIndex,&comlsqr);
    
    
    endTime=MPI_Wtime()-endTime;
    MPI_Barrier(MPI_COMM_WORLD);
    
    ////////  WRITE BIN FILES and FileConstr_GsrSolProps.dat
    
    //////////////////////////// noiter
    if(noIter){
        if(myid==0) printf("\nEnd run: -noiter option givens.\n");
        
        MPI_Finalize();
        exit(EXIT_SUCCESS);
        
    }

    /////////////// MAIN CALL
    
    /////////////
    // This function computes the product of system matrix by precondVect. This avoids to compute the produsct in aprod for each iteration.
    if(myid==0 && debugMode)
    {
        printf("TEST LSQR START nobs=%ld, nunk=%ld, damp=%f, knownTerms[0]=%f, knownTerms[%ld]=%f, atol=%f btol=%f conlim=%f itnlim=%ld systemMatrix[0]=%f, systemMatrix[%ld]=%f   matrixIndex[0]=%ld matrixIndex[%ld]=%ld  instrCol[0]=%d instrCol[%ld]=%d preCondVect[0]=%f preCondVect[%ld]=%f preCondVect[%ld]=%f  nunkSplit=%ld\n",nobs,  nunk, damp, knownTerms[0], mapNoss[myid]-1,knownTerms[mapNoss[myid]-1], atol, btol,conlim,itnlim, systemMatrix[0],mapNcoeff[myid]-1,systemMatrix[mapNcoeff[myid]-1],  matrixIndex[0],mapNoss[myid]*2-1, matrixIndex[mapNoss[myid]*2-1], instrCol[0], mapNoss[myid]*nInstrPSolved -1,instrCol[mapNoss[myid]*nInstrPSolved -1], preCondVect[0],VrIdAstroPDim*nAstroPSolved-1, preCondVect[VrIdAstroPDim*nAstroPSolved-1], VrIdAstroPDim*nAstroPSolved, preCondVect[VrIdAstroPDim*nAstroPSolved], nunkSplit);
    }

    precondSystemMatrix(systemMatrix, preCondVect, matrixIndex,instrCol,comlsqr);

    //systemMatrix is passed to lsqr already conditioned.
    if(myid==0) {
        printf("System built, ready to call LSQR.\nPlease wait. System solution is underway...");
    }
    ////////////////////
    startTime=MPI_Wtime();
    seconds[1] = time(NULL);
    tot_sec[0] = seconds[1] - seconds[0];
    comlsqr.totSec=tot_sec[0];
    if(myid==0) printf("Starting lsqr after sec: %ld\n", tot_sec[0]);
    chdir(wpath);
    if(outputDirOpt) chdir(outputDir); // go to output dir
    lsqr(nobs, nunk, damp, knownTerms, vVect, wVect, xSolution, standardError, atol, btol, conlim,
         (int) itnlim, &istop, &itn, &anorm, &acond, &rnorm, &arnorm,
         &xnorm, systemMatrix, matrixIndex, instrCol, instrConstrIlung,preCondVect,comlsqr);
    
    ///////////////////////////
    MPI_Barrier(MPI_COMM_WORLD);
    endTime=MPI_Wtime();
    double clockTime=MPI_Wtime();
    if(myid==0) printf("Global Time for lsqr =%f sec \n",endTime-startTime);
    seconds[2] = time(NULL);
    tot_sec[1] = seconds[2] - seconds[1];
    
    long thetaCol=0, muthetaCol=0;
    ldiv_t res_ldiv;
    long flagTheta=0, flagMuTheta=0;
    
    if(nAstroPSolved==2) {
        thetaCol=1;
        muthetaCol=0;
    }
    else if(nAstroPSolved==3) {
        thetaCol=2;
        muthetaCol=0;
    }
    else if(nAstroPSolved==4) {
        thetaCol=1;
        muthetaCol=3;
    }
    else if(nAstroPSolved==5) {
        thetaCol=2;
        muthetaCol=4;
    }
    
    double epsilon;
    double localSumX=0, localSumXsq=0, sumX=0, sumXsq=0, average=0, rms=0;
    double dev=0, localMaxDev=0, maxDev=0;
    
    ///////// Each PE runs over the its Astrometric piece
    if(muthetaCol==0)
        for(ii=0; ii<VrIdAstroPDim*nAstroPSolved; ii++) {
            res_ldiv=ldiv((ii-thetaCol),nAstroPSolved);
            flagTheta=res_ldiv.rem;
            if(flagTheta==0)
            {
                xSolution[ii]*=(-preCondVect[ii]);
                if(idtest)
                {
                    epsilon=xSolution[ii]+1.0;
                    localSumX-=epsilon;
                    dev=fabs(epsilon);
                    if(dev>localMaxDev) localMaxDev=dev;
                }
            }
            else
            {
                xSolution[ii]*=preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081)
                if(idtest)
                {
                    epsilon=xSolution[ii]-1.0;
                    localSumX+=epsilon;
                    dev=fabs(epsilon);
                    if(dev>localMaxDev) localMaxDev=dev;
                }
                
            }
            if(idtest) localSumXsq+=epsilon*epsilon;
            standardError[ii]*=preCondVect[ii];
        }
    else
        for(ii=0; ii<VrIdAstroPDim*nAstroPSolved; ii++) {
            res_ldiv=ldiv((ii-thetaCol),nAstroPSolved);
            flagTheta=res_ldiv.rem;
            res_ldiv=ldiv((ii-muthetaCol),nAstroPSolved);
            flagMuTheta=res_ldiv.rem;
            if((flagTheta==0) || (flagMuTheta==0)) {
                xSolution[ii]*=(-preCondVect[ii]);
                if(idtest) {
                    epsilon=xSolution[ii]+1.0;
                    localSumX-=epsilon;
                    dev=fabs(epsilon);
                    if(dev>localMaxDev) localMaxDev=dev;
                }
            } else {
                xSolution[ii]*=preCondVect[ii]; // the corrections in theta are converted for consistency with the naming conventions in the Data Model to corrections in delta by a change of sign (Mantis Issue 0013081)
                if(idtest) {
                    epsilon=xSolution[ii]-1.0;
                    localSumX+=epsilon;
                    dev=fabs(epsilon);
                    if(dev>localMaxDev) localMaxDev=dev;
                }
            }
            if(idtest) localSumXsq+=epsilon*epsilon;
            standardError[ii]*=preCondVect[ii];
        }
    //////////// End of de-preconditioning for the Astrometric unknowns
    
    //////////// Then only PE=0 runs over the shared unknowns (Attitude, Instrument, and Global)
    if(myid==0)
        for(ii=VroffsetAttParam; ii<nunkSplit; ii++) {
            xSolution[ii]*=preCondVect[ii];
            if(idtest) {
                localSumX+=(xSolution[ii]-1.0);
                dev=fabs(1.0-xSolution[ii]);
                if(dev>localMaxDev) localMaxDev=dev;
                localSumXsq+=(xSolution[ii]-1.0)*(xSolution[ii]-1.0);
            }
            standardError[ii]*=preCondVect[ii];
        }
    //////////// End of de-preconditioning for the shared unknowns
    
    
    //////////////// TEST per verificare le somme di idtest.... da commentare/eliminare
    /*
     if(idtest)
     {
     if(myid==0)
     printf("TEST localSum. PE=%d, localSumX=%le, localSumXsq=%le, VrIdAstroPDim*nAstroPSolved+(nunkSplit-VroffsetAttParam)=%ld, localMaxDev=%le\n", myid, localSumX, localSumXsq, VrIdAstroPDim*nAstroPSolved+(nunkSplit-VroffsetAttParam), localMaxDev);
     else
     printf("TEST localSum. PE=%d, localSumX=%le, localSumXsq=%le, VrIdAstroPDim*nAstroPSolved=%ld, localMaxDev=%le\n", myid, localSumX, localSumXsq, VrIdAstroPDim*nAstroPSolved, localMaxDev);
     }
     */
    //////////////// Fine TEST
    
    free(systemMatrix);
    free(matrixIndex);
    free(instrCol);
    free(knownTerms);
    free(vVect);
    free(wVect);
    free(preCondVect);

    if(idtest)
    {
        MPI_Reduce(&localSumX,&sumX,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
        MPI_Reduce(&localSumXsq,&sumXsq,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
        MPI_Reduce(&localMaxDev,&maxDev,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
    }
    
    if(myid==0)
    {
        if(idtest)
        {
            average=sumX/nunk;
            rms=pow(sumXsq/nunk-pow(average,2.0),0.5);
            printf("Average deviation from ID solution: %le.\n", average);
            printf("                               rms: %le.\n", rms);
            printf("Maximum deviation from ID solution: %le.\n", maxDev);
        }
    }
    
    if(istop==1000 && wrsol==0)
    {
        if(myid==0){
            printf("Reached  limit at itn=%d. Execution stopped. CPR written!\n",itn);
            ffcont=fopen("CPR_CONT","w");
            fclose(ffcont);
        }
        MPI_Finalize();
        exit(EXIT_SUCCESS);
    }
    ///////////////
    ///////////// Initialize the output filename
    if(istop==1000){
        sprintf(filenameAstroResults, "%s%d%s", "GsrAstroParamSolution_intermediate_",itn,".bin"); // file storing the Astrometric Parameters
        sprintf(filenameAttResults, "%s%d%s","GsrAttitudeParamSolution_intermediate_",itn,".bin"); // file storing the Attitude Parameters
        sprintf(filenameInstrResults, "%s%d%s","GsrInstrParamSolution_intermediate_",itn,".bin"); // file storing the Instrument Parameters
        sprintf(filenameGlobalResults, "%s%d%s","GsrGlobalParamSolution_intermediate_",itn,".bin"); // file storing the Global Parameters
        sprintf(filenameSolPropsFinal,"GsrFinal_intermediate_%d_%s",itn, argv[argc-1]);
    } else{
        if(myid==0){
            printf("Execution finished END_FILE written!\n");
            chdir(wpath);
            ffcont=fopen("END_FILE","w");
            fclose(ffcont);
            system("rm CPR_CONT");
            if(outputDirOpt) chdir(outputDir); // go to output dir
        }
        sprintf(filenameAstroResults, "%s", "GsrAstroParamSolution.bin"); // file storing the Astrometric Parameters
        sprintf(filenameAttResults, "%s","GsrAttitudeParamSolution.bin"); // file storing the Attitude Parameters
        sprintf(filenameInstrResults, "%s","GsrInstrParamSolution.bin"); // file storing the Instrument Parameters
        sprintf(filenameGlobalResults, "%s","GsrGlobalParamSolution.bin"); // file storing the Global Parameters
        sprintf(filenameSolPropsFinal,"GsrFinal_%s", argv[argc-1]);
    }
    if(debugMode && myid==0) printf("Output %s %s %s %s %s\n",filenameSolProps,filenameAstroResults,	filenameAttResults, filenameInstrResults, filenameGlobalResults);
    /////////////////////////////////////////////
    
    MPI_Status statusMpi;
    if (myid == 0)
    {
        printf("Processing finished.\n");
        fflush(stdout);
        for (ii = 0; ii < 10; ii++)  printf("xSolution[%ld]=%16.12le \t standardError[%ld]=%16.12le\n", ii, xSolution[ii],ii,standardError[ii]);
    }
    
    ////////// Writing final solution
    
    if(myid==0)
    {
        
        ////////// Writing the raw results to the files GsrFinal_GsrSolProps.dat
        
        
        fpSolPropsFinal=fopen(filenameSolPropsFinal,"w");
        if (!fpSolPropsFinal)
        {
            printf("Error while open %s\n",filenameSolPropsFinal);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        fprintf(fpSolPropsFinal, "sphereId= %ld\n",sphereId) ;
        fprintf(fpSolPropsFinal, "nStar= %ld\n",nStar) ;
        fprintf(fpSolPropsFinal, "nAttParam= %d\n",nAttParam) ;
        fprintf(fpSolPropsFinal, "nInstrParam= %d\n",nInstrParam) ;
        fprintf(fpSolPropsFinal, "nInstrParamTot= %d\n",nInstrParamTot) ;
        fprintf(fpSolPropsFinal, "nGlobalParam= %d\n",nGlobalParam) ;
        fprintf(fpSolPropsFinal, "nAstroPSolved= %d\n",nAstroPSolved) ;
        fprintf(fpSolPropsFinal, "nInstrPSolved= %d\n",nInstrPSolved) ;
        fprintf(fpSolPropsFinal, "nFoVs= %d\n",nFoVs);
        fprintf(fpSolPropsFinal, "nCCDs= %d\n",nCCDs);
        fprintf(fpSolPropsFinal, "nPixelColumns= %d\n",nPixelColumns);
        fprintf(fpSolPropsFinal, "nTimeIntervals= %d\n",nTimeIntervals);
        fprintf(fpSolPropsFinal, "lSQRacond= %lf\n",acond) ;
        fprintf(fpSolPropsFinal, "lSQRanorm= %lf\n",anorm) ;
        fprintf(fpSolPropsFinal, "lSQRarnorm= %lf\n",arnorm) ;
        fprintf(fpSolPropsFinal, "lSQRitNumb= %d\n",itn) ;
        fprintf(fpSolPropsFinal, "lSQRrnorm= %lf\n",rnorm) ;
        fprintf(fpSolPropsFinal, "lSQRstopCond= %d\n",istop) ;
        fprintf(fpSolPropsFinal, "lSQRxnorm= %lf\n",xnorm) ;
        fclose(fpSolPropsFinal);
        
        ////////////////////// Writing GsrAstroParamSolution.bin
        if(nAstroPSolved)
        {
            long testOnStars=0;
            fpAstroR=fopen(filenameAstroResults,"wb");
            if (!fpAstroR)
            {
                printf("Error while open %s\n",filenameAstroResults);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fwrite(&sphereId,sizeof(long),1,fpAstroR);
            xAstro = (double *) calloc(VrIdAstroPDimMax*nAstroP, sizeof(double));
            if (!xAstro)
            {
                printf("Error allocating xAstro\n");
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            for(ii=0; ii<VrIdAstroPDim; ii++)
            {
                for(jj=0; jj< nAstroPSolved;jj++)
                {
                    xAstro[ii*nAstroP+mapAstroP[jj]] = xSolution[ii*nAstroPSolved+jj];
                }
            }
            fwrite(xAstro,sizeof(double),VrIdAstroPDim*nAstroP,fpAstroR);
            testOnStars+=VrIdAstroPDim;
            for(kk=1;kk<nproc;kk++)
            {
                offset=0;
                MPI_Recv(&VrIdAstroPDimRecv, 1, MPI_LONG, kk, 10, MPI_COMM_WORLD, &statusMpi);
                mpi_recv(xSolution, VrIdAstroPDimRecv*nAstroPSolved, MPI_DOUBLE, kk, 11,MPI_COMM_WORLD, &statusMpi);
                if(comlsqr.mapStar[kk][0]==comlsqr.mapStar[kk-1][1]) offset=1;
                for(ii=0+offset; ii<VrIdAstroPDimRecv; ii++)
                {