solvergaiaSim.c

/* 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;
    int nthreads=0;
    int ntasks=0;

    ///////////////
    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 =0;// 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 = 5;
        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],"-tasks") == 0)
        {
            ntasks=atoi(argv[i+1]);
            if(ntasks <=0) ntasks=0;
        }
        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(ntasks>0)
            printf(" -tasks %d ",ntasks);
        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);
        }
        printf ("\n");
    }
    /////////////////
    getcwd(wpath, sizePath);
    if (!inputDirOpt)
        strcpy(inputDir, wpath);
    printf("Process %d running on %s\n", myid, processor_name);
////#ifdef OMP
///#pragma omp task
    {
    //    int tid = omp_get_thread_num();
       /*
        const char* s = getenv("NX_SMP_WORKERS");//? "":"1";
        if(s ==NULL)
        	s="1";
        nthreads = atoi(s); //QUI** verifica che la chiamata ritorni il numero di tasks allocati!!
        */
        
        //nthreads sono attuatori
        // ntasks sono la mia suddivisione
        nthreads=omp_get_num_threads();
        if(ntasks==0) 
            ntasks=nthreads;
        printf("PE=%d  on processor %s total num of threads =%d, ntasks=%d\n",myid, processor_name,nthreads,ntasks);
    }
////#pragma omp taskwait
////#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);