42 void epsp(
double * x,
double *root,
int * stepsp,
double * v,
43 int n,
double z,
double lambda0)
47 int rho_1, rho_2, rho, rho_T, rho_S;
48 int V_i_b, V_i_e, V_i;
49 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
50 double s_1, s_2, s, s_T, s_S, v_max, temp;
51 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
57 printf(
"\n z should be nonnegative!");
72 lambda_1=v_max - z; lambda_2=v_max;
87 f_lambda_1=s_1-rho_1* lambda_1 -z;
89 rho_2=0; s_2=0; f_lambda_2=-z;
93 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
99 i=V_i_b; j=V_i_e; rho=0; s=0;
101 while( (i <= V_i_e) && (x[i] <= lambda) ){
104 while( (j>=V_i_b) && (x[j] > lambda) ){
111 temp=x[i]; x[i]=x[j]; x[j]=temp;
116 rho=V_i_e-j; rho+=rho_2; s+=s_2;
117 f_lambda=s-rho*lambda-z;
119 if ( fabs(f_lambda)<
delta ){
124 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
126 V_i_e=j; V_i=V_i_e-V_i_b+1;
129 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
131 V_i_b=i; V_i=V_i_e-V_i_b+1;
153 lambda_T=lambda_1 + f_lambda_1 /rho_1;
155 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
156 lambda_T=lambda_2 + f_lambda_2 /rho_2;
160 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
162 if (fabs(lambda_T-lambda_S) <=
delta){
163 lambda=lambda_T; flag=1;
168 lambda=(lambda_T+lambda_S)/2;
174 while( (i <= V_i_e) && (x[i] <= lambda) ){
180 while( (j>=V_i_b) && (x[j] > lambda) ){
181 if (x[j] > lambda_S){
190 if (x[i] > lambda_S){
201 temp=x[i]; x[i]=x[j]; x[j]=temp;
206 s_S+=s_2; rho_S+=rho_2;
209 f_lambda_S=s_S-rho_S*lambda_S-z;
210 f_lambda=s-rho*lambda-z;
211 f_lambda_T=s_T-rho_T*lambda_T-z;
215 if ( fabs(f_lambda)<
delta ){
220 if ( fabs(f_lambda_S)<
delta ){
222 lambda=lambda_S; flag=1;
225 if ( fabs(f_lambda_T)<
delta ){
227 lambda=lambda_T; flag=1;
238 lambda_2=lambda; s_2=s; rho_2=rho;
241 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
242 f_lambda_1=f_lambda_T;
246 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
249 while( (j>=V_i_b) && (x[j] > lambda_T) ){
257 V_i_b=i; V_i=V_i_e-V_i_b+1;
260 lambda_1=lambda; s_1=s; rho_1=rho;
263 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
264 f_lambda_2=f_lambda_S;
268 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
271 while( (j>=V_i_b) && (x[j] > lambda_S) ){
279 V_i_e=j; V_i=V_i_e-V_i_b+1;
283 lambda=(s - z)/ rho; flag=1;