/* 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 #include #include #include #include #include #include "lsqr.h" #include #include #include "util.h" #ifdef OMP #include #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); 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); 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; comlsqr.nthreads= nthreads; comlsqr.ntasks= ntasks; ///////////////////// 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); }