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
40 /** \file qmarray.c
41 * Provides algorithms for qmarray manipulation
42 */
43
44 #include <stdlib.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 "quasimatrix.h"
55 #include "qmarray.h"
56
57 /* #define ZEROTHRESH 1e2*DBL_EPSILON */
58 /* #define ZERO 1e2*DBL_EPSILON */
59
60 /* #define ZEROTHRESH 1e0*DBL_EPSILON */
61 /* #define ZERO 1e0*DBL_EPSILON */
62
63 #define ZEROTHRESH 1e2*DBL_EPSILON
64 #define ZERO 1e2*DBL_EPSILON
65
66 /* #define ZEROTHRESH 1e-200 */
67 /* #define ZERO 1e-200 */
68
69 #ifndef VQMALU
70 #define VQMALU 0
71 #endif
72
73 #ifndef VQMAMAXVOL
74 #define VQMAMAXVOL 0
75 #endif
76
77 #ifndef VQMAHOUSEHOLDER
78 #define VQMAHOUSEHOLDER 0
79 #endif
80
81 /***********************************************************//**
82 Allocate space for a qmarray
83
84 \param[in] nrows - number of rows of quasimatrix
85 \param[in] ncols - number of cols of quasimatrix
86
87 \return qmarray
88 ***************************************************************/
qmarray_alloc(size_t nrows,size_t ncols)89 struct Qmarray * qmarray_alloc(size_t nrows, size_t ncols){
90
91 struct Qmarray * qm = NULL;
92 if ( NULL == (qm = malloc(sizeof(struct Qmarray)))){
93 fprintf(stderr, "failed to allocate memory for qmarray.\n");
94 exit(1);
95 }
96 qm->nrows = nrows;
97 qm->ncols = ncols;
98 if ( NULL == (qm->funcs=malloc(nrows*ncols*sizeof(struct GenericFunction *)))){
99 fprintf(stderr, "failed to allocate memory for qmarray.\n");
100 exit(1);
101 }
102 size_t ii;
103 for (ii = 0; ii < nrows*ncols; ii++){
104 qm->funcs[ii] = NULL;
105 }
106 return qm;
107 }
108
109 /***********************************************************//**
110 Free memory allocated to qmarray
111
112 \param[in] qm - qmarray
113 ***************************************************************/
qmarray_free(struct Qmarray * qm)114 void qmarray_free(struct Qmarray * qm){
115
116 if (qm != NULL){
117 size_t ii = 0;
118 for (ii = 0; ii < qm->nrows*qm->ncols; ii++){
119 generic_function_free(qm->funcs[ii]);
120 }
121 free(qm->funcs);
122 free(qm);qm=NULL;
123 }
124 }
125
126 /***********************************************************//**
127 (deep)copy qmarray
128
129 \param[in] qm - qmarray
130
131 \return qmo
132 ***************************************************************/
qmarray_copy(const struct Qmarray * qm)133 struct Qmarray * qmarray_copy(const struct Qmarray * qm)
134 {
135 struct Qmarray * qmo = qmarray_alloc(qm->nrows,qm->ncols);
136 size_t ii;
137 for (ii = 0; ii < qm->nrows*qm->ncols; ii++){
138 qmo->funcs[ii] = generic_function_copy(qm->funcs[ii]);
139 }
140
141 return qmo;
142 }
143
144 /***********************************************************//**
145 Perform an elementwise operator on a qmarray
146
147 \param[in] qm - qmarray
148 \param[in] op - operator
149
150 \return qmo
151 ***************************************************************/
qmarray_operate_elements(const struct Qmarray * qm,const struct Operator * op)152 struct Qmarray * qmarray_operate_elements(const struct Qmarray * qm,
153 const struct Operator * op)
154 {
155 struct Qmarray * qmo = qmarray_alloc(qm->nrows,qm->ncols);
156 if (op != NULL){
157 for (size_t ii = 0; ii < qm->nrows*qm->ncols; ii++){
158 qmo->funcs[ii] = op->f(qm->funcs[ii], op->opts);
159 }
160 }
161 else{
162 for (size_t ii = 0; ii < qm->nrows*qm->ncols; ii++){
163 qmo->funcs[ii] = generic_function_copy(qm->funcs[ii]);
164 }
165 }
166
167 return qmo;
168 }
169
170 /***********************************************************//**
171 Serialize a qmarray
172
173 \param[in,out] ser - stream to serialize to
174 \param[in] qma - quasimatrix array
175 \param[in,out] totSizeIn - if NULL then serialize, if not NULL
176 then return size;
177
178 \return ser shifted by number of bytes
179 ***************************************************************/
180 unsigned char *
qmarray_serialize(unsigned char * ser,const struct Qmarray * qma,size_t * totSizeIn)181 qmarray_serialize(unsigned char * ser, const struct Qmarray * qma,
182 size_t *totSizeIn)
183 {
184
185 // nrows -> ncols -> func -> func-> ... -> func
186 size_t ii;
187 size_t totSize = 2*sizeof(size_t);
188 size_t size_temp;
189 for (ii = 0; ii < qma->nrows * qma->ncols; ii++){
190 serialize_generic_function(NULL,qma->funcs[ii],&size_temp);
191 totSize += size_temp;
192 }
193 if (totSizeIn != NULL){
194 *totSizeIn = totSize;
195 return ser;
196 }
197
198 //printf("in serializing qmarray\n");
199 unsigned char * ptr = ser;
200 ptr = serialize_size_t(ptr, qma->nrows);
201 ptr = serialize_size_t(ptr, qma->ncols);
202 for (ii = 0; ii < qma->nrows * qma->ncols; ii++){
203 /* printf("on function number (%zu/%zu)\n",ii+1,qma->nrows * qma->ncols); */
204 // print_generic_function(qma->funcs[ii],3,NULL);
205 ptr = serialize_generic_function(ptr, qma->funcs[ii],NULL);
206 }
207 return ptr;
208 }
209
210 /***********************************************************//**
211 Deserialize a qmarray
212
213 \param[in] ser - serialized qmarray
214 \param[in,out] qma - qmarray
215
216 \return ptr - shifted ser after deserialization
217 ***************************************************************/
218 unsigned char *
qmarray_deserialize(unsigned char * ser,struct Qmarray ** qma)219 qmarray_deserialize(unsigned char * ser, struct Qmarray ** qma)
220 {
221 unsigned char * ptr = ser;
222
223 size_t nrows, ncols;
224 ptr = deserialize_size_t(ptr, &nrows);
225 ptr = deserialize_size_t(ptr, &ncols);
226 /* printf("deserialized rows,cols = (%zu,%zu)\n",nrows,ncols); */
227
228 *qma = qmarray_alloc(nrows, ncols);
229
230 size_t ii;
231 for (ii = 0; ii < nrows * ncols; ii++){
232 ptr = deserialize_generic_function(ptr, &((*qma)->funcs[ii]));
233 }
234
235 return ptr;
236 }
237
238 /***********************************************************//**
239 Save a qmarray in text format
240
241 \param[in] qma - quasimatrix array
242 \param[in,out] fp - stream to which to write
243 \param[in] prec - precision with which to print
244
245 ***************************************************************/
qmarray_savetxt(const struct Qmarray * qma,FILE * fp,size_t prec)246 void qmarray_savetxt(const struct Qmarray * qma, FILE * fp, size_t prec)
247 {
248
249 // nrows -> ncols -> func -> func-> ... -> func
250
251 assert (qma != NULL);
252 fprintf(fp,"%zu ",qma->nrows);
253 fprintf(fp,"%zu ",qma->ncols);
254 for (size_t ii = 0; ii < qma->nrows * qma->ncols; ii++){
255 generic_function_savetxt(qma->funcs[ii],fp,prec);
256 }
257 }
258
259 /***********************************************************//**
260 Load a qmarray saved in text format
261
262 \param[in,out] fp - stream to save to
263
264 \return qmarray
265 ***************************************************************/
qmarray_loadtxt(FILE * fp)266 struct Qmarray * qmarray_loadtxt(FILE * fp)
267 {
268 size_t nrows, ncols;
269 int num = fscanf(fp,"%zu ",&nrows);
270 assert (num == 1);
271 num =fscanf(fp,"%zu ",&ncols);
272 assert(num == 1);
273 struct Qmarray * qma = qmarray_alloc(nrows, ncols);
274
275 size_t ii;
276 for (ii = 0; ii < nrows * ncols; ii++){
277 qma->funcs[ii] = generic_function_loadtxt(fp);
278 }
279
280 return qma;
281 }
282
283 /***********************************************************//**
284 Create a qmarray by approximating 1d functions
285
286 \param[in] nrows - number of rows of quasimatrix
287 \param[in] ncols - number of cols of quasimatrix
288 \param[in] fapp - approximation arguments
289 \param[in] fw - wrapped function
290
291 \return qmarray
292 ***************************************************************/
293 struct Qmarray *
qmarray_approx1d(size_t nrows,size_t ncols,struct OneApproxOpts * fapp,struct Fwrap * fw)294 qmarray_approx1d(size_t nrows, size_t ncols, struct OneApproxOpts * fapp,
295 struct Fwrap * fw)
296 /* double (**funcs)(double, void *), */
297 /* void ** args, */
298 /* enum function_class fc, void * st, double lb, */
299 /* double ub, void * aopts) */
300 {
301 struct Qmarray * qm = qmarray_alloc(nrows,ncols);
302 size_t ii;
303 for (ii = 0; ii < nrows*ncols; ii++){
304 fwrap_set_which_eval(fw,ii);
305 qm->funcs[ii] =
306 generic_function_approximate1d(fapp->fc,fapp->aopts,fw);
307 }
308 return qm;
309 }
310
311 /***********************************************************
312 Create a qmarray from a fiber_cuts array
313
314 \param[in] nrows - number of rows of qmarray
315 \param[in] ncols - number of columns of qmarray
316 \param[in] f - functions
317 \param[in] fcut - array of fiber cuts
318 \param[in] fc - function class of each column
319 \param[in] sub_type - sub_type of each column
320 \param[in] lb - lower bound of inputs to functions
321 \param[in] ub - upper bound of inputs to functions
322 \param[in] aopts - approximation options
323
324 \return qmarray
325 ***************************************************************/
326 /* struct Qmarray * */
327 /* qmarray_from_fiber_cuts(size_t nrows, size_t ncols, */
328 /* double (*f)(double, void *),struct FiberCut ** fcut, */
329 /* enum function_class fc, void * sub_type, double lb, */
330 /* double ub, void * aopts) */
331 /* { */
332 /* struct Qmarray * qm = qmarray_alloc(nrows,ncols); */
333 /* size_t ii; */
334 /* for (ii = 0; ii < nrows*ncols; ii++){ */
335 /* qm->funcs[ii] = generic_function_approximate1d(f, fcut[ii], */
336 /* fc, sub_type, lb, ub, aopts); */
337 /* } */
338 /* return qm; */
339 /* } */
340
341
342
343 /***********************************************************//**
344 Generate a qmarray with orthonormal columns consisting of
345 one dimensional functions of linear elements on a certain grid
346
347 \param[in] nrows - number of rows
348 \param[in] ncols - number of columns
349 \param[in] grid - grid size
350
351 \return qmarray with orthonormal columns
352
353 \note
354 - Not super efficient because of copies
355 ***************************************************************/
356 /* struct Qmarray * */
357 /* qmarray_orth1d_linelm_grid(size_t nrows,size_t ncols, struct c3Vector * grid) */
358 /* { */
359 /* struct Qmarray * qm = qmarray_alloc(nrows,ncols); */
360 /* generic_function_array_orth1d_linelm_columns(qm->funcs,nrows,ncols,grid); */
361 /* return qm; */
362 /* } */
363
364 /***********************************************************//**
365 Generate a qmarray with orthonormal rows
366
367 \param[in] nrows - number of rows
368 \param[in] ncols - number of columns
369 \param[in] opts - approximation optsion
370
371 \return qmarray with orthonormal rows
372
373 \note
374 Not super efficient because of copies
375 ***************************************************************/
376 struct Qmarray *
qmarray_orth1d_rows(size_t nrows,size_t ncols,struct OneApproxOpts * opts)377 qmarray_orth1d_rows(size_t nrows, size_t ncols, struct OneApproxOpts * opts)
378 {
379
380 struct Qmarray * qm = qmarray_alloc(nrows,ncols);
381
382 struct GenericFunction ** funcs = NULL;
383 if ( NULL == (funcs = malloc(nrows*sizeof(struct GenericFunction *)))){
384 fprintf(stderr, "failed to allocate memory for quasimatrix.\n");
385 exit(1);
386 }
387 size_t ii, jj,kk;
388 for (ii = 0; ii < nrows; ii++){
389 funcs[ii] = NULL;
390 }
391 /* printf("wwwhat\n"); */
392 generic_function_array_orth(nrows, funcs, opts->fc, opts->aopts);
393 /* printf("wwwhere\n"); */
394 /* for (size_t ii = 0; ii < nrows; ii++){ */
395 /* printf("after array_orth ii = %zu\n", ii); */
396 /* print_generic_function(funcs[ii], 4, NULL, stdout); */
397 /* } */
398
399 struct GenericFunction * zero = generic_function_constant(0.0,opts->fc,opts->aopts);
400
401 size_t onnon = 0;
402 size_t onorder = 0;
403 for (jj = 0; jj < nrows; jj++){
404 assert (onorder < nrows);
405 qm->funcs[onnon*nrows+jj] = generic_function_copy(funcs[onorder]);
406 double normfunc = generic_function_norm(qm->funcs[onnon*nrows+jj]);
407 /* printf("normfunc rows = %3.15G\n", normfunc); */
408 if (isnan(normfunc)){
409 printf("normfunc rows = %3.15G\n", normfunc);
410 printf("num_rows = %zu\n", nrows);
411 printf("num_cols = %zu\n", ncols);
412 printf("onorder = %zu\n", onorder);
413 print_generic_function(qm->funcs[onnon*nrows+jj], 4, NULL, stdout);
414 exit(1);
415 }
416 for (kk = 0; kk < onnon; kk++){
417 qm->funcs[kk*nrows+jj] = generic_function_copy(zero);
418 }
419 for (kk = onnon+1; kk < ncols; kk++){
420 qm->funcs[kk*nrows+jj] = generic_function_copy(zero);
421 }
422 onnon = onnon+1;
423 if (onnon == ncols){
424 onorder = onorder+1;
425 onnon = 0;
426 }
427 }
428 //printf("max order rows = %zu\n",onorder);
429
430 for (ii = 0; ii < nrows;ii++){
431 generic_function_free(funcs[ii]);
432 funcs[ii] = NULL;
433 }
434 free(funcs); funcs=NULL;
435 generic_function_free(zero); zero = NULL;
436
437 return qm;
438 }
439
440 /***********************************************************//**
441 Generate a qmarray with orthonormal columns consisting of
442 one dimensional functions
443
444 \param[in] nrows - number of rows
445 \param[in] ncols - number of columns
446 \param[in] opts - options
447
448 \return qmarray with orthonormal columns
449
450 \note
451 - Not super efficient because of copies
452 ***************************************************************/
453 struct Qmarray *
qmarray_orth1d_columns(size_t nrows,size_t ncols,struct OneApproxOpts * opts)454 qmarray_orth1d_columns(size_t nrows,size_t ncols,struct OneApproxOpts * opts)
455 {
456
457 struct Qmarray * qm = qmarray_alloc(nrows,ncols);
458
459 struct GenericFunction ** funcs = NULL;
460 if ( NULL == (funcs = malloc(ncols*sizeof(struct GenericFunction *)))){
461 fprintf(stderr, "failed to allocate memory for quasimatrix.\n");
462 exit(1);
463 }
464 /* size_t ii, jj; //,kk; */
465 for (size_t ii = 0; ii < ncols; ii++){
466 funcs[ii] = NULL;
467 }
468
469 struct GenericFunction * zero = generic_function_constant(0.0,opts->fc,opts->aopts);
470 generic_function_array_orth(ncols,funcs,opts->fc,opts->aopts);
471 /* if (nrows >= ncols){ */
472 /* generic_function_array_orth(1,funcs,opts->fc,opts->aopts); */
473 /* for (size_t col = 0; col < ncols; col++){ */
474 /* for (size_t row = 0; row < nrows; row++){ */
475 /* if (row != col){ */
476 /* qm->funcs[col*nrows+row] = generic_function_copy(zero); */
477 /* } */
478 /* else{ */
479 /* qm->funcs[col*nrows+row] = generic_function_copy(funcs[0]); */
480 /* } */
481
482 /* } */
483 /* } */
484 /* } */
485 /* else{ */
486
487 /* size_t ngen = ncols / nrows + 1; */
488 /* generic_function_array_orth(ngen+1,funcs,opts->fc,opts->aopts); */
489 /* size_t oncol = 0; */
490 /* /\* int minnum = 0; *\/ */
491 /* size_t maxnum = 0; */
492 /* int converged = 0; */
493 /* printf("nrows=%zu,ncols=%zu,ngen=%zu\n",nrows,ncols,ngen); */
494 /* while (converged == 0){ */
495 /* size_t count = 0; */
496 /* printf("maxnum = %zu\n",maxnum); */
497 /* /\* print_generic_function(funcs[maxnum],3,NULL); *\/ */
498 /* for (size_t col = oncol; col < oncol+nrows; col++){ */
499 /* printf("col = %zu/%zu\n",col,ncols); */
500 /* if (col == ncols){ */
501 /* printf("conveged!\n"); */
502 /* converged = 1; */
503 /* break; */
504 /* } */
505 /* for (size_t row = 0; row < nrows; row++){ */
506 /* printf("row=%zu/%zu\n",row,nrows); */
507 /* if (row != count){ */
508 /* qm->funcs[col*nrows+row] = generic_function_copy(zero); */
509 /* } */
510 /* else{ */
511 /* qm->funcs[col*nrows+row] = generic_function_copy(funcs[maxnum]); */
512 /* } */
513 /* } */
514 /* count++; */
515 /* } */
516 /* maxnum += 1; */
517 /* oncol += nrows; */
518 /* if (maxnum > ngen){ */
519 /* fprintf(stderr,"Cannot generate enough orthonormal polynomials\n"); */
520 /* assert (1 == 0); */
521 /* } */
522 /* } */
523 /* printf("we are out\n"); */
524 /* } */
525 size_t onnon = 0;
526 size_t onorder = 0;
527 for (size_t jj = 0; jj < ncols; jj++){
528 qm->funcs[jj*nrows+onnon] = generic_function_copy(funcs[onorder]);
529 double normfunc = generic_function_norm(qm->funcs[jj*nrows+onnon]);
530 /* printf("normfunc cols = %3.15G %d \n", normfunc, isnan(normfunc)); */
531 if (isnan(normfunc)){
532 printf("normfunc cols = %3.15G\n", normfunc);
533 printf("num_rows = %zu\n", nrows);
534 printf("num_cols = %zu\n", ncols);
535 print_generic_function(qm->funcs[jj*nrows+onnon], 4, NULL, stdout);
536
537 generic_function_free(qm->funcs[jj*nrows+onnon]);
538 qm->funcs[jj*nrows+onnon] = generic_function_copy(zero);
539 exit(1);
540 }
541
542 for (size_t kk = 0; kk < onnon; kk++){
543 qm->funcs[jj*nrows+kk] = generic_function_copy(zero);
544 }
545 for (size_t kk = onnon+1; kk < nrows; kk++){
546 qm->funcs[jj*nrows+kk] = generic_function_copy(zero);
547 }
548 onnon = onnon+1;
549 if (onnon == nrows){
550 onorder = onorder+1;
551 onnon = 0;
552 }
553 }
554
555 //printf("max order rows = %zu\n",onorder);
556
557 for (size_t ii = 0; ii < ncols; ii++){
558 if (funcs[ii] != NULL){
559 generic_function_free(funcs[ii]);
560 funcs[ii] = NULL;
561 }
562 }
563 free(funcs); funcs=NULL;
564 generic_function_free(zero); zero = NULL;
565 /* printf("returned\n"); */
566 /* return qm; */
567 /* struct Qmarray * qm = qmarray_alloc(nrows,ncols); */
568 /* struct Qmarray * qmtemp = qmarray_alloc(ncols,1); */
569 /* generic_function_array_orth1d_columns(qm->funcs,qmtemp->funcs, */
570 /* fc,st,nrows,ncols,lb,ub); */
571
572 /* qmarray_free(qmtemp); qmtemp = NULL; */
573 return qm;
574 }
575
576 /***********************************************************//**
577 Create a Qmarray of zero functions
578
579 \param[in] nrows - number of rows of quasimatrix
580 \param[in] ncols - number of cols of quasimatrix
581 \param[in] opts - options
582
583
584 \return qmarray
585 ***************************************************************/
586 struct Qmarray *
qmarray_zeros(size_t nrows,size_t ncols,struct OneApproxOpts * opts)587 qmarray_zeros(size_t nrows,size_t ncols, struct OneApproxOpts * opts)
588 {
589 struct Qmarray * qm = qmarray_alloc(nrows,ncols);
590 size_t ii;
591 for (ii = 0; ii < nrows*ncols; ii++){
592 qm->funcs[ii] = generic_function_constant(0.0,opts->fc,opts->aopts);
593 }
594 return qm;
595 }
596
597 /********************************************************//**
598 * Create a qmarray consisting of pseudo-random orth poly expansion
599 *
600 * \param[in] ptype - polynomial type
601 * \param[in] nrows - number of rows
602 * \param[in] ncols - number of columns
603 * \param[in] maxorder - maximum order of the polynomial
604 * \param[in] lower - lower bound of input
605 * \param[in] upper - upper bound of input
606 *
607 * \return qm - qmarray
608 ************************************************************/
609 struct Qmarray *
qmarray_poly_randu(enum poly_type ptype,size_t nrows,size_t ncols,size_t maxorder,double lower,double upper)610 qmarray_poly_randu(enum poly_type ptype,
611 size_t nrows, size_t ncols,
612 size_t maxorder, double lower, double upper)
613 {
614 struct Qmarray * qm = qmarray_alloc(nrows,ncols);
615 size_t ii;
616 for (ii = 0; ii < nrows*ncols; ii++){
617 qm->funcs[ii] = generic_function_poly_randu(ptype, maxorder,
618 lower,upper);
619 }
620 return qm;
621 }
622
623
624
625 // getters and setters
626 /**********************************************************//**
627 Return a typed function
628 **************************************************************/
629 void *
qmarray_get_func_base(const struct Qmarray * qma,size_t row,size_t col)630 qmarray_get_func_base(const struct Qmarray * qma, size_t row, size_t col)
631 {
632 assert (qma != NULL);
633 assert (row < qma->nrows);
634 assert (col < qma->ncols);
635
636 return qma->funcs[col*qma->nrows+row]->f;
637 }
638
639 /**********************************************************//**
640 Return a function
641 **************************************************************/
642 struct GenericFunction *
qmarray_get_func(const struct Qmarray * qma,size_t row,size_t col)643 qmarray_get_func(const struct Qmarray * qma, size_t row, size_t col)
644 {
645 assert (qma != NULL);
646 assert (row < qma->nrows);
647 assert (col < qma->ncols);
648
649 return qma->funcs[col*qma->nrows+row];
650 }
651
652 /**********************************************************//**
653 Set the row of a quasimatrix array to a given quasimatrix by copying
654
655 \param[in,out] qma - qmarray
656 \param[in] row - row to set
657 \param[in] qm - quasimatrix to copy
658 **************************************************************/
659 void
qmarray_set_row(struct Qmarray * qma,size_t row,const struct Quasimatrix * qm)660 qmarray_set_row(struct Qmarray * qma, size_t row, const struct Quasimatrix * qm)
661 {
662 size_t ii;
663 for (ii = 0; ii < qma->ncols; ii++){
664 generic_function_free(qma->funcs[ii*qma->nrows+row]);
665 struct GenericFunction * f = quasimatrix_get_func(qm,ii);
666 qma->funcs[ii*qma->nrows+row] = generic_function_copy(f);
667 }
668 }
669
670 /**********************************************************//**
671 Extract a copy of a quasimatrix from a row of a qmarray
672
673 \param[in] qma - qmarray
674 \param[in] row - row to copy
675
676 \return qm - quasimatrix
677 **************************************************************/
678 struct Quasimatrix *
qmarray_extract_row(const struct Qmarray * qma,size_t row)679 qmarray_extract_row(const struct Qmarray * qma, size_t row)
680 {
681
682 struct Quasimatrix * qm = quasimatrix_alloc(qma->ncols);
683 size_t ii;
684 for (ii = 0; ii < qma->ncols; ii++){
685 quasimatrix_set_func(qm,qma->funcs[ii*qma->nrows+row],ii);
686 /* qm->funcs[ii] = generic_function_copy(qma->funcs[ii*qma->nrows+row]); */
687 }
688 return qm;
689 }
690
691
692 /**********************************************************//**
693 Set the column of a quasimatrix array to a given quasimatrix by copying
694
695 \param[in,out] qma - qmarray to modify
696 \param[in] col - column to set
697 \param[in] qm - quasimatrix to copy
698 **************************************************************/
qmarray_set_column(struct Qmarray * qma,size_t col,const struct Quasimatrix * qm)699 void qmarray_set_column(struct Qmarray * qma, size_t col,
700 const struct Quasimatrix * qm)
701 {
702 size_t ii;
703 for (ii = 0; ii < qma->nrows; ii++){
704 generic_function_free(qma->funcs[col*qma->nrows+ii]);
705 struct GenericFunction * f = quasimatrix_get_func(qm,ii);
706 qma->funcs[col*qma->nrows+ii] = generic_function_copy(f);
707 }
708 }
709
710 /***********************************************************//**
711 Function qmarray_set_column_gf
712
713 Purpose: Set the column of a quasimatrix array
714 to a given array of generic functions
715
716 \param[in,out] qma - qmarray
717 \param[in] col - column to set
718 \param[in] gf - array of generic functions
719 ***************************************************************/
720 void
qmarray_set_column_gf(struct Qmarray * qma,size_t col,struct GenericFunction ** gf)721 qmarray_set_column_gf(struct Qmarray * qma, size_t col,
722 struct GenericFunction ** gf)
723 {
724 size_t ii;
725 for (ii = 0; ii < qma->nrows; ii++){
726 generic_function_free(qma->funcs[col*qma->nrows+ii]);
727 qma->funcs[col*qma->nrows+ii] = generic_function_copy(gf[ii]);
728 }
729 }
730
731 /*********************************************************//**
732 Extract a copy of a quasimatrix from a column of a qmarray
733
734 \param[in] qma - qmarray
735 \param[in] col - column to copy
736
737 \return quasimatrix
738 *************************************************************/
739 struct Quasimatrix *
qmarray_extract_column(const struct Qmarray * qma,size_t col)740 qmarray_extract_column(const struct Qmarray * qma, size_t col)
741 {
742
743 struct Quasimatrix * qm = quasimatrix_alloc(qma->nrows);
744 size_t ii;
745 for (ii = 0; ii < qma->nrows; ii++){
746 quasimatrix_set_func(qm,qma->funcs[col*qma->nrows+ii],ii);
747 /* qm->funcs[ii] = generic_function_copy(qma->funcs[col*qma->nrows+ii]); */
748 }
749 return qm;
750 }
751
752 /***********************************************************//**
753 get number of columns
754 ***************************************************************/
qmarray_get_ncols(const struct Qmarray * qma)755 size_t qmarray_get_ncols(const struct Qmarray * qma)
756 {
757 assert (qma != NULL);
758 return qma->ncols;
759 }
760
761 /***********************************************************//**
762 get number of rows
763 ***************************************************************/
qmarray_get_nrows(const struct Qmarray * qma)764 size_t qmarray_get_nrows(const struct Qmarray * qma)
765 {
766 assert (qma != NULL);
767 return qma->nrows;
768 }
769
770 //////////////////////////////////////////////////////////////////
771
772 /***********************************************************//**
773 Quasimatrix array - matrix multiplication
774
775 \param[in] Q - quasimatrix array
776 \param[in] R - matrix (fortran order)
777 \param[in] b - number of columns of R
778
779 \return B - qmarray
780 ***************************************************************/
qmam(const struct Qmarray * Q,const double * R,size_t b)781 struct Qmarray * qmam(const struct Qmarray * Q, const double * R, size_t b)
782 {
783 size_t nrows = Q->nrows;
784 struct Qmarray * B = qmarray_alloc(nrows,b);
785 size_t ii,jj;
786 for (jj = 0; jj < nrows; jj++){
787 for (ii = 0; ii < b; ii++){
788 B->funcs[ii*nrows+jj] = // Q[jj,:]B[:,ii]
789 generic_function_lin_comb2(Q->ncols,
790 Q->nrows, Q->funcs + jj, 1, R + ii*Q->ncols);
791 }
792 }
793 return B;
794 }
795
796 /***********************************************************//**
797 Transpose Quasimatrix array - matrix multiplication
798
799 \param[in] Q - quasimatrix array
800 \param[in] R - matrix (fortran order)
801 \param[in] b - number of columns of R
802
803 \return B - qmarray
804 ***************************************************************/
805 struct Qmarray *
qmatm(const struct Qmarray * Q,const double * R,size_t b)806 qmatm(const struct Qmarray * Q,const double * R, size_t b)
807 {
808 struct Qmarray * B = qmarray_alloc(Q->ncols,b);
809 size_t ii,jj;
810 for (jj = 0; jj < B->nrows; jj++){
811 for (ii = 0; ii < b; ii++){
812 B->funcs[ii*B->nrows+jj] = // Q[:,jj]^T B[:,ii]
813 generic_function_lin_comb2(Q->nrows, 1,
814 Q->funcs + jj * Q->nrows,
815 1, R + ii*Q->nrows);
816 }
817 }
818 return B;
819 }
820
821 /***********************************************************//**
822 Matrix - Quasimatrix array multiplication
823
824 \param[in] R - matrix (fortran order)
825 \param[in] Q - quasimatrix array
826 \param[in] b - number of rows of R
827
828 \return B - qmarray (b x Q->ncols)
829 ***************************************************************/
830 struct Qmarray *
mqma(const double * R,const struct Qmarray * Q,size_t b)831 mqma(const double * R, const struct Qmarray * Q, size_t b)
832 {
833 struct Qmarray * B = qmarray_alloc(b, Q->ncols);
834 size_t ii,jj;
835 for (jj = 0; jj < Q->ncols; jj++){
836 for (ii = 0; ii < b; ii++){
837 B->funcs[jj*b+ii] = // R[ii,:]Q[:,jj]
838 generic_function_lin_comb2(Q->nrows, 1, Q->funcs + jj*Q->nrows,
839 b, R + ii);
840 }
841 }
842 return B;
843 }
844
845 /***********************************************************//**
846 Qmarray - Qmarray mutliplication
847
848 \param[in] a - qmarray 1
849 \param[in] b - qmarray 2
850
851 \return c - qmarray
852 ***************************************************************/
853 struct Qmarray *
qmaqma(const struct Qmarray * a,const struct Qmarray * b)854 qmaqma(const struct Qmarray * a, const struct Qmarray * b)
855 {
856 struct Qmarray * c = qmarray_alloc(a->nrows, b->ncols);
857 size_t ii,jj;
858 for (jj = 0; jj < b->ncols; jj++){
859 for (ii = 0; ii < a->nrows; ii++){
860 // a[ii,:]b[:,jj]
861 c->funcs[jj*c->nrows+ii] =
862 generic_function_sum_prod(a->ncols, a->nrows,
863 a->funcs + ii, 1,
864 b->funcs + jj*b->nrows);
865 }
866 }
867 return c;
868 }
869
870 /***********************************************************//**
871 Transpose qmarray - qmarray mutliplication
872
873 \param[in] a - qmarray 1
874 \param[in] b - qmarray 2
875
876 \return c - qmarray : c = a^T b
877 ***************************************************************/
878 struct Qmarray *
qmatqma(const struct Qmarray * a,const struct Qmarray * b)879 qmatqma(const struct Qmarray * a, const struct Qmarray * b)
880 {
881 struct Qmarray * c = qmarray_alloc(a->ncols, b->ncols);
882 size_t ii,jj;
883 for (jj = 0; jj < b->ncols; jj++){
884 for (ii = 0; ii < a->ncols; ii++){
885 //printf("computing c[%zu,%zu] \n",ii,jj);
886 // c[jj*a->ncols+ii] = a[:,ii]^T b[:,jj]
887 c->funcs[jj*c->nrows+ii] = generic_function_sum_prod(b->nrows, 1,
888 a->funcs + ii*a->nrows, 1, b->funcs + jj*b->nrows);
889 }
890 }
891 return c;
892 }
893
894 /***********************************************************//**
895 qmarray - transpose qmarray mutliplication
896
897 \param[in] a - qmarray 1
898 \param[in] b - qmarray 2
899
900 \return c - qmarray : c = a b^T
901 ***************************************************************/
902 struct Qmarray *
qmaqmat(const struct Qmarray * a,const struct Qmarray * b)903 qmaqmat(const struct Qmarray * a, const struct Qmarray * b)
904 {
905 struct Qmarray * c = qmarray_alloc(a->nrows, b->nrows);
906 size_t ii,jj;
907 for (jj = 0; jj < b->nrows; jj++){
908 for (ii = 0; ii < a->nrows; ii++){
909 // c[jj*c->ncols+ii] = a[ii,:]^T b[jj,:]^T
910 c->funcs[jj*c->nrows+ii] = generic_function_sum_prod(a->ncols, a->nrows,
911 a->funcs + ii, b->nrows, b->funcs + jj);
912 }
913 }
914 return c;
915 }
916
917
918 /***********************************************************//**
919 Transpose qmarray - transpose qmarray mutliplication
920
921 \param[in] a - qmarray 1
922 \param[in] b - qmarray 2
923
924 \return c - qmarray
925 ***************************************************************/
926 struct Qmarray *
qmatqmat(const struct Qmarray * a,const struct Qmarray * b)927 qmatqmat(const struct Qmarray * a, const struct Qmarray * b)
928 {
929 struct Qmarray * c = qmarray_alloc(a->ncols, b->nrows);
930 size_t ii,jj;
931 for (jj = 0; jj < b->nrows; jj++){
932 for (ii = 0; ii < a->ncols; ii++){
933 // c[jj*a->ncols+ii] = a[:,ii]^T b[jj,:]
934 c->funcs[ii*c->nrows+jj] = generic_function_sum_prod(b->ncols, 1,
935 a->funcs + ii*a->nrows, b->nrows, b->funcs + jj);
936 }
937 }
938 return c;
939 }
940
941 /***********************************************************//**
942 Integral of Transpose qmarray - qmarray mutliplication
943
944 \param[in] a - qmarray 1
945 \param[in] b - qmarray 2
946
947 \return array of integrals of \f$ c = \int a(x)^T b(x) dx \f$
948 ***************************************************************/
949 double *
qmatqma_integrate(const struct Qmarray * a,const struct Qmarray * b)950 qmatqma_integrate(const struct Qmarray * a,const struct Qmarray * b)
951 {
952 double * c = calloc_double(a->ncols*b->ncols);
953 size_t nrows = a->ncols;
954 size_t ncols = b->ncols;
955 size_t ii,jj;
956 for (jj = 0; jj < ncols; jj++){
957 for (ii = 0; ii < nrows; ii++){
958 // c[jj*nrows+ii] = a[:,ii]^T b[jj,:]
959 c[jj*nrows+ii] = generic_function_inner_sum(b->nrows, 1,
960 a->funcs + ii*a->nrows, 1, b->funcs + jj*b->nrows);
961 }
962 }
963 return c;
964 }
965
966 /***********************************************************//**
967 Integral of qmarray - transpose qmarray mutliplication
968
969 \param[in] a - qmarray 1
970 \param[in] b - qmarray 2
971
972 \return array of integrals of \f$ c = \int a(x) b(x)^T dx \f$
973 ***************************************************************/
974 double *
qmaqmat_integrate(const struct Qmarray * a,const struct Qmarray * b)975 qmaqmat_integrate(const struct Qmarray * a,const struct Qmarray * b)
976 {
977 double * c = calloc_double(a->nrows*b->nrows);
978 size_t nrows = a->nrows;
979 size_t ncols = b->nrows;
980 size_t ii,jj;
981 for (jj = 0; jj < ncols; jj++){
982 for (ii = 0; ii < nrows; ii++){
983 // c[jj*nrows+ii] = a[:,ii]^T b[jj,:]
984 c[jj*nrows+ii] =
985 generic_function_inner_sum(a->ncols,
986 a->nrows,
987 a->funcs + ii,
988 b->nrows, b->funcs + jj);
989 }
990 }
991 return c;
992 }
993
994 /***********************************************************//**
995 Integral of Transpose qmarray - transpose qmarray mutliplication
996
997 \param[in] a - qmarray 1
998 \param[in] b - qmarray 2
999
1000 \return array of integrals
1001 ***************************************************************/
1002 double *
qmatqmat_integrate(const struct Qmarray * a,const struct Qmarray * b)1003 qmatqmat_integrate(const struct Qmarray * a,const struct Qmarray * b)
1004 {
1005 double * c = calloc_double(a->ncols*b->nrows);
1006 size_t ii,jj;
1007 for (jj = 0; jj < b->nrows; jj++){
1008 for (ii = 0; ii < a->ncols; ii++){
1009 // c[jj*a->ncols+ii] = a[:,ii]^T b[jj,:]
1010 //c[ii*b->nrows+jj] = generic_function_sum_prod_integrate(b->ncols, 1,
1011 // a->funcs + ii*a->nrows, b->nrows, b->funcs + jj);
1012
1013 c[ii*b->nrows+jj] = generic_function_inner_sum(b->ncols, 1,
1014 a->funcs + ii*a->nrows, b->nrows, b->funcs + jj);
1015 }
1016 }
1017 return c;
1018 }
1019
1020 /***********************************************************//**
1021 (Weighted) Integral of Transpose qmarray - transpose qmarray mutliplication
1022
1023 \param[in] a - qmarray 1
1024 \param[in] b - qmarray 2
1025
1026 \return array of integral values
1027
1028 \note
1029 In order for this function to make sense a(x) and b(x)
1030 shoud have the same underlying basis
1031 ***************************************************************/
1032 double *
qmatqmat_integrate_weighted(const struct Qmarray * a,const struct Qmarray * b)1033 qmatqmat_integrate_weighted(const struct Qmarray * a,const struct Qmarray * b)
1034 {
1035 double * c = calloc_double(a->ncols*b->nrows);
1036 size_t ii,jj;
1037 for (jj = 0; jj < b->nrows; jj++){
1038 for (ii = 0; ii < a->ncols; ii++){
1039 // c[jj*a->ncols+ii] = a[:,ii]^T b[jj,:]
1040 //c[ii*b->nrows+jj] = generic_function_sum_prod_integrate(b->ncols, 1,
1041 // a->funcs + ii*a->nrows, b->nrows, b->funcs + jj);
1042
1043 c[ii*b->nrows+jj] =
1044 generic_function_inner_weighted_sum(b->ncols, 1,
1045 a->funcs + ii*a->nrows, b->nrows, b->funcs + jj);
1046 }
1047 }
1048 return c;
1049 }
1050
1051 /***********************************************************//**
1052 Kronecker product between two qmarrays
1053
1054 \param[in] a - qmarray 1
1055 \param[in] b - qmarray 2
1056
1057 \return c - qmarray kron(a,b)
1058 ***************************************************************/
qmarray_kron(const struct Qmarray * a,const struct Qmarray * b)1059 struct Qmarray * qmarray_kron(const struct Qmarray * a,const struct Qmarray * b)
1060 {
1061 struct Qmarray * c = qmarray_alloc(a->nrows * b->nrows,
1062 a->ncols * b->ncols);
1063 size_t ii,jj,kk,ll,column,row, totrows;
1064
1065 totrows = a->nrows * b->nrows;
1066
1067 for (ii = 0; ii < a->ncols; ii++){
1068 for (jj = 0; jj < a->nrows; jj++){
1069 for (kk = 0; kk < b->ncols; kk++){
1070 for (ll = 0; ll < b->nrows; ll++){
1071 column = ii*b->ncols+kk;
1072 row = jj*b->nrows + ll;
1073 c->funcs[column*totrows + row] = generic_function_prod(
1074 a->funcs[ii*a->nrows+jj],b->funcs[kk*b->nrows+ll]);
1075 }
1076 }
1077 }
1078 }
1079 return c;
1080 }
1081
1082 /***********************************************************//**
1083 Integral of kronecker product between two qmarrays
1084
1085 \param[in] a - qmarray 1
1086 \param[in] b - qmarray 2
1087
1088 \return c - \f$ \int kron(a(x),b(x)) dx \f$
1089 ***************************************************************/
qmarray_kron_integrate(const struct Qmarray * a,const struct Qmarray * b)1090 double * qmarray_kron_integrate(const struct Qmarray * a, const struct Qmarray * b)
1091 {
1092 double * c = calloc_double(a->nrows*b->nrows * a->ncols*b->ncols);
1093
1094 size_t ii,jj,kk,ll,column,row, totrows;
1095
1096 totrows = a->nrows * b->nrows;
1097
1098 for (ii = 0; ii < a->ncols; ii++){
1099 for (jj = 0; jj < a->nrows; jj++){
1100 for (kk = 0; kk < b->ncols; kk++){
1101 for (ll = 0; ll < b->nrows; ll++){
1102 column = ii*b->ncols+kk;
1103 row = jj*b->nrows + ll;
1104 c[column*totrows + row] = generic_function_inner(
1105 a->funcs[ii*a->nrows+jj],b->funcs[kk*b->nrows+ll]);
1106 }
1107 }
1108 }
1109 }
1110
1111 return c;
1112 }
1113
1114 /***********************************************************//**
1115 Weighted integral of kronecker product between two qmarrays
1116
1117 \param[in] a - qmarray 1
1118 \param[in] b - qmarray 2
1119
1120 \return c - \f$ \int kron(a(x),b(x)) w(x) dx \f$
1121
1122 \note
1123 In order for this function to make sense a(x) and b(x)
1124 shoud have the same underlying basis
1125 ***************************************************************/
qmarray_kron_integrate_weighted(const struct Qmarray * a,const struct Qmarray * b)1126 double * qmarray_kron_integrate_weighted(const struct Qmarray * a, const struct Qmarray * b)
1127 {
1128 double * c = calloc_double(a->nrows*b->nrows * a->ncols*b->ncols);
1129
1130 size_t ii,jj,kk,ll,column,row, totrows;
1131
1132 totrows = a->nrows * b->nrows;
1133
1134 for (ii = 0; ii < a->ncols; ii++){
1135 for (jj = 0; jj < a->nrows; jj++){
1136 for (kk = 0; kk < b->ncols; kk++){
1137 for (ll = 0; ll < b->nrows; ll++){
1138 column = ii*b->ncols+kk;
1139 row = jj*b->nrows + ll;
1140 c[column*totrows + row] = generic_function_inner_weighted(
1141 a->funcs[ii*a->nrows+jj],b->funcs[kk*b->nrows+ll]);
1142 }
1143 }
1144 }
1145 }
1146 return c;
1147 }
1148
1149
1150 /***********************************************************//**
1151 If a is vector, b and c are quasimatrices then computes
1152 a^T kron(b,c)
1153
1154 \param[in] a - vector,array (b->ncols * c->ncols);
1155 \param[in] b - qmarray 1
1156 \param[in] c - qmarray 2
1157
1158 \return d - qmarray b->ncols x (c->ncols)
1159 ***************************************************************/
1160 struct Qmarray *
qmarray_vec_kron(const double * a,const struct Qmarray * b,const struct Qmarray * c)1161 qmarray_vec_kron(const double * a,const struct Qmarray * b,const struct Qmarray * c)
1162 {
1163 struct Qmarray * d = NULL;
1164
1165 struct Qmarray * temp = qmatm(b,a,c->nrows); // b->ncols * c->ncols
1166 //printf("got qmatm\n");
1167 //d = qmaqma(temp,c);
1168 d = qmatqmat(c,temp);
1169 qmarray_free(temp); temp = NULL;
1170 return d;
1171 }
1172
1173 /***********************************************************//**
1174 If a is vector, b and c are quasimatrices then computes
1175 \f$ \int a^T kron(b(x),c(x)) dx \f$
1176
1177 \param[in] a - vector,array (b->nrows * c->nrows);
1178 \param[in] b - qmarray 1
1179 \param[in] c - qmarray 2
1180
1181 \return d - b->ncols x (c->ncols)
1182 ***************************************************************/
1183 double *
qmarray_vec_kron_integrate(const double * a,const struct Qmarray * b,const struct Qmarray * c)1184 qmarray_vec_kron_integrate(const double * a, const struct Qmarray * b,const struct Qmarray * c)
1185 {
1186 double * d = NULL;
1187
1188 struct Qmarray * temp = qmatm(b,a,c->nrows); // b->ncols * c->ncols
1189 //printf("got qmatm\n");
1190 //d = qmaqma(temp,c);
1191 d = qmatqmat_integrate(c,temp);
1192 qmarray_free(temp); temp = NULL;
1193 return d;
1194 }
1195
1196 /***********************************************************//**
1197 If a is vector, b and c are quasimatrices then computes
1198 \f$ \int a^T kron(b(x),c(x)) w(x) dx \f$
1199
1200 \param[in] a - vector,array (b->nrows * c->nrows);
1201 \param[in] b - qmarray 1
1202 \param[in] c - qmarray 2
1203
1204 \return d - b->ncols x (c->ncols)
1205
1206 \note
1207 In order for this function to make sense b(x) and c(x)
1208 shoud have the same underlying basis
1209 ***************************************************************/
1210 double *
qmarray_vec_kron_integrate_weighted(const double * a,const struct Qmarray * b,const struct Qmarray * c)1211 qmarray_vec_kron_integrate_weighted(const double * a, const struct Qmarray * b,const struct Qmarray * c)
1212 {
1213 double * d = NULL;
1214
1215 struct Qmarray * temp = qmatm(b,a,c->nrows); // b->ncols * c->ncols
1216 //printf("got qmatm\n");
1217 //d = qmaqma(temp,c);
1218 d = qmatqmat_integrate_weighted(c,temp);
1219 qmarray_free(temp); temp = NULL;
1220 return d;
1221 }
1222
1223 /***********************************************************//**
1224 If a is a matrix, b and c are qmarrays then computes
1225 a kron(b,c) fast
1226
1227 \param[in] r - number of rows of a
1228 \param[in] a - array (k, b->nrows * c->nrows);
1229 \param[in] b - qmarray 1
1230 \param[in] c - qmarray 2
1231
1232 \return d - qmarray (k, b->ncols * c->ncols)
1233
1234 \note
1235 Do not touch this!! Who knows how it happens to be right...
1236 but it is tested (dont touch the test either. and it seems to be right
1237 ***************************************************************/
1238 struct Qmarray *
qmarray_mat_kron(size_t r,const double * a,const struct Qmarray * b,const struct Qmarray * c)1239 qmarray_mat_kron(size_t r, const double * a, const struct Qmarray * b, const struct Qmarray * c)
1240 {
1241 struct Qmarray * temp = qmarray_alloc(r*b->nrows,c->ncols);
1242 generic_function_kronh(1, r, b->nrows, c->nrows, c->ncols, a,
1243 c->funcs,temp->funcs);
1244
1245 struct Qmarray * d = qmarray_alloc(r, b->ncols * c->ncols);
1246 generic_function_kronh2(1,r, c->ncols, b->nrows,b->ncols,
1247 b->funcs,temp->funcs, d->funcs);
1248
1249 qmarray_free(temp); temp = NULL;
1250 return d;
1251 }
1252
1253 /***********************************************************//**
1254 If a is a matrix, b and c are qmarrays then computes
1255 kron(b,c)a fast
1256
1257 \param[in] r - number of rows of a
1258 \param[in] a - array (b->ncols * c->ncols, a);
1259 \param[in] b - qmarray 1
1260 \param[in] c - qmarray 2
1261
1262 \return d - qmarray (b->nrows * c->nrows , k)
1263
1264 \note
1265 Do not touch this!! Who knows how it happens to be right...
1266 but it is tested (dont touch the test either. and it seems to be right
1267 ***************************************************************/
1268 struct Qmarray *
qmarray_kron_mat(size_t r,const double * a,const struct Qmarray * b,const struct Qmarray * c)1269 qmarray_kron_mat(size_t r, const double * a, const struct Qmarray * b,
1270 const struct Qmarray * c)
1271 {
1272 struct Qmarray * temp = qmarray_alloc(c->nrows,b->ncols*r);
1273 generic_function_kronh(0, r, b->ncols, c->nrows, c->ncols, a,
1274 c->funcs,temp->funcs);
1275
1276 struct Qmarray * d = qmarray_alloc(b->nrows * c->nrows,r);
1277 generic_function_kronh2(0,r, c->nrows, b->nrows,b->ncols,
1278 b->funcs,temp->funcs, d->funcs);
1279
1280 qmarray_free(temp); temp = NULL;
1281 return d;
1282 }
1283
1284 // The algorithms below really are not cache-aware and might be slow ...
qmarray_block_kron_mat(char type,int left,size_t nlblocks,struct Qmarray ** lblocks,struct Qmarray * right,size_t r,double * mat,struct Qmarray * out)1285 void qmarray_block_kron_mat(char type, int left, size_t nlblocks,
1286 struct Qmarray ** lblocks, struct Qmarray * right, size_t r,
1287 double * mat, struct Qmarray * out)
1288 {
1289 if (left == 1)
1290 {
1291
1292 size_t ii;
1293 if (type == 'D'){
1294 size_t totrows = 0;
1295 for (ii = 0; ii < nlblocks; ii++){
1296 totrows += lblocks[ii]->nrows;
1297 }
1298 struct GenericFunction ** gfarray =
1299 generic_function_array_alloc(r * right->ncols * totrows);
1300
1301 generic_function_kronh(1, r, totrows, right->nrows, right->ncols,
1302 mat, right->funcs, gfarray);
1303
1304
1305 size_t width = 0;
1306 size_t runrows = 0;
1307 size_t jj,kk,ll;
1308 for (ii = 0; ii < nlblocks; ii++){
1309 struct GenericFunction ** temp =
1310 generic_function_array_alloc(right->ncols * r * lblocks[ii]->nrows);
1311
1312 for (jj = 0; jj < right->ncols; jj++){
1313 for (kk = 0; kk < lblocks[ii]->nrows; kk++){
1314 for (ll = 0; ll < r; ll++){
1315 temp[ll + kk*r + jj * lblocks[ii]->nrows*r] =
1316 gfarray[runrows + ll + kk*r + jj*totrows*r];
1317 }
1318 }
1319 }
1320 generic_function_kronh2(1,r,right->ncols,
1321 lblocks[ii]->nrows, lblocks[ii]->ncols,
1322 lblocks[ii]->funcs, temp,
1323 out->funcs + width * out->nrows);
1324 width += lblocks[ii]->ncols * right->ncols;
1325 runrows += r * lblocks[ii]->nrows;
1326 free(temp); temp = NULL;
1327 }
1328 generic_function_array_free(gfarray, r * right->ncols * totrows);
1329 }
1330 else if (type == 'R'){
1331 size_t row1 = lblocks[0]->nrows;
1332 for (ii = 1; ii < nlblocks; ii++){
1333 assert( row1 == lblocks[ii]->nrows );
1334 }
1335
1336 struct GenericFunction ** gfarray =
1337 generic_function_array_alloc(r * right->ncols * row1);
1338
1339 generic_function_kronh(1,r, row1, right->nrows, right->ncols,
1340 mat, right->funcs, gfarray);
1341
1342
1343 size_t width = 0;
1344 for (ii = 0; ii < nlblocks; ii++){
1345 generic_function_kronh2(1,r, right->ncols, row1,
1346 lblocks[ii]->ncols, lblocks[ii]->funcs, gfarray,
1347 out->funcs + width * out->nrows);
1348 width += lblocks[ii]->ncols * right->ncols;
1349 }
1350 generic_function_array_free(gfarray,r * right->ncols * row1);
1351 }
1352 else if (type == 'C'){
1353
1354 size_t col1 = lblocks[0]->ncols;
1355 size_t totrows = 0;
1356 for (ii = 1; ii < nlblocks; ii++){
1357 assert( col1 == lblocks[ii]->ncols );
1358 totrows += lblocks[ii]->nrows;
1359 }
1360
1361 struct Qmarray * temp = qmarray_alloc(r, col1*nlblocks*right->ncols);
1362 qmarray_block_kron_mat('D',1,nlblocks,lblocks,right,r,mat,temp);
1363
1364 size_t jj,kk,ll;
1365 size_t startelem = 0;
1366 for (ii = 0; ii < col1; ii++){
1367 for (jj = 0; jj < right->ncols; jj++){
1368 for (kk = 0; kk < r; kk++){
1369 out->funcs[kk + jj*r + ii * r*right->ncols] =
1370 generic_function_copy(
1371 temp->funcs[startelem + kk + jj*r + ii * r*right->ncols]);
1372
1373 }
1374 }
1375 }
1376 startelem += r * right->ncols * col1;
1377 for (ll = 1; ll < nlblocks; ll++){
1378 for (ii = 0; ii < col1; ii++){
1379 for (jj = 0; jj < right->ncols; jj++){
1380 for (kk = 0; kk < r; kk++){
1381 generic_function_axpy(1.0,
1382 temp->funcs[startelem + kk + jj*r + ii * r*right->ncols],
1383 out->funcs[kk + jj*r + ii * r*right->ncols]);
1384
1385 }
1386 }
1387 }
1388 startelem += r * right->ncols * col1;
1389 }
1390
1391 qmarray_free(temp); temp = NULL;
1392
1393 }
1394 }
1395 else{
1396 size_t ii;
1397 if (type == 'D'){
1398 size_t totcols = 0;
1399 size_t totrows = 0;
1400 for (ii = 0; ii < nlblocks; ii++){
1401 totcols += lblocks[ii]->ncols;
1402 totrows += lblocks[ii]->nrows;
1403 }
1404
1405
1406 struct GenericFunction ** gfarray =
1407 generic_function_array_alloc(r * right->nrows * totcols);
1408
1409 generic_function_kronh(0, r, totcols, right->nrows, right->ncols,
1410 mat, right->funcs, gfarray);
1411
1412 // printf("number of elements in gfarray is %zu\n",r * right->nrows * totcols);
1413 size_t runcols = 0;
1414 size_t runrows = 0;
1415 size_t jj,kk,ll;
1416 // printf("done with computing kronh1\n");
1417 for (ii = 0; ii < nlblocks; ii++){
1418
1419 size_t st1a = right->nrows;
1420 size_t st1b = lblocks[ii]->ncols;
1421 size_t st = st1a * st1b * r;
1422 struct GenericFunction ** temp = generic_function_array_alloc(st);
1423
1424 size_t rga = right->nrows;;
1425 size_t lastind = right->nrows * totcols;
1426
1427 for (jj = 0; jj < r; jj++){
1428 for (kk = 0; kk < st1b; kk++){
1429 for (ll = 0; ll < st1a; ll++){
1430 temp [ll + kk*st1a + jj*st1a*st1b] =
1431 gfarray[ll + (kk+runcols)*rga + jj*lastind];
1432 }
1433 }
1434 }
1435
1436 size_t st2a = right->nrows;
1437 size_t st2b = lblocks[ii]->nrows;
1438 size_t st2 = st2a * st2b * r;
1439
1440 struct GenericFunction ** temp2 =
1441 generic_function_array_alloc(st2);
1442
1443 generic_function_kronh2(0,r,st2a, st2b, st1b,
1444 lblocks[ii]->funcs, temp, temp2);
1445
1446 for (jj = 0; jj < st2a; jj++){
1447 for (kk = 0; kk < st2b; kk++){
1448 for (ll = 0; ll < r; ll++){
1449 out->funcs[jj + (runrows + kk)*st2a + ll*totrows* st2a] =
1450 temp2[jj + kk*st2a + ll*st2a*st2b];
1451 }
1452 }
1453 }
1454
1455 runcols += lblocks[ii]->ncols;
1456 runrows += lblocks[ii]->nrows;
1457 free(temp); temp = NULL;
1458 free(temp2); temp2 = NULL;
1459 }
1460 generic_function_array_free(gfarray, r*right->nrows*totcols);
1461 }
1462 else if (type == 'R')
1463 {
1464
1465 size_t totcols = lblocks[0]->ncols;
1466 size_t row1 = lblocks[0]->nrows;
1467 for (ii = 1; ii < nlblocks; ii++){
1468 totcols += lblocks[ii]->ncols;
1469 assert ( lblocks[ii]->nrows == row1);
1470 }
1471
1472 struct Qmarray * temp =
1473 qmarray_alloc(row1*nlblocks*right->nrows, r);
1474 qmarray_block_kron_mat('D',0,nlblocks,lblocks,right,r,mat,temp);
1475
1476 size_t jj,kk,ll;
1477 size_t ind1,ind2;
1478 size_t runrows = 0;
1479 for (ii = 0; ii < nlblocks; ii++){
1480 //printf(" on block %zu\n",ii);
1481 if (ii == 0){
1482 for (jj = 0; jj < r; jj++){
1483 for (kk = 0; kk < row1; kk++){
1484 for (ll = 0; ll < right->nrows; ll++){
1485 ind1 = ll + kk * right->nrows +
1486 jj * row1 * right->nrows;
1487 ind2 = ll + kk * right->nrows +
1488 jj * row1*nlblocks * right->nrows;
1489 out->funcs[ind1] =
1490 generic_function_copy(temp->funcs[ind2]);
1491 }
1492 }
1493 }
1494 }
1495 else{
1496 for (jj = 0; jj < r; jj++){
1497 for (kk = 0; kk < row1; kk++){
1498 for (ll = 0; ll < right->nrows; ll++){
1499 ind1 = ll + kk * right->nrows +
1500 jj * row1 * right->nrows;
1501 ind2 = ll + (kk+runrows) * right->nrows +
1502 jj * row1 *nlblocks* right->nrows;
1503 generic_function_axpy(
1504 1.0,
1505 temp->funcs[ind2],
1506 out->funcs[ind1]);
1507 }
1508 }
1509 }
1510 }
1511 runrows += row1;
1512 }
1513 qmarray_free(temp); temp = NULL;
1514 }
1515 else if (type == 'C'){
1516 size_t totrows = lblocks[0]->nrows;
1517 size_t col1 = lblocks[0]->ncols;
1518 for (ii = 1; ii < nlblocks; ii++){
1519 totrows += lblocks[ii]->nrows;
1520 assert ( lblocks[ii]->ncols == col1);
1521 }
1522
1523 struct GenericFunction ** gfarray =
1524 generic_function_array_alloc(r * right->nrows * col1);
1525
1526 generic_function_kronh(0, r, col1, right->nrows, right->ncols,
1527 mat, right->funcs, gfarray);
1528
1529 size_t jj,kk,ll,ind1,ind2;
1530 size_t runrows = 0;
1531 for (ii = 0; ii < nlblocks; ii++){
1532 size_t st2a = right->nrows;
1533 size_t st2b = lblocks[ii]->nrows;
1534 size_t st2 = st2a * st2b * r;
1535
1536 struct GenericFunction ** temp2 =
1537 generic_function_array_alloc(st2);
1538
1539 generic_function_kronh2(0,r,st2a, st2b, lblocks[ii]->ncols,
1540 lblocks[ii]->funcs, gfarray, temp2);
1541
1542 for (jj = 0; jj < r; jj++){
1543 for (kk = 0; kk < st2b; kk++){
1544 for (ll = 0; ll < right->nrows; ll++){
1545 ind2 = ll + kk * right->nrows +
1546 jj * st2b * right->nrows;
1547 ind1 = ll + (kk + runrows)*right->nrows +
1548 jj* totrows * right->nrows;
1549 out->funcs[ind1] = temp2[ind2];
1550 }
1551 }
1552 }
1553 runrows += lblocks[ii]->nrows;
1554 free(temp2); temp2 = NULL;
1555 }
1556
1557 generic_function_array_free(gfarray,
1558 r * right->nrows * col1);
1559
1560 }
1561 }
1562 }
1563
1564
1565 /***********************************************************//**
1566 Integrate all the elements of a qmarray
1567
1568 \param[in] a - qmarray to integrate
1569
1570 \return out - array containing integral of every function in a
1571 ***************************************************************/
qmarray_integrate(const struct Qmarray * a)1572 double * qmarray_integrate(const struct Qmarray * a)
1573 {
1574
1575 double * out = calloc_double(a->nrows*a->ncols);
1576 size_t ii, jj;
1577 for (ii = 0; ii < a->ncols; ii++){
1578 for (jj = 0; jj < a->nrows; jj++){
1579 out[ii*a->nrows + jj] =
1580 generic_function_integrate(a->funcs[ii*a->nrows+jj]);
1581 }
1582 }
1583
1584 return out;
1585 }
1586
1587 /***********************************************************//**
1588 Weighted integration of all the elements of a qmarray
1589
1590 \param[in] a - qmarray to integrate
1591
1592 \return out - array containing integral of every function in a
1593
1594 \note Computes \f$ \int f(x) w(x) dx \f$ for every univariate function
1595 in the qmarray
1596
1597 w(x) depends on underlying parameterization
1598 for example, it is 1/2 for legendre (and default for others),
1599 gauss for hermite,etc
1600 ***************************************************************/
qmarray_integrate_weighted(const struct Qmarray * a)1601 double * qmarray_integrate_weighted(const struct Qmarray * a)
1602 {
1603
1604 double * out = calloc_double(a->nrows*a->ncols);
1605 size_t ii, jj;
1606 for (ii = 0; ii < a->ncols; ii++){
1607 for (jj = 0; jj < a->nrows; jj++){
1608 out[ii*a->nrows + jj] =
1609 generic_function_integrate_weighted(a->funcs[ii*a->nrows+jj]);
1610 }
1611 }
1612
1613 return out;
1614 }
1615
1616 /***********************************************************//**
1617 Norm of a qmarray
1618
1619 \param[in] a - first qmarray
1620
1621 \return L2 norm
1622 ***************************************************************/
qmarray_norm2(const struct Qmarray * a)1623 double qmarray_norm2(const struct Qmarray * a)
1624 {
1625 double out = 0.0;
1626 size_t ii;
1627 for (ii = 0; ii < a->ncols * a->nrows; ii++){
1628 // print_generic_function(a->funcs[ii],0,NULL);
1629 out += generic_function_inner(a->funcs[ii],a->funcs[ii]);
1630 // printf("out in qmarray norm = %G\n",out);
1631 }
1632 if (out < 0){
1633 fprintf(stderr,"Inner product between two qmarrays cannot be negative %G\n",out);
1634 exit(1);
1635 //return 0.0;
1636 }
1637
1638 return sqrt(fabs(out));
1639 }
1640
1641 /***********************************************************//**
1642 Two norm difference between two functions
1643
1644 \param[in] a - first qmarray
1645 \param[in] b - second qmarray
1646
1647 \return out - difference in 2 norm
1648 ***************************************************************/
qmarray_norm2diff(const struct Qmarray * a,const struct Qmarray * b)1649 double qmarray_norm2diff(const struct Qmarray * a,const struct Qmarray * b)
1650 {
1651 struct GenericFunction ** temp = generic_function_array_daxpby(
1652 a->ncols*a->nrows, 1.0,1,a->funcs,-1.0,1,b->funcs);
1653 double out = 0.0;
1654 size_t ii;
1655 for (ii = 0; ii < a->ncols * a->nrows; ii++){
1656 out += generic_function_inner(temp[ii],temp[ii]);
1657 }
1658 if (out < 0){
1659 fprintf(stderr,"Inner product between two qmarrays cannot be negative %G\n",out);
1660 exit(1);
1661 //return 0.0;
1662 }
1663 generic_function_array_free(temp, a->ncols*a->nrows);
1664 //free(temp);
1665 temp = NULL;
1666 return sqrt(out);
1667 }
1668
1669 /***********************************************************//**
1670 Compute \f[ y \leftarrow ax + y \f]
1671
1672 \param[in] a - scale for first qmarray
1673 \param[in] x - first qmarray
1674 \param[in,out] y - second qmarray
1675 ***************************************************************/
qmarray_axpy(double a,const struct Qmarray * x,struct Qmarray * y)1676 void qmarray_axpy(double a, const struct Qmarray * x, struct Qmarray * y)
1677 {
1678 assert (x->nrows == y->nrows );
1679 assert (x->ncols == y->ncols );
1680
1681 size_t ii;
1682 for (ii = 0; ii < x->nrows * x->ncols; ii++){
1683 generic_function_axpy(a,x->funcs[ii],y->funcs[ii]);
1684 }
1685 }
1686
1687 struct Zeros {
1688 double * zeros;
1689 size_t N;
1690 };
1691
eval_zpoly(size_t n,const double * x,double * out,void * args)1692 int eval_zpoly(size_t n, const double * x,double * out, void * args){
1693
1694 struct Zeros * z = args;
1695 for (size_t jj = 0; jj < n; jj++){
1696 out[jj] = 1.0;
1697 size_t ii;
1698 for (ii = 0; ii < z->N; ii++){
1699 out[jj] *= (x[jj] - z->zeros[ii]);
1700 }
1701 }
1702 return 0;
1703 }
1704
create_any_L(struct GenericFunction ** L,size_t nrows,size_t upto,size_t * piv,double * px,struct OneApproxOpts * approxopts,void * optargs)1705 void create_any_L(struct GenericFunction ** L, size_t nrows,
1706 size_t upto,size_t * piv, double * px,
1707 struct OneApproxOpts * approxopts, void * optargs)
1708 /* double lb, double ub, void * optargs) */
1709 {
1710
1711 //create an arbitrary quasimatrix array that has zeros at piv[:upto-1],px[:upto-1]
1712 // and one at piv[upto],piv[upto] less than one every else
1713 size_t ii,jj;
1714 size_t * count = calloc_size_t(nrows);
1715 double ** zeros = malloc( nrows * sizeof(double *));
1716 assert (zeros != NULL);
1717 for (ii = 0; ii < nrows; ii++){
1718 zeros[ii] = calloc_double(upto);
1719 for (jj = 0; jj < upto; jj++){
1720 if (piv[jj] == ii){
1721 zeros[ii][count[ii]] = px[jj];
1722 count[ii]++;
1723 }
1724 }
1725 if (count[ii] == 0){
1726 L[ii] = generic_function_constant(1.0,approxopts->fc,approxopts->aopts);//&ptype,lb,ub,NULL);
1727 }
1728 else{
1729 struct Zeros zz;
1730 zz.zeros = zeros[ii];
1731 zz.N = count[ii];
1732 struct Fwrap * fw = fwrap_create(1,"general-vec");
1733 fwrap_set_fvec(fw,eval_zpoly,&zz);
1734 L[ii] = generic_function_approximate1d(approxopts->fc,approxopts->aopts,fw);
1735 fwrap_destroy(fw);
1736 }
1737 free(zeros[ii]); zeros[ii] = NULL;
1738 }
1739
1740 free(zeros); zeros = NULL;
1741 free(count); count = NULL;
1742
1743 //this line is critical!!
1744 size_t amind;
1745 double xval;
1746 if (VQMALU){
1747 printf("inside of creatin g any L \n");
1748 printf("nrows = %zu \n",nrows);
1749
1750 //for (ii = 0; ii < nrows; ii++){
1751 // print_generic_function(L[ii],3,NULL);
1752 //}
1753 }
1754
1755 //double tval =
1756 generic_function_array_absmax(nrows, 1, L,&amind, &xval,optargs);
1757 px[upto] = xval;
1758 piv[upto] = amind;
1759 double val = generic_function_1d_eval(L[piv[upto]],px[upto]);
1760 generic_function_array_scale(1.0/val,L,nrows);
1761 if (VQMALU){
1762 printf("got new val = %G\n", val);
1763 printf("Values of new array at indices\n ");
1764 for (size_t zz = 0; zz < upto+1; zz++){
1765 double eval = generic_function_1d_eval(L[piv[zz]],px[zz]);
1766 printf("\t ind=%zu x=%G val=%G\n",piv[zz],px[zz],eval);
1767 }
1768 }
1769
1770 assert (fabs(val) > ZEROTHRESH);
1771
1772 }
1773
1774 /* static void create_any_L_nodal(struct GenericFunction ** L, */
1775 /* enum function_class fc, */
1776 /* size_t nrows, */
1777 /* size_t upto,size_t * piv, double * px, */
1778 /* struct OneApproxOpts * approxargs, void * optargs) */
1779 /* { */
1780
1781 /* //create an arbitrary quasimatrix array that has zeros at piv[:upto-1],px[:upto-1] */
1782 /* // and one at piv[upto],piv[upto] less than one every else */
1783 /* // printf("creating any L!\n"); */
1784
1785 /* assert (1 == 0); */
1786 /* double lb,ub; */
1787
1788 /* struct c3Vector * optnodes = NULL; */
1789 /* piv[upto] = 0; */
1790 /* if (optargs != NULL){ */
1791 /* /\* printf("Warning: optargs is not null in create_any_L_linelm so\n"); *\/ */
1792 /* /\* printf(" might not be able to guarantee that the new px\n"); *\/ */
1793 /* /\* printf(" is going to lie on a desired node\n"); *\/ */
1794 /* optnodes = optargs; */
1795 /* assert (optnodes->size > 3); */
1796 /* px[upto] = optnodes->elem[1]; */
1797 /* } */
1798 /* else{ */
1799 /* assert (1 == 0); */
1800 /* /\* px[upto] = lb + ub/(ub-lb) * randu(); *\/ */
1801 /* } */
1802
1803
1804 /* double * zeros = calloc_double(upto); */
1805 /* size_t iz = 0; */
1806 /* for (size_t ii = 0; ii < upto; ii++){ */
1807 /* if (piv[ii] == 0){ */
1808 /* if ( (fabs(px[ii] - lb) > 1e-15) && (fabs(px[ii]-ub) > 1e-15)){ */
1809 /* zeros[iz] = px[ii]; */
1810 /* iz++; */
1811 /* } */
1812 /* } */
1813 /* } */
1814
1815 /* int exists = 0; */
1816 /* size_t count = 2; */
1817 /* do { */
1818 /* exists = 0; */
1819 /* if (optargs != NULL){ */
1820 /* assert (count < optnodes->size); */
1821 /* } */
1822 /* for (size_t ii = 0; ii < iz; ii++){ */
1823 /* if (fabs(px[ii] - px[upto]) < 1e-15){ */
1824 /* exists = 1; */
1825 /* if (optargs != NULL){ */
1826 /* px[upto] = optnodes->elem[count]; */
1827 /* } */
1828 /* else{ */
1829 /* px[upto] = lb + ub/(ub-lb) * randu(); */
1830 /* } */
1831 /* count++; */
1832 /* break; */
1833 /* } */
1834 /* } */
1835 /* }while (exists == 1); */
1836
1837 /* /\* printf("before\n"); *\/ */
1838 /* /\* dprint(upto,px); *\/ */
1839 /* /\* iprint_sz(upto,piv); *\/ */
1840 /* /\* printf("%G\n",px[upto]); *\/ */
1841 /* /\* dprint(iz,zeros); *\/ */
1842 /* /\* L[0] = generic_function_onezero(LINELM,px[upto],iz,zeros,lb,ub); *\/ */
1843 /* L[0] = generic_function_onezero(fc,px[upto],iz,zeros,lb,ub); */
1844
1845 /* for (size_t ii = 1; ii < nrows;ii++){ */
1846 /* L[ii] = generic_function_constant(0.0,fc,approxargs); */
1847 /* } */
1848
1849 /* // for (size_t ii = 0; ii < upto; ii++){ */
1850 /* // double eval = generic_function_1d_eval( */
1851 /* // L[piv[ii]],px[ii]); */
1852 /* // printf("eval should be 0 = %G\n",eval); */
1853 /* // } */
1854 /* // double eval = generic_function_1d_eval(L[piv[upto]],px[upto]); */
1855 /* // printf("eval should be 1 = %G\n", eval); */
1856 /* // printf("after\n"); */
1857
1858 /* free(zeros); zeros = NULL; */
1859 /* } */
1860
1861 /***********************************************************//**
1862 Compute the LU decomposition of a quasimatrix array of 1d functioins
1863
1864 \param[in] A - qmarray to decompose
1865 \param[in,out] L - qmarray representing L factor
1866 \param[in,out] u - allocated space for U factor
1867 \param[in,out] piv - row of pivots
1868 \param[in,out] px - x values of pivots
1869 \param[in] app - approximation arguments
1870 \param[in] optargs - optimization arguments
1871
1872 \return info = 0 full rank <0 low rank ( rank = A->n + info )
1873 ***************************************************************/
qmarray_lu1d(struct Qmarray * A,struct Qmarray * L,double * u,size_t * piv,double * px,struct OneApproxOpts * app,void * optargs)1874 int qmarray_lu1d(struct Qmarray * A, struct Qmarray * L, double * u,
1875 size_t * piv, double * px, struct OneApproxOpts * app,
1876 void * optargs)
1877 {
1878 int info = 0;
1879
1880 size_t ii,kk;
1881 double val, amloc;
1882 size_t amind;
1883 struct GenericFunction ** temp = NULL;
1884
1885 for (kk = 0; kk < A->ncols; kk++){
1886
1887 if (VQMALU){
1888 printf("\n\n\n\n\n");
1889 printf("lu kk=%zu out of %zu \n",kk,A->ncols-1);
1890 printf("norm=%G\n",generic_function_norm(A->funcs[kk*A->nrows]));
1891 }
1892
1893 val = generic_function_array_absmax(A->nrows, 1, A->funcs + kk*A->nrows,
1894 &amind, &amloc,optargs);
1895 //printf("after absmax\n");
1896 piv[kk] = amind;
1897 px[kk] = amloc;
1898 //printf("absmax\n");
1899
1900 val = generic_function_1d_eval(A->funcs[kk*A->nrows+amind], amloc);
1901 if (VQMALU){
1902 printf("locmax=%3.15G\n",amloc);
1903 printf("amindmax=%zu\n",amind);
1904 printf("val of max =%3.15G\n",val);
1905 struct GenericFunction ** At = A->funcs+kk*L->nrows;
1906 double eval2 = generic_function_1d_eval(At[amind],amloc);
1907 printf("eval of thing %G\n",eval2);
1908 //print_generic_function(A->funcs[kk*A->nrows+amind],0,NULL);
1909 }
1910 if (fabs(val) <= ZEROTHRESH) {
1911 if (VQMALU){
1912 printf("creating any L\n");
1913 }
1914
1915 if ((A->funcs[kk*A->nrows]->fc == LINELM) ||
1916 (A->funcs[kk*A->nrows]->fc == CONSTELM)){
1917 //printf("lets go!\n");
1918 /* create_any_L_linelm(L->funcs+kk*L->nrows, */
1919 /* L->nrows,kk,piv,px,app,optargs); */
1920 generic_function_array_onezero(
1921 L->funcs+kk*L->nrows,
1922 L->nrows,
1923 app->fc,
1924 kk,piv,px,app->aopts);
1925 }
1926 else{
1927 create_any_L(L->funcs+kk*L->nrows,L->nrows,
1928 kk,piv,px,app,optargs);
1929 }
1930 amind = piv[kk];
1931 amloc = px[kk];
1932
1933 if (VQMALU){
1934 printf("done creating any L\n");
1935 }
1936 //print_generic_function(L->funcs[kk*L->nrows],0,NULL);
1937 //print_generic_function(L->funcs[kk*L->nrows+1],0,NULL);
1938 //printf("in here\n");
1939 //val = 0.0;
1940 }
1941 else{
1942 generic_function_array_daxpby2(A->nrows,1.0/val, 1,
1943 A->funcs + kk*A->nrows,
1944 0.0, 1, NULL, 1,
1945 L->funcs + kk * L->nrows);
1946 }
1947
1948 if (VQMALU){
1949 // printf("got new val = %G\n", val);
1950 printf("Values of new array at indices\n ");
1951 struct GenericFunction ** Lt = L->funcs+kk*L->nrows;
1952
1953 for (size_t zz = 0; zz < kk; zz++){
1954 double eval = generic_function_1d_eval(Lt[piv[zz]],px[zz]);
1955 printf("\t ind=%zu x=%G val=%G\n",piv[zz],px[zz],eval);
1956 }
1957 printf("-----------------------------------\n");
1958 }
1959
1960 //printf("k start here\n");
1961 for (ii = 0; ii < A->ncols; ii++){
1962 if (VQMALU){
1963 printf("In lu qmarray ii=%zu/%zu \n",ii,A->ncols);
1964 }
1965 if (ii == kk){
1966 u[ii*A->ncols+kk] = val;
1967 }
1968 else{
1969 u[ii*A->ncols+kk] = generic_function_1d_eval(
1970 A->funcs[ii*A->nrows + amind], amloc);
1971 }
1972 if (VQMALU){
1973 printf("u=%3.15G\n",u[ii*A->ncols+kk]);
1974 }
1975 /*
1976 print_generic_function(A->funcs[ii*A->nrows],3,NULL);
1977 print_generic_function(L->funcs[kk*L->nrows],3,NULL);
1978 */
1979 //printf("compute temp\n");
1980 temp = generic_function_array_daxpby(A->nrows, 1.0, 1,
1981 A->funcs + ii*A->nrows,
1982 -u[ii*A->ncols+kk], 1, L->funcs+ kk*L->nrows);
1983 if (VQMALU){
1984 //print_generic_function(A->funcs[ii*A->nrows],0,NULL);
1985 //printf("norm pre=%G\n",generic_function_norm(A->funcs[ii*A->nrows]));
1986 }
1987 //printf("got this daxpby\n");
1988 qmarray_set_column_gf(A,ii,temp);
1989 if (VQMALU){
1990 //printf("norm post=%G\n",generic_function_norm(A->funcs[ii*A->nrows]));
1991 }
1992 generic_function_array_free(temp,A->nrows);
1993 }
1994 }
1995 //printf("done?\n");
1996 // check if matrix is full rank
1997 for (kk = 0; kk < A->ncols; kk++){
1998 if (fabs(u[kk*A->ncols + kk]) < ZEROTHRESH){
1999 info --;
2000 }
2001 }
2002
2003 return info;
2004 }
2005
2006
remove_duplicates(size_t dim,size_t ** pivi,double ** pivx,double lb,double ub)2007 void remove_duplicates(size_t dim, size_t ** pivi, double ** pivx, double lb, double ub)
2008 {
2009
2010 size_t ii,jj;
2011 size_t * piv = *pivi;
2012 double * xiv = *pivx;
2013
2014 //printf("startindices\n");
2015 //dprint(dim,xiv);
2016 int done = 0;
2017 while (done == 0){
2018 done = 1;
2019 //printf("start again\n");
2020 for (ii = 0; ii < dim; ii++){
2021 for (jj = ii+1; jj < dim; jj++){
2022 if (piv[ii] == piv[jj]){
2023 double diff = fabs(xiv[ii] - xiv[jj]);
2024 double difflb = fabs(xiv[jj] - lb);
2025 //double diffub = fabs(xiv[jj] - ub);
2026
2027 if (diff < ZEROTHRESH){
2028 //printf("difflb=%G\n",difflb);
2029 if (difflb > ZEROTHRESH){
2030 xiv[jj] = (xiv[jj] + lb)/2.0;
2031 }
2032 else{
2033 // printf("use upper bound=%G\n",ub);
2034 xiv[jj] = (xiv[jj] + ub)/2.0;
2035 }
2036 done = 0;
2037 //printf("\n ii=%zu, old xiv[%zu]=%G\n",ii,jj,xiv[jj]);
2038 //xiv[jj] = 0.12345;
2039 //printf("new xiv[%zu]=%G\n",jj,xiv[jj]);
2040 break;
2041 }
2042 }
2043 }
2044 if (done == 0){
2045 break;
2046 }
2047 }
2048 //printf("indices\n");
2049 //dprint(dim,xiv);
2050 //done = 1;
2051 }
2052 }
2053
2054 /***********************************************************//**
2055 Compute the LU decomposition of a quasimatrix array of 1d functioins
2056
2057 \param[in] A - qmarray to decompose
2058 \param[in,out] L - qmarray representing L factor
2059 \param[in,out] u - allocated space for U factor
2060 \param[in,out] ps - pivot set
2061 \param[in] app - approximation arguments
2062 \param[in] optargs - optimization arguments
2063
2064 \return info = 0 full rank <0 low rank ( rank = A->n + info )
2065 \note THIS FUNCTION IS NOT READY YET
2066 ***************************************************************/
qmarray_lu1d_piv(struct Qmarray * A,struct Qmarray * L,double * u,struct PivotSet * ps,struct OneApproxOpts * app,void * optargs)2067 int qmarray_lu1d_piv(struct Qmarray * A, struct Qmarray * L, double * u,
2068 struct PivotSet * ps, struct OneApproxOpts * app,
2069 void * optargs)
2070 {
2071 assert (app != NULL);
2072 assert (1 == 0);
2073 int info = 0;
2074
2075 size_t ii,kk;
2076 double val;
2077 struct GenericFunction ** temp = NULL;
2078
2079 for (kk = 0; kk < A->ncols; kk++){
2080
2081 struct Pivot * piv = pivot_set_get_pivot(ps,kk);
2082 if (VQMALU){
2083 printf("\n\n\n\n\n");
2084 printf("lu kk=%zu out of %zu \n",kk,A->ncols-1);
2085 printf("norm=%G\n",generic_function_norm(A->funcs[kk*A->nrows]));
2086 }
2087
2088 val = generic_function_array_absmax_piv(A->nrows,1,A->funcs + kk*A->nrows,
2089 piv,optargs);
2090
2091 val = generic_function_1darray_eval_piv(A->funcs + kk*A->nrows,piv);
2092 if (VQMALU){
2093 printf("val of max =%3.15G\n",val);
2094 printf("eval of thing %G\n",val);
2095 }
2096 if (fabs(val) <= ZEROTHRESH) {
2097 if (VQMALU){
2098 printf("creating any L\n");
2099 }
2100
2101 if (A->funcs[kk*A->nrows]->fc == LINELM){
2102 assert (1 == 0);
2103 /* generic_function_array_onezero( */
2104 /* L->funcs+kk*L->nrows, */
2105 /* L->nrows, */
2106 /* app->fc, */
2107 /* kk,piv,px,app->aopts); */
2108 }
2109 else{
2110 assert (1 == 0);
2111 /* create_any_L(L->funcs+kk*L->nrows,L->nrows, */
2112 /* kk,piv,px,app,optargs); */
2113 }
2114
2115 /* if (VQMALU){ */
2116 /* printf("done creating any L\n"); */
2117 /* } */
2118 }
2119 else{
2120 generic_function_array_daxpby2(A->nrows,1.0/val, 1,
2121 A->funcs + kk*A->nrows,
2122 0.0, 1, NULL, 1,
2123 L->funcs + kk * L->nrows);
2124 }
2125
2126 if (VQMALU){
2127 // printf("got new val = %G\n", val);
2128 printf("Values of new array at indices\n ");
2129 struct GenericFunction ** Lt = L->funcs+kk*L->nrows;
2130
2131 for (size_t zz = 0; zz < kk; zz++){
2132 struct Pivot * pp = pivot_set_get_pivot(ps,zz);
2133 double eval = generic_function_1darray_eval_piv(Lt,pp);
2134 size_t ind = pivot_get_ind(pp);
2135 //printf("\t ind=%zu x=%G val=%G\n",piv[zz],px[zz],eval);
2136 printf("\t ind=%zu val=%G\n",ind,eval);
2137 }
2138 printf("-----------------------------------\n");
2139 }
2140
2141 //printf("k start here\n");
2142 for (ii = 0; ii < A->ncols; ii++){
2143 if (VQMALU){
2144 printf("In lu qmarray ii=%zu/%zu \n",ii,A->ncols);
2145 }
2146 if (ii == kk){
2147 u[ii*A->ncols+kk] = val;
2148 }
2149 else{
2150 u[ii*A->ncols+kk] =
2151 generic_function_1darray_eval_piv(A->funcs + ii*A->nrows,
2152 piv);
2153 }
2154 if (VQMALU){
2155 printf("u=%3.15G\n",u[ii*A->ncols+kk]);
2156 }
2157 /*
2158 print_generic_function(A->funcs[ii*A->nrows],3,NULL);
2159 print_generic_function(L->funcs[kk*L->nrows],3,NULL);
2160 */
2161 //printf("compute temp\n");
2162 temp = generic_function_array_daxpby(A->nrows, 1.0, 1,
2163 A->funcs + ii*A->nrows,
2164 -u[ii*A->ncols+kk], 1, L->funcs+ kk*L->nrows);
2165 if (VQMALU){
2166 //print_generic_function(A->funcs[ii*A->nrows],0,NULL);
2167 //printf("norm pre=%G\n",generic_function_norm(A->funcs[ii*A->nrows]));
2168 }
2169 //printf("got this daxpby\n");
2170 qmarray_set_column_gf(A,ii,temp);
2171 if (VQMALU){
2172 //printf("norm post=%G\n",generic_function_norm(A->funcs[ii*A->nrows]));
2173 }
2174 generic_function_array_free(temp,A->nrows);
2175 }
2176 }
2177 //printf("done?\n");
2178 // check if matrix is full rank
2179 for (kk = 0; kk < A->ncols; kk++){
2180 if (fabs(u[kk*A->ncols + kk]) < ZEROTHRESH){
2181 info --;
2182 }
2183 }
2184
2185 return info;
2186 }
2187
2188
2189 /***********************************************************//**
2190 Perform a greedy maximum volume procedure to find the
2191 maximum volume submatrix of a qmarray
2192
2193 \param[in] A - qmarray
2194 \param[in,out] Asinv - submatrix inv
2195 \param[in,out] pivi - row pivots
2196 \param[in,out] pivx - x value in row pivots
2197 \param[in] appargs - approximation arguments
2198 \param[in] optargs - optimization arguments
2199
2200 \return info =
2201 0 converges,
2202 < 0 no invertible submatrix,
2203 >0 rank may be too high,
2204
2205 \note
2206 naive implementation without rank 1 updates
2207 ***************************************************************/
qmarray_maxvol1d(struct Qmarray * A,double * Asinv,size_t * pivi,double * pivx,struct OneApproxOpts * appargs,void * optargs)2208 int qmarray_maxvol1d(struct Qmarray * A, double * Asinv, size_t * pivi,
2209 double * pivx, struct OneApproxOpts * appargs, void * optargs)
2210 {
2211 //printf("in maxvolqmarray!\n");
2212
2213 int info = 0;
2214
2215 size_t r = A->ncols;
2216
2217 struct Qmarray * L = qmarray_alloc(A->nrows, r);
2218 double * U = calloc_double(r * r);
2219
2220 struct Qmarray * Acopy = qmarray_copy(A);
2221
2222 if (VQMAMAXVOL){
2223 printf("==*=== \n\n\n\n In MAXVOL \n");
2224 /* size_t ll; */
2225 /* for (ll = 0; ll < A->nrows * A->ncols; ll++){ */
2226 /* printf("%G\n", generic_function_norm(Acopy->funcs[ll])); */
2227 /* } */
2228 }
2229
2230 //print_qmarray(Acopy,0,NULL);
2231 info = qmarray_lu1d(Acopy,L,U,pivi,pivx,appargs,optargs);
2232 if (VQMAMAXVOL){
2233 printf("pivot immediate \n");
2234 iprint_sz(A->ncols, pivi);
2235 //pivx[A->ncols-1] = 0.12345;
2236 dprint(A->ncols,pivx);
2237 }
2238 if (info > 0){
2239 //printf("Couldn't find an invertible submatrix for maxvol %d\n",info);
2240 printf("Error in input %d of qmarray_lu1d\n",info);
2241 //printf("Couldn't Error from quasimatrix_lu1d \n");
2242 qmarray_free(L);
2243 qmarray_free(Acopy);
2244 free(U);
2245 return info;
2246 }
2247
2248 size_t ii,jj;
2249 for (ii = 0; ii < r; ii++){
2250 for (jj = 0; jj < r; jj++){
2251 U[jj * r + ii] = generic_function_1d_eval(
2252 A->funcs[jj*A->nrows + pivi[ii]], pivx[ii]);
2253 }
2254 }
2255 /* if (VQMAMAXVOL){ */
2256 /* printf("built U\nn"); */
2257 /* } */
2258
2259 int * ipiv2 = calloc(r, sizeof(int));
2260 size_t lwork = r*r;
2261 double * work = calloc(lwork, sizeof(double));
2262 //int info2;
2263
2264 pinv(r,r,r,U,Asinv,ZEROTHRESH);
2265 /*
2266 dgetrf_(&r,&r, Asinv, &r, ipiv2, &info2);
2267 printf("info2a=%d\n",info2);
2268 dgetri_(&r, Asinv, &r, ipiv2, work, &lwork, &info2); //invert
2269 printf("info2b=%d\n",info2);
2270 */
2271
2272 /* if (VQMAMAXVOL){ */
2273 /* printf("took pseudoinverse \nn"); */
2274 /* } */
2275 struct Qmarray * B = qmam(A,Asinv,r);
2276
2277 /* if (VQMAMAXVOL){ */
2278 /* printf("did qmam \n"); */
2279 /* } */
2280 size_t maxrow, maxcol;
2281 double maxx, maxval, maxval2;
2282
2283 //printf("B size = (%zu,%zu)\n",B->nrows, B->ncols);
2284 // BLAH 2
2285
2286 /* if (VQMAMAXVOL){ */
2287 /* printf("do absmax1d\n"); */
2288 /* } */
2289 qmarray_absmax1d(B, &maxx, &maxrow, &maxcol, &maxval,optargs);
2290 if (VQMAMAXVOL){
2291 printf("maxx=%G,maxrow=%zu maxcol=%zu maxval=%G\n",maxx,maxrow, maxcol, maxval);
2292 }
2293 size_t maxiter = 30;
2294 size_t iter =0;
2295 double delta = 1e-2;
2296 while (maxval > (1.0 + delta)){
2297 pivi[maxcol] = maxrow;
2298 pivx[maxcol] = maxx;
2299
2300 //printf(" update Asinv A size = (%zu,%zu) \n",A->nrows,A->ncols);
2301 for (ii = 0; ii < r; ii++){
2302 //printf("pivx[%zu]=%3.5f pivi=%zu \n",ii,pivx[ii],pivi[ii]);
2303 for (jj = 0; jj < r; jj++){
2304 //printf("jj=%zu\n",jj);
2305 //print_generic_function(A->funcs[jj*A->nrows+pivi[ii]],0,NULL);
2306 U[jj * r + ii] = generic_function_1d_eval(
2307 A->funcs[jj*A->nrows + pivi[ii]], pivx[ii]);
2308 //printf("got jj \n");
2309 }
2310 }
2311
2312 //printf(" before dgetrf \n");
2313
2314 pinv(r,r,r,U,Asinv,ZEROTHRESH);
2315 //dgetrf_(&r,&r, Asinv, &r, ipiv2, &info2);
2316 //dgetri_(&r, Asinv, &r, ipiv2, work, &lwork, &info2); //invert
2317 qmarray_free(B); B = NULL;
2318 B = qmam(A,Asinv,r);
2319 qmarray_absmax1d(B, &maxx, &maxrow, &maxcol, &maxval2,optargs);
2320 if (VQMAMAXVOL){
2321 printf("maxx=%G,maxrow=%zu maxcol=%zu maxval=%G\n",maxx,maxrow, maxcol, maxval2);
2322 }
2323
2324 /* if (fabs(maxval2-maxval)/fabs(maxval) < 1e-2){ */
2325 /* break; */
2326 /* } */
2327 if (iter > maxiter){
2328 if (VQMAMAXVOL){
2329 printf("maximum iterations (%zu) reached\n",maxiter);
2330 }
2331 break;
2332 }
2333 maxval = maxval2;
2334 //printf("maxval=%G\n",maxval);
2335 iter++;
2336 }
2337
2338 if (iter > 0){
2339 pivi[maxcol] = maxrow;
2340 pivx[maxcol] = maxx;
2341 }
2342
2343 if (VQMAMAXVOL){
2344 printf("qmarray_maxvol indices and pivots before removing duplicates\n");
2345 iprint_sz(r,pivi);
2346 dprint(r,pivx);
2347 }
2348 //printf("done\n");
2349 //if ( r < 0){
2350 if ( r > 1){
2351 double lb = generic_function_get_lb(A->funcs[0]);
2352 double ub = generic_function_get_ub(A->funcs[0]);
2353 remove_duplicates(r, &pivi, &pivx,lb,ub);
2354 for (ii = 0; ii < r; ii++){
2355 for (jj = 0; jj < r; jj++){
2356 U[jj * r + ii] = generic_function_1d_eval(
2357 A->funcs[jj*A->nrows + pivi[ii]], pivx[ii]);
2358 }
2359 }
2360 pinv(r,r,r,U,Asinv,ZEROTHRESH);
2361 }
2362
2363 if (VQMAMAXVOL){
2364 printf("qmarray_maxvol indices and pivots after removing duplicates\n");
2365 iprint_sz(r,pivi);
2366 dprint(r,pivx);
2367 }
2368 //printf("done with maxval = %3.4f\n", maxval);
2369
2370 free(work);
2371 free(ipiv2);
2372 qmarray_free(B);
2373 qmarray_free(L);
2374 qmarray_free(Acopy);
2375 free(U);
2376 return info;
2377 }
2378
2379 /***********************************************************//**
2380 Compute the householder triangularization of a quasimatrix array.
2381
2382 \param[in,out] A - qmarray to triangularize (overwritten)
2383 \param[in,out] E - qmarray with orthonormal columns
2384 \param[in,out] V - allocated space for a qmarray (will store reflectors)
2385 \param[in,out] R - allocated space upper triangular factor
2386
2387 \return info - (if != 0 then something is wrong)
2388 ***************************************************************/
2389 int
qmarray_householder(struct Qmarray * A,struct Qmarray * E,struct Qmarray * V,double * R)2390 qmarray_householder(struct Qmarray * A, struct Qmarray * E,
2391 struct Qmarray * V, double * R)
2392 {
2393 //printf("\n\n\n\n\n\n\n\n\nSTARTING QMARRAY_HOUSEHOLDER\n");
2394 size_t ii, jj;
2395 struct Quasimatrix * v = NULL;
2396 struct Quasimatrix * v2 = NULL;
2397 struct Quasimatrix * atemp = NULL;
2398 double rho, sigma, alpha;
2399 double temp1;
2400
2401 //printf("ORTHONORMAL QMARRAY IS \n");
2402 //print_qmarray(E,0,NULL);
2403 //printf("qm array\n");
2404 struct GenericFunction ** vfuncs = NULL;
2405 for (ii = 0; ii < A->ncols; ii++){
2406 if (VQMAHOUSEHOLDER){
2407 printf(" On iter (%zu / %zu) \n", ii, A->ncols);
2408 }
2409 //printf("function is\n");
2410 //print_generic_function(A->funcs[ii*A->nrows],3,NULL);
2411
2412 rho = generic_function_array_norm(A->nrows,1,A->funcs+ii*A->nrows);
2413
2414 if (rho < ZEROTHRESH){
2415 rho = 0.0;
2416 }
2417
2418 if (VQMAHOUSEHOLDER){
2419 printf(" \t ... rho = %3.15G ZEROTHRESH=%3.15G\n",
2420 rho,ZEROTHRESH);
2421 }
2422
2423 //printf(" new E = \n");
2424 //print_generic_function(E->funcs[ii*E->nrows],3,NULL);
2425 R[ii * A->ncols + ii] = rho;
2426 alpha = generic_function_inner_sum(A->nrows,1,E->funcs+ii*E->nrows,
2427 1, A->funcs+ii*A->nrows);
2428
2429 if (fabs(alpha) < ZEROTHRESH)
2430 alpha=0.0;
2431
2432 if (VQMAHOUSEHOLDER){
2433 printf(" \t ... alpha = %3.15G\n",alpha);
2434 }
2435
2436 if (alpha >= ZEROTHRESH){
2437 //printf("flip sign\n");
2438 generic_function_array_flip_sign(E->nrows,1,E->funcs+ii*E->nrows);
2439 }
2440 v = quasimatrix_alloc(A->nrows);
2441 quasimatrix_free_funcs(v);
2442
2443 /* free(v->funcs);v->funcs=NULL; */
2444 //printf("array daxpby nrows = %zu\n",A->nrows);
2445 //printf("epostalpha_qma = \n");
2446 //print_generic_function(E->funcs[ii*E->nrows],3,NULL);
2447 //printf("Apostalpha_qma = \n");
2448 //print_generic_function(A->funcs[ii*A->nrows],3,NULL);
2449
2450 vfuncs = generic_function_array_daxpby(A->nrows,rho,1,
2451 E->funcs+ii*E->nrows,-1.0,
2452 1, A->funcs+ii*A->nrows);
2453 quasimatrix_set_funcs(v,vfuncs);
2454
2455 /* printf("sup!\n"); */
2456 //printf("v update = \n");
2457 //print_generic_function(v->funcs[0],3,NULL);
2458
2459 //printf("QMARRAY MOD IS \n");
2460 //print_qmarray(E,0,NULL);
2461 // improve orthogonality
2462 //*
2463 if (ii > 1){
2464 struct Quasimatrix * tempv = NULL;
2465 for (jj = 0; jj < ii-1; jj++){
2466 //printf("jj=%zu (max is %zu) \n",jj, ii-1);
2467 tempv = quasimatrix_copy(v);
2468 struct GenericFunction ** tvfuncs = NULL;
2469 tvfuncs = quasimatrix_get_funcs(tempv);
2470 //printf("compute temp\n");
2471 //print_quasimatrix(tempv, 0, NULL);
2472 temp1 = generic_function_inner_sum(A->nrows,1,tvfuncs,1,
2473 E->funcs + jj*E->nrows);
2474 /* quasimatrix_set_funcs(tempv,tvfuncs); */
2475
2476 //printf("temp1= %G\n",temp1);
2477 /* generic_function_array_free(vfuncs,A->nrows); */
2478 quasimatrix_free_funcs(v);
2479 /* vfuncs = quasimatrix_get_funcs(v); */
2480 vfuncs = generic_function_array_daxpby(A->nrows,1.0,1,tvfuncs,
2481 -temp1,1,E->funcs+jj*E->nrows);
2482 quasimatrix_set_funcs(v,vfuncs);
2483 //printf("k ok= %G\n",temp1);
2484 quasimatrix_free(tempv);
2485 }
2486 }
2487 //*/
2488
2489 /* printf("compute sigma\n"); */
2490 //print_generic_function(v->funcs[0],3,NULL);
2491 /* sigma = generic_function_array_norm(A->nrows, 1, vfuncs); */
2492 sigma = quasimatrix_norm(v);
2493 /* printf("get it!\n"); */
2494 if (sigma < ZEROTHRESH){
2495 sigma = 0.0;
2496 }
2497
2498 if (VQMAHOUSEHOLDER){
2499 printf(" \t ... sigma = %3.15G ZEROTHRESH=%3.15G\n",
2500 sigma,ZEROTHRESH);
2501 }
2502 //printf(" \t ... sigma = %3.15G ZEROTHRESH=%3.15G\n",
2503 // sigma,ZEROTHRESH);
2504
2505 //printf("v2preqma = \n");
2506 //print_generic_function(v->funcs[0],3,NULL);
2507 double ttol = ZEROTHRESH;
2508 struct GenericFunction ** v2fs;
2509 if (fabs(sigma) <= ttol){
2510 //printf("here sigma too small!\n");
2511 v2 = quasimatrix_alloc(A->nrows);
2512 quasimatrix_free_funcs(v2);
2513 /* free(v2->funcs);v2->funcs=NULL; */
2514 //just a copy
2515 v2fs = generic_function_array_daxpby(A->nrows,1.0, 1,
2516 E->funcs+ii*E->nrows, 0.0, 1, NULL);
2517 quasimatrix_set_funcs(v2,v2fs);
2518 }
2519 else {
2520 //printf("quasi daxpby\n");
2521 //printf("sigma = %G\n",sigma);
2522 /* printf("we are here\n"); */
2523 v2 = quasimatrix_daxpby(1.0/sigma, v, 0.0, NULL);
2524 /* printf("and now there\n"); */
2525 v2fs = quasimatrix_get_funcs(v2);
2526 //printf("quasi got it daxpby\n");
2527 }
2528 quasimatrix_free(v);
2529 qmarray_set_column(V,ii,v2);
2530
2531 //printf("v2qma = \n");
2532 //print_generic_function(v2->funcs[0],3,NULL);
2533 //printf("go into cols \n");
2534 for (jj = ii+1; jj < A->ncols; jj++){
2535
2536 if (VQMAHOUSEHOLDER){
2537 printf("\t On sub-iter %zu\n", jj);
2538 }
2539 //printf("Aqma = \n");
2540 //print_generic_function(A->funcs[jj*A->nrows],3,NULL);
2541 temp1 = generic_function_inner_sum(A->nrows,1,v2fs,1,
2542 A->funcs+jj*A->nrows);
2543
2544 //printf("\t\t temp1 pre = %3.15G DBL_EPS=%G\n", temp1,DBL_EPSILON);
2545 if (fabs(temp1) < ZEROTHRESH){
2546 temp1 = 0.0;
2547 }
2548
2549 atemp = quasimatrix_alloc(A->nrows);
2550 quasimatrix_free_funcs(atemp);
2551 struct GenericFunction ** afuncs = quasimatrix_get_funcs(atemp);
2552 afuncs = generic_function_array_daxpby(A->nrows, 1.0, 1,
2553 A->funcs+jj*A->nrows,
2554 -2.0 * temp1, 1,
2555 v2fs);
2556 quasimatrix_set_funcs(atemp,afuncs);
2557
2558 //printf("atemp = \n");
2559 //print_generic_function(atemp->funcs[0], 0, NULL);
2560 R[jj * A->ncols + ii] = generic_function_inner_sum(E->nrows, 1,
2561 E->funcs + ii*E->nrows, 1, afuncs);
2562
2563
2564 //double aftemp =
2565 // generic_function_array_norm(A->nrows,1,atemp->funcs);
2566
2567 //double eftemp =
2568 // generic_function_array_norm(E->nrows,1,E->funcs+ii*E->nrows);
2569
2570 if (VQMAHOUSEHOLDER){
2571 //printf("\t \t ... temp1 = %3.15G\n",temp1);
2572 //printf("\t\t e^T atemp= %3.15G\n", R[jj*A->ncols+ii]);
2573 }
2574
2575 v = quasimatrix_alloc(A->nrows);
2576 quasimatrix_free_funcs(v);
2577 /* free(v->funcs);v->funcs=NULL; */
2578 vfuncs = generic_function_array_daxpby(A->nrows, 1.0, 1,
2579 afuncs,
2580 -R[jj * A->ncols + ii], 1,
2581 E->funcs + ii*E->nrows);
2582 quasimatrix_set_funcs(v,vfuncs);
2583
2584
2585 if (VQMAHOUSEHOLDER){
2586 //double temprho =
2587 // generic_function_array_norm(A->nrows,1,v->funcs);
2588 //printf("new v func norm = %3.15G\n",temprho);
2589 }
2590
2591 qmarray_set_column(A,jj,v); // overwrites column jj
2592
2593 quasimatrix_free(atemp); atemp=NULL;
2594 quasimatrix_free(v); v = NULL;
2595 }
2596
2597 quasimatrix_free(v2);
2598 }
2599 /* printf("done?\n"); */
2600 //printf("qmarray v after = \n");
2601 //print_qmarray(V,0,NULL);
2602
2603 //printf("qmarray E after = \n");
2604 //print_qmarray(E,0,NULL);
2605 /* printf("\n\n\n\n\n\n\n\n\n Ending QMARRAY_HOUSEHOLDER\n"); */
2606 return 0;
2607 }
2608
2609 /***********************************************************//**
2610 Compute the Q for the QR factor from householder reflectors
2611 for a Quasimatrix Array
2612
2613 \param[in,out] Q - qmarray E obtained afterm qmarray_householder
2614 becomes overwritten
2615 \param[in,out] V - Householder functions, obtained from qmarray_householder
2616
2617 \return info - if != 0 then something is wrong)
2618 ***************************************************************/
qmarray_qhouse(struct Qmarray * Q,struct Qmarray * V)2619 int qmarray_qhouse(struct Qmarray * Q, struct Qmarray * V)
2620 {
2621
2622 int info = 0;
2623 size_t ii, jj;
2624
2625 struct Quasimatrix * q2 = NULL;
2626 double temp1;
2627 size_t counter = 0;
2628 ii = Q->ncols-1;
2629
2630 /*
2631 printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
2632 printf("************Q beforer*************\n");
2633 print_qmarray(Q, 0, NULL);
2634 printf("************V beforer*************\n");
2635 print_qmarray(V, 0, NULL);
2636 */
2637 while (counter < Q->ncols){
2638 /* printf("INCREMENT COUNTER \n"); */
2639 for (jj = ii; jj < Q->ncols; jj++){
2640
2641 temp1 = generic_function_inner_sum(V->nrows,1,V->funcs + ii*V->nrows,
2642 1, Q->funcs + jj*Q->nrows);
2643 //printf("temp1=%G\n",temp1);
2644 q2 = quasimatrix_alloc(Q->nrows);
2645 quasimatrix_free_funcs(q2);
2646 struct GenericFunction ** q2funcs = quasimatrix_get_funcs(q2);
2647 /* printf("lets go\n"); */
2648 q2funcs = generic_function_array_daxpby(Q->nrows, 1.0, 1,
2649 Q->funcs + jj * Q->nrows, -2.0 * temp1, 1,
2650 V->funcs + ii * V->nrows);
2651 /* printf("why?!\n"); */
2652 quasimatrix_set_funcs(q2,q2funcs);
2653 qmarray_set_column(Q,jj, q2);
2654 //printf("Q new = \n");
2655 //print_generic_function(Q->funcs[jj],1,NULL);
2656 quasimatrix_free(q2); q2 = NULL;
2657 }
2658 /*
2659 printf(" *-========================== Q temp ============== \n");
2660 print_qmarray(Q, 0, NULL);
2661 */
2662 ii = ii - 1;
2663 counter = counter + 1;
2664 }
2665 /* printf("done\n"); */
2666 //printf("************Q after*************\n");
2667 //print_qmarray(Q, 0, NULL);
2668 //printf("\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n");
2669
2670 return info;
2671 }
2672
2673 /***********************************************************//**
2674 Compute the householder triangularization of a
2675 quasimatrix array. (LQ= A), Q orthonormal,
2676 L lower triangular
2677
2678 \param[in,out] A - qmarray to triangularize (destroyed
2679 \param[in,out] E - qmarray with orthonormal rows
2680 \param[in,out] V - allocated space for a qmarray (will store reflectors)
2681 \param[in,out] L - allocated space lower triangular factor
2682
2683 \return info (if != 0 then something is wrong)
2684 ***************************************************************/
2685 int
qmarray_householder_rows(struct Qmarray * A,struct Qmarray * E,struct Qmarray * V,double * L)2686 qmarray_householder_rows(struct Qmarray * A, struct Qmarray * E,
2687 struct Qmarray * V, double * L)
2688 {
2689 size_t ii, jj;
2690 struct Quasimatrix * e = NULL;
2691 struct Quasimatrix * x = NULL;
2692 struct Quasimatrix * v = NULL;//quasimatrix_alloc(A->ncols);
2693 struct Quasimatrix * v2 = NULL;
2694 struct Quasimatrix * atemp = NULL;
2695 double rho, sigma, alpha;
2696 double temp1;
2697
2698 struct GenericFunction ** vfuncs = NULL;
2699 for (ii = 0; ii < A->nrows; ii++){
2700 if (VQMAHOUSEHOLDER)
2701 printf(" On iter %zu\n", ii);
2702
2703 rho = generic_function_array_norm(A->ncols,A->nrows,A->funcs+ii);
2704
2705 if (rho < ZEROTHRESH){
2706 rho = 0.0;
2707 }
2708
2709 if (VQMAHOUSEHOLDER){
2710 printf(" \t ... rho = %3.15G ZEROTHRESH=%3.15G\n",
2711 rho,ZEROTHRESH);
2712 }
2713
2714 L[ii * A->nrows + ii] = rho;
2715 alpha = generic_function_inner_sum(A->ncols,E->nrows,E->funcs+ii,
2716 A->nrows, A->funcs+ii);
2717 if (fabs(alpha) < ZEROTHRESH){
2718 alpha = 0.0;
2719 }
2720
2721 if (VQMAHOUSEHOLDER)
2722 printf(" \t ... alpha = %3.15G\n",alpha);
2723
2724 if (alpha >= ZEROTHRESH){
2725 generic_function_array_flip_sign(E->ncols,E->nrows,E->funcs+ii);
2726 }
2727
2728 v = quasimatrix_alloc(A->ncols);
2729 quasimatrix_free_funcs(v);
2730
2731 vfuncs = generic_function_array_daxpby(A->ncols, rho, E->nrows,
2732 E->funcs+ii, -1.0, A->nrows, A->funcs+ii);
2733 quasimatrix_set_funcs(v,vfuncs);
2734
2735 // improve orthogonality
2736 //*
2737 if (ii > 1){
2738 struct Quasimatrix * tempv = NULL;
2739 for (jj = 0; jj < ii-1; jj++){
2740 tempv = quasimatrix_copy(v);
2741 struct GenericFunction ** tvfuncs = quasimatrix_get_funcs(tempv);
2742 temp1 = generic_function_inner_sum(A->ncols,1,tvfuncs,E->nrows,
2743 E->funcs + jj);
2744 quasimatrix_free_funcs(v);
2745 vfuncs = generic_function_array_daxpby(A->ncols,1.0,1,tvfuncs,
2746 -temp1,E->nrows,E->funcs+jj);
2747 quasimatrix_set_funcs(v,vfuncs);
2748 quasimatrix_free(tempv);
2749 }
2750 }
2751 //*/
2752 sigma = quasimatrix_norm(v);//generic_function_array_norm(v->n, 1, v->funcs);
2753 if (sigma < ZEROTHRESH){
2754 sigma = 0.0;
2755 }
2756 //
2757 double ttol = ZEROTHRESH;
2758
2759 if (VQMAHOUSEHOLDER){
2760 printf(" \t ... sigma = %G ttol=%G \n",sigma,ttol);
2761 }
2762 struct GenericFunction ** v2funcs;
2763
2764 if (fabs(sigma) <= ttol){
2765 //printf("HERE SIGMA TOO SMALL\n");
2766 v2 = quasimatrix_alloc(A->ncols);
2767 quasimatrix_free_funcs(v2);
2768 //just a copy
2769 v2funcs = generic_function_array_daxpby(E->ncols,1.0, E->nrows,
2770 E->funcs+ii, 0.0, 1, NULL);
2771 quasimatrix_set_funcs(v2,v2funcs);
2772 }
2773 else {
2774 v2 = quasimatrix_daxpby(1.0/sigma, v, 0.0, NULL);
2775 v2funcs = quasimatrix_get_funcs(v2);
2776 }
2777 quasimatrix_free(v);
2778 qmarray_set_row(V,ii,v2);
2779
2780 for (jj = ii+1; jj < A->nrows; jj++){
2781
2782 if (VQMAHOUSEHOLDER){
2783 //printf("\t On sub-iter %zu\n", jj);
2784 }
2785
2786 temp1 = generic_function_inner_sum(A->ncols,1,v2funcs,
2787 A->nrows, A->funcs+jj);
2788
2789 if (fabs(temp1) < ZEROTHRESH){
2790 temp1 = 0.0;
2791 }
2792 //if (VQMAHOUSEHOLDER)
2793 // printf("\t\t temp1= %3.15G\n", temp1);
2794
2795 atemp = quasimatrix_alloc(A->ncols);
2796 quasimatrix_free_funcs(atemp);
2797 struct GenericFunction ** afuncs = NULL;
2798 /* free(atemp->funcs); */
2799 afuncs = generic_function_array_daxpby(A->ncols, 1.0,
2800 A->nrows, A->funcs+jj,
2801 -2.0 * temp1, 1,
2802 v2funcs);
2803 quasimatrix_set_funcs(atemp,afuncs);
2804
2805 L[ii * A->nrows + jj] = generic_function_inner_sum(E->ncols,
2806 E->nrows,
2807 E->funcs + ii, 1,
2808 afuncs);
2809
2810 //if (VQMAHOUSEHOLDER)
2811 // printf("\t\t e^T atemp= %3.15G\n", L[ii*A->nrows+jj]);
2812
2813 v = quasimatrix_alloc(A->ncols);
2814 quasimatrix_free_funcs(v);
2815 /* free(v->funcs);v->funcs=NULL; */
2816 /* vfuncs = quasimatrix_get_funcs(v); */
2817 vfuncs = generic_function_array_daxpby(A->ncols, 1.0, 1,
2818 afuncs, -L[ii * A->nrows + jj],
2819 E->nrows, E->funcs + ii);
2820 quasimatrix_set_funcs(v,vfuncs);
2821
2822 qmarray_set_row(A,jj,v); // overwrites column jj
2823
2824 quasimatrix_free(atemp); atemp=NULL;
2825 quasimatrix_free(v); v = NULL;
2826 }
2827
2828 quasimatrix_free(x);
2829 quasimatrix_free(e);
2830 quasimatrix_free(v2);
2831 }
2832 return 0;
2833 }
2834
2835 /***********************************************************//**
2836 Compute the Q matrix from householder reflector for LQ decomposition
2837
2838 \param[in,out] Q - qmarray E obtained afterm qmarray_householder_rows (overwrite)
2839 \param[in,out] V - Householder functions, obtained from qmarray_householder
2840
2841 \return info - (if != 0 then something is wrong)
2842 ***************************************************************/
qmarray_qhouse_rows(struct Qmarray * Q,struct Qmarray * V)2843 int qmarray_qhouse_rows(struct Qmarray * Q, struct Qmarray * V)
2844 {
2845 int info = 0;
2846 size_t ii, jj;
2847
2848 struct Quasimatrix * q2 = NULL;
2849 double temp1;
2850 size_t counter = 0;
2851 ii = Q->nrows-1;
2852 while (counter < Q->nrows){
2853 for (jj = ii; jj < Q->nrows; jj++){
2854 temp1 = generic_function_inner_sum(V->ncols, V->nrows, V->funcs + ii,
2855 Q->nrows, Q->funcs + jj);
2856
2857 q2 = quasimatrix_alloc(Q->ncols);
2858 quasimatrix_free_funcs(q2);
2859 struct GenericFunction ** q2funcs = quasimatrix_get_funcs(q2);
2860 q2funcs = generic_function_array_daxpby(Q->ncols, 1.0, Q->nrows,
2861 Q->funcs + jj, -2.0 * temp1, V->nrows,
2862 V->funcs + ii);
2863 quasimatrix_set_funcs(q2,q2funcs);
2864
2865 qmarray_set_row(Q,jj, q2);
2866 quasimatrix_free(q2); q2 = NULL;
2867 }
2868 ii = ii - 1;
2869 counter = counter + 1;
2870 }
2871
2872 return info;
2873 }
2874
2875 /***********************************************************//**
2876 Compute the householder triangularization of a
2877 qmarray. whose elements consist of
2878 one dimensional functions (simple interface)
2879
2880 \param[in] dir - type either "QR" or "LQ"
2881 \param[in,out] A - qmarray to triangularize (destroyed in call)
2882 \param[in,out] R - allocated space upper triangular factor
2883 \param[in] app - approximation arguments
2884
2885 \return Q
2886 ***************************************************************/
2887 struct Qmarray *
qmarray_householder_simple(char * dir,struct Qmarray * A,double * R,struct OneApproxOpts * app)2888 qmarray_householder_simple(char * dir,struct Qmarray * A,double * R,
2889 struct OneApproxOpts * app)
2890 {
2891
2892 size_t ncols = A->ncols;
2893
2894 // generate two quasimatrices needed for householder
2895 struct Qmarray * Q = NULL;
2896 if (strcmp(dir,"QR") == 0){
2897
2898 Q = qmarray_orth1d_columns(A->nrows,ncols,app);
2899
2900
2901 if (app->fc != PIECEWISE){
2902 int out = qmarray_qr(A,&Q,&R,app);
2903 assert (out == 0);
2904 }
2905 else{
2906 struct Qmarray * V = qmarray_alloc(A->nrows,ncols);
2907 int out = qmarray_householder(A,Q,V,R);
2908 assert(out == 0);
2909 out = qmarray_qhouse(Q,V);
2910 assert(out == 0);
2911 qmarray_free(V);
2912 }
2913 }
2914 else if (strcmp(dir, "LQ") == 0){
2915 /* printf("lets go LQ!\n"); */
2916 Q = qmarray_orth1d_rows(A->nrows,ncols,app);
2917 /* printf("got orth1d_rows!\n"); */
2918 if (app->fc != PIECEWISE){
2919 //free(R); R = NULL;
2920 int out = qmarray_lq(A,&Q,&R,app);
2921 assert (out == 0);
2922 }
2923 else{
2924 struct Qmarray * V = qmarray_alloc(A->nrows,ncols);
2925 int out = qmarray_householder_rows(A,Q,V,R);
2926 /* printf("\n\n\n\n"); */
2927 /* printf("right after qmarray_householder_rows \n"); */
2928 /* print_qmarray(Q,0,NULL); */
2929
2930 /* printf("\n\n\n\n"); */
2931 /* printf("VV right after qmarray_householder_rows \n"); */
2932 /* print_qmarray(V,0,NULL); */
2933
2934 assert(out == 0);
2935 out = qmarray_qhouse_rows(Q,V);
2936
2937 /* printf("\n\n\n\n"); */
2938 /* printf("right after qmarray_qhouse_rows \n"); */
2939 /* print_qmarray(Q,0,NULL); */
2940 /* printf("\n\n\n\n"); */
2941
2942 assert(out == 0);
2943 qmarray_free(V); V = NULL;
2944 }
2945 }
2946 else{
2947 fprintf(stderr, "No clear QR/LQ decomposition for type=%s\n",dir);
2948 exit(1);
2949 }
2950 return Q;
2951 }
2952
2953 /***********************************************************
2954 Compute the householder triangularization of a
2955 qmarray. for nodal basis (grid) functions of a fixed
2956 grid
2957
2958 \param[in] dir - type either "QR" or "LQ"
2959 \param[in,out] A - qmarray to triangularize (destroyed in call)
2960 \param[in,out] R - allocated space upper triangular factor
2961 \param[in] grid - grid of locations at which to
2962
2963 \return Q
2964 ***************************************************************/
2965 /* struct Qmarray * */
2966 /* qmarray_householder_simple_grid(char * dir, struct Qmarray * A, double * R, */
2967 /* struct c3Vector * grid) */
2968 /* { */
2969
2970 /* for (size_t ii = 0; ii < A->nrows * A->ncols; ii++){ */
2971 /* assert(A->funcs[ii]->fc == LINELM); */
2972 /* } */
2973
2974 /* struct Qmarray * Q = NULL; */
2975 /* if (strcmp(dir,"QR") == 0){ */
2976 /* Q = qmarray_orth1d_linelm_grid(A->nrows, A->ncols, grid); */
2977 /* int out = qmarray_qr(A,&Q,&R); */
2978 /* assert (out == 0); */
2979 /* } */
2980 /* else if (strcmp(dir, "LQ") == 0){ */
2981 /* struct Qmarray * temp = qmarray_orth1d_linelm_grid(A->nrows,A->ncols,grid); */
2982 /* Q = qmarray_transpose(temp); */
2983 /* int out = qmarray_lq(A,&Q,&R); */
2984 /* assert(out == 0); */
2985 /* qmarray_free(temp); temp = NULL; */
2986 /* } */
2987 /* else{ */
2988 /* fprintf(stderr, "No clear QR/LQ decomposition for type=%s\n",dir); */
2989 /* exit(1); */
2990 /* } */
2991 /* return Q; */
2992 /* } */
2993
2994 /***********************************************************//**
2995 Compute the svd of a quasimatrix array Udiag(lam)vt = A
2996
2997 \param[in,out] A - qmarray to get SVD (destroyed)
2998 \param[in,out] U - qmarray with orthonormal columns
2999 \param[in,out] lam - singular values
3000 \param[in,out] vt - matrix containing right singular vectors
3001 \param[in] opts - options for approximations
3002
3003 \return info - if not == 0 then error
3004 ***************************************************************/
qmarray_svd(struct Qmarray * A,struct Qmarray ** U,double * lam,double * vt,struct OneApproxOpts * opts)3005 int qmarray_svd(struct Qmarray * A, struct Qmarray ** U, double * lam,
3006 double * vt, struct OneApproxOpts * opts)
3007 {
3008
3009 int info = 0;
3010
3011 double * R = calloc_double(A->ncols * A->ncols);
3012 struct Qmarray * temp = qmarray_householder_simple("QR",A,R,opts);
3013
3014 double * u = calloc_double(A->ncols * A->ncols);
3015 svd(A->ncols, A->ncols, A->ncols, R, u,lam,vt);
3016
3017 qmarray_free(*U);
3018 *U = qmam(temp,u,A->ncols);
3019
3020 qmarray_free(temp);
3021 free(u);
3022 free(R);
3023 return info;
3024 }
3025
3026 /***********************************************************//**
3027 Compute the reduced svd of a quasimatrix array
3028 Udiag(lam)vt = A with every singular value greater than delta
3029
3030 \param[in,out] A - qmarray to get SVD (destroyed)
3031 \param[in,out] U - qmarray with orthonormal columns
3032 \param[in,out] lam - singular values
3033 \param[in,out] vt - matrix containing right singular vectors
3034 \param[in,out] delta - threshold
3035 \param[in] opts - approximation options
3036
3037 \return rank - rank of the qmarray
3038
3039 \note
3040 *U*, *lam*, and *vt* are allocated in this function
3041 ***************************************************************/
qmarray_truncated_svd(struct Qmarray * A,struct Qmarray ** U,double ** lam,double ** vt,double delta,struct OneApproxOpts * opts)3042 size_t qmarray_truncated_svd(
3043 struct Qmarray * A, struct Qmarray ** U,
3044 double ** lam, double ** vt, double delta,
3045 struct OneApproxOpts * opts)
3046 {
3047
3048 //int info = 0;
3049 double * R = calloc_double(A->ncols * A->ncols);
3050 struct Qmarray * temp = qmarray_householder_simple("QR",A,R,opts);
3051
3052 double * u = NULL;
3053 //printf("R = \n");
3054 //dprint2d_col(A->ncols,A->ncols,R);
3055 size_t rank = truncated_svd(A->ncols, A->ncols, A->ncols,
3056 R, &u, lam, vt, delta);
3057
3058 qmarray_free(*U);
3059 *U = qmam(temp,u,rank);
3060
3061 qmarray_free(temp);
3062 free(u);
3063 free(R);
3064 return rank;
3065 }
3066
3067 /***********************************************************//**
3068 Compute the location of the (abs) maximum element in aqmarray consisting of 1 dimensional functions
3069
3070 \param[in] qma - qmarray 1
3071 \param[in,out] xmax - x value at maximum
3072 \param[in,out] rowmax - row of maximum
3073 \param[in,out] colmax - column of maximum
3074 \param[in,out] maxval - absolute maximum value
3075 \param[in] optargs - optimization arguments
3076 ***************************************************************/
qmarray_absmax1d(struct Qmarray * qma,double * xmax,size_t * rowmax,size_t * colmax,double * maxval,void * optargs)3077 void qmarray_absmax1d(struct Qmarray * qma, double * xmax, size_t * rowmax,
3078 size_t * colmax, double * maxval, void * optargs)
3079 {
3080 size_t combrowcol;
3081 *maxval = generic_function_array_absmax(qma->nrows * qma->ncols,
3082 1, qma->funcs,
3083 &combrowcol, xmax,optargs);
3084
3085 size_t nsubs = 0;
3086 //printf("combrowcol = %zu \n", combrowcol);
3087 while ((combrowcol+1) > qma->nrows){
3088 nsubs++;
3089 combrowcol -= qma->nrows;
3090 }
3091
3092 *rowmax = combrowcol;
3093 *colmax = nsubs;
3094 //printf("rowmax, colmaxx = (%zu,%zu)\n", *rowmax, *colmax);
3095 }
3096
3097 /***********************************************************//**
3098 Transpose a qmarray
3099
3100 \param a [in] - qmarray
3101
3102 \return b - \f$ b(x) = a(x)^T \f$
3103 ***************************************************************/
3104 struct Qmarray *
qmarray_transpose(struct Qmarray * a)3105 qmarray_transpose(struct Qmarray * a)
3106 {
3107
3108 struct Qmarray * b = qmarray_alloc(a->ncols, a->nrows);
3109 size_t ii, jj;
3110 for (ii = 0; ii < a->nrows; ii++){
3111 for (jj = 0; jj < a->ncols; jj++){
3112 b->funcs[ii*a->ncols+jj] =
3113 generic_function_copy(a->funcs[jj*a->nrows+ii]);
3114 }
3115 }
3116 return b;
3117 }
3118
3119 /***********************************************************//**
3120 Stack two qmarrays horizontally
3121
3122 \param[in] a - qmarray 1
3123 \param[in] b - qmarray 2
3124
3125 \return c - \f$ [a b] \f$
3126 ***************************************************************/
qmarray_stackh(const struct Qmarray * a,const struct Qmarray * b)3127 struct Qmarray * qmarray_stackh(const struct Qmarray * a,const struct Qmarray * b)
3128 {
3129 assert (a != NULL);
3130 assert (b != NULL);
3131 assert (a->nrows == b->nrows);
3132 struct Qmarray * c = qmarray_alloc(a->nrows, a->ncols + b->ncols);
3133 size_t ii;
3134 for (ii = 0; ii < a->nrows * a->ncols; ii++){
3135 c->funcs[ii] = generic_function_copy(a->funcs[ii]);
3136 }
3137 for (ii = 0; ii < b->nrows * b->ncols; ii++){
3138 c->funcs[ii+a->nrows*a->ncols] = generic_function_copy(b->funcs[ii]);
3139 }
3140 return c;
3141 }
3142
3143 /***********************************************************//**
3144 Stack two qmarrays vertically
3145
3146 \param[in] a - qmarray 1
3147 \param[in] b - qmarray 2
3148
3149 \return c : \f$ [a; b] \f$
3150 ***************************************************************/
qmarray_stackv(const struct Qmarray * a,const struct Qmarray * b)3151 struct Qmarray * qmarray_stackv(const struct Qmarray * a, const struct Qmarray * b)
3152 {
3153 assert (a != NULL);
3154 assert (b != NULL);
3155 assert (a->ncols == b->ncols);
3156 struct Qmarray * c = qmarray_alloc(a->nrows+b->nrows, a->ncols);
3157 size_t ii, jj;
3158 for (jj = 0; jj < a->ncols; jj++){
3159 for (ii = 0; ii < a->nrows; ii++){
3160 c->funcs[jj*c->nrows+ii] =
3161 generic_function_copy(a->funcs[jj*a->nrows+ii]);
3162 }
3163 for (ii = 0; ii < b->nrows; ii++){
3164 c->funcs[jj*c->nrows+ii+a->nrows] =
3165 generic_function_copy(b->funcs[jj*b->nrows+ii]);
3166 }
3167 }
3168
3169 return c;
3170 }
3171
3172 /***********************************************************//**
3173 Block diagonal qmarray from two qmarrays
3174
3175 \param a [in] - qmarray 1
3176 \param b [in] - qmarray 2
3177
3178 \return c - \f$ [a 0 ;0 b] \f$
3179 ***************************************************************/
qmarray_blockdiag(struct Qmarray * a,struct Qmarray * b)3180 struct Qmarray * qmarray_blockdiag(struct Qmarray * a, struct Qmarray * b)
3181 {
3182 struct Qmarray * c = qmarray_alloc(a->nrows+b->nrows, a->ncols+b->ncols);
3183 size_t ii, jj;
3184 for (jj = 0; jj < a->ncols; jj++){
3185 for (ii = 0; ii < a->nrows; ii++){
3186 c->funcs[jj*c->nrows+ii] =
3187 generic_function_copy(a->funcs[jj*a->nrows+ii]);
3188 }
3189 for (ii = 0; ii < b->nrows; ii++){
3190 c->funcs[jj*c->nrows+ii+a->nrows] = generic_function_daxpby(0.0,
3191 a->funcs[0],0.0,NULL);
3192 }
3193 }
3194 for (jj = a->ncols; jj < c->ncols; jj++){
3195 for (ii = 0; ii < a->nrows; ii++){
3196 c->funcs[jj*c->nrows+ii] = generic_function_daxpby(0.0,a->funcs[0],0.0,NULL);
3197
3198 }
3199 for (ii = 0; ii < b->nrows; ii++){
3200 c->funcs[jj*c->nrows+ii+a->nrows] =
3201 generic_function_copy(b->funcs[(jj-a->ncols)*b->nrows+ii]);
3202 }
3203 }
3204
3205 return c;
3206 }
3207
3208 /***********************************************************//**
3209 Compute the derivative of every function in the qmarray
3210
3211 \param[in] a - qmarray
3212
3213 \return qmarray of functions representing the derivatives at all values x
3214 ***************************************************************/
qmarray_deriv(struct Qmarray * a)3215 struct Qmarray * qmarray_deriv(struct Qmarray * a)
3216 {
3217
3218 struct Qmarray * b = qmarray_alloc(a->nrows, a->ncols);
3219 size_t ii;
3220 for (ii = 0; ii < a->nrows*a->ncols; ii++){
3221 //printf("function %zu is \n",ii);
3222 //print_generic_function(a->funcs[ii],3,NULL);
3223 b->funcs[ii] = generic_function_deriv(a->funcs[ii]);
3224 //printf("got its deriv\n");
3225 }
3226 return b;
3227 }
3228
3229 /***********************************************************//**
3230 Compute the second derivative of every function in the qmarray
3231
3232 \param[in] a - qmarray
3233
3234 \return qmarray of functions representing the derivatives at all values x
3235 ***************************************************************/
qmarray_dderiv(struct Qmarray * a)3236 struct Qmarray * qmarray_dderiv(struct Qmarray * a)
3237 {
3238
3239 struct Qmarray * b = qmarray_alloc(a->nrows, a->ncols);
3240 size_t ii;
3241 for (ii = 0; ii < a->nrows*a->ncols; ii++){
3242 b->funcs[ii] = generic_function_dderiv(a->funcs[ii]);
3243 }
3244 return b;
3245 }
3246
3247 /***********************************************************//**
3248 Compute the second derivative of every function in the qmarray enforcing periodic boundaries
3249
3250 \param[in] a - qmarray
3251
3252 \return qmarray of functions representing the derivatives at all values x
3253 ***************************************************************/
qmarray_dderiv_periodic(struct Qmarray * a)3254 struct Qmarray * qmarray_dderiv_periodic(struct Qmarray * a)
3255 {
3256
3257 struct Qmarray * b = qmarray_alloc(a->nrows, a->ncols);
3258 size_t ii;
3259 for (ii = 0; ii < a->nrows*a->ncols; ii++){
3260 b->funcs[ii] = generic_function_dderiv_periodic(a->funcs[ii]);
3261 }
3262 return b;
3263 }
3264
3265 /***********************************************************//**
3266 Evaluate the derivative of every function in the qmarray
3267
3268 \param[in] a - qmarray
3269 \param[in] x - location at which to evaluate
3270 \param[in,out] out - evaluations (previously allocated)
3271 ***************************************************************/
qmarray_deriv_eval(const struct Qmarray * a,double x,double * out)3272 void qmarray_deriv_eval(const struct Qmarray * a, double x, double * out)
3273 {
3274
3275
3276 size_t ii;
3277 for (ii = 0; ii < a->nrows*a->ncols; ii++){
3278 out[ii] = generic_function_deriv_eval(a->funcs[ii],x);
3279 }
3280 }
3281
3282
3283 /********************************************************//**
3284 Rounding (thresholding) of quasimatrix array
3285
3286 \param qma [inout] - quasimatrix
3287 \param epsilon [in] - threshold
3288
3289 ***********************************************************/
qmarray_roundt(struct Qmarray ** qma,double epsilon)3290 void qmarray_roundt(struct Qmarray ** qma, double epsilon)
3291 {
3292 size_t ii;
3293 for (ii = 0; ii < (*qma)->nrows * (*qma)->ncols; ii++){
3294 generic_function_roundt(&((*qma)->funcs[ii]),epsilon);
3295 }
3296 }
3297
3298
3299 /***********************************************************//**
3300 Evaluate a qmarray
3301
3302 \param[in] qma - quasimatrix array
3303 \param[in] x - location at which to evaluate
3304 \param[in,out] out - allocated array to store output (qma->nrows * qma->ncols)
3305 ***************************************************************/
qmarray_eval(struct Qmarray * qma,double x,double * out)3306 void qmarray_eval(struct Qmarray * qma, double x, double * out)
3307 {
3308 size_t ii;
3309 for (ii = 0; ii < qma->nrows*qma->ncols; ii++){
3310 out[ii] = generic_function_1d_eval(qma->funcs[ii],x);
3311 }
3312 }
3313
3314 /***********************************************************//**
3315 Get number of parameters describing a particular function
3316 within a quasimatrix array
3317
3318 \param[in] qma - quasimatrix array
3319 \param[in] row - row
3320 \param[in] col - col
3321
3322 \returns number of parameters
3323 ***************************************************************/
qmarray_func_get_nparams(const struct Qmarray * qma,size_t row,size_t col)3324 size_t qmarray_func_get_nparams(const struct Qmarray * qma,
3325 size_t row, size_t col)
3326 {
3327 struct GenericFunction * gf = qma->funcs[row + col * qma->nrows];
3328 return generic_function_get_num_params(gf);
3329 }
3330
3331
3332 /***********************************************************//**
3333 Get number of parameters describing the qmarray
3334
3335 \param[in] qma - quasimatrix array
3336 \param[in,out] nfunc - number of parameters in the
3337 function with most parameters
3338
3339 \returns - total number of parameters
3340 ***************************************************************/
qmarray_get_nparams(const struct Qmarray * qma,size_t * nfunc)3341 size_t qmarray_get_nparams(const struct Qmarray * qma, size_t * nfunc)
3342 {
3343
3344 size_t size = qma->nrows*qma->ncols;
3345 size_t ntot = 0;
3346 size_t maxfun = 0;
3347 for (size_t ii = 0; ii < size; ii++){
3348 size_t nparam = generic_function_get_num_params(qma->funcs[ii]);
3349 ntot += nparam;
3350
3351 if (nparam > maxfun){
3352 maxfun = nparam;
3353 }
3354 }
3355 if (nfunc != NULL){
3356 *nfunc = maxfun;
3357 }
3358
3359 return ntot;
3360 }
3361
3362 /***********************************************************//**
3363 Get Qmarray parameters
3364
3365 \param[in] qma - quasimatrix array whose parameters to update
3366 \param[in,out] param - location to store parameters
3367
3368 \returns total number of parameters
3369 ***************************************************************/
qmarray_get_params(struct Qmarray * qma,double * param)3370 size_t qmarray_get_params(struct Qmarray * qma, double * param)
3371 {
3372 size_t size = qma->nrows*qma->ncols;
3373 size_t onind = 0;
3374 size_t nparam;
3375 for (size_t ii = 0; ii < size; ii++){
3376 if (param == NULL){
3377 nparam = generic_function_get_num_params(qma->funcs[ii]);
3378 }
3379 else{
3380 nparam = generic_function_get_params(qma->funcs[ii],
3381 param+onind);
3382 }
3383 onind += nparam;
3384 }
3385 return onind;
3386 }
3387
3388 /***********************************************************//**
3389 Update parameters of a single generic function
3390
3391 \param[in,out] qma - quasimatrix array whose parameters to update
3392 \param[in] row - row of the function to update
3393 \param[in] col - column of the function
3394 \param[in] nparam - total number of parameters
3395 \param[in] param - parameters with which to update
3396 ***************************************************************/
qmarray_uni_update_params(struct Qmarray * qma,size_t row,size_t col,size_t nparam,const double * param)3397 void qmarray_uni_update_params(struct Qmarray * qma, size_t row, size_t col,
3398 size_t nparam, const double * param)
3399 {
3400
3401
3402 size_t ind = col * qma->nrows + row;
3403 size_t nparam_expect = generic_function_get_num_params(qma->funcs[ind]);
3404 if (nparam_expect != nparam){
3405 fprintf(stderr, "Trying to update a univariate function in\n");
3406 fprintf(stderr, "qmarray_uni_update_params with the wrong number of parameters\n");
3407 exit(1);
3408 }
3409 generic_function_update_params(qma->funcs[ind],nparam,param);
3410 }
3411
3412 /***********************************************************//**
3413 Update Qmarray parameters
3414
3415 \param[in,out] qma - quasimatrix array whose parameters to update
3416 \param[in] nparams - total number of parameters
3417 \param[in] param - parameters with which to update
3418 ***************************************************************/
qmarray_update_params(struct Qmarray * qma,size_t nparams,const double * param)3419 void qmarray_update_params(struct Qmarray * qma, size_t nparams, const double * param)
3420 {
3421
3422 size_t size = qma->nrows*qma->ncols;
3423 size_t onind = 0;
3424 for (size_t ii = 0; ii < size; ii++){
3425 size_t nparam = generic_function_get_num_params(qma->funcs[ii]);
3426 generic_function_update_params(qma->funcs[ii],nparam,param+onind);
3427 onind += nparam;
3428 }
3429 if (onind != nparams){
3430 fprintf(stderr,"Error updating the parameters of a Qmarray\n");
3431 fprintf(stderr,"number of parameters expected is not the\n");
3432 fprintf(stderr,"number of actual parameters\n");
3433 exit(1);
3434 }
3435 }
3436
3437 /***********************************************************//**
3438 Evaluate the gradient of a qmarray with respect to the parameters
3439 of the function at a set of evaluation locations
3440
3441 \param[in] qma - quasimatrix array
3442 \param[in] N - number of locations at which to evaluate
3443 \param[in] x - location at which to evaluate
3444 \param[in] incx - increment between inputs
3445 \param[in,out] out - allocated array to store output (qma->nrows * qma->ncols)
3446 \param[in] incout - increment of output
3447 \param[in,out] grad - compute gradient if not null (store gradient here)
3448 \param[in] incg - incremement of gradient
3449 \param[in,out] workspace - big enough to handle the largest number of params of any one function
3450 ***************************************************************/
qmarray_param_grad_eval(struct Qmarray * qma,size_t N,const double * x,size_t incx,double * out,size_t incout,double * grad,size_t incg,double * workspace)3451 void qmarray_param_grad_eval(struct Qmarray * qma, size_t N,
3452 const double * x, size_t incx,
3453 double * out, size_t incout,
3454 double * grad, size_t incg,
3455 double * workspace)
3456 {
3457 size_t ii;
3458 /* printf("\t cmoooonn\n"); */
3459 size_t size = qma->nrows*qma->ncols;
3460 /* printf("\t evaluate! size=%zu\n",size); */
3461 if (out != NULL){
3462 /* printf("1darray_eval2N\n"); */
3463 generic_function_1darray_eval2N(size,qma->funcs,N,x,incx,out,incout);
3464 /* printf("got 1darray_eval2N incout = %zu\n",incout); */
3465 for (ii = 0; ii < size; ii++){
3466 if (isnan(out[ii])){
3467 fprintf(stderr,"Warning, evaluation in qmarray_param_grad_eval is nan\n");
3468 fprintf(stderr,"ii=%zu,size=%zu,incout=%zu\n",ii,size,incout);
3469 fprintf(stderr,"x = %G\n",x[ii*incx]);
3470 exit(1);
3471 }
3472 else if (isinf(out[ii])){
3473 fprintf(stderr,"Warning, evaluation in qmarray_param_grad_eval inf\n");
3474 exit(1);
3475 }
3476 }
3477 /* printf("finished checking!\n"); */
3478
3479 /* printf("out = "); dprint(size,out); */
3480 }
3481
3482 if (grad != NULL){
3483 for (size_t jj = 0; jj < N; jj++){
3484 size_t onparam = 0;
3485 for (ii = 0; ii < size; ii++){
3486 size_t nparam = generic_function_get_num_params(qma->funcs[ii]);
3487 generic_function_param_grad_eval(qma->funcs[ii],1,x + jj*incx,workspace);
3488 for (size_t kk = 0; kk < nparam; kk++){
3489 // only one element changes so copy everything
3490 for(size_t zz = 0; zz < size; zz++){
3491 grad[onparam*size+zz + jj * incg] = 0.0;
3492 }
3493 // and change the one element
3494 grad[onparam*size+ii + jj * incg] = workspace[kk];
3495 if (isnan(workspace[kk])){
3496 fprintf(stderr,"Warning, gradient in qmarray_param_grad_eval is nan\n");
3497 exit(1);
3498 }
3499 else if (isinf(workspace[kk])){
3500 fprintf(stderr,"Warning, gradient in qmarray_param_grad_eval is inf\n");
3501 exit(1);
3502 }
3503
3504 onparam++;
3505 }
3506 }
3507 }
3508 }
3509 /* printf("done there\n"); */
3510 }
3511
3512 /***********************************************************//**
3513 Evaluate the gradient of a qmarray with respect to the parameters
3514 of the function at a set of evaluation locations
3515
3516 \param[in] qma - quasimatrix array
3517 \param[in] N - number of locations at which to evaluate
3518 \param[in,out] out - allocated array to store output (qma->nrows * qma->ncols)
3519 \param[in] incout - increment of output
3520 \param[in] grad - computed gradient
3521 \param[in] incg - incremement of gradient
3522
3523 ***************************************************************/
qmarray_linparam_eval(struct Qmarray * qma,size_t N,double * out,size_t incout,const double * grad,size_t incg)3524 void qmarray_linparam_eval(struct Qmarray * qma, size_t N,
3525 double * out, size_t incout,
3526 const double * grad, size_t incg)
3527
3528 {
3529
3530 size_t size = qma->nrows*qma->ncols;
3531 size_t nparam,onparam;
3532 size_t ii,jj,indout,indgrad;
3533 double * param;
3534 for (ii = 0; ii < N; ii++){
3535 /* printf("ii = %zu\n",ii); */
3536 onparam=0;
3537 indout = incout*ii;
3538 indgrad = ii*incg;
3539 for (jj = 0; jj < size; jj++){
3540 param = generic_function_get_params_ref(qma->funcs[jj],&nparam);
3541 /* dprint(nparam,grad+ii*incg + onparam); */
3542 out[jj + indout] = cblas_ddot(nparam,param,1,
3543 grad + indgrad + onparam, 1);
3544 onparam += nparam;
3545 }
3546 }
3547 }
3548
3549 /***********************************************************//**
3550 Sum the L2 norms of each function
3551
3552 \param[in] qma - quasimatrix array
3553 \param[in] scale - increment between inputs
3554 \param[in,out] grad - compute gradient if not null (store gradient here)
3555
3556 \returns sum of scale * inner_product of each function
3557 \note
3558 If gradient is not null then it adds scaled values of the new gradients to the
3559 existing gradient. Gradient is stored function first, row second, column third
3560 ***************************************************************/
qmarray_param_grad_sqnorm(struct Qmarray * qma,double scale,double * grad)3561 double qmarray_param_grad_sqnorm(struct Qmarray * qma, double scale, double * grad)
3562 {
3563 /* printf("\t cmoooonn\n"); */
3564 size_t size = qma->nrows*qma->ncols;
3565 double out = 0.0;
3566 if (grad == NULL){
3567 for (size_t ii = 0; ii < size; ii++){
3568 out += scale * generic_function_inner(qma->funcs[ii],qma->funcs[ii]);
3569 }
3570 }
3571 else{
3572 size_t runparam=0;
3573 for (size_t ii = 0; ii < size; ii++){
3574 out += scale * generic_function_inner(qma->funcs[ii],qma->funcs[ii]);
3575 size_t nparam = generic_function_get_num_params(qma->funcs[ii]);
3576 int res = generic_function_squared_norm_param_grad(qma->funcs[ii],
3577 scale, grad+runparam);
3578 assert (res == 0);
3579 runparam += nparam;
3580 }
3581 }
3582
3583 return out;
3584 }
3585
3586 /***********************************************************//**
3587 Evaluate the gradient of a qmarray with respect to the parameters
3588 of the function, store gradient in sparse format, and multiply on the left
3589 by vector
3590
3591 \param[in] qma - quasimatrix array
3592 \param[in] N - number of locations at which to evaluate
3593 \param[in] x - location at which to evaluate
3594 \param[in] incx - increment between inputs
3595 \param[in,out] out - allocated array to store output (qma->nrows * qma->ncols)
3596 \param[in] incout - increment of output
3597 \param[in,out] grad - compute gradient if not null (store gradient here)
3598 \param[in] incg - incremement of gradient between each evaluation
3599 \param[in] left - set of vectors to multiply each gradient
3600 \param[in,out] mult_out - multiplied out values
3601 \param[in] nparam - expected number of parameters
3602
3603
3604 \note
3605 mult_out is a flattened array, increment between values for different
3606 *x* values is nparam * qma->ncols
3607 mult_out[]
3608
3609 ***************************************************************/
qmarray_param_grad_eval_sparse_mult(struct Qmarray * qma,size_t N,const double * x,size_t incx,double * out,size_t incout,double * grad,size_t incg,double * left,double * mult_out,size_t nparam)3610 void qmarray_param_grad_eval_sparse_mult(struct Qmarray * qma, size_t N,
3611 const double * x, size_t incx,
3612 double * out, size_t incout,
3613 double * grad, size_t incg,
3614 double * left, double * mult_out, size_t nparam)
3615 {
3616
3617 /* printf("\t cmoooonn\n"); */
3618 size_t size = qma->nrows*qma->ncols;
3619 /* printf("\t evaluate! size=%zu\n",size); */
3620 if (out != NULL){
3621 generic_function_1darray_eval2N(size,qma->funcs,N,x,incx,out,incout);
3622 }
3623
3624 if (grad != NULL){
3625 size_t inc_out = nparam * qma->ncols;
3626 for (size_t jj = 0; jj < N; jj++){
3627 // zero every output
3628 if (mult_out != NULL){
3629 for (size_t ii = 0; ii < nparam * qma->ncols; ii++){
3630 mult_out[jj * inc_out + ii] = 0.0;
3631 }
3632 }
3633
3634 size_t onnum = 0;
3635 for (size_t ii = 0; ii < qma->ncols; ii++){
3636 for (size_t kk = 0; kk < qma->nrows; kk++){
3637 size_t onfunc = kk + ii * qma->nrows;
3638 size_t nparamf = generic_function_get_num_params(qma->funcs[onfunc]);
3639 generic_function_param_grad_eval(qma->funcs[onfunc],1,
3640 x + jj*incx,
3641 grad + onnum + jj * incg);
3642
3643 if (left != NULL){
3644 /* printf("%zu,%zu,%zu\n",ii,kk,onnum); */
3645 size_t active_left = kk + jj * qma->nrows;
3646
3647
3648 // modifying onnum vector
3649 size_t on_output = onnum * qma->ncols + jj * inc_out;
3650
3651 // only need to update *ith* element
3652 for (size_t ll = 0; ll < nparamf; ll++){
3653 size_t modelem = ll * qma->ncols + ii;
3654 mult_out[on_output + modelem] =
3655 grad[onnum + jj * incg + ll] * left[active_left];
3656 }
3657 }
3658 onnum += nparamf;
3659 }
3660 }
3661 if (mult_out != NULL){
3662 assert (nparam == onnum);
3663 }
3664 }
3665 }
3666 /* printf("done there\n"); */
3667 }
3668
3669
3670 /***********************************************************//**
3671 Interpolate a qmarray using a nodal basis at particular values
3672
3673 \param[in] qma - quasimatrix array
3674 \param[in] N - number of nodes to form nodal basis
3675 \param[in] x - locations of nodal basis
3676
3677 \return Qmarray with interpolated functions at a nodal basis
3678 ***************************************************************/
qmarray_create_nodal(struct Qmarray * qma,size_t N,double * x)3679 struct Qmarray * qmarray_create_nodal(struct Qmarray * qma, size_t N, double * x)
3680 {
3681 struct Qmarray * qmaout = qmarray_alloc(qma->nrows,qma->ncols);
3682 for (size_t ii = 0; ii < qma->nrows * qma->ncols; ii++){
3683 qmaout->funcs[ii] = generic_function_create_nodal(qma->funcs[ii],N,x);
3684 }
3685 return qmaout;
3686 }
3687
3688
3689 /***********************************************************//**
3690 Determine whether kristoffel weighting is active
3691
3692 \param[in] qma - qmarray
3693
3694 \return 1 if active, 0 otherwise
3695
3696 \note It is assumed that if kristoffel is active on one function
3697 it is active on all
3698 ***************************************************************/
qmarray_is_kristoffel_active(const struct Qmarray * qma)3699 int qmarray_is_kristoffel_active(const struct Qmarray * qma)
3700 {
3701
3702 int active = generic_function_is_kristoffel_active(qma->funcs[0]);
3703 for (size_t ii = 1; ii < qma->nrows * qma->ncols; ii++){
3704 if (generic_function_is_kristoffel_active(qma->funcs[ii]) != active){
3705 fprintf(stderr, "Kristoffel weighting cannot be active on some univariate functions\n");
3706 fprintf(stderr, "and inactive on other ones\n");
3707 exit(1);
3708 }
3709 }
3710 return active;
3711 }
3712
qmarray_activate_kristoffel(struct Qmarray * qma)3713 void qmarray_activate_kristoffel(struct Qmarray * qma)
3714 {
3715 for (size_t ii = 0; ii < qma->nrows * qma->ncols; ii++){
3716 generic_function_activate_kristoffel(qma->funcs[ii]);
3717 }
3718 }
3719
qmarray_deactivate_kristoffel(struct Qmarray * qma)3720 void qmarray_deactivate_kristoffel(struct Qmarray * qma)
3721 {
3722 for (size_t ii = 0; ii < qma->nrows * qma->ncols; ii++){
3723 generic_function_deactivate_kristoffel(qma->funcs[ii]);
3724 }
3725 }
3726
3727
3728 /***********************************************************//**
3729 Get the kristoffel normalization factor
3730
3731 \param[in] qma - quasimatrix array
3732 \param[in] x - location at which to obtain normalization
3733
3734 \return normalization factor
3735
3736 \note Assumes that each univariate function has the same weight
3737 ***************************************************************/
qmarray_get_kristoffel_weight(const struct Qmarray * qma,double x)3738 double qmarray_get_kristoffel_weight(const struct Qmarray * qma,
3739 double x)
3740 {
3741 double weight = generic_function_get_kristoffel_weight(qma->funcs[0],x);
3742 /* printf("qmarray kri weight = %G\n",weight); */
3743 return weight;
3744 }
3745
3746
3747
3748
3749 /////////////////////////////////////////////////////////////////////
3750 /////////////////////////////////////////////////////////////////////
3751 /////////////////////////////////////////////////////////////////////
3752 /////////////////////////////////////////////////////////////////////
3753 /////////////////////////////////////////////////////////////////////
3754 /////////////////////////////////////////////////////////////////////
3755 /////////////////////////////////////////////////////////////////////
3756 /////////////////////////////////////////////////////////////////////
3757 /////////////////////////////////////////////////////////////////////
3758
3759
qmarray_qr(struct Qmarray * A,struct Qmarray ** Q,double ** R,struct OneApproxOpts * app)3760 int qmarray_qr(struct Qmarray * A, struct Qmarray ** Q, double ** R,
3761 struct OneApproxOpts * app)
3762 {
3763
3764 size_t ii,kk,ll;
3765
3766 size_t nrows = A->nrows;
3767 size_t ncols = A->ncols;
3768 if ( (*Q) == NULL){
3769 *Q = qmarray_orth1d_columns(nrows,ncols,app);
3770 }
3771
3772 // print_qmarray(*Q,0,NULL);
3773
3774 if ((*R) == NULL){
3775 *R = calloc_double(ncols*ncols);
3776 }
3777
3778 struct Qmarray * V = qmarray_alloc(nrows,ncols);
3779 for (ii = 0; ii < nrows * ncols; ii++){
3780 V->funcs[ii] = generic_function_alloc_base(1);
3781 }
3782
3783 /* printf("\n\n\n ------------------------------ \n"); */
3784 /* printf("in qmarray QR\n"); */
3785 /* for (size_t ii = 0; ii < ncols; ii++){ */
3786 /* struct GenericFunction * gf = qmarray_get_func(*Q,0,ii); */
3787 /* double integral = generic_function_integral(gf); */
3788 /* double norm = generic_function_norm(gf); */
3789 /* printf("\nintegral = %G, norm=%G\n",integral,norm); */
3790 /* /\* print_generic_function(gf,0,NULL); *\/ */
3791 /* } */
3792
3793
3794 double s, rho, alpha, sigma,v;
3795 for (ii = 0; ii < ncols; ii++){
3796 /* printf("ii = %zu \n", ii); */
3797
3798 rho = generic_function_array_norm(nrows,1,A->funcs+ii*nrows);
3799 (*R)[ii*ncols+ii] = rho;
3800 alpha = generic_function_inner_sum(nrows,1,
3801 (*Q)->funcs+ii*nrows,1,
3802 A->funcs+ii*nrows);
3803
3804 /* printf("\t rho=%G, \t alpha=%G\n",rho,alpha); */
3805 if (fabs(alpha) < ZERO){
3806 alpha = 0.0;
3807 s = 0.0;
3808 }
3809 else{
3810 s = -alpha / fabs(alpha);
3811 if (s < 0.0){
3812 generic_function_array_flip_sign(nrows,1,
3813 (*Q)->funcs+ii*nrows);
3814 }
3815 }
3816
3817 // printf("s = %G\n",s);
3818
3819 for (kk = 0; kk < nrows; kk++){
3820 generic_function_weighted_sum_pa(
3821 rho,(*Q)->funcs[ii*nrows+kk],-1.0,
3822 A->funcs[ii*nrows+kk],
3823 &(V->funcs[ii*nrows+kk]));
3824 }
3825
3826 //improve orthogonality
3827 for (size_t zz = 0; zz < ii; zz++){
3828 v = generic_function_inner_sum(nrows,1,
3829 V->funcs+ii*nrows,
3830 1, (*Q)->funcs+zz*nrows);
3831 generic_function_array_axpy(nrows,-v,(*Q)->funcs+zz*nrows,
3832 V->funcs+ii*nrows);
3833 }
3834
3835 sigma = 0.0;
3836 for (kk = 0; kk < nrows; kk++){
3837 sigma += generic_function_inner(V->funcs[ii*nrows+kk],
3838 V->funcs[ii*nrows+kk]);
3839 }
3840
3841 /* printf("\t sigma = %G\n",sigma); */
3842 sigma = sqrt(sigma);
3843 // weird that this line is needed for stabilization
3844 // Need to add reorthoganlization
3845 //sigma = sigma + 1e-16;
3846 if (fabs(sigma) < ZERO){
3847 //printf("less than zero sigma = %G,Zero=%G\n",sigma,ZERO);
3848 //V=E
3849 for (kk = 0; kk < nrows; kk++){
3850 generic_function_free(V->funcs[ii*nrows+kk]);
3851 V->funcs[ii*nrows+kk] = NULL;
3852 V->funcs[ii*nrows+kk] = generic_function_copy(
3853 (*Q)->funcs[ii*nrows+kk]);
3854 //generic_function_copy_pa((*Q)->funcs[ii*nrows+kk],V->funcs[ii*nrows+kk]);
3855 }
3856 }
3857 else{
3858 /* printf("\t scale that bad boy, sigma less? \n"); */
3859 generic_function_array_scale(1.0/sigma,V->funcs+ii*nrows,
3860 nrows);
3861 }
3862
3863 //printf("start inner loop\n");
3864 double ev = generic_function_inner_sum(nrows,1,
3865 (*Q)->funcs+ii*nrows,1,
3866 V->funcs+ii*nrows);
3867 for (kk = ii+1; kk < ncols; kk++){
3868 double temp =
3869 generic_function_inner_sum(nrows,1,
3870 V->funcs+ii*nrows,1,
3871 A->funcs+kk*nrows);
3872 (*R)[kk*ncols+ii]=
3873 generic_function_inner_sum(nrows,1,
3874 (*Q)->funcs+ii*nrows,1,
3875 A->funcs+kk*nrows);
3876 (*R)[kk*ncols+ii] += (-2.0 * ev * temp);
3877 for (ll = 0; ll < nrows; ll++){
3878 int success =
3879 generic_function_axpy(-2.0*temp,
3880 V->funcs[ii*nrows+ll],
3881 A->funcs[kk*nrows+ll]);
3882 if (success == 1){
3883 return 1;
3884 }
3885 //assert (success == 0);
3886 success = generic_function_axpy(
3887 -(*R)[kk*ncols+ii],
3888 (*Q)->funcs[ii*nrows+ll],
3889 A->funcs[kk*nrows+ll]);
3890
3891 if (success == 1){
3892 return 1;
3893 }
3894 // assert (success == 0);
3895 }
3896 }
3897 }
3898
3899 /* printf("\nafter getting the reflectors V = \n"); */
3900 /* for (size_t ii = 0; ii < ncols; ii++){ */
3901 /* struct GenericFunction * gf = qmarray_get_func(V,0,ii); */
3902 /* double integral = generic_function_integral(gf); */
3903 /* double norm = generic_function_norm(gf); */
3904 /* printf("integral = %G, norm=%G\n",integral,norm); */
3905 /* /\* print_generic_function(gf,0,NULL); *\/ */
3906 /* } */
3907
3908 /* printf("\n ------------------------------ \n\n\n"); */
3909 // form Q
3910 size_t counter = 0;
3911 ii = ncols-1;
3912 size_t jj;
3913 while (counter < ncols){
3914
3915 for (jj = ii; jj < ncols; jj++){
3916 double val = generic_function_inner_sum(nrows,1,
3917 (*Q)->funcs+jj*nrows, 1, V->funcs+ii*nrows);
3918
3919 for (ll = 0; ll < nrows; ll++){
3920 generic_function_axpy(-2.0*val,V->funcs[ii*nrows+ll],
3921 (*Q)->funcs[jj*nrows+ll]);
3922
3923 }
3924 }
3925
3926 ii = ii - 1;
3927 counter += 1;
3928 }
3929
3930 qmarray_free(V); V = NULL;
3931 return 0;
3932 }
3933
qmarray_lq(struct Qmarray * A,struct Qmarray ** Q,double ** L,struct OneApproxOpts * app)3934 int qmarray_lq(struct Qmarray * A, struct Qmarray ** Q, double ** L,
3935 struct OneApproxOpts * app)
3936 {
3937
3938 size_t ii,kk,ll;
3939
3940 size_t nrows = A->nrows;
3941 size_t ncols = A->ncols;
3942
3943 if ( (*Q) == NULL){
3944 *Q = qmarray_orth1d_rows(nrows,ncols,app);
3945 }
3946
3947 if ((*L) == NULL){
3948 *L = calloc_double(nrows*nrows);
3949 }
3950
3951 struct Qmarray * V = qmarray_alloc(nrows,ncols);
3952 for (ii = 0; ii < nrows * ncols; ii++){
3953 V->funcs[ii] = generic_function_alloc_base(1);
3954 }
3955
3956 double s, rho, alpha, sigma, v;
3957 for (ii = 0; ii < nrows; ii++){
3958 //printf("ii = %zu \n", ii);
3959
3960 rho = generic_function_array_norm(ncols,nrows,A->funcs+ii);
3961 (*L)[ii*nrows+ii] = rho;
3962 alpha = generic_function_inner_sum(ncols,nrows,(*Q)->funcs+ii,nrows, A->funcs+ii);
3963
3964 //printf("alpha=%G\n",alpha);
3965 if (fabs(alpha) < ZERO){
3966 alpha = 0.0;
3967 s = 0.0;
3968 }
3969 else{
3970 s = -alpha / fabs(alpha);
3971 if (s < 0.0){
3972 generic_function_array_flip_sign(ncols,nrows,(*Q)->funcs+ii);
3973 }
3974 }
3975
3976 for (kk = 0; kk < ncols; kk++){
3977 generic_function_weighted_sum_pa(
3978 rho,(*Q)->funcs[kk*nrows+ii],-1.0,
3979 A->funcs[kk*nrows+ii],
3980 &(V->funcs[kk*nrows+ii]));
3981 }
3982
3983 //improve orthogonality
3984 for (size_t zz = 0; zz < ii; zz++){
3985 v = generic_function_inner_sum(ncols,nrows,
3986 V->funcs+ii,
3987 nrows, (*Q)->funcs+zz);
3988 for (size_t llz = 0; llz < ncols; llz++){
3989 generic_function_axpy(-v,(*Q)->funcs[llz*nrows+zz],
3990 V->funcs[llz*nrows+ii]);
3991 }
3992 }
3993
3994 sigma = 0.0;
3995 for (kk = 0; kk < ncols;kk++){
3996 sigma += generic_function_inner(V->funcs[kk*nrows+ii],
3997 V->funcs[kk*nrows+ii]);
3998 }
3999
4000 sigma = sqrt(sigma);
4001 // weird that this line is needed for stabilization
4002 // Need to add reorthoganlization
4003 //sigma = sigma + ZERO*1e-3;
4004 if (fabs(sigma) < ZERO){
4005 //printf("less than zero sigma = %G,Zero=%G\n",sigma,ZERO);
4006 //V=E
4007 for (kk = 0; kk < ncols; kk++){
4008 generic_function_free(V->funcs[kk*nrows+ii]);
4009 V->funcs[kk*nrows+ii] = NULL;
4010 V->funcs[kk*nrows+ii] = generic_function_copy((*Q)->funcs[kk*nrows+ii]);
4011 //generic_function_copy_pa((*Q)->funcs[ii*nrows+kk],V->funcs[ii*nrows+kk]);
4012 }
4013 }
4014 else{
4015 for (kk = 0; kk < ncols; kk++){
4016 generic_function_scale(1.0/sigma,V->funcs[kk*nrows+ii]);
4017 }
4018 }
4019
4020 //printf("start inner loop\n");
4021 double ev = generic_function_inner_sum(ncols,nrows,(*Q)->funcs+ii,nrows,V->funcs+ii);
4022 for (kk = ii+1; kk < nrows; kk++){
4023 double temp = generic_function_inner_sum(ncols,nrows,V->funcs+ii,nrows,A->funcs+kk);
4024 (*L)[ii*nrows+kk]=generic_function_inner_sum(ncols,nrows,(*Q)->funcs+ii,nrows,A->funcs+kk);
4025 (*L)[ii*nrows+kk] += (-2.0 * ev * temp);
4026 for (ll = 0; ll < ncols; ll++){
4027 int success = generic_function_axpy(-2.0*temp,V->funcs[ll*nrows+ii],A->funcs[ll*nrows+kk]);
4028 assert (success == 0);
4029 success = generic_function_axpy(-(*L)[ii*nrows+kk],
4030 (*Q)->funcs[ll*nrows+ii],A->funcs[ll*nrows+kk]);
4031 assert (success == 0);
4032 }
4033 }
4034 }
4035
4036 // form Q
4037 size_t counter = 0;
4038 ii = nrows-1;
4039 size_t jj;
4040 while (counter < nrows){
4041
4042 for (jj = ii; jj < nrows; jj++){
4043 double val = generic_function_inner_sum(ncols,nrows,
4044 (*Q)->funcs+jj, nrows, V->funcs+ii);
4045
4046 for (ll = 0; ll < ncols; ll++){
4047 generic_function_axpy(-2.0*val,V->funcs[ll*nrows+ii],
4048 (*Q)->funcs[ll*nrows+jj]);
4049
4050 }
4051 }
4052
4053 ii = ii - 1;
4054 counter += 1;
4055 }
4056
4057 qmarray_free(V); V = NULL;
4058 return 0;
4059 }
4060
qmarray_qr_gs(struct Qmarray * A,double ** R)4061 int qmarray_qr_gs(struct Qmarray * A, double ** R)
4062 {
4063 size_t nrows = A->nrows;
4064 size_t ncols = A->ncols;
4065 if ((*R) == NULL){
4066 *R = calloc_double(ncols*ncols);
4067 }
4068
4069 for (size_t ii = 0; ii < ncols; ii++ ){
4070 (*R)[ii*ncols+ii] = generic_function_array_norm(nrows,1,A->funcs + ii*nrows);
4071 if ((*R)[ii*ncols+ii] > ZERO){
4072 generic_function_array_scale(1.0/(*R)[ii*ncols+ii],A->funcs+ii*nrows,nrows);
4073 for (size_t jj = ii+1; jj < ncols; jj++){
4074 (*R)[jj*ncols+ii] = generic_function_inner_sum(nrows,1,A->funcs+ii*nrows,
4075 1, A->funcs + jj*nrows);
4076 generic_function_array_axpy(nrows,-(*R)[jj*ncols+ii],A->funcs+ii*nrows,A->funcs+jj*nrows);
4077 }
4078 }
4079 else{
4080 printf("warning!!\n");
4081 printf("norm = %G\n",(*R)[ii*ncols+ii]);
4082 assert(1 == 0);
4083 }
4084 }
4085 return 0;
4086 }
4087
4088
qmarray_lq_gs(struct Qmarray * A,double ** R)4089 int qmarray_lq_gs(struct Qmarray * A, double ** R)
4090 {
4091 size_t nrows = A->nrows;
4092 size_t ncols = A->ncols;
4093 if ((*R) == NULL){
4094 *R = calloc_double(nrows*nrows);
4095 }
4096
4097 for (size_t ii = 0; ii < nrows; ii++ ){
4098 (*R)[ii*nrows+ii] = generic_function_array_norm(ncols,nrows,A->funcs + ii);
4099 if ((*R)[ii*nrows+ii] > ZERO){
4100 for (size_t jj = 0; jj < ncols; jj++){
4101 generic_function_scale(1.0/(*R)[ii*nrows+ii],A->funcs[jj*nrows+ii]);
4102 }
4103 for (size_t jj = ii+1; jj < nrows; jj++){
4104 (*R)[ii*nrows+jj] = generic_function_inner_sum(ncols,nrows,A->funcs+ii,
4105 nrows, A->funcs + jj);
4106 for (size_t kk = 0; kk < ncols; kk++){
4107 generic_function_axpy(-(*R)[ii*nrows+jj],
4108 A->funcs[kk*nrows+ii],
4109 A->funcs[kk*nrows+jj]);
4110 }
4111 }
4112 }
4113 else{
4114 printf("norm = %G\n",(*R)[ii*nrows+ii]);
4115 assert (1 == 0);
4116 }
4117 }
4118 return 0;
4119 }
4120
4121
print_qmarray(struct Qmarray * qm,size_t prec,void * args,FILE * fp)4122 void print_qmarray(struct Qmarray * qm, size_t prec, void * args, FILE *fp)
4123 {
4124
4125 printf("Quasimatrix Array (%zu,%zu)\n",qm->nrows, qm->ncols);
4126 printf("=========================================\n");
4127 size_t ii,jj;
4128 for (ii = 0; ii < qm->nrows; ii++){
4129 for (jj = 0; jj < qm->ncols; jj++){
4130 printf("(%zu, %zu)\n",ii,jj);
4131 print_generic_function(qm->funcs[jj*qm->nrows+ ii],prec,args,fp);
4132 printf("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
4133 }
4134 }
4135 }
4136