25 #define innerIter 1000
26 #define outerIter 1000
28 void eplb(
double * x,
double *root,
int * steps,
double * v,
int n,
double z,
double lambda0)
32 int rho_1, rho_2, rho, rho_T, rho_S;
33 int V_i_b, V_i_e, V_i;
34 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
35 double s_1, s_2, s, s_T, s_S, v_max, temp;
36 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
44 printf(
"\n z should be nonnegative!");
51 s_1=x[V_i]=v_max=fabs(v[0]);
61 x[V_i]=fabs(v[i]); s_1+= x[V_i]; rho_1++;
80 lambda_1=0; lambda_2=v_max;
83 rho_2=0; s_2=0; f_lambda_2=-z;
87 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
93 i=V_i_b; j=V_i_e; rho=0; s=0;
95 while( (i <= V_i_e) && (x[i] <= lambda) ){
98 while( (j>=V_i_b) && (x[j] > lambda) ){
105 temp=x[i]; x[i]=x[j]; x[j]=temp;
110 rho=V_i_e-j; rho+=rho_2; s+=s_2;
111 f_lambda=s-rho*lambda-z;
113 if ( fabs(f_lambda)<
delta ){
118 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
120 V_i_e=j; V_i=V_i_e-V_i_b+1;
123 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
125 V_i_b=i; V_i=V_i_e-V_i_b+1;
147 lambda_T=lambda_1 + f_lambda_1 /rho_1;
149 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
150 lambda_T=lambda_2 + f_lambda_2 /rho_2;
154 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
156 if (fabs(lambda_T-lambda_S) <=
delta){
157 lambda=lambda_T; flag=1;
162 lambda=(lambda_T+lambda_S)/2;
168 while( (i <= V_i_e) && (x[i] <= lambda) ){
174 while( (j>=V_i_b) && (x[j] > lambda) ){
175 if (x[j] > lambda_S){
184 if (x[i] > lambda_S){
195 temp=x[i]; x[i]=x[j]; x[j]=temp;
200 s_S+=s_2; rho_S+=rho_2;
203 f_lambda_S=s_S-rho_S*lambda_S-z;
204 f_lambda=s-rho*lambda-z;
205 f_lambda_T=s_T-rho_T*lambda_T-z;
209 if ( fabs(f_lambda)<
delta ){
214 if ( fabs(f_lambda_S)<
delta ){
216 lambda=lambda_S; flag=1;
219 if ( fabs(f_lambda_T)<
delta ){
221 lambda=lambda_T; flag=1;
232 lambda_2=lambda; s_2=s; rho_2=rho;
235 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
236 f_lambda_1=f_lambda_T;
240 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
243 while( (j>=V_i_b) && (x[j] > lambda_T) ){
251 V_i_b=i; V_i=V_i_e-V_i_b+1;
254 lambda_1=lambda; s_1=s; rho_1=rho;
257 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
258 f_lambda_2=f_lambda_S;
262 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
265 while( (j>=V_i_b) && (x[j] > lambda_S) ){
273 V_i_e=j; V_i=V_i_e-V_i_b+1;
277 lambda=(s - z)/ rho; flag=1;
297 void epp1(
double *x,
double *v,
int n,
double rho)
316 void epp2(
double *x,
double *v,
int n,
double rho)
340 void eppInf(
double *x,
double * c,
int * iter_step,
double *v,
int n,
double rho,
double c0)
348 eplb(x, c, &steps, v, n, rho, c0);
357 void zerofind(
double *root,
int * iterStep,
double v,
double p,
double c,
double x0)
359 double x, f, fprime, p1=p-1, pp;
364 *root=0; *iterStep=0;
return;
368 *root=v; * iterStep=0;
return;
372 if ( (x0 <v) && (x0>0) )
388 fprime=1 + c* p1 * pp / x;
405 if ( (x<0) || (x>v)){
433 if ( fabs(f) <=
delta){
440 printf(
"\n The number of steps exceed %d, in finding the root for f(x)= x + c x^{p-1} - v, 0< x< v.",
innerIter);
441 printf(
"\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
452 double norm(
double * v,
double p,
int n)
466 return( pow(t, 1/p) );
469 void eppO(
double *x,
double * cc,
int * iter_step,
double *v,
int n,
double rho,
double p)
471 int i, *flag, bisStep, newtonStep=0, totoalStep=0;
472 double vq=0,
epsilon, vmax=0, vmin=1e10;
473 double q=1/(1-1/p), c, c1, c2, root, f, xp;
479 flag=(
int *)malloc(
sizeof(
int)*n);
523 iter_step[0]=iter_step[1]=0;
542 if ( log((1-
epsilon) * vmin) - (p-1) * log(
epsilon* vmin ) >= 709 )
558 eppInf(x, cc, iter_step, v, n, rho, 0);
578 if (fabs(c1-c2) <=
delta){
593 zerofind(&root, &newtonStep, v[i], p, c, x[i]);
595 temp=fabs(root-x[i]);
600 totoalStep+=newtonStep;
605 f=rho * pow(xp, 1-p) - c;
611 if ( (x_diff <=
delta) && (p_n==0) )
618 if ( (x_diff <=
delta) && (p_n==1) )
629 if ( fabs(c1-c2) <=
delta * c2 )
632 printf(
"\n The number of bisection steps exceed %d.",
outerIter);
633 printf(
"\n c1=%e, c2=%e, x_diff=%e, f=%e",c1,c2,x_diff,f);
634 printf(
"\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
660 iter_step[0]=bisStep;
661 iter_step[1]=totoalStep;
664 void epp(
double *x,
double * c,
int * iter_step,
double * v,
int n,
double rho,
double p,
double c0){
666 printf(
"\n rho should be non-negative!");
673 iter_step[0]=iter_step[1]=0;
679 iter_step[0]=iter_step[1]=0;
683 eppInf(x, c, iter_step, v, n, rho, c0);
685 eppO(x, c, iter_step, v, n, rho, p);
double norm(double *v, double p, int n)
void epp2(double *x, double *v, int n, double rho)
void eppInf(double *x, double *c, int *iter_step, double *v, int n, double rho, double c0)
void epp1(double *x, double *v, int n, double rho)
void eppO(double *x, double *cc, int *iter_step, double *v, int n, double rho, double p)
void zerofind(double *root, int *iterStep, double v, double p, double c, double x0)
static const float64_t epsilon
void eplb(double *x, double *root, int *steps, double *v, int n, double z, double lambda0)
void epp(double *x, double *c, int *iter_step, double *v, int n, double rho, double p, double c0)