27 #define innerIter 1000
28 #define outerIter 1000
30 void eplb(
double * x,
double *root,
int * steps,
double * v,
int n,
double z,
double lambda0)
34 int rho_1, rho_2, rho, rho_T, rho_S;
35 int V_i_b, V_i_e, V_i;
36 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
37 double s_1, s_2, s, s_T, s_S, v_max, temp;
38 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
46 printf(
"\n z should be nonnegative!");
53 s_1=x[V_i]=v_max=fabs(v[0]);
63 x[V_i]=fabs(v[i]); s_1+= x[V_i]; rho_1++;
82 lambda_1=0; lambda_2=v_max;
85 rho_2=0; s_2=0; f_lambda_2=-z;
89 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
95 i=V_i_b; j=V_i_e; rho=0; s=0;
97 while( (i <= V_i_e) && (x[i] <= lambda) ){
100 while( (j>=V_i_b) && (x[j] > lambda) ){
107 temp=x[i]; x[i]=x[j]; x[j]=temp;
112 rho=V_i_e-j; rho+=rho_2; s+=s_2;
113 f_lambda=s-rho*lambda-z;
115 if ( fabs(f_lambda)< delta ){
120 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
122 V_i_e=j; V_i=V_i_e-V_i_b+1;
125 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
127 V_i_b=i; V_i=V_i_e-V_i_b+1;
149 lambda_T=lambda_1 + f_lambda_1 /rho_1;
151 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
152 lambda_T=lambda_2 + f_lambda_2 /rho_2;
156 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
158 if (fabs(lambda_T-lambda_S) <= delta){
159 lambda=lambda_T; flag=1;
164 lambda=(lambda_T+lambda_S)/2;
170 while( (i <= V_i_e) && (x[i] <= lambda) ){
176 while( (j>=V_i_b) && (x[j] > lambda) ){
177 if (x[j] > lambda_S){
186 if (x[i] > lambda_S){
197 temp=x[i]; x[i]=x[j]; x[j]=temp;
202 s_S+=s_2; rho_S+=rho_2;
205 f_lambda_S=s_S-rho_S*lambda_S-z;
206 f_lambda=s-rho*lambda-z;
207 f_lambda_T=s_T-rho_T*lambda_T-z;
211 if ( fabs(f_lambda)< delta ){
216 if ( fabs(f_lambda_S)< delta ){
218 lambda=lambda_S; flag=1;
221 if ( fabs(f_lambda_T)< delta ){
223 lambda=lambda_T; flag=1;
234 lambda_2=lambda; s_2=s; rho_2=rho;
237 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
238 f_lambda_1=f_lambda_T;
242 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
245 while( (j>=V_i_b) && (x[j] > lambda_T) ){
253 V_i_b=i; V_i=V_i_e-V_i_b+1;
256 lambda_1=lambda; s_1=s; rho_1=rho;
259 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
260 f_lambda_2=f_lambda_S;
264 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
267 while( (j>=V_i_b) && (x[j] > lambda_S) ){
275 V_i_e=j; V_i=V_i_e-V_i_b+1;
279 lambda=(s - z)/ rho; flag=1;
299 void epp1(
double *x,
double *v,
int n,
double rho)
318 void epp2(
double *x,
double *v,
int n,
double rho)
342 void eppInf(
double *x,
double * c,
int * iter_step,
double *v,
int n,
double rho,
double c0)
350 eplb(x, c, &steps, v, n, rho, c0);
359 void zerofind(
double *root,
int * iterStep,
double v,
double p,
double c,
double x0)
361 double x, f, fprime, p1=p-1, pp;
366 *root=0; *iterStep=0;
return;
370 *root=v; * iterStep=0;
return;
374 if ( (x0 <v) && (x0>0) )
390 fprime=1 + c* p1 * pp / x;
407 if ( (x<0) || (x>v)){
435 if ( fabs(f) <= delta){
441 if (step>=innerIter){
442 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);
443 printf(
"\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
454 double norm(
double * v,
double p,
int n)
468 return( pow(t, 1/p) );
471 void eppO(
double *x,
double * cc,
int * iter_step,
double *v,
int n,
double rho,
double p)
473 int i, *flag, bisStep, newtonStep=0, totoalStep=0;
474 double vq=0, epsilon, vmax=0, vmin=1e10;
475 double q=1/(1-1/p), c, c1, c2, root, f, xp;
481 flag=(
int *)malloc(
sizeof(
int)*n);
525 iter_step[0]=iter_step[1]=0;
541 epsilon=(vq -rho)/ vq;
544 if ( log((1-epsilon) * vmin) - (p-1) * log( epsilon* vmin ) >= 709 )
560 eppInf(x, cc, iter_step, v, n, rho, 0);
566 c1= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
567 c2= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
571 c2= (1-epsilon) * vmax / pow(epsilon* vmax, p-1);
572 c1= (1-epsilon) * vmin / pow(epsilon* vmin, p-1);
580 if (fabs(c1-c2) <= delta){
595 zerofind(&root, &newtonStep, v[i], p, c, x[i]);
597 temp=fabs(root-x[i]);
602 totoalStep+=newtonStep;
607 f=rho * pow(xp, 1-p) - c;
609 if ( fabs(f)<=delta || fabs(c1-c2)<=delta )
613 if ( (x_diff <=delta) && (p_n==0) )
620 if ( (x_diff <=delta) && (p_n==1) )
628 if (bisStep>=outerIter){
631 if ( fabs(c1-c2) <=delta * c2 )
634 printf(
"\n The number of bisection steps exceed %d.", outerIter);
635 printf(
"\n c1=%e, c2=%e, x_diff=%e, f=%e",c1,c2,x_diff,f);
636 printf(
"\n If you meet with this problem, please contact Jun Liu (j.liu@asu.edu). Thanks!");
662 iter_step[0]=bisStep;
663 iter_step[1]=totoalStep;
666 void epp(
double *x,
double * c,
int * iter_step,
double * v,
int n,
double rho,
double p,
double c0){
668 printf(
"\n rho should be non-negative!");
675 iter_step[0]=iter_step[1]=0;
681 iter_step[0]=iter_step[1]=0;
685 eppInf(x, c, iter_step, v, n, rho, c0);
687 eppO(x, c, iter_step, v, n, rho, p);
690 #endif //USE_GPL_SHOGUN