#include #include #include "util.h" #include #include #ifdef OMP #include #endif #define __MSGSIZ_MAX 100000 /* The function instrIndexIdToColIndexGlobal gets *instrIndexPointer and *instrConst, and computes *matrixIndexIntsr, which contains the columns of the 6 instrument parameters in the design matrix. The input pointers refers to: a) which FoV, CCD, PixelColumn, and TimeInterval each observation has occurred according to the following schema instrIndexPointer[0]=FoV instrIndexPointer[1]=CCD instrIndexPointer[2]=PixelColumn instrIndexPointer[3]=TimeInterval b) the characteristics of the instrument according to the following schema instrConst[0]=nFoVs instrConst[1]=nCCDs instrConst[2]=nPixelColumns instrConst[3]=nTimeIntervals For performance reasons, it accepts some (five) pre-calculated constant offsets, i.e.: 1) offsetCMag = nCCDs = instrConst[1] 2) offsetCnu = nCCDs*(1+nFoVs) = instrConst[1]*(1+instrConst[0]) 3) offsetCdelta_eta = nCCDs*(1+nFoVs+nPixelColumns) = instrConst[1]*(1+instrConst[0]+instrConst[2]) 4) offsetCDelta_eta_1 = nCCDs*(1+nFoVs*(1+nTimeIntervals)+nPixelColumns) = instrConst[1]*(1+instrConst[0]*(1+instrConst[3])+instrConst[2]) 5) offsetCDelta_eta_2 = nCCDs*(1+nFoVs*(1+2*nTimeIntervals)+nPixelColumns) = instrConst[1]*(1+instrConst[0]*(1+2*instrConst[3])+instrConst[2]) */ void instrIndexIdToColIndexGlobal(int* instrIndexPointer, int* instrConst,int totrows ,struct comData comlsqr, int* instrCols) { long relPos_ls; int nInstrParam=comlsqr.nInstrParam; short nInstrPSolved=comlsqr.nInstrPSolved; int maInstrFlag=comlsqr.maInstrFlag; int nuInstrFlag=comlsqr.nuInstrFlag; int ssInstrFlag=comlsqr.ssInstrFlag; int lsInstrFlag=comlsqr.lsInstrFlag; int nFoVs =1+instrConst[0]; int nCCDs = instrConst[1]; int nPixCols = instrConst[2]; int nTInts = instrConst[3]; for(int jj=0;jj=nInstrParam ||instrCols[jj*nInstrPSolved+k]<0 ){ printf("SEVERE ERROR on instrCols[%d]=%d > nInstrparam=%d\n",jj*nInstrPSolved+k,instrCols[jj*nInstrPSolved+k],nInstrParam); MPI_Abort(MPI_COMM_WORLD, 1); exit(EXIT_FAILURE); } } } return; } void ColIndexToinstrIndexIdGlobal(int* instrIndexPointer, int* instrConst,int totrows ,struct comData comlsqr, int* instrCols) { long relPos_Delta_eta; int testConstr=0; for(int jj=0;jjd_name, "dpccu3dpctavugsrgsrsystemrow", 28) == 0) ? 1 : 0); } int selextConstrStar(const struct dirent *a) { return ((strstr(a->d_name, "nullspaceastrofit") != NULL) ? 1 : 0); } int selextConstrAtt(const struct dirent *a) { return ((strstr(a->d_name, "nullspaceattitudefit") != NULL) ? 1 : 0); } int selbarConstrStar(const struct dirent *a) { return ((strstr(a->d_name, "barconstrastrofit") != NULL) ? 1 : 0); } int selSM(const struct dirent *a) { return ((strstr(a->d_name, "_SM.bin") != NULL) ? 1 : 0); } int selKT(const struct dirent *a) { return ((strstr(a->d_name, "_KT.bin") != NULL) ? 1 : 0); } int selII(const struct dirent *a) { return ((strstr(a->d_name, "_II.bin") != NULL) ? 1 : 0); } int selMI(const struct dirent *a) { return ((strstr(a->d_name, "_MI.bin") != NULL) ? 1 : 0); } int selAll(const struct dirent *a) { return ((strstr(a->d_name, ".bin") != NULL) ? 1 : 0); } int selGSB(const struct dirent *a) { return ((strstr(a->d_name, ".gsb") != NULL) ? 1 : 0); } int selLastGSB(const struct dirent *a) { return ((strstr(a->d_name, "CPRLast") != NULL) ? 1 : 0); } /* This function returns the values associated to the FoV, the CCD, the PixelColumn and the TimeInterval coded in instrOutput, storing them in instrOutput[i], with i=0...3 respectively. NB: the information about the FoV is coded as 0 or 1 to use a single bit in the mask, but it should be 1 or 2 respectively, so one has to add 1 to the demasking result in order to obtain the correct value in instrOutput[0]. */ void instrDeMask(long instrInput, int acSc ,int* instrOutput) { int CCD_OFFSET = 1; int PIXEL_OFFSET = 9; int TIME_OFFSET = 20; int FOV_MASK = 0x01; int CCD_MASK = 0xFF; int PIXEL_MASK = 0x7FF; int TIME_MASK = 0x7FF; instrOutput[0]= (int) ((instrInput & FOV_MASK)); instrOutput[1]= (int) ((instrInput >> CCD_OFFSET) & CCD_MASK); instrOutput[2]= (int) ((instrInput >> PIXEL_OFFSET) & PIXEL_MASK); instrOutput[3]= (int) ((instrInput >> TIME_OFFSET) & TIME_MASK); instrOutput[4]= acSc; return; } void restartSetup(int *itn, double *knownTerms, double *beta, double *alpha, double *vVect, double *anorm, double *rhobar, double *phibar, double *wVect, double *xSolution, double *standardError, double *dnorm, double *sn2, double *cs2, double *z, double *xnorm1, double *res2, int *nstop, struct comData comlsqr) { FILE *fChekPointPtr; int rank, status,size,cksize; int existCPR; int globalCPR; char lastBinName[80]; char rankStr[8]; long int * mapNoss, * mapNcoeff, nunkSplit; long dummyLong[1]; int dummyInt[1]; struct dirent **namelistGSB; long replicaKT=comlsqr.nEqExtConstr+comlsqr.nEqBarConstr+comlsqr.nOfInstrConstr; long localVect=comlsqr.VrIdAstroPDimMax*comlsqr.nAstroPSolved; long replicaVect=comlsqr.nAttParam+comlsqr.nInstrParam+comlsqr.nGlobalParam; mapNcoeff=comlsqr.mapNcoeff; mapNoss=comlsqr.mapNoss; nunkSplit=comlsqr.nunkSplit; MPI_Comm_size(MPI_COMM_WORLD, &size); MPI_Comm_rank(MPI_COMM_WORLD, &rank); sprintf (rankStr, "%d", rank); int nFilesLastGsb= scandir(".", &namelistGSB, selLastGSB, alphasort); cksize=size; sprintf(lastBinName,"GaiaGsrCPRLast_%06d.gsb", rank); existCPR=0; globalCPR=0; fChekPointPtr=NULL; if( (fChekPointPtr=fopen(lastBinName,"rb")) !=NULL ) // If GaiaGsrCPRLast_#PE.gsb does not exist... { existCPR=1; fread(&cksize,sizeof(int),1,fChekPointPtr); } MPI_Allreduce(&existCPR, &globalCPR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); if( (globalCPR!=0 && globalCPR!=size) || size!=cksize) { status=1; if(rank==0) printf("PE=%d SEVERE ERROR on CPR files. MPI Abort\n",rank); // ...we have MPI_Abort(MPI_COMM_WORLD,status); } if(existCPR) { int nFilesLastGsb= scandir(".", &namelistGSB, selLastGSB, alphasort); if (nFilesLastGsb != size){ status=1; if(rank==0) printf("PE=%d SEVERE ERROR on CPR files, gsb files =%d not equal to size=%d. MPI Abort\n",rank,nFilesLastGsb,size); // ...we have MPI_Abort(MPI_COMM_WORLD,status); } fread(dummyLong,sizeof(long),1,fChekPointPtr); if(dummyLong[0] != mapNoss[rank]){ printf("PE=%d CPR Severe error: Invalid read mapNoss=%ld but must to be=%ld\n",rank,dummyLong[0],mapNoss[rank]); status=1; MPI_Abort(MPI_COMM_WORLD,status); } fread(dummyLong,sizeof(long),1,fChekPointPtr); if(dummyLong[0] != replicaKT){ printf("PE=%d CPR Severe Error. Read value replicaKT=%ld is different from computed replicaKt=%ld\n",rank,dummyLong[0],replicaKT); status=1; MPI_Abort(MPI_COMM_WORLD,status); } fread(dummyLong,sizeof(long),1,fChekPointPtr); if(dummyLong[0] !=localVect){ printf("PE=%d CPR Severe Error. Read value localVect=%ld is different from computed localVect=%ld\n",rank,dummyLong[0],localVect); status=1; MPI_Abort(MPI_COMM_WORLD,status); } fread(dummyInt,sizeof(int),1,fChekPointPtr); if(dummyInt[0] !=replicaVect){ printf("PE=%d CPR Severe Error. Read value replicaVect=%ld is different from computed replicaVect=%ld/n",rank,dummyLong[0],replicaVect); status=1; MPI_Abort(MPI_COMM_WORLD,status); } fread(itn, sizeof(int),1,fChekPointPtr); fread(knownTerms, sizeof(double),mapNoss[rank]+comlsqr.nEqExtConstr+comlsqr.nEqBarConstr+comlsqr.nOfInstrConstr,fChekPointPtr); fread(beta, sizeof(double),1,fChekPointPtr); fread(alpha, sizeof(double),1,fChekPointPtr); fread(vVect, sizeof(double),nunkSplit,fChekPointPtr); fread(anorm, sizeof(double),1,fChekPointPtr); fread(rhobar, sizeof(double),1,fChekPointPtr); fread(phibar, sizeof(double),1,fChekPointPtr); fread(wVect, sizeof(double),nunkSplit,fChekPointPtr); fread(xSolution, sizeof(double),nunkSplit,fChekPointPtr); fread(standardError, sizeof(double),nunkSplit,fChekPointPtr); fread(dnorm, sizeof(double),1,fChekPointPtr); fread(sn2, sizeof(double),1,fChekPointPtr); fread(cs2, sizeof(double),1,fChekPointPtr); fread(z, sizeof(double),1,fChekPointPtr); fread(xnorm1, sizeof(double),1,fChekPointPtr); fread(res2, sizeof(double),1,fChekPointPtr); fread(nstop, sizeof(int),1,fChekPointPtr); // printf("PE=%d RCPR itn =%d\n",rank,*itn); // printf("PE=%d RCPR knownTerms ==> %f %f\n",rank,knownTerms[0], knownTerms[mapNoss[rank]-1]); // printf("PE=%d RCPR beta =%f\n",rank,*beta); // printf("PE=%d RCPR alpha =%f\n",rank,*alpha); // printf("PE=%d RCPR vVect ==> %f %f\n",rank,vVect[0], vVect[nunkSplit-1]); // printf("PE=%d RCPR anorm =%f\n",rank,*anorm); // printf("PE=%d RCPR rhobar =%f\n",rank,*rhobar); // printf("PE=%d RCPR phibar =%f\n",rank,*phibar); // printf("PE=%d RCPR wVect ==> %f %f\n",rank,wVect[0], wVect[nunkSplit-1]); // printf("PE=%d RCPR xSolution ==> %f %f\n",rank,xSolution[0], xSolution[nunkSplit-1]); // printf("PE=%d RCPR standardError ==> %f %f\n",rank,standardError[0], standardError[nunkSplit-1]); // printf("PE=%d RCPR dnorm =%f\n",rank,*dnorm); // printf("PE=%d RCPR sn2 =%f\n",rank,*sn2); // printf("PE=%d RCPR cs2 =%f\n",rank,*cs2); // printf("PE=%d RCPR z =%f\n",rank,*z); // printf("PE=%d RCPR xnorm1 =%f\n",rank,*xnorm1); // printf("PE=%d RCPR res2 =%f\n",rank,*res2); // printf("PE=%d RCPR nstop =%d\n",rank,*nstop); fclose(fChekPointPtr); } } void writeCheckPoint(int itn, double *knownTerms, double beta, double alpha, double *vVect, double anorm, double rhobar, double phibar, double *wVect, double *xSolution, double *standardError, double dnorm, double sn2, double cs2, double z, double xnorm1, double res2, int nstop, struct comData comlsqr) { FILE *fChekPointPtr; int rank, size,status; int noCPR; int globalCPR; long replicaKT; long VrIdAstroPDimMax; long localVect; int replicaVect; char prevBinName[80]; char lastBinName[80]; char rankStr[8]; long int * mapNoss, * mapNcoeff, nunkSplit; mapNcoeff=comlsqr.mapNcoeff; mapNoss=comlsqr.mapNoss; nunkSplit=comlsqr.nunkSplit; VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; localVect=VrIdAstroPDimMax*comlsqr.nAstroPSolved; replicaVect=comlsqr.nAttParam+comlsqr.nInstrParam+comlsqr.nGlobalParam; MPI_Comm_size(MPI_COMM_WORLD, &size); MPI_Comm_rank(MPI_COMM_WORLD, &rank); sprintf (rankStr, "%d", rank); sprintf(prevBinName,"GaiaGsrCPRPrev_%06d.gsb", rank); sprintf(lastBinName,"GaiaGsrCPRLast_%06d.gsb", rank); noCPR=0; globalCPR=0; fChekPointPtr=NULL; replicaKT= comlsqr.nEqExtConstr+comlsqr.nEqBarConstr+comlsqr.nOfInstrConstr; if((fChekPointPtr=fopen(lastBinName,"rb"))==NULL) // If GaiaGsrCPRLast_#PE.gsb does not exist... { if(rank==0) printf("PE=%d No checkpoint yet. Writing a new one:\n",rank); // ...we have to write it from scratch noCPR=1; } else fclose(fChekPointPtr); MPI_Allreduce(&noCPR, &globalCPR,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); if(globalCPR!=0) { if((fChekPointPtr=fopen(lastBinName,"wb"))==NULL) { status=1; if(rank==0) printf("PE=%d SEVERE ERROR No CPR file. MPI Abort\n",rank); // ...we have MPI_Abort(MPI_COMM_WORLD,status); } } else { // else... // Remove penultimate checkpoint // Move last checkpoint files to penultimate if(rename(lastBinName,prevBinName) == -1) { status=1; printf("PE=%d SEVERE ERROR Cannot rename previous CPR file %s. MPI Abort.\n",rank, lastBinName); MPI_Abort(MPI_COMM_WORLD,status); } } // Can now write new last checkpoint if((fChekPointPtr=fopen(lastBinName,"wb"))==NULL) // If GaiaGsrCPRLast.gst does not exist... { status=1; printf("PE=%d Severe Warning no Checkpoint file will be produced:\n",rank); // ...we have to write it from scratch MPI_Abort(MPI_COMM_WORLD,status); } else { // else... fwrite(&size, sizeof(int),1,fChekPointPtr); fwrite(&mapNoss[rank],sizeof(long),1,fChekPointPtr); fwrite(&replicaKT,sizeof(long),1,fChekPointPtr); fwrite(&localVect, sizeof(long),1,fChekPointPtr); fwrite(&replicaVect, sizeof(int),1,fChekPointPtr); fwrite(&itn, sizeof(int),1,fChekPointPtr); fwrite(knownTerms, sizeof(double),mapNoss[rank]+comlsqr.nEqExtConstr+comlsqr.nEqBarConstr+comlsqr.nOfInstrConstr,fChekPointPtr); fwrite(&beta, sizeof(double),1,fChekPointPtr); fwrite(&alpha, sizeof(double),1,fChekPointPtr); fwrite(vVect, sizeof(double),nunkSplit,fChekPointPtr); fwrite(&anorm, sizeof(double),1,fChekPointPtr); fwrite(&rhobar, sizeof(double),1,fChekPointPtr); fwrite(&phibar, sizeof(double),1,fChekPointPtr); fwrite(wVect, sizeof(double),nunkSplit,fChekPointPtr); fwrite(xSolution, sizeof(double),nunkSplit,fChekPointPtr); fwrite(standardError, sizeof(double),nunkSplit,fChekPointPtr); fwrite(&dnorm, sizeof(double),1,fChekPointPtr); fwrite(&sn2, sizeof(double),1,fChekPointPtr); fwrite(&cs2, sizeof(double),1,fChekPointPtr); fwrite(&z, sizeof(double),1,fChekPointPtr); fwrite(&xnorm1, sizeof(double),1,fChekPointPtr); fwrite(&res2, sizeof(double),1,fChekPointPtr); fwrite(&nstop, sizeof(int),1,fChekPointPtr); // printf("PE=%d WCPR size =%d\n",rank,size); // printf("PE=%d WCPR itn =%d\n",rank,itn); // printf("PE=%d WCPR knownTerms ==> %f %f\n",rank,knownTerms[0], knownTerms[mapNoss[rank]-1]); // printf("PE=%d WCPR beta =%f\n",rank,beta); // printf("PE=%d WCPR alpha =%f\n",rank,alpha); // printf("PE=%d WCPR vVect ==> %f %f\n",rank,vVect[0], vVect[nunkSplit-1]); // printf("PE=%d WCPR anorm =%f\n",rank,anorm); // printf("PE=%d WCPR rhobar =%f\n",rank,rhobar); // printf("PE=%d WCPR phibar =%f\n",rank,phibar); // printf("PE=%d WCPR wVect ==> %f %f\n",rank,wVect[0], wVect[nunkSplit-1]); // printf("PE=%d WCPR xSolution ==> %f %f\n",rank,xSolution[0], xSolution[nunkSplit-1]); // printf("PE=%d WCPR standardError ==> %f %f\n",rank,standardError[0], standardError[nunkSplit-1]); // printf("PE=%d WCPR dnorm =%f\n",rank,dnorm); // printf("PE=%d WCPR sn2 =%f\n",rank,sn2); // printf("PE=%d WCPR cs2 =%f\n",rank,cs2); // printf("PE=%d WCPR z =%f\n",rank,z); // printf("PE=%d WCPR xnorm1 =%f\n",rank,xnorm1); // printf("PE=%d WCPR res2 =%f\n",rank,res2); // printf("PE=%d WCPR nstop =%d\n",rank,nstop); fclose(fChekPointPtr); if(rank==0)printf("Checkpoint writing ended successfully.\n"); } } void SumCirc(double *vectToSum,struct comData comlsqr) { int rank, size, npeSend, npeRecv; MPI_Status status; MPI_Request req2,req3; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); int nMov=2; if(size==2) nMov=1; if(size==1) return; double *tempSendBuf, *tempRecvBuf; int tempSendIdBuf[2],tempRecvIdBuf[2]; tempSendBuf=(double *) calloc(comlsqr.multMI*comlsqr.nAstroPSolved,sizeof(double)); tempRecvBuf=(double *) calloc(comlsqr.multMI*comlsqr.nAstroPSolved,sizeof(double)); npeSend=rank+1; if(npeSend==size) npeSend=0; npeRecv=rank-1; if(npeRecv<0) npeRecv=size-1; tempSendIdBuf[0]=comlsqr.mapStar[rank][0]; //strating star tempSendIdBuf[1]=comlsqr.mapStar[rank][1]; //ending star for(int i=0;imyid; comlsqr->nthreads=1; #ifdef OMP comlsqr->nthreads = omp_get_max_threads(); #endif int nthreads=comlsqr->nthreads; /// Prepare the structure for the division of the for cycle in aprod mode=2 comlsqr->mapForThread=(long **) calloc(nthreads,sizeof(long *)); for(int n=0;nmapForThread[n]=(long *) calloc(3,sizeof(long)); int nElements=comlsqr->mapNoss[myid]/nthreads; comlsqr->mapForThread[0][0]=0; comlsqr->mapForThread[0][1]=nElements/2; comlsqr->mapForThread[0][2]=nElements; if(comlsqr->mapNoss[myid]%nthreads>0) comlsqr->mapForThread[0][2]++; for(int n=1;nmapForThread[n][0]=comlsqr->mapForThread[n-1][2]; comlsqr->mapForThread[n][1]=comlsqr->mapForThread[n][0]+nElements/2; comlsqr->mapForThread[n][2]=comlsqr->mapForThread[n][0]+nElements; if(comlsqr->mapNoss[myid]%nthreads>n) comlsqr->mapForThread[n][2]++; } comlsqr->mapForThread[nthreads-1][2]=comlsqr->mapNoss[myid]; ////////////////////////////////// // Check for the NOT super-imposed stars at half cycle if(comlsqr->nAstroPSolved>0){ int smpFailed=0; for(int n=1;nmapForThread[n-1][2]-1)]==matrixIndex[2*comlsqr->mapForThread[n][0]]) { if(comlsqr->mapForThread[n][0]==comlsqr->mapForThread[n][2]) { smpFailed=1; printf("PE=%d. SEVERE WARNING. Smp not applicable. mapForThread[%d][0] =%ld and mapForThread[%d][2]=%ld\n",myid, n,comlsqr->mapForThread[n][0],n,comlsqr->mapForThread[n][1]); break; } comlsqr->mapForThread[n][0]++; comlsqr->mapForThread[n-1][2]++; if(smpFailed) break; } } if(smpFailed) { printf("UTIL: SEVERE WARNING PE=%d smpFailed\n",myid); comlsqr->mapForThread[0][0]=0; comlsqr->mapForThread[0][1]=comlsqr->mapNoss[myid]; comlsqr->mapForThread[0][2]=comlsqr->mapForThread[0][1]; for(int n=1;nmapForThread[n][0]=comlsqr->mapNoss[myid]; comlsqr->mapForThread[n][1]=comlsqr->mapNoss[myid]; comlsqr->mapForThread[n][2]=comlsqr->mapNoss[myid]; } } } ///// if(comlsqr->myid==0) printf("\n\nRunning with OMP: nthreads=%d\n\n",nthreads); } void blocksAttIndep(long int *matrixIndex,struct comData *comlsqr) { //ATTENZIONE NON E? CORRETTA E NON VA MAI USATA SE nAstroPSolved=0 int myid=comlsqr->myid; comlsqr->nthreads=1; #ifdef OMP comlsqr->nthreads = omp_get_num_threads(); #endif int nthreads=comlsqr->nthreads; comlsqr->nSubsetAtt=1; comlsqr->NOnSubsetAtt=1; if(nthreads==1) return; int dependancyFound; //Ogni sottointervallo è diviso in nthreads. Ogni volta che trovo una dipenedenza di indice nel singolo sottointervallo (ne basta uno) moltiplico per due i sottointervalli in modo da ridurre la probabilità di dipendenza. L'indipendenza dell'indice va cercata nel singolo sottointervallo while(1){ // ogni volta su questo while moltiplico x 2 i sottointervalli dependancyFound=0; for(int i=0;inSubsetAtt;i++) //ciclo su tutti i sotto intervall { long totalEleInBlock=comlsqr->mapNoss[myid]/comlsqr->nSubsetAtt; //numero di elementi totali inclusi in tutte le thread nel blocco TBV long totalEleInBlockThread= (comlsqr->mapNoss[myid]/comlsqr->nSubsetAtt)/nthreads; //elementi nella singola thread TBV for(int nsubBlock=0;nsubBlocknAttParAxis || indexFound>kindexFound+comlsqr->nAttParAxis)) { dependancyFound=1; break; } if(dependancyFound) break; } //for k if(dependancyFound) break; }// for j if(dependancyFound) break; }// for nsubBlock if(dependancyFound) break; }// for i if(dependancyFound) { comlsqr->nSubsetAtt=comlsqr->nSubsetAtt*2; if(comlsqr->nSubsetAtt>256) break; } else break; }// while if(dependancyFound) { printf("PE=%d WARNING impossible to find on 256 subSet index independancy for Attitude Parameters\n",myid); comlsqr->NOnSubsetAtt=1; //variabile che indica che MAI si ha indipendenza indici fino a 256 sottointervalli } else { printf("PE=%d Attitude index independancy with %d nSubsetAtt\n",myid,comlsqr->nSubsetAtt); comlsqr->NOnSubsetAtt=0; //variabile che indica che ABBIAMO indipendenza indici fino a 256 sottointervalli con comlsrq->nSubsetAtt sottointervalli } } // This function computes the product of system matrix by precondVect. This avoids to compute the produsct in aprod for each iteration. void precondSystemMatrix(double *systemMatrix, double *preCondVect, long int *matrixIndex,int *instrCol,struct comData comlsqr) { int myid; long int *mapNoss, *mapNcoeff; long int j, l=0; int ii; int setBound[4]; myid=comlsqr.myid; mapNcoeff=comlsqr.mapNcoeff; mapNoss=comlsqr.mapNoss; short nAstroPSolved=comlsqr.nAstroPSolved; short nInstrPSolved=comlsqr.nInstrPSolved; long nparam=comlsqr.parOss; int multMI=comlsqr.multMI; short nAttParAxis=comlsqr.nAttParAxis; long counterAxis=0, counterInstr=0; long nDegFredoomAtt=comlsqr.nDegFreedomAtt; long VrIdAstroPDimMax=comlsqr.VrIdAstroPDimMax; long offsetAttParam=comlsqr.offsetAttParam; long offsetInstrParam=comlsqr.offsetInstrParam; long offsetGlobParam=comlsqr.offsetGlobParam; int extConstraint=comlsqr.extConstraint; int nEqExtConstr=comlsqr.nEqExtConstr; int numOfExtStar=comlsqr.numOfExtStar; int barConstraint=comlsqr.barConstraint; int nEqBarConstr=comlsqr.nEqBarConstr; int numOfBarStar=comlsqr.numOfBarStar; int numOfExtAttCol=comlsqr.numOfExtAttCol; int startingAttColExtConstr=comlsqr.startingAttColExtConstr; short nAttAxes=comlsqr.nAttAxes; int nElemIC=comlsqr.nElemIC; long VroffsetAttParam=comlsqr.VroffsetAttParam; setBound[0]=comlsqr.setBound[0]; setBound[1]=comlsqr.setBound[1]; setBound[2]=comlsqr.setBound[2]; setBound[3]=comlsqr.setBound[3]; for(long i=0;i=setBound[0] && ii< setBound[1]) { if(ii==setBound[0]) { long numOfStarPos=matrixIndex[i*multMI]/nAstroPSolved; j=(numOfStarPos-comlsqr.mapStar[myid][0])*nAstroPSolved; } else j++; } if(ii>=setBound[1] && ii< setBound[2]) { if(((ii-setBound[1]) % nAttParAxis)==0) { j=matrixIndex[multMI*i+multMI-1]+counterAxis*nDegFredoomAtt+(VrIdAstroPDimMax*nAstroPSolved-offsetAttParam); counterAxis++; } else j++; } if(ii>=setBound[2] && ii< setBound[3]) { j=offsetInstrParam+instrCol[i*nInstrPSolved+counterInstr]+(VrIdAstroPDimMax*nAstroPSolved-offsetAttParam); counterInstr++; } if(ii>=setBound[3]) { if(ii==comlsqr.setBound[3]) j=offsetGlobParam+(VrIdAstroPDimMax*nAstroPSolved-offsetAttParam); else j++; } systemMatrix[l]=systemMatrix[l]*preCondVect[j]; l++; }//for(ii }//for i if(extConstraint){ for(int i=0;i0){ for(int i=0;i0) { if(lcount<=count) { count=lcount; lcount=0; } else lcount=lcount - count; MPI_Allreduce(&source[ncyc*chunch], &dest[ncyc*chunch], count, datatype, op, comm); ncyc++; } } void mpi_recv(double *source, long int lcount, MPI_Datatype datatype, int peSource, int tag, MPI_Comm comm,MPI_Status *status) { int ncyc=0; int count=__MSGSIZ_MAX; int chunch=__MSGSIZ_MAX; int localTag=tag; while (lcount>0) { if(lcount<=count) { count=lcount; lcount=0; } else lcount=lcount - count; MPI_Recv(&source[ncyc*chunch], count, datatype, peSource, localTag,comm,status); ncyc++; localTag+=100; } } void mpi_send(double *source, long int lcount, MPI_Datatype datatype, int peDest, int tag, MPI_Comm comm) { int ncyc=0; int count=__MSGSIZ_MAX; int chunch=__MSGSIZ_MAX; int localTag=tag; while (lcount>0) { if(lcount<=count) { count=lcount; lcount=0; } else lcount=lcount - count; MPI_Send(&source[ncyc*chunch], count, datatype, peDest, localTag,comm); ncyc++; localTag+=100; } } /* Generates a pseudo-random number having a gaussian distribution * with mean ave e rms sigma. * The init2 parameter is used only when the the pseudo-random * extractor is the ran2() from Numerical Recipes instead of the * standard rand() system function. */ double gauss(double ave, double sigma, long init2) { int i; double rnd; rnd=0.0; for(i=1; i<=12; i++) // comment the following line and uncomment the next one // to use the system rountine for random numbers rnd += ran2(&init2); // rnd += ((double) rand()/RAND_MAX); rnd -= 6.0; rnd = ave+sigma*rnd; return rnd; } /* From "Numerical Recipes in C". Generates random numbers. * Requires a pointer to long as seed and gives a double as the result. */ double ran2(long *idum) /* Long period (> 2 . 10 18 ) random number generator of L'Ecuyer with Bays-Durham shu.e and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter idum between successive deviates in a sequence. RNMX should approximate the largest oating value that is less than 1. */ { int j; long k; static long idum2=123456789; static long iy=0; static long iv[NTAB]; double temp; if (*idum <= 0) { // Initialize. if (-(*idum) < 1) *idum=1; // Be sure to prevent idum = 0. else *idum = -(*idum); idum2=(*idum); for (j=NTAB+7;j>=0;j--) { // Load the shu.e table (after 8 warm-ups). k=(*idum)/IQ1; *idum=IA1*(*idum-k*IQ1)-k*IR1; if (*idum < 0) *idum += IM1; if (j < NTAB) iv[j] = *idum; } iy=iv[0]; } k=(*idum)/IQ1; // Start here when not initializing. *idum=IA1*(*idum-k*IQ1)-k*IR1; // Compute idum=(IA1*idum) % IM1 without // over ows by Schrage's method. if (*idum < 0) *idum += IM1; k=idum2/IQ2; idum2=IA2*(idum2-k*IQ2)-k*IR2; // Compute idum2=(IA2*idum) % IM2 likewise. if (idum2 < 0) idum2 += IM2; j=iy/NDIV; // Will be in the range 0..NTAB-1. iy=iv[j]-idum2; // Here idum is shu.ed, idum and idum2 are // combined to generate output. iv[j] = *idum; if (iy < 1) iy += IMM1; if ((temp=AM*iy) > RNMX) return RNMX; // Because users don't expect endpoint values. else return temp; } void ByteSwap(unsigned char * b, int n) { register int i = 0; register int j = n-1; char temp; while (i5) return 1; if(nAstroPSolved==2){ inv[0]=0; inv[1]=1; } if(nAstroPSolved==3){ inv[0]=1; inv[1]=2; } if(nAstroPSolved==4){ inv[0]=0; inv[1]=1; inv[2]=2; inv[3]=3; } if(nAstroPSolved==5){ inv[0]=1; inv[1]=2; inv[2]=3; inv[3]=4; } return 0; } void writeBinFiles(double* systemMatrix,long* matrixIndex,int* instrCol,double* knownTerms,char* wrfileDir,char* wpath, struct comData comlsqr, int debugMode){ int nproc,myid; FILE *fpSM,*fpMI,*fpII,*fpKT,*fpBar; int nparam; MPI_Status status; char actpath[1024]; size_t sizePath=1020; int extConstraint; int barConstraint; MPI_Comm_size(MPI_COMM_WORLD, &nproc); MPI_Comm_rank(MPI_COMM_WORLD, &myid); nparam=comlsqr.nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrPSolved+comlsqr.nGlobP; getcwd(actpath,sizePath); chdir(wpath); if(!(chdir(wrfileDir)==0)) { printf("wrfile directory does not exist. Aborting\n"); MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } long int rowCounter=1; long int rowWritten=0; long int stillOpen=-1; int nOfFile=1; // file numerati da 1 in avanti. Il corrente file da scrivere è nOfFile int lastStargsrId=-1; rowWritten=0; rowCounter=0; long rowCounterPrev=0; extConstraint=comlsqr.extConstraint; barConstraint=comlsqr.barConstraint; if(myid>0){ MPI_Recv(&lastStargsrId, 1, MPI_INT, myid-1, 0, MPI_COMM_WORLD, &status); MPI_Recv(&rowCounterPrev, 1, MPI_LONG, myid-1, 1, MPI_COMM_WORLD, &status); //righe scritte nell'ultimo file MPI_Recv(&nOfFile, 1, MPI_INT, myid-1, 2, MPI_COMM_WORLD, &status); // printf("TP0 PE=%d ricevo da %d lastStargsrId=%d rowCounterPrev=%ld nOfFile=%ld\n",myid,myid-1,lastStargsrId,rowCounterPrev,nOfFile); } if(lastStargsrId==1000*nOfFile-1) // sono in fase di chiusura file { if(lastStargsrId!=matrixIndex[0]/comlsqr.nAstroPSolved ){ //rinomino il file precedente era alla fine, io sono in stella successiva! char filenameCoeffBinMatrixIndex[512]; strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); char tmp1[20],tmp2[20],tmp3[20]; sprintf(tmp1,"%09d",1000*(nOfFile-1)); sprintf(tmp2,"%09d",lastStargsrId); sprintf(tmp3,"%07ld",rowCounterPrev); rowCounterPrev=0; strcat(filenameCoeffBinMatrixIndex, tmp1); strcat(filenameCoeffBinMatrixIndex, "_"); strcat(filenameCoeffBinMatrixIndex, tmp2); strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-"); strcat(filenameCoeffBinMatrixIndex, tmp3); strcat(filenameCoeffBinMatrixIndex, "_MI.bin"); char filenameCoeffBinSystemMatrix[512]; strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinSystemMatrix, tmp1); strcat(filenameCoeffBinSystemMatrix, "_"); strcat(filenameCoeffBinSystemMatrix, tmp2); strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-"); strcat(filenameCoeffBinSystemMatrix, tmp3); strcat(filenameCoeffBinSystemMatrix, "_SM.bin"); char filenameCoeffBinIstrIndex[512]; strcpy(filenameCoeffBinIstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinIstrIndex, tmp1); strcat(filenameCoeffBinIstrIndex, "_"); strcat(filenameCoeffBinIstrIndex, tmp2); strcat(filenameCoeffBinIstrIndex, "_000000_nrows-"); strcat(filenameCoeffBinIstrIndex, tmp3); strcat(filenameCoeffBinIstrIndex, "_II.bin"); char filenameCoeffBinKnownTerms[512]; strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinKnownTerms, tmp1); strcat(filenameCoeffBinKnownTerms, "_"); strcat(filenameCoeffBinKnownTerms, tmp2); strcat(filenameCoeffBinKnownTerms, "_000000_nrows-"); strcat(filenameCoeffBinKnownTerms, tmp3); strcat(filenameCoeffBinKnownTerms, "_KT.bin"); stillOpen=0; rename("MI.bin", filenameCoeffBinMatrixIndex); rename("SM.bin", filenameCoeffBinSystemMatrix); rename("II.bin", filenameCoeffBinIstrIndex); rename("KT.bin", filenameCoeffBinKnownTerms); nOfFile++; // printf("TP3 PE=%d rinominati file %s %s %s %s\n",myid,filenameCoeffBinMatrixIndex, filenameCoeffBinSystemMatrix,filenameCoeffBinIstrIndex,filenameCoeffBinKnownTerms); } // if(lastStargsrId!=matrixIndex[0]/comlsqr.nAstroPSolved } // if(lastStargsrId==1000*nOfFile-1 fpMI=fopen("MI.bin","ab"); fpSM=fopen("SM.bin","ab"); fpII=fopen("II.bin","ab"); fpKT=fopen("KT.bin","ab"); // printf("TP4 PE=%d aperti in append files MI.bin ecc\n",myid); for(long int i=0;i0) fwrite(&instrCol[rowWritten*comlsqr.nInstrPSolved],sizeof(int),rowCounter*comlsqr.nInstrPSolved,fpII); fwrite(&knownTerms[rowWritten],sizeof(double),rowCounter,fpKT); fclose(fpMI); fclose(fpSM); fclose(fpII); fclose(fpKT); // printf("TP6 PE=%d i=%ld, rowCounter=%ld rowWritten=%ld rowCounterPrev=%ld\n",myid,i,rowCounter,rowWritten,rowCounterPrev); if(myid==nproc-1){ lastStargsrId=matrixIndex[comlsqr.mapNoss[myid]*2-comlsqr.multMI]/comlsqr.nAstroPSolved; // printf("TP7 PE=%d i=%ld lastStargsrId=%d\n",myid,i,lastStargsrId); char filenameCoeffBinMatrixIndex[512]; strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); char tmp1[20],tmp2[20],tmp3[20]; sprintf(tmp1,"%09d",1000*(nOfFile-1)); sprintf(tmp2,"%09d",lastStargsrId); sprintf(tmp3,"%07ld",rowCounter+rowCounterPrev); strcat(filenameCoeffBinMatrixIndex, tmp1); strcat(filenameCoeffBinMatrixIndex, "_"); strcat(filenameCoeffBinMatrixIndex, tmp2); strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-"); strcat(filenameCoeffBinMatrixIndex, tmp3); strcat(filenameCoeffBinMatrixIndex, "_MI.bin"); char filenameCoeffBinSystemMatrix[512]; strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinSystemMatrix, tmp1); strcat(filenameCoeffBinSystemMatrix, "_"); strcat(filenameCoeffBinSystemMatrix, tmp2); strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-"); strcat(filenameCoeffBinSystemMatrix, tmp3); strcat(filenameCoeffBinSystemMatrix, "_SM.bin"); char filenameCoeffBinIstrIndex[512]; strcpy(filenameCoeffBinIstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinIstrIndex, tmp1); strcat(filenameCoeffBinIstrIndex, "_"); strcat(filenameCoeffBinIstrIndex, tmp2); strcat(filenameCoeffBinIstrIndex, "_000000_nrows-"); strcat(filenameCoeffBinIstrIndex, tmp3); strcat(filenameCoeffBinIstrIndex, "_II.bin"); char filenameCoeffBinKnownTerms[512]; strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcat(filenameCoeffBinKnownTerms, tmp1); strcat(filenameCoeffBinKnownTerms, "_"); strcat(filenameCoeffBinKnownTerms, tmp2); strcat(filenameCoeffBinKnownTerms, "_000000_nrows-"); strcat(filenameCoeffBinKnownTerms, tmp3); strcat(filenameCoeffBinKnownTerms, "_KT.bin"); rename("MI.bin", filenameCoeffBinMatrixIndex); rename("SM.bin", filenameCoeffBinSystemMatrix); rename("II.bin", filenameCoeffBinIstrIndex); rename("KT.bin", filenameCoeffBinKnownTerms); rowCounterPrev=rowCounter; //inutile // printf("TP7.1 PE=%d i=%ld, rowCounter=%ld rowWritten=%ld rinominati i files MI.bin ecc FINITO!!\n",myid,i,rowCounter,rowWritten); break; } //if(myid==nproc-1) rowCounterPrev+=rowCounter; // printf("TP7.2 PE=%d i=%ld, rowCounter=%ld rowWritten=%ld rowCounterPrev=%ld PROX PE!!\n",myid,i,rowCounter,rowWritten,rowCounterPrev); break; //esce dal ciclo for }//if(i==comlsqr.mapNoss[myid]-1 ) if ( (matrixIndex[i*comlsqr.multMI]/comlsqr.nAstroPSolved) < 1000*nOfFile ) { rowCounter++; continue; } else { fwrite(&matrixIndex[rowWritten*comlsqr.multMI],sizeof(long int),rowCounter*comlsqr.multMI,fpMI); fclose(fpMI); for(int q=0;q0) fwrite(&instrCol[rowWritten*comlsqr.nInstrPSolved],sizeof(int),rowCounter*comlsqr.nInstrPSolved,fpII); fclose(fpII); fwrite(&knownTerms[rowWritten],sizeof(double),rowCounter,fpKT); fclose(fpKT); char filenameCoeffBinMatrixIndex[512]; char filenameCoeffBinSystemMatrix[512]; char filenameCoeffBinInstrIndex[512]; char filenameCoeffBinKnownTerms[512]; strcpy(filenameCoeffBinMatrixIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcpy(filenameCoeffBinSystemMatrix,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcpy(filenameCoeffBinInstrIndex,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); strcpy(filenameCoeffBinKnownTerms,"Gsr_dpccu3dpctavugsrgsrsystemrow_0000_"); char tmp1[20],tmp2[20],tmp3[20]; sprintf(tmp1,"%09d",1000*(nOfFile-1)); sprintf(tmp2,"%09d",1000*nOfFile-1); sprintf(tmp3,"%07ld",rowCounter+rowCounterPrev); strcat(filenameCoeffBinMatrixIndex, tmp1); strcat(filenameCoeffBinMatrixIndex, "_"); strcat(filenameCoeffBinMatrixIndex, tmp2); strcat(filenameCoeffBinMatrixIndex, "_000000_nrows-"); strcat(filenameCoeffBinMatrixIndex, tmp3); strcat(filenameCoeffBinMatrixIndex, "_MI.bin"); strcat(filenameCoeffBinSystemMatrix, tmp1); strcat(filenameCoeffBinSystemMatrix, "_"); strcat(filenameCoeffBinSystemMatrix, tmp2); strcat(filenameCoeffBinSystemMatrix, "_000000_nrows-"); strcat(filenameCoeffBinSystemMatrix, tmp3); strcat(filenameCoeffBinSystemMatrix, "_SM.bin"); strcat(filenameCoeffBinInstrIndex, tmp1); strcat(filenameCoeffBinInstrIndex, "_"); strcat(filenameCoeffBinInstrIndex, tmp2); strcat(filenameCoeffBinInstrIndex, "_000000_nrows-"); strcat(filenameCoeffBinInstrIndex, tmp3); strcat(filenameCoeffBinInstrIndex, "_II.bin"); strcat(filenameCoeffBinKnownTerms, tmp1); strcat(filenameCoeffBinKnownTerms, "_"); strcat(filenameCoeffBinKnownTerms, tmp2); strcat(filenameCoeffBinKnownTerms, "_000000_nrows-"); strcat(filenameCoeffBinKnownTerms, tmp3); strcat(filenameCoeffBinKnownTerms, "_KT.bin"); rename("MI.bin", filenameCoeffBinMatrixIndex); rename("SM.bin", filenameCoeffBinSystemMatrix); rename("II.bin", filenameCoeffBinInstrIndex); rename("KT.bin", filenameCoeffBinKnownTerms); // printf("TP9 PE=%d i=%ld Rinominati sui file GSR_... MI.bin ecc %s \n",myid,i,filenameCoeffBinMatrixIndex ); fpMI=fopen("MI.bin","wb"); fpSM=fopen("SM.bin","wb"); fpII=fopen("II.bin","wb"); fpKT=fopen("KT.bin","wb"); rowWritten+=rowCounter; rowCounter=1; // sono già sulla successiva stella // printf("TP10 PE=%d aperti in wb files MI.bin ecc\n",myid); rowCounterPrev=0; nOfFile++; }//else //Scrivere su file } if(myidnumOfExtStarinFile+nwr)numOfExtStartoWrite=numOfExtStarinFile; else numOfExtStartoWrite=nStar-nwr; startStar=nwr; endStar=startStar+numOfExtStartoWrite-1; double *buffArray; buffArray=(double *)calloc(numOfExtStartoWrite*nAstroPSolved*nEqExtConstr,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.; buffArray[localCounter]=randVal; localCounter++; } char filenameExtConstrAstro[512]; char tmp1[20],tmp2[20],tmp3[20]; sprintf(tmp1,"%09d",startStar); sprintf(tmp2,"%09d",endStar); sprintf(tmp3,"%09d",nEqExtConstr); strcpy(filenameExtConstrAstro,"Gsr_nullspaceastrofit_0000_"); strcat(filenameExtConstrAstro, tmp1); strcat(filenameExtConstrAstro, "_"); strcat(filenameExtConstrAstro, tmp2); strcat(filenameExtConstrAstro, "_000000_nrows-"); strcat(filenameExtConstrAstro, tmp3); strcat(filenameExtConstrAstro,".bin"); fpSM=fopen(filenameExtConstrAstro,"wb"); fwrite(buffArray,sizeof(double),nEqExtConstr*numOfExtStartoWrite*nAstroPSolved,fpSM); fclose(fpSM); free(buffArray); } // for(in nwr=.. /// write extConstr for Att double *buffArray; buffArray=(double *)calloc(nDegFreedomAtt,sizeof(double)); for(int i=0;inumOfBarStarinFile+nwr)numOfBarStartoWrite=numOfBarStarinFile; else numOfBarStartoWrite=nStar-nwr; startStar=nwr; endStar=startStar+numOfBarStartoWrite-1; double *buffArray; buffArray=(double *)calloc(numOfBarStartoWrite*nAstroPSolved*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.; buffArray[localCounter]=randVal; localCounter++; } char filenameBarConstrAstro[512]; char tmp1[20],tmp2[20],tmp3[20]; sprintf(tmp1,"%09d",startStar); sprintf(tmp2,"%09d",endStar); sprintf(tmp3,"%09d",nEqBarConstr); strcpy(filenameBarConstrAstro,"Gsr_barconstrastrofit_0000_"); strcat(filenameBarConstrAstro, tmp1); strcat(filenameBarConstrAstro, "_"); strcat(filenameBarConstrAstro, tmp2); strcat(filenameBarConstrAstro, "_000000_nrows-"); strcat(filenameBarConstrAstro, tmp3); strcat(filenameBarConstrAstro,".bin"); fpBar=fopen(filenameBarConstrAstro,"wb"); fwrite(buffArray,sizeof(double),nEqBarConstr*numOfBarStartoWrite*nAstroPSolved,fpBar); fclose(fpBar); free(buffArray); } // for(in nwr=.. }//if (barConstraint MPI_Barrier(MPI_COMM_WORLD); chdir(wpath); chdir(actpath); } int cmpfunc (const void * a, const void * b) { return ( *(int*)a - *(int*)b ); } int randint(int max) { return (int) (max*(rand()/(RAND_MAX+1.0))); } // This function generates a psuedo-random integer n, where min <= n <= max int randint1(int min, int max) { if (min>max) { return max+(int)((min-max+1)*(rand()/(RAND_MAX+1.0))); } else { return min+(int)((max-min+1)*(rand()/(RAND_MAX+1.0))); } } // This function generates a psuedo-random long integer n, where 0 <= n < max long randlong(long max) { return (long) (max*(rand()/(RAND_MAX+1.0))); } // This function generates a psuedo-random long integer n, where min <= n <= max long randlong1(long min, long max) { if (min>max) { return max+(long)((min-max+1)*(rand()/(RAND_MAX+1.0))); } else { return min+(long)((max-min+1)*(rand()/(RAND_MAX+1.0))); } } // This function computes the instrument hash given the number of FoV, CCD, PixelColumn and TimeInterval. // Since these numbers must fill in 1, 8, 11, and 11 bits respectively, they are checked against appropriate // maximum values, and if their values are greater, the function returns -1. // If the function is called with FoV, CCD, PixelColumn and TimeInterval it returns the instrId hash, // if it is called with nFoVs, nCCDs, nPixelColumns and nTimeIntervals it returns the instrSetUp hash. // NB: the FoV is treated differently since it can be 1 or 2, so what is actually hashed is nFoVs-1 to keep it into 1 bit. long instr_hash(int FoV, int CCD, int PixelColumn, int TimeInterval) { short CCDOffset = 1; short PixelOffset = 9; short TimeOffset = 20; short FoVMaxValue = ((short) pow(2.0, CCDOffset)); short CCDMaxValue = ((short) pow(2.0, PixelOffset-CCDOffset)); short PixelMaxValue = ((short) pow(2.0, TimeOffset-PixelOffset)); short TimeMaxValue = ((short) pow(2.0, TimeOffset-PixelOffset)); if(FoV>FoVMaxValue || CCD>CCDMaxValue || PixelColumn>PixelMaxValue || TimeInterval>TimeMaxValue) return -1; // Some parameter is greater than its maximum possible value. The hash cannot be calculated. else return (((long) FoV-1) | (((long) CCD) << (CCDOffset)) | (((long) PixelColumn) << (PixelOffset)) | (((long) TimeInterval) << (TimeOffset))); } // This function extracts an integer in the range pos_min..pos_max which points to a position of the array . // This implies that values is an array with pos_max+1 positions. void swap(long *i, long *j) { long t = *i; *i = *j; *j = t; } int fill_extract(long *values, long *pos_min, long pos_max, long *number) { long position; if(*pos_min==pos_max-1) position=*pos_min; else position=randlong1(*pos_min, pos_max-1); *number=values[position]; swap(&values[*pos_min], &values[position]); (*pos_min)++; if(*pos_min>=pos_max) return 1; return 0; } struct nullSpace cknullSpace(double * systemMatrix,long * matrixIndex,double *attNS,struct comData comlsqr){ struct nullSpace results; int nproc, myid; long nunkSplitNS; double * nullSpaceVect; double * prodNS; long nElements, nStar; int nEqExtConstr; int nparam; int nLocalStar; int firstStarConstr,lastStarConstr; int nOfElextObs,numOfExtStar,numOfExtAttCol,startingAttColExtConstr; short nAstroPSolved,nAttParAxis,nAttAxes; double *nullSpaceFPN; double *productNS; int npeSend,npeRecv; double sum,extConstrW; long sumVer; int setBound[4]; long int l1,k,l,j; long int nDegFreedomAtt,localAstroMax,offsetAttParam; long int *mapNoss, *mapNcoeff; time_t seconds[2], tot_sec; int **mapStar; MPI_Status status; MPI_Request req1; MPI_Comm_size(MPI_COMM_WORLD, &nproc); MPI_Comm_rank(MPI_COMM_WORLD, &myid); nEqExtConstr=comlsqr.nEqExtConstr; firstStarConstr=comlsqr.firstStarConstr; lastStarConstr=comlsqr.lastStarConstr; nAstroPSolved=comlsqr.nAstroPSolved; nOfElextObs=comlsqr.nOfElextObs; nDegFreedomAtt=comlsqr.nDegFreedomAtt; nAttParAxis=comlsqr.nAttParAxis; localAstroMax=comlsqr.VrIdAstroPDimMax*nAstroPSolved; offsetAttParam=comlsqr.offsetAttParam; numOfExtStar=comlsqr.numOfExtStar; nAttAxes=comlsqr.nAttAxes; mapNcoeff=comlsqr.mapNcoeff; mapNoss=comlsqr.mapNoss; numOfExtAttCol=comlsqr.numOfExtAttCol; startingAttColExtConstr=comlsqr.startingAttColExtConstr; setBound[0]=comlsqr.setBound[0]; setBound[1]=comlsqr.setBound[1]; setBound[2]=comlsqr.setBound[2]; setBound[3]=comlsqr.setBound[3]; extConstrW=comlsqr.extConstrW; int multMI=comlsqr.multMI; mapStar=comlsqr.mapStar; nStar=comlsqr.nStar; int nAttParam=comlsqr.nAttParam; int nInstrParam=comlsqr.nInstrParam; int nGlobalParam=comlsqr.nGlobalParam; // errore ==> nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrP+comlsqr.nGlobP; nparam=nAstroPSolved+comlsqr.nAttP+comlsqr.nInstrPSolved+comlsqr.nGlobP; nLocalStar=comlsqr.mapStar[myid][1]-comlsqr.mapStar[myid][0]+1; nElements = mapNoss[myid]+nEqExtConstr; nullSpaceFPN = (double *) calloc(nAstroPSolved*nEqExtConstr, sizeof(double)); if (!nullSpaceFPN) exit(err_malloc("nullSpaceFPN",myid)); productNS = (double *) calloc(nElements, sizeof(double)); if (!productNS) exit(err_malloc("productNS",myid)); npeSend=myid-1; npeRecv=myid+1; for(int i=0; i0){ MPI_Isend(&systemMatrix[mapNoss[myid]*nparam+i*nOfElextObs],nAstroPSolved, MPI_DOUBLE, npeSend, 1,MPI_COMM_WORLD, &req1); } if(myid0) MPI_Wait(&req1,&status); MPI_Barrier(MPI_COMM_WORLD); } prodNS = (double *) calloc(nElements, sizeof(double)); if (!prodNS) exit(err_malloc("prodNS",myid)); nunkSplitNS=localAstroMax + nAttParam+nInstrParam+nGlobalParam; nullSpaceVect = (double *) calloc(nunkSplitNS, sizeof(double)); if (!nullSpaceVect) exit(err_malloc("nullSpaceVect",myid)); for(int ic=0; iclastStarConstr){ for(int j=0;j=mapNoss[myid]*multMI){ printf("ERROR: PE=%d i=%d lset=%d multMI*i=%d greather than mapNoss[myid]*multMI=%ld\n",myid,i,lset,multMI*i,mapNoss[myid]*multMI); continue; }*/ long numOfStarPos=matrixIndex[multMI*i]/nAstroPSolved; j=(numOfStarPos-comlsqr.mapStar[myid][0])*nAstroPSolved; } else j++; } if(lset>=setBound[1] && lset=mapNoss[myid]*multMI){ printf("ERROR: PE=%d i=%d lset=%d multMI*i+(multMI-1)=%d greather than mapNoss[myid]*multMI=%ld\n",myid,i,lset,multMI*i+(multMI-1),mapNoss[myid]*multMI); continue; }*/ j=matrixIndex[multMI*i+(multMI-1)]+((lset-setBound[1])/nAttParAxis)*nDegFreedomAtt+(localAstroMax-offsetAttParam); } else j++; } sum=sum+systemMatrix[l]*nullSpaceVect[j]; double NSVal; if(lset=setBound[1]+4 && lset=setBound[1]+8) NSVal=(matrixIndex[multMI*i+1]-nStar*nAstroPSolved)+ic*nDegFreedomAtt+lset-(setBound[1]+8); else NSVal=0; } sumVer=sumVer+1.0*NSVal; } lset++; }//for(l /* if(i>=nElements){ printf("ERROR: PE=%d i=%d greather than nElements=%ld\n",myid,i ,nElements); continue; }*/ if(sumVer != chkSumVer){ printf("ERROR: PE=%d NullSapce Equation ic=%d, sumVer[%d]=%d and chkSumVer=%ld are not equal\n",myid,ic,i,sumVer, chkSumVer); MPI_Abort(MPI_COMM_WORLD, 0); exit(1); } productNS[i]=sum; }//for i //////////// in case of ExtConstr /* TO BE DECIDED orthogonalty of nullspace base vectors for(i2=0;i2productNS[j]) localMin=productNS[j]; if(localMax=bins){ printf("Error while computing local myid=0 ix=%d on cknullSpace at element %d on equation %d. prodNS=%lf greater than max=%lf\n",ix,j,ic,productNS[j],results.compMax[ic]); MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } isto[ix]++; } for(int k=1;k=bins){ printf("Error while computing remote peid=%d ix=%d on cknullSpace at element %d on equation %d. prodNS=%lf greater than max=%lf\n",k,ix,j,ic,productNS[j],results.compMax[ic]); MPI_Abort(MPI_COMM_WORLD,1); exit(EXIT_FAILURE); } isto[ix]++; } } fprintf(fpNS, "%le\n%le\n%le\n%le\n%le\n%d\n%le\n",results.compMin[ic],results.compMax[ic],results.compAvg[ic],results.compVar[ic], results.vectNorm[ic],bins,binWidth); for(int j=0;j input_len) { return NULL; } strncpy (dest, input + offset, len); return dest; } /*--------------------------------------------------------------------------*/ // Compute the value of the shifted Legendre polinomial of degree deg. Here we // just need deg<3. double legendre(int deg, double x) { double res; switch(deg) { case 0: res = 1.0; break; case 1: res = 2.0*(x-0.5); break; case 2: res = 6.0*(x-0.5)*(x-0.5)-0.5; break; default: res = -1.0; break; } return res; } //void computeInstrContr (int lsInstrFlag, int ssInstrFlag, long nFoVs,long nCCDs, long nPixelColumns,long nTimeIntervals,double * instrCoeffConstr,int * instrColsConstr,int * instrConstrIlung) int computeInstrConstr (struct comData comlsqr,double * instrCoeffConstr,int * instrColsConstr,int * instrConstrIlung) { ////////////////// Writing instrConstrRows_xxxx.bin file // 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. // MUST BE ALLOCATED THE FOLLOWING VECTORS: //instrCoeffConstr=(double *) calloc(nElemIC, sizeof(double)); //instrColsConstr=(int *) calloc(nElemIC, sizeof(int)); //instrConstrIlung=(int *) calloc(nOfInstrConstr, sizeof(int)); int lsInstrFlag=comlsqr.lsInstrFlag; int ssInstrFlag=comlsqr.ssInstrFlag; long nFoVs=1+comlsqr.instrConst[0]; long nCCDs=comlsqr.instrConst[1]; long nPixelColumns=comlsqr.instrConst[2]; long nTimeIntervals=comlsqr.instrConst[3]; long cCDLSAACZP=comlsqr.cCDLSAACZP; int nElemICLSAL =comlsqr.nElemICLSAL; int nElemICLSAC =comlsqr.nElemICLSAC; int nElemICSS = comlsqr.nElemICSS; int nOfInstrConstr = comlsqr.nOfInstrConstr; int nElemIC = comlsqr.nElemIC; int counterElem=0; int counterEqs=0; int elemAcc=0; if(lsInstrFlag){ // generate large scale constraint eq. AL for(int i=0; i