util.c

#include <stdlib.h>
#include <stdio.h>
#include "util.h"
#include <mpi.h>
#include <math.h>
#ifdef OMP
#include <omp.h>
#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<totrows;jj++){
    
    if(instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)]==-1){ //is a Constraint
		for(int kk=0;kk<nInstrPSolved;kk++)
            instrCols[jj*nInstrPSolved+kk]=0;
        continue;
    }
	// FoV = instrIndexPointer[jj*DEFAULT_NINSTRINDEXES], CCD = instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]
	// PixelColumn = instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+2], TimeInterval = instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+3]
	int FoV = instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)];
	int CCD = instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)+1];
    int PixCol = instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)+2];
    int TInt = instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)+3];
    int ACALFlag = instrIndexPointer[jj*(DEFAULT_NINSTRINDEXES+1)+4];

	int counter=0;
		if(maInstrFlag) {
			instrCols[jj*nInstrPSolved+counter] = CCD-1; // Index_CMag = CCD-1
			counter++;
		}
		if(nuInstrFlag) {
			// Index_Cnu = offsetCMag + (FoV−1)*nCCDs + (CCD−1)
			instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCMag + FoV*nCCDs + (CCD-1);
			counter++;
		}
		if(ssInstrFlag) {
			if(ACALFlag) {
				// Index Cdelta_zeta = offsetCDelta_eta_3 + (CCD−1)*nPixelColumns + (PixelColumn−1)
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCDelta_eta_3 + (CCD-1)*nPixCols + (PixCol-1);
				counter++;
			} else {
				// Index Cdelta_eta = offsetCnu + (CCD−1)*nPixelColumns + (PixelColumn−1)
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCnu + (CCD-1)*nPixCols + (PixCol-1);
				counter++;
			}
		}
	// For performance reasons, compute the relative Delta_eta index only once
    // relPos_Delta_eta   = (FoV-1)*nCCDs*nTimeIntervals+(CCD-1)*nTimeIntervals+(TimeInterval-1)
    //                    = (instrIndexPointer[0]-1)*instrConst[1]*instrConst[3]+(instrIndexPointer[1]-1)*instrConst[3]+(instrIndexPointer[3]-1)
		if(lsInstrFlag) {
			relPos_ls = FoV*nCCDs*nTInts+(CCD-1)*nTInts+(TInt-1);
			if(ACALFlag) {
				// Index CDelta_zeta_1 = offsetCdelta_zeta + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCdelta_zeta + relPos_ls;
				counter++;
				// Index CDelta_zeta_2 = offsetCDelta_zeta_1 + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCDelta_zeta_1 + relPos_ls;
				counter++;
				// Index CDelta_zeta_3 = offsetCDelta_zeta_2 + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCDelta_zeta_2 + relPos_ls;
				counter++;
			} else {
				// Index CDelta_eta_1 = offsetCdelta_eta + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCdelta_eta + relPos_ls;
				counter++;
				// Index CDelta_eta_2 = offsetCDelta_eta_1 + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCDelta_eta_1 + relPos_ls;
				counter++;
				// Index CDelta_eta_3 = offsetCDelta_eta_2 + relPos_Delta_zeta
				instrCols[jj*nInstrPSolved+counter] = comlsqr.offsetCDelta_eta_2 + relPos_ls;
				counter++;
			}
		}
//	}
	if(counter!=nInstrPSolved) {
		printf("SEVERE ERROR PE=%d counter=%d != nInstrPSolved=%d on row #%d\n",comlsqr.myid, counter, nInstrPSolved, jj);
		MPI_Abort(MPI_COMM_WORLD, 1);
		exit(EXIT_FAILURE);

	}
    for(int k=0;k<nInstrPSolved;k++){
        if(instrCols[jj*nInstrPSolved+k]>=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;jj<totrows;jj++){
        testConstr=0;
         for(int k=0;k<6;k++){
            testConstr+=instrCols[jj*DEFAULT_NINSTRVALUES+k];
        }
        if(testConstr==0){ //is a Constraint
            instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+0]=0;
            instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]=0;
            instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+2]=0;
            instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+3]=0;
            continue;
        }
        // Index_CMag = CCD-1
        instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]=instrCols[jj*DEFAULT_NINSTRVALUES+0]+1;
        // Index_Cnu = offsetCMag + (FoV−1)*nCCDs + (CCD−1)
        instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+0]=(instrCols[jj*DEFAULT_NINSTRVALUES+1]-instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]+1-comlsqr.offsetCMag)/instrConst[1] + 1;
        // Index Cdelta_eta = offsetCnu + (CCD−1)*nPixelColumns + (PixelColumn−1)
        instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+2]=instrCols[jj*DEFAULT_NINSTRVALUES+2]-comlsqr.offsetCnu-(instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]-1)*instrConst[2]+1;
        
        relPos_Delta_eta=instrCols[jj*DEFAULT_NINSTRVALUES+3]-comlsqr.offsetCdelta_eta;
        instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+3]=relPos_Delta_eta -(instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+0]-1)*instrConst[1]*instrConst[3]-(instrIndexPointer[jj*DEFAULT_NINSTRINDEXES+1]-1)*instrConst[3]+1;
    }
    return;
}

/*--------------------------------------------------------------------------*/
void printerror(int status) {
	/*****************************************************/
	/* Print out cfitsio error messages and exit program */
	/*****************************************************/

	if (status) {
//		fits_report_error(stderr, status); /* print error report */
		MPI_Abort(MPI_COMM_WORLD, status);
		exit(status); /* terminate the program, returning error status */
	}
	return;
}

/*--------------------------------------------------------------------------*/
void printerrorsingle(int status) {
	/*****************************************************/
	/* Print out cfitsio error messages and exit program */
	/*****************************************************/

	if (status) {
//		fits_report_error(stderr, status); /* print error report */
		exit(status); /* terminate the program, returning error status */
	}
	return;
}

int err_malloc(const char *s,int id) {
	printf("out of memory while allocating %s on PE=%d.\n", s, id);
	MPI_Abort(MPI_COMM_WORLD, 1);
	return 1;
}

int sel(const struct dirent *a) {
	return ((strncmp(a->d_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;i<comlsqr.nAstroPSolved;i++)
		tempSendBuf[i]=vectToSum[i];
	for(int i=0;i<comlsqr.nAstroPSolved;i++)
		tempSendBuf[i+comlsqr.nAstroPSolved]=vectToSum[(comlsqr.VrIdAstroPDim-1)*comlsqr.nAstroPSolved+i];

	for(int i=0;i<nMov;i++)
	{
		if(i==0) //forward propagation!
		{
			npeSend=rank+1;
			if(npeSend==size) npeSend=0;
			npeRecv=rank-1;
			if(npeRecv<0) npeRecv=size-1;
		}
		if(i==1) //backward propagation!
		{
			npeSend=rank-1;
			if(npeSend<0) npeSend=size-1;
			npeRecv=rank+1;
			if(npeRecv==size) npeRecv=0;
		}
		MPI_Isend(tempSendIdBuf, comlsqr.multMI, MPI_INT, npeSend, 1,MPI_COMM_WORLD, &req2);
		MPI_Isend(tempSendBuf, comlsqr.multMI*comlsqr.nAstroPSolved, MPI_DOUBLE, npeSend, 2,
              MPI_COMM_WORLD, &req3);

		MPI_Recv(tempRecvIdBuf, comlsqr.multMI, MPI_INT, npeRecv, 1,MPI_COMM_WORLD, &status);
		MPI_Recv(tempRecvBuf, comlsqr.multMI*comlsqr.nAstroPSolved, MPI_DOUBLE, npeRecv, 2,
              MPI_COMM_WORLD, &status);

		MPI_Wait(&req2,&status);
		MPI_Wait(&req3,&status);

		MPI_Barrier(MPI_COMM_WORLD);					
		
		
		int okupd=0;
		if(tempRecvIdBuf[1]==comlsqr.mapStar[rank][0])
		{
		   for(int ns=0;ns<comlsqr.nAstroPSolved;ns++)
					vectToSum[ns]+=tempRecvBuf[comlsqr.nAstroPSolved+ns];
		   okupd=1;
		}
		if(tempRecvIdBuf[1]==comlsqr.mapStar[rank][1] && okupd==0)
		{
		   for(int ns=0;ns<comlsqr.nAstroPSolved;ns++)
					vectToSum[(comlsqr.VrIdAstroPDim-1)*comlsqr.nAstroPSolved+ns]+=tempRecvBuf[comlsqr.nAstroPSolved+ns];
		}
		
	    okupd=0;
		
		if(tempRecvIdBuf[0]!=tempRecvIdBuf[1])
		{
		
		  if(tempRecvIdBuf[0]==comlsqr.mapStar[rank][1] )
 		  {
			for(int ns=0;ns<comlsqr.nAstroPSolved;ns++)
				vectToSum[(comlsqr.VrIdAstroPDim-1)*comlsqr.nAstroPSolved+ns]+=tempRecvBuf[ns];
			okupd=1;
		}
		if(tempRecvIdBuf[0]==comlsqr.mapStar[rank][0])
		{
		   for(int ns=0;ns<comlsqr.nAstroPSolved;ns++)
					vectToSum[ns]+=tempRecvBuf[ns];
		}
		
		}//iftempRecvIdBuf[0]!=tempRecvIdBuf[1]
		
		
		
	} // next for
		
		

		free(tempSendBuf); 
		free(tempRecvBuf);
	
}

void initThread(long int  *matrixIndex,struct comData *comlsqr)
{
int myid=comlsqr->myid;

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;n<nthreads;n++)
	comlsqr->mapForThread[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;n<nthreads;n++)
{
	comlsqr->mapForThread[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;n<nthreads;n++)
    {

		while(matrixIndex[2*(comlsqr->mapForThread[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;n<nthreads;n++)
        {
            comlsqr->mapForThread[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;i<comlsqr->nSubsetAtt;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;nsubBlock<nthreads;nsubBlock++) //nel sottointervallo cerco l'indipendenza degli indici nelle nthreads del sistema a partire dal blocco della prima tread che confronto con le successive e poi la seconda thread che confronto con la terza e seguenti  ecc. ecc.
	{
	for(long  j=totalEleInBlockThread*i;j<totalEleInBlockThread*(i+1)+1;j++) //j spazzola tutti gli elementi nel blocco della thread "nsubBlok" per poi confrontarlo con tutti gli elelementi delle thread seguenti  
	{
		int indexFound=matrixIndex[j*2+1];
		for(long k=totalEleInBlockThread*(nsubBlock+1)+1;k<totalEleInBlock;k++) //k spazzola a partire dal primo elemento della thread seguente (nsubBlock) fino a tutti gli elementi del sottointervallo 
		{
			int kindexFound=matrixIndex[k*2+1];
			if(!(indexFound<=kindexFound-comlsqr->nAttParAxis || 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<comlsqr.mapNoss[myid];i++){

        counterAxis=0;
        counterInstr=0;


        for(ii=0;ii<nparam;ii++){
            if(ii>=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;i<nEqExtConstr;i++){
            for(int ns=0;ns<nAstroPSolved*numOfExtStar;ns++){
                systemMatrix[l]=systemMatrix[l]*preCondVect[ns];
                l++;
            }
            for(int naxis=0;naxis<nAttAxes;naxis++){
                for(int j=0;j<numOfExtAttCol;j++){
                    int ncolumn = VrIdAstroPDimMax*nAstroPSolved+startingAttColExtConstr+j+naxis*nDegFredoomAtt;
                    systemMatrix[l]=systemMatrix[l]*preCondVect[ncolumn];
                    l++;
                }
            }

        }
    }

    if(barConstraint){
        for(int i=0;i<nEqBarConstr;i++){
            for(int ns=0;ns<nAstroPSolved*numOfBarStar;ns++){
                systemMatrix[l]=systemMatrix[l]*preCondVect[ns];
                l++;
            }
            
        }
    }

    if(nElemIC>0){
        for(int i=0;i<nElemIC;i++){
            int ncolumn=offsetInstrParam+(VroffsetAttParam-offsetAttParam)+instrCol[mapNoss[myid]*nInstrPSolved+i];
            systemMatrix[l]=systemMatrix[l]*preCondVect[ncolumn];
            l++;
        }
    }

}    
void mpi_allreduce(double *source, double *dest,  long int lcount, 
                   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
{
	int ncyc=0;
	int count=__MSGSIZ_MAX;
	int chunch=__MSGSIZ_MAX;
	while (lcount>0)
	{
		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.
   */
{