1 #include <stdio.h>
2 #include <string.h>
3 #include <stdlib.h>
4 #include <math.h>
5 #include <assert.h>
6 #include <time.h>
7 
8 #include "array.h"
9 #include "linalg.h"
10 #include "lib_clinalg.h"
11 #include "lib_funcs.h"
12 
main(int argc,char * argv[])13 int main(int argc, char * argv[])
14 {
15     srand(time(NULL));
16     if ((argc != 2)){
17        printf("Correct function call = ./dmrgprodbench type\n");
18        printf("Options for type name include: \n");
19        printf("\t \"0\": time vs r comparison between dmrg and standard prod-round (DOESN't WORK!!!!!)\n");
20        printf("\t \"1\": run for profile purposes\n");
21        printf("\t \"2\": time vs r comparison for just dmrg returns (rank,exact_rank,rounded_rank,als_time,,direct_time)\n");
22        return 0;
23     }
24 
25 
26     size_t dim = 4;
27     enum poly_type ptype = LEGENDRE;
28     struct OpeOpts * opts = ope_opts_alloc(ptype);
29     struct OneApproxOpts * qmopts = one_approx_opts_alloc(POLYNOMIAL,opts);
30     struct MultiApproxOpts * fopts = multi_approx_opts_alloc(dim);
31     multi_approx_opts_set_all_same(fopts,qmopts);
32     struct BoundingBox * bds = bounding_box_init_std(dim);
33 
34     if (strcmp(argv[1],"0") == 0){
35         printf("DOESNT WORK YET!");
36         size_t ranks[5] = {1,4,4,4,1};
37         size_t maxorder = 40;
38 
39         size_t nrepeat = 1;
40         size_t nranks = 40;
41         size_t jj,kk;
42 
43         double * results = calloc_double(nranks*4);
44         for (jj = 0; jj < nranks; jj++){
45             for (kk = 1; kk < dim; kk++){
46                 ranks[kk] = 2 + 1*jj;
47             }
48 
49             results[jj] = (double) ranks[1];
50 
51             printf("base ranks are ");
52             iprint_sz(dim+1,ranks);
53 
54             double time = 0.0;
55             double time2 = 0.0;
56             double mrank = 0.0;
57             size_t ii;
58 
59             for (ii = 0; ii < nrepeat; ii++){
60 
61                 struct FunctionTrain * a = NULL;
62                 a = function_train_poly_randu(ptype,bds,ranks,maxorder);
63 
64                 struct FunctionTrain * start = NULL;
65                 a = function_train_constant(1.0,fopts);
66 
67                 struct FunctionTrain * acopy = function_train_copy(a);
68 
69                 clock_t tic = clock();
70                 struct FunctionTrain * at = function_train_product(a,a);
71                 struct FunctionTrain * p = function_train_copy(at);
72                 struct FunctionTrain * pr = function_train_round(p,1e-10,fopts);
73                 clock_t toc = clock();
74                 time += (double)(toc - tic) / CLOCKS_PER_SEC;
75 
76                 printf("Actual ranks are ");
77                 iprint_sz(dim+1,function_train_get_ranks(at));
78                 mrank += pr->ranks[2];
79 
80                 tic = clock();
81                 struct FunctionTrain * finish = dmrg_product(start,acopy,acopy,1e-5,10,1e-10,0,fopts);
82                 toc = clock();
83                 time2 += (double)(toc - tic) / CLOCKS_PER_SEC;
84 
85                 printf("finish ranks "); iprint_sz(dim+1,function_train_get_ranks(finish));
86                 //double diff = function_train_relnorm2diff(pr,finish);
87                 double diff = function_train_relnorm2diff(at,finish);
88                 printf("diff = %G\n",diff);
89                 assert (diff*diff < 1e-10);
90                 function_train_free(p); p = NULL;
91                 function_train_free(pr); p = NULL;
92                 function_train_free(finish); finish=NULL;
93 
94                 function_train_free(a); a = NULL;
95                 function_train_free(at); at = NULL;
96                 function_train_free(start); start = NULL;
97             }
98             time /= nrepeat;
99             time2 /= nrepeat;
100             mrank /= nrepeat;
101             printf("Average times are (%G,%G) max rank is %G \n", time,time2,mrank);
102 
103             results[nranks+jj] = time;
104             results[2*nranks+jj] = time2;
105             results[3*nranks+jj] = mrank;
106 
107         }
108         darray_save(nranks,4,results,"time_vs_rank.dat",1);
109 
110     }
111     else{
112         size_t ranks[5] = {1,15,15,15,1};
113         size_t maxorder = 10;
114 
115 
116         struct FunctionTrain * a = function_train_poly_randu(ptype,bds,ranks,
117                                                              maxorder);
118         struct FunctionTrain * start = function_train_constant(1.0,fopts);
119 
120         clock_t tic = clock();
121         struct FunctionTrain * finish = dmrg_product(start,a,a,1e-5,10,1e-10,0,fopts);
122         clock_t toc = clock();
123         printf("timing is %G\n", (double)(toc - tic) / CLOCKS_PER_SEC);
124 
125         function_train_free(a); a = NULL;
126         function_train_free(start); start = NULL;
127         function_train_free(finish); finish = NULL;
128     }
129     if (strcmp(argv[1],"2") == 0){
130         size_t ranks[5] = {1,4,4,4,1};
131         size_t maxorder = 10;
132 
133         size_t nrepeat = 5;
134         size_t nranks = 20;
135         size_t jj,kk;
136 
137         // storage (rank, rank_exact, rank_rounded, time_als, time_rounded)
138         double * results = calloc_double(nranks*5);
139         for (jj = 0; jj < nranks; jj++){
140             for (kk = 1; kk < dim; kk++){
141                 ranks[kk] = 2 + 1*jj;
142             }
143 
144             results[jj] = (double) ranks[1];
145 
146             printf("base ranks are ");
147             iprint_sz(dim+1,ranks);
148 
149             double time2 = 0.0;
150             double time = 0.0;
151             double mrank = 0.0;
152             size_t ii;
153             for (ii = 0; ii < nrepeat; ii++){
154 
155                 struct FunctionTrain * a = NULL;
156                 a = function_train_poly_randu(ptype,bds,ranks,maxorder);
157                 struct FunctionTrain * start = NULL;
158                 start = function_train_constant(1.0,fopts);
159 
160                 clock_t tic = clock();
161                 struct FunctionTrain * finish = dmrg_product(start,a,a,1e-5,10,1e-10,0,fopts);
162                 clock_t toc = clock();
163                 time2 += (double)(toc - tic) / CLOCKS_PER_SEC;
164 
165 
166                 tic = clock();
167                 struct FunctionTrain * doubled = function_train_product(a,a);
168                 struct FunctionTrain * rounded = function_train_round(doubled,1e-10,fopts);
169                 toc = clock();
170                 time += (double)(toc - tic) / CLOCKS_PER_SEC;
171 
172                 results[nranks+jj] = 1.0*function_train_get_maxrank(doubled);
173                 mrank += function_train_get_maxrank(rounded);
174 
175                 /* printf("finish ranks"); */
176                 /* iprint_sz(dim+1,function_train_get_ranks(finish)); */
177 
178                 /* printf("rounded ranks"); */
179                 /* iprint_sz(dim+1,function_train_get_ranks(rounded)); */
180 
181                 double diff = function_train_relnorm2diff(rounded,finish);
182                 printf("diff = %G\n",diff);
183 
184                 function_train_free(doubled); doubled = NULL;
185                 function_train_free(rounded); rounded = NULL;
186 
187                 function_train_free(finish); finish=NULL;
188                 function_train_free(a); a = NULL;
189                 function_train_free(start); start = NULL;
190             }
191             time2 /= nrepeat;
192             time /= nrepeat;
193             mrank /= nrepeat;
194             results[2*nranks+jj] = mrank;
195             printf("Average time (ALS/DIRECT) (%G/%G) max rank is %G \n",time2,time,mrank);
196 
197             results[3*nranks+jj] = time2;
198             results[4*nranks+jj] = time;
199 
200         }
201         darray_save(nranks,5,results,"time_vs_rank_dmrg.dat",1);
202     }
203 
204 
205     bounding_box_free(bds); bds = NULL;
206     multi_approx_opts_free(fopts);
207     one_approx_opts_free_deep(&qmopts);
208     return 0;
209 }
210