44 void epsp(
double * x,
double *root,
int * stepsp,
double * v,
45 int n,
double z,
double lambda0)
49 int rho_1, rho_2, rho, rho_T, rho_S;
50 int V_i_b, V_i_e, V_i;
51 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
52 double s_1, s_2, s, s_T, s_S, v_max, temp;
53 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
59 printf(
"\n z should be nonnegative!");
74 lambda_1=v_max - z; lambda_2=v_max;
89 f_lambda_1=s_1-rho_1* lambda_1 -z;
91 rho_2=0; s_2=0; f_lambda_2=-z;
95 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
101 i=V_i_b; j=V_i_e; rho=0; s=0;
103 while( (i <= V_i_e) && (x[i] <= lambda) ){
106 while( (j>=V_i_b) && (x[j] > lambda) ){
113 temp=x[i]; x[i]=x[j]; x[j]=temp;
118 rho=V_i_e-j; rho+=rho_2; s+=s_2;
119 f_lambda=s-rho*lambda-z;
121 if ( fabs(f_lambda)<
delta ){
126 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
128 V_i_e=j; V_i=V_i_e-V_i_b+1;
131 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
133 V_i_b=i; V_i=V_i_e-V_i_b+1;
155 lambda_T=lambda_1 + f_lambda_1 /rho_1;
157 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
158 lambda_T=lambda_2 + f_lambda_2 /rho_2;
162 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
164 if (fabs(lambda_T-lambda_S) <=
delta){
165 lambda=lambda_T; flag=1;
170 lambda=(lambda_T+lambda_S)/2;
176 while( (i <= V_i_e) && (x[i] <= lambda) ){
182 while( (j>=V_i_b) && (x[j] > lambda) ){
183 if (x[j] > lambda_S){
192 if (x[i] > lambda_S){
203 temp=x[i]; x[i]=x[j]; x[j]=temp;
208 s_S+=s_2; rho_S+=rho_2;
211 f_lambda_S=s_S-rho_S*lambda_S-z;
212 f_lambda=s-rho*lambda-z;
213 f_lambda_T=s_T-rho_T*lambda_T-z;
217 if ( fabs(f_lambda)<
delta ){
222 if ( fabs(f_lambda_S)<
delta ){
224 lambda=lambda_S; flag=1;
227 if ( fabs(f_lambda_T)<
delta ){
229 lambda=lambda_T; flag=1;
240 lambda_2=lambda; s_2=s; rho_2=rho;
243 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
244 f_lambda_1=f_lambda_T;
248 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
251 while( (j>=V_i_b) && (x[j] > lambda_T) ){
259 V_i_b=i; V_i=V_i_e-V_i_b+1;
262 lambda_1=lambda; s_1=s; rho_1=rho;
265 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
266 f_lambda_2=f_lambda_S;
270 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
273 while( (j>=V_i_b) && (x[j] > lambda_S) ){
281 V_i_e=j; V_i=V_i_e-V_i_b+1;
285 lambda=(s - z)/ rho; flag=1;
void epsp(double *x, double *root, int *stepsp, double *v, int n, double z, double lambda0)