/* This program solves the relativistic sphere with an
 iterative method (conjugate gradients). Dynamic case.
 It is the c translation of the original fortran program
 solvergaia.for
 M.G. Lattanzi, A. Vecchiato, B. Bucciarelli (23 May 1996)
 A. Vecchiato, R. Morbidelli for the Fortran 90 version with
 dynamical memory allocation (21 March 2005)
 
 A. Vecchiato, June 16, 2008
 
 Version 2.1 , Feb 8, 2012
 Version history:
 - version 2.1, Feb  8 2012 - Added debug mode
 - version 2.0, Jan 27 2012 - Added Checkpoint & Restart
 - version 1.2, Jan 24 2012 -
 - version 5.0 - May 2013 from Ugo Becciani and Alberto Vecchiato
 */



//#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <dirent.h>
#include <math.h>
#include "lsqr.h"
#include <time.h>
#include <mpi.h>
#include "util.h"
#ifdef OMP
#include <omp.h>
#endif

#define MAX_CONSTR 1000


long instr_hash(int FoV, int CCD, int PixelColumn, int TimeInterval);
long randlong(long max);
long randlong1(long min, long max);
int randint(int max);
int randint1(int min, int max);
int fill_extract(long *values, long *pos_min, long pos_max, long *number);



/* Start of main program */
int main(int argc, char **argv) {
    int debugMode,wrsol;
    int i;
    
    int  idtest=0, precond=1;
    double srIDtest, pert;
    
    int inputDirOpt=0, outputDirOpt=0, inputDirLen;
    char inputDir[1024]="", outputDir[1024]="",wrfileDir[1024]="";
    char wpath[1024];
    size_t sizePath=1020;
    
    char filenameSolProps[150];
    char filenameSolPropsFinal[150];
    char filenameAstroResults[150]; /* file storing the Astrometric Parameters */
    char filenameAttResults[150]; /* file storing the Attitude Parameters */
    char filenameInstrResults[150]; /* file storing the Instrument Parameters */
    char filenameGlobalResults[150]; /* file storing the Global Parameters */
    
    long ii, jj, kk;
    long sphereId; // Id of the sphere to be solved
    long idum; // variable to initialize the random number generator aggiunta la variabile per inizializzare ran2()
    
    // LSQR input parameters
    long itnlim; // maximum number of iteration allowed for the LSQR convergence
    double damp, atol, btol, conlim; // other LSQR input parameters
    double aTol; //read by command line, overrides atol if >=0
    
    // LSQR output parameters
    int istop; // LSQR stopping condition
    int itn; // LSQR iteration number
    double anorm, acond, rnorm, arnorm, xnorm; // other LSQR output parameters
    
    // Properties of the system to be solved
    // Astrometric parameters:
    long nStar;      // number of stars
    short nAstroP;   // number of astrometric parameters for each observation
    
    short nAstroPSolved,nInstrPSolved;	// number of astrometric and instrumental parameters taken as unknowns in the system of equations
    int *mapAstroP;	 // indeces of the solved astrometric parameters
    int lsInstrFlag,ssInstrFlag,nuInstrFlag,maInstrFlag;
    int nElemIC=0,nOfInstrConstr=0;
    int nOfInstrConstrLSAL=0, nElemICLSAL=0, nOfInstrConstrLSAC=0, nElemICLSAC=0 , nOfInstrConstrSS=0 ,nElemICSS=0;

    // Attitude parameters:
    
    
    long nDegFreedomAtt=0; 	        // number of degrees of freedom for each attitude axis
    long cCDLSAACZP=14;                 // CCD light sensitive area AC zero point
    short nAttAxes=0, nAttParAxis;  // number of attitude axes to be reconstructed, number of non-zero attitude coefficients per obs. & axis
    short nAttP=0;     // number of attitude parameters for each observation
    // Instrument parameters:
    long instrSetUp; // number coding the set up of the instrument in terms of: # of FoVs, # of CCDs, # of pixel columns, # of time intervals
    short nInstrP;   // number of instrument parameters for each observation
    // Global parameters:
    short nGlobP;    // number of global parameters
    
    long nobs; // number of observations
    long nConstrLong, nConstrLat, nConstrMuLong, nConstrMuLat; // number of constraints for the astrometric parameters
    long *constrLongId, *constrLatId, *constrMuLongId, *constrMuLatId; // arrays with the gsrIDs of the constrained sources
    double *constrLongW, *constrLatW, *constrMuLongW, *constrMuLatW; // arrays with the weights of the constraints
    double extConstrW; // weight of the null space constraint equations
    double barConstrW; // weight of the baricentric constraint equations
    const double attExtConstrFact=-0.5; // the coefficients of the attitude part of the null space constraint equations must be multiplied by a factor -1/2
    
    // Variables to be read from the GsrSystemRows
    
    // Internal variables
    long nAstroParam;
    int nAttParam, nInstrParam, nGlobalParam,nInstrParamTot; // total number of unknowns for the Astrometric, Attitude, Instrument, and Global parameters respectively
    int nparam; // total number of unknowns
    long nunk, nunkSplit, nConstr;
    long ncoeff, ielem, global_ielem, precnofrows,global_precnofrows, ncolumn, nrowsToRead;
    long offsetAttParam, offsetInstrParam, offsetGlobParam, VroffsetAttParam; // offests for the Attitude, Instrument, and Global columns of the coefficient matrix
    unsigned long int totmem, nElements; // total required memory and temporary variable to store the memory being allocated
    long VrIdAstroPDimMax=0; //
    long VrIdAstroPDim=0; //
    long VrIdAstroPDimRecv;
    int offset;
    long nObsxStar, nobsOri;
    int nfileProc=3;
    int addElementAtt=0;
    int addElementextStar=0;
    int addElementbarStar=0;
   int nOfElextObs=0;
    int nOfElBarObs=0;
    double *attNS;

    
    // Array arguments of the LSQR function
    double *knownTerms, *vVect, *wVect, *xSolution, *standardError; // arrays containing the known terms (knownterms), the bidiagonalization (wVect, wVect), the unknowns (xSolution), and the estimated variances (standardError)
    long int  *matrixIndex; // ia[i] contains the column number of the coefficient systemMatrix[i]
    int *instrConst;
   	int *instrCol;  // columns on each observation for instumental Parameters
    int * instrConstrIlung; // vector containing the length of each instrConstr eq.
    double *systemMatrix, *preCondVect; // array of the non-zero coefficients of the system and of the column normalization respectively
    int wgInstrCoeff=1;
    
    /////////////////////
    // Arrays containing the solution coming from the LSQR
    double *xAstro, *standardErrorAstro; // solution and standard errors for the Astrometric parameters
    double *xAtt, *standardErrorAtt; // solution and standard errors for the Attitude parameters
    double *xInstr, *standardErrorInstr; // solution and standard errors for the Instrument parameters
    double *xGlobal, *standardErrorGlobal; // solution and standard errors for the Global parameters
    
    
    // file pointers
    FILE *fpSolProps,*fpSolPropsFinal,*fpSolPropsFileBin;
    FILE *fpMI, *fpSM, *fpII, *fpKT, *fpCPR;
    FILE *fpAstroR, *fpAttR, *fpInstrR, *fpGlobR, *ffcont;
    
    time_t seconds[10], tot_sec[10];
    seconds[0] = time(NULL);
    double timeToReadFiles;
    
    //MPI definitions
    int  nproc, myid, namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];
    double startTime,endTime;
    
    // Distributed array maps
    long int * mapNcoeff, * mapNoss;
    long int   mapNcoeffBefore, mapNossBefore;
    long int   mapNcoeffAfter, mapNossAfter;
    // struct to communicate with lsqr
    int timeCPR, timeLimit,itnCPR,itnCPRstop=1000000;
    int itnLimit; //read by command line, overrides itnlim if >0
    struct comData comlsqr;
    // serach for map
    long lastStar=-1;
    long firstStar=-1;
    long lastStarConstr=-1;
    long firstStarConstr=-1;
    int seqStar=1;
    int withFile=1;
    int noConstr=0;
    int   zeroAtt=0, zeroInstr=0, zeroGlob=0, wrFilebin=0;
    int constraintFound[MAX_CONSTR][2]; // first: number of file of the constraint; second: row in file
    int rowInFile=0;
    int noIter=0;
    int noCPR=0;
    int extConstraint=0, nEqExtConstr=0;
    int barConstraint=0, nEqBarConstr=0;
    int startingAttColExtConstr,endingAttColExtConstr,numOfExtAttCol,starOverlap=0,numOfExtStar,numOfBarStar;
    char constranitFoundFileName[MAX_CONSTR][256];
    struct dirent  **namelistMI;
    struct nullSpace nullSpaceCk;
    double nullSpaceAttfact;
    int autoRun=0;
    float memGlobal=0;
    float memGB;
    ///////////////
    for (int i=1;i<argc-1;i++)
    {
        if(strcmp (argv[i],"-testXCode") == 0)
           sleep(600);
    }

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &nproc);
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Get_processor_name(processor_name,&namelen);
    startTime= MPI_Wtime();
    // Initialize the random number generator for each PE
    idum=-((long)startTime+myid);
    
    
    if (argc == 1) {
        if(myid==0){
            printf("Usage:  solvergaiaSim [-auto] [-memGlobal value]  [-IDtest [value]   -noFile -noConstr -numFilexproc nfileProc -Precond [on|off] -timeCPR hours -timelimit hours -itnCPR numberOfIterations -noCPR -itnlimit numberOfIterations   -atol value  -inputDir inputdir -outputDir outputdir      -zeroAtt -zeroInstr -zeroGlob -wgic value]] -wrfilebin writedir -wrsol -noiter  -extConstr weight -barConstr weight filename\n\n");
        }
        MPI_Finalize();
        exit(EXIT_FAILURE);
    }
    debugMode=0;
    wrsol=0;
    int testTime=0;
    short testSolved=-1;
    nAstroPSolved=-1;
    idtest=0;
    
    //idtest=1;  // IDTest forced for the simulator
    
    timeCPR=DEFAULT_TIMECPR;
    timeLimit=DEFAULT_TIMELIMIT;
    itnCPR=DEFAULT_ITNCPR;
    itnLimit=0;
    aTol=-1.0;
    if(strcmp (argv[1],"-auto") == 0){
        printf(" solvergaiaSim -auto ");
        idtest=1;
        srIDtest=1;
        withFile=0;
        noConstr=1;
        noCPR=1;
        itnLimit=20000;
        wgInstrCoeff=100;
        barConstraint=1;
        nEqBarConstr=DEFAULT_BARCONSTROWS;
        barConstrW=10.0;
        autoRun=1;
        for (int i=2;i<argc-1;i++){
            if(strcmp (argv[i],"-memGlobal") == 0)
            {
                memGlobal=atof(argv[i+1]);
                 printf("-memGlobal %f",memGlobal);
            }
 
        }
        printf("\n\n");
    }
    
     for (int i=1;i<argc-1;i++)
     {
        if(strcmp (argv[i],"-IDtest") == 0)
        {
            idtest=1;
            srIDtest=atof(argv[i+1]);
        }
        if(strcmp (argv[i],"-Precond") == 0)  // aggiunto precond on/off
            if(strcmp (argv[i+1],"off") == 0)  // aggiunto precond on/off
                precond=0;
        
        if(strcmp (argv[i],"-noFile") == 0) withFile=0;
        if(strcmp (argv[i],"-noConstr") == 0) noConstr=1;
        
        
        if(strcmp (argv[i],"-numFilexproc") == 0)
        {
            nfileProc=atoi(argv[i+1]);
            if(nfileProc <=0) nfileProc=3;
        }
        if(strcmp (argv[i],"-timeCPR") == 0)
        {
            testTime=atoi(argv[i+1]);
            if(testTime >0) timeCPR=testTime;
        }
        if(strcmp (argv[i],"-timelimit") == 0)
        {
            testTime=atoi(argv[i+1]);
            if(testTime >0) timeLimit=testTime;
        }
        
        if(strcmp (argv[i],"-itnCPR") == 0)
        {
            testTime=atoi(argv[i+1]);
            if(testTime >0) itnCPR=testTime;
        }
        if(strcmp (argv[i],"-noCPR") == 0)
        {
            noCPR=1;
        }
        if(strcmp (argv[i],"-itnlimit") == 0)
        {
            testTime=atoi(argv[i+1]);
            if(testTime >0) itnLimit=testTime;
        }
        if(strcmp (argv[i],"-atol") == 0)
        {
            double dummy=atof(argv[i+1]);
            if(dummy >=0) aTol=dummy;
        }
        if(strcmp (argv[i],"-zeroAtt") == 0)
        {
            zeroAtt=1;
        }
        if(strcmp (argv[i],"-zeroInstr") == 0)
        {
            zeroInstr=1;
        }
        if(strcmp (argv[i],"-zeroGlob") == 0)
        {
            zeroGlob=1;
        }
        if(strcmp (argv[i],"-wgic") == 0)
        {
            wgInstrCoeff=atoi(argv[i+1]);
            if(wgInstrCoeff<=0) wgInstrCoeff=1;
        }

        // inputDir e outputDir
        
        if(strcmp (argv[i],"-inputDir") == 0)
        {
            sprintf(inputDir, "%s", argv[i+1]);
            inputDirOpt=1;
        }
        
        if(strcmp (argv[i],"-outputDir") == 0)
        {
            sprintf(outputDir, "%s", argv[i+1]);
            outputDirOpt=1;
        }
        
        // Fine input e output Dir
        
        if(strcmp (argv[i],"-debug") == 0)
            debugMode=1;
        
        if(strcmp (argv[i],"-wrsol") == 0){
            wrsol=1;
        }
        
        if(strcmp (argv[i],"-noiter") == 0){
            noIter=1;
        }
        
        if(strcmp (argv[i],"-wrfilebin") == 0)
        {
            sprintf(wrfileDir, "%s", argv[i+1]);
            wrFilebin=1;
        }
        if(strcmp (argv[i],"-extConstr") == 0)
        {
            extConstraint=1; //extendet external constraint
            nEqExtConstr=DEFAULT_EXTCONSTROWS;
            extConstrW=atof(argv[i+1]);
            if(extConstrW==0.0){
                printf("ERROR: PE=%d -extConstr option given with no value or zero value ==> %le. Aborting\n",myid,extConstrW);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        }
        if(strcmp (argv[i],"-barConstr") == 0)
        {
            barConstraint=1; //extendet external constraint
            nEqBarConstr=DEFAULT_BARCONSTROWS;
            barConstrW=atof(argv[i+1]);
            if(barConstrW==0.0){
                printf("ERROR: PE=%d -barConstr option given with no value or zero value ==> %le. Aborting\n",myid,barConstrW);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        }
        
    }
    if (extConstraint) noConstr=2;
    if (barConstraint) noConstr=2;
   
    
    if(myid==0)
    {
        printf("Execution of solvergaia Simulator version %s on %d mpi-tasks\n solvergaiaSim ",VER,nproc); 
        if(idtest)printf ("-IDtest %le ", srIDtest);
        if(precond)printf ("-Precond on");
        else printf ("-Precond off");
        if(nfileProc!=3) printf ("-numFilexproc %d ",nfileProc);
       if(wrFilebin) printf(" -wrfilebin %s ",wrfileDir);
        if(!withFile) printf(" -noFile ");
        if(itnLimit>0) printf(" -itnlimit %d ",itnLimit);
        if(noCPR>0) printf(" -noCPR ");
        if(itnCPR!= DEFAULT_ITNCPR) printf(" -itnCPR %d ", itnCPR);
        if(zeroAtt) printf(" -zeroAtt ");
        if(noConstr==1) printf(" -noConstr ");
        if(zeroInstr) printf(" -zeroInstr ");
        if(zeroGlob) printf(" -zeroGlob ");
        if(inputDirOpt) printf(" -inputDir %s  ", inputDir);
        if(outputDirOpt) printf(" -outputDir %s ", outputDir);
        if(wrsol) printf(" -wrsol ");
        if(noIter) printf(" -noiter ");
        if(debugMode) printf(" -debug ");
        if(extConstraint)printf("-extConstr %le ",extConstrW);
        if(barConstraint)printf("-barConstr %le ",barConstrW);
		printf("-wgic %d", wgInstrCoeff);
        if(!autoRun) printf(" %s\n",argv[argc-1]);
        if(extConstraint && barConstraint) {
            printf("Error: baricentric anld null space constraints are mutually exclusive. Aborting\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin && strcmp(inputDir, wrfileDir)==0){
            printf("inputDir and wrfilebinDir are the same. Execution Aborted.\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
            
        }
    }
    /////////////////
    getcwd(wpath,sizePath);
    if(!inputDirOpt)
        strcpy(inputDir, wpath);
    printf("Process %d running on %s\n",myid,processor_name);
#ifdef OMP
#pragma omp parallel
    {
        int tid = omp_get_thread_num();
        int nthreads = omp_get_num_threads();
        printf("PE=%d  on processor %s total num of threads =%d   ID thread =%d \n",myid, processor_name,nthreads,tid);
       	
    }
#endif
    nullSpaceAttfact=1.0*attExtConstrFact*extConstrW;
    comlsqr.extConstrW=extConstrW;
    comlsqr.nullSpaceAttfact=nullSpaceAttfact;
    comlsqr.barConstrW=barConstrW;
   
    
    if(inputDirOpt)
    {
        if(myid==0) printf("Checking input directory %s ... ", inputDir);
        if(!(chdir(inputDir)==0))  {
            printf("Input directory does not exist. Aborting\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        chdir(wpath);
        if(myid==0) printf("success.\n");
    }
    inputDirLen=strlen(inputDir);
    sprintf(filenameSolProps, "%s", argv[argc-1]);
    sprintf(filenameSolPropsFinal,"GsrFinal_%s", argv[argc-1]);
    
    if(outputDirOpt)
    {
        if(myid==0) printf("Checking output directory %s ...", outputDir);
        if(!(chdir(outputDir)==0)) {
            printf("Output directory does not exist on PE %d. Aborting\n", myid);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(myid==0) printf("success.\n");
    }
  
    ///////////// Initialize the output filename
    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
    
    
//    instrConst=(int *) calloc(DEFAULT_NINSTRINDEXES+1 , sizeof(int));
    instrConst=(int *) calloc(DEFAULT_NINSTRINDEXES , sizeof(int));
    nAttParAxis=4;
    //START READING filenameSolProps
    //START READING filenameSolProps GsrSolProps.dat
    if(myid==0 && wrFilebin){
        chdir(wpath);
        chdir(wrfileDir);
        fpSolPropsFileBin=fopen(filenameSolProps,"w");
    }
    chdir(wpath);
    if(inputDirOpt) chdir(inputDir);
    if(myid==0)
    {
        if (autoRun){
            sphereId=0;
            atol= 0.000000;
            btol= 0.000000;
            conlim= 10000000000000.000000;
            itnlim= 2000;
            if(itnLimit>0)
                itnlim=itnLimit;
            damp= 0.000000;
            nStar= 10000;
            nAstroP= 5;
            nAstroPSolved= 5;
            nConstrLat= 0;
            nConstrLong= 0;
            nConstrMuLat= 0;
            nConstrMuLong= 0;
            nDegFreedomAtt= 230489;
            nAttAxes= 3;
            instrConst[0]= 1;
            instrConst[1]= 62;
            instrConst[2]= 1966;
            instrConst[3]= 22;
            nInstrP= 6;
            nInstrPSolved= 6;
            lsInstrFlag= 1;
            ssInstrFlag= 1;
            nuInstrFlag= 1;
            maInstrFlag= 1;
            nOfInstrConstr= -1;
            nElemIC= -1;
            nGlobP= 0;
            nobs= 2400000;

            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;
            if(nInstrP==0) nInstrParamTot=0;
            nGlobalParam = nGlobP;
            
            nparam = nAstroPSolved + nAttP + nInstrPSolved + nGlobP;
            if(memGlobal!=0){
                memGB=simfullram(nStar,nobs,memGlobal,nparam,nAttParam,nInstrParam);
                printf("Running with memory %f GB, nStar=%d nobs=%ld\n",memGB, nStar,nobs);
            }
            printf("sphereId= %7ld\n", sphereId);
            printf("atol= %18.15lf\n", atol);
            printf("btol= %18.15lf\n", btol);
            printf("conlim= %18.15le\n", conlim);
            printf("itnlim= %7ld\n", itnlim);
            printf("damp= %18.15lf\n", damp);
            printf("nStar= %7ld\n", nStar);
            printf("nAstroP= %7hd\n", nAstroP);
            printf("nAstroPSolved= %hd\n", nAstroPSolved);
            printf("nConstrLat= %5ld\n", nConstrLat);
            printf("nConstrLong= %5ld\n", nConstrLong);
            printf("nConstrMuLat= %5ld\n", nConstrMuLat);
            printf("nConstrMuLong= %5ld\n", nConstrMuLong);
            printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt);
            printf("nAttAxes= %7hd\n", nAttAxes);
            printf("nFovs= %7d ", instrConst[0]+1);
            printf("(instrConst[0])= %7d\n", instrConst[0]);
            printf("nCCDs= %7d\n", instrConst[1]);
            printf("nPixelColumns= %7d\n", instrConst[2]);
            printf("nTimeIntervals= %7d\n", instrConst[3]);
            printf("nInstrP= %7hd\n", nInstrP);
            printf("nInstrPSolved= %7hd\n", nInstrPSolved);
            printf("lsInstrFlag= %7d\n", lsInstrFlag);
            printf("ssInstrFlag= %7d\n", ssInstrFlag);
            printf("nuInstrFlag= %7d\n", nuInstrFlag);
            printf("maInstrFlag= %7d\n", maInstrFlag);
            printf("nOfInstrConstr= %7d\n", nOfInstrConstr);
            printf("nElemIC= %7d\n", nElemIC);
            printf("nGlobP= %7hd\n", nGlobP);
            printf("nObs= %10ld\n", nobs);
            if(wrFilebin){
                fprintf(fpSolPropsFileBin,"sphereId= %ld\n", sphereId);
                fprintf(fpSolPropsFileBin,"atol= %lf\n", atol);
                fprintf(fpSolPropsFileBin,"btol= %lf\n", btol);
                fprintf(fpSolPropsFileBin,"conlim= %lf\n", conlim);
                fprintf(fpSolPropsFileBin,"itnlim= %ld\n", itnlim);
                fprintf(fpSolPropsFileBin,"damp= %lf\n", damp);
                fprintf(fpSolPropsFileBin,"nStar= %ld\n", nStar);
                fprintf(fpSolPropsFileBin,"nAstroP= %hd\n", nAstroP);
                fprintf(fpSolPropsFileBin,"nAstroPSolved= %hd\n", nAstroPSolved);
                fprintf(fpSolPropsFileBin,"nConstrLat= 0\n");
                fprintf(fpSolPropsFileBin,"nConstrLong= 0\n");
                fprintf(fpSolPropsFileBin,"nConstrMuLat= 0\n");
                fprintf(fpSolPropsFileBin,"nConstrMuLong= 0\n");
                fprintf(fpSolPropsFileBin,"nDegFreedomAtt= %ld\n", nDegFreedomAtt);
                fprintf(fpSolPropsFileBin,"nAttAxes= %hd\n", nAttAxes);
                fprintf(fpSolPropsFileBin,"nFoVs= %d\n", instrConst[0]);
                fprintf(fpSolPropsFileBin,"nCCDs= %d\n", instrConst[1]);
                fprintf(fpSolPropsFileBin,"nPixelColumns= %d\n", instrConst[2]);
                fprintf(fpSolPropsFileBin,"nTimeIntervals= %d\n", instrConst[3]);
                fprintf(fpSolPropsFileBin,"nInstrP= %hd\n", nInstrP);
                fprintf(fpSolPropsFileBin,"nInstrPSolved= %hd\n", nInstrPSolved);
                fprintf(fpSolPropsFileBin,"lsInstrFlag= %hd\n", lsInstrFlag);
                fprintf(fpSolPropsFileBin,"ssInstrFlag= %hd\n", ssInstrFlag);
                fprintf(fpSolPropsFileBin,"nuInstrFlag= %hd\n", nuInstrFlag);
                fprintf(fpSolPropsFileBin,"maInstrFlag= %hd\n", maInstrFlag);
                fprintf(fpSolPropsFileBin,"nOfInstrConstr= %d\n", nOfInstrConstr);
                fprintf(fpSolPropsFileBin,"nElemIC= %d\n", nElemIC);
                fprintf(fpSolPropsFileBin,"nGlobP= %hd\n", nGlobP);
                fprintf(fpSolPropsFileBin,"nObs= %ld\n", nobs);

                fclose(fpSolPropsFileBin);
            }
        }else{
        fpSolProps=fopen(filenameSolProps,"r");
        if (!fpSolProps)
        {
            printf("Error while open %s\n",filenameSolProps);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        
        if(fscanf(fpSolProps, "sphereId= %ld\n",&sphereId) != 1) {
            printf("Error reading sphereId  %ld \n",sphereId);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("sphereId= %7ld\n", sphereId);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"sphereId= %ld\n", sphereId);
        
        if(fscanf(fpSolProps, "atol= %lf\n",&atol) != 1) {
            printf("Error reading atol  %lf \n",atol);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(aTol >=0)
            atol=aTol;
        printf("atol= %18.15lf\n", atol);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"atol= %lf\n", atol);
        
        if(fscanf(fpSolProps, "btol= %lf\n",&btol) != 1) {
            printf("Error reading btol  %lf \n",btol);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("btol= %18.15lf\n", btol);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"btol= %lf\n", btol);
        
        if(fscanf(fpSolProps, "conlim= %lf\n",&conlim) != 1) {
            printf("Error reading conlim  %lf \n",conlim);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("conlim= %18.15le\n", conlim);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"conlim= %lf\n", conlim);
        
        if(fscanf(fpSolProps, "itnlim= %ld\n",&itnlim) != 1) {
            printf("Error reading itnlim  %ld \n",itnlim);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(itnLimit>0)
            itnlim=itnLimit;
        printf("itnlim= %7ld\n", itnlim);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"itnlim= %ld\n", itnlim);
        
        if(fscanf(fpSolProps, "damp= %lf\n",&damp) != 1) {
            printf("Error reading damp  %lf \n",damp);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("damp= %18.15lf\n", damp);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"damp= %lf\n", damp);
        
        if(fscanf(fpSolProps, "nStar= %ld\n",&nStar) != 1) {
            printf("Error reading nStar  %ld \n",nStar);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nStar= %7ld\n", nStar);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nStar= %ld\n", nStar);
        
        if(fscanf(fpSolProps, "nAstroP= %hd\n",&nAstroP) != 1) {
            printf("Error reading nAstroP  %hd \n",nAstroP);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nAstroP= %7hd\n", nAstroP);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nAstroP= %hd\n", nAstroP);
        
        if(fscanf(fpSolProps, "nAstroPSolved= %hd\n",&nAstroPSolved) != 1) {
            printf("Error reading nAstroPSolved  %hd \n",nAstroPSolved);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nAstroPSolved= %hd\n", nAstroPSolved);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nAstroPSolved= %hd\n", nAstroPSolved);
        
        if(fscanf(fpSolProps, "nConstrLat= %ld\n",&nConstrLat) != 1) {
            printf("Error reading nConstrLat  %ld \n",nConstrLat);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nConstrLat= %5ld\n", nConstrLat);
        if(wrFilebin){
            if(noConstr) fprintf(fpSolPropsFileBin,"nConstrLat= 0\n");
            else fprintf(fpSolPropsFileBin,"nConstrLat= %ld\n", nConstrLat);
        }
        if(nConstrLat>0)
        {
            constrLatId = (long *) calloc(nConstrLat, sizeof(long));
            for(i=0;i<nConstrLat-1;i++) {
                if(fscanf(fpSolProps, "%ld ",&constrLatId[i]) != 1) {
                    printf("Error reading constrLatId[%d]  %ld \n",i,constrLatId[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrLatId[%d]=%7ld ",i,constrLatId[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld ",constrLatId[i]);
            }
            if(fscanf(fpSolProps, "%ld\n",&constrLatId[nConstrLat-1]) != 1) {
                printf("Error reading constrLatId[nConstrLat-1] %ld \n", constrLatId[nConstrLat-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrLatId[nConstrLat-1]=%7ld\n",constrLatId[nConstrLat-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld\n",constrLatId[nConstrLat-1]);
            
            constrLatW = (double *) calloc(nConstrLat, sizeof(double));
            for(i=0;i<nConstrLat-1;i++) {
                if(fscanf(fpSolProps, "%lf ",&constrLatW[i]) != 1) {
                    printf("Error reading constrLatW[%d]  %lf \n",i,constrLatW[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrLatW[%d]=%18.15le ",i,constrLatW[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrLatW[i]);
            }
            if(fscanf(fpSolProps, "%lf\n",&constrLatW[nConstrLat-1]) != 1) {
                printf("Error reading constrLatW[nConstrLat-1]  %lf \n", constrLatW[nConstrLat-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrLatW[nConstrLat-1]=%18.15lf\n",constrLatW[nConstrLat-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrLatW[nConstrLat-1]);
        }
        
        if(fscanf(fpSolProps, "nConstrLong= %ld\n",&nConstrLong) != 1) {
            printf("Error reading nConstrLong  %ld \n",nConstrLong);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nConstrLong= %5ld\n", nConstrLong);
        if(wrFilebin){
            if(noConstr) fprintf(fpSolPropsFileBin,"nConstrLong= 0\n");
            else fprintf(fpSolPropsFileBin,"nConstrLong= %ld\n", nConstrLong);
        }
        if(nConstrLong>0)
        {
            constrLongId = (long *) calloc(nConstrLong, sizeof(long));
            for(i=0;i<nConstrLong-1;i++) {
                if(fscanf(fpSolProps, "%ld ",&constrLongId[i]) != 1) {
                    printf("Error reading constrLongId[%d]  %ld \n",i,constrLongId[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrLongId[%d]=%7ld ",i,constrLongId[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld ",constrLongId[i]);
            }
            if(fscanf(fpSolProps, "%ld\n",&constrLongId[nConstrLong-1]) != 1) {
                printf("Error reading constrLongId[nConstrLong-1]  %ld \n", constrLongId[nConstrLong-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrLongId[nConstrLong-1]=%7ld\n",constrLongId[nConstrLong-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld\n",constrLongId[nConstrLong-1]);
            
            constrLongW = (double *) calloc(nConstrLong, sizeof(double));
            for(i=0;i<nConstrLong-1;i++) {
                if(fscanf(fpSolProps, "%lf ",&constrLongW[i]) != 1) {
                    printf("Error reading constrLongW[%d]  %lf \n",i,constrLongW[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrLongW[%d]=%18.15le ",i,constrLongW[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrLongW[i]);
            }
            if(fscanf(fpSolProps, "%lf\n",&constrLongW[nConstrLong-1]) != 1) {
                printf("Error reading constrLongW[nConstrLong-1]  %lf \n", constrLongW[nConstrLong-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrLongW[nConstrLong-1]=%18.15lf\n",constrLongW[nConstrLong-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrLongW[nConstrLong-1]);
        }
        
        if(fscanf(fpSolProps, "nConstrMuLat= %ld\n",&nConstrMuLat) != 1) {
            printf("Error reading nConstrMuLat  %ld \n",nConstrMuLat);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nConstrMuLat= %5ld\n", nConstrMuLat);
        if(wrFilebin){
            if(noConstr) fprintf(fpSolPropsFileBin,"nConstrMuLat= 0\n");
            else fprintf(fpSolPropsFileBin,"nConstrMuLat= %ld\n", nConstrMuLat);
        }
        
        if(nConstrMuLat>0)
        {
            constrMuLatId = (long *) calloc(nConstrMuLat, sizeof(long));
            for(i=0;i<nConstrMuLat-1;i++) {
                if(fscanf(fpSolProps, "%ld ",&constrMuLatId[i]) != 1) {
                    printf("Error reading constrMuLatId[%d]  %ld \n",i,constrMuLatId[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrMuLatId[%d]=%7ld ",i,constrMuLatId[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld ",constrMuLatId[i]);
            }
            if(fscanf(fpSolProps, "%ld\n",&constrMuLatId[nConstrMuLat-1]) != 1) {
                printf("Error reading constrMuLatId[nConstrMuLat-1]  %ld \n", constrMuLatId[nConstrMuLat-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrMuLatId[nConstrMuLat-1]=%7ld\n",constrMuLatId[nConstrMuLat-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld\n",constrMuLatId[nConstrMuLat-1]);
            
            constrMuLatW = (double *) calloc(nConstrMuLat, sizeof(double));
            for(i=0;i<nConstrMuLat-1;i++) {
                if(fscanf(fpSolProps, "%lf ",&constrMuLatW[i]) != 1) {
                    printf("Error reading constrMuLatW[%d]  %lf \n",i,constrMuLatW[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrMuLatW[%d]=%18.15le ",i,constrMuLatW[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrMuLatW[i]);
            }
            if(fscanf(fpSolProps, "%lf\n",&constrMuLatW[nConstrMuLat-1]) != 1) {
                printf("Error reading constrMuLatW[nConstrMuLat-1]  %lf \n", constrMuLatW[nConstrMuLat-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrMuLatW[nConstrMuLat-1]=%18.15lf\n",constrMuLatW[nConstrMuLat-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrMuLatW[nConstrMuLat-1]);
            
        }
        
        if(fscanf(fpSolProps, "nConstrMuLong= %ld\n",&nConstrMuLong) != 1) {
            printf("Error reading nConstrMuLong  %ld \n",nConstrMuLong);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nConstrMuLong= %5ld\n", nConstrMuLong);
        if(wrFilebin){
            if(noConstr) fprintf(fpSolPropsFileBin,"nConstrMuLong= 0\n");
            else fprintf(fpSolPropsFileBin,"nConstrMuLong= %ld\n", nConstrMuLong);
        }
        if(nConstrMuLong>0)
        {
            constrMuLongId = (long *) calloc(nConstrMuLong, sizeof(long));
            for(i=0;i<nConstrMuLong-1;i++) {
                if(fscanf(fpSolProps, "%ld ",&constrMuLongId[i]) != 1) {
                    printf("Error reading constrMuLongId[%d]  %ld \n",i,constrMuLongId[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrMuLongId[%d]=%7ld ",i,constrMuLongId[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%7ld ",constrMuLongId[i]);
            }
            if(fscanf(fpSolProps, "%ld\n",&constrMuLongId[nConstrMuLong-1]) != 1) {
                printf("Error reading constrMuLongId[nConstrMuLong-1]  %ld \n", constrMuLongId[nConstrMuLong-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrMuLongId[nConstrMuLong-1]=%7ld\n",constrMuLongId[nConstrMuLong-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%ld\n",constrMuLongId[nConstrMuLong-1]);
            
            constrMuLongW = (double *) calloc(nConstrMuLong, sizeof(double));
            for(i=0;i<nConstrMuLong-1;i++) {
                if(fscanf(fpSolProps, "%lf ",&constrMuLongW[i]) != 1) {
                    printf("Error reading constrMuLongW[%d]  %lf \n",i,constrMuLongW[i]);
                    MPI_Abort(MPI_COMM_WORLD,1);
                    exit(EXIT_FAILURE);
                }
                if(debugMode) printf("constrMuLongW[%d]=%18.15le ",i,constrMuLongW[i]);
                if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf ",constrMuLongW[i]);
            }
            if(fscanf(fpSolProps, "%lf\n",&constrMuLongW[nConstrMuLong-1]) != 1) {
                printf("Error reading constrMuLongW[nConstrMuLong-1]  %lf \n", constrMuLongW[nConstrMuLong-1]);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            if(debugMode) printf("constrMuLongW[nConstrMuLong-1]=%18.15lf\n",constrMuLongW[nConstrMuLong-1]);
            if(wrFilebin && !noConstr) fprintf(fpSolPropsFileBin,"%lf\n",constrMuLongW[nConstrMuLong-1]);
        }
        
        if(fscanf(fpSolProps, "nDegFreedomAtt= %ld\n",&nDegFreedomAtt) != 1) {
            printf("Error reading nDegFreedomAtt  %ld \n",nDegFreedomAtt);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroAtt) nDegFreedomAtt=0;
        printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nDegFreedomAtt= %ld\n", nDegFreedomAtt);
        
        if(fscanf(fpSolProps, "nAttAxes= %hd\n",&nAttAxes) != 1) {
            printf("Error reading nAttAxes  %hd \n",nAttAxes);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroAtt) nAttAxes=0;
        printf("nAttAxes= %7hd\n", nAttAxes);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nAttAxes= %hd\n", nAttAxes);
        
        if(fscanf(fpSolProps, "nFoVs= %d\n",&instrConst[0]) != 1) {
            printf("Error reading nFoVs  %d \n",instrConst[0]);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroInstr) instrConst[0]=0;
        // instrConst[0] must be  0 or 1 because nFoVs = 2 for Gaia and nFoVs=1+instrConst[0]
        if(instrConst[0]<0 || instrConst[0]>1){
            printf("Execution aborted  with invalid nFovs=%d\n",instrConst[0]+1);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
            
        }
        printf("nFovs= %7d ", instrConst[0]+1);
        printf("(instrConst[0])= %7d\n", instrConst[0]);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nFoVs= %d\n", instrConst[0]);
        
        if(fscanf(fpSolProps, "nCCDs= %d\n",&instrConst[1]) != 1) {
            printf("Error reading nCCDs  %d \n",instrConst[1]);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroInstr) instrConst[1]=0;
        printf("nCCDs= %7d\n", instrConst[1]);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nCCDs= %d\n", instrConst[1]);
        
        if(fscanf(fpSolProps, "nPixelColumns= %d\n",&instrConst[2]) != 1) {
            printf("Error reading nPixelColumns  %d \n",instrConst[2]);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroInstr) instrConst[2]=0;
        printf("nPixelColumns= %7d\n", instrConst[2]);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nPixelColumns= %d\n", instrConst[2]);
        
        if(fscanf(fpSolProps, "nTimeIntervals= %d\n",&instrConst[3]) != 1) {
            printf("Error reading nTimeIntervals  %d \n",instrConst[3]);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroInstr) instrConst[3]=0;
        printf("nTimeIntervals= %7d\n", instrConst[3]);
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nTimeIntervals= %d\n", instrConst[3]);
        
        if(fscanf(fpSolProps, "nInstrP= %hd\n",&nInstrP) != 1) {
            printf("Error reading nInstrP  %hd \n",nInstrP);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(zeroInstr) nInstrP=0;
        printf("nInstrP= %7hd\n", nInstrP);
        if(nInstrP!=0 && nInstrP !=6){
            printf("Execution aborted  with invalid nInstrP\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nInstrP= %hd\n", nInstrP);
        if(fscanf(fpSolProps, "nInstrPSolved= %hd\n",&nInstrPSolved) != 1)
        {
            printf("Error reading nInstrPSolved  %hd \n",nInstrPSolved);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("nInstrPSolved= %7hd\n", nInstrPSolved);
        
        if(nInstrPSolved<0 || nInstrPSolved >6){
            printf("Execution aborted  with invalid nInstrPSolved\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"nInstrPSolved= %hd\n", nInstrPSolved);
        
        if(fscanf(fpSolProps, "lsInstrFlag= %d\n",&lsInstrFlag) != 1)
        {
            printf("Error reading lsInstrFlag  %d \n",lsInstrFlag);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        printf("lsInstrFlag= %7d\n", lsInstrFlag);
        
        if(lsInstrFlag<0 || lsInstrFlag >1){
            printf("Execution aborted  with invalid lsInstrFlag\n");
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        if(wrFilebin) fprintf(fpSolPropsFileBin,"lsInstrFlag= %hd\n", lsInstrFlag);
        if(fscanf(fpSolProps, "ssInstrFlag= %d\n",&ssInstrFlag) != 1)
        {
            printf("Error reading ssInstrFlag  %d \n",ssInstrFlag);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        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);
    } //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++)
                {
                    for(jj=0; jj< nAstroPSolved;jj++)
                    {
                        xAstro[(ii-offset)*nAstroP+mapAstroP[jj]] = xSolution[ii*nAstroPSolved+jj];
                    }
                }
                fwrite(xAstro,sizeof(double),(VrIdAstroPDimRecv-offset)*nAstroP,fpAstroR);
                testOnStars+=VrIdAstroPDimRecv-offset;
            }
            free(xAstro);
            
            if(testOnStars!=nStar)
            {
                printf("Error on xAstro testOnStar =%ld not equal to nmStar=%ld\n",testOnStars,nStar);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            testOnStars=0;
            standardErrorAstro = (double *) calloc(VrIdAstroPDimMax*nAstroP, sizeof(double));
            if (!standardErrorAstro)
            {
                printf("Error allocating standardErrorAstro\n");
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            for(ii=0; ii<VrIdAstroPDim; ii++)
            {
                for(jj=0; jj< nAstroPSolved;jj++)
                {
                    standardErrorAstro[ii*nAstroP+mapAstroP[jj]] =standardError[ii*nAstroPSolved+jj];
                }
            }
            fwrite(standardErrorAstro,sizeof(double),VrIdAstroPDim*nAstroP,fpAstroR);
            testOnStars+=VrIdAstroPDim;
            for(kk=1;kk<nproc;kk++)
            {
                offset=0;
                MPI_Recv(&VrIdAstroPDimRecv, 1, MPI_LONG, kk, 12, MPI_COMM_WORLD, &statusMpi);
                MPI_Recv(standardError, VrIdAstroPDimRecv*nAstroPSolved , MPI_DOUBLE, kk, 13,MPI_COMM_WORLD, &statusMpi);
                if(comlsqr.mapStar[kk][0]==comlsqr.mapStar[kk-1][1]) offset=1;
                for(ii=0+offset; ii<VrIdAstroPDimRecv; ii++)
                {
                    for(jj=0; jj< nAstroPSolved;jj++)
                    {
                        standardErrorAstro[(ii-offset)*nAstroP+mapAstroP[jj]] = standardError[ii*nAstroPSolved+jj];
                    }
                }
                fwrite(standardErrorAstro,sizeof(double),(VrIdAstroPDimRecv-offset)*nAstroP,fpAstroR);
                testOnStars+=VrIdAstroPDimRecv-offset;
            }
            free(standardErrorAstro);
            if(testOnStars!=nStar)
            {
                printf("Error on standardErrorAstro testOnStar =%ld not equal to nmStar=%ld\n",testOnStars,nStar);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fclose(fpAstroR);
            //////////////////////////////////////////////////////////
        }
        ////////////////////// writing  GsrAttitudeParamSolution.bin
        fpAttR=fopen(filenameAttResults,"wb");
        if (!fpAttR)
        {
            printf("Error while open %s\n",filenameAttResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        fwrite(&sphereId,sizeof(long),1,fpAttR);
        fwrite(&xSolution[VroffsetAttParam],sizeof(double),nAttParam,fpAttR);
        fwrite(&standardError[VroffsetAttParam],sizeof(double),nAttParam,fpAttR);
        fclose(fpAttR);
        //////////////////////////////////////////////////////////
        
        ////////////////////// writing GsrInstrParamSolution.bin
        fpInstrR=fopen(filenameInstrResults,"wb");
        if (!fpInstrR)
        {
            printf("Error while open %s\n",filenameInstrResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        double *xInstr,*seInstr;
        xInstr=(double *) calloc(nInstrParamTot, sizeof(double));
        seInstr=(double *) calloc(nInstrParamTot, sizeof(double));
        if(nInstrPSolved!=nInstrP)
        {
            int fixedOffsetCMag = nCCDs;
            int fixedOffsetCnu = fixedOffsetCMag+nFoVs*nCCDs;
            int fixedOffsetCdelta_eta = fixedOffsetCnu+nCCDs*nPixelColumns;
            int fixedOffsetCDelta_eta = fixedOffsetCdelta_eta+3*nFoVs*nCCDs*nTimeIntervals;
            int fixedOffsetCdelta_zeta = fixedOffsetCDelta_eta+nCCDs*nPixelColumns;
            
            int counterInstr=0;
            if(maInstrFlag){
                for(int k=0;k<fixedOffsetCMag;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(nuInstrFlag){
                for(int k=fixedOffsetCMag;k<fixedOffsetCnu;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(ssInstrFlag){
                for(int k=fixedOffsetCnu;k<fixedOffsetCdelta_eta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(lsInstrFlag){
                for(int k=fixedOffsetCdelta_eta;k<fixedOffsetCDelta_eta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(ssInstrFlag){
                for(int k=fixedOffsetCDelta_eta;k<fixedOffsetCdelta_zeta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(lsInstrFlag){
                for(int k=fixedOffsetCdelta_zeta;k<nInstrParamTot;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(counterInstr!=nInstrParam) {
                printf("SEVERE ERROR: counterInstr=%ld does not match nInstrParam=%ld while filling in xInstr and seInstr.\n", counterInstr, nInstrParam);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        }
        fwrite(&sphereId,sizeof(long),1,fpInstrR);
        fwrite(xInstr,sizeof(double),nInstrParamTot,fpInstrR);
        fwrite(seInstr,sizeof(double),nInstrParamTot,fpInstrR);
        fclose(fpInstrR);
        free(xInstr);
        free(seInstr);
        //////////////////////////////////////////////////////////
        
        ////////////////////// writing GsrGlobalParamSolution.bin
        fpGlobR=fopen(filenameGlobalResults,"wb");
        if (!fpGlobR)
        {
            printf("Error while open %s\n",filenameGlobalResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        fwrite(&sphereId,sizeof(long),1,fpGlobR);
        fwrite(&xSolution[VroffsetAttParam+nAttParam+nInstrParam],sizeof(double),nGlobalParam,fpGlobR);
        fwrite(&standardError[VroffsetAttParam+nAttParam+nInstrParam],sizeof(double),nGlobalParam,fpGlobR);
        fclose(fpGlobR);
        //////////////////////////////////////////////////////////
        
    } //if (myid==0)
    else
    {
        //////////////// send for AstroResults.bin
        if(nAstroPSolved)
        {
            MPI_Send(&VrIdAstroPDim, 1, MPI_LONG, 0, 10, MPI_COMM_WORLD);
            mpi_send(xSolution, VrIdAstroPDim*nAstroPSolved, MPI_DOUBLE, 0, 11,MPI_COMM_WORLD);
            MPI_Send(&VrIdAstroPDim, 1, MPI_LONG, 0, 12, MPI_COMM_WORLD);
            MPI_Send(standardError, VrIdAstroPDim*nAstroPSolved , MPI_DOUBLE, 0, 13,MPI_COMM_WORLD);
        }
        //////////////////////////////////////////////////////////
    }
    MPI_Barrier(MPI_COMM_WORLD);
    //  tolti i free() mancanti
    free(xSolution);
    free(standardError);
    
    if(myid==0) printf("\nEnd run.\n");
    seconds[3] = time(NULL);
    tot_sec[2] = seconds[3] - seconds[2];
    tot_sec[3] = seconds[3] - seconds[0];
    
    if(myid==0) printf("time before lsqr in sec: %ld\ntime for lsqr: %ld\n", tot_sec[0],
                       tot_sec[1]);
    if(myid==0) printf("time after lsqr in sec: %ld\ntime TOTAL: %ld\n", tot_sec[2],
                       tot_sec[3]);
    
    MPI_Finalize();
    exit(EXIT_SUCCESS);
}