1 /*
2 * gretl -- Gnu Regression, Econometrics and Time-series Library
3 * Copyright (C) 2001 Allin Cottrell and Riccardo "Jack" Lucchetti
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 *
18 */
19
20 #include "libgretl.h"
21 #include "gretl_matrix.h"
22 #include "matrix_extra.h"
23 #include "gretl_cmatrix.h"
24 #include "gretl_normal.h"
25 #include "usermat.h"
26 #include "genparse.h"
27 #include "uservar.h"
28
29 #define MDEBUG 0
30 #define CONTIG_DEBUG 0
31
32 /* mspec convenience macros */
33
34 #define mspec_get_offset(m) (m->lsel.range[0])
35 #define mspec_get_n_elem(m) (m->lsel.range[1])
36
37 #define mspec_set_offset(m,o) (m->lsel.range[0] = o)
38 #define mspec_set_n_elem(m,n) (m->lsel.range[1] = n)
39
40 #define mspec_set_element(m,i) (m->lsel.range[0] = i)
41
42 #define spec_is_single(s) (s->ltype==SEL_SINGLE || s->rtype==SEL_SINGLE)
43 #define spec_single_val(s) (s->ltype==SEL_SINGLE ? s->lsel.range[0] : \
44 s->rsel.range[0])
45
46 #define singleton_left_range(s) (s->ltype == SEL_RANGE && \
47 s->lsel.range[0] == s->lsel.range[1])
48 #define singleton_right_range(s) (s->rtype == SEL_RANGE && \
49 s->rsel.range[0] == s->rsel.range[1])
50
51 #define all_or_null(t) (t == SEL_ALL || t == SEL_NULL)
52
53 #define lhs_is_scalar(s,m) (s->ltype == SEL_ELEMENT || \
54 (m->rows == 1 && s->ltype == SEL_ALL) || \
55 singleton_left_range(s))
56 #define rhs_is_scalar(s,m) (s->rtype == SEL_ELEMENT || \
57 (m->cols == 1 && all_or_null(s->rtype)) || \
58 singleton_right_range(s))
59
60 #define rowmax(s,m) (s->range[1] == MSEL_MAX ? m->rows : s->range[1])
61 #define colmax(s,m) (s->range[1] == MSEL_MAX ? m->cols : s->range[1])
62
63 /* end mspec convenience macros */
64
65 /**
66 * get_matrix_by_name:
67 * @name: name of the matrix.
68 *
69 * Looks up a user-defined matrix by name.
70 *
71 * Returns: pointer to matrix, or %NULL if not found.
72 */
73
get_matrix_by_name(const char * name)74 gretl_matrix *get_matrix_by_name (const char *name)
75 {
76 gretl_matrix *ret = NULL;
77
78 if (name != NULL && *name != '\0') {
79 user_var *u;
80
81 u = get_user_var_of_type_by_name(name, GRETL_TYPE_MATRIX);
82 if (u != NULL) {
83 ret = user_var_get_value(u);
84 }
85 }
86
87 return ret;
88 }
89
90 /**
91 * steal_matrix_by_name:
92 * @name: name of the matrix.
93 *
94 * Looks up a user-defined matrix by name and if found,
95 * grabs the matrix, leaving the matrix pointer on the
96 * named matrix as %NULL.
97 *
98 * Returns: pointer to matrix, or %NULL if not found.
99 */
100
steal_matrix_by_name(const char * name)101 gretl_matrix *steal_matrix_by_name (const char *name)
102 {
103 gretl_matrix *ret = NULL;
104
105 if (name != NULL && *name != '\0') {
106 user_var *u;
107
108 u = get_user_var_of_type_by_name(name, GRETL_TYPE_MATRIX);
109
110 if (u != NULL) {
111 ret = user_var_steal_value(u);
112 }
113 }
114
115 return ret;
116 }
117
118 /**
119 * get_matrix_copy_by_name:
120 * @name: name of the matrix.
121 * @err: location to receive error code;
122 *
123 * Looks up a user-defined matrix by name.
124 *
125 * Returns: a copy of the named matrix, or %NULL if not found.
126 */
127
get_matrix_copy_by_name(const char * name,int * err)128 gretl_matrix *get_matrix_copy_by_name (const char *name, int *err)
129 {
130 gretl_matrix *m = get_matrix_by_name(name);
131 gretl_matrix *ret = NULL;
132
133 if (m == NULL) {
134 *err = E_UNKVAR;
135 } else {
136 ret = gretl_matrix_copy(m);
137 if (ret == NULL) {
138 *err = E_ALLOC;
139 }
140 }
141
142 return ret;
143 }
144
145 /* @s is the user-supplied selection vector.
146 @n is the max possible value of row or col.
147 @pslice is location to receive positive selection
148 list.
149 */
150
handle_vector_exclusion(const gretl_vector * s,int n,int ** pslice)151 static int handle_vector_exclusion (const gretl_vector *s,
152 int n, int **pslice)
153 {
154 int slen = gretl_vector_get_length(s);
155 int i, j, k, nsel = n;
156 int err = 0;
157
158 /* do we need this check here? */
159 for (i=0; i<slen; i++) {
160 k = (int) fabs(s->val[i]);
161 if (k > n) {
162 gretl_errmsg_sprintf(_("Index value %d is out of bounds"), k);
163 return E_DATA;
164 }
165 }
166
167 /* how many of the @n values are selected? */
168 for (i=1; i<=n; i++) {
169 for (j=0; j<slen; j++) {
170 if (s->val[j] == -i) {
171 /* element i is dropped */
172 nsel--;
173 break;
174 }
175 }
176 }
177
178 if (nsel == 0) {
179 *pslice = gretl_null_list();
180 } else {
181 /* compose positive selection list */
182 int *slice = gretl_list_new(nsel);
183 int excluded, k = 1;
184
185 if (slice != NULL) {
186 for (i=1; i<=n; i++) {
187 excluded = 0;
188 for (j=0; j<slen; j++) {
189 if (s->val[j] == -i) {
190 excluded = 1;
191 break;
192 }
193 }
194 if (!excluded) {
195 slice[k++] = i;
196 }
197 }
198 *pslice = slice;
199 }
200 }
201
202 if (*pslice == NULL) {
203 err = E_ALLOC;
204 }
205
206 return err;
207 }
208
bad_sel_vector(const gretl_vector * v,int n)209 static int bad_sel_vector (const gretl_vector *v, int n)
210 {
211 int i, k, len = gretl_vector_get_length(v);
212 int exclude = v->val[0] < 0;
213
214 for (i=0; i<len; i++) {
215 k = exclude ? -v->val[i] : v->val[i];
216 if (k < 1 || k > n) {
217 gretl_errmsg_sprintf(_("Index value %d is out of bounds"), k);
218 return E_INVARG;
219 }
220 }
221
222 return 0;
223 }
224
bad_sel_range(int * range,int n)225 static int bad_sel_range (int *range, int n)
226 {
227 int i, k, err = 0;
228
229 for (i=0; i<2; i++) {
230 k = range[i];
231 if (k != MSEL_MAX && (k < 1 || k > n)) {
232 err = E_INVARG;
233 if (k > n && k > 2000000000) {
234 k -= IDX_TBD - n;
235 if (k > 0) {
236 range[i] = k;
237 err = 0;
238 }
239 }
240 if (err) {
241 gretl_errmsg_sprintf(_("Index value %d is out of bounds"), k);
242 break;
243 }
244 }
245 }
246
247 return err;
248 }
249
bad_sel_single(int * pk,int n)250 static int bad_sel_single (int *pk, int n)
251 {
252 int err = 0;
253
254 if (*pk != MSEL_MAX && (*pk < 1 || *pk > n)) {
255 err = E_INVARG;
256 if (*pk > n && *pk > 2000000000) {
257 int k = *pk + (n - IDX_TBD);
258
259 *pk = k;
260 if (k > 0) {
261 err = 0;
262 }
263 }
264 if (err) {
265 gretl_errmsg_sprintf(_("Index value %d is out of bounds"), *pk);
266 }
267 }
268
269 return err;
270 }
271
272 /* convert a matrix subspec component into list of rows
273 or columns: n is the maximal dimension of selectable items
274 */
275
mspec_make_list(int type,union msel * sel,int n,int * err)276 int *mspec_make_list (int type, union msel *sel, int n,
277 int *err)
278 {
279 int *slice = NULL;
280 int single_exclude = 0;
281 int i, ns = 0;
282
283 if (type == SEL_ALL || type == SEL_NULL) {
284 return NULL;
285 }
286
287 if (type == SEL_MATRIX) {
288 if (sel->m == NULL) {
289 gretl_errmsg_set(_("Range is non-positive!"));
290 *err = E_DATA;
291 } else if (sel->m->val[0] < 0) {
292 *err = handle_vector_exclusion(sel->m, n, &slice);
293 return slice; /* we're done */
294 } else {
295 ns = gretl_vector_get_length(sel->m);
296 }
297 } else {
298 /* range or single exclusion */
299 int sr0 = sel->range[0];
300 int sr1 = sel->range[1];
301
302 if (sr1 == MSEL_MAX) {
303 sr1 = sel->range[1] = n;
304 }
305 if (sr0 < 0 && sr1 == sr0) {
306 /* excluding a single row or column? */
307 sr0 = -sr0;
308 if (sr0 > n) {
309 gretl_errmsg_sprintf(_("Index value %d is out of bounds"),
310 sr0);
311 *err = E_DATA;
312 } else {
313 ns = n - 1;
314 single_exclude = sr0;
315 }
316 } else {
317 *err = bad_sel_range(sel->range, n);
318 if (!*err) {
319 ns = sel->range[1] - sel->range[0] + 1;
320 if (ns <= 0) {
321 gretl_errmsg_sprintf(_("Range %d to %d is non-positive!"),
322 sel->range[0], sel->range[1]);
323 *err = E_DATA;
324 }
325 }
326 }
327 }
328
329 if (!*err) {
330 slice = gretl_list_new(ns);
331 if (slice == NULL) {
332 *err = E_ALLOC;
333 }
334 }
335
336 if (!*err) {
337 /* compose and check slice */
338 if (single_exclude > 0) {
339 int k = 1;
340
341 for (i=1; i<=slice[0]; i++) {
342 if (i == single_exclude) {
343 k++;
344 }
345 slice[i] = k++;
346 }
347 } else {
348 for (i=0; i<slice[0]; i++) {
349 if (type == SEL_MATRIX) {
350 slice[i+1] = sel->m->val[i];
351 } else {
352 slice[i+1] = sel->range[0] + i;
353 }
354 }
355 }
356
357 for (i=1; i<=slice[0] && !*err; i++) {
358 if (slice[i] < 1 || slice[i] > n) {
359 gretl_errmsg_sprintf(_("Index value %d is out of bounds"),
360 slice[i]);
361 *err = 1;
362 }
363 }
364 }
365
366 if (*err) {
367 free(slice);
368 slice = NULL;
369 }
370
371 return slice;
372 }
373
set_element_index(matrix_subspec * spec,const gretl_matrix * m)374 static int set_element_index (matrix_subspec *spec,
375 const gretl_matrix *m)
376
377 {
378 int i = spec->lsel.range[0];
379 int j = spec->rsel.range[0];
380 int k;
381
382 if (spec->ltype == SEL_ALL && spec->rtype == SEL_ALL) {
383 k = 0;
384 } else if (spec->ltype == SEL_ELEMENT ||
385 (spec->ltype == SEL_RANGE && spec->rtype == SEL_RANGE)) {
386 k = (j-1) * m->rows + (i-1);
387 } else {
388 if (spec->ltype == SEL_CONTIG) {
389 fprintf(stderr, "CONTIG on left\n");
390 }
391 if (spec->rtype == SEL_CONTIG) {
392 fprintf(stderr, "CONTIG on right\n");
393 }
394 k = i > j ? (i-1) : (j-1);
395 }
396
397 #if CONTIG_DEBUG
398 fprintf(stderr, "set_element_index: i=%d, j=%d -> k=%d\n",
399 i, j, k);
400 #endif
401
402 /* record single (zero-based) index */
403 mspec_set_offset(spec, k);
404 mspec_set_n_elem(spec, 1);
405
406 return 0;
407 }
408
spec_check_dimensions(matrix_subspec * spec,const gretl_matrix * m)409 static int spec_check_dimensions (matrix_subspec *spec,
410 const gretl_matrix *m)
411 {
412 int lt = spec->ltype;
413 int rt = spec->rtype;
414 int err = 0;
415
416 /* check spec->lsel against rows */
417 if (lt == SEL_MATRIX) {
418 err = bad_sel_vector(spec->lsel.m, m->rows);
419 } else if (lt == SEL_RANGE || lt == SEL_ELEMENT) {
420 err = bad_sel_range(spec->lsel.range, m->rows);
421 }
422 if (!err) {
423 /* check spec->rsel against cols */
424 if (rt == SEL_MATRIX) {
425 err = bad_sel_vector(spec->rsel.m, m->cols);
426 } else if (rt == SEL_RANGE || rt == SEL_ELEMENT) {
427 err = bad_sel_range(spec->rsel.range, m->cols);
428 }
429 }
430
431 return err;
432 }
433
434 #if CONTIG_DEBUG
435 # define IN_USERMAT
436 # include "mspec_debug.c"
437 #endif
438
439 /* When we come here we've already screened out DIAG
440 selection and single-element selection.
441 */
442
submatrix_contig(matrix_subspec * spec,const gretl_matrix * m)443 static int submatrix_contig (matrix_subspec *spec, const gretl_matrix *m)
444 {
445 int contig = 0;
446
447 if (spec->ltype == SEL_MATRIX || spec->rtype == SEL_MATRIX ||
448 spec->ltype == SEL_EXCL || spec->rtype == SEL_EXCL) {
449 ; /* cannot treat as contig */
450 } else if (m->rows == 1 || m->cols == 1) {
451 /* must be contig */
452 contig = 1;
453 } else if (singleton_right_range(spec)) {
454 /* single column selected */
455 contig = 1;
456 } else if ((spec->rtype == SEL_ALL || spec->rtype == SEL_RANGE) &&
457 spec->ltype == SEL_ALL) {
458 /* range of columns, all rows */
459 contig = 1;
460 }
461
462 return contig;
463 }
464
465 /* Determine the offset from the start of m->val, and the
466 number of elements, in a matrix slice that has been found
467 to consist of contiguous data.
468 */
469
get_offset_nelem(matrix_subspec * spec,const gretl_matrix * m,int * poff,int * pn)470 static int get_offset_nelem (matrix_subspec *spec,
471 const gretl_matrix *m,
472 int *poff, int *pn)
473 {
474 SelType stype = spec->ltype;
475 union msel *sel = &spec->lsel;
476 int offset = 0, nelem = 0;
477
478 if (m->cols == 1) {
479 /* column vector */
480 if (stype == SEL_ALL) {
481 nelem = m->rows;
482 } else {
483 int rmax = rowmax(sel, m);
484
485 offset = sel->range[0] - 1;
486 nelem = rmax - sel->range[0] + 1;
487 }
488 } else if (m->rows == 1) {
489 /* row vector */
490 if (spec->rtype != SEL_NULL) {
491 stype = spec->rtype;
492 sel = &spec->rsel;
493 }
494 if (stype == SEL_ALL) {
495 nelem = m->cols;
496 } else {
497 int cmax = colmax(sel, m);
498
499 offset = sel->range[0] - 1;
500 nelem = cmax - sel->range[0] + 1;
501 }
502 } else if (singleton_right_range(spec)) {
503 /* a single matrix column */
504 if (stype == SEL_ALL) {
505 nelem = m->rows;
506 } else {
507 int rmax = rowmax(sel, m);
508
509 offset = sel->range[0] - 1;
510 nelem = rmax - sel->range[0] + 1;
511 }
512 offset += m->rows * (spec->rsel.range[0] - 1);
513 } else {
514 /* multiple columns, all rows */
515 if (spec->rtype == SEL_ALL) {
516 nelem = m->rows * m->cols;
517 } else {
518 int cmax, nc;
519
520 sel = &spec->rsel;
521 cmax = colmax(sel, m);
522 nc = cmax - sel->range[0] + 1;
523 offset = m->rows * (sel->range[0] - 1);
524 nelem = m->rows * nc;
525 }
526 }
527
528 *poff = offset;
529 *pn = nelem;
530
531 return 0;
532 }
533
534 /* to make clear that spec->lsel pertains to columns, move
535 it to the right position */
536
commute_selectors(matrix_subspec * spec)537 static void commute_selectors (matrix_subspec *spec)
538 {
539 spec->rtype = spec->ltype;
540 spec->rsel = spec->lsel;
541 memset(&spec->lsel, 0, sizeof spec->lsel);
542 spec->ltype = SEL_ALL;
543 }
544
545 /* Catch the case of an implicit column or row specification for
546 a sub-matrix of an (n x 1) or (1 x m) matrix; also catch the
547 error of giving just one row/col spec for a matrix that has
548 more than one row and more than one column.
549 */
550
check_matrix_subspec(matrix_subspec * spec,const gretl_matrix * m)551 int check_matrix_subspec (matrix_subspec *spec, const gretl_matrix *m)
552 {
553 int veclen, err = 0;
554
555 if (spec->checked) {
556 return 0;
557 } else if (spec->ltype == SEL_REAL && !m->is_complex) {
558 spec->ltype = spec->rtype = SEL_ALL;
559 goto check_done;
560 } else if (is_sel_dummy(spec->ltype)) {
561 /* SEL_DIAG, SEL_UPPER, SEL_LOWER, SEL_IMAG:
562 nothing to be done here?
563 */
564 goto check_done;
565 }
566
567 #if CONTIG_DEBUG
568 subspec_debug_print(spec, m);
569 #endif
570
571 if (spec->rtype == SEL_NULL && m->rows == 1) {
572 /* row vector: transfer spec to column dimension */
573 commute_selectors(spec);
574 }
575
576 veclen = gretl_vector_get_length(m);
577
578 if (spec->rtype == SEL_NULL && m->rows > 1 && m->cols > 1) {
579 gretl_errmsg_set(_("Ambiguous matrix index"));
580 return E_DATA;
581 }
582
583 if (spec_is_single(spec)) {
584 int k = spec_single_val(spec);
585
586 err = bad_sel_single(&k, veclen);
587 if (!err) {
588 spec->ltype = spec->rtype = SEL_ELEMENT;
589 mspec_set_offset(spec, k - 1);
590 mspec_set_n_elem(spec, 1);
591 }
592 /* nothing more to do */
593 goto check_done;
594 }
595
596 err = spec_check_dimensions(spec, m);
597 if (err) {
598 return err;
599 }
600
601 if (m->rows == 0 || m->cols == 0) {
602 /* nothing more to do */
603 goto check_done;
604 }
605
606 if (lhs_is_scalar(spec, m) && rhs_is_scalar(spec, m)) {
607 /* we're looking at just one element */
608 set_element_index(spec, m);
609 spec->ltype = spec->rtype = SEL_ELEMENT;
610 goto check_done;
611 }
612
613 if (submatrix_contig(spec, m)) {
614 int offset, nelem;
615
616 get_offset_nelem(spec, m, &offset, &nelem);
617 spec->ltype = SEL_CONTIG;
618 mspec_set_offset(spec, offset);
619 mspec_set_n_elem(spec, nelem);
620 if (offset < 0 || nelem <= 0) {
621 fprintf(stderr, "*** offset = %d, nelem = %d (m: %dx%d) ***\n",
622 offset, nelem, m->rows, m->cols);
623 }
624 }
625
626 check_done:
627
628 if (!err) {
629 spec->checked = 1;
630 }
631
632 return err;
633 }
634
mspec_get_string(matrix_subspec * spec,int i)635 const char *mspec_get_string (matrix_subspec *spec, int i)
636 {
637 if (spec == NULL || i < 0 || i > 1) {
638 return NULL;
639 } else if (i == 0) {
640 return spec->ltype != SEL_STR ? NULL : spec->lsel.str;
641 } else {
642 return spec->rtype != SEL_STR ? NULL : spec->rsel.str;
643 }
644 }
645
get_slices(matrix_subspec * spec,const gretl_matrix * M)646 static int get_slices (matrix_subspec *spec,
647 const gretl_matrix *M)
648 {
649 int err = 0;
650
651 spec->rslice = mspec_make_list(spec->ltype, &spec->lsel,
652 M->rows, &err);
653 if (!err) {
654 spec->cslice = mspec_make_list(spec->rtype, &spec->rsel,
655 M->cols, &err);
656 }
657
658 #if MDEBUG
659 if (err) {
660 fprintf(stderr, "usermat: get_slices, err=%d\n", err);
661 }
662 #endif
663
664 return err;
665 }
666
assign_scalar_to_submatrix(gretl_matrix * M,const gretl_matrix * S,double x,matrix_subspec * spec)667 int assign_scalar_to_submatrix (gretl_matrix *M,
668 const gretl_matrix *S,
669 double x,
670 matrix_subspec *spec)
671 {
672 double complex z = NADBL;
673 int mr = M->rows;
674 int mc = M->cols;
675 int i, err = 0;
676
677 if (spec == NULL) {
678 fprintf(stderr, "matrix_replace_submatrix: spec is NULL!\n");
679 return E_DATA;
680 }
681
682 if (M->is_complex) {
683 z = (S != NULL)? S->z[0] : x;
684 }
685
686 if (spec->ltype == SEL_CONTIG) {
687 int ini = mspec_get_offset(spec);
688 int fin = ini + mspec_get_n_elem(spec);
689
690 for (i=ini; i<fin; i++) {
691 if (M->is_complex) {
692 M->z[i] = z;
693 } else {
694 M->val[i] = x;
695 }
696 }
697 return 0;
698 }
699
700 if (is_sel_dummy(spec->ltype)) {
701 if (M->is_complex) {
702 err = gretl_matrix_set_part(M, S, x, spec->ltype);
703 } else {
704 err = gretl_matrix_set_part(M, NULL, x, spec->ltype);
705 }
706 return err; /* we're done */
707 }
708
709 if (spec->rslice == NULL && spec->cslice == NULL) {
710 /* parse @spec into lists of affected rows and columns */
711 err = get_slices(spec, M);
712 }
713
714 if (!err) {
715 int sr = (spec->rslice == NULL)? mr : spec->rslice[0];
716 int sc = (spec->cslice == NULL)? mc : spec->cslice[0];
717 int j, l, k = 0;
718 int mi, mj;
719
720 for (i=0; i<sr; i++) {
721 mi = (spec->rslice == NULL)? k++ : spec->rslice[i+1] - 1;
722 l = 0;
723 for (j=0; j<sc; j++) {
724 mj = (spec->cslice == NULL)? l++ : spec->cslice[j+1] - 1;
725 if (M->is_complex) {
726 gretl_cmatrix_set(M, mi, mj, z);
727 } else {
728 gretl_matrix_set(M, mi, mj, x);
729 }
730 }
731 }
732 }
733
734 return err;
735 }
736
matrix_subspec_new(void)737 matrix_subspec *matrix_subspec_new (void)
738 {
739 matrix_subspec *spec = calloc(1, sizeof *spec);
740
741 if (spec != NULL) {
742 spec->rslice = spec->cslice = NULL;
743 }
744
745 return spec;
746 }
747
contig_cols(matrix_subspec * spec,const gretl_matrix * m)748 static int contig_cols (matrix_subspec *spec,
749 const gretl_matrix *m)
750 {
751 if (spec->rtype == SEL_ALL) {
752 return m->cols;
753 } else if (spec->rtype == SEL_RANGE) {
754 int r1 = spec->rsel.range[1];
755 int cmax = r1 == MSEL_MAX ? m->cols : r1;
756
757 return cmax - spec->rsel.range[0] + 1;
758 } else {
759 return 1;
760 }
761 }
762
transcribe_cols_8(gretl_matrix * M,const gretl_matrix * S,const int * cslice,int sscalar)763 static void transcribe_cols_8 (gretl_matrix *M,
764 const gretl_matrix *S,
765 const int *cslice,
766 int sscalar)
767 {
768 int i, j, mcol, nr = M->rows;
769
770 if (sscalar) {
771 /* write RHS scalar into all rows of each selected col */
772 for (j=1; j<=cslice[0]; j++) {
773 mcol = cslice[j] - 1;
774 for (i=0; i<nr; i++) {
775 gretl_matrix_set(M, i, mcol, S->val[0]);
776 }
777 }
778 } else {
779 /* zap from (cols of) S into selected cols of M */
780 double *xtarg, *xsrc = S->val;
781 size_t colsize = nr * sizeof *M->val;
782
783 for (j=1; j<=cslice[0]; j++) {
784 mcol = cslice[j] - 1;
785 xtarg = M->val + mcol * nr;
786 memcpy(xtarg, xsrc, colsize);
787 xsrc += nr;
788 }
789 }
790 }
791
transcribe_cols_16(gretl_matrix * M,const gretl_matrix * S,const int * cslice,int sscalar)792 static void transcribe_cols_16 (gretl_matrix *M,
793 const gretl_matrix *S,
794 const int *cslice,
795 int sscalar)
796 {
797 int i, j, mcol, nr = M->rows;
798
799 if (sscalar) {
800 /* write RHS scalar into all rows of each selected col */
801 double complex z = S->is_complex ? S->z[0] : S->val[0];
802
803 for (j=1; j<=cslice[0]; j++) {
804 mcol = cslice[j] - 1;
805 for (i=0; i<nr; i++) {
806 gretl_cmatrix_set(M, i, mcol, z);
807 }
808 }
809 } else {
810 /* zap from (cols of) S into selected cols of M */
811 double complex *ztarg, *zsrc = S->z;
812 double *xsrc = S->val;
813 size_t colsize = nr * sizeof *M->z;
814
815 for (j=1; j<=cslice[0]; j++) {
816 mcol = cslice[j] - 1;
817 ztarg = M->z + mcol * nr;
818 if (!S->is_complex) {
819 for (i=0; i<nr; i++) {
820 ztarg[i] = xsrc[i];
821 }
822 xsrc += nr;
823 } else {
824 memcpy(ztarg, zsrc, colsize);
825 zsrc += nr;
826 }
827 }
828 }
829 }
830
831 /* @M is the target for partial replacement, @S is the source to
832 substitute, and @spec tells how/where to make the
833 substitution.
834 */
835
matrix_replace_submatrix(gretl_matrix * M,const gretl_matrix * S,matrix_subspec * spec)836 int matrix_replace_submatrix (gretl_matrix *M,
837 const gretl_matrix *S,
838 matrix_subspec *spec)
839 {
840 int mr = gretl_matrix_rows(M);
841 int mc = gretl_matrix_cols(M);
842 int sr = gretl_matrix_rows(S);
843 int sc = gretl_matrix_cols(S);
844 int sscalar = 0;
845 int i, j;
846 int err = 0;
847
848 if (!M->is_complex && S->is_complex) {
849 fputs("matrix_replace_submatrix: M is real but S is complex\n", stderr);
850 return E_MIXED;
851 }
852
853 #if MDEBUG
854 fprintf(stderr, "\nmatrix_replace_submatrix\n");
855 #endif
856
857 if (spec == NULL) {
858 fprintf(stderr, "matrix_replace_submatrix: spec is NULL!\n");
859 return E_DATA;
860 }
861
862 if (is_sel_dummy(spec->ltype)) {
863 return gretl_matrix_set_part(M, S, 0, spec->ltype);
864 }
865
866 if (spec->ltype == SEL_CONTIG) {
867 int ini = mspec_get_offset(spec);
868 int n = mspec_get_n_elem(spec);
869 int ccols = contig_cols(spec, M);
870
871 if (M->is_complex && S->rows == 1 && S->cols == 1) {
872 for (i=0; i<n; i++) {
873 M->z[ini + i] = S->z[0];
874 }
875 } else if (S->rows * S->cols != n) {
876 err = E_NONCONF;
877 } else if (ccols > 1 && M->rows != S->rows) {
878 err = E_NONCONF;
879 } else if (M->is_complex && !S->is_complex) {
880 /* can't use memcpy here! */
881 for (i=0; i<n; i++) {
882 M->z[ini + i] = S->val[i];
883 }
884 } else if (M->is_complex) {
885 memcpy(M->z + ini, S->z, n * sizeof *M->z);
886 } else {
887 memcpy(M->val + ini, S->val, n * sizeof *M->val);
888 }
889 return err;
890 }
891
892 if (sr > mr || sc > mc) {
893 /* the replacement matrix won't fit into M */
894 fprintf(stderr, "matrix_replace_submatrix: target is %d x %d but "
895 "replacement part is %d x %d\n", mr, mc, sr, sc);
896 return E_NONCONF;
897 }
898
899 if (spec->rslice == NULL && spec->cslice == NULL) {
900 /* parse mspec into lists of affected rows and columns */
901 #if MDEBUG
902 fprintf(stderr, "calling get_slices\n");
903 #endif
904 err = get_slices(spec, M);
905 if (err) {
906 return err;
907 }
908 }
909
910 #if MDEBUG
911 printlist(spec->rslice, "rslice (rows list)");
912 printlist(spec->cslice, "cslice (cols list)");
913 fprintf(stderr, "orig M = %d x %d, S = %d x %d\n", mr, mc, sr, sc);
914 #endif
915
916 if (sr == 1 && sc == 1) {
917 /* the replacement is a scalar */
918 sscalar = 1;
919 sr = (spec->rslice == NULL)? mr : spec->rslice[0];
920 sc = (spec->cslice == NULL)? mc : spec->cslice[0];
921 } else if (spec->rslice != NULL && spec->rslice[0] != sr) {
922 fprintf(stderr, "mspec has %d rows but substitute matrix has %d\n",
923 spec->rslice[0], sr);
924 err = E_NONCONF;
925 } else if (spec->cslice != NULL && spec->cslice[0] != sc) {
926 fprintf(stderr, "mspec has %d cols but substitute matrix has %d\n",
927 spec->cslice[0], sc);
928 err = E_NONCONF;
929 }
930
931 if (!err && spec->rslice == NULL && spec->cslice != NULL) {
932 /* the target is just specified by column(s) */
933 if (M->is_complex) {
934 transcribe_cols_16(M, S, spec->cslice, sscalar);
935 } else {
936 transcribe_cols_8(M, S, spec->cslice, sscalar);
937 }
938 } else if (!err) {
939 /* the general case, no special shortcuts */
940 double complex z = 0;
941 double x = 0;
942 int l, k = 0;
943 int mi, mj;
944
945 if (sscalar && S->is_complex) {
946 z = S->z[0];
947 } else if (sscalar) {
948 x = S->val[0];
949 z = x; /* in case M is complex */
950 }
951
952 for (j=0; j<sc; j++) {
953 mj = (spec->cslice == NULL)? k++ : spec->cslice[j+1] - 1;
954 l = 0;
955 for (i=0; i<sr; i++) {
956 mi = (spec->rslice == NULL)? l++ : spec->rslice[i+1] - 1;
957 if (!sscalar) {
958 if (S->is_complex) {
959 z = gretl_cmatrix_get(S, i, j);
960 } else {
961 x = gretl_matrix_get(S, i, j);
962 }
963 }
964 if (M->is_complex) {
965 gretl_cmatrix_set(M, mi, mj, z);
966 } else {
967 gretl_matrix_set(M, mi, mj, x);
968 }
969 }
970 }
971 }
972
973 return err;
974 }
975
matrix_get_submatrix(const gretl_matrix * M,matrix_subspec * spec,int prechecked,int * err)976 gretl_matrix *matrix_get_submatrix (const gretl_matrix *M,
977 matrix_subspec *spec,
978 int prechecked,
979 int *err)
980 {
981 gretl_matrix *S = NULL;
982 int r, c;
983
984 if (M == NULL || spec == NULL) {
985 *err = E_DATA;
986 return NULL;
987 }
988
989 if (!prechecked) {
990 *err = check_matrix_subspec(spec, M);
991 if (*err) {
992 return NULL;
993 }
994 }
995
996 if (spec->ltype == SEL_DIAG) {
997 return gretl_matrix_get_diagonal(M, err);
998 } else if (spec->ltype == SEL_UPPER || spec->ltype == SEL_LOWER) {
999 int upper = (spec->ltype == SEL_UPPER);
1000
1001 return gretl_matrix_get_triangle(M, upper, err);
1002 } else if (spec->ltype == SEL_REAL || spec->ltype == SEL_IMAG) {
1003 int im = (spec->ltype == SEL_IMAG);
1004
1005 return gretl_cmatrix_extract(M, im, err);
1006 } else if (spec->ltype == SEL_CONTIG) {
1007 return matrix_get_chunk(M, spec, err);
1008 } else if (spec->ltype == SEL_ELEMENT) {
1009 int i = mspec_get_element(spec);
1010
1011 if (M->is_complex) {
1012 S = gretl_cmatrix_from_scalar(M->z[i], err);
1013 } else {
1014 S = gretl_matrix_from_scalar(M->val[i]);
1015 }
1016 if (S == NULL) {
1017 *err = E_ALLOC;
1018 }
1019 return S;
1020 }
1021
1022 if (spec->rslice == NULL && spec->cslice == NULL) {
1023 *err = get_slices(spec, M);
1024 if (*err) {
1025 return NULL;
1026 }
1027 }
1028
1029 #if MDEBUG
1030 printlist(spec->rslice, "rslice");
1031 printlist(spec->cslice, "cslice");
1032 fprintf(stderr, "M = %d x %d\n", M->rows, M->cols);
1033 #endif
1034
1035 r = (spec->rslice == NULL)? M->rows : spec->rslice[0];
1036 c = (spec->cslice == NULL)? M->cols : spec->cslice[0];
1037
1038 S = gretl_matching_matrix_new(r, c, M);
1039
1040 if (S == NULL) {
1041 *err = E_ALLOC;
1042 } else {
1043 int j, mj;
1044
1045 if (spec->cslice != NULL && spec->rslice == NULL) {
1046 /* copying entire columns */
1047 if (M->is_complex) {
1048 double complex *dest = S->z;
1049 size_t csize = r * sizeof *dest;
1050
1051 for (j=0; j<c; j++) {
1052 mj = spec->cslice[j+1] - 1;
1053 memcpy(dest, M->z + mj * r, csize);
1054 dest += r;
1055 }
1056 } else {
1057 double *dest = S->val;
1058 size_t csize = r * sizeof *dest;
1059
1060 for (j=0; j<c; j++) {
1061 mj = spec->cslice[j+1] - 1;
1062 memcpy(dest, M->val + mj * r, csize);
1063 dest += r;
1064 }
1065 }
1066 } else {
1067 int i, mi, l, k = 0;
1068 double complex z;
1069 double x;
1070
1071 for (j=0; j<c; j++) {
1072 mj = (spec->cslice == NULL)? k++ : spec->cslice[j+1] - 1;
1073 l = 0;
1074 for (i=0; i<r; i++) {
1075 mi = (spec->rslice == NULL)? l++ : spec->rslice[i+1] - 1;
1076 if (M->is_complex) {
1077 z = gretl_cmatrix_get(M, mi, mj);
1078 gretl_cmatrix_set(S, i, j, z);
1079 } else {
1080 x = gretl_matrix_get(M, mi, mj);
1081 gretl_matrix_set(S, i, j, x);
1082 }
1083 }
1084 }
1085 }
1086 }
1087
1088 if (S != NULL) {
1089 /* try transcribing metadata on @M if applicable */
1090 if (S->rows == M->rows && gretl_matrix_is_dated(M)) {
1091 gretl_matrix_transcribe_obs_info(S, M);
1092 }
1093 if (S->cols == M->cols) {
1094 const char **cnames = gretl_matrix_get_colnames(M);
1095
1096 if (cnames != NULL) {
1097 char **cpy = strings_array_dup((char **) cnames, M->cols);
1098
1099 if (cpy != NULL) {
1100 gretl_matrix_set_colnames(S, cpy);
1101 }
1102 }
1103 }
1104 }
1105
1106 return S;
1107 }
1108
1109 /* Return the element of @M that occupies 0-based position @i in
1110 its internal vec representation.
1111 */
1112
matrix_get_element(const gretl_matrix * M,int i,int * err)1113 double matrix_get_element (const gretl_matrix *M, int i, int *err)
1114 {
1115 double x = NADBL;
1116
1117 if (M == NULL) {
1118 *err = E_DATA;
1119 } else if (i < 0 || i >= M->rows * M->cols) {
1120 gretl_errmsg_sprintf(_("Index value %d is out of bounds"), i+1);
1121 *err = E_INVARG;
1122 } else {
1123 x = M->val[i];
1124 }
1125
1126 return x;
1127 }
1128
1129 /* Copy a contiguous chunk of data out of @M, the offset and
1130 extent of which are given by @spec.
1131 */
1132
matrix_get_chunk(const gretl_matrix * M,matrix_subspec * spec,int * err)1133 gretl_matrix *matrix_get_chunk (const gretl_matrix *M,
1134 matrix_subspec *spec,
1135 int *err)
1136 {
1137 int offset = spec->lsel.range[0];
1138 int nelem = spec->lsel.range[1];
1139 gretl_matrix *ret;
1140 int rows, cols;
1141
1142 if (offset < 0) {
1143 fprintf(stderr, "matrix_get_chunk: offset = %d\n", offset);
1144 *err = E_DATA;
1145 return NULL;
1146 }
1147
1148 if (M->cols > 1 && M->rows > 1) {
1149 cols = contig_cols(spec, M);
1150 rows = nelem / cols;
1151 } else if (M->rows == 1) {
1152 rows = 1;
1153 cols = nelem;
1154 } else {
1155 rows = nelem;
1156 cols = 1;
1157 }
1158
1159 ret = gretl_matching_matrix_new(rows, cols, M);
1160
1161 if (ret == NULL) {
1162 *err = E_ALLOC;
1163 } else {
1164 size_t sz;
1165
1166 if (M->is_complex) {
1167 sz = nelem * sizeof *M->z;
1168 memcpy(ret->z, M->z + offset, sz);
1169 } else {
1170 sz = nelem * sizeof *M->val;
1171 memcpy(ret->val, M->val + offset, sz);
1172 }
1173 if (M->rows > 1 && rows == M->rows && offset % rows == 0 &&
1174 gretl_matrix_is_dated(M)) {
1175 gretl_matrix_transcribe_obs_info(ret, M);
1176 }
1177 }
1178
1179 return ret;
1180 }
1181
1182 /* Handle the case where we got a single string as argument
1183 to colnames() or rownames(), for a matrix with more than
1184 one column or row: construct specific names by appending
1185 a column or row index.
1186 */
1187
expand_names(const char * s,int n,int * err)1188 static char **expand_names (const char *s, int n, int *err)
1189 {
1190 char **S = NULL;
1191
1192 if (!gretl_is_ascii(s)) {
1193 *err = E_INVARG;
1194 return NULL;
1195 }
1196
1197 S = strings_array_new(n);
1198
1199 if (S == NULL) {
1200 *err = E_ALLOC;
1201 } else {
1202 char tmp[12];
1203 int i, m;
1204
1205 for (i=0; i<n && !*err; i++) {
1206 sprintf(tmp, "%d", i+1);
1207 m = strlen(tmp);
1208 sprintf(tmp, "%.*s%d", 9-m, s, i+1);
1209 S[i] = gretl_strdup(tmp);
1210 if (S[i] == NULL) {
1211 *err = E_ALLOC;
1212 }
1213 }
1214 if (*err) {
1215 strings_array_free(S, n);
1216 S = NULL;
1217 }
1218 }
1219
1220 return S;
1221 }
1222
umatrix_set_names_from_string(gretl_matrix * M,const char * s,int byrow)1223 int umatrix_set_names_from_string (gretl_matrix *M,
1224 const char *s,
1225 int byrow)
1226 {
1227 int n, err = 0;
1228
1229 n = (byrow)? M->rows : M->cols;
1230
1231 if (s == NULL || *s == '\0') {
1232 if (byrow) {
1233 gretl_matrix_set_rownames(M, NULL);
1234 } else {
1235 gretl_matrix_set_colnames(M, NULL);
1236 }
1237 } else {
1238 char **S;
1239 int ns;
1240
1241 S = gretl_string_split(s, &ns, " \n\t");
1242
1243 if (S == NULL) {
1244 err = E_ALLOC;
1245 } else if (ns == 1 && n > 1) {
1246 strings_array_free(S, ns);
1247 S = expand_names(s, n, &err);
1248 } else if (ns != n) {
1249 err = E_NONCONF;
1250 strings_array_free(S, ns);
1251 }
1252
1253 if (!err) {
1254 if (byrow) {
1255 gretl_matrix_set_rownames(M, S);
1256 } else {
1257 gretl_matrix_set_colnames(M, S);
1258 }
1259 }
1260 }
1261
1262 return err;
1263 }
1264
umatrix_set_names_from_array(gretl_matrix * M,void * data,int byrow)1265 int umatrix_set_names_from_array (gretl_matrix *M,
1266 void *data,
1267 int byrow)
1268 {
1269 gretl_array *A = data;
1270 int n, err = 0;
1271
1272 n = (byrow)? M->rows : M->cols;
1273
1274 if (A == NULL || gretl_array_get_length(A) == 0) {
1275 if (byrow) {
1276 gretl_matrix_set_rownames(M, NULL);
1277 } else {
1278 gretl_matrix_set_colnames(M, NULL);
1279 }
1280 } else {
1281 char **AS;
1282 int i, ns;
1283
1284 AS = gretl_array_get_strings(A, &ns);
1285
1286 if (ns != n) {
1287 err = E_NONCONF;
1288 } else {
1289 for (i=0; i<ns && !err; i++) {
1290 if (AS[i] == NULL || AS[i][0] == '\0') {
1291 fprintf(stderr, "Missing string in colnames/rownames\n");
1292 err = E_INVARG;
1293 }
1294 }
1295 }
1296
1297 if (!err) {
1298 char **S = strings_array_dup(AS, ns);
1299
1300 if (S == NULL) {
1301 err = E_ALLOC;
1302 } else if (byrow) {
1303 gretl_matrix_set_rownames(M, S);
1304 } else {
1305 gretl_matrix_set_colnames(M, S);
1306 }
1307 }
1308 }
1309
1310 return err;
1311 }
1312
umatrix_set_names_from_list(gretl_matrix * M,const int * list,const DATASET * dset,int byrow)1313 int umatrix_set_names_from_list (gretl_matrix *M,
1314 const int *list,
1315 const DATASET *dset,
1316 int byrow)
1317 {
1318 int i, n, err = 0;
1319
1320 n = (byrow)? M->rows : M->cols;
1321
1322 if (list == NULL || list[0] == 0) {
1323 if (byrow) {
1324 gretl_matrix_set_rownames(M, NULL);
1325 } else {
1326 gretl_matrix_set_colnames(M, NULL);
1327 }
1328 } else if (list[0] != n) {
1329 err = E_NONCONF;
1330 } else {
1331 char **S = strings_array_new(n);
1332
1333 if (S == NULL) {
1334 err = E_ALLOC;
1335 }
1336
1337 for (i=0; i<n && !err; i++) {
1338 S[i] = gretl_strndup(dset->varname[list[i+1]], 12);
1339 if (S[i] == NULL) {
1340 err = E_ALLOC;
1341 }
1342 }
1343
1344 if (err) {
1345 strings_array_free(S, n);
1346 } else if (byrow) {
1347 gretl_matrix_set_rownames(M, S);
1348 } else {
1349 gretl_matrix_set_colnames(M, S);
1350 }
1351 }
1352
1353 return err;
1354 }
1355
user_matrix_get_column_name(const gretl_matrix * M,int col,int * err)1356 char *user_matrix_get_column_name (const gretl_matrix *M, int col,
1357 int *err)
1358 {
1359 char *ret = NULL;
1360
1361 if (M == NULL || col < 1 || col > M->cols) {
1362 *err = E_DATA;
1363 } else {
1364 const char **S = gretl_matrix_get_colnames(M);
1365
1366 if (S == NULL) {
1367 ret = gretl_strdup("");
1368 } else {
1369 ret = gretl_strdup(S[col-1]);
1370 }
1371 if (ret == NULL) {
1372 *err = E_ALLOC;
1373 }
1374 }
1375
1376 return ret;
1377 }
1378
user_matrix_get_row_name(const gretl_matrix * M,int row,int * err)1379 char *user_matrix_get_row_name (const gretl_matrix *M, int row,
1380 int *err)
1381 {
1382 char *ret = NULL;
1383
1384 if (M == NULL || row < 1 || row > M->rows) {
1385 *err = E_DATA;
1386 } else {
1387 const char **S = gretl_matrix_get_rownames(M);
1388
1389 if (S == NULL) {
1390 ret = gretl_strdup("");
1391 } else {
1392 ret = gretl_strdup(S[row-1]);
1393 }
1394 if (ret == NULL) {
1395 *err = E_ALLOC;
1396 }
1397 }
1398
1399 return ret;
1400 }
1401
1402 double
user_matrix_get_determinant(gretl_matrix * m,int tmpmat,int f,int * err)1403 user_matrix_get_determinant (gretl_matrix *m, int tmpmat,
1404 int f, int *err)
1405 {
1406 gretl_matrix *R = NULL;
1407 double d = NADBL;
1408
1409 if (gretl_is_null_matrix(m)) {
1410 return d;
1411 } else if (tmpmat) {
1412 /* it's OK to overwrite @m */
1413 R = m;
1414 } else {
1415 /* @m should not be over-written! */
1416 R = gretl_matrix_copy(m);
1417 }
1418
1419 if (R != NULL) {
1420 if (f == F_LDET) {
1421 d = gretl_matrix_log_determinant(R, err);
1422 } else {
1423 d = gretl_matrix_determinant(R, err);
1424 }
1425 if (R != m) {
1426 gretl_matrix_free(R);
1427 }
1428 }
1429
1430 return d;
1431 }
1432
maybe_replace_content(gretl_matrix * targ,gretl_matrix * tmp,int err)1433 static void maybe_replace_content (gretl_matrix *targ,
1434 gretl_matrix *tmp,
1435 int err)
1436 {
1437 if (!err) {
1438 gretl_matrix_replace_content(targ, tmp);
1439 }
1440 gretl_matrix_free(tmp);
1441 }
1442
matrix_invert_in_place(gretl_matrix * m)1443 int matrix_invert_in_place (gretl_matrix *m)
1444 {
1445 gretl_matrix *tmp = gretl_matrix_copy(m);
1446 int err = 0;
1447
1448 if (tmp == NULL) {
1449 err = E_ALLOC;
1450 } else {
1451 err = gretl_invert_matrix(tmp);
1452 maybe_replace_content(m, tmp, err);
1453 }
1454
1455 return err;
1456 }
1457
matrix_cholesky_in_place(gretl_matrix * m)1458 int matrix_cholesky_in_place (gretl_matrix *m)
1459 {
1460 gretl_matrix *tmp = gretl_matrix_copy(m);
1461 int err = 0;
1462
1463 if (tmp == NULL) {
1464 err = E_ALLOC;
1465 } else {
1466 err = gretl_matrix_cholesky_decomp(tmp);
1467 maybe_replace_content(m, tmp, err);
1468 }
1469
1470 return err;
1471 }
1472
matrix_transpose_in_place(gretl_matrix * m)1473 int matrix_transpose_in_place (gretl_matrix *m)
1474 {
1475 gretl_matrix *tmp = gretl_matrix_copy_transpose(m);
1476 int err = 0;
1477
1478 if (tmp == NULL) {
1479 err = E_ALLOC;
1480 } else {
1481 maybe_replace_content(m, tmp, 0);
1482 }
1483
1484 return err;
1485 }
1486
matrix_XTX_in_place(gretl_matrix * m)1487 int matrix_XTX_in_place (gretl_matrix *m)
1488 {
1489 gretl_matrix *tmp = gretl_matrix_alloc(m->cols, m->cols);
1490 int err;
1491
1492 if (tmp == NULL) {
1493 err = E_ALLOC;
1494 } else {
1495 err = gretl_matrix_multiply_mod(m, GRETL_MOD_TRANSPOSE,
1496 m, GRETL_MOD_NONE,
1497 tmp, GRETL_MOD_NONE);
1498 maybe_replace_content(m, tmp, err);
1499 }
1500
1501 return err;
1502 }
1503
user_matrix_vec(const gretl_matrix * m,int * err)1504 gretl_matrix *user_matrix_vec (const gretl_matrix *m, int *err)
1505 {
1506 gretl_matrix *R = NULL;
1507
1508 if (gretl_is_null_matrix(m)) {
1509 R = gretl_null_matrix_new();
1510 } else {
1511 R = gretl_matching_matrix_new(m->rows * m->cols, 1, m);
1512 if (R != NULL) {
1513 gretl_matrix_vectorize(R, m);
1514 }
1515 }
1516
1517 if (R == NULL) {
1518 *err = E_ALLOC;
1519 }
1520
1521 return R;
1522 }
1523
user_matrix_vech(const gretl_matrix * m,int * err)1524 gretl_matrix *user_matrix_vech (const gretl_matrix *m, int *err)
1525 {
1526 gretl_matrix *R = NULL;
1527
1528 if (gretl_is_null_matrix(m)) {
1529 R = gretl_null_matrix_new();
1530 } else if (m->rows != m->cols) {
1531 *err = E_NONCONF;
1532 } else {
1533 int n = m->rows;
1534 int k = n * (n + 1) / 2;
1535
1536 R = gretl_matching_matrix_new(k, 1, m);
1537 if (R != NULL) {
1538 *err = gretl_matrix_vectorize_h(R, m);
1539 }
1540 }
1541
1542 if (R == NULL && !*err) {
1543 *err = E_ALLOC;
1544 }
1545
1546 return R;
1547 }
1548
user_matrix_unvech(const gretl_matrix * m,int * err)1549 gretl_matrix *user_matrix_unvech (const gretl_matrix *m, int *err)
1550 {
1551 gretl_matrix *R = NULL;
1552
1553 if (gretl_is_null_matrix(m)) {
1554 R = gretl_null_matrix_new();
1555 } else if (m->cols != 1) {
1556 *err = E_NONCONF;
1557 } else {
1558 int n = (int) ((sqrt(1.0 + 8.0 * m->rows) - 1.0) / 2.0);
1559
1560 R = gretl_matching_matrix_new(n, n, m);
1561 if (R != NULL) {
1562 *err = gretl_matrix_unvectorize_h(R, m);
1563 }
1564 }
1565
1566 if (R == NULL && !*err) {
1567 *err = E_ALLOC;
1568 }
1569
1570 return R;
1571 }
1572
1573 static int
real_user_matrix_QR_decomp(const gretl_matrix * m,gretl_matrix ** Q,gretl_matrix ** R)1574 real_user_matrix_QR_decomp (const gretl_matrix *m, gretl_matrix **Q,
1575 gretl_matrix **R)
1576 {
1577 int mc = gretl_matrix_cols(m);
1578 int err = 0;
1579
1580 *Q = gretl_matrix_copy(m);
1581
1582 if (*Q == NULL) {
1583 err = E_ALLOC;
1584 } else if (R != NULL) {
1585 *R = gretl_matrix_alloc(mc, mc);
1586 if (*R == NULL) {
1587 err = E_ALLOC;
1588 }
1589 }
1590
1591 if (!err) {
1592 err = gretl_matrix_QR_decomp(*Q, (R == NULL)? NULL : *R);
1593 }
1594
1595 if (err) {
1596 gretl_errmsg_set(_("Matrix decomposition failed"));
1597 gretl_matrix_free(*Q);
1598 *Q = NULL;
1599 if (R != NULL) {
1600 gretl_matrix_free(*R);
1601 *R = NULL;
1602 }
1603 }
1604
1605 return err;
1606 }
1607
user_matrix_QR_decomp(const gretl_matrix * m,gretl_matrix * R,int * err)1608 gretl_matrix *user_matrix_QR_decomp (const gretl_matrix *m,
1609 gretl_matrix *R,
1610 int *err)
1611 {
1612 gretl_matrix *Q = NULL;
1613 gretl_matrix *Rtmp = NULL;
1614 gretl_matrix **pR;
1615
1616 if (gretl_is_null_matrix(m)) {
1617 *err = E_DATA;
1618 return NULL;
1619 }
1620
1621 pR = (R != NULL)? &Rtmp : NULL;
1622
1623 *err = real_user_matrix_QR_decomp(m, &Q, pR);
1624 if (Rtmp != NULL) {
1625 maybe_replace_content(R, Rtmp, *err);
1626 }
1627
1628 return Q;
1629 }
1630
user_matrix_SVD(const gretl_matrix * m,gretl_matrix * U,gretl_matrix * V,int * err)1631 gretl_matrix *user_matrix_SVD (const gretl_matrix *m,
1632 gretl_matrix *U,
1633 gretl_matrix *V,
1634 int *err)
1635 {
1636 gretl_matrix *S = NULL;
1637 gretl_matrix *Utmp = NULL;
1638 gretl_matrix *Vtmp = NULL;
1639 gretl_matrix **pU, **pV;
1640
1641 if (gretl_is_null_matrix(m)) {
1642 *err = E_DATA;
1643 return NULL;
1644 }
1645
1646 pU = (U != NULL)? &Utmp : NULL;
1647 pV = (V != NULL)? &Vtmp : NULL;
1648
1649 if (m->is_complex) {
1650 *err = gretl_cmatrix_SVD(m, pU, &S, pV, 0);
1651 } else {
1652 *err = gretl_matrix_SVD(m, pU, &S, pV, 0);
1653 }
1654 if (!*err) {
1655 if (Utmp != NULL) {
1656 maybe_replace_content(U, Utmp, *err);
1657 }
1658 if (Vtmp != NULL) {
1659 maybe_replace_content(V, Vtmp, *err);
1660 }
1661 }
1662
1663 return S;
1664 }
1665
null_OLS(const gretl_matrix * Y,gretl_matrix * U,gretl_matrix * V,int * err)1666 static gretl_matrix *null_OLS (const gretl_matrix *Y,
1667 gretl_matrix *U,
1668 gretl_matrix *V,
1669 int *err)
1670 {
1671 gretl_matrix *B = NULL;
1672
1673 if (U != NULL) {
1674 /* residuals = Y */
1675 gretl_matrix *Ycpy;
1676
1677 Ycpy = gretl_matrix_copy(Y);
1678 if (Ycpy == NULL) {
1679 *err = E_ALLOC;
1680 } else {
1681 gretl_matrix_replace_content(U, Ycpy);
1682 gretl_matrix_free(Ycpy);
1683 }
1684 }
1685
1686 if (!*err && V != NULL) {
1687 /* variance is empty */
1688 gretl_matrix *nullV;
1689
1690 nullV = gretl_null_matrix_new();
1691 gretl_matrix_replace_content(V, nullV);
1692 gretl_matrix_free(nullV);
1693 }
1694
1695 if (!*err) {
1696 B = gretl_matrix_alloc(0, Y->cols);
1697 }
1698
1699 return B;
1700 }
1701
user_matrix_ols(const gretl_matrix * Y,const gretl_matrix * X,gretl_matrix * U,gretl_matrix * V,gretlopt opt,int * err)1702 gretl_matrix *user_matrix_ols (const gretl_matrix *Y,
1703 const gretl_matrix *X,
1704 gretl_matrix *U,
1705 gretl_matrix *V,
1706 gretlopt opt,
1707 int *err)
1708 {
1709 gretl_matrix *B = NULL;
1710 int g, k, T;
1711
1712 if (gretl_is_null_matrix(Y) || X == NULL) {
1713 *err = E_DATA;
1714 return NULL;
1715 }
1716
1717 if (gretl_is_complex(Y) || gretl_is_complex(X) ||
1718 gretl_is_complex(U) || gretl_is_complex(V)) {
1719 fprintf(stderr, "E_CMPLX in user_matrix_ols\n");
1720 *err = E_CMPLX;
1721 return NULL;
1722 }
1723
1724 T = Y->rows;
1725 k = X->cols;
1726 g = Y->cols;
1727
1728 if (X->rows != T) {
1729 *err = E_NONCONF;
1730 return NULL;
1731 }
1732
1733 if (X->cols == 0) {
1734 /* handle the null model case */
1735 return null_OLS(Y, U, V, err);
1736 }
1737
1738 if (g > 1 && (opt & OPT_M)) {
1739 /* multiple precision: we accept only one y var */
1740 *err = E_DATA;
1741 return NULL;
1742 }
1743
1744 if (U != NULL) {
1745 *err = gretl_matrix_realloc(U, T, g);
1746 if (*err) {
1747 return NULL;
1748 }
1749 }
1750
1751 if (V != NULL && g == 1) {
1752 /* single regressand case */
1753 *err = gretl_matrix_realloc(V, k, k);
1754 }
1755
1756 if (!*err) {
1757 B = gretl_matrix_alloc(k, g);
1758 if (B == NULL) {
1759 *err = E_ALLOC;
1760 }
1761 }
1762
1763 if (!*err) {
1764 if (g == 1) {
1765 /* single regressand */
1766 double s2 = 0;
1767 double *ps2 = (V != NULL)? &s2 : NULL;
1768
1769 if (opt & OPT_M) {
1770 /* use multiple precision */
1771 *err = gretl_matrix_mp_ols(Y, X, B, V, U, ps2);
1772 } else {
1773 *err = gretl_matrix_ols(Y, X, B, V, U, ps2);
1774 }
1775 } else {
1776 /* multiple regressands: @V will actually be (X'X)^{-1} */
1777 gretl_matrix *Vtmp = NULL;
1778 gretl_matrix **Vp = (V != NULL)? &Vtmp : NULL;
1779
1780 *err = gretl_matrix_multi_ols(Y, X, B, U, Vp);
1781 if (Vtmp != NULL) {
1782 maybe_replace_content(V, Vtmp, *err);
1783 }
1784 }
1785 }
1786
1787 if (*err) {
1788 gretl_matrix_free(B);
1789 B = NULL;
1790 }
1791
1792 return B;
1793 }
1794
user_matrix_rls(const gretl_matrix * Y,const gretl_matrix * X,const gretl_matrix * R,const gretl_matrix * Q,gretl_matrix * U,gretl_matrix * V,int * err)1795 gretl_matrix *user_matrix_rls (const gretl_matrix *Y,
1796 const gretl_matrix *X,
1797 const gretl_matrix *R,
1798 const gretl_matrix *Q,
1799 gretl_matrix *U,
1800 gretl_matrix *V,
1801 int *err)
1802 {
1803 gretl_matrix *B = NULL;
1804 int g, k, T;
1805
1806 if (gretl_is_null_matrix(Y) || gretl_is_null_matrix(X)) {
1807 *err = E_DATA;
1808 return NULL;
1809 }
1810
1811 if (gretl_is_complex(Y) || gretl_is_complex(X) ||
1812 gretl_is_complex(R) || gretl_is_complex(Q)) {
1813 fprintf(stderr, "E_CMPLX in user_matrix_rls\n");
1814 *err = E_CMPLX;
1815 return NULL;
1816 }
1817
1818 T = Y->rows;
1819 k = X->cols;
1820 g = Y->cols;
1821
1822 if (X->rows != T) {
1823 *err = E_NONCONF;
1824 return NULL;
1825 }
1826
1827 if (U != NULL) {
1828 *err = gretl_matrix_realloc(U, T, g);
1829 if (*err) {
1830 return NULL;
1831 }
1832 }
1833
1834 if (!*err) {
1835 B = gretl_matrix_alloc(k, g);
1836 if (B == NULL) {
1837 *err = E_ALLOC;
1838 }
1839 }
1840
1841 if (!*err) {
1842 /* note: &V will actually be M (X'X)^{-1} M' */
1843 gretl_matrix *Vtmp = NULL;
1844 gretl_matrix **Vp = (V != NULL)? &Vtmp : NULL;
1845
1846 *err = gretl_matrix_restricted_multi_ols(Y, X, R, Q, B,
1847 U, Vp);
1848 if (Vtmp != NULL) {
1849 maybe_replace_content(V, Vtmp, *err);
1850 }
1851 }
1852
1853 if (*err) {
1854 gretl_matrix_free(B);
1855 B = NULL;
1856 }
1857
1858 return B;
1859 }
1860
user_matrix_GHK(const gretl_matrix * C,const gretl_matrix * A,const gretl_matrix * B,const gretl_matrix * U,gretl_matrix * dP,int * err)1861 gretl_matrix *user_matrix_GHK (const gretl_matrix *C,
1862 const gretl_matrix *A,
1863 const gretl_matrix *B,
1864 const gretl_matrix *U,
1865 gretl_matrix *dP,
1866 int *err)
1867 {
1868 gretl_matrix *P = NULL;
1869 int n, m, npar;
1870
1871 if (gretl_is_null_matrix(A) || gretl_is_null_matrix(C)) {
1872 *err = E_DATA;
1873 return NULL;
1874 }
1875
1876 n = A->rows;
1877 m = C->rows;
1878 npar = m + m + m*(m+1)/2;
1879
1880 if (dP != NULL) {
1881 *err = gretl_matrix_realloc(dP, n, npar);
1882 }
1883
1884 if (!*err) {
1885 P = gretl_GHK2(C, A, B, U, dP, err);
1886 }
1887
1888 return P;
1889 }
1890
user_matrix_eigensym(const gretl_matrix * m,gretl_matrix * R,int * err)1891 gretl_matrix *user_matrix_eigensym (const gretl_matrix *m,
1892 gretl_matrix *R,
1893 int *err)
1894 {
1895 gretl_matrix *C = NULL;
1896 gretl_matrix *E = NULL;
1897 int vecs = 0;
1898
1899 if (gretl_is_null_matrix(m)) {
1900 *err = E_DATA;
1901 return NULL;
1902 }
1903
1904 if (gretl_matrix_na_check(m)) {
1905 *err = E_NAN;
1906 return NULL;
1907 }
1908
1909 if (R != NULL) {
1910 /* computing (right) eigenvectors */
1911 vecs = 1;
1912 }
1913
1914 /* @m would be destroyed by this operation */
1915 C = gretl_matrix_copy(m);
1916 if (C == NULL) {
1917 *err = E_ALLOC;
1918 }
1919
1920 if (!*err) {
1921 if (m->is_complex) {
1922 E = gretl_zheev(C, vecs, err);
1923 } else {
1924 E = gretl_symmetric_matrix_eigenvals(C, vecs, err);
1925 }
1926 }
1927
1928 if (!*err && vecs) {
1929 maybe_replace_content(R, C, 0);
1930 }
1931
1932 if (!vecs) {
1933 gretl_matrix_free(C);
1934 }
1935
1936 return E;
1937 }
1938
user_gensymm_eigenvals(const gretl_matrix * A,const gretl_matrix * B,gretl_matrix * V,int * err)1939 gretl_matrix *user_gensymm_eigenvals (const gretl_matrix *A,
1940 const gretl_matrix *B,
1941 gretl_matrix *V,
1942 int *err)
1943 {
1944 gretl_matrix *E = NULL;
1945
1946 if (gretl_is_null_matrix(A) || gretl_is_null_matrix(B)) {
1947 *err = E_DATA;
1948 return NULL;
1949 }
1950
1951 if (gretl_matrix_na_check(A) || gretl_matrix_na_check(B)) {
1952 *err = E_NAN;
1953 return NULL;
1954 }
1955
1956 if (V != NULL) {
1957 *err = gretl_matrix_realloc(V, B->cols, A->rows);
1958 }
1959
1960 if (!*err) {
1961 E = gretl_gensymm_eigenvals(A, B, V, err);
1962 }
1963
1964 return E;
1965 }
1966
gretl_matrix_set_part(gretl_matrix * targ,const gretl_matrix * src,double x,SelType sel)1967 int gretl_matrix_set_part (gretl_matrix *targ,
1968 const gretl_matrix *src,
1969 double x, SelType sel)
1970 {
1971 int err = 0;
1972
1973 if (sel == SEL_DIAG) {
1974 if (targ->is_complex) {
1975 err = gretl_cmatrix_set_diagonal(targ, src, x);
1976 } else {
1977 err = gretl_matrix_set_diagonal(targ, src, x);
1978 }
1979 } else if (sel == SEL_LOWER || sel == SEL_UPPER) {
1980 int upper = (sel == SEL_UPPER);
1981
1982 if (targ->is_complex) {
1983 err = gretl_cmatrix_set_triangle(targ, src, x, upper);
1984 } else {
1985 err = gretl_matrix_set_triangle(targ, src, x, upper);
1986 }
1987 } else if (sel == SEL_REAL || sel == SEL_IMAG) {
1988 if (targ->is_complex) {
1989 int im = (sel == SEL_IMAG);
1990
1991 err = gretl_cmatrix_set_part(targ, src, x, im);
1992 } else {
1993 err = E_TYPES;
1994 }
1995 }
1996
1997 return err;
1998 }
1999