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