util.c

    if(ssInstrFlag){
        // generate small scale constraint eq. AL
		double x;
		for(int i=0; i<nCCDs; i++) {
            for(int j=0; j<3; j++) { // each CCD generates 3 constraint equations, one for each order of the legendre polinomials
                instrConstrIlung[counterEqs] = nElemICSS;
                elemAcc+=nElemICSS;
                counterEqs++;
                for(int k=0; k<nPixelColumns; k++) {
					x=(k+0.5)/nPixelColumns;
                    instrCoeffConstr[counterElem] = legendre(j,x);
					if(instrCoeffConstr[counterElem]==-1.0) {
						printf("Error from legendre function when i=%d, j=%d, k=%d\n", i, j, k);
						return 0;
					}
                    instrColsConstr[counterElem] = comlsqr.offsetCnu + i*nPixelColumns + k;
                    counterElem++;
                }
            }
        }
        // generate small scale constraint eq. AC
        for(int i=0; i<nCCDs; i++) {
            for(int j=0; j<3; j++) { // each CCD generates 3 constraint equations, one for each order of the legendre polinomials
                instrConstrIlung[counterEqs] = nElemICSS;
                elemAcc+=nElemICSS;
                counterEqs++;
                for(int k=0; k<nPixelColumns; k++) {
					x=(k+0.5)/nPixelColumns;
                    instrCoeffConstr[counterElem] = legendre(j,x);
					if(instrCoeffConstr[counterElem]==-1.0) {
						printf("Error from legendre function when i=%d, j=%d, k=%d\n", i, j, k);
						return 0;
					}
                    instrColsConstr[counterElem] = comlsqr.offsetCDelta_eta_3 + i*nPixelColumns + k;
                    counterElem++;
                }
            }
        }
    }
}
if(counterEqs!=nOfInstrConstr) {
    printf("SEVERE ERROR  counterEqs =%d does not coincide with nOfInstrConstr=%d\n", counterEqs, nOfInstrConstr);
    return 0;
}
if(counterElem!=nElemIC) {
    printf("SEVERE ERROR  counterElem =%d does not coincide with nElemIC=%d\n", counterElem, nElemIC);
    return 0;
}
if(elemAcc!=nElemIC) {
    printf("SEVERE ERROR   elemAcc =%d does not coincide with nElemIC=%d\n", elemAcc, nElemIC);
    return 0;
}
    return 1;
}
void swapInstrCoeff(double * instrCoeff, long repeat, long nrows){
    double * tmp;
    tmp = (double *) calloc(repeat * nrows , sizeof(double));
    for (int j=0;j<nrows;j++){
        tmp[j*repeat+0]=instrCoeff[j*repeat+4];
        tmp[j*repeat+1]=instrCoeff[j*repeat+5];
        tmp[j*repeat+2]=instrCoeff[j*repeat+0];
        tmp[j*repeat+3]=instrCoeff[j*repeat+1];
        tmp[j*repeat+4]=instrCoeff[j*repeat+2];
        tmp[j*repeat+5]=instrCoeff[j*repeat+3];
    }
    for (int j=0;j<nrows*repeat;j++)
        instrCoeff[j]=tmp[j];
    free(tmp);
    return;
}
float simfullram(long &nStar, long &nobs, float memGlobal, int nparam, int nAttParam, int nInstrParam){
    float smGB=0., ktGB=0., miGB=0.,iiGB=0.,auxGB=0., memGB=0., prevmemGB;
    long prevnStar, prevnobs;
    long gigaByte=1024*1024*1024;
    long ncoeff;

        
        ncoeff = nparam * nobs; // total number of non-zero coefficients of the system
        smGB=(float)(ncoeff)*8/(gigaByte);  //systemMatrix
        ktGB=(float)(nobs)*8/(gigaByte);     //knownTerms
        miGB=(float)(nobs*2)*8/(gigaByte);   //matrixIndex
        iiGB=(float)(nobs*6)*4/(gigaByte);   //InstrIndex
        auxGB=(float)(nStar*5+nAttParam+nInstrParam+0)*8/(gigaByte); //precondVect+vVect+wVect+xSolution+standardError
        memGB=smGB+miGB+ktGB+iiGB+5*auxGB;
        if(memGlobal < memGB){
            return memGlobal;
        }
        
        while(memGB < memGlobal){
            prevnStar=nStar;
            prevnobs=nobs;
            prevmemGB=memGB;
            nStar*=2;
            nobs*=3;
            ncoeff = nparam * nobs;
            smGB=(float)(ncoeff)*8/(gigaByte);  //systemMatrix
            ktGB=(float)(nobs)*8/(gigaByte);     //knownTerms
            miGB=(float)(nobs*2)*8/(gigaByte);   //matrixIndex
            iiGB=(float)(nobs*6)*4/(gigaByte);   //InstrIndex
            auxGB=(float)(nStar*5+nAttParam+nInstrParam+0)*8/(gigaByte); //precondVect+vVect+wVect+xSolution+standardError
            memGB=smGB+miGB+ktGB+iiGB+5*auxGB;
        }
    nStar=prevnStar;
    nobs=prevnobs;
    memGB=prevmemGB;
    while(memGB < memGlobal){
        prevnStar=nStar;
        prevnobs=nobs;
        prevmemGB=memGB;
        nobs+=10000;
        ncoeff = nparam * nobs;
        smGB=(float)(ncoeff)*8/(gigaByte);  //systemMatrix
        ktGB=(float)(nobs)*8/(gigaByte);     //knownTerms
        miGB=(float)(nobs*2)*8/(gigaByte);   //matrixIndex
        iiGB=(float)(nobs*6)*4/(gigaByte);   //InstrIndex
        auxGB=(float)(nStar*5+nAttParam+nInstrParam+0)*8/(gigaByte); //precondVect+vVect+wVect+xSolution+standardError
        memGB=smGB+miGB+ktGB+iiGB+5*auxGB;
    }
    nStar=prevnStar;
    nobs=prevnobs;
    memGB=prevmemGB;

    return prevmemGB;
}