diff --git a/GaiaGsrParSim b/GaiaGsrParSim deleted file mode 100755 index 0a6cb5c6896514f89dbcac26744e306710081b61..0000000000000000000000000000000000000000 Binary files a/GaiaGsrParSim and /dev/null differ diff --git a/Makefile b/Makefile index 60a1e3831f0ccf4352609e401332e5cdffbfaf13..aeaed93e7f42726f7f4004e09eb4e2943bd0d3bb 100644 --- a/Makefile +++ b/Makefile @@ -1,89 +1,39 @@ # Gaia GSRPar Makefile -COMPILER = /opt/ompss/bin/mcc -##COMPILER = /opt/openmpi_icc/bin/mpic++ -COMPILERCPP = /opt/ompss/bin/mcxx --ompss -DOMP -##COMPILERCPP = /opt/openmpi_icc/bin/mpic++ +COMPILER = mpiicpc -qopenmp -DOMP +COMPILERCPP = mpiicpc -qopenmp -DOMP -CC = $(COMPILERCPP) +CC = $(COMPILER) CPP= $(COMPILERCPP) -CFITSIOLIB=/opt/cfitsio3_1_0/lib +CFITSIOLIB=/opt/cluster/cfitsio/intel/2017.1/lib/ +GAIAINC=/u/ubecciani/EuroExa/CoreApp7.2/src - - -GAIAINC=/home-volume/ube/Gaia/ParallelCodeV7/src7.2 - - - -GAIAGSR= aprod.o lsqrblas.o lsqr.o solvergaia.o util.o GAIAGSRSIM= aprod.o lsqrblas.o lsqr.o solvergaiaSim.o util.o -GAIAGSRFits2Bin= fits2bin.o util.o lsqrblas.o -GAIAGSRBin2Fits= bin2fits.o -GAIAGSRBin2Reduced= bin2reduced.o util.o lsqrblas.o -GAIAGSRCkEmptyCols= util.o ckemptycols.o lsqrblas.o -GAIAGSRBin2Asc= bin2asc.o -GAIAGSRChTask= changeTask.o -GAIAGSRRepairTask= ripristino.o MEMREQ= memRequest.o - -INCLUDE = -I$(GAIAINC) -I/opt/local/include -I/opt/openmpi_icc/include/ -I/opt/cfitsio3_1_0/include/ -I/opt/openmpi/include -pthread +##INCLUDE = -I$(GAIAINC) -I/opt/local/include -I/opt/local/include/mpich2 -I/opt/cluster/cfitsio/intel/2017.1/include/ +INCLUDE = -I$(GAIAINC) -I/opt/cluster/cfitsio/intel/2017.1/include/ #INCLUDE = -I$(GAIAINC) #CFLAGS= $(INCLUDE) -std=c99 -CPPFLAGS= $(INCLUDE) -g -fsanitize=address -mllvm -asan-stack -CPPFLAGS= $(INCLUDE) -LIB = -pthread -L/usr/local/pbspro/lib -Wl,-rpath -Wl,/usr/local/pbspro/lib -Wl,-rpath -Wl,/opt/openmpi/lib -Wl,--enable-new-dtags -L/opt/openmpi/lib -lmpi -L$(CFITSIOLIB) -lcfitsio -lm +#CPPFLAGS= $(INCLUDE) -g -fsanitize=address -mllvm -asan-stack +#CPPFLAGS= $(INCLUDE) -lirc -limf -lsvml +CPPFLAGS= $(INCLUDE) +LIB = -L$(CFITSIOLIB) -lcfitsio -lm -all: GaiaGsrParSim -####all: GaiaGsrPar MemReq GaiaGsrParSim GaiaFits2Bin GaiaBin2Fits GaiaBin2Reduced GaiaCkEmptyCols GaiaBin2Asc GaiaChTask GaiaRepTask +all: MemReq GaiaGsrParSim ###all: GaiaGsrPar MemReq GaiaGsrParTest GaiaFits2Bin GaiaBin2Fits GaiaBin2Reduced -ckemptycols.o: ckemptycols.cpp - $(CPP) $(CPPFLAGS) -c ckemptycols.cpp - -bin2fits.o: bin2fits.cpp - $(CPP) $(CPPFLAGS) -c bin2fits.cpp - -bin2asc.o: bin2asc.cpp - $(CPP) $(CPPFLAGS) -c bin2asc.cpp - - - -GaiaGsrPar: $(GAIAGSR) - $(CPP) $(CPPFLAGS) -o GaiaGsrPar $(GAIAGSR) $(INCLUDE) $(LIB) - -MemReq: $(MEMREQ) - $(CPP) $(CPPFLAGS) -o MemReq $(MEMREQ) $(INCLUDE) $(LIB) - GaiaGsrParSim: $(GAIAGSRSIM) $(CPP) $(CPPFLAGS) -o GaiaGsrParSim $(GAIAGSRSIM) $(INCLUDE) $(LIB) -GaiaFits2Bin: $(GAIAGSRFits2Bin) - $(CPP) $(CPPFLAGS) -o GaiaFits2Bin $(GAIAGSRFits2Bin) $(INCLUDE) $(LIB) - -GaiaBin2Fits: $(GAIAGSRBin2Fits) - $(CPP) $(CPPFLAGS) -o GaiaBin2Fits $(GAIAGSRBin2Fits) $(INCLUDE) $(LIB) - -GaiaBin2Reduced: $(GAIAGSRBin2Reduced) - $(CPP) $(CPPFLAGS) -o GaiaBin2Reduced $(GAIAGSRBin2Reduced) $(INCLUDE) $(LIB) - -GaiaBin2Asc: $(GAIAGSRBin2Asc) - $(CPP) $(CPPFLAGS) -o GaiaBin2Asc $(GAIAGSRBin2Asc) $(INCLUDE) $(LIB) - - -GaiaCkEmptyCols: $(GAIAGSRCkEmptyCols) - $(CPP) $(CPPFLAGS) -o GaiaCkEmptyCols $(GAIAGSRCkEmptyCols) $(LIB) - -GaiaChTask: $(GAIAGSRChTask) - $(CPP) -o GaiaChTask $(GAIAGSRChTask) $(LIB) +MemReq: $(MEMREQ) + $(CPP) $(CPPFLAGS) -o MemReq $(MEMREQ) $(INCLUDE) $(LIB) -GaiaRepTask: $(GAIAGSRRepairTask) - $(CPP) $(CPPFLAGS) -o GaiaRepairTask $(GAIAGSRRepairTask) $(LIB) - clean: rm -f *.o core + diff --git a/Makefile.marconi_new b/Makefile.marconi_new deleted file mode 100644 index d5ce443830e3162c2f6bcecc2e1d968f1e8cec4b..0000000000000000000000000000000000000000 --- a/Makefile.marconi_new +++ /dev/null @@ -1,38 +0,0 @@ -# Gaia GSRPar Makefile - -COMPILER = mpiicpc -qopenmp -DOMP -COMPILERCPP = mpiicpc -qopenmp -DOMP - -CC = $(COMPILER) -CPP= $(COMPILERCPP) - -CFITSIOLIB=$(CFITSIO_LIB) - -GAIAINC=. - -GAIAGSRSIM= aprod.o lsqrblas.o lsqr.o solvergaiaSim.o util.o -MEMREQ= memRequest.o - - -##INCLUDE = -I$(GAIAINC) -I/opt/local/include -I/opt/local/include/mpich2 -I/opt/cluster/cfitsio/intel/2017.1/include/ -INCLUDE = -I$(GAIAINC) -I$(CFITSIO_INC) -#INCLUDE = -I$(GAIAINC) -#CFLAGS= $(INCLUDE) -std=c99 -#CPPFLAGS= $(INCLUDE) -g -fsanitize=address -mllvm -asan-stack -#CPPFLAGS= $(INCLUDE) -lirc -limf -lsvml -CPPFLAGS= $(INCLUDE) -DOMP= -LIB = -L$(CFITSIOLIB) -lcfitsio -lm - - -all: GaiaGsrParSim -###all: GaiaGsrPar MemReq GaiaGsrParTest GaiaFits2Bin GaiaBin2Fits GaiaBin2Reduced - -GaiaGsrParSim: $(GAIAGSRSIM) - $(CPP) $(CPPFLAGS) -o GaiaGsrParSim $(GAIAGSRSIM) $(INCLUDE) $(LIB) - -MemReq: $(MEMREQ) - $(CPP) $(CPPFLAGS) -o MemReq $(MEMREQ) $(INCLUDE) $(LIB) - -clean: - rm -f *.o core - diff --git a/Makefile_marconi b/Makefile_marconi deleted file mode 100644 index 9e08b9fd3fbdc326de6582ebc376fa7926d01d5a..0000000000000000000000000000000000000000 --- a/Makefile_marconi +++ /dev/null @@ -1,81 +0,0 @@ -# Gaia GSRPar Makefile - -##COMPILER = mpiicpc -g -traceback -O3 -qopenmp -##COMPILERCPP = mpiicpc -g -traceback -O3 -qopenmp -COMPILER = mpiicpc -O3 -qopenmp -COMPILERCPP = mpiicpc -O3 -qopenmp - -CC = $(COMPILER) -CPP = $(COMPILERCPP) - - -# CFITSIOLIB=/opt/cfitsio/lib -#CFITSIOLIB=/cineca/prod/opt/libraries/cfitsio/3.390/intel--pe-xe-2016--binary/lib/ -CFITSIOLIB=$(CFITSIO_LIB) - -GAIAINC=. - - -GAIAGSR= aprod.o lsqrblas.o lsqr.o solvergaia.o util.o -GAIAGSRSIM= aprod.o lsqrblas.o lsqr.o solvergaiaSim.o util.o -GAIAGSRFits2Bin= fits2bin.o util.o lsqrblas.o -GAIAGSRBin2Fits= bin2fits.o -GAIAGSRBin2Reduced= bin2reduced.o util.o lsqrblas.o -GAIAGSRCkEmptyCols= util.o ckemptycols.o lsqrblas.o -GAIAGSRBin2Asc= bin2asc.o -GAIAGSRChTask= changeTask.o -GAIAGSRRepairTask= ripristino.o -MEMREQ= memRequest.o - - -# INCLUDE = -I$(GAIAINC) -I/opt/cfitsio/include -#INCLUDE = -I$(GAIAINC) -I/cineca/prod/opt/libraries/cfitsio/3.390/intel--pe-xe-2016--binary/include/ -INCLUDE = -I$(GAIAINC) -I$(CFITSIO_INC) -CPPFLAGS= $(INCLUDE) -DOMP= -LIB = -L$(CFITSIOLIB) -lcfitsio -lm - -all: GaiaGsrPar MemReq GaiaGsrParSim GaiaFits2Bin GaiaBin2Fits GaiaBin2Reduced GaiaCkEmptyCols GaiaBin2Asc GaiaChTask GaiaRepTask -###all: GaiaGsrPar MemReq GaiaGsrParTest GaiaFits2Bin GaiaBin2Fits GaiaBin2Reduced - -ckemptycols.o: ckemptycols.cpp - $(CPP) $(CPPFLAGS) -c ckemptycols.cpp - -bin2fits.o: bin2fits.cpp - $(CPP) $(CPPFLAGS) -c bin2fits.cpp - -bin2asc.o: bin2asc.cpp - $(CPP) $(CPPFLAGS) -c bin2asc.cpp - -GaiaGsrPar: $(GAIAGSR) - $(CPP) -o GaiaGsrPar $(GAIAGSR) $(INCLUDE) $(LIB) - -MemReq: $(MEMREQ) - $(CPP) -o MemReq $(MEMREQ) $(INCLUDE) $(LIB) - -GaiaGsrParSim: $(GAIAGSRSIM) - $(CPP) -o GaiaGsrParSim $(GAIAGSRSIM) $(INCLUDE) $(LIB) - -GaiaFits2Bin: $(GAIAGSRFits2Bin) - $(CPP) -o GaiaFits2Bin $(GAIAGSRFits2Bin) $(INCLUDE) $(LIB) - -GaiaBin2Fits: $(GAIAGSRBin2Fits) - $(CPP) -o GaiaBin2Fits $(GAIAGSRBin2Fits) $(INCLUDE) $(LIB) - -GaiaBin2Reduced: $(GAIAGSRBin2Reduced) - $(CPP) -o GaiaBin2Reduced $(GAIAGSRBin2Reduced) $(INCLUDE) $(LIB) - -GaiaBin2Asc: $(GAIAGSRBin2Asc) - $(CPP) -o GaiaBin2Asc $(GAIAGSRBin2Asc) $(INCLUDE) $(LIB) - -GaiaChTask: $(GAIAGSRChTask) - $(CPP) -o GaiaChTask $(GAIAGSRChTask) $(LIB) - -GaiaCkEmptyCols: $(GAIAGSRCkEmptyCols) - $(CPP) -o GaiaCkEmptyCols $(GAIAGSRCkEmptyCols) $(LIB) - -GaiaRepTask: $(GAIAGSRRepairTask) - $(CPP) $(CPPFLAGS) -o GaiaRepairTask $(GAIAGSRRepairTask) $(LIB) - -clean: - rm -f *.o core - diff --git a/a b/a deleted file mode 100644 index 8243688b002ac7237467477fdcecfc22f883fccf..0000000000000000000000000000000000000000 --- a/a +++ /dev/null @@ -1 +0,0 @@ -scp *.c *.cpp *.h ubecciani@ssh.hca.bsc.es:/home/ubecciani/AVU-GSR/OMP/src/ diff --git a/aprod.c b/aprod.c index fa389f8e14a934b2106eccc39cd4397b9aedd565..1e1507b99a8fa4a9dde4c0bf9e86345b453d4c1b 100644 --- a/aprod.c +++ b/aprod.c @@ -7,724 +7,458 @@ //#include "pardef.h" #include "util.h" void aprod(int mode, long int m, long int n, double *vVect, double *knownTerms, - double *systemMatrix, long int *matrixIndex, int *instrCol, int *instrConstrIlung, struct comData comlsqr, time_t *ompSec) + double *systemMatrix,long int *matrixIndex, int *instrCol,int *instrConstrIlung, struct comData comlsqr,time_t *ompSec) { - // Parallel definitions - int myid, nproc; +// Parallel definitions + int myid,nproc; long int *mapNoss, *mapNcoeff; - int nthreads, tid, ntasks; - long **mapForThread; - /// struct comData *comlsqr; - // + int nthreads, tid; + long **mapForThread; +/// struct comData *comlsqr; +// - FILE *fk, *fk0; + FILE *fk,*fk0; - double zero = 0.0; - double sum, yi; - long int l1, l2; - long int i, i1; - long int l, j, k; - int i2 = 0, j2 = 0, j3 = 0, na = 0; + double zero=0.0; + double sum, yi; + long int l1,l2; + long int i, i1; + long int l, j, k; + int i2=0,j2=0,j3=0,na=0; int setBound[4]; + + double localSum; + short nAstroPSolved=comlsqr.nAstroPSolved; - double localSum; - short nAstroPSolved = comlsqr.nAstroPSolved; - - long localAstro = comlsqr.VrIdAstroPDim * nAstroPSolved; - long localAstroMax = comlsqr.VrIdAstroPDimMax * nAstroPSolved; - // Initialize. - myid = comlsqr.myid; - nproc = comlsqr.nproc; - mapNcoeff = comlsqr.mapNcoeff; - mapNoss = comlsqr.mapNoss; - nthreads = comlsqr.nthreads; - ntasks= comlsqr.ntasks; + long localAstro= comlsqr.VrIdAstroPDim*nAstroPSolved; + long localAstroMax=comlsqr.VrIdAstroPDimMax*nAstroPSolved; +// Initialize. +myid=comlsqr.myid; +nproc=comlsqr.nproc; +mapNcoeff=comlsqr.mapNcoeff; +mapNoss=comlsqr.mapNoss; +nthreads=comlsqr.nthreads; +mapForThread=comlsqr.mapForThread; +int multMI=comlsqr.multMI; +long nparam=comlsqr.parOss; +short nAttAxes=comlsqr.nAttAxes; +int numOfExtStar=comlsqr.numOfExtStar; +int numOfBarStar=comlsqr.numOfBarStar; +int numOfExtAttCol=comlsqr.numOfExtAttCol; +long VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; +int startingAttColExtConstr=comlsqr.startingAttColExtConstr; +int nOfElextObs=comlsqr.nOfElextObs; +int nEqExtConstr=comlsqr.nEqExtConstr; +int nOfElBarObs=comlsqr.nOfElBarObs; +int nEqBarConstr=comlsqr.nEqBarConstr; +int debugMode=comlsqr.debugMode; +short nInstrPSolved=comlsqr.nInstrPSolved; +int nOfInstrConstr=comlsqr.nOfInstrConstr; +int nElemIC=comlsqr.nElemIC; +short nAttP=comlsqr.nAttP; +short nGlobP=comlsqr.nGlobP; + +setBound[0]=comlsqr.setBound[0]; + setBound[1]=comlsqr.setBound[1]; + setBound[2]=comlsqr.setBound[2]; + setBound[3]=comlsqr.setBound[3]; + long nDegFreedomAtt=comlsqr.nDegFreedomAtt; + short nAttParAxis=comlsqr.nAttParAxis; + long offsetAttParam=comlsqr.offsetAttParam; + long offsetInstrParam=comlsqr.offsetInstrParam; + long offsetGlobParam=comlsqr.offsetGlobParam; + +nthreads=1; +tid=0; + FILE *fp1,*fp2; +// fp1=fopen("test1_aprod","w"); +// fp2=fopen("test2_aprod","w"); - mapForThread = comlsqr.mapForThread; - int multMI = comlsqr.multMI; - long nparam = comlsqr.parOss; - short nAttAxes = comlsqr.nAttAxes; - int numOfExtStar = comlsqr.numOfExtStar; - int numOfBarStar = comlsqr.numOfBarStar; - int numOfExtAttCol = comlsqr.numOfExtAttCol; - long VrIdAstroPDimMax = comlsqr.VrIdAstroPDimMax; - int startingAttColExtConstr = comlsqr.startingAttColExtConstr; - int nOfElextObs = comlsqr.nOfElextObs; - int nEqExtConstr = comlsqr.nEqExtConstr; - int nOfElBarObs = comlsqr.nOfElBarObs; - int nEqBarConstr = comlsqr.nEqBarConstr; - int debugMode = comlsqr.debugMode; - short nInstrPSolved = comlsqr.nInstrPSolved; - int nOfInstrConstr = comlsqr.nOfInstrConstr; - int nElemIC = comlsqr.nElemIC; - short nAttP = comlsqr.nAttP; - short nGlobP = comlsqr.nGlobP; - - setBound[0] = comlsqr.setBound[0]; - setBound[1] = comlsqr.setBound[1]; - setBound[2] = comlsqr.setBound[2]; - setBound[3] = comlsqr.setBound[3]; - long nDegFreedomAtt = comlsqr.nDegFreedomAtt; - short nAttParAxis = comlsqr.nAttParAxis; - long offsetAttParam = comlsqr.offsetAttParam; - long offsetInstrParam = comlsqr.offsetInstrParam; - long offsetGlobParam = comlsqr.offsetGlobParam; - //nthreads = 1; - tid = 0; - FILE *fp1, *fp2; - // fp1=fopen("test1_aprod","w"); - // fp2=fopen("test2_aprod","w"); - if (mode != 1 && mode != 2) - { - printf("ERROR: Invalid mode=%d in aprod function\n", mode); - exit(1); - } - l1 = 0; - l2 = 0; - myid = comlsqr.myid; - if (mode == 1) +if(mode!=1 && mode !=2) +{ + printf("ERROR: Invalid mode=%d in aprod function\n",mode); + exit(1); +} +l1=0; +l2=0; +myid=comlsqr.myid; +if(mode==1) +{ + time_t startTime=time(NULL); +#pragma omp parallel private(myid, sum, k, l1, l2, l, j,tid,nthreads,i2,na) shared(mapNoss,instrCol,comlsqr,vVect,systemMatrix,matrixIndex,knownTerms,j2) { - time_t startTime = time(NULL); -//// #pragma omp parallel private(myid, sum, k, l1, l2, l, j, tid, nthreads, i2, na) shared(mapNoss, instrCol, comlsqr, vVect, systemMatrix, matrixIndex, knownTerms, j2) + myid=comlsqr.myid; + if(comlsqr.itn==1) { - myid = comlsqr.myid; - /* - //FV_ EDIT ompSs - if (comlsqr.itn == 1) - { - #ifdef OMP - tid = omp_get_thread_num(); - nthreads = omp_get_num_threads(); - #endif - }*/ - if (comlsqr.itn == 1 && debugMode) - printf("PE=%d Aprod1 OpenMP num of threads =%d from thread =%d icycle=%ld comlsqr.itn=%d\n", myid, nthreads, tid, i, comlsqr.itn); - long miValAstro = 0; - long miValAtt = 0; - long jstartAtt = 0; - long jstartAstro = 0; - long lset = 0; - long offLocalAstro = 0; - long offLocalAtt = 0; - long offLocalInstr = 0; //Offset on Instruments - long ixInstr = 0; - int nInstrVal = 0; - offLocalInstr = offsetInstrParam + (localAstroMax - offsetAttParam); //Offset on Instruments - nInstrVal = nAstroPSolved + nAttP; - offLocalAstro = comlsqr.mapStar[myid][0] * nAstroPSolved; //Offset on my mstars - offLocalAtt = localAstroMax - offsetAttParam; //Offset on attitude - long offLocalGlob = offsetGlobParam + (localAstroMax - offsetAttParam); //Offset on GlobP - int nGlobVal = nAstroPSolved + nAttP + nInstrPSolved; - jstartAstro = miValAstro - offLocalAstro; - - //FV_ EDIT ompSs - - - for(int nt=0; nt < ntasks; nt++ ) - { - #pragma omp task //shared(matrixIndex,mapForThread,vVect,systemMatrix,knownTerms) in(nt,multMI,offLocalAstro, lset) - { - for (long ix = mapForThread[nt][0]; ix < mapForThread[nt][2]; ix++) - { - sum = 0.; - ///////////////////////////////////////////////////// - /// Mode 1 Astrometric Sect - if (nAstroPSolved) - { - - lset = ix * nparam; - if (matrixIndex[multMI * ix] != miValAstro) - { - miValAstro = matrixIndex[multMI * ix]; - jstartAstro = miValAstro - offLocalAstro; - } - for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++) - { - sum = sum + systemMatrix[lset] * vVect[jx]; - lset++; - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Attitude Sect - if (nAttP) - { - lset = ix * nparam + nAstroPSolved; - miValAtt = matrixIndex[multMI * ix + (multMI - 1)]; - for (int nax = 0; nax < nAttAxes; nax++) - { - jstartAtt = miValAtt + offLocalAtt + nax * nDegFreedomAtt; - for (long inpax = jstartAtt; inpax < jstartAtt + nAttParAxis; inpax++) - { - sum = sum + systemMatrix[lset] * vVect[inpax]; - lset++; - } - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Instrument Sect - if (nInstrPSolved) - { - - lset = ix * nparam + nInstrVal; - long iiVal = ix * nInstrPSolved; - for (int inInstr = 0; inInstr < nInstrPSolved; inInstr++) - { - ixInstr = offLocalInstr + instrCol[iiVal + inInstr]; - sum = sum + systemMatrix[lset] * vVect[ixInstr]; - lset++; - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Global sect - if (nGlobP) - { - lset = ix * nparam + nGlobVal; - for (long inGlob = offLocalGlob; inGlob < offLocalGlob + nGlobP; inGlob++) - { - sum = sum + systemMatrix[lset] * vVect[inGlob]; - lset++; - } - } - ////////////////////////////////////////////////////// - knownTerms[ix] += sum; - } //for ix - } - #pragma omp taskwait - - } - - - - - - - - /* - //#pragma omp for - for (long ix = 0; ix < mapNoss[myid]; ix++) - { - sum = 0.; - ///////////////////////////////////////////////////// - /// Mode 1 Astrometric Sect - if (nAstroPSolved) - { - - lset = ix * nparam; - if (matrixIndex[multMI * ix] != miValAstro) - { - miValAstro = matrixIndex[multMI * ix]; - jstartAstro = miValAstro - offLocalAstro; - } - for (long jx = jstartAstro; jx < jstartAstro + nAstroPSolved; jx++) - { - sum = sum + systemMatrix[lset] * vVect[jx]; - lset++; - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Attitude Sect - if (nAttP) - { - lset = ix * nparam + nAstroPSolved; - miValAtt = matrixIndex[multMI * ix + (multMI - 1)]; - for (int nax = 0; nax < nAttAxes; nax++) - { - jstartAtt = miValAtt + offLocalAtt + nax * nDegFreedomAtt; - for (long inpax = jstartAtt; inpax < jstartAtt + nAttParAxis; inpax++) - { - sum = sum + systemMatrix[lset] * vVect[inpax]; - lset++; - } - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Instrument Sect - if (nInstrPSolved) - { +#ifdef OMP + tid = omp_get_thread_num(); + nthreads = omp_get_num_threads(); +#endif + } + if(comlsqr.itn==1 && debugMode) + printf("PE=%d Aprod1 OpenMP num of threads =%d from thread =%d icycle=%ld comlsqr.itn=%d\n",myid, nthreads,tid,i,comlsqr.itn); + long miValAstro=0; + long miValAtt=0; + long jstartAtt=0; + long jstartAstro=0; + long lset=0; + long offLocalAstro=0; + long offLocalAtt=0; + long offLocalInstr=0; //Offset on Instruments + long ixInstr=0; + int nInstrVal=0; + offLocalInstr=offsetInstrParam+(localAstroMax-offsetAttParam); //Offset on Instruments + nInstrVal=nAstroPSolved+nAttP; + offLocalAstro=comlsqr.mapStar[myid][0]*nAstroPSolved; //Offset on my mstars + offLocalAtt=localAstroMax-offsetAttParam; //Offset on attitude + long offLocalGlob=offsetGlobParam+(localAstroMax-offsetAttParam); //Offset on GlobP + int nGlobVal=nAstroPSolved+nAttP+nInstrPSolved; + jstartAstro=miValAstro-offLocalAstro; - lset = ix * nparam + nInstrVal; - long iiVal = ix * nInstrPSolved; - for (int inInstr = 0; inInstr < nInstrPSolved; inInstr++) - { - ixInstr = offLocalInstr + instrCol[iiVal + inInstr]; - sum = sum + systemMatrix[lset] * vVect[ixInstr]; - lset++; - } - } - ////////////////////////////////////////////////////// - /// Mode 1 Global sect - if (nGlobP) - { - lset = ix * nparam + nGlobVal; - for (long inGlob = offLocalGlob; inGlob < offLocalGlob + nGlobP; inGlob++) - { - sum = sum + systemMatrix[lset] * vVect[inGlob]; - lset++; - } - } - ////////////////////////////////////////////////////// - knownTerms[ix] += sum; - } //for ix - */ - - - - - /// Mode 1 ExtConstr - if (nEqExtConstr) - { - long offExtAtt; - long offExtAttConstr = VrIdAstroPDimMax * nAstroPSolved + startingAttColExtConstr; - long vVIx; - long ktIx = mapNoss[myid]; - long offExtConstr; - for (int iexc = 0; iexc < nEqExtConstr; iexc++) - { - sum = 0.0; - offExtConstr = mapNcoeff[myid] + iexc * nOfElextObs; - //FV_ EDIT ompSs - //#pragma omp for - for (int j3 = 0; j3 < numOfExtStar * nAstroPSolved; j3++) - sum += systemMatrix[offExtConstr + j3] * vVect[j3]; - for (int nax = 0; nax < nAttAxes; nax++) - { - offExtAtt = offExtConstr + numOfExtStar * nAstroPSolved + nax * numOfExtAttCol; - vVIx = offExtAttConstr + nax * nDegFreedomAtt; - //FV_ EDIT ompSs - //#pragma omp for - for (int j3 = 0; j3 < numOfExtAttCol; j3++) - sum += systemMatrix[offExtAtt + j3] * vVect[vVIx + j3]; - } - //FV_ EDIT ompSs - //#pragma omp atomic - knownTerms[ktIx + iexc] += sum; - } //for iexc - } - ////////////////////////////////////////////////////// - /// Mode 1 BarConstr - if (nEqBarConstr) - { - long offBarConstr = mapNcoeff[myid] + nOfElextObs * nEqExtConstr; - long offBarConstrIx; - long ktIx = mapNoss[myid] + nEqExtConstr; - for (int iexc = 0; iexc < nEqBarConstr; iexc++) - { - sum = 0.0; - offBarConstrIx = offBarConstr + iexc * nOfElBarObs; - //FV_ EDIT ompSs - //#pragma omp for - for (int j3 = 0; j3 < numOfBarStar * nAstroPSolved; j3++) - sum += systemMatrix[offBarConstrIx + j3] * vVect[j3]; - //FV_ EDIT ompSs - //#pragma omp atomic - knownTerms[ktIx + iexc] += sum; - } //for iexc - } - ////////////////////////////////////////////////////// - /// Mode 1 InstrConstr - if (nOfInstrConstr) +#pragma omp for + for(long ix=0;ix ZERO, - wantse = standardError != NULL; -long int i; -int - maxdx, - nconv, nstop; -double - alfopt, - alpha, arnorm0, beta, bnorm, - cs, cs1, cs2, ctol, - delta, dknorm, dknormSum, dnorm, dxk, dxmax, - gamma, gambar, phi, phibar, psi, - res2, rho, rhobar, rhbar1, - rhs, rtol, sn, sn1, sn2, - t, tau, temp, test1, test2, test3, - theta, t1, t2, t3, xnorm1, z, zbar; -char - enter[] = "Enter LSQR. ", - exitLsqr[] = "Exit LSQR. ", - msg[6][100] = + const bool + extra = false, // true for extra printing below. + damped = damp > ZERO, + wantse = standardError != NULL; + long int i; + int + maxdx, nconv, nstop; + double + alfopt, alpha, arnorm0, beta, bnorm, + cs, cs1, cs2, ctol, + delta, dknorm, dknormSum, dnorm, dxk, dxmax, + gamma, gambar, phi, phibar, psi, + res2, rho, rhobar, rhbar1, + rhs, rtol, sn, sn1, sn2, + t, tau, temp, test1, test2, test3, + theta, t1, t2, t3, xnorm1, z, zbar; + char + enter[] = "Enter LSQR. ", + exitLsqr[] = "Exit LSQR. ", + msg[6][100] = { {"The exact solution is x = 0"}, {"A solution to Ax = b was found, given atol, btol"}, {"A least-squares solution was found, given atol"}, {"A damped least-squares solution was found, given atol"}, {"Cond(Abar) seems to be too large, given conlim"}, - {"The iteration limit was reached"}}; -char lsqrOut[] = "lsqr-output"; // (buffer flush problem) -char wpath[1024]; -size_t sizePath = 1020; - + {"The iteration limit was reached"} + }; + char lsqrOut[] = "lsqr-output"; // (buffer flush problem) + char wpath[1024]; + size_t sizePath=1020; + //----------------------------------------------------------------------- // Format strings. -char fmt_1000[] = - " %s Least-squares solution of Ax = b\n" - " The matrix A has %ld rows and %ld columns\n" - " damp = %-22.2e wantse = %10i\n" - " atol = %-22.2e conlim = %10.2e\n" - " btol = %-22.2e itnlim = %10d\n\n"; -char fmt_1200[] = - " Itn x(1) Function" - " Compatible LS Norm A Cond A\n"; -char fmt_1300[] = - " Itn x(1) Function" - " Compatible LS Norm Abar Cond Abar\n"; -char fmt_1400[] = - " phi dknorm dxk alfa_opt\n"; -char fmt_1500_extra[] = - " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e %8.1e %7.1e %7.1e %7.1e\n"; -char fmt_1500[] = - " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e\n"; -char fmt_1550[] = - " %6d %16.9e %16.9e %9.2e %9.2e\n"; -char fmt_1600[] = - "\n"; -char fmt_2000[] = - "\n" - " %s istop = %-10d itn = %-10d\n" - " %s anorm = %11.5e acond = %11.5e\n" - " %s vnorm = %11.5e xnorm = %11.5e\n" - " %s rnorm = %11.5e arnorm = %11.5e\n"; -char fmt_2100[] = - " %s max dx = %7.1e occured at itn %-9d\n" - " %s = %7.1e*xnorm\n"; -char fmt_3000[] = - " %s %s\n"; + char fmt_1000[] = + " %s Least-squares solution of Ax = b\n" + " The matrix A has %ld rows and %ld columns\n" + " damp = %-22.2e wantse = %10i\n" + " atol = %-22.2e conlim = %10.2e\n" + " btol = %-22.2e itnlim = %10d\n\n"; + char fmt_1200[] = + " Itn x(1) Function" + " Compatible LS Norm A Cond A\n"; + char fmt_1300[] = + " Itn x(1) Function" + " Compatible LS Norm Abar Cond Abar\n"; + char fmt_1400[] = + " phi dknorm dxk alfa_opt\n"; + char fmt_1500_extra[] = + " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e %8.1e %7.1e %7.1e %7.1e\n"; + char fmt_1500[] = + " %6d %16.9e %16.9e %9.2e %9.2e %8.1e %8.1e\n"; + char fmt_1550[] = + " %6d %16.9e %16.9e %9.2e %9.2e\n"; + char fmt_1600[] = + "\n"; + char fmt_2000[] = + "\n" + " %s istop = %-10d itn = %-10d\n" + " %s anorm = %11.5e acond = %11.5e\n" + " %s vnorm = %11.5e xnorm = %11.5e\n" + " %s rnorm = %11.5e arnorm = %11.5e\n"; + char fmt_2100[] = + " %s max dx = %7.1e occured at itn %-9d\n" + " %s = %7.1e*xnorm\n"; + char fmt_3000[] = + " %s %s\n"; ///////////// Specific definitions -time_t startCycleTime, endCycleTime, totTime, partialTime, CPRtimeStart, CPRtimeEnd, ompSec = 0; -int writeCPR; //=0 no write CPR, =1 write CPR and continue, =2 write CPR and return -int noCPR; -int myid, nproc; -long int *mapNoss, *mapNcoeff; -MPI_Status status; -long int nunkSplit, localAstro, localAstroMax, other; -int nAstroElements; -int debugMode; -long VrIdAstroPDim, VrIdAstroPDimMax; -long nDegFreedomAtt; -int nAttAxes; -int nAttParam, nInstrParam, nGlobalParam, itnTest, nAttP, numOfExtStar, numOfBarStar, numOfExtAttCol, nobs; -long mapNossBefore; -short nAstroPSolved, nInstrPsolved; -double cycleStartMpiTime, cycleEndMpiTime; -double *vAuxVect; -int CPRCount; -FILE *fpCPRend, *nout, *fpItnLimitNew; -int nEqExtConstr, nEqBarConstr; -int nElemIC, nOfInstrConstr; -double *kAuxcopy; -double *kcopy; -int ntasks; -int nthreads; -//////////////////////////////// + time_t startCycleTime,endCycleTime,totTime,partialTime,CPRtimeStart,CPRtimeEnd,ompSec=0; + int writeCPR; //=0 no write CPR, =1 write CPR and continue, =2 write CPR and return + int noCPR; + int myid,nproc; + long int *mapNoss, *mapNcoeff; + MPI_Status status; + long int nunkSplit, localAstro, localAstroMax, other; + int nAstroElements; + int debugMode; + long VrIdAstroPDim, VrIdAstroPDimMax; + long nDegFreedomAtt; + int nAttAxes; + int nAttParam,nInstrParam,nGlobalParam,itnTest,nAttP,numOfExtStar,numOfBarStar,numOfExtAttCol,nobs; + long mapNossBefore; + short nAstroPSolved,nInstrPsolved; + double cycleStartMpiTime, cycleEndMpiTime; + double *vAuxVect; + int CPRCount; + FILE *fpCPRend, *nout, *fpItnLimitNew; + int nEqExtConstr,nEqBarConstr; + int nElemIC,nOfInstrConstr; + double *kAuxcopy; + double *kcopy; +//////////////////////////////// // Initialize. -myid = comlsqr.myid; -nproc = comlsqr.nproc; -mapNcoeff = comlsqr.mapNcoeff; -mapNoss = comlsqr.mapNoss; -nAttParam = comlsqr.nAttParam; -nInstrParam = comlsqr.nInstrParam; -nGlobalParam = comlsqr.nGlobalParam; -nunkSplit = comlsqr.nunkSplit; -VrIdAstroPDim = comlsqr.VrIdAstroPDim; -VrIdAstroPDimMax = comlsqr.VrIdAstroPDimMax; -nAstroPSolved = comlsqr.nAstroPSolved; -nEqExtConstr = comlsqr.nEqExtConstr; -nEqBarConstr = comlsqr.nEqBarConstr; -nDegFreedomAtt = comlsqr.nDegFreedomAtt; -nAttAxes = comlsqr.nAttAxes; -localAstro = VrIdAstroPDim * nAstroPSolved; -localAstroMax = VrIdAstroPDimMax * nAstroPSolved; -numOfExtStar = comlsqr.numOfExtStar; -numOfBarStar = comlsqr.numOfBarStar; -nAttP = comlsqr.nAttP; -numOfExtAttCol = comlsqr.numOfExtAttCol; -mapNossBefore = comlsqr.mapNossBefore; -nobs = comlsqr.nobs; -debugMode = comlsqr.debugMode; -noCPR = comlsqr.noCPR; -nElemIC = comlsqr.nElemIC; -nOfInstrConstr = comlsqr.nOfInstrConstr; -nInstrPsolved = comlsqr.nInstrPSolved; -ntasks= comlsqr.ntasks; -nthreads= comlsqr.nthreads; - - other = (long)nAttParam + nInstrParam + nGlobalParam; -if (nAstroPSolved) - vAuxVect = (double *)calloc(localAstroMax, sizeof(double)); - -comlsqr.itn = itn; -totTime = comlsqr.totSec; -partialTime = comlsqr.totSec; - -if (debugMode) -{ - for (long i = 0; i < mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr; i++) - { - printf("PE=%d Knowterms[%d]=%e\n", myid, i, knownTerms[i]); +myid=comlsqr.myid; +nproc=comlsqr.nproc; +mapNcoeff=comlsqr.mapNcoeff; +mapNoss=comlsqr.mapNoss; +nAttParam=comlsqr.nAttParam; +nInstrParam=comlsqr.nInstrParam; +nGlobalParam=comlsqr.nGlobalParam; +nunkSplit=comlsqr.nunkSplit; +VrIdAstroPDim=comlsqr.VrIdAstroPDim; +VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; +nAstroPSolved=comlsqr.nAstroPSolved; +nEqExtConstr=comlsqr.nEqExtConstr; +nEqBarConstr=comlsqr.nEqBarConstr; +nDegFreedomAtt=comlsqr.nDegFreedomAtt; +nAttAxes=comlsqr.nAttAxes; +localAstro= VrIdAstroPDim*nAstroPSolved; +localAstroMax=VrIdAstroPDimMax*nAstroPSolved; +numOfExtStar=comlsqr.numOfExtStar; +numOfBarStar=comlsqr.numOfBarStar; +nAttP=comlsqr.nAttP; +numOfExtAttCol=comlsqr.numOfExtAttCol; +mapNossBefore=comlsqr.mapNossBefore; +nobs=comlsqr.nobs; +debugMode=comlsqr.debugMode; +noCPR=comlsqr.noCPR; +nElemIC=comlsqr.nElemIC; +nOfInstrConstr=comlsqr.nOfInstrConstr; +nInstrPsolved=comlsqr.nInstrPSolved; + +other=(long)nAttParam + nInstrParam + nGlobalParam; +if(nAstroPSolved) vAuxVect=(double *) calloc(localAstroMax,sizeof(double)); + +comlsqr.itn=itn; +totTime=comlsqr.totSec; +partialTime=comlsqr.totSec; + + if(debugMode){ + for(long i=0; i ZERO) - ctol = ONE / conlim; -anorm = ZERO; -acond = ZERO; -dnorm = ZERO; -dxmax = ZERO; -res2 = ZERO; -psi = ZERO; -xnorm = ZERO; -xnorm1 = ZERO; -cs2 = -ONE; -sn2 = ZERO; -z = ZERO; - -CPRCount = 0; +itn = 0; +istop = 0; +nstop = 0; +maxdx = 0; +ctol = ZERO; +if (conlim > ZERO) ctol = ONE / conlim; +anorm = ZERO; +acond = ZERO; +dnorm = ZERO; +dxmax = ZERO; +res2 = ZERO; +psi = ZERO; +xnorm = ZERO; +xnorm1 = ZERO; +cs2 = - ONE; +sn2 = ZERO; +z = ZERO; + +CPRCount=0; // ------------------------------------------------------------------ // Set up the first vectors u and v for the bidiagonalization. // These satisfy beta*u = b, alpha*v = A(transpose)*u. // ------------------------------------------------------------------ -dload(nunkSplit, 0.0, vVect); -dload(nunkSplit, 0.0, xSolution); -if (wantse) - dload(nunkSplit, 0.0, standardError); -alpha = ZERO; +dload( nunkSplit, 0.0, vVect ); +dload( nunkSplit, 0.0, xSolution ); +if ( wantse ) + dload( nunkSplit, 0.0, standardError ); +alpha = ZERO; // Find the maximum value on u -double betaLoc; -if (myid == 0) - betaLoc = cblas_dnrm2(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, knownTerms, 1); -else - betaLoc = cblas_dnrm2(mapNoss[myid], knownTerms, 1); + double betaLoc; + if(myid==0) betaLoc=cblas_dnrm2(mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr,knownTerms,1); + else betaLoc=cblas_dnrm2(mapNoss[myid],knownTerms,1); -double betaLoc2 = betaLoc * betaLoc; +double betaLoc2=betaLoc*betaLoc; double betaGlob; -MPI_Allreduce(&betaLoc2, &betaGlob, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -beta = sqrt(betaGlob); -if (beta > ZERO) +MPI_Allreduce(&betaLoc2, &betaGlob,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); +beta=sqrt(betaGlob); +if (beta > ZERO) { - cblas_dscal(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, (ONE / beta), knownTerms, 1); - MPI_Barrier(MPI_COMM_WORLD); + cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (ONE / beta), knownTerms, 1 ); + MPI_Barrier(MPI_COMM_WORLD); //only processor 0 accumulate and distribute the v array -//// - - // FV EDIT - //#pragma omp for - for (i = 0; i < localAstroMax; i++) - { - vAuxVect[i] = vVect[i]; - vVect[i] = 0; - } - if (myid != 0) +#pragma omp for + for(i=0;i ZERO) { - double alphaOther = cblas_dnrm2(comlsqr.nunkSplit - localAstroMax, &vVect[localAstroMax], 1); - double alphaOther2 = alphaOther * alphaOther; - alphaLoc2 = alphaLoc2 + alphaOther2; + cblas_dscal ( nunkSplit, (ONE / alpha), vVect, 1 ); + cblas_dcopy ( nunkSplit, vVect, 1, wVect, 1 ); } - double alphaGlob2 = 0; - MPI_Allreduce(&alphaLoc2, &alphaGlob2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - alpha = sqrt(alphaGlob2); -} - -if (alpha > ZERO) -{ - cblas_dscal(nunkSplit, (ONE / alpha), vVect, 1); - cblas_dcopy(nunkSplit, vVect, 1, wVect, 1); -} // printf("LSQR T4: PE=%d alpha=%15.12lf beta=%15.12lf arnorm=%15.12lf\n",myid,alpha,beta,arnorm); -arnorm = alpha * beta; -arnorm0 = arnorm; + arnorm = alpha * beta; + arnorm0 = arnorm; // printf("LSQR T5: PE=%d alpha=%15.12lf beta=%15.12lf arnorm=%15.12lf\n",myid,alpha,beta,arnorm); -if (arnorm == ZERO) - goto goto_800; -rhobar = alpha; -phibar = beta; -bnorm = beta; -rnorm = beta; + if (arnorm == ZERO) goto goto_800; + rhobar = alpha; + phibar = beta; + bnorm = beta; + rnorm = beta; -if (myid == 0) -{ - nout = fopen(lsqrOut, "a"); - if (nout != NULL) - { - if (damped) - fprintf(nout, fmt_1300); - else - fprintf(nout, fmt_1200); - } - else + if (myid==0) { - printf("Error while opening %s (2)\n", lsqrOut); - MPI_Abort(MPI_COMM_WORLD, 1); - exit(EXIT_FAILURE); + nout=fopen(lsqrOut,"a"); + if ( nout != NULL) { + if ( damped ) + fprintf(nout, fmt_1300); + else + fprintf(nout, fmt_1200); + } else { + printf("Error while opening %s (2)\n",lsqrOut); + MPI_Abort(MPI_COMM_WORLD,1); + exit(EXIT_FAILURE); + } + + test1 = ONE; + test2 = alpha / beta; + + if ( extra ) + fprintf(nout, fmt_1400); + + fprintf(nout, fmt_1550, itn, xSolution[0], rnorm, test1, test2); + fprintf(nout, fmt_1600); + fclose(nout); } - test1 = ONE; - test2 = alpha / beta; - - if (extra) - fprintf(nout, fmt_1400); - - fprintf(nout, fmt_1550, itn, xSolution[0], rnorm, test1, test2); - fprintf(nout, fmt_1600); - fclose(nout); -} - // ================================================================== // CheckPointRestart // ================================================================== - -if (!noCPR) - if (myid == 0) - printf("PE=%d call restart setup\n", myid); -restartSetup(&itn, - knownTerms, - &beta, - &alpha, - vVect, - &anorm, - &rhobar, - &phibar, - wVect, - xSolution, - standardError, - &dnorm, - &sn2, - &cs2, - &z, - &xnorm1, - &res2, - &nstop, - comlsqr); -if (myid == 0) - printf("PE=%d end restart setup\n", myid); + + if(!noCPR) + if(myid==0) printf("PE=%d call restart setup\n",myid); + restartSetup( &itn, + knownTerms, + &beta, + &alpha, + vVect, + &anorm, + &rhobar, + &phibar, + wVect, + xSolution, + standardError, + &dnorm, + &sn2, + &cs2, + &z, + &xnorm1, + &res2, + &nstop, + comlsqr); + if(myid==0) printf("PE=%d end restart setup\n",myid); // ================================================================== // Main iteration loop. // ================================================================== -if (myid == 0) -{ - startCycleTime = time(NULL); -} -//////////////////////// START ITERATIONS -while (1) -{ - writeCPR = 0; // No write CPR - cycleStartMpiTime = MPI_Wtime(); - itnTest = -1; - if (myid == 0) - { - fpItnLimitNew = NULL; - fpItnLimitNew = fopen("ITNLIMIT_NEW", "r"); - if (fpItnLimitNew != NULL) - { - fscanf(fpItnLimitNew, "%d\n", &itnTest); - if (itnTest > 0) + if(myid==0) + { + startCycleTime=time(NULL); + } +//////////////////////// START ITERATIONS + while (1) { + writeCPR=0; // No write CPR + cycleStartMpiTime=MPI_Wtime(); + itnTest=-1; + if(myid==0){ + fpItnLimitNew=NULL; + fpItnLimitNew=fopen("ITNLIMIT_NEW","r"); + if (fpItnLimitNew!=NULL) + { + fscanf(fpItnLimitNew, "%d\n",&itnTest); + if(itnTest>0) { - comlsqr.itnLimit = itnTest; - if (itnlim != itnTest) - { - itnlim = itnTest; - if (myid == 0) - printf("itnLimit forced with ITNLIMIT_NEW file to value %d\n", itnlim); + comlsqr.itnLimit=itnTest; + if(itnlim != itnTest){ + itnlim=itnTest; + if(myid==0) printf("itnLimit forced with ITNLIMIT_NEW file to value %d\n",itnlim); } } fclose(fpItnLimitNew); + } } - } - - if (myid == 0) - { - endCycleTime = time(NULL) - startCycleTime; - startCycleTime = time(NULL); - totTime = totTime + endCycleTime; - partialTime = partialTime + endCycleTime; - if (!noCPR) + + if(myid==0) { - - if (partialTime > comlsqr.timeCPR * 60) + endCycleTime=time(NULL)-startCycleTime; + startCycleTime=time(NULL); + totTime=totTime+endCycleTime; + partialTime=partialTime+endCycleTime; + if(!noCPR){ + + if(partialTime>comlsqr.timeCPR*60) + { + writeCPR=1; + partialTime=0; + } + if(totTime+endCycleTime*2>comlsqr.timeLimit*60) writeCPR=2; + + if(CPRCount>0 && CPRCount%comlsqr.itnCPR==0) writeCPR=1; + if(CPRCount>0 && CPRCount==comlsqr.itnCPRstop) writeCPR=2; + if(CPRCount==itnlim) writeCPR=2; + fpCPRend=NULL; + fpCPRend=fopen("CPR_END","r"); + if (fpCPRend!=NULL) { - writeCPR = 1; - partialTime = 0; - } - if (totTime + endCycleTime * 2 > comlsqr.timeLimit * 60) - writeCPR = 2; - - if (CPRCount > 0 && CPRCount % comlsqr.itnCPR == 0) - writeCPR = 1; - if (CPRCount > 0 && CPRCount == comlsqr.itnCPRstop) - writeCPR = 2; - if (CPRCount == itnlim) - writeCPR = 2; - fpCPRend = NULL; - fpCPRend = fopen("CPR_END", "r"); - if (fpCPRend != NULL) - { - itnTest = -1; - fscanf(fpCPRend, "%d\n", &itnTest); - if (itnTest > 0) + itnTest=-1; + fscanf(fpCPRend, "%d\n",&itnTest); + if(itnTest>0) { - if (comlsqr.itnCPRend != itnTest) - { - comlsqr.itnCPRend = itnTest; - printf("itnCPRend forced with CPR_END file to value %d\n", comlsqr.itnCPRend); + if(comlsqr.itnCPRend != itnTest){ + comlsqr.itnCPRend=itnTest; + printf("itnCPRend forced with CPR_END file to value %d\n",comlsqr.itnCPRend); } } fclose(fpCPRend); } - if (comlsqr.itnCPRend > 0 && comlsqr.itnCPRend <= itn) + if(comlsqr.itnCPRend>0 && comlsqr.itnCPRend<=itn) { - writeCPR = 2; - printf("itnCPRend condition triggered to value %d. writeCPR set to 2.\n", comlsqr.itnCPRend); + writeCPR=2; + printf("itnCPRend condition triggered to value %d. writeCPR set to 2.\n",comlsqr.itnCPRend); } - } //if(!noCPR) - } //if(myid==0) - MPI_Barrier(MPI_COMM_WORLD); - MPI_Bcast(&writeCPR, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&itnlim, 1, MPI_INT, 0, MPI_COMM_WORLD); - MPI_Bcast(&comlsqr.itnLimit, 1, MPI_INT, 0, MPI_COMM_WORLD); - if (myid == 0) - { - printf("lsqr: Iteration number %d. Iteration seconds %ld. OmpSec %d Global Seconds %ld \n", itn, endCycleTime, ompSec, totTime); - if (writeCPR == 1) - printf("... writing CPR files\n"); - if (writeCPR == 2) - printf("... writing CPR files and return\n"); - } - - if (writeCPR > 0) - { - if (myid == 0) - CPRtimeStart = time(NULL); - writeCheckPoint(itn, - knownTerms, - beta, - alpha, - vVect, - anorm, - rhobar, - phibar, - wVect, - xSolution, - standardError, - dnorm, - sn2, - cs2, - z, - xnorm1, - res2, - nstop, - comlsqr); - - if (myid == 0) - { - CPRtimeEnd = time(NULL) - CPRtimeStart; - totTime += CPRtimeEnd; - partialTime += CPRtimeEnd; - printf("CPR: itn=%d writing CPR seconds %ld. Global Seconds %ld \n", itn, CPRtimeEnd, totTime); + }//if(!noCPR) + }//if(myid==0) + MPI_Barrier(MPI_COMM_WORLD); + MPI_Bcast( &writeCPR, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &itnlim, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &comlsqr.itnLimit, 1, MPI_INT, 0, MPI_COMM_WORLD); + if(myid==0) + { + printf("lsqr: Iteration number %d. Iteration seconds %ld. OmpSec %d Global Seconds %ld \n",itn,endCycleTime,ompSec,totTime); + if(writeCPR==1) printf("... writing CPR files\n"); + if(writeCPR==2) printf("... writing CPR files and return\n"); + } + + + if(writeCPR>0) + { + if(myid==0) + CPRtimeStart=time(NULL); + writeCheckPoint(itn, + knownTerms, + beta, + alpha, + vVect, + anorm, + rhobar, + phibar, + wVect, + xSolution, + standardError, + dnorm, + sn2, + cs2, + z, + xnorm1, + res2, + nstop, + comlsqr); + + if(myid==0){ + CPRtimeEnd=time(NULL)-CPRtimeStart; + totTime+=CPRtimeEnd; + partialTime+=CPRtimeEnd; + printf("CPR: itn=%d writing CPR seconds %ld. Global Seconds %ld \n",itn,CPRtimeEnd, totTime); + } + + if(writeCPR==2) + { + *istop_out = 1000; + *itn_out = itn; + if(nAstroPSolved) free(vAuxVect); + return; + } + if(myid==0) startCycleTime=time(NULL); } - if (writeCPR == 2) + itn = itn + 1; + comlsqr.itn=itn; + CPRCount++; +// ------------------------------------------------------------------ +// Perform the next step of the bidiagonalization to obtain the +// next beta, u, alpha, v. These satisfy the relations +// beta*u = A*v - alpha*u, +// alpha*v = A(transpose)*u - beta*v. +// ------------------------------------------------------------------ + cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (- alpha), knownTerms, 1 ); + kAuxcopy=(double *) calloc(nEqExtConstr+nEqBarConstr+nOfInstrConstr, sizeof(double)); + for (int kk=0;kk ZERO) - { - cblas_dscal(mapNoss[myid] + nEqExtConstr + nEqBarConstr + nOfInstrConstr, (ONE / beta), knownTerms, 1); - cblas_dscal(nunkSplit, (-beta), vVect, 1); -//only processor 0 accumulate and distribute the v array - // FV EDIT - //#pragma omp for - for (i = 0; i < localAstroMax; i++) - { - vAuxVect[i] = vVect[i]; - vVect[i] = 0; - } - if (myid != 0) - { - // FV EDIT - //#pragma omp for - for (i = localAstroMax; i < nunkSplit; i++) - vVect[i] = 0; - } + temp = d2norm( alpha, beta ); + temp = d2norm( temp , damp ); + anorm = d2norm( anorm, temp ); - aprod(2, m, n, vVect, knownTerms, systemMatrix, matrixIndex, instrCol, instrConstrIlung, comlsqr, &ompSec); - MPI_Barrier(MPI_COMM_WORLD); + if (beta > ZERO) { + cblas_dscal ( mapNoss[myid]+nEqExtConstr+nEqBarConstr+nOfInstrConstr, (ONE / beta), knownTerms, 1 ); + cblas_dscal ( nunkSplit, (- beta), vVect, 1 ); +//only processor 0 accumulate and distribute the v array +#pragma omp for + for(i=0;i ZERO) { + cblas_dscal ( nunkSplit, (ONE / alpha), vVect, 1 ); + } } - double alphaLoc = 0; - if (nAstroPSolved) - alphaLoc = cblas_dnrm2(nAstroElements * comlsqr.nAstroPSolved, vVect, 1); - double alphaLoc2 = alphaLoc * alphaLoc; - if (myid == 0) - { - double alphaOther = cblas_dnrm2(comlsqr.nunkSplit - localAstroMax, &vVect[localAstroMax], 1); - double alphaOther2 = alphaOther * alphaOther; - alphaLoc2 = alphaLoc2 + alphaOther2; +// ------------------------------------------------------------------ +// Use a plane rotation to eliminate the damping parameter. +// This alters the diagonal (rhobar) of the lower-bidiagonal matrix. +// ------------------------------------------------------------------ + rhbar1 = rhobar; + if ( damped ) { + rhbar1 = d2norm( rhobar, damp ); + cs1 = rhobar / rhbar1; + sn1 = damp / rhbar1; + psi = sn1 * phibar; + phibar = cs1 * phibar; } - double alphaGlob2 = 0; - MPI_Allreduce(&alphaLoc2, &alphaGlob2, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - alpha = sqrt(alphaGlob2); - if (alpha > ZERO) - { - cblas_dscal(nunkSplit, (ONE / alpha), vVect, 1); +// ------------------------------------------------------------------ +// Use a plane rotation to eliminate the subdiagonal element (beta) +// of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. +// ------------------------------------------------------------------ + rho = d2norm( rhbar1, beta ); + cs = rhbar1 / rho; + sn = beta / rho; + theta = sn * alpha; + rhobar = - cs * alpha; + phi = cs * phibar; + phibar = sn * phibar; + tau = sn * phi; + +// ------------------------------------------------------------------ +// Update x, w and (perhaps) the standard error estimates. +// ------------------------------------------------------------------ + t1 = phi / rho; + t2 = - theta / rho; + t3 = ONE / rho; + dknorm = ZERO; + + for (i = 0; i < nAstroElements*comlsqr.nAstroPSolved; i++) { + t = wVect[i]; + t = (t3*t)*(t3*t); + dknorm = t + dknorm; } - } - - // ------------------------------------------------------------------ - // Use a plane rotation to eliminate the damping parameter. - // This alters the diagonal (rhobar) of the lower-bidiagonal matrix. - // ------------------------------------------------------------------ - rhbar1 = rhobar; - if (damped) - { - rhbar1 = d2norm(rhobar, damp); - cs1 = rhobar / rhbar1; - sn1 = damp / rhbar1; - psi = sn1 * phibar; - phibar = cs1 * phibar; - } - - // ------------------------------------------------------------------ - // Use a plane rotation to eliminate the subdiagonal element (beta) - // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. - // ------------------------------------------------------------------ - rho = d2norm(rhbar1, beta); - cs = rhbar1 / rho; - sn = beta / rho; - theta = sn * alpha; - rhobar = -cs * alpha; - phi = cs * phibar; - phibar = sn * phibar; - tau = sn * phi; - - // ------------------------------------------------------------------ - // Update x, w and (perhaps) the standard error estimates. - // ------------------------------------------------------------------ - t1 = phi / rho; - t2 = -theta / rho; - t3 = ONE / rho; - dknorm = ZERO; - - for (i = 0; i < nAstroElements * comlsqr.nAstroPSolved; i++) - { - t = wVect[i]; - t = (t3 * t) * (t3 * t); - dknorm = t + dknorm; - } - dknormSum = 0; - MPI_Allreduce(&dknorm, &dknormSum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); - dknorm = dknormSum; - - if (wantse) - { - // FV EDIT - //#pragma omp for - for (i = 0; i < localAstro; i++) - { - t = wVect[i]; - xSolution[i] = t1 * t + xSolution[i]; - wVect[i] = t2 * t + vVect[i]; - t = (t3 * t) * (t3 * t); - standardError[i] = t + standardError[i]; + dknormSum=0; + MPI_Allreduce(&dknorm,&dknormSum,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); + dknorm=dknormSum; + + if ( wantse ) { +#pragma omp for + for (i = 0; i < localAstro; i++) { + t = wVect[i]; + xSolution[i] = t1*t + xSolution[i]; + wVect[i] = t2*t + vVect[i]; + t = (t3*t)*(t3*t); + standardError[i] = t + standardError[i]; + } +#pragma omp for + for (i = localAstroMax; i < localAstroMax+other; i++) { + t = wVect[i]; + xSolution[i] = t1*t + xSolution[i]; + wVect[i] = t2*t + vVect[i]; + t = (t3*t)*(t3*t); + standardError[i] = t + standardError[i]; + } + for (i = localAstroMax; i < localAstroMax+other; i++) { + t = wVect[i]; + t = (t3*t)*(t3*t); + dknorm = t + dknorm; + } } - // FV EDIT - //#pragma omp for - for (i = localAstroMax; i < localAstroMax + other; i++) - { - t = wVect[i]; - xSolution[i] = t1 * t + xSolution[i]; - wVect[i] = t2 * t + vVect[i]; - t = (t3 * t) * (t3 * t); - standardError[i] = t + standardError[i]; + else { +#pragma omp for + for (i = 0; i < localAstro; i++) { + t = wVect[i]; + xSolution[i] = t1*t + xSolution[i]; + wVect[i] = t2*t + vVect[i]; + } +#pragma omp for + for (i = localAstroMax; i < localAstroMax+other; i++) { + t = wVect[i]; + xSolution[i] = t1*t + xSolution[i]; + wVect[i] = t2*t + vVect[i]; + } + for (i = localAstroMax; i < localAstroMax+other; i++) { + t = wVect[i]; + dknorm = (t3*t)*(t3*t) + dknorm; + } } - for (i = localAstroMax; i < localAstroMax + other; i++) - { - t = wVect[i]; - t = (t3 * t) * (t3 * t); - dknorm = t + dknorm; + +// ------------------------------------------------------------------ +// Monitor the norm of d_k, the update to x. +// dknorm = norm( d_k ) +// dnorm = norm( D_k ), where D_k = (d_1, d_2, ..., d_k ) +// dxk = norm( phi_k d_k ), where new x = x_k + phi_k d_k. +// ------------------------------------------------------------------ + dknorm = sqrt( dknorm ); + dnorm = d2norm( dnorm, dknorm ); + dxk = fabs( phi * dknorm ); + if (dxmax < dxk ) { + dxmax = dxk; + maxdx = itn; } - } - else - { - // FV EDIT - //#pragma omp for - for (i = 0; i < localAstro; i++) - { - t = wVect[i]; - xSolution[i] = t1 * t + xSolution[i]; - wVect[i] = t2 * t + vVect[i]; + +// ------------------------------------------------------------------ +// Use a plane rotation on the right to eliminate the +// super-diagonal element (theta) of the upper-bidiagonal matrix. +// Then use the result to estimate norm(x). +// ------------------------------------------------------------------ + delta = sn2 * rho; + gambar = - cs2 * rho; + rhs = phi - delta * z; + zbar = rhs / gambar; + xnorm = d2norm( xnorm1, zbar ); + gamma = d2norm( gambar, theta ); + cs2 = gambar / gamma; + sn2 = theta / gamma; + z = rhs / gamma; + xnorm1 = d2norm( xnorm1, z ); + +// ------------------------------------------------------------------ +// Test for convergence. +// First, estimate the norm and condition of the matrix Abar, +// and the norms of rbar and Abar(transpose)*rbar. +// ------------------------------------------------------------------ + acond = anorm * dnorm; + res2 = d2norm( res2 , psi ); + rnorm = d2norm( res2 , phibar ); + arnorm = alpha * fabs( tau ); + +// Now use these norms to estimate certain other quantities, +// some of which will be small near a solution. + + alfopt = sqrt( rnorm / (dnorm * xnorm) ); + test1 = rnorm / bnorm; + test2 = ZERO; + if (rnorm > ZERO) test2 = arnorm / (anorm * rnorm); +// if (arnorm0 > ZERO) test2 = arnorm / arnorm0; //(Michael Friedlander's modification) + test3 = ONE / acond; + t1 = test1 / (ONE + anorm * xnorm / bnorm); + rtol = btol + atol * anorm * xnorm / bnorm; + +// The following tests guard against extremely small values of +// atol, btol or ctol. (The user may have set any or all of +// the parameters atol, btol, conlim to zero.) +// The effect is equivalent to the normal tests using +// atol = relpr, btol = relpr, conlim = 1/relpr. + + t3 = ONE + test3; + t2 = ONE + test2; + t1 = ONE + t1; +// printf("LSQR TP8 PE=%d itn=%d test1=%18.16lf, rnorm=%lf bnorm=%lf anorm=%lf xnorm=%lf rtol=%lf t1=%18.16lf\n",myid,itn,test1,rnorm,bnorm,anorm,xnorm,rtol,t1); + + if (itn >= itnlim) istop = 5; + if (t3 <= ONE ) istop = 4; + if (t2 <= ONE ) istop = 2; + if (t1 <= ONE ) istop = 1; +// if (t1 <= ONE ) printf("PE=%d t1=%lf\n",myid,t1); + +// Allow for tolerances set by the user. + + if (test3 <= ctol) istop = 4; + if (test2 <= atol) istop = 2; + if (test1 <= rtol) istop = 1; //(Michael Friedlander had this commented out) +// if (test1 <= rtol) printf("PE=%d test1=%lf\n",myid,test1); + +// ------------------------------------------------------------------ +// See if it is time to print something. +// ------------------------------------------------------------------ +// if (nout == NULL ) goto goto_600; // Delete for buffer flush modification??? TBV + if (n <= 40 ) goto goto_400; + if (itn <= 10 ) goto goto_400; + if (itn >= itnlim-10) goto goto_400; + if (itn % 10 == 0 ) goto goto_400; + if (test3 <= 2.0*ctol) goto goto_400; + if (test2 <= 10.0*atol) goto goto_400; + if (test1 <= 10.0*rtol) goto goto_400; + if (istop != 0 ) goto goto_400; + goto goto_600; + +// Print a line for this iteration. +// "extra" is for experimental purposes. + + goto_400: + if(!comlsqr.Test) + if(myid==0) { + nout=fopen(lsqrOut,"a"); + if ( nout != NULL) { + if ( extra ) { + fprintf(nout, fmt_1500_extra, + itn, xSolution[0], rnorm, test1, test2, anorm, + acond, phi, dknorm, dxk, alfopt); + } else { + fprintf(nout, fmt_1500, + itn, xSolution[0], rnorm, test1, test2, anorm, acond); + } + if (itn % 10 == 0) fprintf(nout, fmt_1600); + } else { + printf("Error while opening %s (3)\n",lsqrOut); + MPI_Abort(MPI_COMM_WORLD,1); + exit(EXIT_FAILURE); + } + fclose(nout); } - // FV EDIT - //#pragma omp for - for (i = localAstroMax; i < localAstroMax + other; i++) - { - t = wVect[i]; - xSolution[i] = t1 * t + xSolution[i]; - wVect[i] = t2 * t + vVect[i]; +// ------------------------------------------------------------------ +// Stop if appropriate. +// The convergence criteria are required to be met on nconv +// consecutive iterations, where nconv is set below. +// Suggested value: nconv = 1, 2 or 3. +// ------------------------------------------------------------------ + goto_600: + if (istop == 0) { + nstop = 0; } - for (i = localAstroMax; i < localAstroMax + other; i++) - { - t = wVect[i]; - dknorm = (t3 * t) * (t3 * t) + dknorm; + else { + nconv = 1; + nstop = nstop + 1; + if (nstop < nconv && itn < itnlim) istop = 0; } - } - - // ------------------------------------------------------------------ - // Monitor the norm of d_k, the update to x. - // dknorm = norm( d_k ) - // dnorm = norm( D_k ), where D_k = (d_1, d_2, ..., d_k ) - // dxk = norm( phi_k d_k ), where new x = x_k + phi_k d_k. - // ------------------------------------------------------------------ - dknorm = sqrt(dknorm); - dnorm = d2norm(dnorm, dknorm); - dxk = fabs(phi * dknorm); - if (dxmax < dxk) - { - dxmax = dxk; - maxdx = itn; - } - - // ------------------------------------------------------------------ - // Use a plane rotation on the right to eliminate the - // super-diagonal element (theta) of the upper-bidiagonal matrix. - // Then use the result to estimate norm(x). - // ------------------------------------------------------------------ - delta = sn2 * rho; - gambar = -cs2 * rho; - rhs = phi - delta * z; - zbar = rhs / gambar; - xnorm = d2norm(xnorm1, zbar); - gamma = d2norm(gambar, theta); - cs2 = gambar / gamma; - sn2 = theta / gamma; - z = rhs / gamma; - xnorm1 = d2norm(xnorm1, z); - - // ------------------------------------------------------------------ - // Test for convergence. - // First, estimate the norm and condition of the matrix Abar, - // and the norms of rbar and Abar(transpose)*rbar. - // ------------------------------------------------------------------ - acond = anorm * dnorm; - res2 = d2norm(res2, psi); - rnorm = d2norm(res2, phibar); - arnorm = alpha * fabs(tau); - - // Now use these norms to estimate certain other quantities, - // some of which will be small near a solution. - - alfopt = sqrt(rnorm / (dnorm * xnorm)); - test1 = rnorm / bnorm; - test2 = ZERO; - if (rnorm > ZERO) - test2 = arnorm / (anorm * rnorm); - // if (arnorm0 > ZERO) test2 = arnorm / arnorm0; //(Michael Friedlander's modification) - test3 = ONE / acond; - t1 = test1 / (ONE + anorm * xnorm / bnorm); - rtol = btol + atol * anorm * xnorm / bnorm; - - // The following tests guard against extremely small values of - // atol, btol or ctol. (The user may have set any or all of - // the parameters atol, btol, conlim to zero.) - // The effect is equivalent to the normal tests using - // atol = relpr, btol = relpr, conlim = 1/relpr. - - t3 = ONE + test3; - t2 = ONE + test2; - t1 = ONE + t1; - // printf("LSQR TP8 PE=%d itn=%d test1=%18.16lf, rnorm=%lf bnorm=%lf anorm=%lf xnorm=%lf rtol=%lf t1=%18.16lf\n",myid,itn,test1,rnorm,bnorm,anorm,xnorm,rtol,t1); - - if (itn >= itnlim) - istop = 5; - if (t3 <= ONE) - istop = 4; - if (t2 <= ONE) - istop = 2; - if (t1 <= ONE) - istop = 1; - // if (t1 <= ONE ) printf("PE=%d t1=%lf\n",myid,t1); - - // Allow for tolerances set by the user. - - if (test3 <= ctol) - istop = 4; - if (test2 <= atol) - istop = 2; - if (test1 <= rtol) - istop = 1; //(Michael Friedlander had this commented out) - // if (test1 <= rtol) printf("PE=%d test1=%lf\n",myid,test1); - - // ------------------------------------------------------------------ - // See if it is time to print something. - // ------------------------------------------------------------------ - // if (nout == NULL ) goto goto_600; // Delete for buffer flush modification??? TBV - if (n <= 40) - goto goto_400; - if (itn <= 10) - goto goto_400; - if (itn >= itnlim - 10) - goto goto_400; - if (itn % 10 == 0) - goto goto_400; - if (test3 <= 2.0 * ctol) - goto goto_400; - if (test2 <= 10.0 * atol) - goto goto_400; - if (test1 <= 10.0 * rtol) - goto goto_400; - if (istop != 0) - goto goto_400; - goto goto_600; - - // Print a line for this iteration. - // "extra" is for experimental purposes. - -goto_400: - if (!comlsqr.Test) - if (myid == 0) + + if(comlsqr.Test) { - nout = fopen(lsqrOut, "a"); - if (nout != NULL) - { - if (extra) - { - fprintf(nout, fmt_1500_extra, - itn, xSolution[0], rnorm, test1, test2, anorm, - acond, phi, dknorm, dxk, alfopt); - } - else - { - fprintf(nout, fmt_1500, - itn, xSolution[0], rnorm, test1, test2, anorm, acond); - } - if (itn % 10 == 0) - fprintf(nout, fmt_1600); - } - else - { - printf("Error while opening %s (3)\n", lsqrOut); - MPI_Abort(MPI_COMM_WORLD, 1); - exit(EXIT_FAILURE); - } - fclose(nout); - } - // ------------------------------------------------------------------ - // Stop if appropriate. - // The convergence criteria are required to be met on nconv - // consecutive iterations, where nconv is set below. - // Suggested value: nconv = 1, 2 or 3. - // ------------------------------------------------------------------ -goto_600: - if (istop == 0) - { - nstop = 0; + if(itn n) - t = m - n; - if (damped) - t = m; - t = rnorm / sqrt(t); - - for (i = 0; i < nunkSplit; i++) - standardError[i] = t * sqrt(standardError[i]); -} + if ( wantse ) { + t = ONE; + if (m > n) t = m - n; + if ( damped ) t = m; + t = rnorm / sqrt( t ); + + for (i = 0; i < nunkSplit; i++) + standardError[i] = t * sqrt( standardError[i] ); + + } // Decide if istop = 2 or 3. // Print the stopping condition. -goto_800 : if (damped && istop == 2) istop = 3; -if (myid == 0) -{ - nout = fopen(lsqrOut, "a"); - if (nout != NULL) - { - fprintf(nout, fmt_2000, - exitLsqr, istop, itn, - exitLsqr, anorm, acond, - exitLsqr, bnorm, xnorm, - exitLsqr, rnorm, arnorm); - fprintf(nout, fmt_2100, - exitLsqr, dxmax, maxdx, - exitLsqr, dxmax / (xnorm + 1.0e-20)); - fprintf(nout, fmt_3000, - exitLsqr, msg[istop]); - } - else - { - printf("Error while opening %s (4)\n", lsqrOut); - MPI_Abort(MPI_COMM_WORLD, 1); - exit(EXIT_FAILURE); - } - fclose(nout); -} + goto_800: + if (damped && istop == 2) istop = 3; + if (myid == 0) { + nout=fopen(lsqrOut,"a"); + if (nout != NULL) { + fprintf(nout, fmt_2000, + exitLsqr, istop, itn, + exitLsqr, anorm, acond, + exitLsqr, bnorm, xnorm, + exitLsqr, rnorm, arnorm); + fprintf(nout, fmt_2100, + exitLsqr, dxmax, maxdx, + exitLsqr, dxmax/(xnorm + 1.0e-20)); + fprintf(nout, fmt_3000, + exitLsqr, msg[istop]); + } else { + printf("Error while opening %s (4)\n",lsqrOut); + MPI_Abort(MPI_COMM_WORLD,1); + exit(EXIT_FAILURE); + } + fclose(nout); + } // Assign output variables from local copies. -*istop_out = istop; -*itn_out = itn; -*anorm_out = anorm; -*acond_out = acond; -*rnorm_out = rnorm; -*arnorm_out = test2; -*xnorm_out = xnorm; - -if (nAstroPSolved) - free(vAuxVect); -return; + *istop_out = istop; + *itn_out = itn; + *anorm_out = anorm; + *acond_out = acond; + *rnorm_out = rnorm; + *arnorm_out = test2; + *xnorm_out = xnorm; + + if(nAstroPSolved) free(vAuxVect); + return; } + + diff --git a/memRequest.c b/memRequest.c old mode 100644 new mode 100755 diff --git a/run_ppl.cmd b/run_ppl.cmd new file mode 100755 index 0000000000000000000000000000000000000000..636fd0f8866ad4cc9a9601660ed8ef548bdf13e4 --- /dev/null +++ b/run_ppl.cmd @@ -0,0 +1,28 @@ +#!/bin/bash +#SBATCH --job-name=AVUGSR_OMP +#SBATCH --time=0:30:00 +#SBATCH --out=out.log +#SBATCH --err=err.log +###SBATCH --nodes=2 --ntasks-per-node=2 --ntasks-per-socket=2 +#SBATCH --nodes=1 +#SBATCH --ntasks=2 +#SBATCH --cpus-per-task=4 +#SBATCH --mem=100000 +#SBATCH --partition=merlin +####SBATCH --account=INAF_GAIA14_3 +# +# +# +###module load profile/astro +##module load autoload intelmpi cfitsio +export OMP_NUM_THREADS=4 +export KMP_AFFINITY=compact +# +cd /home/ubecciani/AVU-GSR/OMP/src/run + +srun ../GaiaGsrParSim -auto + +exit 0 + +# +# diff --git a/solvergaiaSim.c b/solvergaiaSim.c index 30b8f0aa6fb9ec457d3b6fa7da402c193f6bdacf..d001f9ad8618aa8f4c38c26d998eaa5d63086274 100644 --- a/solvergaiaSim.c +++ b/solvergaiaSim.c @@ -16,6 +16,8 @@ - version 5.0 - May 2013 from Ugo Becciani and Alberto Vecchiato */ + + //#include #include #include @@ -32,6 +34,7 @@ #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); @@ -39,570 +42,505 @@ 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 main(int argc, char **argv) { + int debugMode,wrsol; int i; - - int idtest = 0, precond = 1; + + int idtest=0, precond=1; double srIDtest, pert; - - int inputDirOpt = 0, outputDirOpt = 0, inputDirLen; - char inputDir[1024] = "", outputDir[1024] = "", wrfileDir[1024] = ""; + + int inputDirOpt=0, outputDirOpt=0, inputDirLen; + char inputDir[1024]="", outputDir[1024]="",wrfileDir[1024]=""; char wpath[1024]; - size_t sizePath = 1020; - + 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 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() - + 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 + 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 - + double aTol; //read by command line, overrides atol if >=0 + // LSQR output parameters - int istop; // LSQR stopping condition - int itn; // LSQR iteration number + 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; + 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 + + + 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 + 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 - + 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 + 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 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; // + 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; + 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] + 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. + 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; - + 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 *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 *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; + int nproc, myid, namelen; char processor_name[MPI_MAX_PROCESSOR_NAME]; - double startTime, endTime; - + double startTime,endTime; + // Distributed array maps - long int *mapNcoeff, *mapNoss; - long int mapNcoeffBefore, mapNossBefore; - long int mapNcoeffAfter, mapNossAfter; + long int * mapNcoeff, * mapNoss; + long int mapNcoeffBefore, mapNossBefore; + long int mapNcoeffAfter, mapNossAfter; // struct to communicate with lsqr - int timeCPR, timeLimit, itnCPR, itnCPRstop = 1000000; + 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; + 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; + 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 dirent **namelistMI; struct nullSpace nullSpaceCk; double nullSpaceAttfact; - int autoRun = 0; - float memGlobal = 0; + int autoRun=0; + float memGlobal=0; float memGB; - int nthreads=0; - int ntasks=0; - /////////////// - for (int i = 1; i < argc - 1; i++) + for (int i=1;i 0) - timeCPR = testTime; + nfileProc=atoi(argv[i+1]); + if(nfileProc <=0) nfileProc=3; } - if (strcmp(argv[i], "-timelimit") == 0) + if(strcmp (argv[i],"-timeCPR") == 0) { - testTime = atoi(argv[i + 1]); - if (testTime > 0) - timeLimit = testTime; + testTime=atoi(argv[i+1]); + if(testTime >0) timeCPR=testTime; } - - if (strcmp(argv[i], "-itnCPR") == 0) + if(strcmp (argv[i],"-timelimit") == 0) { - testTime = atoi(argv[i + 1]); - if (testTime > 0) - itnCPR = testTime; + testTime=atoi(argv[i+1]); + if(testTime >0) timeLimit=testTime; } - if (strcmp(argv[i], "-noCPR") == 0) + + if(strcmp (argv[i],"-itnCPR") == 0) { - noCPR = 1; + testTime=atoi(argv[i+1]); + if(testTime >0) itnCPR=testTime; } - if (strcmp(argv[i], "-itnlimit") == 0) + if(strcmp (argv[i],"-noCPR") == 0) { - testTime = atoi(argv[i + 1]); - if (testTime > 0) - itnLimit = testTime; + noCPR=1; } - if(strcmp (argv[i],"-tasks") == 0) + if(strcmp (argv[i],"-itnlimit") == 0) { - ntasks=atoi(argv[i+1]); - if(ntasks <=0) ntasks=0; + testTime=atoi(argv[i+1]); + if(testTime >0) itnLimit=testTime; } - if (strcmp(argv[i], "-atol") == 0) + if(strcmp (argv[i],"-atol") == 0) { - double dummy = atof(argv[i + 1]); - if (dummy >= 0) - aTol = dummy; + double dummy=atof(argv[i+1]); + if(dummy >=0) aTol=dummy; } - if (strcmp(argv[i], "-zeroAtt") == 0) + if(strcmp (argv[i],"-zeroAtt") == 0) { - zeroAtt = 1; + zeroAtt=1; } - if (strcmp(argv[i], "-zeroInstr") == 0) + if(strcmp (argv[i],"-zeroInstr") == 0) { - zeroInstr = 1; + zeroInstr=1; } - if (strcmp(argv[i], "-zeroGlob") == 0) + if(strcmp (argv[i],"-zeroGlob") == 0) { - zeroGlob = 1; + zeroGlob=1; } - if (strcmp(argv[i], "-wgic") == 0) + if(strcmp (argv[i],"-wgic") == 0) { - wgInstrCoeff = atoi(argv[i + 1]); - if (wgInstrCoeff <= 0) - wgInstrCoeff = 1; + wgInstrCoeff=atoi(argv[i+1]); + if(wgInstrCoeff<=0) wgInstrCoeff=1; } // inputDir e outputDir - - if (strcmp(argv[i], "-inputDir") == 0) + + if(strcmp (argv[i],"-inputDir") == 0) { - sprintf(inputDir, "%s", argv[i + 1]); - inputDirOpt = 1; + sprintf(inputDir, "%s", argv[i+1]); + inputDirOpt=1; } - - if (strcmp(argv[i], "-outputDir") == 0) + + if(strcmp (argv[i],"-outputDir") == 0) { - sprintf(outputDir, "%s", argv[i + 1]); - outputDirOpt = 1; + 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],"-debug") == 0) + debugMode=1; + + if(strcmp (argv[i],"-wrsol") == 0){ + wrsol=1; } - - if (strcmp(argv[i], "-noiter") == 0) - { - noIter = 1; + + if(strcmp (argv[i],"-noiter") == 0){ + noIter=1; } - - if (strcmp(argv[i], "-wrfilebin") == 0) + + if(strcmp (argv[i],"-wrfilebin") == 0) { - sprintf(wrfileDir, "%s", argv[i + 1]); - wrFilebin = 1; + sprintf(wrfileDir, "%s", argv[i+1]); + wrFilebin=1; } - if (strcmp(argv[i], "-extConstr") == 0) + 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); + 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) + 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); + 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) - { + if (extConstraint) noConstr=2; + if (barConstraint) noConstr=2; + + + if(myid==0) + { + printf("Execution of solvergaia Simulator version %s on %d mpi-tasks\n solvergaiaSim ",VER,nproc); + if(idtest)printf ("-IDtest %le ", srIDtest); + if(precond)printf ("-Precond on"); + else printf ("-Precond off"); + if(nfileProc!=3) printf ("-numFilexproc %d ",nfileProc); + if(wrFilebin) printf(" -wrfilebin %s ",wrfileDir); + if(!withFile) printf(" -noFile "); + if(itnLimit>0) printf(" -itnlimit %d ",itnLimit); + if(noCPR>0) printf(" -noCPR "); + if(itnCPR!= DEFAULT_ITNCPR) printf(" -itnCPR %d ", itnCPR); + if(zeroAtt) printf(" -zeroAtt "); + if(noConstr==1) printf(" -noConstr "); + if(zeroInstr) printf(" -zeroInstr "); + if(zeroGlob) printf(" -zeroGlob "); + if(inputDirOpt) printf(" -inputDir %s ", inputDir); + if(outputDirOpt) printf(" -outputDir %s ", outputDir); + if(wrsol) printf(" -wrsol "); + if(noIter) printf(" -noiter "); + if(debugMode) printf(" -debug "); + if(extConstraint)printf("-extConstr %le ",extConstrW); + if(barConstraint)printf("-barConstr %le ",barConstrW); + printf("-wgic %d", wgInstrCoeff); + if(!autoRun) printf(" %s\n",argv[argc-1]); + if(extConstraint && barConstraint) { printf("Error: baricentric anld null space constraints are mutually exclusive. Aborting\n"); - MPI_Abort(MPI_COMM_WORLD, 1); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } - if (wrFilebin && strcmp(inputDir, wrfileDir) == 0) - { + if(wrFilebin && strcmp(inputDir, wrfileDir)==0){ printf("inputDir and wrfilebinDir are the same. Execution Aborted.\n"); - MPI_Abort(MPI_COMM_WORLD, 1); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); + } - printf ("\n"); } ///////////////// - getcwd(wpath, sizePath); - if (!inputDirOpt) + getcwd(wpath,sizePath); + if(!inputDirOpt) strcpy(inputDir, wpath); - printf("Process %d running on %s\n", myid, processor_name); -////#ifdef OMP -///#pragma omp task + printf("Process %d running on %s\n",myid,processor_name); +#ifdef OMP +#pragma omp parallel { - // 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); + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + printf("PE=%d on processor %s total num of threads =%d ID thread =%d \n",myid, processor_name,nthreads,tid); + } -////#pragma omp taskwait -////#endif - nullSpaceAttfact = 1.0 * attExtConstrFact * extConstrW; - comlsqr.extConstrW = extConstrW; - comlsqr.nullSpaceAttfact = nullSpaceAttfact; - comlsqr.barConstrW = barConstrW; - - if (inputDirOpt) +#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)) - { + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } chdir(wpath); - if (myid == 0) - printf("success.\n"); + if(myid==0) printf("success.\n"); } - inputDirLen = strlen(inputDir); - sprintf(filenameSolProps, "%s", argv[argc - 1]); - sprintf(filenameSolPropsFinal, "GsrFinal_%s", argv[argc - 1]); - - if (outputDirOpt) + 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)) - { + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } - if (myid == 0) - printf("success.\n"); + 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; + 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) - { + if(myid==0 && wrFilebin){ chdir(wpath); chdir(wrfileDir); - fpSolPropsFileBin = fopen(filenameSolProps, "w"); + 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; + 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]; + 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; + 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); + 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); @@ -619,7 +557,7 @@ int main(int argc, char **argv) printf("nConstrMuLong= %5ld\n", nConstrMuLong); printf("nDegFreedomAtt= %7ld\n", nDegFreedomAtt); printf("nAttAxes= %7hd\n", nAttAxes); - printf("nFovs= %7d ", instrConst[0] + 1); + 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]); @@ -634,750 +572,616 @@ int main(int argc, char **argv) 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); + 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); } - else + 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) { - 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); + constrLatId = (long *) calloc(nConstrLat, sizeof(long)); + for(i=0;i 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); + 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 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); + 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 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); + 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 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 - } + 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); - endTime = MPI_Wtime(); - if ((nDegFreedomAtt == 0 && nAttAxes > 0) || (nDegFreedomAtt > 0 && nAttAxes == 0)) - { - if (myid == 0) + 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("inconsistent values for nDegFreedomAtt=%ld and nAttaxes=%d\n", nDegFreedomAtt, nAttAxes); - MPI_Abort(MPI_COMM_WORLD, 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 } - - 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)); + + 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); + } } - 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; + + if(myid==0) printf("Time to read initial parameters =%f sec.\n",endTime-startTime); + // Start section for parameter broadcast + MPI_Bcast( &atol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( &btol, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( &conlim, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( &damp, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( &nAstroP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nAstroPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nInstrP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nInstrPSolved, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &lsInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &ssInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nuInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &maInstrFlag, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nElemIC, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nOfInstrConstr, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast( &nGlobP, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &itnlim, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nStar, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nDegFreedomAtt, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nAttAxes, 1, MPI_SHORT, 0, MPI_COMM_WORLD); + MPI_Bcast( &instrSetUp, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( instrConst, DEFAULT_NINSTRINDEXES , MPI_INT, 0, MPI_COMM_WORLD); // errore + MPI_Bcast( &nobs, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nConstrLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nConstrLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nConstrMuLat, 1, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( &nConstrMuLong, 1, MPI_LONG, 0, MPI_COMM_WORLD); + if(myid !=0) + { + constrLatId=(long *) calloc(nConstrLat , sizeof(long)); + constrLatW=(double *) calloc(nConstrLat , sizeof(double)); + constrLongId=(long *) calloc(nConstrLong , sizeof(long)); + constrLongW=(double *) calloc(nConstrLong , sizeof(double)); + constrMuLatId=(long *) calloc(nConstrMuLat , sizeof(long)); + constrMuLatW=(double *) calloc(nConstrMuLat , sizeof(double)); + constrMuLongId=(long *) calloc(nConstrMuLong , sizeof(long)); + constrMuLongW=(double *) calloc(nConstrMuLong , sizeof(double)); + } + MPI_Bcast( constrLatId, nConstrLat, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( constrLatW, nConstrLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( constrLongId, nConstrLong, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( constrLongW, nConstrLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( constrMuLatId, nConstrMuLat, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( constrMuLatW, nConstrMuLat, MPI_DOUBLE, 0, MPI_COMM_WORLD); + MPI_Bcast( constrMuLongId, nConstrMuLong, MPI_LONG, 0, MPI_COMM_WORLD); + MPI_Bcast( constrMuLongW, nConstrMuLong, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + if(nInstrPSolved==0) + zeroInstr=1; /////////// Multiplicity of matrixIndex - + int multMI; - - if (nAstroPSolved) - multMI = 2; - else - { - multMI = 1; - extConstraint = 0; - barConstraint = 0; + + 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]; + 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; + 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; + + 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 + 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; + 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 @@ -1388,1012 +1192,909 @@ int main(int argc, char **argv) // 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); - } + + 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) + + nObsxStar=nobs/nStar; + nobsOri=nobs; + if(noConstr) { - nConstr = 0; - nConstrLong = 0; - nConstrLat = 0; - nConstrMuLong = 0; - nConstrMuLat = 0; - } - else + nConstr=0; + nConstrLong=0; + nConstrLat=0; + nConstrMuLong=0; + nConstrMuLat=0; + } else { - for (int q = 0; q < nConstrLat; q++) - if (constrLatId[q] >= nStar) + for(int q=0;q=nStar) { - printf("PE=%d Error invalit Lat Constraint %ld\n", myid, constrLatId[q]); - MPI_Abort(MPI_COMM_WORLD, 1); + 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) + + for(int q=0;q=nStar) { - printf("PE=%d Error invalit Long Constraint %ld\n", myid, constrLongId[q]); - MPI_Abort(MPI_COMM_WORLD, 1); + 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) + for(int q=0;q=nStar) { - printf("PE=%d Error invalit MuLat Constraint %ld\n", myid, constrMuLatId[q]); - MPI_Abort(MPI_COMM_WORLD, 1); + 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) + + for(int q=0;q=nStar) { - printf("PE=%d Error invalit MuLat Constraint %ld\n", myid, constrMuLongId[q]); - MPI_Abort(MPI_COMM_WORLD, 1); + 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); + + + 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); + + 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=(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= i + 1) - mapNoss[i]++; - mapNcoeff[i] = mapNoss[i] * nparam; - if (i < myid) + mapNoss[i]=(nobs)/nproc; + if(nobs % nproc >=i+1) mapNoss[i]++; + mapNcoeff[i]=mapNoss[i]*nparam; + if(i myid) + if(i>myid) { - mapNossAfter += mapNoss[i]; - mapNcoeffAfter += mapNcoeff[i]; + 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) - { + if(extConstraint){ + int * sumNObsxStar; + sumNObsxStar=(int *) calloc(nStar, sizeof(int)); + int irest=nobs % nStar; + for(int i=0;i 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; + + } + long counterObsxStar=0; + for(int i=0;imapNossBefore && 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 + 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 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; + + } + int counterObsxStar=0; + for(int i=0;imapNossBefore && 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 + 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; + + + comlsqr.nStar=nStar; + comlsqr.nAstroP=nAstroP; + comlsqr.nAstroPSolved=nAstroPSolved; + comlsqr.nAttP=nAttP; + comlsqr.nInstrP=nInstrP; + comlsqr.nInstrPSolved=nInstrPSolved; + comlsqr.nGlobP=nGlobP; + comlsqr.mapNossBefore=mapNossBefore; + comlsqr.mapNossAfter=mapNossAfter; + comlsqr.myid=myid; + comlsqr.nproc=nproc; + comlsqr.mapNoss=mapNoss; + comlsqr.mapNcoeff=mapNcoeff; + comlsqr.multMI=multMI; + comlsqr.debugMode=debugMode; + comlsqr.noCPR=noCPR; + comlsqr.nAttParam=nAttParam; + comlsqr.extConstraint=extConstraint; + comlsqr.nEqExtConstr=nEqExtConstr; + comlsqr.numOfExtStar=numOfExtStar; + comlsqr.barConstraint=barConstraint; + comlsqr.nEqBarConstr=nEqBarConstr; + comlsqr.numOfBarStar=numOfBarStar; + comlsqr.firstStarConstr=firstStarConstr; + comlsqr.lastStarConstr=lastStarConstr; + comlsqr.numOfExtAttCol=numOfExtAttCol; + comlsqr.startingAttColExtConstr=startingAttColExtConstr; + comlsqr.setBound[0]=0; + comlsqr.setBound[1]=nAstroPSolved; + comlsqr.setBound[2]=nAstroPSolved+nAttP; + comlsqr.setBound[3]=nAstroPSolved+nAttP+nInstrPSolved; + comlsqr.nDegFreedomAtt=nDegFreedomAtt; + comlsqr.nAttParAxis=nAttParAxis; + comlsqr.nAttAxes=nAttAxes; + comlsqr.nobs=nobs; + comlsqr.lsInstrFlag=lsInstrFlag; + comlsqr.ssInstrFlag=ssInstrFlag; + comlsqr.nuInstrFlag=nuInstrFlag; + comlsqr.maInstrFlag=maInstrFlag; + comlsqr.myid=myid; + comlsqr.cCDLSAACZP=cCDLSAACZP; + comlsqr.nOfInstrConstr=nOfInstrConstr; + comlsqr.nElemIC=nElemIC; + comlsqr.nElemICLSAL=nElemICLSAL; + comlsqr.nElemICLSAC=nElemICLSAC; + comlsqr.nElemICSS=nElemICSS; + comlsqr.instrConst[0]=instrConst[0]; + comlsqr.instrConst[1]=instrConst[1]; + comlsqr.instrConst[2]=instrConst[2]; + comlsqr.instrConst[3]=instrConst[3]; + comlsqr.nInstrParam=nInstrParam; + comlsqr.nGlobalParam=nGlobalParam; ///////////////////// end buidl map distributed arrays - + // Allocate the memory for the vectors and compute the total memory allocated totmem = 0; nElements = mapNcoeff[myid]; - if (extConstraint) - { - addElementextStar = (lastStarConstr - firstStarConstr + 1) * nAstroPSolved; - addElementAtt = numOfExtAttCol * nAttAxes; + if (extConstraint){ + addElementextStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved; + addElementAtt=numOfExtAttCol*nAttAxes; } - if (barConstraint) - { - addElementbarStar = (lastStarConstr - firstStarConstr + 1) * nAstroPSolved; + if (barConstraint){ + addElementbarStar= (lastStarConstr-firstStarConstr+1)*nAstroPSolved; } - nOfElextObs = addElementextStar + addElementAtt; - nOfElBarObs = addElementbarStar; - comlsqr.nOfElextObs = nOfElextObs; - comlsqr.nOfElBarObs = nOfElBarObs; - - nElements += nOfElextObs * nEqExtConstr + nOfElBarObs * nEqBarConstr + nElemIC; + nOfElextObs=addElementextStar+addElementAtt; + nOfElBarObs=addElementbarStar; + comlsqr.nOfElextObs=nOfElextObs; + comlsqr.nOfElBarObs=nOfElBarObs; - systemMatrix = (double *)calloc(nElements, sizeof(double)); + 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)); + 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)); + 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); + 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 + 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); - + 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(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); - + 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"); - + 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); + 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); + 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); + 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); + 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(myid==0) printf( + "Found %d SM data files in the directory. Reading the coefficients...\n", + (nfileProc)*nproc); + timeToReadFiles=MPI_Wtime(); + + for(ii=0;ii 0) + obsTotal+=nObsxStar; + if(startingStar0) { - 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++; + for(int q=0;q 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++; + if(debugMode) + printf("PE=%d instrStartFreedom=%ld instrEndFreedom=%ld nInstrParam=%d\n",myid,instrStartFreedom,instrEndFreedom,nInstrParam); + if(obsStar==0) + { + obsStar=nObsxStar; + if(currentStar0) + { + for(int q=0;q 0) + 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++; + for(int q=0;q= nDegFreedomAtt - nAttParAxis) - lastFreedom = nDegFreedomAtt - nAttParAxis; - matrixIndex[ii * multMI + (multMI - 1)] = offsetAttParam + lastFreedom; - if (lastFreedom >= endFreedom || lastFreedom >= nDegFreedomAtt - nAttParAxis) - freedomReached = 1; - lastFreedom += nAttParAxis; + 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 + } 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; + lastFreedom=0; + constraintFound[counterConstr][0]=currentStar/1000; + constraintFound[counterConstr][1]=rowInFile; counterConstr++; } - matrixIndex[ii * multMI + (multMI - 1)] = offsetAttParam + lastFreedom; + 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; + + if(!instrFreedomReached && nInstrPSolved) + { + if((obsStar-counterStarObs)<=isConstraint) + { //constraint + for(int kk=0;kkinstrEndFreedom) instrLastFreedom=instrEndFreedom; + instrCol[ii*nInstrPSolved]=instrLastFreedom; + for(int kk=1;kk 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) - { + if((obsStar-counterStarObs)>isConstraint ) + { + for(int q=0;q0) + { + if(ii!=0) offsetConstraint=0; + int foundedConstraint=(obsStar-counterStarObs)+offsetConstraint; + int itis=0; + for(int q=0;q 0) - { - for (int q = 0; q < nConstrLat; q++) - if (constrLatId[q] == currentStar) + changedStar=1; + isConstraint=0; + obsStar=nObsxStar; + if(currentStar0) + { + for(int q=0;q= 2) - randVal = 0.; - if (nAstroPSolved == 5 && i % nAstroPSolved == 0) - randVal = 0.; - if (nAstroPSolved == 5 && i % nAstroPSolved > 2 && j < 3) - randVal = 0.; + exit(err_malloc("attNS",myid)); + if(myid==0){ + for(int i=0;i=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(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==1 || j==4){ + if(i=addElementextStar+2*addElementAtt/nAttAxes) randVal=0.0; } - if (j == 2 || j == 5) - { - if (i < addElementextStar + 2 * addElementAtt / nAttAxes) - randVal = 0.0; + if(j==2 || j==5){ + if(i= 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); + accumulator=(double *) calloc(nEqBarConstr,sizeof(double)); + + for(int j=0;j=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)); + 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 (!computeInstrConstr(comlsqr, instrCoeffConstr, instrColsConstr, instrConstrIlung)) + + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } ////////////////////////// - for (int k = 0; k < nElemIC; k++) - instrCoeffConstr[k] = instrCoeffConstr[k] * wgInstrCoeff; + for(int k=0;k 0) + if(seqStar<=1 && nAstroPSolved>0) { - printf("ERROR PE=%d Only %d star Run not allowed with this PE numbers .n", myid, seqStar); + 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; + comlsqr.VrIdAstroPDim=seqStar; + long tempDimMax=comlsqr.VrIdAstroPDim; long tempVrIdAstroPDimMax; - MPI_Allreduce(&tempDimMax, &tempVrIdAstroPDimMax, 1, MPI_LONG, MPI_MAX, MPI_COMM_WORLD); - - comlsqr.VrIdAstroPDimMax = 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)); - + tempStarSend=(int **) calloc(nproc,sizeof(int *)); + for(int i=0;i 0) + + 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); + 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) - { + if(myid==0) + for(int i=0;id_name); + if(debugMode) + for(ii=0;iid_name); } - if (nFiles <= 0) + if(nFiles <=0) { - printf("error on scandir n=%d\n", nFiles); - MPI_Abort(MPI_COMM_WORLD, 1); + printf("error on scandir n=%d\n",nFiles); + MPI_Abort(MPI_COMM_WORLD,1); } - + ////////////////// Writing FileConstr_GsrSolProps.dat - char fileConstr[512] = ""; + 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]); + fpFilePosConstr=fopen(fileConstr,"w"); + for(int k=0;kd_name,constraintFound[k][1]); + } + if(debugMode) + for(int k=0;kd_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]); + for(int np=1;npd_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]); + if(debugMode) + for(int k=0;kd_name,constraintFound[k][1]); } } fclose(fpFilePosConstr); - } - else - { + } 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); + 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)); + + 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); - + 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)); + vVect = (double *) calloc(nElements, sizeof(double)); if (!vVect) - exit(err_malloc("vVect", myid)); - totmem += nElements * sizeof(double); - - wVect = (double *)calloc(nElements, sizeof(double)); + 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)); + 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)); + 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 - + 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) + if(myid==0) { printf("LOCAL %ld MB of memory allocated on each task.\n", totmem); - printf("TOTAL MB memory allocated= %ld\n", nproc * 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 int aIndex=0; long instrLongIndex[6]; - long iIelem = 0; - - for (ii = 0; ii < mapNoss[myid] * multMI; ii = ii + multMI) + 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) + 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); + 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++) + for(int ns=0;ns= nunkSplit || ncolumn < 0) - { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if (preCondVect[ncolumn] == 0.0) + 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++) + for(int naxis=0;naxis= 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); + 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 + 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) + + + if(nInstrPSolved>0) { - for (int ns = 0; ns < nInstrPSolved; ns++) + for(int ns=0;ns= nunkSplit || ncolumn < 0) + 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("ERROR. PE=%d ii=%ld ",myid, ii); + for(int ke=0;ke= nunkSplit || ncolumn < 0) - { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii+2]=%ld\n", myid, ncolumn, ii, + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if (preCondVect[ncolumn] == 0) + 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]); + 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) + for(int i=0;i0) 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); + 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++) + for(int ns=0;ns= nunkSplit || ncolumn < 0) - { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } preCondVect[ncolumn] += systemMatrix[aIndex] * systemMatrix[aIndex]; - if (preCondVect[ncolumn] == 0.0) + 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); + for(int naxis=0;naxis= 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++; + 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) + if(barConstraint){ + for(int i=0;i0) 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); + 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++) + for(int ns=0;ns= nunkSplit || ncolumn < 0) - { - printf("ERROR. PE=%d ncolumn=%ld ii=%ld matrixIndex[ii]=%ld\n", myid, ncolumn, ii, + 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); + 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); + if(preCondVect[ncolumn]==0.0) + printf("Astrometric: preCondVect[%ld]=0.0\n", ncolumn); aIndex++; } // @@ -2877,328 +2558,291 @@ int main(int argc, char **argv) ///// 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); + if(nElemIC>0){ + for(int i=0;i= 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]); + 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++) + 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 && 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); + if(ii-VrIdAstroPDimMax*nAstroPSolvednAttParam && ii-VrIdAstroPDimMax*nAstroPSolvednAttParam+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) + 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); + 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) + } 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); + 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) + 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); + if(ii-VrIdAstroPDimMax*nAstroPSolvednAttParam && ii-VrIdAstroPDimMax*nAstroPSolvednAttParam+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) + + 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; + 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"); - + 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) + 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); + 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); + precondSystemMatrix(systemMatrix, preCondVect, matrixIndex,instrCol,comlsqr); //systemMatrix is passed to lsqr already conditioned. - if (myid == 0) - { + if(myid==0) { printf("System built, ready to call LSQR.\nPlease wait. System solution is underway..."); } //////////////////// - startTime = MPI_Wtime(); + 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]); + 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 + 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); - + (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); + 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; + + long thetaCol=0, muthetaCol=0; ldiv_t res_ldiv; - long flagTheta = 0, flagMuTheta = 0; - - if (nAstroPSolved == 2) - { - thetaCol = 1; - muthetaCol = 0; + long flagTheta=0, flagMuTheta=0; + + if(nAstroPSolved==2) { + thetaCol=1; + muthetaCol=0; } - else if (nAstroPSolved == 3) - { - thetaCol = 2; - muthetaCol = 0; + else if(nAstroPSolved==3) { + thetaCol=2; + muthetaCol=0; } - else if (nAstroPSolved == 4) - { - thetaCol = 1; - muthetaCol = 3; + else if(nAstroPSolved==4) { + thetaCol=1; + muthetaCol=3; } - else if (nAstroPSolved == 5) - { - thetaCol = 2; - muthetaCol = 4; + 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; - + 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) + if(muthetaCol==0) + for(ii=0; ii localMaxDev) - localMaxDev = dev; + 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) + 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; + epsilon=xSolution[ii]-1.0; + localSumX+=epsilon; + dev=fabs(epsilon); + if(dev>localMaxDev) localMaxDev=dev; } + } - if (idtest) - localSumXsq += epsilon * epsilon; - standardError[ii] *= preCondVect[ii]; + 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; + for(ii=0; iilocalMaxDev) 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; + } 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]; + 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); + if(myid==0) + for(ii=VroffsetAttParam; iilocalMaxDev) localMaxDev=dev; + localSumXsq+=(xSolution[ii]-1.0)*(xSolution[ii]-1.0); } - standardError[ii] *= preCondVect[ii]; + standardError[ii]*=preCondVect[ii]; } //////////// End of de-preconditioning for the shared unknowns - + + //////////////// TEST per verificare le somme di idtest.... da commentare/eliminare /* if(idtest) @@ -3210,7 +2854,7 @@ int main(int argc, char **argv) } */ //////////////// Fine TEST - + free(systemMatrix); free(matrixIndex); free(instrCol); @@ -3219,31 +2863,30 @@ int main(int argc, char **argv) free(wVect); free(preCondVect); - if (idtest) + 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); + 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(myid==0) { - if (idtest) + if(idtest) { - average = sumX / nunk; - rms = pow(sumXsq / nunk - pow(average, 2.0), 0.5); + 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(istop==1000 && wrsol==0) { - if (myid == 0) - { - printf("Reached limit at itn=%d. Execution stopped. CPR written!\n", itn); - ffcont = fopen("CPR_CONT", "w"); + if(myid==0){ + printf("Reached limit at itn=%d. Execution stopped. CPR written!\n",itn); + ffcont=fopen("CPR_CONT","w"); fclose(ffcont); } MPI_Finalize(); @@ -3251,304 +2894,283 @@ int main(int argc, char **argv) } /////////////// ///////////// 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) - { + 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"); + 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(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); + 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]); + 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) + + if(myid==0) { - + ////////// Writing the raw results to the files GsrFinal_GsrSolProps.dat - - fpSolPropsFinal = fopen(filenameSolPropsFinal, "w"); + + + fpSolPropsFinal=fopen(filenameSolPropsFinal,"w"); if (!fpSolPropsFinal) { - printf("Error while open %s\n", filenameSolPropsFinal); - MPI_Abort(MPI_COMM_WORLD, 1); + 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); + 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) + if(nAstroPSolved) { - long testOnStars = 0; - fpAstroR = fopen(filenameAstroResults, "wb"); + long testOnStars=0; + fpAstroR=fopen(filenameAstroResults,"wb"); if (!fpAstroR) { - printf("Error while open %s\n", filenameAstroResults); - MPI_Abort(MPI_COMM_WORLD, 1); + 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)); + 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); + MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } - for (ii = 0; ii < VrIdAstroPDim; ii++) + for(ii=0; ii