1 // Copyright (c) 2015-2016, Massachusetts Institute of Technology
2 // Copyright (c) 2016-2017 Sandia Corporation
3 // Copyright (c) 2017 NTESS, LLC.
4
5 // This file is part of the Compressed Continuous Computation (C3) Library
6 // Author: Alex A. Gorodetsky
7 // Contact: alex@alexgorodetsky.com
8
9 // All rights reserved.
10
11 // Redistribution and use in source and binary forms, with or without modification,
12 // are permitted provided that the following conditions are met:
13
14 // 1. Redistributions of source code must retain the above copyright notice,
15 // this list of conditions and the following disclaimer.
16
17 // 2. Redistributions in binary form must reproduce the above copyright notice,
18 // this list of conditions and the following disclaimer in the documentation
19 // and/or other materials provided with the distribution.
20
21 // 3. Neither the name of the copyright holder nor the names of its contributors
22 // may be used to endorse or promote products derived from this software
23 // without specific prior written permission.
24
25 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
28 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
29 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
31 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
32 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
33 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
34 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36 //Code
37
38
39 /** \file ft.c
40 * Provides algorithms for function train manipulation
41 */
42
43 #include <stdlib.h>
44 #include <stdio.h>
45 #include <assert.h>
46 #include <math.h>
47 #include <string.h>
48
49 #include "stringmanip.h"
50 #include "array.h"
51 #include "linalg.h"
52 #include "lib_funcs.h"
53
54 #include "ft.h"
55
56 #ifndef ZEROTHRESH
57 /// @private
58 #define ZEROTHRESH 1e0* DBL_EPSILON
59 /* #define ZEROTHRESH 1e0*DBL_EPSILON */
60 /* #define ZEROTHRESH 1e-200 */
61 #endif
62
63 /* #define ZEROTHRESH 1e-3*DBL_EPSILON */
64
65 #ifndef VPREPCORE
66 /// @private
67 #define VPREPCORE 0
68 #endif
69
70 #ifndef VFTCROSS
71 /// @private
72 #define VFTCROSS 0
73 #endif
74
75
76 /// @private
77 struct FTMemSpace
78 {
79 size_t num_spaces;
80 size_t space_size; // increment 1
81 double * vals; // num_spaces x space_size
82 };
83
84 /// @private
ft_mem_space_alloc(size_t num_spaces,size_t space_size)85 struct FTMemSpace * ft_mem_space_alloc(size_t num_spaces, size_t space_size)
86 {
87 struct FTMemSpace * ftmem = malloc(sizeof(struct FTMemSpace));
88 if (ftmem == NULL){
89 fprintf(stderr, "Failure to allocate memory in FT\n");
90 return NULL;
91 }
92
93 ftmem->num_spaces = num_spaces;
94 ftmem->space_size = space_size;
95 ftmem->vals = calloc_double(num_spaces * space_size);
96
97 return ftmem;
98 }
99
100 /// @private
ft_mem_space_arr_alloc(size_t dim,size_t num_spaces,size_t space_size)101 struct FTMemSpace ** ft_mem_space_arr_alloc(size_t dim, size_t num_spaces, size_t space_size)
102 {
103 struct FTMemSpace ** ftmem = malloc(dim * sizeof(struct FTMemSpace * ));
104 if (ftmem == NULL){
105 fprintf(stderr, "Failure to allocate memory in FT\n");
106 return NULL;
107 }
108 for (size_t ii = 0; ii < dim; ii++){
109 ftmem[ii] = ft_mem_space_alloc(num_spaces,space_size);
110 }
111
112 return ftmem;
113 }
114
115 /// @private
ft_mem_space_free(struct FTMemSpace * ftmem)116 void ft_mem_space_free(struct FTMemSpace * ftmem)
117 {
118 if (ftmem != NULL){
119 free(ftmem->vals); ftmem->vals = NULL;
120 free(ftmem); ftmem = NULL;
121 }
122 }
123
124 /// @private
ft_mem_space_arr_free(size_t dim,struct FTMemSpace ** ftmem)125 void ft_mem_space_arr_free(size_t dim, struct FTMemSpace ** ftmem)
126 {
127 if (ftmem != NULL){
128 for (size_t ii = 0; ii < dim; ii++){
129 ft_mem_space_free(ftmem[ii]);
130 ftmem[ii] = NULL;
131 }
132 free(ftmem); ftmem = NULL;
133 }
134 }
135
136 /// @private
ft_mem_space_check(struct FTMemSpace * ftmem,size_t num_spaces,size_t space_size)137 int ft_mem_space_check(struct FTMemSpace * ftmem,
138 size_t num_spaces,
139 size_t space_size)
140 {
141 if (ftmem == NULL){
142 return 1;
143 }
144 if ( (ftmem->num_spaces < num_spaces) ||
145 (ftmem->space_size != space_size))
146 {
147 ftmem->num_spaces = num_spaces;
148 ftmem->space_size = space_size;
149 free(ftmem->vals); ftmem->vals = NULL;
150 ftmem->vals = calloc_double(num_spaces * space_size);
151 }
152
153 return 0;
154
155 }
156
157 /// @private
ft_mem_space_get_incr(struct FTMemSpace * ftmem)158 size_t ft_mem_space_get_incr(struct FTMemSpace * ftmem)
159 {
160 return ftmem->space_size;
161 }
162
163 /***********************************************************//**
164 Allocate space for a function_train
165
166 \param[in] dim - dimension of function train
167
168 \return ft - function train
169 ***************************************************************/
function_train_alloc(size_t dim)170 struct FunctionTrain * function_train_alloc(size_t dim)
171 {
172 /* printf("Allocating Function Train\n"); */
173 struct FunctionTrain * ft = NULL;
174 if ( NULL == (ft = malloc(sizeof(struct FunctionTrain)))){
175 fprintf(stderr, "failed to allocate memory for function train.\n");
176 exit(1);
177 }
178 ft->dim = dim;
179 ft->ranks = calloc_size_t(dim+1);
180 if ( NULL == (ft->cores = malloc(dim * sizeof(struct Qmarray *)))){
181 fprintf(stderr, "failed to allocate memory for function train.\n");
182 exit(1);
183 }
184 size_t ii;
185 for (ii = 0; ii < dim; ii++){
186 ft->cores[ii] = NULL;
187 }
188 ft->evalspace1 = NULL;
189 ft->evalspace2 = NULL;
190 ft->evalspace3 = NULL;
191 ft->evalspace4 = NULL;
192 ft->evaldd1 = NULL;
193 ft->evaldd2 = NULL;
194 ft->evaldd3 = NULL;
195 ft->evaldd4 = NULL;
196
197 /* printf("Done Allocating Function Train\n"); */
198 return ft;
199 }
200
201 /**********************************************************//**
202 Copy a function train
203
204 \param[in] a - function train to copy
205
206 \return b - function train
207 **************************************************************/
function_train_copy(const struct FunctionTrain * a)208 struct FunctionTrain * function_train_copy(const struct FunctionTrain * a)
209 {
210
211 if (a == NULL){
212 return NULL;
213 }
214 struct FunctionTrain * b = function_train_alloc(a->dim);
215
216 size_t ii;
217 for (ii = 0; ii < a->dim; ii++){
218 b->ranks[ii] = a->ranks[ii];
219 b->cores[ii] = qmarray_copy(a->cores[ii]);
220 }
221 b->ranks[a->dim] = a->ranks[a->dim];
222 return b;
223 }
224
225 /**********************************************************//**
226 Free memory allocated to a function train
227
228 \param[in,out] ft - function train to free
229 **************************************************************/
function_train_free(struct FunctionTrain * ft)230 void function_train_free(struct FunctionTrain * ft)
231 {
232 /* printf("\t Freeing Function Train\n"); */
233 if (ft != NULL){
234 /* printf("\t\t FT is not null\n"); */
235 free(ft->ranks);
236 ft->ranks = NULL;
237 size_t ii;
238 for (ii = 0; ii < ft->dim; ii++){
239 qmarray_free(ft->cores[ii]);ft->cores[ii] = NULL;
240 }
241 free(ft->cores); ft->cores=NULL;
242
243 ft_mem_space_free(ft->evalspace1); ft->evalspace1 = NULL;
244 ft_mem_space_free(ft->evalspace2); ft->evalspace2 = NULL;
245 ft_mem_space_free(ft->evalspace3); ft->evalspace3 = NULL;
246
247 ft_mem_space_arr_free(ft->dim, ft->evalspace4); ft->evalspace4 = NULL;
248
249 free_dd(ft->dim,ft->evaldd1); ft->evaldd1 = NULL;
250 /* free_dd(2*ft->dim,ft->evaldd2); ft->evaldd2 = NULL; */
251 free_dd(ft->dim,ft->evaldd3); ft->evaldd3 = NULL;
252 free_dd(ft->dim,ft->evaldd4); ft->evaldd4 = NULL;
253 free(ft); ft = NULL;
254 }
255 /* printf("\t\tDone Freeing Function Train\n"); */
256 }
257
258 struct FunctionTrain *
function_train_operate(const struct FunctionTrain * ft,struct Operator ** op)259 function_train_operate(const struct FunctionTrain * ft,
260 struct Operator ** op)
261 {
262
263 struct FunctionTrain * ft_out = function_train_alloc(ft->dim);
264 ft_out->ranks[1] = ft->ranks[1];
265 for (size_t ii = 0; ii < ft_out->dim; ii++){
266 ft_out->cores[ii] = qmarray_operate_elements(ft->cores[ii],
267 op[ii]);
268 ft_out->ranks[ii+1] = ft->ranks[ii+1];
269 }
270
271 return ft_out;
272 }
273
274
275 /***********************************************************//**
276 Serialize a function_train
277
278 \param[in,out] ser - stream to serialize to
279 \param[in] ft - function train
280 \param[in,out] totSizeIn - if NULL then serialize, if not NULL then return size
281
282 \return ptr - ser shifted by number of ytes
283 ***************************************************************/
284 unsigned char *
function_train_serialize(unsigned char * ser,struct FunctionTrain * ft,size_t * totSizeIn)285 function_train_serialize(unsigned char * ser, struct FunctionTrain * ft,
286 size_t *totSizeIn)
287 {
288
289 // dim -> nranks -> core1 -> core2 -> ... -> cored
290
291 size_t ii;
292 //dim + ranks
293 size_t totSize = sizeof(size_t) + (ft->dim+1)*sizeof(size_t);
294 size_t size_temp;
295 for (ii = 0; ii < ft->dim; ii++){
296 qmarray_serialize(NULL,ft->cores[ii],&size_temp);
297 totSize += size_temp;
298 }
299 if (totSizeIn != NULL){
300 *totSizeIn = totSize;
301 return ser;
302 }
303
304 unsigned char * ptr = ser;
305 ptr = serialize_size_t(ptr, ft->dim);
306 //printf("serializing ranks\n");
307 for (ii = 0; ii < ft->dim+1; ii++){
308 ptr = serialize_size_t(ptr,ft->ranks[ii]);
309 }
310 //printf("done serializing ranks dim\n");
311 for (ii = 0; ii < ft->dim; ii++){
312 /* printf("Serializing core (%zu/%zu)\n",ii+1,ft->dim); */
313 ptr = qmarray_serialize(ptr, ft->cores[ii],NULL);
314 }
315 return ptr;
316 }
317
318 /***********************************************************//**
319 Deserialize a function train
320
321 \param[in,out] ser - serialized function train
322 \param[in,out] ft - function_train
323
324 \return ptr - shifted ser after deserialization
325 ***************************************************************/
326 unsigned char *
function_train_deserialize(unsigned char * ser,struct FunctionTrain ** ft)327 function_train_deserialize(unsigned char * ser, struct FunctionTrain ** ft)
328 {
329 unsigned char * ptr = ser;
330
331 size_t dim;
332 ptr = deserialize_size_t(ptr, &dim);
333 /* printf("deserialized dim=%zu\n",dim); */
334 *ft = function_train_alloc(dim);
335
336 size_t ii;
337 for (ii = 0; ii < dim+1; ii++){
338 ptr = deserialize_size_t(ptr, &((*ft)->ranks[ii]));
339 }
340 /* printf("deserialized ranks = \n"); iprint_sz(dim+1,(*ft)->ranks); */
341 for (ii = 0; ii < dim; ii++){
342 ptr = qmarray_deserialize(ptr, &((*ft)->cores[ii]));
343 }
344
345 return ptr;
346 }
347
348 /***********************************************************//**
349 Save a function train to text file for portability
350
351 \param[in] ft - function train
352 \param[in,out] stream - stream to write to
353 \param[in] prec - precision with which to save
354
355 ***************************************************************/
function_train_savetxt(struct FunctionTrain * ft,FILE * stream,size_t prec)356 void function_train_savetxt(struct FunctionTrain * ft, FILE * stream,
357 size_t prec)
358 {
359
360 // dim -> nranks -> core1 -> core2 -> ... -> cored
361
362 size_t ii;
363 //dim + ranks
364
365 fprintf(stream,"%zu ",ft->dim);
366 for (ii = 0; ii < ft->dim+1; ii++){
367 fprintf(stream,"%zu ",ft->ranks[ii]);
368 }
369
370 for (ii = 0; ii < ft->dim; ii++){
371 qmarray_savetxt(ft->cores[ii],stream,prec);
372 }
373 }
374
375 /***********************************************************//**
376 Load a function train stored as a text file
377
378 \param[in,out] fp - stream from which to load
379
380 \return function train
381 ***************************************************************/
function_train_loadtxt(FILE * fp)382 struct FunctionTrain * function_train_loadtxt(FILE * fp)
383 {
384 size_t dim;
385 int num = fscanf(fp,"%zu ",&dim);
386 assert (num == 1);
387
388 struct FunctionTrain * ft = function_train_alloc(dim);
389
390 size_t ii;
391 for (ii = 0; ii < dim+1; ii++){
392 num = fscanf(fp,"%zu ",ft->ranks + ii);
393 assert (num == 1);
394 }
395 for (ii = 0; ii < dim; ii++){
396 ft->cores[ii] = qmarray_loadtxt(fp);
397 }
398
399 return ft;
400 }
401
402 /***********************************************************//**
403 Save a function train to file
404
405 \param[in] ft - function train to save
406 \param[in] filename - name of file to save to
407
408 \return success (1) or failure (0) of opening the file
409 ***************************************************************/
function_train_save(struct FunctionTrain * ft,char * filename)410 int function_train_save(struct FunctionTrain * ft, char * filename)
411 {
412
413 FILE *fp;
414 fp = fopen(filename, "w");
415 if (fp == NULL){
416 fprintf(stderr, "cat: can't open %s\n", filename);
417 return 0;
418 }
419
420 size_t totsize;
421 function_train_serialize(NULL,ft, &totsize);
422
423 unsigned char * data = malloc(totsize+sizeof(size_t));
424 if (data == NULL){
425 fprintf(stderr, "can't allocate space for saving density\n");
426 return 0;
427 }
428
429 // serialize size first!
430 unsigned char * ptr = serialize_size_t(data,totsize);
431 ptr = function_train_serialize(ptr,ft,NULL);
432
433 fwrite(data,sizeof(unsigned char),totsize+sizeof(size_t),fp);
434
435 free(data); data = NULL;
436 fclose(fp);
437 return 1;
438 }
439
440 /***********************************************************//**
441 Load a function train from a file
442
443 \param[in] filename - filename of file to load
444
445 \return ft if successfull NULL otherwise
446 ***************************************************************/
function_train_load(char * filename)447 struct FunctionTrain * function_train_load(char * filename)
448 {
449 FILE *fp;
450 fp = fopen(filename, "r");
451 if (fp == NULL){
452 fprintf(stderr, "cat: can't open %s\n", filename);
453 return NULL;
454 }
455
456 size_t totsize;
457 size_t k = fread(&totsize,sizeof(size_t),1,fp);
458 if ( k != 1){
459 printf("error reading file %s\n",filename);
460 return NULL;
461 }
462
463 unsigned char * data = malloc(totsize);
464 if (data == NULL){
465 fprintf(stderr, "can't allocate space for loading density\n");
466 return NULL;
467 }
468
469 k = fread(data,sizeof(unsigned char),totsize,fp);
470
471 struct FunctionTrain * ft = NULL;
472 function_train_deserialize(data,&ft);
473
474 free(data); data = NULL;
475 fclose(fp);
476 return ft;
477 }
478
479 /***********************************************************//**
480 Get dim
481 ***************************************************************/
function_train_get_dim(const struct FunctionTrain * ft)482 size_t function_train_get_dim(const struct FunctionTrain * ft)
483 {
484 assert(ft != NULL);
485 return ft->dim;
486 }
487
488 /***********************************************************//**
489 Get core
490 ***************************************************************/
function_train_get_core(const struct FunctionTrain * ft,size_t ind)491 struct Qmarray* function_train_get_core(const struct FunctionTrain * ft, size_t ind)
492 {
493 assert(ft != NULL);
494 assert (ind < ft->dim);
495 return ft->cores[ind];
496 }
497
498 /***********************************************************//**
499 Get rank
500 ***************************************************************/
function_train_get_rank(const struct FunctionTrain * ft,size_t index)501 size_t function_train_get_rank(const struct FunctionTrain * ft, size_t index)
502 {
503 assert (ft != NULL);
504 assert (ft->ranks != NULL);
505 assert (index < ft->dim+1);
506 return ft->ranks[index];
507 }
508
509 /***********************************************************//**
510 Get ranks
511 ***************************************************************/
function_train_get_ranks(const struct FunctionTrain * ft)512 size_t * function_train_get_ranks(const struct FunctionTrain * ft)
513 {
514 assert (ft != NULL);
515 assert (ft->ranks != NULL);
516 return ft->ranks;
517 }
518
519 /***********************************************************//**
520 Computes the maximum rank of a FT
521
522 \param[in] ft - function train
523
524 \return maxrank
525 ***************************************************************/
function_train_get_maxrank(const struct FunctionTrain * ft)526 size_t function_train_get_maxrank(const struct FunctionTrain * ft)
527 {
528 assert (ft != NULL);
529 assert (ft->ranks != NULL);
530
531 size_t maxrank = 1;
532 size_t ii;
533 for (ii = 0; ii < ft->dim+1; ii++){
534 if (ft->ranks[ii] > maxrank){
535 maxrank = ft->ranks[ii];
536 }
537 }
538
539 return maxrank;
540 }
541
542 /***********************************************************//**
543 Computes the average rank of a FT. Doesn't cound first and last ranks
544
545 \param[in] ft - function train
546
547 \return avgrank
548 ***************************************************************/
function_train_get_avgrank(const struct FunctionTrain * ft)549 double function_train_get_avgrank(const struct FunctionTrain * ft)
550 {
551
552 double avg = 0;
553 if (ft->dim == 1){
554 return 1.0;
555 }
556 else if (ft->dim == 2){
557 return (double) ft->ranks[1];
558 }
559 else{
560 size_t ii;
561 for (ii = 1; ii < ft->dim; ii++){
562 avg += (double)ft->ranks[ii];
563 }
564 }
565 return (avg / (ft->dim - 1.0));
566 }
567
568 /***********************************************************//**
569 Evaluate a new core and multiply against the previous evaluations (Left-Right)
570
571 \param[in] ft - function train
572 \param[in] core - which core to evaluate
573 \param[in] N - number of evaluations
574 \param[in] x - evaluation location
575 \param[in] incx - increment of x
576 \param[in] prev_eval - result of evaluating previous cores
577 \param[in] inc_pr - increment for previous evaluation
578 \param[in,out] new_eval_space - space for the evaluation of a core
579 \param[in] inc_new_space - increment for new evaluation
580 \param[in,out] new_eval - output location for prev_eval * new_eval_space
581 \param[in] inc_new - increment for final multiplication
582
583 \note Increments are between storage locations for each evaluation
584 ***************************************************************/
function_train_eval_next_core(struct FunctionTrain * ft,size_t core,size_t N,const double * x,size_t incx,double * prev_eval,size_t inc_pr,double * new_eval_space,size_t inc_new_space,double * new_eval,size_t inc_new)585 void function_train_eval_next_core(struct FunctionTrain * ft, size_t core, size_t N,
586 const double * x, size_t incx,
587 double * prev_eval, size_t inc_pr,
588 double * new_eval_space, size_t inc_new_space,
589 double * new_eval, size_t inc_new)
590 {
591 assert (core > 0);
592
593 size_t r1 = ft->cores[core]->nrows;
594 size_t r2 = ft->cores[core]->ncols;
595 generic_function_1darray_eval2N(
596 r1*r2, ft->cores[core]->funcs, N, x, incx, new_eval_space, inc_new_space);
597
598 for (size_t ii = 0; ii < N; ii++){
599 cblas_dgemv(CblasColMajor,CblasTrans,
600 r1,r2, 1.0,
601 new_eval_space + ii * inc_new_space, r1,
602 prev_eval + ii*inc_pr, 1, 0.0, new_eval + ii*inc_new, 1);
603 }
604
605 }
606
607 /***********************************************************//**
608 Evaluate a new core and multiply against the previous evaluations (Right-Left)
609
610 \param[in] ft - function train
611 \param[in] core - which core to evaluate
612 \param[in] N - number of evaluations
613 \param[in] x - evaluation location
614 \param[in] incx - increment of x
615 \param[in] prev_eval - result of evaluating previous cores
616 \param[in] inc_pr - increment for previous evaluation
617 \param[in,out] new_eval_space - space for the evaluation of a core
618 \param[in] inc_new_space - increment for new evaluation
619 \param[in,out] new_eval - output location for prev_eval * new_eval_space
620 \param[in] inc_new - increment for final multiplication
621
622 \note Increments are between storage locations for each evaluation
623 ***************************************************************/
function_train_eval_prev_core(struct FunctionTrain * ft,size_t core,size_t N,const double * x,size_t incx,double * prev_eval,size_t inc_pr,double * new_eval_space,size_t inc_new_space,double * new_eval,size_t inc_new)624 void function_train_eval_prev_core(struct FunctionTrain * ft, size_t core, size_t N,
625 const double * x, size_t incx,
626 double * prev_eval, size_t inc_pr,
627 double * new_eval_space, size_t inc_new_space,
628 double * new_eval, size_t inc_new)
629 {
630 assert (core > 0);
631
632 size_t r1 = ft->cores[core]->nrows;
633 size_t r2 = ft->cores[core]->ncols;
634 generic_function_1darray_eval2N(
635 r1*r2, ft->cores[core]->funcs, N, x, incx, new_eval_space, inc_new_space);
636
637 for (size_t ii = 0; ii < N; ii++){
638 cblas_dgemv(CblasColMajor,CblasNoTrans,
639 r1, r2, 1.0,
640 new_eval_space + ii * inc_new_space, r1,
641 prev_eval + ii * inc_pr, 1, 0.0, new_eval + ii*inc_new, 1);
642 }
643
644 }
645
646 /***********************************************************//**
647 Some overhead for function train evaluation functions
648 ***************************************************************/
function_train_pre_eval(struct FunctionTrain * ft,size_t * maxrank,double ** t1,double ** t2,double ** t3,size_t N)649 void function_train_pre_eval(struct FunctionTrain * ft, size_t * maxrank,
650 double ** t1,double ** t2,double ** t3, size_t N)
651 {
652
653 *maxrank = function_train_get_maxrank(ft);
654 size_t mr2 = *maxrank * *maxrank;
655 int ok = ft_mem_space_check(ft->evalspace1, N, mr2);
656 if (ok == 1){
657 ft->evalspace1 = ft_mem_space_alloc(N, mr2);
658 }
659 ok = ft_mem_space_check(ft->evalspace2, N, mr2);
660 if (ok == 1){
661 ft->evalspace2 = ft_mem_space_alloc(N, mr2);
662 }
663 ok = ft_mem_space_check(ft->evalspace3, N, mr2);
664 if (ok == 1){
665 ft->evalspace3 = ft_mem_space_alloc(N, mr2);
666 }
667
668 *t1 = ft->evalspace1->vals;
669 *t2 = ft->evalspace2->vals;
670 *t3 = ft->evalspace3->vals;
671 }
672
673 /***********************************************************//**
674 Evaluate a function train up to, and including core k
675
676 \param[in] ft - function train
677 \param[in] core - last core to evaluate
678 \param[in] x - location at which to evaluate
679 \param[in,out] out - evaluation, with enough space
680 \param[in,out] out_dim - size of out
681
682 \note NOT TESTED
683 ***************************************************************/
function_train_eval_up_to_core(struct FunctionTrain * ft,size_t core,const double * x,double * out,size_t * out_dim)684 void function_train_eval_up_to_core(struct FunctionTrain * ft, size_t core, const double * x, double * out, size_t *out_dim)
685 {
686 size_t dim = ft->dim;
687 assert(ft->ranks[0] == 1);
688 assert(ft->ranks[dim] == 1);
689
690 size_t maxrank,ii = 0;
691 double *t1=NULL, *t2=NULL, *t3=NULL;
692 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
693
694 generic_function_1darray_eval2(
695 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
696 ft->cores[ii]->funcs, x[ii],t1);
697
698 int onsol = 1;
699 for (ii = 1; ii < core; ii++){
700 if (ii%2 == 1){
701 function_train_eval_next_core(ft,ii,1,x+ii,0,t1,0,t2,0,t3,0); // zeros because only 1 data point
702 onsol = 2;
703 }
704 else{
705 function_train_eval_next_core(ft,ii,1,x+ii,0,t3,0,t2,0,t1,0);
706 onsol=1;
707 }
708 }
709
710 *out_dim = ft->cores[core]->ncols;
711 if (onsol == 2){
712 memmove(out,t3,*out_dim*sizeof(double));
713 }
714 else if ( onsol == 1){
715 memmove(out,t1,*out_dim*sizeof(double));
716 }
717 else{
718 fprintf(stderr,"Weird error in function_train_val\n");
719 exit(1);
720 }
721 }
722
723 /***********************************************************//**
724 Evaluate a function train
725
726 \param[in] ft - function train
727 \param[in] x - location at which to evaluate
728
729 \return val - value of the function train
730 ***************************************************************/
function_train_eval(struct FunctionTrain * ft,const double * x)731 double function_train_eval(struct FunctionTrain * ft, const double * x)
732 {
733
734 size_t dim = ft->dim;
735 assert(ft->ranks[0] == 1);
736 assert(ft->ranks[dim] == 1);
737
738 size_t maxrank;
739 double *t1=NULL, *t2=NULL, *t3=NULL;
740 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
741
742 size_t ii = 0;
743 generic_function_1darray_eval2(
744 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
745 ft->cores[ii]->funcs, x[ii],t1);
746 int onsol = 1;
747 /* printf("Start running eval: "); */
748 /* dprint(ft->ranks[ii+1],t1); */
749 for (ii = 1; ii < dim; ii++){
750 if (ii%2 == 1){
751 function_train_eval_next_core(ft,ii,1,x+ii,0,t1,0,t2,0,t3,0); // zero because 1 data pint
752 onsol = 2;
753 /* printf("new eval = \n"); */
754 /* dprint2d_col(ft->ranks[ii],ft->ranks[ii+1],t2); */
755 /* printf("\t running eval: "); */
756 /* dprint(ft->ranks[ii+1],t3); */
757 }
758 else{
759 function_train_eval_next_core(ft,ii,1,x+ii,0,t3,0,t2,0,t1,0);
760 onsol=1;
761 /* printf("new eval = \n"); */
762 /* dprint2d_col(ft->ranks[ii],ft->ranks[ii+1],t2); */
763 /* printf("\t running eval: "); */
764 /* dprint(ft->ranks[ii+1],t1); */
765 }
766
767 /* if (ii == dim-1){ */
768 /* if (onsol == 2){ */
769 /* printf("pre_eval = "); dprint(2,t1); */
770 /* } */
771 /* else{ */
772 /* printf("pre_eval = "); dprint(2,t3); */
773 /* } */
774 /* printf("last_core_eval = "); dprint(2,t2); */
775 /* } */
776 }
777
778 double out;// = 0.123456789;
779 if (onsol == 2){
780 out = t3[0];
781 // free(t3); t3 = NULL;
782 }
783 else if ( onsol == 1){
784 out = t1[0];
785 // free(t1); t1 = NULL;
786 }
787 else{
788 fprintf(stderr,"Weird error in function_train_val\n");
789 exit(1);
790 }
791
792 return out;
793 }
794
795
796 /***********************************************************//**
797 Get number of parameters of a particular univariate function
798
799 \param[in] ft - functoin train
800 \param[in] core - which core parameters to count
801 \param[in] ii - row
802 \param[in] jj - column
803
804 \returns number of parameters in the core
805 ***************************************************************/
function_train_func_get_nparams(const struct FunctionTrain * ft,size_t core,size_t ii,size_t jj)806 size_t function_train_func_get_nparams(const struct FunctionTrain * ft,
807 size_t core, size_t ii, size_t jj)
808 {
809 size_t nparam =
810 qmarray_func_get_nparams(ft->cores[core],ii,jj);
811 return nparam;
812 }
813
814 /***********************************************************//**
815 Get number of parameters describing a core
816
817 \param[in] ft - functoin train
818 \param[in] core - which core parameters to count
819 \param[in,out] maxparamfunc - number of parameters in the function
820 with most parameters within the core
821
822 \returns number of parameters in the core
823 ***************************************************************/
function_train_core_get_nparams(const struct FunctionTrain * ft,size_t core,size_t * maxparamfunc)824 size_t function_train_core_get_nparams(const struct FunctionTrain * ft,
825 size_t core,
826 size_t * maxparamfunc)
827 {
828 size_t ncoreparams = qmarray_get_nparams(ft->cores[core],maxparamfunc);
829 return ncoreparams;
830 }
831
832 /***********************************************************//**
833 Get number of parameters describing a ft
834
835 \param[in] ft - function train
836
837 \returns number of parameters in the FT
838 ***************************************************************/
function_train_get_nparams(const struct FunctionTrain * ft)839 size_t function_train_get_nparams(const struct FunctionTrain * ft)
840 {
841
842 size_t ncore;
843 size_t maxp;
844 size_t tot = 0;
845 for (size_t ii = 0; ii < ft->dim; ii++){
846 ncore = qmarray_get_nparams(ft->cores[ii],&maxp);
847 tot += ncore;
848 }
849
850 return tot;
851 }
852
853 /***********************************************************//**
854 Get core parameters
855 ***************************************************************/
function_train_core_get_params(const struct FunctionTrain * ft,size_t core,double * param)856 size_t function_train_core_get_params(const struct FunctionTrain * ft,
857 size_t core,
858 double * param)
859 {
860 size_t nparam = qmarray_get_params(ft->cores[core],param);
861 return nparam;
862 }
863
864 /***********************************************************//**
865 Get function paramaters
866 ***************************************************************/
function_train_get_uni(const struct FunctionTrain * ft,size_t core,size_t row,size_t col)867 void * function_train_get_uni(const struct FunctionTrain * ft,
868 size_t core, size_t row, size_t col)
869 {
870 return qmarray_get_func_base(ft->cores[core], row, col);
871 }
872
873 /***********************************************************//**
874 Get function paramaters
875 ***************************************************************/
876 struct GenericFunction *
function_train_get_gfuni(const struct FunctionTrain * ft,size_t core,size_t row,size_t col)877 function_train_get_gfuni(const struct FunctionTrain * ft,
878 size_t core, size_t row, size_t col)
879 {
880 return qmarray_get_func(ft->cores[core], row, col);
881 }
882
883 /***********************************************************//**
884 Get parameters
885 ***************************************************************/
function_train_get_params(const struct FunctionTrain * ft,double * param)886 size_t function_train_get_params(const struct FunctionTrain * ft,
887 double * param)
888 {
889 size_t running = 0;
890 for (size_t ii = 0; ii < ft->dim; ii++){
891 size_t nparam = qmarray_get_params(ft->cores[ii],param+running);
892 running += nparam;
893 }
894 return running;
895 }
896
897
898 /***********************************************************//**
899 Update the parameters of a single univariate function
900 ***************************************************************/
function_train_uni_update_params(struct FunctionTrain * ft,size_t core,size_t row,size_t col,size_t nparam,const double * param)901 void function_train_uni_update_params(struct FunctionTrain * ft,
902 size_t core,
903 size_t row,
904 size_t col,
905 size_t nparam,
906 const double * param)
907 {
908 qmarray_uni_update_params(ft->cores[core], row, col, nparam, param);
909 }
910
911 /***********************************************************//**
912 Update the parameters defining a function train core
913 ***************************************************************/
function_train_core_update_params(struct FunctionTrain * ft,size_t core,size_t nparam,const double * param)914 void function_train_core_update_params(struct FunctionTrain * ft,
915 size_t core,
916 size_t nparam,
917 const double * param)
918 {
919 qmarray_update_params(ft->cores[core],nparam,param);
920 }
921
922 /***********************************************************//**
923 Update the parameters of a function train
924 ***************************************************************/
function_train_update_params(struct FunctionTrain * ft,const double * param)925 size_t function_train_update_params(struct FunctionTrain * ft, const double * param)
926 {
927 size_t running = 0;
928 size_t incore = 0;
929 for (size_t ii = 0; ii < ft->dim; ii++){
930 incore = qmarray_get_nparams(ft->cores[ii],NULL);
931 qmarray_update_params(ft->cores[ii],incore,param+running);
932 running += incore;
933 }
934 return running;
935 }
936
937 /***********************************************************//**
938 Sum the L2 norms of each function of each core
939
940 \param[in] ft - quasimatrix array
941 \param[in] weights - weights for the sum of each column
942 \param[in,out] grad - compute gradient if not null (store gradient here)
943
944 \returns sum of squared L2 norms of each univariate function
945
946 \note
947 If gradient is not null then it adds scaled values of the new gradients to the
948 existing gradient. Gradient is stored function first, row second, column third, core fourth
949 ***************************************************************/
function_train_param_grad_sqnorm(struct FunctionTrain * ft,double * weights,double * grad)950 double function_train_param_grad_sqnorm(struct FunctionTrain * ft, double * weights, double * grad)
951 {
952
953 size_t running = 0;
954 size_t incr;
955 size_t nfmax; // not used
956 double out = 0.0;
957 for (size_t ii = 0; ii < ft->dim; ii++){
958 incr = qmarray_get_nparams(ft->cores[ii], &nfmax);
959 if (grad != NULL){
960 out += qmarray_param_grad_sqnorm(ft->cores[ii],weights[ii],grad+running);
961 }
962 else{
963 out += qmarray_param_grad_sqnorm(ft->cores[ii],weights[ii],NULL);
964 }
965 running+=incr;
966 /* printf("out = %3.15G\n",out); */
967 }
968
969 return out;
970 }
971
972
973 /***********************************************************//**
974 Evaluate a function train at particular indices which
975 should be valid for underlying univariate functions
976 that are nodal basis functions
977
978 \param[in] ft - function train
979 \param[in] ind - location at which to evaluate
980
981 \return val - value of the function train
982 ***************************************************************/
function_train_eval_ind(struct FunctionTrain * ft,const size_t * ind)983 double function_train_eval_ind(struct FunctionTrain * ft, const size_t * ind)
984 {
985
986 size_t dim = ft->dim;
987 assert(ft->ranks[0] == 1);
988 assert(ft->ranks[dim] == 1);
989
990 size_t maxrank;
991 double *t1=NULL, *t2=NULL, *t3=NULL;
992 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
993
994 size_t ii = 0;
995 generic_function_1darray_eval2_ind(
996 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
997 ft->cores[ii]->funcs, ind[ii],t1);
998
999 int onsol = 1;
1000 for (ii = 1; ii < dim; ii++){
1001 generic_function_1darray_eval2_ind(
1002 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1003 ft->cores[ii]->funcs, ind[ii],t2);
1004
1005 if (ii%2 == 1){
1006 // previous times new core
1007 cblas_dgemv(CblasColMajor,CblasTrans,
1008 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1009 t2, ft->ranks[ii],
1010 t1, 1, 0.0, t3, 1);
1011 onsol = 2;
1012 }
1013 else {
1014 // cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, 1,
1015 // ft->ranks[ii+1], ft->ranks[ii], 1.0, t3, 1, t2,
1016 // ft->ranks[ii], 0.0, t1, 1)
1017 cblas_dgemv(CblasColMajor,CblasTrans,
1018 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1019 t2, ft->ranks[ii],
1020 t3,1,0.0,t1,1);
1021 onsol = 1;
1022 }
1023 }
1024
1025 double out;// = 0.123456789;
1026 if (onsol == 2){
1027 out = t3[0];
1028 // free(t3); t3 = NULL;
1029 }
1030 else if ( onsol == 1){
1031 out = t1[0];
1032 // free(t1); t1 = NULL;
1033 }
1034 else{
1035 fprintf(stderr,"Weird error in function_train_val\n");
1036 exit(1);
1037 }
1038
1039 return out;
1040 }
1041
1042 /***********************************************************//**
1043 Evaluate a fiber of the function-train at particular indices
1044 when each fiber is represented in a (fixed) nodal basis
1045
1046 \param[in] ft - function train
1047 \param[in] fixed_ind - indices of fixed_values
1048 \param[in] N - number of evaluations
1049 \param[in] ind - locations of evaluations
1050 \param[in] dim_vary - dimension that is varying
1051 \param[in,out] out - output values
1052
1053 ***************************************************************/
function_train_eval_fiber_ind(struct FunctionTrain * ft,const size_t * fixed_ind,size_t N,const size_t * ind,size_t dim_vary,double * out)1054 void function_train_eval_fiber_ind(struct FunctionTrain * ft,
1055 const size_t * fixed_ind,
1056 size_t N, const size_t * ind,
1057 size_t dim_vary,
1058 double * out)
1059 {
1060
1061 size_t dim = ft->dim;
1062 assert(ft->ranks[0] == 1);
1063 assert(ft->ranks[dim] == 1);
1064
1065 size_t maxrank;
1066 double *t1=NULL, *t2=NULL, *t3=NULL;
1067 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
1068
1069 double * temp_tensor = NULL;
1070 double * temp_mat = NULL;
1071
1072 size_t nrows, ncols;
1073
1074 struct GenericFunction ** arr = ft->cores[dim_vary]->funcs;
1075 nrows = ft->cores[dim_vary]->nrows;
1076 ncols = ft->cores[dim_vary]->ncols;
1077 temp_tensor = calloc_double(nrows * ncols * N);
1078 temp_mat = calloc_double(ncols * N);
1079 for (size_t jj = 0; jj < N; jj++){
1080 generic_function_1darray_eval2_ind(nrows * ncols,arr,ind[jj],
1081 temp_tensor + jj * nrows * ncols);
1082 }
1083
1084 if (dim_vary == 0){
1085 assert (nrows == 1);
1086 memmove(temp_mat, temp_tensor, ncols * N * sizeof(double));
1087 }
1088 else{ // (dim_vary != 0){
1089 arr = ft->cores[0]->funcs;
1090 nrows = ft->cores[0]->nrows;
1091 ncols = ft->cores[0]->ncols;
1092 generic_function_1darray_eval2_ind(nrows*ncols,arr,fixed_ind[0],t1);
1093
1094 /* printf("first core is \n \t"); */
1095 /* dprint(ncols,t1); */
1096 int onsol = 1;
1097 for (size_t ii = 1; ii < dim_vary; ii++){
1098 nrows = ft->cores[ii]->nrows;
1099 ncols = ft->cores[ii]->ncols;
1100 arr = ft->cores[ii]->funcs;
1101
1102 generic_function_1darray_eval2_ind(nrows * ncols, arr,fixed_ind[ii],t2);
1103 /* printf("\t\t t2 = \n"); */
1104 /* dprint2d_col(nrows,ncols,t2); */
1105 if (onsol == 1){
1106 // previous times new core
1107 /* cblas_dgemv(CblasColMajor,CblasTrans, */
1108 /* nrows, ncols, 1.0, t2, nrows, */
1109 /* t1, 1, 0.0, t3, 1); */
1110 cblas_dgemv(CblasColMajor,CblasTrans,
1111 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1112 t2, ft->ranks[ii],
1113 t1, 1, 0.0, t3, 1);
1114
1115 onsol = 2;
1116 /* printf("\t next is t3 \n \t"); */
1117 /* dprint(ncols,t3); */
1118 }
1119 else {
1120 cblas_dgemv(CblasColMajor,CblasTrans,
1121 nrows, ncols, 1.0,
1122 t2, nrows,
1123 t3,1,0.0,t1,1);
1124 onsol = 1;
1125 /* printf("\t next is t1 \n \t"); */
1126 /* dprint(ncols,t1); */
1127 }
1128 }
1129
1130 nrows = ft->cores[dim_vary]->nrows;
1131 ncols = ft->cores[dim_vary]->ncols;
1132 if (onsol == 1){
1133 /* printf("t1 before multiplying by matrix is \n \t"); */
1134 /* dprint(nrows,t1); */
1135 /* dprint2d_col(nrows,ncols,temp_tensor); */
1136 for (size_t jj = 0; jj < N; jj++){
1137 cblas_dgemv(CblasColMajor,CblasTrans,
1138 nrows, ncols, 1.0, temp_tensor + jj*nrows*ncols,
1139 nrows, t1, 1, 0.0, temp_mat + jj*ncols, 1);
1140 }
1141 }
1142 else{
1143 for (size_t jj = 0; jj < N; jj++){
1144 cblas_dgemv(CblasColMajor,CblasTrans,
1145 nrows, ncols, 1.0, temp_tensor + jj*nrows*ncols,
1146 nrows, t3, 1, 0.0, temp_mat + jj*ncols, 1);
1147 }
1148 }
1149 }
1150
1151 // got everything before and including and stored in temp_mat
1152 if (dim_vary == dim-1){ // done b/c ncols = 1
1153 /* printf("we are here\n"); */
1154 memmove(out,temp_mat, N * sizeof(double));
1155 /* printf("out = "); dprint(N,out); */
1156 }
1157 else{ // come from the back
1158
1159 nrows = ft->cores[dim-1]->nrows;
1160 ncols = ft->cores[dim-1]->ncols;
1161 arr = ft->cores[dim-1]->funcs;
1162 generic_function_1darray_eval2_ind(nrows * ncols,arr,fixed_ind[dim-1],t1);
1163 int onsol = 1;
1164 for (size_t ii = dim-2; ii > dim_vary; ii--){
1165 nrows = ft->cores[ii]->nrows;
1166 ncols = ft->cores[ii]->ncols;
1167 arr = ft->cores[ii]->funcs;
1168 generic_function_1darray_eval2_ind(nrows*ncols,arr,fixed_ind[ii],t2);
1169 if (onsol == 1){
1170 // previous times new core
1171 cblas_dgemv(CblasColMajor,CblasNoTrans,
1172 nrows, ncols, 1.0, t2, nrows,
1173 t1, 1, 0.0, t3, 1);
1174 onsol = 2;
1175 }
1176 else {
1177 cblas_dgemv(CblasColMajor,CblasNoTrans,
1178 nrows, ncols, 1.0,
1179 t2, nrows,
1180 t3,1,0.0,t1,1);
1181 onsol = 1;
1182 }
1183 }
1184
1185 // now combine
1186 ncols = ft->cores[dim_vary]->ncols;
1187 if (onsol == 1){
1188 cblas_dgemv(CblasColMajor,CblasTrans,
1189 ncols, N, 1.0, temp_mat,
1190 ncols, t1, 1, 0.0, out, 1);
1191 }
1192 else{
1193 cblas_dgemv(CblasColMajor,CblasTrans,
1194 ncols, N, 1.0, temp_mat,
1195 ncols, t3, 1, 0.0, out, 1);
1196 }
1197
1198 }
1199
1200 free(temp_mat); temp_mat = NULL;
1201 free(temp_tensor); temp_tensor = NULL;
1202 }
1203
1204 /***********************************************************//**
1205 Evaluate a function train with perturbations
1206 in every coordinate
1207
1208 \param[in] ft - function train
1209 \param[in] x - location at which to evaluate
1210 \param[in] pert - perturbed vals (order of storage is )
1211 (dim 1) - +
1212 (dim 2) - +
1213 (dim 3) - +
1214 (dim 4) - +
1215 ... for a total of 2d values
1216 \param[in,out] vals - values at perturbation points
1217
1218 \return val - value of the function train at x
1219 ***************************************************************/
function_train_eval_co_perturb(struct FunctionTrain * ft,const double * x,const double * pert,double * vals)1220 double function_train_eval_co_perturb(struct FunctionTrain * ft,
1221 const double * x, const double * pert,
1222 double * vals)
1223 {
1224
1225 size_t dim = ft->dim;
1226 assert(ft->ranks[0] == 1);
1227 assert(ft->ranks[dim] == 1);
1228
1229
1230 size_t maxrank;
1231 double *t1=NULL, *t2=NULL, *t3=NULL;
1232 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
1233
1234 if (ft->evaldd1 == NULL){
1235 ft->evaldd1 = malloc_dd(dim);
1236 for (size_t ii = 0; ii < dim; ii++){
1237 ft->evaldd1[ii] = calloc_double(maxrank * maxrank);
1238 }
1239 }
1240 if (ft->evaldd3 == NULL){
1241 ft->evaldd3 = malloc_dd(dim);
1242 for (size_t ii = 0; ii < dim; ii++){
1243 ft->evaldd3[ii] = calloc_double(maxrank);
1244 }
1245 }
1246 if (ft->evaldd4 == NULL){
1247 ft->evaldd4 = malloc_dd(dim);
1248 for (size_t ii = 0; ii < dim; ii++){
1249 ft->evaldd4[ii] = calloc_double(maxrank);
1250 }
1251 }
1252
1253 // center cores
1254 // double ** cores_center = malloc_dd(dim);
1255
1256 /* double ** cores_neigh = malloc_dd(2*dim); // cores for neighbors */
1257 /* double ** bprod = malloc_dd(dim); */
1258 /* double ** fprod = malloc_dd(dim); */
1259
1260 double ** cores_center = ft->evaldd1;
1261 double ** bprod = ft->evaldd3;
1262 double ** fprod = ft->evaldd4;
1263
1264 // evaluate the cores for the center
1265 for (size_t ii = 0; ii < dim; ii++){
1266 /* fprod[ii] = calloc_double(ft->ranks[ii+1]); */
1267 /* cores_center[ii] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1268 if (ii == 0){
1269 generic_function_1darray_eval2(
1270 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1271 ft->cores[ii]->funcs, x[ii], fprod[ii]);
1272 memmove(cores_center[ii],fprod[ii],ft->ranks[1]*sizeof(double));
1273 }
1274 else{
1275 generic_function_1darray_eval2(
1276 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1277 ft->cores[ii]->funcs, x[ii], cores_center[ii]);
1278 cblas_dgemv(CblasColMajor,CblasTrans,
1279 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1280 cores_center[ii], ft->ranks[ii],
1281 fprod[ii-1], 1, 0.0, fprod[ii], 1);
1282 }
1283
1284 /* cores_neigh[2*ii] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1285 /* cores_neigh[2*ii+1] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1286
1287 /* generic_function_1darray_eval2( */
1288 /* ft->cores[ii]->nrows * ft->cores[ii]->ncols, */
1289 /* ft->cores[ii]->funcs, pert[2*ii], cores_neigh[2*ii]); */
1290 /* generic_function_1darray_eval2( */
1291 /* ft->cores[ii]->nrows * ft->cores[ii]->ncols, */
1292 /* ft->cores[ii]->funcs, pert[2*ii+1], cores_neigh[2*ii+1]); */
1293 }
1294
1295 for (int ii = dim-1; ii >= 0; ii--){
1296 /* bprod[ii] = calloc_double(ft->ranks[ii]); */
1297 if (ii == (int)dim-1){
1298 generic_function_1darray_eval2(
1299 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1300 ft->cores[ii]->funcs, x[ii], bprod[ii]);
1301 }
1302 else{
1303 cblas_dgemv(CblasColMajor,CblasNoTrans,
1304 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1305 cores_center[ii], ft->ranks[ii],
1306 bprod[ii+1], 1, 0.0, bprod[ii], 1);
1307 }
1308 }
1309
1310 for (size_t ii = 0; ii < dim; ii++){
1311
1312 generic_function_1darray_eval2(
1313 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1314 ft->cores[ii]->funcs, pert[2*ii],t1);
1315 generic_function_1darray_eval2(
1316 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1317 ft->cores[ii]->funcs, pert[2*ii+1], t2);
1318 if (ii == 0){
1319 vals[ii] = cblas_ddot(ft->ranks[1],t1,1,bprod[1],1);
1320 vals[ii+1] = cblas_ddot(ft->ranks[1],t2,1,bprod[1],1);
1321 }
1322 else if (ii == dim-1){
1323 vals[2*ii] = cblas_ddot(ft->ranks[ii],t1,1,
1324 fprod[dim-2],1);
1325 vals[2*ii+1] = cblas_ddot(ft->ranks[ii],t2,
1326 1,fprod[dim-2],1);
1327 }
1328 else{
1329 /* double * temp = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1330 cblas_dgemv(CblasColMajor,CblasNoTrans,
1331 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1332 t1, ft->ranks[ii],
1333 bprod[ii+1], 1, 0.0, t3, 1);
1334 vals[2*ii] = cblas_ddot(ft->ranks[ii],t3,1,
1335 fprod[ii-1],1);
1336
1337 cblas_dgemv(CblasColMajor,CblasNoTrans,
1338 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1339 t2, ft->ranks[ii],
1340 bprod[ii+1], 1, 0.0, t3, 1);
1341
1342 vals[2*ii+1] = cblas_ddot(ft->ranks[ii],t3,1,
1343 fprod[ii-1],1);
1344 /* free(temp); temp = NULL; */
1345 }
1346
1347 }
1348
1349 double out = fprod[dim-1][0];
1350 /* printf("out = %G\n",out); */
1351 /* printf("should equal = %G\n",bprod[0][0]); */
1352 /* free_dd(dim,cores_center); */
1353 /* free_dd(2*dim,cores_neigh); */
1354 /* free_dd(dim,bprod); bprod = NULL; */
1355 /* free_dd(dim,fprod); fprod = NULL; */
1356
1357 return out;
1358 }
1359
1360 /***********************************************************//**
1361 Evaluate a function train with perturbations
1362 in every coordinate for nodal basis function
1363
1364 \param[in] ft - function train
1365 \param[in] x - location at which to evaluate
1366 \param[in] pert - perturbed vals (order of storage is )
1367 (dim 1) - +
1368 (dim 2) - +
1369 (dim 3) - +
1370 (dim 4) - +
1371 ... for a total of 2d values
1372 \param[in,out] vals - values at perturbation points
1373
1374 \return val - value of the function train at x
1375 ***************************************************************/
function_train_eval_co_perturb_ind(struct FunctionTrain * ft,const size_t * x,const size_t * pert,double * vals)1376 double function_train_eval_co_perturb_ind(struct FunctionTrain * ft,
1377 const size_t * x, const size_t * pert,
1378 double * vals)
1379 {
1380
1381 size_t dim = ft->dim;
1382 assert(ft->ranks[0] == 1);
1383 assert(ft->ranks[dim] == 1);
1384
1385 size_t maxrank;
1386 double *t1=NULL, *t2=NULL, *t3=NULL;
1387 function_train_pre_eval(ft,&maxrank,&t1,&t2,&t3,1);
1388 if (ft->evaldd1 == NULL){
1389 ft->evaldd1 = malloc_dd(dim);
1390 for (size_t ii = 0; ii < dim; ii++){
1391 ft->evaldd1[ii] = calloc_double(maxrank * maxrank);
1392 }
1393 }
1394 if (ft->evaldd3 == NULL){
1395 ft->evaldd3 = malloc_dd(dim);
1396 for (size_t ii = 0; ii < dim; ii++){
1397 ft->evaldd3[ii] = calloc_double(maxrank);
1398 }
1399 }
1400 if (ft->evaldd4 == NULL){
1401 ft->evaldd4 = malloc_dd(dim);
1402 for (size_t ii = 0; ii < dim; ii++){
1403 ft->evaldd4[ii] = calloc_double(maxrank);
1404 }
1405 }
1406
1407 // center cores
1408 // double ** cores_center = malloc_dd(dim);
1409
1410 /* double ** cores_neigh = malloc_dd(2*dim); // cores for neighbors */
1411 /* double ** bprod = malloc_dd(dim); */
1412 /* double ** fprod = malloc_dd(dim); */
1413
1414 double ** cores_center = ft->evaldd1;
1415 double ** bprod = ft->evaldd3;
1416 double ** fprod = ft->evaldd4;
1417
1418 // evaluate the cores for the center
1419 for (size_t ii = 0; ii < dim; ii++){
1420 /* fprod[ii] = calloc_double(ft->ranks[ii+1]); */
1421 /* cores_center[ii] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1422 if (ii == 0){
1423 generic_function_1darray_eval2_ind(
1424 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1425 ft->cores[ii]->funcs, x[ii], fprod[ii]);
1426 memmove(cores_center[ii],fprod[ii],ft->ranks[1]*sizeof(double));
1427 }
1428 else{
1429 generic_function_1darray_eval2_ind(
1430 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1431 ft->cores[ii]->funcs, x[ii], cores_center[ii]);
1432 cblas_dgemv(CblasColMajor,CblasTrans,
1433 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1434 cores_center[ii], ft->ranks[ii],
1435 fprod[ii-1], 1, 0.0, fprod[ii], 1);
1436 }
1437
1438 /* cores_neigh[2*ii] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1439 /* cores_neigh[2*ii+1] = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1440
1441 /* generic_function_1darray_eval2( */
1442 /* ft->cores[ii]->nrows * ft->cores[ii]->ncols, */
1443 /* ft->cores[ii]->funcs, pert[2*ii], cores_neigh[2*ii]); */
1444 /* generic_function_1darray_eval2( */
1445 /* ft->cores[ii]->nrows * ft->cores[ii]->ncols, */
1446 /* ft->cores[ii]->funcs, pert[2*ii+1], cores_neigh[2*ii+1]); */
1447 }
1448
1449 for (int ii = dim-1; ii >= 0; ii--){
1450 /* bprod[ii] = calloc_double(ft->ranks[ii]); */
1451 if (ii == (int)dim-1){
1452 generic_function_1darray_eval2_ind(
1453 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1454 ft->cores[ii]->funcs, x[ii], bprod[ii]);
1455 }
1456 else{
1457 cblas_dgemv(CblasColMajor,CblasNoTrans,
1458 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1459 cores_center[ii], ft->ranks[ii],
1460 bprod[ii+1], 1, 0.0, bprod[ii], 1);
1461 }
1462 }
1463
1464 for (size_t ii = 0; ii < dim; ii++){
1465
1466 generic_function_1darray_eval2_ind(
1467 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1468 ft->cores[ii]->funcs, pert[2*ii], t1);
1469 generic_function_1darray_eval2_ind(
1470 ft->cores[ii]->nrows * ft->cores[ii]->ncols,
1471 ft->cores[ii]->funcs, pert[2*ii+1], t2);
1472 if (ii == 0){
1473 vals[ii] = cblas_ddot(ft->ranks[1],t1,1,bprod[1],1);
1474 vals[ii+1] = cblas_ddot(ft->ranks[1],t2,1,bprod[1],1);
1475 }
1476 else if (ii == dim-1){
1477 vals[2*ii] = cblas_ddot(ft->ranks[ii],t1,1,
1478 fprod[dim-2],1);
1479 vals[2*ii+1] = cblas_ddot(ft->ranks[ii],t2,
1480 1,fprod[dim-2],1);
1481 }
1482 else{
1483 /* double * temp = calloc_double(ft->ranks[ii] * ft->ranks[ii+1]); */
1484 cblas_dgemv(CblasColMajor,CblasNoTrans,
1485 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1486 t1, ft->ranks[ii],
1487 bprod[ii+1], 1, 0.0, t3, 1);
1488 vals[2*ii] = cblas_ddot(ft->ranks[ii],t3,1,
1489 fprod[ii-1],1);
1490
1491 cblas_dgemv(CblasColMajor,CblasNoTrans,
1492 ft->ranks[ii], ft->ranks[ii+1], 1.0,
1493 t2, ft->ranks[ii],
1494 bprod[ii+1], 1, 0.0, t3, 1);
1495
1496 vals[2*ii+1] = cblas_ddot(ft->ranks[ii],t3,1,
1497 fprod[ii-1],1);
1498 /* free(temp); temp = NULL; */
1499 }
1500
1501 }
1502
1503 double out = fprod[dim-1][0];
1504 /* printf("out = %G\n",out); */
1505 /* printf("should equal = %G\n",bprod[0][0]); */
1506 /* free_dd(dim,cores_center); */
1507 /* free_dd(2*dim,cores_neigh); */
1508 /* free_dd(dim,bprod); bprod = NULL; */
1509 /* free_dd(dim,fprod); fprod = NULL; */
1510
1511 return out;
1512 }
1513
1514 /********************************************************//**
1515 * Create a function train consisting of pseudo-random orth poly expansion
1516 *
1517 * \param[in] ptype - polynomial type
1518 * \param[in] bds - boundaries
1519 * \param[in] ranks - (dim+1,1) array of ranks
1520 * \param[in] maxorder - maximum order of the polynomial
1521 *
1522 * \return function train
1523 ************************************************************/
1524 struct FunctionTrain *
function_train_poly_randu(enum poly_type ptype,struct BoundingBox * bds,size_t * ranks,size_t maxorder)1525 function_train_poly_randu(enum poly_type ptype,struct BoundingBox * bds,
1526 size_t * ranks, size_t maxorder)
1527 {
1528 size_t dim = bounding_box_get_dim(bds);
1529 double * lb = bounding_box_get_lb(bds);
1530 double * ub = bounding_box_get_ub(bds);
1531 struct FunctionTrain * ft = function_train_alloc(dim);
1532 memmove(ft->ranks,ranks, (dim+1)*sizeof(size_t));
1533
1534 size_t ii;
1535 for (ii = 0; ii < dim; ii++){
1536 ft->cores[ii] =
1537 qmarray_poly_randu(ptype,ranks[ii],ranks[ii+1],maxorder,
1538 lb[ii],ub[ii]);
1539 }
1540 return ft;
1541 }
1542
1543
1544 /***********************************************************//**
1545 Compute a function train representation of \f$ f1(x1)*f2(x2)*...* fd(xd)\f$
1546
1547 \param[in] ftargs - approximation options
1548 \param[in] fw - wrapped functions
1549
1550 \return function train
1551 ***************************************************************/
1552 struct FunctionTrain *
function_train_rankone(struct MultiApproxOpts * ftargs,struct Fwrap * fw)1553 function_train_rankone(struct MultiApproxOpts * ftargs, struct Fwrap * fw)
1554 {
1555 assert (ftargs != NULL);
1556 assert (fw != NULL);
1557 size_t dim = multi_approx_opts_get_dim(ftargs);
1558 struct OneApproxOpts * aopt;
1559
1560 struct FunctionTrain * ft = function_train_alloc(dim);
1561 size_t onDim;
1562 for (onDim = 0; onDim < dim; onDim++){
1563
1564 aopt = multi_approx_opts_get_aopts(ftargs,onDim);
1565
1566 fwrap_set_which_eval(fw,onDim);
1567
1568 ft->ranks[onDim] = 1;
1569 ft->cores[onDim] = qmarray_alloc(1,1);
1570 ft->cores[onDim]->funcs[0] =
1571 generic_function_approximate1d(aopt->fc,aopt->aopts,fw);
1572 }
1573 ft->ranks[dim] = 1;
1574 return ft;
1575 }
1576
1577 /***********************************************************//**
1578 Compute a function train representation of \f$ f1(x_1) + f2(x_2) + ... + fd(x_d) \f$
1579
1580 \param[in] ftargs - approximation options
1581 \param[in] fw - wrapped functions
1582
1583 \return ft - function train
1584 ***************************************************************/
1585 struct FunctionTrain *
function_train_initsum(struct MultiApproxOpts * ftargs,struct Fwrap * fw)1586 function_train_initsum(struct MultiApproxOpts * ftargs, struct Fwrap * fw)
1587 {
1588 assert (ftargs != NULL);
1589 assert (fw != NULL);
1590 size_t dim = multi_approx_opts_get_dim(ftargs);
1591 struct OneApproxOpts * aopt;
1592
1593 struct FunctionTrain * ft = function_train_alloc(dim);
1594
1595 size_t onDim = 0;
1596 fwrap_set_which_eval(fw,onDim);
1597 aopt = multi_approx_opts_get_aopts(ftargs, onDim);
1598 ft->ranks[onDim] = 1;
1599 if (dim == 1){
1600 ft->cores[onDim] = qmarray_alloc(1,1);
1601 }
1602 else{
1603 ft->cores[onDim] = qmarray_alloc(1,2);
1604 }
1605
1606 ft->cores[onDim]->funcs[0] =
1607 generic_function_approximate1d(aopt->fc,aopt->aopts,fw);
1608 ft->cores[onDim]->funcs[1] = generic_function_constant(1.0, aopt->fc,
1609 aopt->aopts);
1610
1611 if (dim > 1){
1612 for (onDim = 1; onDim < dim-1; onDim++){
1613
1614 fwrap_set_which_eval(fw,onDim);
1615 aopt = multi_approx_opts_get_aopts(ftargs, onDim);
1616 ft->ranks[onDim] = 2;
1617 ft->cores[onDim] = qmarray_alloc(2,2);
1618
1619 ft->cores[onDim]->funcs[0] =
1620 generic_function_constant(1.0,aopt->fc,aopt->aopts);
1621
1622 ft->cores[onDim]->funcs[1] =
1623 generic_function_approximate1d(aopt->fc,aopt->aopts,fw);
1624
1625 ft->cores[onDim]->funcs[2] =
1626 generic_function_constant(0.0, aopt->fc, aopt->aopts);
1627
1628 ft->cores[onDim]->funcs[3] =
1629 generic_function_constant(1.0, aopt->fc, aopt->aopts);
1630 }
1631
1632 onDim = dim-1;
1633 fwrap_set_which_eval(fw,onDim);
1634 aopt = multi_approx_opts_get_aopts(ftargs, onDim);
1635
1636 ft->ranks[onDim] = 2;
1637 ft->ranks[onDim+1] = 1;
1638 ft->cores[onDim] = qmarray_alloc(2,1);
1639
1640 ft->cores[onDim]->funcs[0] =
1641 generic_function_constant(1.0, aopt->fc, aopt->aopts);
1642
1643 ft->cores[onDim]->funcs[1] =
1644 generic_function_approximate1d(aopt->fc,aopt->aopts,fw);
1645 }
1646
1647 return ft;
1648 }
1649
1650 /***********************************************************//**
1651 Create a zeros function train
1652
1653 \param[in] aopts - approximation options
1654 \param[in] ranks - ranks
1655
1656 \return function train
1657
1658 \note
1659 Puts the constant into the first core
1660 ***************************************************************/
1661 struct FunctionTrain *
function_train_zeros(struct MultiApproxOpts * aopts,const size_t * ranks)1662 function_train_zeros(struct MultiApproxOpts * aopts, const size_t * ranks)
1663 {
1664 assert (aopts != NULL);
1665
1666 size_t dim = multi_approx_opts_get_dim(aopts);
1667 struct FunctionTrain * ft = function_train_alloc(dim);
1668 memmove(ft->ranks,ranks,(dim+1) * sizeof(size_t));
1669
1670 struct OneApproxOpts * oneopt = NULL;
1671 for (size_t onDim = 0; onDim < dim; onDim++){
1672 oneopt = multi_approx_opts_get_aopts(aopts,onDim);
1673 ft->cores[onDim] = qmarray_alloc(ranks[onDim],ranks[onDim+1]);
1674 for (size_t jj = 0; jj < ranks[onDim] * ranks[onDim+1]; jj++){
1675 ft->cores[onDim]->funcs[jj] =
1676 generic_function_zero(oneopt->fc,oneopt->aopts,1);
1677 }
1678 }
1679 return ft;
1680 }
1681
1682 /***********************************************************//**
1683 Compute a function train representation of \f$ a \f$
1684
1685 \param[in] a - value
1686 \param[in] aopts - approximation options
1687
1688 \return function train
1689
1690 \note
1691 Puts the constant into the first core
1692 ***************************************************************/
1693 struct FunctionTrain *
function_train_constant(double a,struct MultiApproxOpts * aopts)1694 function_train_constant(double a, struct MultiApproxOpts * aopts)
1695 {
1696 assert (aopts != NULL);
1697 size_t dim = multi_approx_opts_get_dim(aopts);
1698 struct FunctionTrain * ft = function_train_alloc(dim);
1699
1700 size_t onDim = 0;
1701
1702
1703 double sign = 1;
1704 if (a < 0){
1705 sign = -1;
1706 }
1707 double apow = pow(fabs(a),1.0/dim);
1708 /* printf("a = %3.15G apow = %3.15G\n",a,apow); */
1709 struct OneApproxOpts * oneopt = multi_approx_opts_get_aopts(aopts,onDim);
1710 ft->ranks[onDim] = 1;
1711 ft->cores[onDim] = qmarray_alloc(1,1);
1712 ft->cores[onDim]->funcs[0] =
1713 generic_function_constant(sign*apow,oneopt->fc,oneopt->aopts);
1714
1715 for (onDim = 1; onDim < dim; onDim++){
1716 oneopt = multi_approx_opts_get_aopts(aopts,onDim);
1717 ft->ranks[onDim] = 1;
1718 ft->cores[onDim] = qmarray_alloc(1,1);
1719 ft->cores[onDim]->funcs[0] =
1720 generic_function_constant(apow,oneopt->fc,oneopt->aopts);
1721 }
1722 ft->ranks[dim] = 1;
1723
1724 return ft;
1725 }
1726
1727 /***********************************************************//**
1728 Compute a tensor train representation of
1729
1730 \f$ (x_1c_1+a_1) + (x_2c_2+a_2) + .... + (x_dc_d+a_d) \f$
1731
1732 \param[in] c - slope of the function in each dimension (ldc x dim)
1733 \param[in] ldc - stride of c
1734 \param[in] a - value
1735 \param[in] lda - stride of a
1736 \param[in] aopts - approximation options
1737
1738 \return function train
1739 ***************************************************************/
1740 struct FunctionTrain *
function_train_linear(const double * c,size_t ldc,const double * a,size_t lda,struct MultiApproxOpts * aopts)1741 function_train_linear(const double * c, size_t ldc,
1742 const double * a, size_t lda,
1743 struct MultiApproxOpts * aopts)
1744 {
1745
1746 assert (aopts != NULL);
1747 size_t dim = multi_approx_opts_get_dim(aopts);
1748 struct FunctionTrain * ft = function_train_alloc(dim);
1749
1750 size_t onDim = 0;
1751 ft->ranks[onDim] = 1;
1752 if (dim == 1){
1753 ft->cores[onDim] = qmarray_alloc(1,1);
1754 }
1755 else{
1756 ft->cores[onDim] = qmarray_alloc(1,2);
1757 }
1758
1759 struct OneApproxOpts * o = NULL;
1760 o = multi_approx_opts_get_aopts(aopts,onDim);
1761
1762 ft->cores[onDim]->funcs[0] =
1763 generic_function_linear(c[onDim*ldc], a[onDim*lda],o->fc,o->aopts);
1764
1765 if (dim > 1){
1766
1767 ft->cores[onDim]->funcs[1] =
1768 generic_function_constant(1.0, o->fc,o->aopts);
1769
1770 for (onDim = 1; onDim < dim-1; onDim++){
1771 o = multi_approx_opts_get_aopts(aopts,onDim);
1772 ft->ranks[onDim] = 2;
1773 ft->cores[onDim] = qmarray_alloc(2,2);
1774
1775 ft->cores[onDim]->funcs[0] =
1776 generic_function_constant(1.0, o->fc,o->aopts);
1777
1778 ft->cores[onDim]->funcs[1] =
1779 generic_function_linear(c[onDim*ldc], a[onDim*lda],
1780 o->fc,o->aopts);
1781
1782 ft->cores[onDim]->funcs[2] =
1783 generic_function_constant(0.0, o->fc, o->aopts);
1784
1785 ft->cores[onDim]->funcs[3] =
1786 generic_function_constant(1.0, o->fc, o->aopts);
1787 }
1788
1789 onDim = dim-1;
1790 o = multi_approx_opts_get_aopts(aopts,onDim);
1791 ft->ranks[onDim] = 2;
1792 ft->ranks[onDim+1] = 1;
1793 ft->cores[onDim] = qmarray_alloc(2,1);
1794
1795 ft->cores[onDim]->funcs[0] =
1796 generic_function_constant(1.0, o->fc, o->aopts);
1797
1798 ft->cores[onDim]->funcs[1] =
1799 generic_function_linear(c[onDim*ldc], a[onDim*lda],
1800 o->fc, o->aopts);
1801 }
1802
1803 return ft;
1804 }
1805
1806 /***********************************************************//**
1807 Update a function-train representation of
1808
1809 \f$ (x_1c_1+a_1) + (x_2c_2+a_2) + .... + (x_dc_d+a_d) \f$
1810
1811 \param[in] ft - existing linear ft, formed by function_train_linear
1812 \param[in] c - slope of the function in each dimension (ldc x dim)
1813 \param[in] ldc - stride of c
1814 \param[in] a - value
1815 \param[in] lda - stride of a
1816
1817 \return 0 if successfull, 1 otherwise
1818
1819 \note This assumes that we are updating an FT created (and not modified)
1820 with function_train_linear.
1821 ***************************************************************/
1822 int
function_train_linear_update(struct FunctionTrain * ft,const double * c,size_t ldc,const double * a,size_t lda)1823 function_train_linear_update(
1824 struct FunctionTrain * ft,
1825 const double * c, size_t ldc,
1826 const double * a, size_t lda)
1827 {
1828
1829 size_t dim = ft->dim;
1830 size_t onDim = 0;
1831 generic_function_linear_update(ft->cores[onDim]->funcs[0],
1832 c[onDim*ldc], a[onDim*lda]);
1833
1834
1835 if (dim > 1){
1836 for (onDim = 1; onDim < dim-1; onDim++){
1837 generic_function_linear_update(ft->cores[onDim]->funcs[1],
1838 c[onDim*ldc],
1839 a[onDim*lda]);
1840 }
1841
1842 onDim = dim-1;
1843 generic_function_linear_update(ft->cores[onDim]->funcs[1],
1844 c[onDim*ldc],
1845 a[onDim*lda]);
1846 }
1847
1848 return 0;
1849 }
1850
1851 /***********************************************************//**
1852 Compute a function train representation of \f$ \prod_{i=1}^d(a_ix_i) \f$
1853
1854 \param[in] coeffs - coefficients
1855 \param[in] aopts - approximation options
1856
1857 \returns function train
1858
1859 ***************************************************************/
1860 struct FunctionTrain *
function_train_rankone_prod(const double * coeffs,struct MultiApproxOpts * aopts)1861 function_train_rankone_prod(const double * coeffs,
1862 struct MultiApproxOpts * aopts)
1863 {
1864 assert (aopts != NULL);
1865 size_t dim = multi_approx_opts_get_dim(aopts);
1866 assert (dim > 1); //
1867
1868 struct FunctionTrain * ft = function_train_alloc(dim);
1869 struct OneApproxOpts * o = NULL;
1870
1871 size_t onDim = 0;
1872 ft->ranks[onDim] = 1;
1873 ft->cores[onDim] = qmarray_alloc(1,1);
1874
1875 o = multi_approx_opts_get_aopts(aopts,onDim);
1876
1877 // should be quadratic
1878 ft->cores[onDim]->funcs[0] =
1879 generic_function_linear(coeffs[onDim],0.0,o->fc,o->aopts);
1880
1881 for (onDim = 1; onDim < dim-1; onDim++){
1882 //printf("on dimension (%zu/%zu)\n",onDim+1,dim);
1883
1884 o = multi_approx_opts_get_aopts(aopts,onDim);
1885 ft->ranks[onDim] = 1;
1886 ft->cores[onDim] = qmarray_alloc(1,1);
1887 ft->cores[onDim]->funcs[0] =
1888 generic_function_linear(coeffs[onDim],0.0,o->fc,o->aopts);
1889 }
1890
1891 onDim = dim-1;
1892 o = multi_approx_opts_get_aopts(aopts,onDim);
1893 ft->ranks[onDim] = 1;
1894 ft->cores[onDim] = qmarray_alloc(1,1);
1895
1896 ft->cores[onDim]->funcs[0] =
1897 generic_function_linear(coeffs[onDim],0.0,o->fc,o->aopts);
1898
1899 ft->ranks[dim] = 1;
1900 return ft;
1901 }
1902
1903 /***********************************************************//**
1904 Compute a function train representation of \f$ (x-m)^T Q (x-m) \f$
1905
1906 \param[in] coeffs - Q matrix
1907 \param[in] m - m (dim,)
1908 \param[in] aopts - approximation options
1909
1910 \returns function train
1911
1912 \note
1913 Could be more efficient with a better distribution of ranks
1914 ***************************************************************/
1915 struct FunctionTrain *
function_train_quadratic(const double * coeffs,const double * m,struct MultiApproxOpts * aopts)1916 function_train_quadratic(const double * coeffs,
1917 const double * m, struct MultiApproxOpts * aopts)
1918 {
1919 assert (aopts != NULL);
1920 size_t dim = multi_approx_opts_get_dim(aopts);
1921 assert (dim > 1); //
1922
1923 struct FunctionTrain * ft = function_train_alloc(dim);
1924 struct OneApproxOpts * o = NULL;
1925
1926 double temp;
1927 size_t kk,ll;
1928 size_t onDim = 0;
1929
1930 ft->ranks[onDim] = 1;
1931 ft->cores[onDim] = qmarray_alloc(1,dim+1);
1932
1933 o = multi_approx_opts_get_aopts(aopts,onDim);
1934
1935 // should be quadratic
1936 ft->cores[onDim]->funcs[0] =
1937 generic_function_quadratic(coeffs[onDim*dim+onDim], m[onDim],
1938 o->fc,o->aopts);
1939
1940 for (kk = 1; kk < dim; kk++)
1941 {
1942 //printf("on dimension 1 (%zu/%zu)\n",kk,dim);
1943 temp = coeffs[onDim*dim+kk] + coeffs[kk*dim + onDim];
1944 ft->cores[onDim]->funcs[kk] =
1945 generic_function_linear(temp,-temp*m[onDim], o->fc,o->aopts);
1946 }
1947
1948 ft->cores[onDim]->funcs[dim] = generic_function_constant(1.0,o->fc,o->aopts);
1949
1950 for (onDim = 1; onDim < dim-1; onDim++){
1951 //printf("on dimension (%zu/%zu)\n",onDim+1,dim);
1952
1953 o = multi_approx_opts_get_aopts(aopts,onDim);
1954 ft->ranks[onDim] = dim-onDim+2;
1955 ft->cores[onDim] = qmarray_alloc(dim-onDim+2,dim-onDim+1);
1956
1957 for (kk = 0; kk < (dim-onDim+1); kk++){ // loop over columns
1958 for (ll = 0; ll < (dim - onDim + 2); ll++){ //rows;
1959 if ( (ll == 0) && (kk == 0)){ // upper left corner
1960 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1961 generic_function_constant(1.0, o->fc,o->aopts);
1962 }
1963 else if ( (kk == 0) && (ll == 1) ){ // first element of lower diagonal
1964 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1965 generic_function_linear(1.0,-m[onDim], o->fc,
1966 o->aopts);
1967 }
1968 else if (ll == (kk+1)){ // lower diagonal
1969 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1970 generic_function_constant(1.0, o->fc,o->aopts);
1971 }
1972 else if ( (ll == (dim-onDim+1)) && (kk == 0)){ //lower left corner
1973 //quadratic
1974 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1975 generic_function_quadratic(
1976 coeffs[onDim*dim+onDim], m[onDim],
1977 o->fc, o->aopts);
1978 }
1979 else if ( ll == (dim-onDim+1) ){ // rest of bottom row
1980 temp = coeffs[onDim*dim+onDim+kk] +
1981 coeffs[(onDim+kk)*dim + onDim];
1982
1983 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1984 generic_function_linear(temp,-temp*m[onDim],o->fc,o->aopts);
1985 }
1986 else{ // zeros
1987 ft->cores[onDim]->funcs[kk*ft->cores[onDim]->nrows+ll] =
1988 generic_function_constant(0.0,o->fc,o->aopts);
1989 }
1990 }
1991 }
1992 }
1993
1994 onDim = dim-1;
1995 o = multi_approx_opts_get_aopts(aopts,onDim);
1996 ft->ranks[onDim] = dim-onDim+2;
1997 ft->cores[onDim] = qmarray_alloc(dim-onDim+2,1);
1998
1999 ft->cores[onDim]->funcs[0] = generic_function_constant(1.0, o->fc,o->aopts);
2000 ft->cores[onDim]->funcs[1] = generic_function_linear(1.0, -m[onDim],
2001 o->fc,o->aopts);
2002 ft->cores[onDim]->funcs[2] = generic_function_quadratic(
2003 coeffs[onDim*dim+onDim], m[onDim],
2004 o->fc, o->aopts);
2005 ft->ranks[dim] = 1;
2006 return ft;
2007 }
2008
2009 /***********************************************************//**
2010 Compute a tensor train representation of \f[ (x_1-m_1)^2c_1 + (x_2-m_2)^2c_2 + .... + (x_d-m_d)^2c_d \f]
2011
2012 \param[in] coeffs - coefficients for each dimension
2013 \param[in] m - offset in each dimension
2014 \param[in] aopts - approximation arguments
2015
2016 \return a function train
2017 ***************************************************************/
2018 struct FunctionTrain *
function_train_quadratic_aligned(const double * coeffs,const double * m,struct MultiApproxOpts * aopts)2019 function_train_quadratic_aligned(const double * coeffs, const double * m,
2020 struct MultiApproxOpts * aopts)
2021 {
2022
2023 assert (aopts != NULL);
2024 size_t dim = multi_approx_opts_get_dim(aopts);
2025
2026
2027 struct FunctionTrain * ft = function_train_alloc(dim);
2028 struct OneApproxOpts * o = NULL;
2029
2030 size_t onDim = 0;
2031 ft->ranks[onDim] = 1;
2032 if (dim == 1){
2033 ft->cores[onDim] = qmarray_alloc(1,1);
2034 }
2035 else{
2036 ft->cores[onDim] = qmarray_alloc(1,2);
2037 }
2038
2039 if (dim > 1){
2040 o = multi_approx_opts_get_aopts(aopts,onDim);
2041
2042 ft->cores[onDim]->funcs[0] =
2043 generic_function_quadratic(coeffs[onDim], m[onDim],o->fc,o->aopts);
2044
2045 ft->cores[onDim]->funcs[1] =
2046 generic_function_constant(1.0, o->fc, o->aopts);
2047
2048 for (onDim = 1; onDim < dim-1; onDim++){
2049 o = multi_approx_opts_get_aopts(aopts,onDim);
2050 ft->ranks[onDim] = 2;
2051 ft->cores[onDim] = qmarray_alloc(2,2);
2052
2053 ft->cores[onDim]->funcs[0] = generic_function_constant(1.0,
2054 o->fc,o->aopts);
2055
2056 ft->cores[onDim]->funcs[1] =
2057 generic_function_quadratic(coeffs[onDim], m[onDim],
2058 o->fc, o->aopts);
2059
2060
2061 ft->cores[onDim]->funcs[2] = generic_function_constant(0.0, o->fc,
2062 o->aopts);
2063
2064 ft->cores[onDim]->funcs[3] = generic_function_constant(1.0, o->fc,
2065 o->aopts);
2066 }
2067
2068 onDim = dim-1;
2069
2070 ft->ranks[onDim] = 2;
2071 ft->ranks[onDim+1] = 1;
2072 ft->cores[onDim] = qmarray_alloc(2,1);
2073
2074 ft->cores[onDim]->funcs[0] = generic_function_constant(1.0, o->fc,
2075 o->aopts);
2076
2077 ft->cores[onDim]->funcs[1] =
2078 generic_function_quadratic(coeffs[onDim], m[onDim],o->fc,o->aopts);
2079 }
2080 else{
2081 o = multi_approx_opts_get_aopts(aopts,onDim);
2082 ft->cores[onDim]->funcs[0] =
2083 generic_function_quadratic(coeffs[onDim], m[onDim],o->fc,o->aopts);
2084 }
2085
2086 return ft;
2087 }
2088
2089 /********************************************************//**
2090 Interpolate a function-train onto a particular grid forming
2091 another function_train with a nodela basis
2092
2093 \param[in] fta - Function train array to evaluate
2094 \param[in] N - number of nodes in each dimension
2095 \param[in] x - nodes in each dimension
2096
2097 \return new function train
2098 ***********************************************************/
2099 struct FunctionTrain *
function_train_create_nodal(const struct FunctionTrain * fta,size_t * N,double ** x)2100 function_train_create_nodal(const struct FunctionTrain * fta, size_t * N, double ** x)
2101 {
2102 struct FunctionTrain * newft = function_train_alloc(fta->dim);
2103 memmove(newft->ranks,fta->ranks, (fta->dim+1)*sizeof(size_t));
2104 for (size_t ii = 0; ii < newft->dim; ii ++){
2105 newft->cores[ii] = qmarray_create_nodal(fta->cores[ii],N[ii],x[ii]);
2106 }
2107 return newft;
2108 }
2109
2110
2111 /********************************************************//**
2112 Integrate a function in function train format
2113
2114 \param[in] ft - Function train 1
2115
2116 \return Integral \f$ \int f(x) dx \f$
2117 ***********************************************************/
2118 double
function_train_integrate(const struct FunctionTrain * ft)2119 function_train_integrate(const struct FunctionTrain * ft)
2120 {
2121 size_t dim = ft->dim;
2122 size_t ii = 0;
2123 double * t1 = qmarray_integrate(ft->cores[ii]);
2124
2125 double * t2 = NULL;
2126 double * t3 = NULL;
2127 for (ii = 1; ii < dim; ii++){
2128 t2 = qmarray_integrate(ft->cores[ii]);
2129 if (ii%2 == 1){
2130 // previous times new core
2131 t3 = calloc_double(ft->ranks[ii+1]);
2132 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, 1,
2133 ft->ranks[ii+1], ft->ranks[ii], 1.0, t1, 1, t2,
2134 ft->ranks[ii], 0.0, t3, 1);
2135 free(t2); t2 = NULL;
2136 free(t1); t1 = NULL;
2137 }
2138 else {
2139 t1 = calloc_double(ft->ranks[ii+1]);
2140 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, 1,
2141 ft->ranks[ii+1], ft->ranks[ii], 1.0, t3, 1, t2,
2142 ft->ranks[ii], 0.0, t1, 1);
2143 free(t2); t2 = NULL;
2144 free(t3); t3 = NULL;
2145 }
2146 }
2147
2148 double out = 0.123456789;
2149 if (t1 == NULL){
2150 out = t3[0];
2151 free(t3); t3 = NULL;
2152 }
2153 else if ( t3 == NULL){
2154 out = t1[0];
2155 free(t1); t1 = NULL;
2156 }
2157
2158 return out;
2159 }
2160
2161 /********************************************************//**
2162 Integrate a function in function train format with respect
2163 to an appropriately weight
2164
2165 \param[in] ft - Function train 1
2166
2167 \return Integral \f$ \int f(x) w(x) dx \f$
2168
2169 \note w(x) depends on underlying parameterization
2170 for example, it is 1/2 for legendre (and default for others),
2171 gauss for hermite,etc
2172 ***********************************************************/
2173 double
function_train_integrate_weighted(const struct FunctionTrain * ft)2174 function_train_integrate_weighted(const struct FunctionTrain * ft)
2175 {
2176 size_t dim = ft->dim;
2177 size_t ii = 0;
2178 double * t1 = qmarray_integrate_weighted(ft->cores[ii]);
2179
2180 double * t2 = NULL;
2181 double * t3 = NULL;
2182 for (ii = 1; ii < dim; ii++){
2183 t2 = qmarray_integrate_weighted(ft->cores[ii]);
2184 if (ii%2 == 1){
2185 // previous times new core
2186 t3 = calloc_double(ft->ranks[ii+1]);
2187 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, 1,
2188 ft->ranks[ii+1], ft->ranks[ii], 1.0, t1, 1, t2,
2189 ft->ranks[ii], 0.0, t3, 1);
2190 free(t2); t2 = NULL;
2191 free(t1); t1 = NULL;
2192 }
2193 else {
2194 t1 = calloc_double(ft->ranks[ii+1]);
2195 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, 1,
2196 ft->ranks[ii+1], ft->ranks[ii], 1.0, t3, 1, t2,
2197 ft->ranks[ii], 0.0, t1, 1);
2198 free(t2); t2 = NULL;
2199 free(t3); t3 = NULL;
2200 }
2201 }
2202
2203 double out = 0.123456789;
2204 if (t1 == NULL){
2205 out = t3[0];
2206 free(t3); t3 = NULL;
2207 }
2208 else if ( t3 == NULL){
2209 out = t1[0];
2210 free(t1); t1 = NULL;
2211 }
2212
2213 return out;
2214 }
2215
2216
2217 /********************************************************//**
2218 Contract the function along certain dimensions with respect
2219 to an appropriate reight
2220
2221 \param[in] ft - Function train 1
2222 \param[in] ndim_contract - number of dimensions to contract
2223 \param[in] dims_contract - dimensions over which to contract (increasing order)
2224
2225 \return Evaluates \f$ g(x_{\neq c}) = \int f(x_c) w(x_c) dx_c \f$
2226
2227 \note w(x) depends on underlying parameterization
2228 for example, it is 1/2 for legendre (and default for others),
2229 gauss for hermite,etc
2230 ***********************************************************/
2231 struct FunctionTrain *
function_train_integrate_weighted_subset(const struct FunctionTrain * ft,size_t ndim_contract,size_t * dims_contract)2232 function_train_integrate_weighted_subset(
2233 const struct FunctionTrain * ft,
2234 size_t ndim_contract,
2235 size_t * dims_contract)
2236 {
2237
2238
2239 size_t newdim = ft->dim - ndim_contract;
2240 assert (newdim > 0);
2241 /* printf("newdim = %zu\n",newdim); */
2242 struct FunctionTrain * newft = function_train_alloc(newdim);
2243 newft->ranks[0] = 1;
2244 newft->ranks[newdim] = 1;
2245
2246 double * left_mult = NULL;
2247 size_t ncols = 0;
2248
2249 size_t on_contract = 0;
2250 size_t ii = 0;
2251
2252 while ((dims_contract[on_contract] == ii ) && (on_contract < ndim_contract)){
2253 double * mat = qmarray_integrate_weighted(ft->cores[ii]);
2254
2255 if (left_mult == NULL){
2256 ncols = ft->cores[ii]->ncols;
2257 left_mult = calloc_double(ncols);
2258 memmove(left_mult,mat,ncols * sizeof(double));
2259 }
2260 else{
2261 double * new_left_mult = calloc_double(ft->cores[ii]->ncols);
2262 c3linalg_multiple_vec_mat(1,ft->cores[ii]->nrows,ft->cores[ii]->ncols,
2263 left_mult,ncols,
2264 mat,ft->cores[ii]->nrows*ft->cores[ii]->ncols,
2265 new_left_mult,ft->cores[ii]->ncols);
2266 free(left_mult); left_mult = NULL;
2267
2268 ncols = ft->cores[ii]->ncols;
2269 left_mult = calloc_double(ncols);
2270 memmove(left_mult,new_left_mult,ncols*sizeof(double));
2271 free(new_left_mult); new_left_mult = NULL;
2272 }
2273
2274 free(mat); mat = NULL;
2275 ii++;
2276 on_contract++;
2277
2278 if (on_contract == ndim_contract){
2279 break;
2280 }
2281 }
2282
2283 size_t on_new_core = 0;
2284 int done = 0;
2285 if (on_contract == ndim_contract){
2286 on_contract--;
2287 done = 1;
2288 }
2289 for (; ii < ft->dim; ii++){
2290 if ((dims_contract[on_contract] == ii) && (done == 0)){
2291 double * mat = qmarray_integrate_weighted(ft->cores[ii]);
2292
2293 struct Qmarray * newcore = qmam(newft->cores[on_new_core-1],mat,ft->cores[ii]->ncols);
2294 qmarray_free(newft->cores[on_new_core-1]);
2295 newft->cores[on_new_core-1] = qmarray_copy(newcore);
2296
2297 on_contract++;
2298 if (on_contract == ndim_contract){
2299 done = 1;
2300 on_contract--;
2301 }
2302 free(mat); mat = NULL;
2303 qmarray_free(newcore); newcore = NULL;
2304 }
2305 else{
2306 if (left_mult != NULL){
2307 newft->cores[on_new_core] = mqma(left_mult,ft->cores[ii],1);
2308 free(left_mult); left_mult = NULL;
2309 }
2310 else{
2311 newft->cores[on_new_core] = qmarray_copy(ft->cores[ii]);
2312 }
2313 on_new_core++;
2314 }
2315 }
2316
2317
2318 for (ii = 0; ii < newdim; ii++){
2319 newft->ranks[ii] = newft->cores[ii]->nrows;
2320 newft->ranks[ii+1] = newft->cores[ii]->ncols;
2321 }
2322
2323 return newft;
2324 }
2325
2326
2327 /********************************************************//**
2328 Inner product between two functions in FT form
2329
2330 \param[in] a - Function train 1
2331 \param[in] b - Function train 2
2332
2333 \return Inner product \f$ \int a(x)b(x) dx \f$
2334
2335 ***********************************************************/
function_train_inner(const struct FunctionTrain * a,const struct FunctionTrain * b)2336 double function_train_inner(const struct FunctionTrain * a,
2337 const struct FunctionTrain * b)
2338 {
2339 double out = 0.123456789;
2340
2341 double * temp = qmarray_kron_integrate(b->cores[0],a->cores[0]);
2342 double * temp2 = NULL;
2343
2344 for (size_t ii = 1; ii < a->dim; ii++){
2345 temp2 = qmarray_vec_kron_integrate(temp, a->cores[ii],b->cores[ii]);
2346 size_t stemp = a->cores[ii]->ncols * b->cores[ii]->ncols;
2347 free(temp);temp=NULL;
2348 temp = calloc_double(stemp);
2349 memmove(temp, temp2,stemp*sizeof(double));
2350
2351 free(temp2); temp2 = NULL;
2352 }
2353 out = temp[0];
2354 free(temp); temp=NULL;
2355 return out;
2356 }
2357
2358 /********************************************************//**
2359 Inner product between two functions in FT form with weighted
2360 space w(x)
2361
2362 \param[in] a - Function train 1
2363 \param[in] b - Function train 2
2364
2365 \return Inner product \f$ \int a(x)b(x) w(x) dx \f$
2366
2367 \note
2368 In order for this function to make sense a(x) and b(x)
2369 shoud have the same underlying basis
2370 ***********************************************************/
function_train_inner_weighted(const struct FunctionTrain * a,const struct FunctionTrain * b)2371 double function_train_inner_weighted(const struct FunctionTrain * a,
2372 const struct FunctionTrain * b)
2373 {
2374 double out = 0.123456789;
2375 size_t ii;
2376 double * temp = qmarray_kron_integrate_weighted(b->cores[0],a->cores[0]);
2377 double * temp2 = NULL;
2378
2379 //size_t ii;
2380 for (ii = 1; ii < a->dim; ii++){
2381 temp2 = qmarray_vec_kron_integrate_weighted(temp, a->cores[ii],b->cores[ii]);
2382 size_t stemp = a->cores[ii]->ncols * b->cores[ii]->ncols;
2383 free(temp);temp=NULL;
2384 temp = calloc_double(stemp);
2385 memmove(temp, temp2,stemp*sizeof(double));
2386
2387 free(temp2); temp2 = NULL;
2388 }
2389
2390 out = temp[0];
2391 free(temp); temp=NULL;
2392
2393 return out;
2394 }
2395
2396 /********************************************************//**
2397 Compute the L2 norm of a function in FT format
2398
2399 \param[in] a - Function train
2400
2401 \return L2 Norm \f$ \sqrt{int a^2(x) dx} \f$
2402 ***********************************************************/
function_train_norm2(const struct FunctionTrain * a)2403 double function_train_norm2(const struct FunctionTrain * a)
2404 {
2405 //printf("in norm2\n");
2406 double out = function_train_inner(a,a);
2407 if (out < -ZEROTHRESH){
2408 if (out * out > ZEROTHRESH){
2409 /* fprintf(stderr, "inner product of FT with itself should not be neg %G \n",out); */
2410 /* exit(1); */
2411 }
2412 }
2413 return sqrt(fabs(out));
2414 }
2415
2416
2417 /********************************************************//**
2418 Right orthogonalize the cores (except the first one) of the function train
2419
2420 \param[in,out] a - FT (overwritten)
2421 \param[in] aopts - approximation options to QR
2422
2423 \return ftrl - new ft with orthogonalized cores
2424 ***********************************************************/
2425 struct FunctionTrain *
function_train_orthor(struct FunctionTrain * a,struct MultiApproxOpts * aopts)2426 function_train_orthor(struct FunctionTrain * a,
2427 struct MultiApproxOpts * aopts)
2428 {
2429 struct OneApproxOpts * o = NULL;
2430 //printf("orthor\n");
2431 //right left sweep
2432 struct FunctionTrain * ftrl = function_train_alloc(a->dim);
2433 double * L = NULL;
2434 struct Qmarray * temp = NULL;
2435 size_t ii = 1;
2436 size_t core = a->dim-ii;
2437 L = calloc_double(a->cores[core]->nrows * a->cores[core]->nrows);
2438 memmove(ftrl->ranks,a->ranks,(a->dim+1)*sizeof(size_t));
2439
2440
2441 // update last core
2442 // printf("dim = %zu\n",a->dim);
2443 o = multi_approx_opts_get_aopts(aopts,core);
2444 /* printf("go orthor, core norm = %3.15G\n", qmarray_norm2(a->cores[core])); */
2445 ftrl->cores[core] = qmarray_householder_simple("LQ",a->cores[core],L,o);
2446 /* printf("go orthor, core norm = %3.15G\n", qmarray_norm2(ftrl->cores[core])); */
2447 /* ftrl->cores[core] = mqma(L, ftrl->cores[core], ftrl->cores[core]->nrows); */
2448 /* printf("done with LQ\n"); */
2449 for (ii = 2; ii < a->dim; ii++){
2450 core = a->dim-ii;
2451 //printf("on core %zu\n",core);
2452 temp = qmam(a->cores[core],L,ftrl->ranks[core+1]);
2453 free(L); L = NULL;
2454 L = calloc_double(ftrl->ranks[core]*ftrl->ranks[core]);
2455 o = multi_approx_opts_get_aopts(aopts,core);
2456 ftrl->cores[core] = qmarray_householder_simple("LQ",temp,L,o);
2457 /* printf("go orthor, core norm = %3.15G\n", qmarray_norm2(ftrl->cores[core])); */
2458 qmarray_free(temp); temp = NULL;
2459
2460 }
2461 ftrl->cores[0] = qmam(a->cores[0],L,ftrl->ranks[1]);
2462 /* ftrl->cores[0] = qmarray_copy(a->cores[0]); */
2463 /* printf("go orthor, core norm = %3.15G\n", qmarray_norm2(ftrl->cores[core])); */
2464 free(L); L = NULL;
2465 return ftrl;
2466 }
2467
2468 /********************************************************//**
2469 Rounding of a function train
2470
2471 \param[in] ain - FT
2472 \param[in] epsilon - threshold
2473 \param[in] aopts - approximation options to QR
2474
2475 \return ft - rounded function train
2476 ***********************************************************/
2477 struct FunctionTrain *
function_train_round(struct FunctionTrain * ain,double epsilon,struct MultiApproxOpts * aopts)2478 function_train_round(struct FunctionTrain * ain, double epsilon,
2479 struct MultiApproxOpts * aopts)
2480 {
2481 double normain = function_train_norm2(ain);
2482 /* printf("normain = %3.15G\n", normain); */
2483 if (normain < ZEROTHRESH){
2484 return function_train_constant(0.0, aopts);
2485 }
2486 struct OneApproxOpts * o = NULL;
2487 struct FunctionTrain * a = function_train_copy(ain);
2488 // size_t ii;
2489 size_t core;
2490 struct Qmarray * temp = NULL;
2491
2492 double delta = function_train_norm2(a);
2493 delta = delta * epsilon / sqrt(a->dim-1);
2494 //double delta = epsilon;
2495
2496 /* printf("Rounding Starting norm = %G\n",function_train_norm2(a)); */
2497 // printf("begin orho\n");
2498 struct FunctionTrain * ftrl = function_train_orthor(a,aopts);
2499
2500 /* printf("Rounding Ortho norm %G\n",function_train_norm2(ftrl)); */
2501 // printf("ortho gonalized\n");
2502 struct FunctionTrain * ft = function_train_alloc(a->dim);
2503 //struct FunctionTrain * ft = function_train_copy(ftrl);
2504 double * vt = NULL;
2505 double * s = NULL;
2506 double * sdiag = NULL;
2507 double * svt = NULL;
2508 size_t rank;
2509 size_t sizer;
2510 core = 0;
2511 ft->ranks[core] = 1;
2512 sizer = ftrl->cores[core]->ncols;
2513
2514 //*
2515 //printf("rank before = %zu\n",ftrl->ranks[core+1]);
2516 o = multi_approx_opts_get_aopts(aopts,core);
2517 rank = qmarray_truncated_svd(ftrl->cores[core], &(ft->cores[core]),
2518 &s, &vt, delta,o);
2519 //printf("rankdone\n");
2520
2521 //printf("rank after = %zu\n",rank);
2522 ft->ranks[core+1] = rank;
2523 sdiag = diag(rank,s);
2524 svt = calloc_double(rank * sizer);
2525 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rank,
2526 sizer, rank, 1.0, sdiag, rank, vt, rank, 0.0,
2527 svt, rank);
2528 //printf("rank\n");
2529 temp = mqma(svt,ftrl->cores[core+1],rank);
2530 //*/
2531
2532 free(vt); vt = NULL;
2533 free(s); s = NULL;
2534 free(sdiag); sdiag = NULL;
2535 free(svt); svt = NULL;
2536 for (core = 1; core < a->dim-1; core++){
2537 o = multi_approx_opts_get_aopts(aopts,core);
2538 rank = qmarray_truncated_svd(temp, &(ft->cores[core]), &s, &vt, delta,o);
2539 qmarray_free(temp); temp = NULL;
2540
2541 ft->ranks[core+1] = rank;
2542 sdiag = diag(rank,s);
2543 svt = calloc_double(rank * a->ranks[core+1]);
2544 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans, rank,
2545 a->ranks[core+1], rank, 1.0, sdiag, rank, vt, rank, 0.0,
2546 svt, rank);
2547 temp = mqma(svt,ftrl->cores[core+1],rank);
2548 free(vt); vt = NULL;
2549 free(s); s = NULL;
2550 free(sdiag); sdiag = NULL;
2551 free(svt); svt = NULL;
2552 }
2553 core = a->dim-1;
2554 ft->cores[core] = temp;
2555 ft->ranks[a->dim] = 1;
2556
2557 function_train_free(ftrl);
2558 /*
2559 for (ii = 0; ii < ft->dim; ii++){
2560 qmarray_roundt(&ft->cores[ii], epsilon);
2561 }
2562 */
2563
2564 function_train_free(a); a = NULL;
2565 return ft;
2566 }
2567
2568 /********************************************************//**
2569 Truncate functions in FT
2570
2571 \param[in,out] ft - FT
2572 \param[in] epsilon - threshold
2573 ***********************************************************/
function_train_truncate(struct FunctionTrain * ft,double epsilon)2574 void function_train_truncate(struct FunctionTrain * ft, double epsilon)
2575 {
2576 for (size_t ii = 0; ii < ft->dim; ii++){
2577 qmarray_roundt(&ft->cores[ii], epsilon);
2578 }
2579 }
2580
2581 /********************************************************//**
2582 Reapproximate a function train with new approximation options
2583 in a nonadaptive way. MultiApproxOpts should have all info needed
2584
2585 \param[in,out] ft - function train to reapproximate
2586 \param[in] aopts - approximation options
2587
2588 ***********************************************************/
2589 /* void function_train_reapprox_nonadapt(struct FunctionTrain * ft, struct MultiApproxOpts * aopts) */
2590 /* { */
2591 /* assert (ft != NULL); */
2592 /* assert (aopts != NULL); */
2593 /* struct OneApproxOpts * o = NULL; */
2594 /* for (size_t ii = 0; ii < ft->dim; ii++){ */
2595 /* o = multi_approx_opts_get_aopts(aopts,ii); */
2596 /* } */
2597 /* } */
2598
2599 /********************************************************//**
2600 Addition of two functions in FT format
2601
2602 \param[in] a - FT 1
2603 \param[in] b - FT 2
2604
2605 \return ft - function representing a+b
2606 ***********************************************************/
2607 struct FunctionTrain *
function_train_sum(const struct FunctionTrain * a,const struct FunctionTrain * b)2608 function_train_sum(const struct FunctionTrain * a,
2609 const struct FunctionTrain * b)
2610 {
2611 struct FunctionTrain * ft = function_train_alloc(a->dim);
2612
2613 // first core
2614 ft->cores[0] = qmarray_stackh(a->cores[0], b->cores[0]);
2615 ft->ranks[0] = a->ranks[0];
2616 // middle cores
2617
2618 size_t ii;
2619 for (ii = 1; ii < ft->dim-1; ii++){
2620 ft->cores[ii] = qmarray_blockdiag(a->cores[ii], b->cores[ii]);
2621 ft->ranks[ii] = a->ranks[ii] + b->ranks[ii];
2622 }
2623 ft->ranks[ft->dim-1] = a->ranks[a->dim-1] + b->ranks[b->dim-1];
2624
2625 // last core
2626 ft->cores[ft->dim-1] = qmarray_stackv(a->cores[a->dim-1],
2627 b->cores[b->dim-1]);
2628 ft->ranks[ft->dim] = a->ranks[ft->dim];
2629
2630 return ft;
2631 }
2632
2633 /********************************************************//**
2634 Scale a function train representation
2635
2636 \param[in,out] x - Function train
2637 \param[in] a - scaling factor
2638 ***********************************************************/
function_train_scale(struct FunctionTrain * x,double a)2639 void function_train_scale(struct FunctionTrain * x, double a)
2640 {
2641 struct GenericFunction ** temp = generic_function_array_daxpby(
2642 x->cores[0]->nrows*x->cores[0]->ncols,a,1,x->cores[0]->funcs,0.0,1,NULL);
2643 generic_function_array_free(x->cores[0]->funcs,
2644 x->cores[0]->nrows * x->cores[0]->ncols);
2645 x->cores[0]->funcs = temp;
2646 }
2647
2648 /* /\********************************************************\//\** */
2649 /* \f[ af(x) + b \f] */
2650
2651 /* \param[in] a - scaling factor */
2652 /* \param[in] b - offset */
2653 /* \param[in] f - object to scale */
2654 /* \param[in] epsilon - rounding tolerance */
2655
2656 /* \return ft - function representing a+b */
2657 /* ***********************************************************\/ */
2658 /* struct FunctionTrain * */
2659 /* function_train_afpb(double a, double b, */
2660 /* const struct FunctionTrain * f, double epsilon) */
2661 /* { */
2662 /* struct BoundingBox * bds = function_train_bds(f); */
2663 /* struct FunctionTrain * off = NULL; */
2664 /* if (f->cores[0]->funcs[0]->fc != LINELM){ */
2665 /* enum poly_type ptype = LEGENDRE; */
2666 /* off = function_train_constant(POLYNOMIAL,&ptype, */
2667 /* f->dim,b,bds,NULL); */
2668 /* } */
2669 /* else{ */
2670 /* off = function_train_constant(LINELM,NULL, */
2671 /* f->dim,b,bds,NULL); */
2672 /* } */
2673 /* struct FunctionTrain * af = function_train_copy(f); */
2674 /* function_train_scale(af,a); */
2675
2676 /* struct FunctionTrain * temp = function_train_sum(af,off); */
2677
2678
2679 /* //printf("round it \n"); */
2680 /* //printf("temp ranks \n"); */
2681 /* //iprint_sz(temp->dim+1,temp->ranks); */
2682
2683 /* struct FunctionTrain * ft = function_train_round(temp,epsilon); */
2684 /* //printf("got it \n"); */
2685
2686 /* function_train_free(off); off = NULL; */
2687 /* function_train_free(af); af = NULL; */
2688 /* function_train_free(temp); temp = NULL; */
2689 /* bounding_box_free(bds); bds = NULL; */
2690
2691 /* return ft; */
2692 /* } */
2693
2694 /********************************************************//**
2695 Product of two functions in function train form
2696
2697 \param[in] a - Function train 1
2698 \param[in] b - Function train 2
2699
2700 \return Product \f$ f(x) = a(x)b(x) \f$
2701 ***********************************************************/
2702 struct FunctionTrain *
function_train_product(const struct FunctionTrain * a,const struct FunctionTrain * b)2703 function_train_product(const struct FunctionTrain * a,
2704 const struct FunctionTrain * b)
2705 {
2706 struct FunctionTrain * ft = function_train_alloc(a->dim);
2707
2708 size_t ii;
2709 for (ii = 0; ii < ft->dim; ii++){
2710 ft->ranks[ii] = a->ranks[ii] * b->ranks[ii];
2711 ft->cores[ii] = qmarray_kron(a->cores[ii], b->cores[ii]);
2712 }
2713
2714 ft->ranks[ft->dim] = a->ranks[ft->dim] * b->ranks[ft->dim];
2715 return ft;
2716 }
2717
2718 /********************************************************//**
2719 Compute the L2 norm of the difference between two functions
2720
2721 \param[in] a - function train
2722 \param[in] b - function train 2
2723
2724 \return L2 difference \f$ \sqrt{ \int (a(x)-b(x))^2 dx } \f$
2725 ***********************************************************/
function_train_norm2diff(const struct FunctionTrain * a,const struct FunctionTrain * b)2726 double function_train_norm2diff(const struct FunctionTrain * a,
2727 const struct FunctionTrain * b)
2728 {
2729
2730 struct FunctionTrain * c = function_train_copy(b);
2731 function_train_scale(c,-1.0);
2732 struct FunctionTrain * d = function_train_sum(a,c);
2733 //printf("in function_train_norm2diff\n");
2734 double val = function_train_norm2(d);
2735 function_train_free(c);
2736 function_train_free(d);
2737 return val;
2738 }
2739
2740 /********************************************************//**
2741 Compute the relative L2 norm of the difference between two functions
2742
2743 \param[in] a - function train
2744 \param[in] b - function train 2
2745
2746 \return Relative L2 difference
2747 \f$ \sqrt{ \int (a(x)-b(x))^2 dx } / \lVert b(x) \rVert \f$
2748 ***********************************************************/
function_train_relnorm2diff(const struct FunctionTrain * a,const struct FunctionTrain * b)2749 double function_train_relnorm2diff(const struct FunctionTrain * a,
2750 const struct FunctionTrain * b)
2751 {
2752
2753 struct FunctionTrain * c = function_train_copy(b);
2754 function_train_scale(c,-1.0);
2755
2756 double den = function_train_inner(c,c);
2757
2758 //printf("compute sum \n");
2759 struct FunctionTrain * d = function_train_sum(a,c);
2760 //printf("computed \n");
2761 double num = function_train_inner(d,d);
2762
2763 double val = num;
2764 if (fabs(den) > ZEROTHRESH){
2765 val /= den;
2766 }
2767 if (val < -ZEROTHRESH){
2768 // fprintf(stderr, "relative error between two FT should not be neg %G<-%G \n",val,-ZEROTHRESH);
2769 // exit(1);
2770 }
2771 val = sqrt(fabs(val));
2772
2773 function_train_free(c); c = NULL;
2774 function_train_free(d); d = NULL;
2775 return val;
2776 }
2777
2778
2779 //////////////////////////////////////////////////////////////////////////
2780 //////////////////////////////////////////////////////////////////////////
2781 //////////////////////////////////////////////////////////////////////////
2782 //////////////////////////////////////////////////////////////////////////
2783 //////////////////////////////////////////////////////////////////////////
2784 //////////////////////////////////////////////////////////////////////////
2785 ///////////////////// ///////////////////////////
2786 ///////////////////// Cross Approximation ///////////////////////////
2787 ///////////////////// ///////////////////////////
2788 //////////////////////////////////////////////////////////////////////////
2789 //////////////////////////////////////////////////////////////////////////
2790 //////////////////////////////////////////////////////////////////////////
2791 //////////////////////////////////////////////////////////////////////////
2792 //////////////////////////////////////////////////////////////////////////
2793 //////////////////////////////////////////////////////////////////////////
2794 //////////////////////////////////////////////////////////////////////////
2795 //////////////////////////////////////////////////////////////////////////
2796
2797 /** \struct FtCrossArgs
2798 * \brief Arguments for function-train cross approximation
2799 * \var FtCrossArgs::dim
2800 * dimension of function
2801 * \var FtCrossArgs::ranks
2802 * (dim+1,) array of ranks
2803 * \var FtCrossArgs::epsilon
2804 * cross-approximation convergence criteria
2805 * \var FtCrossArgs::maxiter
2806 * maximum number of iteration for cross approximation
2807 * \var FtCrossArgs::epsround
2808 * rounding tolerance for adaptive rank cross approximation
2809 * \var FtCrossArgs::kickrank
2810 * size of rank increase for adaptation
2811 * \var FtCrossArgs::maxranks
2812 * maximum rank to go to during adaptation (dim-1,1)
2813 * \var FtCrossArgs::verbose
2814 * verbosity level (0,1,2)
2815 * */
2816 struct FtCrossArgs
2817 {
2818 size_t dim;
2819 size_t * ranks;
2820 double epsilon;
2821 size_t maxiter;
2822
2823 // adaptation parameters
2824 int adapt;
2825 double epsround;
2826 size_t kickrank;
2827 size_t * maxranks; //maxiteradapt;
2828
2829 int verbose;
2830 };
2831
2832 /***********************************************************//**
2833 Allocate space for cross approximation arguments
2834 Set elements to a default value
2835 ***************************************************************/
ft_cross_args_alloc(size_t dim,size_t start_rank)2836 struct FtCrossArgs * ft_cross_args_alloc(size_t dim, size_t start_rank)
2837 {
2838
2839 struct FtCrossArgs * ftc = malloc(sizeof(struct FtCrossArgs));
2840 if (ftc == NULL){
2841 fprintf(stderr,"Failure allocating FtCrossArgs\n");
2842 exit(1);
2843 }
2844
2845 ftc->dim = dim;
2846 ftc->ranks = calloc_size_t(dim+1);
2847 ftc->ranks[0] = 1;
2848 for (size_t ii = 1; ii < dim; ii++){
2849 ftc->ranks[ii] = start_rank;
2850 }
2851 ftc->ranks[dim] = 1;
2852 ftc->epsilon = 1e-10;
2853 ftc->maxiter = 5;
2854
2855 ftc->adapt = 1;
2856 ftc->epsround = 1e-10;
2857 ftc->kickrank = 5;
2858 // ftc->maxiteradapt = 5;
2859 ftc->maxranks = calloc_size_t(dim-1);
2860 for (size_t ii = 0; ii < dim-1;ii++){
2861 // 4 adaptation steps
2862 ftc->maxranks[ii] = ftc->ranks[ii+1] + ftc->kickrank*4;
2863 }
2864
2865 ftc->verbose = 0;
2866
2867 return ftc;
2868 }
2869
2870 /***********************************************************//**
2871 Set the rank for a particular index;
2872 ***************************************************************/
ft_cross_args_set_rank(struct FtCrossArgs * fca,size_t ind,size_t rank)2873 void ft_cross_args_set_rank(struct FtCrossArgs * fca, size_t ind, size_t rank)
2874 {
2875 fca->ranks[ind] = rank;
2876 }
2877
2878
2879 /***********************************************************//**
2880 Set the rounding tolerance
2881 ***************************************************************/
ft_cross_args_set_round_tol(struct FtCrossArgs * fca,double epsround)2882 void ft_cross_args_set_round_tol(struct FtCrossArgs * fca, double epsround)
2883 {
2884 fca->epsround = epsround;
2885 }
2886
2887 /***********************************************************//**
2888 Set the kickrank
2889 ***************************************************************/
ft_cross_args_set_kickrank(struct FtCrossArgs * fca,size_t kickrank)2890 void ft_cross_args_set_kickrank(struct FtCrossArgs * fca, size_t kickrank)
2891 {
2892 fca->kickrank = kickrank;
2893 }
2894
2895 /***********************************************************//**
2896 Set the maxiter
2897 ***************************************************************/
ft_cross_args_set_maxiter(struct FtCrossArgs * fca,size_t maxiter)2898 void ft_cross_args_set_maxiter(struct FtCrossArgs * fca, size_t maxiter)
2899 {
2900 fca->maxiter = maxiter;
2901 }
2902
2903 /***********************************************************//**
2904 Turn off adaptation
2905 ***************************************************************/
ft_cross_args_set_no_adaptation(struct FtCrossArgs * fca)2906 void ft_cross_args_set_no_adaptation(struct FtCrossArgs * fca)
2907 {
2908 fca->adapt = 0;
2909 }
2910
2911 /***********************************************************//**
2912 Turn onn adaptation
2913 ***************************************************************/
ft_cross_args_set_adaptation(struct FtCrossArgs * fca)2914 void ft_cross_args_set_adaptation(struct FtCrossArgs * fca)
2915 {
2916 fca->adapt = 1;
2917 }
2918
2919 /***********************************************************//**
2920 Set maximum ranks for adaptation
2921 ***************************************************************/
ft_cross_args_set_maxrank_all(struct FtCrossArgs * fca,size_t maxrank)2922 void ft_cross_args_set_maxrank_all(struct FtCrossArgs * fca, size_t maxrank)
2923 {
2924 // fca->maxiteradapt = maxiteradapt;
2925 for (size_t ii = 0; ii < fca->dim-1; ii++){
2926 fca->maxranks[ii] = maxrank;
2927 }
2928 }
2929
2930 /***********************************************************//**
2931 Set maximum ranks for adaptation per dimension
2932 ***************************************************************/
2933 void
ft_cross_args_set_maxrank_ind(struct FtCrossArgs * fca,size_t maxrank,size_t ind)2934 ft_cross_args_set_maxrank_ind(struct FtCrossArgs * fca, size_t maxrank, size_t ind)
2935 {
2936 // fca->maxiteradapt = maxiteradapt;
2937 assert (ind < (fca->dim-1));
2938 fca->maxranks[ind] = maxrank;
2939 }
2940
2941 /***********************************************************//**
2942 Set the cross approximation tolerance
2943 ***************************************************************/
ft_cross_args_set_cross_tol(struct FtCrossArgs * fca,double tol)2944 void ft_cross_args_set_cross_tol(struct FtCrossArgs * fca, double tol)
2945 {
2946 fca->epsilon = tol;
2947 }
2948
2949 /***********************************************************//**
2950 Set the verbosity level
2951 ***************************************************************/
ft_cross_args_set_verbose(struct FtCrossArgs * fca,int verbose)2952 void ft_cross_args_set_verbose(struct FtCrossArgs * fca, int verbose)
2953 {
2954 fca->verbose = verbose;
2955 }
2956
2957 /***********************************************************//**
2958 Get the ranks
2959 ***************************************************************/
ft_cross_args_get_ranks(const struct FtCrossArgs * fca)2960 size_t * ft_cross_args_get_ranks(const struct FtCrossArgs * fca)
2961 {
2962 assert (fca != NULL);
2963 return fca->ranks;
2964 }
2965
2966
2967
2968 /***********************************************************//**
2969 Copy cross approximation arguments
2970 ***************************************************************/
ft_cross_args_copy(const struct FtCrossArgs * fca)2971 struct FtCrossArgs * ft_cross_args_copy(const struct FtCrossArgs * fca)
2972 {
2973 if (fca == NULL){
2974 return NULL;
2975 }
2976 else{
2977 struct FtCrossArgs * f2 = malloc(sizeof(struct FtCrossArgs));
2978 assert (f2 != NULL);
2979 f2->dim = fca->dim;
2980 f2->ranks = calloc_size_t(fca->dim+1);
2981 memmove(f2->ranks,fca->ranks,(fca->dim+1) * sizeof(size_t));
2982 f2->epsilon = fca->epsilon;
2983 f2->maxiter = fca->maxiter;
2984 f2->adapt = fca->adapt;
2985 f2->epsround = fca->epsround;
2986 f2->kickrank = fca->kickrank;
2987 f2->maxranks = calloc_size_t(fca->dim-1);
2988 memmove(f2->maxranks,fca->maxranks, (fca->dim-1)*sizeof(size_t));
2989 f2->verbose = fca->verbose;
2990
2991 return f2;
2992 }
2993 }
2994
2995 /***********************************************************//**
2996 free cross approximation arguments
2997 ***************************************************************/
ft_cross_args_free(struct FtCrossArgs * fca)2998 void ft_cross_args_free(struct FtCrossArgs * fca)
2999 {
3000 if (fca != NULL){
3001 free(fca->ranks); fca->ranks = NULL;
3002 free(fca->maxranks); fca->maxranks = NULL;
3003 free(fca); fca = NULL;
3004 }
3005 }
3006
3007 struct Qmarray *
prepCore(size_t ii,size_t nrows,size_t ncols,struct Fwrap * fw,struct OneApproxOpts * o,struct CrossIndex ** left_ind,struct CrossIndex ** right_ind)3008 prepCore(size_t ii, size_t nrows, size_t ncols,
3009 struct Fwrap * fw, struct OneApproxOpts * o,
3010 struct CrossIndex ** left_ind,struct CrossIndex ** right_ind)
3011 {
3012 assert (fw != NULL);
3013 assert (o != NULL);
3014 assert (left_ind != NULL);
3015 assert (right_ind != NULL);
3016
3017 if (VPREPCORE) printf("in prepCore \n");
3018 size_t ncuts = nrows * ncols;
3019 fwrap_initialize_fiber_approx(fw,ii,ncuts);
3020 if (VPREPCORE) printf("initialized fiber approx ncuts=%zu \n",ncuts);
3021
3022 double * left = NULL, *right = NULL;
3023 size_t nl, nr;
3024 for (size_t jj = 0; jj < ncols; jj++){
3025 nr = 0;
3026 right = cross_index_get_node_value(right_ind[ii],jj,&nr);
3027 for (size_t kk = 0; kk < nrows; kk++){
3028 nl = 0;
3029 left = cross_index_get_node_value(left_ind[ii],kk,&nl);
3030 if (VPREPCORE){
3031 printf("left = "); dprint(nl,left);
3032 printf("right = "); dprint(nr,right);
3033 }
3034
3035 fwrap_add_fiber(fw,jj*nrows+kk,nl,left,nr,right);
3036 }
3037 }
3038
3039 if (VPREPCORE) printf("compute from fibercuts \n");
3040
3041 struct Qmarray * temp = qmarray_alloc(nrows,ncols);
3042 for (size_t jj = 0; jj < ncols; jj++){
3043 for (size_t kk = 0; kk < nrows; kk++){
3044 fwrap_set_which_fiber(fw,jj*nrows+kk);
3045 temp->funcs[jj*nrows+kk] =
3046 generic_function_approximate1d(o->fc,o->aopts,fw);
3047 }
3048 }
3049
3050 if (VPREPCORE) printf("computed!\n");
3051
3052
3053 fwrap_clean_fiber_approx(fw);
3054 return temp;
3055 }
3056
3057
3058
3059 /***********************************************************//**
3060 Initialize an FT from left and right indices
3061 *UNTESTED* USE AT YOUR OWN RISK
3062 *UNTESTED* USE AT YOUR OWN RISK
3063 \param[in] fw - wrapped function
3064 \param[in] ranks - ranks to use
3065 \param[in] left_ind - left indices
3066 \param[in] right_ind - right indices
3067 \param[in] apargs - approximation arguments
3068 *UNTESTED* USE AT YOUR OWN RISK
3069 \return function train decomposition of \f$ f \f$
3070 *UNTESTED* USE AT YOUR OWN RISK
3071 ***************************************************************/
3072 struct FunctionTrain *
function_train_init_from_fibers(struct Fwrap * fw,size_t * ranks,struct CrossIndex ** left_ind,struct CrossIndex ** right_ind,struct MultiApproxOpts * apargs)3073 function_train_init_from_fibers(struct Fwrap * fw,
3074 size_t * ranks,
3075 struct CrossIndex ** left_ind,
3076 struct CrossIndex ** right_ind,
3077 struct MultiApproxOpts * apargs)
3078 {
3079 assert (fw != NULL);
3080 assert (left_ind != NULL);
3081 assert (right_ind != NULL);
3082 assert (apargs != NULL);
3083
3084 size_t d = multi_approx_opts_get_dim(apargs);
3085 struct FunctionTrain * ft = function_train_alloc(d);
3086 memmove(ft->ranks, ranks, (d+1)*sizeof(size_t));
3087
3088
3089 struct OneApproxOpts * o = NULL;
3090 for (size_t ii = 0; ii < d; ii++){
3091 o = multi_approx_opts_get_aopts(apargs,ii);
3092 ft->cores[ii] = prepCore(ii,ranks[ii],ranks[ii+1],fw,o,
3093 left_ind,right_ind);
3094
3095 }
3096 struct FunctionTrain * ft2 = function_train_round(ft,ZEROTHRESH,apargs);
3097 function_train_free(ft); ft = NULL;
3098
3099 return ft2;
3100 }
3101
3102
3103
3104
3105 /***********************************************************//**
3106 Cross approximation of a of a dim-dimensional function
3107 (with adaptation)
3108
3109 \param[in] fw - wrapped function
3110 \param[in] cargs - cross approximation arguments
3111 \param[in,out] left_ind - left indices
3112 \param[in,out] right_ind - right indices
3113 \param[in] apargs - approximation arguments
3114 \param[in] optargs - fiber optimizationa arguments
3115 \param[in] ftref - reference ft for first error approximation
3116
3117 \return function train decomposition of \f$ f \f$
3118
3119 \note
3120 both left and right indices are nested
3121 ***************************************************************/
3122 struct FunctionTrain *
ftapprox_cross(struct Fwrap * fw,struct FtCrossArgs * cargs,struct CrossIndex ** left_ind,struct CrossIndex ** right_ind,struct MultiApproxOpts * apargs,struct FiberOptArgs * optargs,struct FunctionTrain * ftref)3123 ftapprox_cross(struct Fwrap * fw,
3124 struct FtCrossArgs * cargs,
3125 struct CrossIndex ** left_ind,
3126 struct CrossIndex ** right_ind,
3127 struct MultiApproxOpts * apargs,
3128 struct FiberOptArgs * optargs,
3129 struct FunctionTrain * ftref)
3130 {
3131 assert (fw != NULL);
3132 assert (cargs != NULL);
3133 assert (left_ind != NULL);
3134 assert (right_ind != NULL);
3135 assert (apargs != NULL);
3136 assert (optargs != NULL);
3137
3138 size_t dim = multi_approx_opts_get_dim(apargs);
3139 struct OneApproxOpts * o = NULL;
3140 void * opt = NULL;
3141 int info;
3142 /* size_t nrows, ii, oncore; */
3143 size_t ii,oncore;
3144 struct Qmarray * temp = NULL;
3145 struct Qmarray * Q = NULL;
3146 struct Qmarray * Qt = NULL;
3147 double * R = NULL;
3148 size_t * pivind = NULL;
3149 double * pivx = NULL;
3150 double diff, diff2, den;
3151
3152 struct FunctionTrain * ft = function_train_alloc(dim);
3153 memmove(ft->ranks, cargs->ranks, (dim+1)*sizeof(size_t));
3154 size_t * ranks = function_train_get_ranks(ft);
3155
3156 struct FunctionTrain * fti = function_train_copy(ftref);
3157 struct FunctionTrain * fti2 = NULL;
3158
3159 int done = 0;
3160 size_t iter = 0;
3161
3162 while (done == 0){
3163 if (cargs->verbose > 0)
3164 printf("cross iter=%zu \n",iter);
3165
3166 // left right sweep;
3167 // nrows = 1;
3168 for (ii = 0; ii < dim-1; ii++){
3169 o = multi_approx_opts_get_aopts(apargs,ii);
3170 opt = fiber_opt_args_get_opts(optargs,ii);
3171
3172 // printf("sub_type ftcross= %d\n",
3173 // *(int *)ft_approx_args_getst(apargs,ii));
3174 if (cargs->verbose > 1){
3175 printf(" ............. on left-right sweep (%zu/%zu)\n",ii,dim-1);
3176 }
3177 //printf("ii=%zu\n",ii);
3178 pivind = calloc_size_t(ft->ranks[ii+1]);
3179 pivx = calloc_double(ft->ranks[ii+1]);
3180
3181 if (VFTCROSS){
3182 printf("=======================================\n\n\n\n");
3183 printf( "prepCore \n");
3184 printf( "left index set = \n");
3185 print_cross_index(left_ind[ii]);
3186 printf( "right index set = \n");
3187 print_cross_index(right_ind[ii]);
3188 }
3189
3190 /* printf("prepCore\n"); */
3191 temp = prepCore(ii,ranks[ii],ranks[ii+1],fw,o,
3192 left_ind,right_ind);
3193 /* printf("prepped\n"); */
3194
3195 if (VFTCROSS == 2){
3196 printf ("got it \n");
3197 //print_qmarray(temp,0,NULL);
3198 struct Qmarray * tempp = qmarray_copy(temp);
3199 printf("core is \n");
3200 //print_qmarray(tempp,0,NULL);
3201 R = calloc_double(temp->ncols * temp->ncols);
3202 Q = qmarray_householder_simple("QR",temp,R,o);
3203 printf("R=\n");
3204 dprint2d_col(temp->ncols, temp->ncols, R);
3205 // print_qmarray(Q,0,NULL);
3206
3207 struct Qmarray * mult = qmam(Q,R,temp->ncols);
3208 //print_qmarray(Q,3,NULL);
3209 double difftemp = qmarray_norm2diff(mult,tempp);
3210 printf("difftemp = %3.15G\n",difftemp);
3211 qmarray_free(tempp);
3212 qmarray_free(mult);
3213 }
3214 else{
3215 R = calloc_double(temp->ncols * temp->ncols);
3216 Q = qmarray_householder_simple("QR", temp,R,o);
3217 }
3218
3219
3220 info = qmarray_maxvol1d(Q,R,pivind,pivx,o,opt);
3221
3222 if (VFTCROSS){
3223 printf( " got info=%d\n",info);
3224 printf("indices and pivots\n");
3225 iprint_sz(ft->ranks[ii+1],pivind);
3226 dprint(ft->ranks[ii+1],pivx);
3227 }
3228
3229 if (info < 0){
3230 fprintf(stderr, "no invertible submatrix in maxvol in cross\n");
3231 }
3232 if (info > 0){
3233 fprintf(stderr, " error in qmarray_maxvol1d \n");
3234 exit(1);
3235 }
3236
3237 cross_index_free(left_ind[ii+1]);
3238 if (ii > 0){
3239 left_ind[ii+1] =
3240 cross_index_create_nested_ind(0,ft->ranks[ii+1],pivind,
3241 pivx,left_ind[ii]);
3242 }
3243 else{
3244 left_ind[ii+1] = cross_index_alloc(1);
3245 for (size_t zz = 0; zz < ft->ranks[1]; zz++){
3246 cross_index_add_index(left_ind[1],1,&(pivx[zz]),sizeof(double));
3247 }
3248 }
3249
3250 qmarray_free(ft->cores[ii]); ft->cores[ii]=NULL;
3251 ft->cores[ii] = qmam(Q,R, temp->ncols);
3252 // nrows = left_ind[ii+1]->n;
3253
3254 qmarray_free(temp); temp = NULL;
3255 qmarray_free(Q); Q = NULL;
3256 free(pivind); pivind =NULL;
3257 free(pivx); pivx = NULL;
3258 free(R); R=NULL;
3259
3260 }
3261 ii = dim-1;
3262 if (cargs->verbose > 1){
3263 printf(" ............. on left-right sweep (%zu/%zu)\n",ii,dim-1);
3264 }
3265 qmarray_free(ft->cores[ii]); ft->cores[ii] = NULL;
3266 /* ft->cores[ii] = prepCore(ii,cargs->ranks[ii],f,args,bd, */
3267 /* left_ind,right_ind,cargs,apargs,1); */
3268 o = multi_approx_opts_get_aopts(apargs,ii);
3269 ft->cores[ii] = prepCore(ii,ranks[ii],ranks[ii+1],fw,o,
3270 left_ind,right_ind);
3271
3272 if (VFTCROSS == 2){
3273 printf ("got it \n");
3274 //print_qmarray(ft->cores[ii],0,NULL);
3275 printf("integral = %G\n",function_train_integrate(ft));
3276 struct FunctionTrain * tprod = function_train_product(ft,ft);
3277 printf("prod integral = %G\n",function_train_integrate(tprod));
3278 printf("norm2 = %G\n",function_train_norm2(ft));
3279 print_qmarray(tprod->cores[0],0,NULL,stdout);
3280 //print_qmarray(tprod->cores[1],0,NULL);
3281 function_train_free(tprod);
3282 }
3283
3284 if (VFTCROSS){
3285 printf("\n\n\n Index sets after Left-Right cross\n");
3286 for (ii = 0; ii < dim; ii++){
3287 printf("ii = %zu\n",ii);
3288 printf( "left index set = \n");
3289 print_cross_index(left_ind[ii]);
3290 printf( "right index set = \n");
3291 print_cross_index(right_ind[ii]);
3292 }
3293 printf("\n\n\n");
3294 }
3295
3296
3297 //printf("compute difference \n");
3298 //printf("norm fti = %G\n",function_train_norm2(fti));
3299 //printf("norm ft = %G\n",function_train_norm2(ft));
3300
3301 //printf("diff = %G\n",diff);
3302
3303 /* diff = function_train_norm2diff(ft,fti); */
3304 /* if (den > ZEROTHRESH){ */
3305 /* diff /= den; */
3306 /* } */
3307
3308 /* printf("compute norm\n"); */
3309
3310 diff = function_train_norm2diff(ft,fti);
3311 den = function_train_norm2(ft);
3312
3313 if (den > 1e0){
3314 diff = diff / den;
3315 }
3316
3317 // uncomment below to have error checking after a L/R sweep
3318 /* diff = function_train_relnorm2diff(ft,fti); */
3319 /* den = function_train_norm2(ft); */
3320 if (cargs->verbose > 0){
3321 den = function_train_norm2(ft);
3322 printf("...... New FT norm L/R Sweep = %3.15E\n",den);
3323 printf("...... Error L/R Sweep = %E\n",diff);
3324 }
3325
3326
3327 /* if (diff < cargs->epsilon){ */
3328 /* done = 1; */
3329 /* break; */
3330 /* } */
3331
3332
3333 // Can't do this because rank-adaptation doesn't care about
3334 // left_ind so need to end on a right-left sweep
3335 /* if (diff < cargs->epsilon){ */
3336 /* done = 1; */
3337 /* break; */
3338 /* } */
3339
3340
3341 /* function_train_free(fti); fti=NULL; */
3342 /* fti = function_train_copy(ft); */
3343
3344 //Up to here
3345 ////
3346
3347
3348 //printf("copied \n");
3349 //printf("copy diff= %G\n", function_train_norm2diff(ft,fti));
3350
3351 ///////////////////////////////////////////////////////
3352 // right-left sweep
3353 for (oncore = 1; oncore < dim; oncore++){
3354
3355 ii = dim-oncore;
3356 o = multi_approx_opts_get_aopts(apargs,ii);
3357 opt = fiber_opt_args_get_opts(optargs,ii);
3358
3359
3360 if (cargs->verbose > 1){
3361 printf(" ............. on right_left sweep (%zu/%zu)\n",ii,dim-1);
3362 }
3363
3364 // nrows = ft->ranks[ii];
3365
3366 if (VFTCROSS){
3367 printf("do prep\n");
3368 printf( "left index set = \n");
3369 print_cross_index(left_ind[ii]);
3370 printf( "right index set = \n");
3371 print_cross_index(right_ind[ii]);
3372 }
3373 //printf("prep core\n");
3374 /* temp = prepCore(ii,nrows,f,args,bd,left_ind,right_ind,cargs,apargs,0); */
3375 temp = prepCore(ii,ranks[ii],ranks[ii+1],fw,o,left_ind,right_ind);
3376
3377 /* printf("right after prepping core\n"); */
3378 /* print_qmarray(temp,0,NULL); */
3379 //printf("prepped core\n");
3380
3381 R = calloc_double(temp->nrows * temp->nrows);
3382 Q = qmarray_householder_simple("LQ", temp,R,o);
3383
3384 /* printf("right after taking QR of prepped core\n"); */
3385 /* print_qmarray(Q,0,NULL); */
3386
3387 Qt = qmarray_transpose(Q);
3388 pivind = calloc_size_t(ft->ranks[ii]);
3389 pivx = calloc_double(ft->ranks[ii]);
3390 info = qmarray_maxvol1d(Qt,R,pivind,pivx,o,opt);
3391
3392 if (VFTCROSS){
3393 printf("got info=%d\n",info);
3394 printf("indices and pivots\n");
3395 iprint_sz(ft->ranks[ii],pivind);
3396 dprint(ft->ranks[ii],pivx);
3397 }
3398
3399 //printf("got maxvol\n");
3400 if (info < 0){
3401 fprintf(stderr, "noinvertible submatrix in maxvol in rl cross\n");
3402 }
3403
3404 if (info > 0){
3405 fprintf(stderr, " error in qmarray_maxvol1d \n");
3406 exit(1);
3407 }
3408 qmarray_free(Q); Q = NULL;
3409
3410 //printf("pivx \n");
3411 //dprint(ft->ranks[ii], pivx);
3412
3413 Q = qmam(Qt,R, temp->nrows);
3414
3415 qmarray_free(ft->cores[ii]); ft->cores[ii] = NULL;
3416 ft->cores[ii] = qmarray_transpose(Q);
3417
3418 cross_index_free(right_ind[ii-1]);
3419 if (ii < dim-1){
3420 //printf("are we really here? oncore=%zu,ii=%zu\n",oncore,ii);
3421 right_ind[ii-1] =
3422 cross_index_create_nested_ind(1,ft->ranks[ii],pivind,
3423 pivx,right_ind[ii]);
3424 }
3425 else{
3426 //printf("lets update the cross index ii=%zu\n",ii);
3427 right_ind[ii-1] = cross_index_alloc(1);
3428 for (size_t zz = 0; zz < ft->ranks[ii]; zz++){
3429 cross_index_add_index(right_ind[ii-1],1,&(pivx[zz]),sizeof(double));
3430 }
3431 //printf("updated\n");
3432 }
3433
3434 qmarray_free(temp); temp = NULL;
3435 qmarray_free(Q); Q = NULL;
3436 qmarray_free(Qt); Qt = NULL;
3437 free(pivind);
3438 free(pivx);
3439 free(R); R=NULL;
3440
3441 /* break; //AG REMOVE THIS 5/18 */
3442
3443 }
3444
3445 ii = 0;
3446 qmarray_free(ft->cores[ii]); ft->cores[ii] = NULL;
3447
3448 if (cargs->verbose > 1)
3449 printf(" ............. on right_left sweep (%zu/%zu)\n",ii,dim-1);
3450 /* ft->cores[ii] = prepCore(ii,1,f,args,bd,left_ind,right_ind,cargs,apargs,-1); */
3451 o = multi_approx_opts_get_aopts(apargs,ii);
3452 ft->cores[ii] = prepCore(ii,ranks[ii],ranks[ii+1],fw,o,left_ind,right_ind);
3453 if (cargs->verbose > 1)
3454 printf(" ............. done with right left sweep\n");
3455
3456
3457 if (VFTCROSS){
3458 printf("\n\n\n Index sets after Right-left cross\n");
3459 for (ii = 0; ii < dim; ii++){
3460 printf("ii = %zu\n",ii);
3461 printf( "left index set = \n");
3462 print_cross_index(left_ind[ii]);
3463 printf( "right index set = \n");
3464 print_cross_index(right_ind[ii]);
3465 }
3466 printf("\n\n\n");
3467 }
3468
3469 diff = function_train_relnorm2diff(ft,fti);
3470 if (fti2 != NULL){
3471 diff2 = function_train_relnorm2diff(ft,fti2);
3472 }
3473 else{
3474 diff2 = diff;
3475 }
3476
3477 //den = function_train_norm2(ft);
3478 //diff = function_train_norm2diff(ft,fti);
3479 //if (den > ZEROTHRESH){
3480 // diff /= den;
3481 // }
3482
3483 if (cargs->verbose > 0){
3484 den = function_train_norm2(ft);
3485 printf("...... New FT norm R/L Sweep = %3.15E\n",den);
3486 printf("...... Error R/L Sweep = %E,%E\n",diff,diff2);
3487 }
3488
3489 if ( (diff2 < cargs->epsilon) && (diff < cargs->epsilon)){
3490 /* if ( (diff2 < cargs->epsilon) || (diff < cargs->epsilon)){ */
3491 /* if ( diff < cargs->epsilon){ */
3492 done = 1;
3493 break;
3494 }
3495
3496 function_train_free(fti2); fti2 = NULL;
3497 fti2 = function_train_copy(fti);
3498 function_train_free(fti); fti=NULL;
3499 fti = function_train_copy(ft);
3500
3501 iter++;
3502 if (iter == cargs->maxiter){
3503 done = 1;
3504 break;
3505 }
3506 }
3507
3508 function_train_free(fti); fti=NULL;
3509 function_train_free(fti2); fti2=NULL;
3510 return ft;
3511 }
3512
3513
3514 /***********************************************************//**
3515 Cross approximation of a of a dim-dimensional function with rank adaptation
3516
3517 \param[in] fw - wrapped function
3518 \param[in] cargs - cross approximation arguments
3519 \param[in,out] isl - left indices
3520 \param[in,out] isr - right indices
3521 \param[in] apargs - approximation arguments
3522 \param[in] optargs - fiber optimizationa arguments
3523 \param[in] ftref - reference ft for first error approximation
3524
3525 \return function train decomposition of f
3526
3527 \note
3528 both left and right indices are nested
3529 ***************************************************************/
3530 struct FunctionTrain *
ftapprox_cross_rankadapt(struct Fwrap * fw,struct FtCrossArgs * cargs,struct CrossIndex ** isl,struct CrossIndex ** isr,struct MultiApproxOpts * apargs,struct FiberOptArgs * optargs,struct FunctionTrain * ftref)3531 ftapprox_cross_rankadapt(struct Fwrap * fw,
3532 struct FtCrossArgs * cargs,
3533 struct CrossIndex ** isl,
3534 struct CrossIndex ** isr,
3535 struct MultiApproxOpts * apargs,
3536 struct FiberOptArgs * optargs,
3537 struct FunctionTrain * ftref)
3538 {
3539
3540
3541 assert (fw != NULL);
3542 assert (cargs != NULL);
3543 assert (isl != NULL);
3544 assert (isr != NULL);
3545 assert (apargs != NULL);
3546 assert (optargs != NULL);
3547
3548 size_t dim = multi_approx_opts_get_dim(apargs);
3549 double eps = cargs->epsround;
3550 size_t kickrank = cargs->kickrank;
3551
3552 struct FunctionTrain * ft = NULL;
3553
3554 ft = ftapprox_cross(fw,cargs,isl,isr,apargs,optargs,ftref);
3555 if (cargs->adapt == 0){
3556 return ft;
3557 }
3558 // printf("found left index\n");
3559 // print_cross_index(isl[1]);
3560 // printf("found right index\n");
3561 //print_cross_index(isr[0]);
3562
3563 //return ft;
3564 if (cargs->verbose > 0){
3565 printf("done with first cross... rounding\n");
3566 }
3567 size_t * ranks_found = calloc_size_t(dim+1);
3568 memmove(ranks_found,cargs->ranks,(dim+1)*sizeof(size_t));
3569
3570 struct FunctionTrain * ftc = function_train_copy(ft);
3571 struct FunctionTrain * ftr = function_train_round(ft,eps,apargs);
3572 if (cargs->verbose > 0){
3573 printf("rounded ranks = "); iprint_sz(dim+1,ftr->ranks);
3574 }
3575 //struct FunctionTrain * ftr = function_train_copy(ft);
3576 //return ftr;
3577 //printf("DOOONNTT FORGET MEEE HEERREEEE \n");
3578 int adapt = 0;
3579 size_t ii;
3580 for (ii = 1; ii < dim; ii++){
3581 /* printf("%zu == %zu\n",ranks_found[ii],ftr->ranks[ii]); */
3582 if (ranks_found[ii] == ftr->ranks[ii]){
3583 if (cargs->ranks[ii] < cargs->maxranks[ii-1]){
3584 adapt = 1;
3585 size_t kicksize;
3586 if ( (ranks_found[ii]+kickrank) <= cargs->maxranks[ii-1]){
3587 kicksize = kickrank;
3588 }
3589 else{
3590 kicksize = cargs->maxranks[ii-1] -ranks_found[ii];
3591 }
3592 cargs->ranks[ii] = ranks_found[ii] + kicksize;
3593
3594 // simply repeat the last nodes
3595 // this is not efficient but I don't have a better
3596 // idea. Could do it using random locations but
3597 // I don't want to.
3598 // I don't need to modify isl because it is rewritten on the
3599 // first left-right sweep anyways
3600 cross_index_copylast(isr[ii-1],kicksize);
3601
3602
3603 ranks_found[ii] = cargs->ranks[ii];
3604 }
3605 }
3606 }
3607
3608 //printf("adapt here! adapt=%zu\n",adapt);
3609
3610 size_t iter = 0;
3611 // double * xhelp = NULL;
3612 while ( adapt == 1 )
3613 {
3614 adapt = 0;
3615 if (cargs->verbose > 0){
3616 printf("adapting \n");
3617 printf("Increasing rank to \n");
3618 iprint_sz(ft->dim+1,cargs->ranks);
3619 }
3620
3621 function_train_free(ft); ft = NULL;
3622
3623 ft = ftapprox_cross(fw,cargs,isl,isr,apargs,optargs,ftref);
3624
3625 function_train_free(ftc); ftc = NULL;
3626 function_train_free(ftr); ftr = NULL;
3627
3628 //printf("copying\n");
3629 function_train_free(ftc); ftc = NULL;
3630 ftc = function_train_copy(ft);
3631 //printf("rounding\n");
3632 function_train_free(ftr); ftr = NULL;
3633 //ftr = function_train_copy(ft);//, eps);
3634 ftr = function_train_round(ft,eps,apargs);
3635 if (cargs->verbose > 0){
3636 printf("rounded ranks = "); iprint_sz(dim+1,ftr->ranks);
3637 }
3638 //printf("done rounding\n");
3639 for (ii = 1; ii < dim; ii++){
3640 if (ranks_found[ii] == ftr->ranks[ii]){
3641 if (cargs->ranks[ii] < cargs->maxranks[ii-1]){
3642 adapt = 1;
3643 size_t kicksize;
3644 if ( (ranks_found[ii]+kickrank) <= cargs->maxranks[ii-1]){
3645 kicksize = kickrank;
3646 }
3647 else{
3648 kicksize = cargs->maxranks[ii-1] -ranks_found[ii];
3649 }
3650 cargs->ranks[ii] = ranks_found[ii] + kicksize;
3651
3652
3653 // simply repeat the last nodes
3654 // this is not efficient but I don't have a better
3655 // idea. Could do it using random locations but
3656 // I don't want to.
3657 cross_index_copylast(isr[ii-1],kicksize);
3658
3659 ranks_found[ii] = cargs->ranks[ii];
3660 }
3661 }
3662 }
3663
3664 iter++;
3665 /* if (iter == fca->maxiteradapt) { */
3666 /* adapt = 0; */
3667 /* } */
3668 //adapt = 0;
3669
3670 }
3671
3672 if (cargs->verbose > 0){
3673 printf("Final norm = %G\n",function_train_norm2(ftr));
3674 printf("Final ranks: ");
3675 iprint_sz(ftr->dim+1,ftr->ranks);
3676 }
3677
3678 function_train_free(ft); ft = NULL;
3679 function_train_free(ftc); ftc = NULL;
3680 free(ranks_found);
3681 return ftr;
3682 }
3683
3684
3685 /***********************************************************//**
3686 Determine whether kristoffel weighting is active
3687
3688 \param[in] ft - function_train
3689
3690 \return 1 if active, 0 otherwise
3691
3692 \note It is assumed that if kristoffel is active on 1
3693 core it is active on all cores
3694 ***************************************************************/
function_train_is_kristoffel_active(const struct FunctionTrain * ft)3695 int function_train_is_kristoffel_active(const struct FunctionTrain * ft)
3696 {
3697 int active = qmarray_is_kristoffel_active(ft->cores[0]);
3698 for (size_t ii = 1; ii < ft->dim; ii++){
3699
3700 if (qmarray_is_kristoffel_active(ft->cores[ii]) != active){
3701 fprintf(stderr, "Kristoffel weighting cannot be active on some cores\n");
3702 fprintf(stderr, "and inactive on other cores\n");
3703 exit(1);
3704 }
3705 }
3706
3707 return active;
3708 }
3709
function_train_activate_kristoffel(struct FunctionTrain * ft)3710 void function_train_activate_kristoffel(struct FunctionTrain *ft)
3711 {
3712 for (size_t ii = 0; ii < ft->dim; ii++){
3713 qmarray_activate_kristoffel(ft->cores[ii]);
3714 }
3715 }
3716
function_train_deactivate_kristoffel(struct FunctionTrain * ft)3717 void function_train_deactivate_kristoffel(struct FunctionTrain *ft)
3718 {
3719 for (size_t ii = 0; ii < ft->dim; ii++){
3720 qmarray_deactivate_kristoffel(ft->cores[ii]);
3721 }
3722 }
3723
3724 /***********************************************************//**
3725 Get the kristoffel normalization factor
3726
3727 \param[in] ft - function train
3728 \param[in] x - location at which to obtain normalization
3729
3730 \return normalization factor
3731 ***************************************************************/
function_train_get_kristoffel_weight(const struct FunctionTrain * ft,const double * x)3732 double function_train_get_kristoffel_weight(const struct FunctionTrain * ft,
3733 const double * x)
3734 {
3735 double weight = 1.0;
3736 for (size_t ii = 0; ii < ft->dim; ii++){
3737 weight *= qmarray_get_kristoffel_weight(ft->cores[ii],x[ii]);
3738 }
3739
3740 return weight;
3741 }
3742
3743
3744
3745
3746
3747 //////////////////////////////////////////////////////////////////////////
3748 //////////////////////////////////////////////////////////////////////////
3749 //////////////////////////////////////////////////////////////////////////
3750 //////////////////////////////////////////////////////////////////////////
3751 //////////////////////////////////////////////////////////////////////////
3752 //////////////////////////////////////////////////////////////////////////
3753 ///////////////////// ///////////////////////////
3754 ///////////////////// Starting FT1D Array ///////////////////////////
3755 ///////////////////// ///////////////////////////
3756 //////////////////////////////////////////////////////////////////////////
3757 //////////////////////////////////////////////////////////////////////////
3758 //////////////////////////////////////////////////////////////////////////
3759 //////////////////////////////////////////////////////////////////////////
3760 //////////////////////////////////////////////////////////////////////////
3761 //////////////////////////////////////////////////////////////////////////
3762 //////////////////////////////////////////////////////////////////////////
3763 //////////////////////////////////////////////////////////////////////////
3764
3765 /***********************************************************//**
3766 Allocate a 1d array of function trains
3767
3768 \param[in] dimout - number of function trains
3769
3770 \return function train array
3771
3772 \note
3773 Each ft of the array is set to NULL;
3774 ***************************************************************/
ft1d_array_alloc(size_t dimout)3775 struct FT1DArray * ft1d_array_alloc(size_t dimout)
3776 {
3777 struct FT1DArray * fta = malloc(sizeof(struct FT1DArray));
3778 if (fta == NULL){
3779 fprintf(stderr, "Error allocating 1d function train array");
3780 exit(1);
3781 }
3782
3783 fta->size = dimout;
3784
3785 fta->ft = malloc(dimout * sizeof(struct FunctionTrain *));
3786 if (fta == NULL){
3787 fprintf(stderr, "Error allocating 1d function train array");
3788 exit(1);
3789 }
3790
3791 size_t ii;
3792 for (ii = 0; ii < dimout; ii ++){
3793 fta->ft[ii] = NULL;
3794 }
3795
3796 return fta;
3797 }
3798
3799 /***********************************************************//**
3800 Serialize a function train array
3801
3802 \param[in,out] ser - stream to serialize to
3803 \param[in] ft - function train array
3804 \param[in,out] totSizeIn - if NULL then serialize, if not NULL then return size
3805
3806 \return ptr - ser shifted by number of bytes
3807 ***************************************************************/
3808 unsigned char *
ft1d_array_serialize(unsigned char * ser,struct FT1DArray * ft,size_t * totSizeIn)3809 ft1d_array_serialize(unsigned char * ser, struct FT1DArray * ft,
3810 size_t *totSizeIn)
3811 {
3812
3813 // size -> ft1 -> ft2 -> ... ->ft[size]
3814
3815 // size
3816 size_t totSize = sizeof(size_t);
3817 size_t size_temp;
3818 for (size_t ii = 0; ii < ft->size; ii++){
3819 function_train_serialize(NULL,ft->ft[ii],&size_temp);
3820 totSize += size_temp;
3821 }
3822 if (totSizeIn != NULL){
3823 *totSizeIn = totSize;
3824 return ser;
3825 }
3826
3827 // serialize the size
3828 unsigned char * ptr = ser;
3829 ptr = serialize_size_t(ptr, ft->size);
3830
3831 // serialize each function train
3832 for (size_t ii = 0; ii < ft->size; ii++){
3833 ptr = function_train_serialize(ptr,ft->ft[ii],NULL);
3834 }
3835
3836 return ptr;
3837 }
3838
3839 /***********************************************************//**
3840 Deserialize a function train array
3841
3842 \param[in,out] ser - serialized function train array
3843 \param[in,out] ft - function_train array
3844
3845 \return ptr - shifted ser after deserialization
3846 ***************************************************************/
3847 unsigned char *
ft1d_array_deserialize(unsigned char * ser,struct FT1DArray ** ft)3848 ft1d_array_deserialize(unsigned char * ser, struct FT1DArray ** ft)
3849 {
3850 unsigned char * ptr = ser;
3851
3852 // deserialize the number of fts in the array
3853 size_t size;
3854 ptr = deserialize_size_t(ptr, &size);
3855 *ft = ft1d_array_alloc(size);
3856
3857 // deserialize each function train
3858 for (size_t ii = 0; ii < size; ii++){
3859 ptr = function_train_deserialize(ptr, &((*ft)->ft[ii]));
3860 }
3861 return ptr;
3862 }
3863
3864 /***********************************************************//**
3865 Save a function train array to file
3866
3867 \param[in] ft - function train array to save
3868 \param[in] filename - name of file to save to
3869
3870 \return success (1) or failure (0) of opening the file
3871 ***************************************************************/
ft1d_array_save(struct FT1DArray * ft,char * filename)3872 int ft1d_array_save(struct FT1DArray * ft, char * filename)
3873 {
3874
3875 FILE *fp;
3876 fp = fopen(filename, "w");
3877 if (fp == NULL){
3878 fprintf(stderr, "cat: can't open %s\n", filename);
3879 return 0;
3880 }
3881
3882 size_t totsize;
3883 ft1d_array_serialize(NULL,ft, &totsize);
3884
3885 unsigned char * data = malloc(totsize+sizeof(size_t));
3886 if (data == NULL){
3887 fprintf(stderr, "can't allocate space for saving density\n");
3888 return 0;
3889 }
3890
3891 // serialize size first!
3892 unsigned char * ptr = serialize_size_t(data,totsize);
3893 ptr = ft1d_array_serialize(ptr,ft,NULL);
3894
3895 fwrite(data,sizeof(unsigned char),totsize+sizeof(size_t),fp);
3896
3897 free(data); data = NULL;
3898 fclose(fp);
3899 return 1;
3900 }
3901
3902 /***********************************************************//**
3903 Load a function train array from a file
3904
3905 \param[in] filename - filename of file to load
3906
3907 \return ft if successfull NULL otherwise
3908 ***************************************************************/
ft1d_array_load(char * filename)3909 struct FT1DArray * ft1d_array_load(char * filename)
3910 {
3911 FILE *fp;
3912 fp = fopen(filename, "r");
3913 if (fp == NULL){
3914 fprintf(stderr, "cat: can't open %s\n", filename);
3915 return NULL;
3916 }
3917
3918 size_t totsize;
3919 size_t k = fread(&totsize,sizeof(size_t),1,fp);
3920 if ( k != 1){
3921 printf("error reading file %s\n",filename);
3922 return NULL;
3923 }
3924
3925 unsigned char * data = malloc(totsize);
3926 if (data == NULL){
3927 fprintf(stderr, "can't allocate space for loading density\n");
3928 return NULL;
3929 }
3930
3931 k = fread(data,sizeof(unsigned char),totsize,fp);
3932
3933 struct FT1DArray * ft = NULL;
3934 ft1d_array_deserialize(data,&ft);
3935
3936 free(data); data = NULL;
3937 fclose(fp);
3938 return ft;
3939 }
3940
3941 /***********************************************************//**
3942 Copy an array of function trains
3943
3944 \param[in] fta - array to coppy
3945
3946 \return ftb - copied array
3947 ***************************************************************/
ft1d_array_copy(const struct FT1DArray * fta)3948 struct FT1DArray * ft1d_array_copy(const struct FT1DArray * fta)
3949 {
3950 struct FT1DArray * ftb = ft1d_array_alloc(fta->size);
3951 size_t ii;
3952 for (ii = 0; ii < fta->size; ii++){
3953 ftb->ft[ii] = function_train_copy(fta->ft[ii]);
3954 }
3955 return ftb;
3956 }
3957
3958 /***********************************************************//**
3959 Free a 1d array of function trains
3960
3961 \param[in,out] fta - function train array to free
3962 ***************************************************************/
ft1d_array_free(struct FT1DArray * fta)3963 void ft1d_array_free(struct FT1DArray * fta)
3964 {
3965 if (fta != NULL){
3966 size_t ii = 0;
3967 for (ii = 0; ii < fta->size; ii++){
3968 function_train_free(fta->ft[ii]);
3969 fta->ft[ii] = NULL;
3970 }
3971 free(fta->ft);
3972 fta->ft = NULL;
3973 free(fta);
3974 fta = NULL;
3975 }
3976 }
3977
3978 /********************************************************//**
3979 Compute the gradient of a function train
3980
3981 \param[in] ft - Function train
3982
3983 \return gradient
3984 ***********************************************************/
function_train_gradient(const struct FunctionTrain * ft)3985 struct FT1DArray * function_train_gradient(const struct FunctionTrain * ft)
3986 {
3987 struct FT1DArray * ftg = ft1d_array_alloc(ft->dim);
3988 size_t ii;
3989 for (ii = 0; ii < ft->dim; ii++){
3990 ftg->ft[ii] = function_train_copy(ft);
3991 qmarray_free(ftg->ft[ii]->cores[ii]);
3992 ftg->ft[ii]->cores[ii] = NULL;
3993 ftg->ft[ii]->cores[ii] = qmarray_deriv(ft->cores[ii]);
3994 }
3995
3996 return ftg;
3997 }
3998
3999
4000
4001 /********************************************************//**
4002 Evaluate the gradient of a function train
4003
4004 \param[in] ft - Function train
4005 \param[in] x - location at which to evaluate the gradient
4006 \param[in,out] out - allocate space for the gradient
4007 ***********************************************************/
function_train_gradient_eval(const struct FunctionTrain * ft,const double * x,double * out)4008 void function_train_gradient_eval(const struct FunctionTrain * ft,
4009 const double * x,
4010 double * out)
4011 {
4012 /* (void)ft; */
4013 /* (void)x; */
4014 /* (void)out; */
4015 /* assert (1 == 0); */
4016
4017 size_t dim = ft->dim;
4018 size_t maxrank = function_train_get_maxrank(ft);
4019 size_t * ranks = ft->ranks;
4020
4021 size_t mr2 = maxrank*maxrank;
4022 double * evals_lr = calloc_double(maxrank * dim + maxrank);
4023 double * evals_rl = calloc_double(maxrank * dim + maxrank);
4024 double * grads_lr = calloc_double(maxrank * dim + maxrank);
4025 double * evals = calloc_double(mr2 * dim);
4026 double * grads = calloc_double(mr2 * dim);
4027 size_t onuni = 0;
4028 size_t indmem = 0;
4029
4030 size_t kk = 0;
4031 for (size_t col = 0; col < ranks[kk+1]; col++){
4032 evals_lr[col] = generic_function_1d_eval(ft->cores[kk]->funcs[col],x[kk]);
4033 grads_lr[col] = generic_function_deriv_eval(ft->cores[kk]->funcs[col],x[kk]);
4034 onuni++;
4035 }
4036
4037 for (kk = 1; kk < dim; kk++){
4038 for (size_t col = 0; col < ranks[kk+1]; col++){
4039 evals_lr[indmem+ranks[kk]+col] = 0.0;
4040 grads_lr[indmem+ranks[kk]+col] = 0.0;
4041 for (size_t row = 0; row < ranks[kk]; row++){
4042 evals[onuni] = generic_function_1d_eval(ft->cores[kk]->funcs[row + col * ranks[kk]],x[kk]);
4043 grads[onuni] = generic_function_deriv_eval(ft->cores[kk]->funcs[row + col * ranks[kk]],x[kk]);
4044 evals_lr[indmem+ranks[kk]+col] += evals_lr[indmem+row] * evals[onuni];
4045 grads_lr[indmem+ranks[kk]+col] += evals_lr[indmem+row] * grads[onuni];
4046 onuni++;
4047 }
4048 }
4049 indmem += ranks[kk];
4050 }
4051
4052 out[dim-1] = grads_lr[indmem];
4053
4054 size_t ind_eval = onuni-ranks[dim-1];
4055 size_t ind_rl = 0,r1,r2,backind;
4056 memmove(evals_rl, evals + ind_eval, ranks[dim-1] * sizeof(double));
4057
4058 for (size_t zz = 1; zz < dim-1; zz++)
4059 {
4060 backind = dim-1-zz;
4061 r1 = ranks[backind];
4062 r2 = ranks[backind+1];
4063
4064 indmem -= r2;
4065
4066 out[backind] = cblas_ddot(r2,grads_lr + indmem,1,evals_rl + ind_rl, 1);
4067
4068
4069 ind_eval -= ranks[backind]*ranks[backind+1];
4070
4071 cblas_dgemv(CblasColMajor, CblasNoTrans,
4072 r1,r2,1.0,
4073 evals + ind_eval, r1,
4074 evals_rl + ind_rl,1,
4075 0.0,evals_rl + ind_rl + r2, 1);
4076 ind_rl += r2;
4077 }
4078
4079 backind = 0;
4080 out[backind] = cblas_ddot(ranks[1],grads_lr,1,evals_rl + ind_rl, 1);
4081
4082 free(evals_lr); evals_lr = NULL;
4083 free(evals_rl); evals_rl = NULL;
4084 free(grads_lr); grads_lr = NULL;
4085 free(evals); evals = NULL;
4086 free(grads); grads = NULL;
4087 }
4088
4089
4090 /********************************************************//**
4091 Compute the Jacobian of a Function Train 1darray
4092
4093 \param[in] fta - Function train array
4094
4095 \return jacobian
4096 ***********************************************************/
ft1d_array_jacobian(const struct FT1DArray * fta)4097 struct FT1DArray * ft1d_array_jacobian(const struct FT1DArray * fta)
4098 {
4099 struct FT1DArray * jac = ft1d_array_alloc(fta->size * fta->ft[0]->dim);
4100 size_t ii,jj;
4101 for (ii = 0; ii < fta->ft[0]->dim; ii++){
4102 for (jj = 0; jj < fta->size; jj++){
4103 jac->ft[ii*fta->size+jj] = function_train_copy(fta->ft[jj]);
4104 qmarray_free(jac->ft[ii*fta->size+jj]->cores[ii]);
4105 jac->ft[ii*fta->size+jj]->cores[ii] = NULL;
4106 jac->ft[ii*fta->size+jj]->cores[ii] =
4107 qmarray_deriv(fta->ft[jj]->cores[ii]);
4108
4109 }
4110 }
4111 return jac;
4112 }
4113
4114
4115
4116 /********************************************************//**
4117 Compute the hessian of a function train
4118
4119 \param[in] fta - Function train
4120
4121 \return hessian of a function train
4122 ***********************************************************/
function_train_hessian(const struct FunctionTrain * fta)4123 struct FT1DArray * function_train_hessian(const struct FunctionTrain * fta)
4124 {
4125 struct FT1DArray * ftg = function_train_gradient(fta);
4126
4127 struct FT1DArray * fth = ft1d_array_jacobian(ftg);
4128
4129 ft1d_array_free(ftg); ftg = NULL;
4130 return fth;
4131 }
4132
4133 /********************************************************//**
4134 Compute the diagpnal of the hessian a function train
4135
4136 \param[in] ft - Function train
4137
4138 \return gradient
4139 ***********************************************************/
function_train_hessian_diag(const struct FunctionTrain * ft)4140 struct FT1DArray * function_train_hessian_diag(const struct FunctionTrain * ft)
4141 {
4142 struct FT1DArray * ftg = ft1d_array_alloc(ft->dim);
4143 size_t ii;
4144 for (ii = 0; ii < ft->dim; ii++){
4145 ftg->ft[ii] = function_train_copy(ft);
4146 qmarray_free(ftg->ft[ii]->cores[ii]);
4147 ftg->ft[ii]->cores[ii] = NULL;
4148 ftg->ft[ii]->cores[ii] = qmarray_dderiv(ft->cores[ii]);
4149 }
4150
4151 return ftg;
4152 }
4153
4154 /********************************************************//**
4155 Compute the diagpnal of the hessian a function train
4156 enforce periodic boundary conditions
4157
4158 \param[in] ft - Function train
4159
4160 \return gradient
4161 ***********************************************************/
function_train_hessian_diag_periodic(const struct FunctionTrain * ft)4162 struct FT1DArray * function_train_hessian_diag_periodic(const struct FunctionTrain * ft)
4163 {
4164 struct FT1DArray * ftg = ft1d_array_alloc(ft->dim);
4165 size_t ii;
4166 for (ii = 0; ii < ft->dim; ii++){
4167 ftg->ft[ii] = function_train_copy(ft);
4168 qmarray_free(ftg->ft[ii]->cores[ii]);
4169 ftg->ft[ii]->cores[ii] = NULL;
4170 ftg->ft[ii]->cores[ii] = qmarray_dderiv_periodic(ft->cores[ii]);
4171 }
4172
4173 return ftg;
4174 }
4175
4176 /********************************************************//**
4177 Scale a function train array
4178
4179 \param[in,out] fta - function train array
4180 \param[in] n - number of elements in the array to scale
4181 \param[in] inc - increment between elements of array
4182 \param[in] scale - value by which to scale
4183 ***********************************************************/
ft1d_array_scale(struct FT1DArray * fta,size_t n,size_t inc,double scale)4184 void ft1d_array_scale(struct FT1DArray * fta, size_t n, size_t inc, double scale)
4185 {
4186 size_t ii;
4187 for (ii = 0; ii < n; ii++){
4188 function_train_scale(fta->ft[ii*inc],scale);
4189 }
4190 }
4191
4192 /********************************************************//**
4193 Evaluate a function train 1darray
4194
4195 \param[in] fta - Function train array to evaluate
4196 \param[in] x - location at which to obtain evaluations
4197
4198 \return evaluation
4199 ***********************************************************/
ft1d_array_eval(const struct FT1DArray * fta,const double * x)4200 double * ft1d_array_eval(const struct FT1DArray * fta, const double * x)
4201 {
4202 double * out = calloc_double(fta->size);
4203 size_t ii;
4204 for (ii = 0; ii < fta->size; ii++){
4205 out[ii] = function_train_eval(fta->ft[ii], x);
4206 }
4207 return out;
4208 }
4209
4210 /********************************************************//**
4211 Evaluate a function train 1darray
4212
4213 \param[in] fta - Function train array to evaluate
4214 \param[in] x - location at which to obtain evaluations
4215 \param[in,out] out - evaluation
4216
4217 ***********************************************************/
ft1d_array_eval2(const struct FT1DArray * fta,const double * x,double * out)4218 void ft1d_array_eval2(const struct FT1DArray * fta, const double * x, double * out)
4219 {
4220 size_t ii;
4221 for (ii = 0; ii < fta->size; ii++){
4222 out[ii] = function_train_eval(fta->ft[ii], x);
4223 }
4224 }
4225
4226
4227 /********************************************************
4228 Multiply together and sum the elements of two function train arrays
4229 \f[
4230 out(x) = \sum_{i=1}^{N} coeff[i] f_i(x) g_i(x)
4231 \f]
4232
4233 \param[in] N - number of function trains in each array
4234 \param[in] coeff - coefficients to multiply each element
4235 \param[in] f - first array
4236 \param[in] g - second array
4237 \param[in] epsilon - rounding accuracy
4238
4239 \return function train
4240
4241 \note
4242 this needs the approximation options because it does rounding
4243 this is a bit unfortunate -- commenting it out for now
4244 ***********************************************************/
4245 /* struct FunctionTrain * */
4246 /* ft1d_array_sum_prod(size_t N, double * coeff, */
4247 /* const struct FT1DArray * f, const struct FT1DArray * g, */
4248 /* double epsilon) */
4249 /* { */
4250
4251 /* struct FunctionTrain * ft1 = NULL; */
4252 /* struct FunctionTrain * ft2 = NULL; */
4253 /* struct FunctionTrain * out = NULL; */
4254 /* struct FunctionTrain * temp = NULL; */
4255
4256 /* // printf("compute product\n"); */
4257 /* temp = function_train_product(f->ft[0],g->ft[0]); */
4258 /* function_train_scale(temp,coeff[0]); */
4259 /* // printf("scale N = %zu\n",N); */
4260 /* out = function_train_round(temp,epsilon); */
4261 /* // printf("rounded\n"); */
4262 /* function_train_free(temp); temp = NULL; */
4263 /* size_t ii; */
4264 /* for (ii = 1; ii < N; ii++){ */
4265 /* temp = function_train_product(f->ft[ii],g->ft[ii]); */
4266 /* function_train_scale(temp,coeff[ii]); */
4267
4268 /* ft2 = function_train_round(temp,epsilon); */
4269 /* ft1 = function_train_sum(out, ft2); */
4270
4271 /* function_train_free(temp); temp = NULL; */
4272 /* function_train_free(ft2); ft2 = NULL; */
4273 /* function_train_free(out); out = NULL; */
4274
4275 /* out = function_train_round(ft1,epsilon); */
4276
4277 /* function_train_free(ft1); ft1 = NULL; */
4278 /* } */
4279
4280 /* return out; */
4281 /* } */
4282
4283
4284 //////////////////////////////////////////////////////////////////////////
4285 //////////////////////////////////////////////////////////////////////////
4286 //////////////////////////////////////////////////////////////////////////
4287 //////////////////////////////////////////////////////////////////////////
4288 //////////////////////////////////////////////////////////////////////////
4289 //////////////////////////////////////////////////////////////////////////
4290 ///////////////////// ///////////////////////////
4291 ///////////////////// BLAS TYPE INTERFACE ///////////////////////////
4292 ///////////////////// ///////////////////////////
4293 //////////////////////////////////////////////////////////////////////////
4294 //////////////////////////////////////////////////////////////////////////
4295 //////////////////////////////////////////////////////////////////////////
4296 //////////////////////////////////////////////////////////////////////////
4297 //////////////////////////////////////////////////////////////////////////
4298 //////////////////////////////////////////////////////////////////////////
4299 //////////////////////////////////////////////////////////////////////////
4300 //////////////////////////////////////////////////////////////////////////
4301
4302
4303 /***********************************************************//**
4304 Computes
4305 \f[
4306 y \leftarrow \texttt{round}(a x + y, epsilon)
4307 \f]
4308
4309 \param[in] a - scaling factor
4310 \param[in] x - first function train
4311 \param[in,out] y - second function train
4312 \param[in] epsilon - rounding accuracy (0 for exact)
4313
4314 \note
4315 putting NULL it to round will cause an error
4316 ***************************************************************/
c3axpy(double a,struct FunctionTrain * x,struct FunctionTrain ** y,double epsilon,struct MultiApproxOpts * opts)4317 void c3axpy(double a, struct FunctionTrain * x, struct FunctionTrain ** y,
4318 double epsilon, struct MultiApproxOpts * opts)
4319 {
4320
4321 struct FunctionTrain * temp = function_train_copy(x);
4322 function_train_scale(temp,a);
4323 struct FunctionTrain * z = function_train_sum(temp,*y);
4324 function_train_free(*y); *y = NULL;
4325 if (epsilon > 0){
4326 *y = function_train_round(z,epsilon,opts);
4327 function_train_free(z); z = NULL;
4328 }
4329 else{
4330 *y = z;
4331 }
4332 function_train_free(temp); temp = NULL;
4333 }
4334
4335 /***********************************************************//**
4336 Computes
4337 \f$
4338 \langle x,y \rangle
4339 \f$
4340
4341 \param[in] x - first function train
4342 \param[in] y - second function train
4343
4344 \return out - inner product between two function trains
4345 ***************************************************************/
c3dot(struct FunctionTrain * x,struct FunctionTrain * y)4346 double c3dot(struct FunctionTrain * x, struct FunctionTrain * y)
4347 {
4348 double out = function_train_inner(x,y);
4349 return out;
4350 }
4351
4352
4353
4354 /***********************************************************//**
4355 Computes
4356 \f[
4357 y \leftarrow alpha \sum_{i=1}^n \texttt{round}(\texttt{product}(A[i*inca],x),epsilon) +
4358 beta y
4359 \f]
4360 \f[
4361 y \leftarrow \texttt{round}(y,epsilon)
4362 \f]
4363
4364 \note
4365 Also rounds after every summation. Will cause an error because putting a NULL into it
4366 ***************************************************************/
c3gemv(double alpha,size_t n,struct FT1DArray * A,size_t inca,struct FunctionTrain * x,double beta,struct FunctionTrain * y,double epsilon,struct MultiApproxOpts * opts)4367 void c3gemv(double alpha, size_t n, struct FT1DArray * A, size_t inca,
4368 struct FunctionTrain * x,double beta, struct FunctionTrain * y,
4369 double epsilon, struct MultiApproxOpts * opts)
4370 {
4371
4372 size_t ii;
4373 assert (y != NULL);
4374 function_train_scale(y,beta);
4375
4376
4377 if (epsilon > 0){
4378 struct FunctionTrain * runinit = function_train_product(A->ft[0],x);
4379 struct FunctionTrain * run = function_train_round(runinit,epsilon,opts);
4380 for (ii = 1; ii < n; ii++){
4381 struct FunctionTrain * temp =
4382 function_train_product(A->ft[ii*inca],x);
4383 struct FunctionTrain * tempround = function_train_round(temp,epsilon,opts);
4384 c3axpy(1.0,tempround,&run,epsilon,opts);
4385 function_train_free(temp); temp = NULL;
4386 function_train_free(tempround); tempround = NULL;
4387 }
4388 c3axpy(alpha,run,&y,epsilon,opts);
4389 function_train_free(run); run = NULL;
4390 function_train_free(runinit); runinit = NULL;
4391 }
4392 else{
4393 struct FunctionTrain * run = function_train_product(A->ft[0],x);
4394 for (ii = 1; ii < n; ii++){
4395 struct FunctionTrain * temp =
4396 function_train_product(A->ft[ii*inca],x);
4397 c3axpy(1.0,temp,&run,epsilon,opts);
4398 function_train_free(temp); temp = NULL;
4399 }
4400 c3axpy(alpha,run,&y,epsilon,opts);
4401 function_train_free(run); run = NULL;
4402 }
4403
4404 }
4405
4406 /***********************************************************//**
4407 Computes for \f$ i=1 \ldots n\f$
4408 \f[
4409 y[i*incy] \leftarrow \texttt{round}(a * x[i*incx] + y[i*incy],epsilon)
4410 \f]
4411
4412 ***************************************************************/
c3vaxpy(size_t n,double a,struct FT1DArray * x,size_t incx,struct FT1DArray ** y,size_t incy,double epsilon,struct MultiApproxOpts * opts)4413 void c3vaxpy(size_t n, double a, struct FT1DArray * x, size_t incx,
4414 struct FT1DArray ** y, size_t incy, double epsilon,
4415 struct MultiApproxOpts * opts)
4416 {
4417 size_t ii;
4418 for (ii = 0; ii < n; ii++){
4419 c3axpy(a,x->ft[ii*incx],&((*y)->ft[ii*incy]),epsilon, opts);
4420 }
4421 }
4422
4423 /***********************************************************//**
4424 Computes for \f$ i=1 \ldots n\f$
4425 \f[
4426 z \leftarrow alpha\sum_{i=1}^n y[incy*ii]*x[ii*incx] + beta * z
4427 \f]
4428
4429 ***************************************************************/
c3vaxpy_arr(size_t n,double alpha,struct FT1DArray * x,size_t incx,double * y,size_t incy,double beta,struct FunctionTrain ** z,double epsilon,struct MultiApproxOpts * opts)4430 void c3vaxpy_arr(size_t n, double alpha, struct FT1DArray * x,
4431 size_t incx, double * y, size_t incy, double beta,
4432 struct FunctionTrain ** z, double epsilon,
4433 struct MultiApproxOpts * opts)
4434 {
4435 assert (*z != NULL);
4436 function_train_scale(*z,beta);
4437
4438 size_t ii;
4439 for (ii = 0; ii < n; ii++){
4440 c3axpy(alpha*y[incy*ii],x->ft[ii*incx], z,epsilon, opts);
4441 }
4442 }
4443
4444 /***********************************************************//**
4445 Computes
4446 \f[
4447 z \leftarrow \texttt{round}(a\sum_{i=1}^n \texttt{round}(x[i*incx]*y[i*incy],epsilon) + beta * z,epsilon)
4448
4449 \f]
4450
4451 \note
4452 Also rounds after every summation.
4453
4454 ***************************************************************/
c3vprodsum(size_t n,double a,struct FT1DArray * x,size_t incx,struct FT1DArray * y,size_t incy,double beta,struct FunctionTrain ** z,double epsilon,struct MultiApproxOpts * opts)4455 void c3vprodsum(size_t n, double a, struct FT1DArray * x, size_t incx,
4456 struct FT1DArray * y, size_t incy, double beta,
4457 struct FunctionTrain ** z, double epsilon, struct MultiApproxOpts * opts)
4458 {
4459 assert (*z != NULL);
4460 function_train_scale(*z,beta);
4461
4462 size_t ii;
4463
4464 if (epsilon > 0){
4465 struct FunctionTrain * runinit =
4466 function_train_product(x->ft[0],y->ft[0]);
4467 struct FunctionTrain * run = function_train_round(runinit,epsilon,opts);
4468 for (ii = 1; ii < n; ii++){
4469 struct FunctionTrain * tempinit =
4470 function_train_product(x->ft[ii*incx],y->ft[ii*incy]);
4471 struct FunctionTrain * temp = function_train_round(tempinit,epsilon,opts);
4472 c3axpy(1.0,temp,&run,epsilon, opts);
4473 function_train_free(temp); temp = NULL;
4474 function_train_free(tempinit); tempinit = NULL;
4475 }
4476 c3axpy(a,run,z,epsilon, opts);
4477 function_train_free(run); run = NULL;
4478 function_train_free(runinit); runinit = NULL;
4479 }
4480 else{
4481 struct FunctionTrain * run =
4482 function_train_product(x->ft[0],y->ft[0]);
4483 for (ii = 1; ii < n; ii++){
4484 struct FunctionTrain * temp =
4485 function_train_product(x->ft[ii*incx],y->ft[ii*incy]);
4486 c3axpy(1.0,temp,&run,epsilon, opts);
4487 function_train_free(temp); temp = NULL;
4488 }
4489 c3axpy(a,run,z,epsilon, opts);
4490 function_train_free(run); run = NULL;
4491 }
4492 }
4493
4494 /***********************************************************//**
4495 Computes for \f$ i = 1 \ldots m \f$
4496 \f[
4497 y[i*incy] \leftarrow alpha*\sum_{j=1}^n \texttt{product}(A[i,j*lda],x[j*incx]) + beta * y[i*incy]
4498 \f]
4499
4500 \note
4501 Rounds with tolerance epsilon after summation and multiplication. Not shown to avoid clutter.
4502
4503 ***************************************************************/
c3vgemv(size_t m,size_t n,double alpha,struct FT1DArray * A,size_t lda,struct FT1DArray * x,size_t incx,double beta,struct FT1DArray ** y,size_t incy,double epsilon,struct MultiApproxOpts * opts)4504 void c3vgemv(size_t m, size_t n, double alpha, struct FT1DArray * A, size_t lda,
4505 struct FT1DArray * x, size_t incx, double beta, struct FT1DArray ** y,
4506 size_t incy, double epsilon, struct MultiApproxOpts * opts)
4507 {
4508
4509 size_t ii;
4510 assert (*y != NULL);
4511 ft1d_array_scale(*y,m,incy,beta);
4512
4513 struct FT1DArray * run = ft1d_array_alloc(m);
4514 for (ii = 0; ii < m; ii++){
4515 struct FT1DArray ftatemp;
4516 ftatemp.ft = A->ft+ii;
4517 c3vprodsum(n,1.0,&ftatemp,lda,x,incx,0.0,&(run->ft[ii*incy]),epsilon, opts);
4518 }
4519 c3vaxpy(m,alpha,run,1,y,incy,epsilon, opts);
4520
4521 ft1d_array_free(run); run = NULL;
4522 }
4523
4524 /***********************************************************//**
4525 Computes for \f$ i = 1 \ldots m \f$
4526 \f[
4527 y[i*incy] \leftarrow alpha*\sum_{j=1}^n \texttt{product}(A[i,j],B[j*incb]) + beta * y[i*incy]
4528 \f]
4529
4530 \note
4531 Rounds with tolerance epsilon after summation and multiplication. Not shown to avoid clutter.
4532 trans = 0 means not to transpose A
4533 trans = 1 means to transpose A
4534 ***************************************************************/
c3vgemv_arr(int trans,size_t m,size_t n,double alpha,double * A,size_t lda,struct FT1DArray * B,size_t incb,double beta,struct FT1DArray ** y,size_t incy,double epsilon,struct MultiApproxOpts * opts)4535 void c3vgemv_arr(int trans, size_t m, size_t n, double alpha, double * A,
4536 size_t lda, struct FT1DArray * B, size_t incb, double beta,
4537 struct FT1DArray ** y, size_t incy, double epsilon, struct MultiApproxOpts * opts)
4538 {
4539 size_t ii;
4540 assert (*y != NULL);
4541 ft1d_array_scale(*y,m,incy,beta);
4542
4543 if (trans == 0){
4544 for (ii = 0; ii < m; ii++){
4545 c3vaxpy_arr(n,alpha,B,incb,A+ii,lda,beta,&((*y)->ft[ii*incy]),epsilon, opts);
4546 }
4547 }
4548 else if (trans == 1){
4549 for (ii = 0; ii < m; ii++){
4550 c3vaxpy_arr(n,alpha,B,incb,A+ii*lda,1,beta,&((*y)->ft[ii*incy]),epsilon, opts);
4551 }
4552 }
4553 }
4554