Skip to content
solvergaiaSim.c 133 KiB
Newer Older
Fabio Roberto Vitello's avatar
Fabio Roberto Vitello committed
                    for(jj=0; jj< nAstroPSolved;jj++)
                    {
                        xAstro[(ii-offset)*nAstroP+mapAstroP[jj]] = xSolution[ii*nAstroPSolved+jj];
                    }
                }
                fwrite(xAstro,sizeof(double),(VrIdAstroPDimRecv-offset)*nAstroP,fpAstroR);
                testOnStars+=VrIdAstroPDimRecv-offset;
            }
            free(xAstro);
            
            if(testOnStars!=nStar)
            {
                printf("Error on xAstro testOnStar =%ld not equal to nmStar=%ld\n",testOnStars,nStar);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            testOnStars=0;
            standardErrorAstro = (double *) calloc(VrIdAstroPDimMax*nAstroP, sizeof(double));
            if (!standardErrorAstro)
            {
                printf("Error allocating standardErrorAstro\n");
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            for(ii=0; ii<VrIdAstroPDim; ii++)
            {
                for(jj=0; jj< nAstroPSolved;jj++)
                {
                    standardErrorAstro[ii*nAstroP+mapAstroP[jj]] =standardError[ii*nAstroPSolved+jj];
                }
            }
            fwrite(standardErrorAstro,sizeof(double),VrIdAstroPDim*nAstroP,fpAstroR);
            testOnStars+=VrIdAstroPDim;
            for(kk=1;kk<nproc;kk++)
            {
                offset=0;
                MPI_Recv(&VrIdAstroPDimRecv, 1, MPI_LONG, kk, 12, MPI_COMM_WORLD, &statusMpi);
                MPI_Recv(standardError, VrIdAstroPDimRecv*nAstroPSolved , MPI_DOUBLE, kk, 13,MPI_COMM_WORLD, &statusMpi);
                if(comlsqr.mapStar[kk][0]==comlsqr.mapStar[kk-1][1]) offset=1;
                for(ii=0+offset; ii<VrIdAstroPDimRecv; ii++)
                {
                    for(jj=0; jj< nAstroPSolved;jj++)
                    {
                        standardErrorAstro[(ii-offset)*nAstroP+mapAstroP[jj]] = standardError[ii*nAstroPSolved+jj];
                    }
                }
                fwrite(standardErrorAstro,sizeof(double),(VrIdAstroPDimRecv-offset)*nAstroP,fpAstroR);
                testOnStars+=VrIdAstroPDimRecv-offset;
            }
            free(standardErrorAstro);
            if(testOnStars!=nStar)
            {
                printf("Error on standardErrorAstro testOnStar =%ld not equal to nmStar=%ld\n",testOnStars,nStar);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
            fclose(fpAstroR);
            //////////////////////////////////////////////////////////
        }
        ////////////////////// writing  GsrAttitudeParamSolution.bin
        fpAttR=fopen(filenameAttResults,"wb");
        if (!fpAttR)
        {
            printf("Error while open %s\n",filenameAttResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        fwrite(&sphereId,sizeof(long),1,fpAttR);
        fwrite(&xSolution[VroffsetAttParam],sizeof(double),nAttParam,fpAttR);
        fwrite(&standardError[VroffsetAttParam],sizeof(double),nAttParam,fpAttR);
        fclose(fpAttR);
        //////////////////////////////////////////////////////////
        
        ////////////////////// writing GsrInstrParamSolution.bin
        fpInstrR=fopen(filenameInstrResults,"wb");
        if (!fpInstrR)
        {
            printf("Error while open %s\n",filenameInstrResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        double *xInstr,*seInstr;
        xInstr=(double *) calloc(nInstrParamTot, sizeof(double));
        seInstr=(double *) calloc(nInstrParamTot, sizeof(double));
        if(nInstrPSolved!=nInstrP)
        {
            int fixedOffsetCMag = nCCDs;
            int fixedOffsetCnu = fixedOffsetCMag+nFoVs*nCCDs;
            int fixedOffsetCdelta_eta = fixedOffsetCnu+nCCDs*nPixelColumns;
            int fixedOffsetCDelta_eta = fixedOffsetCdelta_eta+3*nFoVs*nCCDs*nTimeIntervals;
            int fixedOffsetCdelta_zeta = fixedOffsetCDelta_eta+nCCDs*nPixelColumns;
            
            int counterInstr=0;
            if(maInstrFlag){
                for(int k=0;k<fixedOffsetCMag;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(nuInstrFlag){
                for(int k=fixedOffsetCMag;k<fixedOffsetCnu;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(ssInstrFlag){
                for(int k=fixedOffsetCnu;k<fixedOffsetCdelta_eta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(lsInstrFlag){
                for(int k=fixedOffsetCdelta_eta;k<fixedOffsetCDelta_eta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(ssInstrFlag){
                for(int k=fixedOffsetCDelta_eta;k<fixedOffsetCdelta_zeta;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(lsInstrFlag){
                for(int k=fixedOffsetCdelta_zeta;k<nInstrParamTot;k++){
                    xInstr[k]=xSolution[VroffsetAttParam+nAttParam+counterInstr];
                    seInstr[k]=standardError[VroffsetAttParam+nAttParam+counterInstr];
                    counterInstr++;
                }
            }
            if(counterInstr!=nInstrParam) {
                printf("SEVERE ERROR: counterInstr=%ld does not match nInstrParam=%ld while filling in xInstr and seInstr.\n", counterInstr, nInstrParam);
                MPI_Abort(MPI_COMM_WORLD,1);
                exit(EXIT_FAILURE);
            }
        }
        fwrite(&sphereId,sizeof(long),1,fpInstrR);
        fwrite(xInstr,sizeof(double),nInstrParamTot,fpInstrR);
        fwrite(seInstr,sizeof(double),nInstrParamTot,fpInstrR);
        fclose(fpInstrR);
        free(xInstr);
        free(seInstr);
        //////////////////////////////////////////////////////////
        
        ////////////////////// writing GsrGlobalParamSolution.bin
        fpGlobR=fopen(filenameGlobalResults,"wb");
        if (!fpGlobR)
        {
            printf("Error while open %s\n",filenameGlobalResults);
            MPI_Abort(MPI_COMM_WORLD,1);
            exit(EXIT_FAILURE);
        }
        fwrite(&sphereId,sizeof(long),1,fpGlobR);
        fwrite(&xSolution[VroffsetAttParam+nAttParam+nInstrParam],sizeof(double),nGlobalParam,fpGlobR);
        fwrite(&standardError[VroffsetAttParam+nAttParam+nInstrParam],sizeof(double),nGlobalParam,fpGlobR);
        fclose(fpGlobR);
        //////////////////////////////////////////////////////////
        
    } //if (myid==0)
    else
    {
        //////////////// send for AstroResults.bin
        if(nAstroPSolved)
        {
            MPI_Send(&VrIdAstroPDim, 1, MPI_LONG, 0, 10, MPI_COMM_WORLD);
            mpi_send(xSolution, VrIdAstroPDim*nAstroPSolved, MPI_DOUBLE, 0, 11,MPI_COMM_WORLD);
            MPI_Send(&VrIdAstroPDim, 1, MPI_LONG, 0, 12, MPI_COMM_WORLD);
            MPI_Send(standardError, VrIdAstroPDim*nAstroPSolved , MPI_DOUBLE, 0, 13,MPI_COMM_WORLD);
        }
        //////////////////////////////////////////////////////////
    }
    MPI_Barrier(MPI_COMM_WORLD);
    //  tolti i free() mancanti
    free(xSolution);
    free(standardError);
    
    if(myid==0) printf("\nEnd run.\n");
    seconds[3] = time(NULL);
    tot_sec[2] = seconds[3] - seconds[2];
    tot_sec[3] = seconds[3] - seconds[0];
    
    if(myid==0) printf("time before lsqr in sec: %ld\ntime for lsqr: %ld\n", tot_sec[0],
                       tot_sec[1]);
    if(myid==0) printf("time after lsqr in sec: %ld\ntime TOTAL: %ld\n", tot_sec[2],
                       tot_sec[3]);
    
    MPI_Finalize();
    exit(EXIT_SUCCESS);
}