Commit 5cdccaeb authored by ofgeorg's avatar ofgeorg

added BBOB 2010

parent 1c023ac3
......@@ -35,10 +35,6 @@ can be passed around
struct paramStruct{
/* all chars have a fixed length, to ease memory management */
/* row ou col (defaut) - format of the design variable data when in arrays of vectors
We might forbid the "row" format in C ... */
char inputFormat[4];
/* a string for the name of the optimiser (used in the post-processing)
If not specified, PARAMS.algName = 'not-specified'; */
char algName[DefaultStringLength];
......@@ -101,8 +97,8 @@ typedef struct paramStruct ParamStruct;
BestFEval
**********************************/
struct lastEvalStruct {
/* Number of evaluations */
unsigned int num;
/* Number of evaluations : double to achieve larger precision */
double num;
/* fitness value */
double F;
......@@ -127,13 +123,14 @@ typedef struct lastEvalStruct LastEvalStruct;
***********************************/
double fgeneric_initialize(ParamStruct PARAMS);
double fgeneric_finalize();
double fgeneric_finalize();
double fgeneric_evaluate(double * X);
double fgeneric_ftarget();
unsigned int fgeneric_maxevals(unsigned int DIM);
unsigned int fgeneric_evaluations();
double fgeneric_maxevals(unsigned int DIM);
double fgeneric_evaluations();
double fgeneric_best();
void fgeneric_noiseseed(unsigned int seed);
int fgeneric_exist(unsigned int FUNC_ID);
void fgeneric_noiseseed(unsigned long seed);
ParamStruct fgeneric_getDefaultPARAMS();
#endif
......@@ -29,8 +29,8 @@ Loosely inspired by fgeneric.m, the matlab version
#include <float.h>
#include <string.h>
//int strcasecmp (const char *s1, const char *s2);
//int strncasecmp (const char *s1, const char *s2, size_t n);
int strcasecmp (const char *s1, const char *s2);
int strncasecmp (const char *s1, const char *s2, size_t n);
/* the specific includes for BBOB */
#include "bbobStructures.h"
......@@ -42,9 +42,9 @@ Loosely inspired by fgeneric.m, the matlab version
void writeNewIndexEntry(ParamStruct PARAMS);
void addIndexEntry(ParamStruct PARAMS);
void addDatIndexEntry(ParamStruct PARAMS);
void writeFinalData(ParamStruct PARAMS, LastEvalStruct BestFEval, unsigned int lastWriteEval, LastEvalStruct LastEval, double Fopt);
void writeBestF(char * dataFile, unsigned int BestFEval, double Fopt);
void sprintData(FILE* fout, unsigned int evals, double F, double bestF, double Fnoisy, double bestFnoisy, double * x, double Fopt);
void writeFinalData(ParamStruct PARAMS, LastEvalStruct BestFEval, double lastWriteEval, LastEvalStruct LastEval, double Fopt);
void writeBestF(char * dataFile, double BestFEval, double Fopt);
void sprintData(FILE* fout, double evals, double F, double bestF, double Fnoisy, double bestFnoisy, double * x, double Fopt);
void writeDataHeader(char * dataFile, double Fopt);
ParamStruct setNextDataFile(ParamStruct PARAMS, unsigned int isAllParamsMatching);
......@@ -53,7 +53,8 @@ ParamStruct setNextDataFile(ParamStruct PARAMS, unsigned int isAllParamsMatching
unsigned int isInitDone;
unsigned int trialid;
double Fopt, fTrigger;
unsigned int evalsTrigger, idxEvalsTrigger, idxDIMEvalsTrigger;
double evalsTrigger;
unsigned int idxEvalsTrigger, idxDIMEvalsTrigger;
int idxFTrigger;
double maxFunEvalsFactor = 1e6;
......@@ -63,7 +64,7 @@ ParamStruct setNextDataFile(ParamStruct PARAMS, unsigned int isAllParamsMatching
bbobFunction actFunc = NULL;
LastEvalStruct LastEval, BestFEval;
unsigned int lastWriteEval;
double lastWriteEval;
/* The default values for the LastEvalStruct structures */
LastEvalStruct lastEvalInit = {0, DBL_MAX, DBL_MAX, DBL_MAX, {0}, 0};
......@@ -93,7 +94,6 @@ LastEvalStruct lastEvalInit = {0, DBL_MAX, DBL_MAX, DBL_MAX,
/* The default values for the ParamStruct structures */
/* UGLY cut-and-paste, don't forget to modify all if you modify one */
ParamStruct DefaultParam = {
"row", /* char inputFormat[4]; */
"not-specified", /* char algName[DefaultStringLength]; */
"", /* char comments[LongDefaultStringLength]; */
".", /* char dataPath[DefaultStringLength]; */
......@@ -115,7 +115,6 @@ ParamStruct DefaultParam = {
};
ParamStruct PreviousPARAMS = {
"row", /* char inputFormat[4]; */
"not-specified", /* char algName[DefaultStringLength]; */
"", /* char comments[LongDefaultStringLength]; */
".", /* char dataPath[DefaultStringLength]; */
......@@ -137,10 +136,9 @@ ParamStruct PreviousPARAMS = {
};
ParamStruct CurrentPARAMS = {
"row", /* char inputFormat[4]; */
"not-specified", /* char algName[DefaultStringLength]; */
"", /* char comments[LongDefaultStringLength]; */
"", /* char dataPath[DefaultStringLength]; */
".", /* char dataPath[DefaultStringLength]; */
"bbobexp", /* char filePrefix[DefaultStringLength]; */
0, /* unsigned int DIM; */
1e-8, /* double precision; AKA deltaFTarget */
......@@ -221,8 +219,7 @@ double fgeneric_initialize_(ParamStruct PARAMS)
/* checks if no current run is still up */
if (actFunc != NULL) {
WARNING("Calling fgeneric_initialize while an experiment is still running (DIM %d, function %d, instance %d)\nCalling fgeneric_finalize to close it properly", CurrentPARAMS.DIM, CurrentPARAMS.funcId, CurrentPARAMS.instanceId);
fgeneric_finalize_();
//FGENERIC_FINALIZE();
fgeneric_finalize();
}
/* OK, here actFunc is NULL, and we can start a new experiment */
/* First check the dimenstion, and (re)allocate memory if necessary */
......@@ -282,10 +279,6 @@ double fgeneric_initialize_(ParamStruct PARAMS)
PARAMS.runCounter = 1;
if ( strcasecmp(PARAMS.inputFormat, "row") )
// if ( strcmpi(PARAMS.inputFormat, "row") )
ERROR("Only ROW format allowed in the C version");
/* now the paths, file names, and different headers to write */
/************************************************************/
if ( strlen(PARAMS.dataPath) == 0 )
......@@ -352,8 +345,7 @@ double fgeneric_initialize_(ParamStruct PARAMS)
% and of the final function evaluation. It closes the data files.
% RETURN the best true fitness value ever obtained.
*/
double fgeneric_finalize_()
//FGENERIC_FINALIZE()
double fgeneric_finalize()
{
if (actFunc == NULL)
ERROR("Finalization process of fgeneric is called before the initialization process has occurred.");
......@@ -373,7 +365,7 @@ double fgeneric_finalize_()
/* returns the number of function
evaluations since FGENERIC('initialize', ...) was called.
*/
unsigned int fgeneric_evaluations()
double fgeneric_evaluations()
{
if (actFunc == NULL)
WARNING("fgeneric_evaluations: fgeneric_initialized has not been called.");
......@@ -407,15 +399,15 @@ double fgeneric_ftarget()
/* returns the maximum number
of function evaluations.
*/
unsigned int fgeneric_maxevals(unsigned int DIM)
double fgeneric_maxevals(unsigned int DIM)
{
return maxFunEvalsFactor * DIM;
return (double)maxFunEvalsFactor * DIM;
}
/* Sets the seed of the noise (in benchmarkshelper) with a modulo 1e9 (constraint
* on the seed provided to benchmarkshelper, this can happen when using time(NULL)).
*/
void fgeneric_noiseseed(unsigned int seed)
void fgeneric_noiseseed(unsigned long seed)
{
seed = seed % 1000000000;
setNoiseSeed(seed, seed);
......@@ -440,8 +432,7 @@ void fgeneric_evaluate_vector(double * XX, unsigned int howMany, double * result
/* a simple loop over the arrays */
for (i=0; i<howMany; i++)
{
*resultTmp = fgeneric_evaluate_(XXtmp);
//FGENERIC_EVALUATE(XXtmp);
*resultTmp = fgeneric_evaluate(XXtmp);
resultTmp++;
XXtmp += DIM;
}
......@@ -450,12 +441,11 @@ void fgeneric_evaluate_vector(double * XX, unsigned int howMany, double * result
/* returns the fitness value of X,
% X being the decision variable vector
*/
double fgeneric_evaluate_(double * X)
//FGENERIC_EVALUATE(double * X)
double fgeneric_evaluate(double * X)
{
unsigned int i, boolImprovement = 0;
double Fvalue, Ftrue;
unsigned int evalsj;
double evalsj;
FILE * dataFileId;
FILE * hdataFileId;
TwoDoubles res;
......@@ -632,7 +622,7 @@ void addDatIndexEntry(ParamStruct PARAMS)
/* --------------------------------------------------------------- */
/* complete the data file with unwritten information */
void writeFinalData(ParamStruct PARAMS, LastEvalStruct BestFEval, unsigned int lastWriteEval, LastEvalStruct LastEval, double Fopt)
void writeFinalData(ParamStruct PARAMS, LastEvalStruct BestFEval, double lastWriteEval, LastEvalStruct LastEval, double Fopt)
{
FILE * dataFileId; /* historical name */
FILE * indexFileId; /* historical name */
......@@ -660,16 +650,16 @@ void writeFinalData(ParamStruct PARAMS, LastEvalStruct BestFEval, unsigned int
/* now the index file */
indexFileId = bbobOpenFile(PARAMS.indexFile);
fprintf(indexFileId, ":%d|%.1e", LastEval.num, BestFEval.F - Fopt - PARAMS.precision);
fprintf(indexFileId, ":%.0f|%.1e", LastEval.num, BestFEval.F - Fopt - PARAMS.precision);
fclose(indexFileId);
}
/* ------------------------------------------------------------------ */
/* rewrite the data file with the information about the best F. */
void writeBestF(char * dataFile, unsigned int BestFEval, double Fopt)
void writeBestF(char * dataFile, double BestFEval, double Fopt)
{
printf("Calling writeBestF with %s, %d, %g\n", dataFile, BestFEval, Fopt);
printf("Calling writeBestF with %s, %.0f, %g\n", dataFile, BestFEval, Fopt);
}
/* %TODO: this function lacks efficiency.
......@@ -729,13 +719,14 @@ buffertmp = '';
/* ------------------------------------------------------------------- */
/* write a formatted line into a data file */
void sprintData(FILE* fout, unsigned int evals, double F, double bestF, double Fnoisy, double bestFnoisy, double * x, double Fopt)
void sprintData(FILE* fout, double evals, double F, double bestF, double Fnoisy, double bestFnoisy, double * x, double Fopt)
{
unsigned int i;
fprintf(fout, "%d",evals);
fprintf(fout, "%.0f", evals);
fprintf(fout, " %+10.9e %+10.9e %+10.9e %+10.9e", F-Fopt, bestF-Fopt, Fnoisy, bestFnoisy);
for (i=0; i<DIM; i++)
fprintf(fout, " %+5.4e",x[i]);
if (DIM < 22)
for (i=0; i<DIM; i++)
fprintf(fout, " %+5.4e",x[i]);
fprintf(fout, "\n");
}
......
......@@ -7,7 +7,7 @@
#include "bbobStructures.h" /* Include all declarations for BBOB calls */
void init_bbob_ (int *dim,int *fun,char *output_path)
void init_bbob (int *dim,int *fun,char *output_path)
{
/* declare your own internal stuff */
......
#------------------------------------------------------------------------
# pCMALib: a parallel fortran 90 library for the evolution strategy with
# covariance matrix adaptation
# Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
# MOSAIC group, ETH Zurich, Switzerland
#-------------------------------------------------------------------------
# output data is saved into this folder (relative to workdir)
OUTPUT_FOLDER = LJ_with_BFGS
# Dimension of the problem
DIMENSIONS = 108
# if special bounds should be used for initalization of population
USE_INIT_BOUNDS = true
# Upper bounds for initialization
INIT_UBOUNDS = 1.5
# Lower bounds for initialization
INIT_LBOUNDS = -1.5
# Upper bounds on all dimensions
ALLDIM_UBOUNDS = 4
# Lower bounds on all dimensions
ALLDIM_LBOUNDS = -4
#the global optimum
GLOBAL_MIN = -173.92
#use the Lennard Jones potential with compression as target function
USE_LJ = true
#size of the inital sigma relative to the bounds - only used if ABS_SIGMA is not set
REL_SIGMA = 0.5
# the maximum fold increase the population is allowed to grow to compared to
# the inital population size
MAXINCSIZE = 10
# factor by which the population size is increased every restart
INCPOPSIZE = 1.1
# Successful run if global_min -f(x) < accuracy
ACCURACY = 1.E-3
# if multi restart CMA (IPOP) should be used
RESTART_CMA = true
# if BFGS should be used to assist CMA. This is still in development!
BFGS_USE = true
# 1 = replace X values by local minimum X\\&& 2 = replace F values with F
# values at local minimum
BFGS_POSITION = 1
# (0) restart randomly within bounds, (1) restart from point of
# convergence, (2) restart from same startpoint all the time
RESTART_TYPE = 0
#------------------------------------------------------------------------
# pCMALib: a parallel fortran 90 library for the evolution strategy with
# covariance matrix adaptation
# Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
# MOSAIC group, ETH Zurich, Switzerland
#-------------------------------------------------------------------------
# output data is saved into this folder (relative to workdir)
OUTPUT_FOLDER = LJ_with_BFGS
# Dimension of the problem
DIMENSIONS = 108
# if special bounds should be used for initalization of population
USE_INIT_BOUNDS = true
# Upper bounds for initialization
INIT_UBOUNDS = 1.5
# Lower bounds for initialization
INIT_LBOUNDS = -1.5
# Upper bounds on all dimensions
ALLDIM_UBOUNDS = 4
# Lower bounds on all dimensions
ALLDIM_LBOUNDS = -4
#the global optimum
GLOBAL_MIN = -173.92
#use the Lennard Jones potential with compression as target function
USE_LJ = true
#size of the inital sigma relative to the bounds - only used if ABS_SIGMA is not set
REL_SIGMA = 0.5
# the maximum fold increase the population is allowed to grow to compared to
# the inital population size
MAXINCSIZE = 10
# factor by which the population size is increased every restart
INCPOPSIZE = 1.1
# Successful run if global_min -f(x) < accuracy
ACCURACY = 1.E-3
# if multi restart CMA (IPOP) should be used
RESTART_CMA = true
# if BFGS should be used to assist CMA. This is still in development!
BFGS_USE = true
# 1 = replace X values by local minimum X\\&& 2 = replace F values with F
# values at local minimum
BFGS_POSITION = 1
# (0) restart randomly within bounds, (1) restart from point of
# convergence, (2) restart from same startpoint all the time
RESTART_TYPE = 0
#------------------------------------------------------------------------
# pCMALib: a parallel fortran 90 library for the evolution strategy with
# covariance matrix adaptation
# Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
# MOSAIC group, ETH Zurich, Switzerland
#-------------------------------------------------------------------------
# output data is saved into this folder (relative to workdir)
OUTPUT_FOLDER = water_pscma
# Dimension of the problem
DIMENSIONS = 48
# Upper bounds on all dimensions
ALLDIM_UBOUNDS = 6
# Lower bounds on all dimensions
ALLDIM_LBOUNDS = -6
#the global optimum
GLOBAL_MIN = -305.5190
# use the TIP4P water potential as target function
USE_TIP = true
# the maximum fold increase the population is allowed to grow to compared to
# the inital population size
ACCURACY = 1.E-3
#Switch PS-CMA-ES on or off
PSCMA = true
# Intervall length between PSO updates
PSOFREQ = 1000
# how much the PSO update influences the covariance
PSOWEIGHT = 0.8
#------------------------------------------------------------------------
# pCMALib: a parallel fortran 90 library for the evolution strategy with
# covariance matrix adaptation
# Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
# MOSAIC group, ETH Zurich, Switzerland
#-------------------------------------------------------------------------
# output data is saved into this folder (relative to workdir)
OUTPUT_FOLDER = water_pscma
# Dimension of the problem
DIMENSIONS = 48
# Upper bounds on all dimensions
ALLDIM_UBOUNDS = 6
# Lower bounds on all dimensions
ALLDIM_LBOUNDS = -6
#the global optimum
GLOBAL_MIN = -305.5190
# use the TIP4P water potential as target function
USE_TIP = true
# the maximum fold increase the population is allowed to grow to compared to
# the inital population size
ACCURACY = 1.E-3
#Switch PS-CMA-ES on or off
PSCMA = true
# Intervall length between PSO updates
PSOFREQ = 1000
# how much the PSO update influences the covariance
PSOWEIGHT = 0.8
......@@ -49,7 +49,7 @@
#ifdef __BBOB
EXTERNAL :: benchmark_bobb
EXTERNAL :: benchmark_bbob
#endif
EXTERNAL :: LJ
EXTERNAL :: TIP4P
......@@ -90,14 +90,14 @@
CALL cmaes_start(DoubleFunnel,options%xstart,options%insigma)
ELSEIF (options%use_CEC .EQ. .TRUE.) THEN
CALL cmaes_start(benchmark,options%xstart,options%insigma)
ELSEIF (options%use_BOBB .EQ. .TRUE.) THEN
ELSEIF (options%use_BBOB .EQ. .TRUE.) THEN
#ifdef __BBOB
CALL cmaes_start(benchmark_bobb,options%xstart,options%insigma)
CALL cmaes_start(benchmark_bbob,options%xstart,options%insigma)
#else
WRITE(*,*) 'You compiled CMA without BBOB and therefore cant call it'
#endif
ELSEIF (options%use_MATFUNC .EQ. .TRUE.) THEN
#ifdef __HAVE_MATLAB__
#ifdef __HAVE_MATLAB__
CALL cmaes_start(matlab_objfunc,options%xstart,options%insigma)
#else
WRITE(*,*) 'You compiled CMA without MATLAB support and therefore cant call it'
......
!-------------------------------------------------------------------------
! Function : cmaes_funcwrap
!-------------------------------------------------------------------------
!
! Purpose : This is a wrapper function
!
! Input : x(:) (R) Vector to evaluate fitfun for
! N (I) Dimension of x(:)
! fitfun External Subroutine defined in cma.f90
!
! Output : Function value at x(:)
!
! Remarks :
!
! References : cma.f90
!
! Revisions :
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Function : cmaes_funcwrap
!-------------------------------------------------------------------------
!
! Purpose : This is a wrapper function
!
! Input : x(:) (R) Vector to evaluate fitfun for
! N (I) Dimension of x(:)
! fitfun External Subroutine defined in cma.f90
!
! Output : Function value at x(:)
!
! Remarks :
!
! References : cma.f90
!
! Revisions :
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! pCMALib: a parallel fortran 90 library for the evolution strategy with
! covariance matrix adaptation
! Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
! MOSAIC group, ETH Zurich, Switzerland
!-------------------------------------------------------------------------
#ifdef __SP
REAL FUNCTION cmaes_funcwrap(x,N,fitfun)
#else
DOUBLE PRECISION FUNCTION cmaes_funcwrap(x,N,fitfun)
#endif
USE cmaes_param_mod
USE cmaes_opts_mod
USE cmaes_mod
IMPLICIT NONE
!-------------------------------------------------------------------------
! Parameters
!-------------------------------------------------------------------------
INTEGER,INTENT(in) :: N
REAL(MK),DIMENSION(N),INTENT(in) :: x
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
EXTERNAL :: fitfun
!-------------------------------------------------------------------------
! Local Variables
!-------------------------------------------------------------------------
REAL(MK),DIMENSION(1) :: temp1
REAL(MK),DIMENSION(N,1) :: temp2
temp2(:,1) = x
temp1(1) = 0.0_MK
countEval = countEval + 1
CALL fitfun(temp1,temp2,N,1,options%LBounds,options%UBounds)
!CALL fitfun(temp1,temp2,N,1)
cmaes_funcwrap = temp1(1)
RETURN
! covariance matrix adaptation
! Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
! MOSAIC group, ETH Zurich, Switzerland
!-------------------------------------------------------------------------
#ifdef __SP
REAL FUNCTION cmaes_funcwrap(x,N,fitfun)
#else
DOUBLE PRECISION FUNCTION cmaes_funcwrap(x,N,fitfun)
#endif
USE cmaes_param_mod
USE cmaes_opts_mod
USE cmaes_mod
IMPLICIT NONE
!-------------------------------------------------------------------------
! Parameters
!-------------------------------------------------------------------------
INTEGER,INTENT(in) :: N
REAL(MK),DIMENSION(N),INTENT(in) :: x
!-------------------------------------------------------------------------
! Externals
!-------------------------------------------------------------------------
EXTERNAL :: fitfun
!-------------------------------------------------------------------------
! Local Variables
!-------------------------------------------------------------------------
REAL(MK),DIMENSION(1) :: temp1
REAL(MK),DIMENSION(N,1) :: temp2
temp2(:,1) = x
temp1(1) = 0.0_MK
countEval = countEval + 1
CALL fitfun(temp1,temp2,N,1,options%LBounds,options%UBounds)
!CALL fitfun(temp1,temp2,N,1)
cmaes_funcwrap = temp1(1)
RETURN
END FUNCTION
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
!-------------------------------------------------------------------------
! Function : cmaes_myprctile
!-------------------------------------------------------------------------
!
! Purpose : This function computes the percentiles in vector perc
! from vector inar
!
! Input : inar(:) (R) Vector containing data to compute
! percentiles from
! perc(:) (R) Vector containing percentiles to compute
!
! N,rN (I) Dimension of vecor inar(N) and perc(rN)
!
! Output : res(:) (R) Vector containing computed percentiles
!
! Remarks :
!
! References : this function is a Fortran translation of cmaes_myprctile in
! cmaes.m, lines 1618ff.
!
! Revisions :
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! Function : cmaes_myprctile
!-------------------------------------------------------------------------
!
! Purpose : This function computes the percentiles in vector perc
! from vector inar
!
! Input : inar(:) (R) Vector containing data to compute
! percentiles from
! perc(:) (R) Vector containing percentiles to compute
!
! N,rN (I) Dimension of vecor inar(N) and perc(rN)
!
! Output : res(:) (R) Vector containing computed percentiles
!
! Remarks :
!
! References : this function is a Fortran translation of cmaes_myprctile in
! cmaes.m, lines 1618ff.
!
! Revisions :
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
! pCMALib: a parallel fortran 90 library for the evolution strategy with
! covariance matrix adaptation
! Christian L. Mueller, Benedikt Baumgartner, Georg Ofenbeck
! MOSAIC group, ETH Zurich, Switzerland
!-------------------------------------------------------------------------
SUBROUTINE cmaes_myprctile(inar, N, perc, rN, res)
USE cmaes_param_mod
USE tool_mrgrnk_mod
IMPLICIT NONE
!-------------------------------------------------------------------------
! Arguments
!-------------------------------------------------------------------------
INTEGER,INTENT(in) :: N
INTEGER,INTENT(in) :: rN
REAL(MK),INTENT(in) :: inar(N)
REAL(MK),INTENT(in) :: perc(rN)
REAL(MK),DIMENSION(rN),INTENT(out) :: res
!-------------------------------------------------------------------------
! Local Variables
!-------------------------------------------------------------------------
REAL(MK),DIMENSION(N) :: sar, availablepercentiles
INTEGER,DIMENSION(N) :: idx ! Used to call SORTC correctly (l.52)
INTEGER :: i,k
availablepercentiles = 0.0_MK
!-------------------------------------------------------------------------
! Sort inar in ascending order
!-------------------------------------------------------------------------
CALL mrgrnk(inar,idx)
DO i = 1, N
sar(i) = inar(idx(i))
END DO
DO i = 1, rN
IF(perc(i) .LE. (100.0_MK*0.5_MK/real(N))) THEN
res(i) = sar(1)
ELSE IF(perc(i) .GE. (100.0_MK*(real(N)-0.5_MK)/real(N)) ) THEN
res(i) = sar(N)
ELSE
!-------------------------------------------------------------------
! Find largest index smaller than required percentile
!-------------------------------------------------------------------
DO k = 1,N
availablepercentiles(k) = 100.0 * (real(k)-0.5) / real(N)
END DO
find:DO k = 1, N
IF(availablepercentiles(k) .GE. perc(i)) EXIT find
END DO find
k = k - 1
!-------------------------------------------------------------------
! Interpolate linearly
!-------------------------------------------------------------------
res(i) = sar(k) + (sar(k+1)-sar(k)) * (perc(i)