00001
00009 #include <stdio.h>
00010 #include <float.h>
00011 #include <math.h>
00012 #include "rlp.h"
00013
00018 #define POS(L,C,NC) ((L)*(NC)+(C))
00019
00030 static void fmprint(double* matrix,int nrows,int ncols,FILE* file)
00031 {
00032 int i,j;
00033
00034 fprintf(file,"\n +");
00035 for(i=0;i<ncols;i++)
00036 fprintf(file,"--------------+");
00037
00038 fputs("\n",file);
00039
00040 for(i=0;i<nrows;i++)
00041 {
00042 fprintf(file," |");
00043 for(j=0;j<ncols;j++)
00044 fprintf(file," %12.4f |",matrix[POS(i,j,ncols)]);
00045
00046 fputs("\n",file);
00047 }
00048
00049 fprintf(file," +");
00050 for(i=0;i<ncols;i++)
00051 fprintf(file,"--------------+");
00052
00053 fputs("\n\n",file);
00054 }
00055
00056
00057
00070 static double minimumc(double* matrix,int nrows,int ncols,int col,int* row)
00071 {
00072 int i;
00073 double res=DBL_MAX;
00074
00075 for(i=0;i<nrows;i++)
00076 if(matrix[POS(i,col,ncols)]<res)
00077 {
00078 res=matrix[POS(i,col,ncols)];
00079 *row=i;
00080 }
00081
00082 return res;
00083 }
00084
00085
00086
00098 static double minimumr(double* matrix,int ncols,int row,int* col)
00099 {
00100 int i;
00101 double res=DBL_MAX;
00102
00103 for(i=0;i<ncols;i++)
00104 if(matrix[POS(row,i,ncols)]<res)
00105 {
00106 res=matrix[POS(row,i,ncols)];
00107 *col=i;
00108 }
00109
00110 return res;
00111 }
00112
00113
00114
00115 int simplex(double* a,int n,int m,FILE* file)
00116 {
00117 double aux;
00118 int pos,nm2,err;
00119
00120 nm2=n+m+2;
00121
00122 aux=minimumc(&a[nm2],m,nm2,n+m,&pos);
00123 err=0;
00124
00125 if(file) fmprint(a,m+1,nm2,file);
00126
00127 if(aux<0) err=simplexd(a,n,m,pos+1,file);
00128
00129 if(err) return err;
00130 else
00131 {
00132 aux=minimumr(a,nm2,0,&pos);
00133
00134 if(aux<0) err=simplexp(a,n,m,pos,file);
00135
00136 if(err) return err;
00137 else return 0;
00138 }
00139 }
00140
00141
00142
00143 int simplexp(double* a,int n,int m,int pos,FILE* file)
00144 {
00145 int nm,pivotc,pivotr,i,j,k;
00146 double pivot,coef,aux,r;
00147
00148 nm=n+m;
00149 k=nm+2;
00150 aux=-1;
00151 pivotc=pos;
00152
00153 while(aux<0)
00154 {
00155 pivotr=-1;
00156 aux=DBL_MAX;
00157
00158 for(i=1;i<m+1;i++)
00159 if(a[POS(i,pivotc,k)]>0)
00160 {
00161 r=a[POS(i,nm,k)]/a[POS(i,pivotc,k)];
00162 if(r<aux)
00163 {
00164 aux=r;
00165 pivotr=i;
00166 }
00167 }
00168
00169 if(pivotr==-1) return 1;
00170
00171 pivot=a[POS(pivotr,pivotc,k)];
00172
00173 for(i=0;i<nm+1;i++)
00174 a[POS(pivotr,i,k)]/=pivot;
00175
00176 for(i=1;i<m+1;i++)
00177 if(i!=pivotr)
00178 {
00179 coef=-a[POS(i,pivotc,k)];
00180 for(j=0;j<nm+1;j++)
00181 a[POS(i,j,k)]+=coef*a[POS(pivotr,j,k)];
00182 }
00183
00184 aux=DBL_MAX;
00185 coef=-a[POS(0,pivotc,k)];
00186
00187 a[POS(pivotr,nm+1,k)]=pivotc;
00188
00189 for(i=0;i<nm+1;i++)
00190 {
00191 a[POS(0,i,k)]+=coef*a[POS(pivotr,i,k)];
00192 if(a[POS(0,i,k)]<aux)
00193 {
00194 aux=a[POS(0,i,k)];
00195 pivotc=i;
00196 }
00197 }
00198
00199 if(file) fmprint(a,m+1,k,file);
00200 }
00201
00202 return 0;
00203 }
00204
00205
00206
00207 int simplexd(double* a,int n,int m,int pos,FILE* file)
00208 {
00209 int pivotc,pivotr,nm,i,j,k;
00210 double aux,pivot,r,coef;
00211
00212 nm=n+m;
00213 k=nm+2;
00214
00215 aux=-1;
00216 pivotr=pos;
00217
00218 while(aux<0)
00219 {
00220 pivotc=-1;
00221 aux=DBL_MAX;
00222
00223 for(i=0;i<nm;i++)
00224 if(a[POS(pivotr,i,k)]<0)
00225 {
00226 r=a[POS(0,i,k)]/a[POS(pivotr,i,k)];
00227 if(fabs(r)<aux)
00228 {
00229 aux=r;
00230 pivotc=i;
00231 }
00232 }
00233
00234 if(pivotc==-1) return 1;
00235
00236 pivot=a[POS(pivotr,pivotc,k)];
00237
00238 for(i=0;i<nm+1;i++)
00239 a[POS(pivotr,i,k)]/=pivot;
00240
00241 aux=a[POS(pivotr,nm,k)];
00242
00243 for(i=0;i<m+1;i++)
00244 if(i!=pivotr)
00245 {
00246 coef=a[POS(i,pivotc,k)]>0?a[POS(i,pivotc,k)]:-a[POS(i,pivotc,k)];
00247 for(j=0;j<nm+1;j++)
00248 {
00249 a[POS(i,j,k)]+=coef*a[POS(pivotr,j,k)];
00250 }
00251
00252 if(a[POS(i,nm,k)]<aux)
00253 {
00254 aux=a[POS(i,nm,k)];
00255 pivotr=i;
00256 }
00257 }
00258
00259 a[POS(pivotr,nm+1,k)]=pivotc;
00260
00261 if(file) fmprint(a,m+1,k,file);
00262 }
00263
00264 return 0;
00265 }