Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #ifndef EPSP_SLEP
00018 #define EPSP_SLEP
00019
00020 #include <stdlib.h>
00021 #include <stdio.h>
00022 #include <time.h>
00023 #include <math.h>
00024
00025 #define delta 1e-12
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 void epsp(double * x, double *root, int * stepsp, double * v,
00043 int n, double z, double lambda0)
00044 {
00045
00046 int i, j, flag=0;
00047 int rho_1, rho_2, rho, rho_T, rho_S;
00048 int V_i_b, V_i_e, V_i;
00049 double lambda_1, lambda_2, lambda_T, lambda_S, lambda;
00050 double s_1, s_2, s, s_T, s_S, v_max, temp;
00051 double f_lambda_1, f_lambda_2, f_lambda, f_lambda_T, f_lambda_S;
00052 int iter_step=0;
00053
00054
00055
00056 if (z< 0){
00057 printf("\n z should be nonnegative!");
00058 exit(-1);
00059 }
00060
00061
00062
00063
00064
00065 v_max=v[0];
00066 for (i=1;i<n; i++){
00067 if (v[i] > v_max)
00068 v_max=v[i];
00069 }
00070
00071
00072 lambda_1=v_max - z; lambda_2=v_max;
00073
00074
00075
00076
00077 V_i=0;s_1=0; rho_1=0;
00078 for (i=0;i<n;i++){
00079 if (v[i] > lambda_1){
00080 x[V_i]=v[i];
00081
00082 s_1+=x[V_i]; rho_1++;
00083
00084 V_i++;
00085 }
00086 }
00087 f_lambda_1=s_1-rho_1* lambda_1 -z;
00088
00089 rho_2=0; s_2=0; f_lambda_2=-z;
00090 V_i_b=0; V_i_e=V_i-1;
00091
00092 lambda=lambda0;
00093 if ( (lambda<lambda_2) && (lambda> lambda_1) ){
00094
00095
00096
00097
00098
00099 i=V_i_b; j=V_i_e; rho=0; s=0;
00100 while (i <= j){
00101 while( (i <= V_i_e) && (x[i] <= lambda) ){
00102 i++;
00103 }
00104 while( (j>=V_i_b) && (x[j] > lambda) ){
00105 s+=x[j];
00106 j--;
00107 }
00108 if (i<j){
00109 s+=x[i];
00110
00111 temp=x[i]; x[i]=x[j]; x[j]=temp;
00112 i++; j--;
00113 }
00114 }
00115
00116 rho=V_i_e-j; rho+=rho_2; s+=s_2;
00117 f_lambda=s-rho*lambda-z;
00118
00119 if ( fabs(f_lambda)< delta ){
00120 flag=1;
00121 }
00122
00123 if (f_lambda <0){
00124 lambda_2=lambda; s_2=s; rho_2=rho; f_lambda_2=f_lambda;
00125
00126 V_i_e=j; V_i=V_i_e-V_i_b+1;
00127 }
00128 else{
00129 lambda_1=lambda; rho_1=rho; s_1=s; f_lambda_1=f_lambda;
00130
00131 V_i_b=i; V_i=V_i_e-V_i_b+1;
00132 }
00133
00134 if (V_i==0){
00135
00136
00137
00138
00139 lambda=(s - z)/ rho;
00140 flag=1;
00141 }
00142
00143
00144
00145
00146
00147 }
00148
00149 while (!flag){
00150 iter_step++;
00151
00152
00153 lambda_T=lambda_1 + f_lambda_1 /rho_1;
00154 if(rho_2 !=0){
00155 if (lambda_2 + f_lambda_2 /rho_2 > lambda_T)
00156 lambda_T=lambda_2 + f_lambda_2 /rho_2;
00157 }
00158
00159
00160 lambda_S=lambda_2 - f_lambda_2 *(lambda_2-lambda_1)/(f_lambda_2-f_lambda_1);
00161
00162 if (fabs(lambda_T-lambda_S) <= delta){
00163 lambda=lambda_T; flag=1;
00164 break;
00165 }
00166
00167
00168 lambda=(lambda_T+lambda_S)/2;
00169
00170 s_T=s_S=s=0;
00171 rho_T=rho_S=rho=0;
00172 i=V_i_b; j=V_i_e;
00173 while (i <= j){
00174 while( (i <= V_i_e) && (x[i] <= lambda) ){
00175 if (x[i]> lambda_T){
00176 s_T+=x[i]; rho_T++;
00177 }
00178 i++;
00179 }
00180 while( (j>=V_i_b) && (x[j] > lambda) ){
00181 if (x[j] > lambda_S){
00182 s_S+=x[j]; rho_S++;
00183 }
00184 else{
00185 s+=x[j]; rho++;
00186 }
00187 j--;
00188 }
00189 if (i<j){
00190 if (x[i] > lambda_S){
00191 s_S+=x[i]; rho_S++;
00192 }
00193 else{
00194 s+=x[i]; rho++;
00195 }
00196
00197 if (x[j]> lambda_T){
00198 s_T+=x[j]; rho_T++;
00199 }
00200
00201 temp=x[i]; x[i]=x[j]; x[j]=temp;
00202 i++; j--;
00203 }
00204 }
00205
00206 s_S+=s_2; rho_S+=rho_2;
00207 s+=s_S; rho+=rho_S;
00208 s_T+=s; rho_T+=rho;
00209 f_lambda_S=s_S-rho_S*lambda_S-z;
00210 f_lambda=s-rho*lambda-z;
00211 f_lambda_T=s_T-rho_T*lambda_T-z;
00212
00213
00214
00215 if ( fabs(f_lambda)< delta ){
00216
00217 flag=1;
00218 break;
00219 }
00220 if ( fabs(f_lambda_S)< delta ){
00221
00222 lambda=lambda_S; flag=1;
00223 break;
00224 }
00225 if ( fabs(f_lambda_T)< delta ){
00226
00227 lambda=lambda_T; flag=1;
00228 break;
00229 }
00230
00231
00232
00233
00234
00235
00236
00237 if (f_lambda <0){
00238 lambda_2=lambda; s_2=s; rho_2=rho;
00239 f_lambda_2=f_lambda;
00240
00241 lambda_1=lambda_T; s_1=s_T; rho_1=rho_T;
00242 f_lambda_1=f_lambda_T;
00243
00244 V_i_e=j; i=V_i_b;
00245 while (i <= j){
00246 while( (i <= V_i_e) && (x[i] <= lambda_T) ){
00247 i++;
00248 }
00249 while( (j>=V_i_b) && (x[j] > lambda_T) ){
00250 j--;
00251 }
00252 if (i<j){
00253 x[j]=x[i];
00254 i++; j--;
00255 }
00256 }
00257 V_i_b=i; V_i=V_i_e-V_i_b+1;
00258 }
00259 else{
00260 lambda_1=lambda; s_1=s; rho_1=rho;
00261 f_lambda_1=f_lambda;
00262
00263 lambda_2=lambda_S; s_2=s_S; rho_2=rho_S;
00264 f_lambda_2=f_lambda_S;
00265
00266 V_i_b=i; j=V_i_e;
00267 while (i <= j){
00268 while( (i <= V_i_e) && (x[i] <= lambda_S) ){
00269 i++;
00270 }
00271 while( (j>=V_i_b) && (x[j] > lambda_S) ){
00272 j--;
00273 }
00274 if (i<j){
00275 x[i]=x[j];
00276 i++; j--;
00277 }
00278 }
00279 V_i_e=j; V_i=V_i_e-V_i_b+1;
00280 }
00281
00282 if (V_i==0){
00283 lambda=(s - z)/ rho; flag=1;
00284
00285 break;
00286 }
00287 }
00288
00289
00290 for(i=0;i<n;i++){
00291 if (v[i] > lambda)
00292 x[i]=v[i]-lambda;
00293 else
00294 x[i]=0;
00295 }
00296 *root=lambda;
00297 *stepsp=iter_step;
00298 }
00299 #endif
00300