00001 #include "math.h"
00002 #include "stdlib.h"
00003 #include "stdio.h"
00004
00005 #include "adolc.h"
00006 #include "sparse/sparsedrivers.h"
00007
00008 #define repnum 10
00009
00010
00011
00012
00013 #include <sys/time.h>
00014 #include "ColPackHeaders.h"
00015
00016 using namespace ColPack;
00017
00018 double k_getTime() {
00019 struct timeval v;
00020 struct timezone z;
00021 gettimeofday(&v, &z);
00022 return ((double)v.tv_sec)+((double)v.tv_usec)/1000000;
00023 }
00024
00025
00026
00027
00028 using namespace std;
00029
00030 #include <list>
00031 #include <map>
00032 #include <string>
00033 #include <vector>
00034
00035
00036
00037
00038 #define tag_f 1
00039 #define tag_red 2
00040 #define tag_HP 3
00041
00042 #define tag_c 4
00043
00044
00045 void init_dim(int *n, int *m);
00046 void init_startpoint(double *x, int n);
00047 double feval(double *x, int n);
00048 adouble feval_ad(double *x, int n);
00049 void ceval(double *x, double *c, int n);
00050 void ceval_ad(double *x, adouble *c, int n);
00051 adouble feval_ad_mod(double *x, int n);
00052 adouble feval_ad_modHP(double *x, int n);
00053
00054 void printmat(char* kette, int n, int m, double** M);
00055 void printmatint(char* kette, int n, int m, int** M);
00056 void printmatint_c(char* kette, int m,unsigned int** M);
00057
00058
00059 int main()
00060 {
00061 int i, j, k, l, sum, n, m, nnz, direct = 1, found;
00062 double f;
00063 double *x, *c;
00064 adouble fad, *xad, *cad;
00065
00066
00067 double** J;
00068
00069
00070 int tape_stats[11];
00071 int num;
00072 FILE *fp_JP;
00073
00074 double **Seed_J;
00075 double **Jc;
00076 int p_J;
00077
00078 int recover = 1;
00079 int jac_vec = 1;
00080 int compute_full = 1;
00081 int output_sparsity_pattern_J = 0;
00082
00083
00084
00085
00086
00087
00088 double t_f_1, t_f_2, div_f=0, div_c=0, div_JP=0, div_JP2=0, div_Seed=0, div_Seed_C=0, div_Jc=0, div_Jc_C=0, div_rec=0, div_rec_C=0, div_J=0;
00089
00090 unsigned int *rind;
00091 unsigned int *cind;
00092 double *values;
00093
00094
00095
00096
00097
00098
00099
00100
00101 init_dim(&n,&m);
00102
00103 printf(" n = %d m = %d\n",n,m);
00104
00105 x = (double*) malloc(n*sizeof(double));
00106 c = (double*) malloc(m*sizeof(double));
00107 cad = new adouble[m];
00108
00109 init_startpoint(x,n);
00110
00111 t_f_1 = k_getTime();
00112 for(i=0;i<repnum;i++)
00113 f = feval(x,n);
00114 t_f_2 = k_getTime();
00115 div_f = (t_f_2 - t_f_1)*1.0/repnum;
00116 printf("XXX The time needed for function evaluation: %10.6f \n \n", div_f);
00117
00118
00119 t_f_1 = k_getTime();
00120 for(i=0;i<repnum;i++)
00121 ceval(x,c,n);
00122 t_f_2 = k_getTime();
00123 div_c = (t_f_2 - t_f_1)*1.0/repnum;
00124 printf("XXX The time needed for constraint evaluation: %10.6f \n \n", div_c);
00125
00126
00127 trace_on(tag_f);
00128
00129 fad = feval_ad(x, n);
00130
00131 fad >>= f;
00132
00133 trace_off();
00134
00135 trace_on(tag_c);
00136
00137 ceval_ad(x, cad, n);
00138
00139 for(i=0;i<m;i++)
00140 cad[i] >>= f;
00141
00142 trace_off();
00143
00144
00145 tapestats(tag_c,tape_stats);
00146 printf("\n independents %d\n",tape_stats[0]);
00147 printf(" dependents %d\n",tape_stats[1]);
00148 printf(" operations %d\n",tape_stats[5]);
00149 printf(" buffer size %d\n",tape_stats[4]);
00150 printf(" maxlive %d\n",tape_stats[2]);
00151 printf(" valstack size %d\n\n",tape_stats[3]);
00152
00153
00154
00155
00156
00157 div_J = -1;
00158
00159 if(compute_full == 1)
00160 {
00161 J = myalloc2(m,n);
00162
00163 t_f_1 = k_getTime();
00164 jacobian(tag_c,m,n,x,J);
00165 t_f_2 = k_getTime();
00166 div_J = (t_f_2 - t_f_1);
00167
00168 printf("XXX The time needed for full Jacobian: %10.6f \n \n", div_J);
00169 printf("XXX runtime ratio: %10.6f \n \n", div_J/div_c);
00170
00171
00172
00173 fp_JP = fopen("jac_full.mtx","w");
00174
00175 fprintf(fp_JP,"%d %d\n",m,n);
00176
00177 for (i=0;i<m;i++)
00178 {
00179 for (j=0;j<n;j++)
00180 if(J[i][j]!=0.0) fprintf(fp_JP,"%d %d %10.6f\n",i,j,J[i][j] );
00181 }
00182 fclose(fp_JP);
00183 }
00184
00185
00186 printf("XXX THE 4 STEP TO COMPUTE SPARSE MATRICES USING ColPack \n \n");
00187
00188
00189 unsigned int *rb=NULL;
00190 unsigned int *cb=NULL;
00191 unsigned int **JP=NULL;
00192 int ctrl[2];
00193
00194 JP = (unsigned int **) malloc(m*sizeof(unsigned int*));
00195 ctrl[0] = 0; ctrl[1] = 0;
00196
00197
00198 t_f_1 = k_getTime();
00199 jac_pat(tag_c, m, n, x, JP, ctrl);
00200 t_f_2 = k_getTime();
00201 div_JP = (t_f_2 - t_f_1);
00202
00203 printf("XXX STEP 1: The time needed for Jacobian pattern: %10.6f \n \n", div_JP);
00204 printf("XXX STEP 1: runtime ratio: %10.6f \n \n", div_JP/div_c);
00205
00206
00207 nnz = 0;
00208 for (i=0;i<m;i++)
00209 nnz += JP[i][0];
00210
00211 printf(" nnz %d \n",nnz);
00212 printf(" hier 1a\n");
00213
00214
00215
00216
00217
00218 double tg_C;
00219 int dummy;
00220
00221 t_f_1 = k_getTime();
00222
00223 BipartiteGraphPartialColoringInterface * gGraph = new BipartiteGraphPartialColoringInterface();
00224 gGraph->GenerateSeedJacobian(JP, m, n, &Seed_J, &dummy, &p_J,
00225 "RIGHT_PARTIAL_DISTANCE_TWO", "SMALLEST_LAST");
00226
00227 t_f_2 = k_getTime();
00228 tg_C = t_f_2 - t_f_1;
00229
00230
00231 printf("XXX STEP 2: The time needed for Seed generation: %10.6f \n \n", tg_C);
00232 printf("XXX STEP 2: runtime ratio: %10.6f \n \n", tg_C/div_c);
00233
00234
00235
00236
00237
00238
00239
00240
00241 if (jac_vec == 1)
00242 {
00243
00244 Jc = myalloc2(m,p_J);
00245 t_f_1 = k_getTime();
00246 printf(" hier 1\n");
00247 fov_forward(tag_c,m,n,p_J,x,Seed_J,c,Jc);
00248 printf(" hier 2\n");
00249 t_f_2 = k_getTime();
00250 div_Jc = (t_f_2 - t_f_1);
00251
00252 printf("XXX STEP 3: The time needed for Jacobian-matrix product: %10.6f \n \n", div_Jc);
00253 printf("XXX STEP 3: runtime ratio: %10.6f \n \n", div_Jc/div_c);
00254
00255
00256 }
00257
00258
00259
00260
00261
00262 if (recover == 1)
00263 {
00264
00265
00266 JacobianRecovery1D jr1d;
00267
00268 printf("m = %d, n = %d, p_J = %d \n",m,n,p_J);
00269
00270
00271
00272 t_f_1 = k_getTime();
00273 jr1d.RecoverD2Cln_CoordinateFormat (gGraph, Jc, JP, &rind, &cind, &values);
00274 t_f_2 = k_getTime();
00275 div_rec_C = (t_f_2 - t_f_1);
00276
00277 printf("XXX STEP 4: The time needed for Recovery: %10.6f \n \n", div_rec_C);
00278 printf("XXX STEP 4: runtime ratio: %10.6f \n \n", div_rec_C/div_c);
00279
00280
00281
00282 fp_JP = fopen("jac_recovered.mtx","w");
00283
00284 fprintf(fp_JP,"%d %d %d\n",m,n,nnz);
00285
00286 for (i=0;i<nnz;i++)
00287 {
00288 fprintf(fp_JP,"%d %d %10.6f\n",rind[i],cind[i],values[i] );
00289 }
00290 fclose(fp_JP);
00291
00292 }
00293
00294
00295
00296
00297
00298
00299 free(JP);
00300 delete[] cad;
00301 free(c);
00302 free(x);
00303 delete gGraph;
00304 if(jac_vec == 1) {
00305 myfree2(Jc);
00306 }
00307 if(compute_full == 1)
00308 {
00309 myfree2(J);
00310 }
00311
00312 return 0;
00313
00314 }
00315
00316
00317
00318
00319
00320
00321 void printmat(char* kette, int n, int m, double** M)
00322 { int i,j;
00323
00324 printf("%s \n",kette);
00325 for(i=0; i<n ;i++)
00326 {
00327 printf("\n %d: ",i);
00328 for(j=0;j<m ;j++)
00329 printf(" %10.4f ", M[i][j]);
00330 }
00331 printf("\n");
00332 }
00333
00334 void printmatint(char* kette, int n, int m, int** M)
00335 { int i,j;
00336
00337 printf("%s \n",kette);
00338 for(i=0; i<n ;i++)
00339 {
00340 printf("\n %d: ",i);
00341 for(j=0;j<m ;j++)
00342 printf(" %d ", M[i][j]);
00343 }
00344 printf("\n");
00345 }
00346
00347
00348 void printmatint_c(char* kette, int m,unsigned int** M)
00349 { int i,j;
00350
00351 printf("%s \n",kette);
00352 for (i=0;i<m;i++)
00353 {
00354 printf("\n %d (%d): ",i,M[i][0]);
00355 for (j=1;j<=M[i][0];j++)
00356 printf("\t%d ",M[i][j]);
00357 }
00358 printf("\n");
00359 }