1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: matrix.c *
8 * *
9 * Routines for computations with square matrices. *
10 * Matrices are stored as double array (i.e.: double *matrix). *
11 * The rows of the matrix have to be stored consecutively in this array, *
12 * i.e. the entry with index [i,j] can be entered via (i*dim+j), where *
13 * dim is the dimension of the dim x dim - matrix. *
14 * *
15 *****************************************************************************
16 * *
17 * Routines are mainly adapted from the GSL (GNU Scientifiy Library) *
18 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough *
19 * *
20 * adapted by Wolfgang Hoermann and Josef Leydold *
21 * Department of Statistics and Mathematics, WU Wien, Austria *
22 * *
23 * This program is free software; you can redistribute it and/or modify *
24 * it under the terms of the GNU General Public License as published by *
25 * the Free Software Foundation; either version 2 of the License, or *
26 * (at your option) any later version. *
27 * *
28 * This program is distributed in the hope that it will be useful, *
29 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
30 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
31 * GNU General Public License for more details. *
32 * *
33 * You should have received a copy of the GNU General Public License *
34 * along with this program; if not, write to the *
35 * Free Software Foundation, Inc., *
36 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
37 * *
38 *****************************************************************************/
39
40 /*--------------------------------------------------------------------------*/
41
42 #include <unur_source.h>
43 #include "matrix_source.h"
44
45 /*---------------------------------------------------------------------------*/
46
47 static int _unur_matrix_swap_rows (int dim, double *A, int i, int j);
48 /* Swap rows i and j in square matrix A of dimension dim. */
49
50 static int _unur_matrix_permutation_swap (int dim, int *p, int i, int j);
51 /* Swap entries i and j of the integer array p. */
52
53 static int _unur_matrix_LU_decomp (int dim, double *A, int *P, int *signum);
54 /* Factorise a general dim x dim matrix A into P A = L U. */
55
56 static int _unur_matrix_backsubstitution_dtrsv(int dim, double *LU, double *X);
57 /* Backsubstitution used for inversion alg. _unur_matrix_LU_invert(). */
58
59 static int _unur_matrix_forwardsubstitution_dtrsv(int dim, double *LU, double *X);
60 /* Forwardsubstitution used for inversion alg. _unur_matrix_LU_invert(). */
61
62 static int _unur_matrix_LU_invert (int dim, double *LU, int *p, double *inverse);
63 /* Compute inverse of matrix with given LU decomposition. */
64
65 /*---------------------------------------------------------------------------*/
66
67 int
_unur_matrix_transform_diagonal(int dim,const double * M,const double * D,double * res)68 _unur_matrix_transform_diagonal (int dim, const double *M, const double *D, double *res)
69 /*----------------------------------------------------------------------*/
70 /* Computes the transformation M^t . D . M of diagonal matrix D */
71 /* and stores it in matrix 'res'. */
72 /* */
73 /* input: */
74 /* dim ... number of columns and rows of m and diag */
75 /* M ... square matrix */
76 /* D ... diagonal entries of a diagonal matrix */
77 /* */
78 /* output: */
79 /* res ... square matrix to store result */
80 /* */
81 /* return: */
82 /* UNUR_SUCCESS on success */
83 /* error code otherwise */
84 /*----------------------------------------------------------------------*/
85 {
86 #define idx(a,b) ((a)*dim+(b))
87
88 int i,j,k;
89 double sum;
90
91 /* check arguments */
92 CHECK_NULL(M,UNUR_ERR_NULL);
93 CHECK_NULL(D,UNUR_ERR_NULL);
94 CHECK_NULL(res,UNUR_ERR_NULL);
95
96 for (i=0; i<dim; i++)
97 for (j=0; j<dim; j++) {
98 for (sum=0., k=0; k<dim; k++)
99 sum += D[k] * M[idx(k,i)] * M[idx(k,j)];
100 res[idx(i,j)] = sum;
101 }
102
103 return UNUR_SUCCESS;
104
105 #undef idx
106 } /* end of _unur_matrix_transform_diagonal() */
107
108 /*---------------------------------------------------------------------------*/
109
110 int
_unur_matrix_swap_rows(int dim,double * A,int i,int j)111 _unur_matrix_swap_rows (int dim, double *A, int i, int j)
112 /*----------------------------------------------------------------------*/
113 /* Swap rows i and j in square matrix A of dimension dim. */
114 /* The resulting matrix is again stored in A. */
115 /* */
116 /* input: */
117 /* dim ... number of columns and rows of */
118 /* A ... dim x dim -matrix */
119 /* i,j ... row numbers that should be swapped */
120 /* */
121 /* output: */
122 /* the given matrix A contains the new matrix with swapped rows */
123 /* */
124 /* return: */
125 /* UNUR_SUCCESS on success */
126 /* error code otherwise */
127 /*----------------------------------------------------------------------*/
128 {
129 /* check arguments */
130 CHECK_NULL(A,UNUR_ERR_NULL);
131
132 if (i != j)
133 {
134 double *row1 = A + i * dim;
135 double *row2 = A + j * dim;
136
137 int k;
138
139 for (k = 0; k < dim; k++)
140 {
141 double tmp = row1[k] ;
142 row1[k] = row2[k] ;
143 row2[k] = tmp ;
144 }
145 }
146
147 return UNUR_SUCCESS;
148 } /* end of _unur_matrix_swap_rows() */
149
150 /*---------------------------------------------------------------------------*/
151
152 int
_unur_matrix_permutation_swap(int dim ATTRIBUTE__UNUSED,int * p,int i,int j)153 _unur_matrix_permutation_swap (int dim ATTRIBUTE__UNUSED, int *p, int i, int j)
154 /*----------------------------------------------------------------------*/
155 /* Swap entries i and j of the integer array p. */
156 /* */
157 /* input: */
158 /* dim ... number of elements */
159 /* p ... permutation vector of length dim */
160 /* (contains all integers between 0 and dim-1) */
161 /* i,j ... indexes that should be swapped */
162 /* */
163 /* output: */
164 /* permutation p with swapped elements */
165 /*----------------------------------------------------------------------*/
166 {
167 if (i != j)
168 {
169 int tmp = p[i];
170 p[i] = p[j];
171 p[j] = tmp;
172 }
173
174 return UNUR_SUCCESS;
175 } /* end of _unur_matrix_permutation_swap() */
176
177 /*---------------------------------------------------------------------------*/
178
179 int
_unur_matrix_LU_decomp(int dim,double * A,int * p,int * signum)180 _unur_matrix_LU_decomp (int dim, double *A, int *p, int *signum)
181 /*----------------------------------------------------------------------*/
182 /* Factorise a general dim x dim matrix A into, */
183 /* */
184 /* P A = L U */
185 /* */
186 /* where P is a permutation matrix, L is unit lower triangular and U */
187 /* is upper triangular. */
188 /* */
189 /* L is stored in the strict lower triangular part of the input */
190 /* matrix. The diagonal elements of L are unity and are not stored. */
191 /* */
192 /* U is stored in the diagonal and upper triangular part of the */
193 /* input matrix. */
194 /* */
195 /* P is stored in the permutation p. Column j of P is column k of the */
196 /* identity matrix, where k = permutation->data[j] */
197 /* */
198 /* signum gives the sign of the permutation, (-1)^n, where n is the */
199 /* number of interchanges in the permutation. */
200 /* */
201 /* See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss */
202 /* Elimination with Partial Pivoting). */
203 /* */
204 /* input: */
205 /* dim ... number of columns and rows of */
206 /* A ... dim x dim matrix */
207 /* p ... pointer to array where permutation should be stored */
208 /* signum ... pointer where sign of permutation should be stored */
209 /* */
210 /* output: */
211 /* A ... L (strict lower triangular part) and U (upper part) */
212 /* p ... permutation vector of length dim */
213 /* (contains all integers between 0 and dim-1) */
214 /* signum ... sign of the permutation (1 or -1) */
215 /* */
216 /* return: */
217 /* UNUR_SUCCESS on success */
218 /* error code otherwise */
219 /*----------------------------------------------------------------------*/
220 {
221 #define idx(a,b) ((a)*dim+(b))
222 int i, j, k;
223
224 /* check arguments */
225 CHECK_NULL(A,UNUR_ERR_NULL);
226 CHECK_NULL(p,UNUR_ERR_NULL);
227 CHECK_NULL(signum,UNUR_ERR_NULL);
228
229 /* initialize permutation vector */
230 *signum = 1;
231 for(i=0;i<dim;i++) p[i]=i;
232
233 for (j = 0; j < dim - 1; j++){
234 /* Find maximum in the j-th column */
235
236 double ajj, max = fabs (A[idx(j,j)]);
237 int i_pivot = j;
238 for (i = j + 1; i < dim; i++){
239 double aij = fabs (A[idx(i,j)]);
240 if (aij > max){
241 max = aij;
242 i_pivot = i;
243 }
244 }
245 if (i_pivot != j){
246 _unur_matrix_swap_rows (dim, A, j, i_pivot);
247 _unur_matrix_permutation_swap (dim, p, j, i_pivot);
248 *signum = -(*signum);
249 }
250
251 ajj = A[idx(j,j)];
252
253 if ( !_unur_iszero(ajj) ){
254 for (i = j + 1; i < dim; i++){
255 double aij = A[idx(i,j)] / ajj;
256 A[idx(i,j)] = aij;
257 for (k = j + 1; k < dim; k++){
258 double aik = A[idx(i,k)];
259 double ajk = A[idx(j,k)];
260 A[idx(i,k)]= aik - aij * ajk;
261 }
262 }
263 }
264 }
265
266 return UNUR_SUCCESS;
267
268 #undef idx
269 } /* end of _unur_matrix_LU_decomp() */
270
271 /*---------------------------------------------------------------------------*/
272
273 int
_unur_matrix_backsubstitution_dtrsv(int dim,double * LU,double * X)274 _unur_matrix_backsubstitution_dtrsv(int dim, double *LU, double *X)
275 /*----------------------------------------------------------------------*/
276 /* Backsubstitution used for inversion alg. _unur_matrix_LU_invert() */
277 /* */
278 /* input: */
279 /* dim ... number of columns and rows of */
280 /* LU ... dim x dim -matrix that contains LU decomposition of matrix */
281 /* X ... vector */
282 /* */
283 /* output: */
284 /* X */
285 /* */
286 /* return: */
287 /* UNUR_SUCCESS on success */
288 /* error code otherwise */
289 /*----------------------------------------------------------------------*/
290 {
291 #define idx(a,b) ((a)*dim+(b))
292
293 int ix,jx,i,j;
294
295 /* check arguments */
296 CHECK_NULL(LU,UNUR_ERR_NULL);
297 CHECK_NULL(X,UNUR_ERR_NULL);
298
299 /* backsubstitution */
300 ix = (dim - 1);
301
302 X[ix] = X[ix] / LU[idx(ix,ix)];
303
304 ix--;
305 for (i = dim - 1; i > 0 && i--;) {
306 double tmp = X[ix];
307 jx = ix + 1;
308 for (j = i + 1; j < dim; j++) {
309 tmp -= LU[idx(i,j)] * X[jx];
310 jx ++;
311 }
312
313 X[ix] = tmp / LU[idx(i,i)];
314 ix --;
315 }
316
317 return UNUR_SUCCESS;
318
319 #undef idx
320 } /* end of _unur_matrix_backsubstitution_dtrsv() */
321
322 /*---------------------------------------------------------------------------*/
323
324 int
_unur_matrix_forwardsubstitution_dtrsv(int dim,double * LU,double * X)325 _unur_matrix_forwardsubstitution_dtrsv(int dim, double *LU, double *X)
326 /*----------------------------------------------------------------------*/
327 /* Forwardsubstitution used for inversion alg. _unur_matrix_LU_invert() */
328 /* */
329 /* input: */
330 /* dim ... number of columns and rows of */
331 /* LU ... dim x dim -matrix that contains LU decomposition of matrix */
332 /* X ... vector */
333 /* */
334 /* output: */
335 /* X */
336 /* */
337 /* return: */
338 /* UNUR_SUCCESS on success */
339 /* error code otherwise */
340 /*----------------------------------------------------------------------*/
341 {
342 #define idx(a,b) ((a)*dim+(b))
343
344 int ix,jx,i,j;
345
346 /* check arguments */
347 CHECK_NULL(LU,UNUR_ERR_NULL);
348 CHECK_NULL(X,UNUR_ERR_NULL);
349
350 /* forward substitution */
351 ix = 0;
352 ix++;
353 for (i = 1; i < dim; i++) {
354 double tmp = X[ix];
355 jx = 0;
356 for (j = 0; j < i; j++) {
357 tmp -= LU[idx(i,j)] * X[jx];
358 jx += 1;
359 }
360 X[ix] = tmp;
361 ix ++;
362 }
363
364 return UNUR_SUCCESS;
365
366 #undef idx
367 } /* end of _unur_matrix_forwardsubstitution_dtrsv() */
368
369 /*---------------------------------------------------------------------------*/
370
371 int
_unur_matrix_LU_invert(int dim,double * LU,int * p,double * inverse)372 _unur_matrix_LU_invert (int dim, double *LU, int *p, double *inverse)
373 /*----------------------------------------------------------------------*/
374 /* Compute inverse of matrix with given LU decomposition. */
375 /* */
376 /* input: */
377 /* dim ... number of columns and rows of */
378 /* LU ... dim x dim -matrix, LU decomposition */
379 /* p ... permutation vector of length dim */
380 /* (contains all integers between 0 and dim-1) */
381 /* inverse ... pointer to array where inverse should be stored */
382 /* */
383 /* output */
384 /* inverse ... the inverse matrix */
385 /* */
386 /* return: */
387 /* UNUR_SUCCESS on success */
388 /* error code otherwise */
389 /*----------------------------------------------------------------------*/
390 {
391 #define idx(a,b) ((a)*dim+(b))
392
393 double *vector;
394 int i,j;
395
396 /* check arguments */
397 CHECK_NULL(LU,UNUR_ERR_NULL);
398 CHECK_NULL(p,UNUR_ERR_NULL);
399 CHECK_NULL(inverse,UNUR_ERR_NULL);
400
401 /* allocate working array */
402 vector = _unur_xmalloc(dim*sizeof(double));
403
404 for (i = 0; i < dim; i++){
405 for(j=0;j<dim;j++) vector[j] = 0.;
406 vector[i] = 1.;
407 /* Solve for c using forward-substitution, L c = P b */
408 _unur_matrix_forwardsubstitution_dtrsv (dim, LU, vector);
409 /* Perform back-substitution, U x = c */
410 _unur_matrix_backsubstitution_dtrsv (dim, LU, vector);
411
412 for(j=0;j<dim;j++){
413 inverse[idx(j,p[i])] = vector[j]; /* set column vector of inverse matrix */
414 }
415 }
416
417 /* free working arrays */
418 free(vector);
419
420 return UNUR_SUCCESS;
421
422 #undef idx
423 } /* end of _unur_matrix_LU_invert() */
424
425 /*---------------------------------------------------------------------------*/
426
427 int
_unur_matrix_invert_matrix(int dim,const double * A,double * Ainv,double * det)428 _unur_matrix_invert_matrix(int dim, const double *A, double *Ainv, double *det)
429 /*----------------------------------------------------------------------*/
430 /* Calculates the inverse matrix (by means of LU decomposition). */
431 /* As a side effect det(A) is comuted. */
432 /* */
433 /* input: */
434 /* dim ... dimension of the square matrix A */
435 /* A ... dim x dim -matrix */
436 /* Ainv ... pointer to array where inverse matrix should be stored */
437 /* det ... pointer where det(A) should be stored */
438 /* */
439 /* output: */
440 /* Ainv ... inverse matrix of A */
441 /* det ... determinant of A */
442 /* */
443 /* return: */
444 /* UNUR_SUCCESS on success */
445 /* other error code, otherwise */
446 /*----------------------------------------------------------------------*/
447 {
448 #define idx(a,b) ((a)*dim+(b))
449
450 int *p, s, i;
451 double *LU; /* array for storing LU decomposition of matrix A */
452
453 /* check arguments */
454 CHECK_NULL(A,UNUR_ERR_NULL);
455 CHECK_NULL(Ainv,UNUR_ERR_NULL);
456 CHECK_NULL(det,UNUR_ERR_NULL);
457 if (dim<1) {
458 _unur_error("matrix",UNUR_ERR_GENERIC,"dimension < 1");
459 return UNUR_ERR_GENERIC;
460 }
461
462 /* allocate working space */
463 p = _unur_xmalloc(dim*sizeof(int));
464 LU = _unur_xmalloc(dim*dim*sizeof(double));
465
466 /* compute LU decomposition */
467 memcpy(LU, A, dim*dim*sizeof(double));
468 _unur_matrix_LU_decomp(dim, LU, p, &s);
469
470 /* compute determinant */
471 *det = s;
472 for(i=0;i<dim;i++)
473 *det *= LU[idx(i,i)];
474
475 /* compute inverse by means of LU factors */
476 _unur_matrix_LU_invert(dim, LU, p, Ainv);
477
478 /* free working space */
479 free(LU);
480 free(p);
481
482 return UNUR_SUCCESS;
483
484 #undef idx
485 } /* end of _unur_matrix_invert_matrix() */
486
487 /*---------------------------------------------------------------------------*/
488
489 int
_unur_matrix_multiplication(int dim,const double * A,const double * B,double * AB)490 _unur_matrix_multiplication(int dim, const double *A, const double *B, double *AB)
491 /*----------------------------------------------------------------------*/
492 /* Calculates the matrix multiplication of two matrices A and B */
493 /* */
494 /* input: */
495 /* dim ... dimension of the square matrix A */
496 /* A ... dim x dim -matrix */
497 /* B ... dim x dim -matrix */
498 /* */
499 /* output: */
500 /* AB ... A*B */
501 /* */
502 /* return: */
503 /* UNUR_SUCCESS on success */
504 /* other error code, otherwise */
505 /*----------------------------------------------------------------------*/
506 {
507 #define idx(a,b) ((a)*dim+(b))
508
509 int i, j, k;
510
511 /* check arguments */
512 CHECK_NULL(A,UNUR_ERR_NULL);
513 CHECK_NULL(B,UNUR_ERR_NULL);
514 CHECK_NULL(AB,UNUR_ERR_NULL);
515
516 if (dim<1) {
517 _unur_error("matrix",UNUR_ERR_GENERIC,"dimension < 1");
518 return UNUR_ERR_GENERIC;
519 }
520
521 /* compute product */
522 for(i=0;i<dim;i++)
523 for(j=0;j<dim;j++) {
524 AB[idx(i,j)]=0.;
525 for(k=0;k<dim;k++) {
526 AB[idx(i,j)] += A[idx(i,k)]*B[idx(k,j)];
527 }
528 }
529
530 return UNUR_SUCCESS;
531
532 #undef idx
533 } /* end of _unur_matrix_multiplication() */
534
535 /*---------------------------------------------------------------------------*/
536
537 double
_unur_matrix_determinant(int dim,const double * A)538 _unur_matrix_determinant ( int dim, const double *A )
539 /*----------------------------------------------------------------------*/
540 /* Calculates the determinant of the matrix A */
541 /* (by means of LU decomposition). */
542 /* input: */
543 /* dim ... dimension of the square matrix A */
544 /* A ... dim x dim -matrix */
545 /* */
546 /* return: */
547 /* determinant of A */
548 /* */
549 /* error: */
550 /* return INFINITY */
551 /*----------------------------------------------------------------------*/
552 {
553 #define idx(a,b) ((a)*dim+(b))
554
555 int *p, s, i;
556 double *LU; /* array for storing LU decomposition of matrix A */
557 double det;
558
559 /* check arguments */
560 CHECK_NULL(A, INFINITY);
561
562 /* one-dimensional case */
563 if (dim==1) return A[0];
564
565 /* allocate working space */
566 p = _unur_xmalloc(dim*sizeof(int));
567 LU = _unur_xmalloc(dim*dim*sizeof(double));
568
569 /* compute LU decomposition */
570 memcpy(LU, A, dim*dim*sizeof(double));
571 _unur_matrix_LU_decomp(dim, LU, p, &s);
572
573 /* compute determinant */
574 det = s;
575 for(i=0;i<dim;i++)
576 det *= LU[idx(i,i)];
577
578 /* free working space */
579 free(LU);
580 free(p);
581
582 return det;
583
584 #undef idx
585 } /* end of _unur_matrix_determinant() */
586
587 /*---------------------------------------------------------------------------*/
588
589 double
_unur_matrix_qf(int dim,double * x,double * A)590 _unur_matrix_qf (int dim, double *x, double *A)
591 /*----------------------------------------------------------------------*/
592 /* Compute quadratic form x'Ax */
593 /* */
594 /* input: */
595 /* dim ... number of columns and rows of */
596 /* x ... vector */
597 /* A ... dim x dim matrix */
598 /* */
599 /* return: */
600 /* returns the result of x'Ax */
601 /* */
602 /* error: */
603 /* return INFINITY */
604 /*----------------------------------------------------------------------*/
605 {
606 #define idx(a,b) ((a)*dim+(b))
607
608 int i,j;
609 double sum,outersum;
610
611 /* check arguments */
612 CHECK_NULL(x,INFINITY);
613 CHECK_NULL(A,INFINITY);
614 if (dim<1) {
615 _unur_error("matrix",UNUR_ERR_GENERIC,"dimension < 1");
616 return INFINITY;
617 }
618
619 outersum=0.;
620 for(i=0;i<dim;i++){
621 sum=0.;
622 for(j=0;j<dim;j++)
623 sum+=A[idx(i,j)]*x[j];
624 outersum+=sum*x[i];
625 }
626
627 return outersum;
628
629 #undef idx
630 } /* end of _unur_matrix_qf() */
631
632 /*---------------------------------------------------------------------------*/
633
634 int
_unur_matrix_cholesky_decomposition(int dim,const double * S,double * L)635 _unur_matrix_cholesky_decomposition (int dim, const double *S, double *L )
636 /*----------------------------------------------------------------------*/
637 /* The Colesky factor L of a variance-covariance matrix S is computed: */
638 /* S = LL' */
639 /* */
640 /* parameters: */
641 /* dim ... dimension of covariance matrix S */
642 /* S ... variance-covariance matrix */
643 /* L ... pointer to array where cholesky factor should be stored */
644 /* */
645 /* return: */
646 /* UNUR_SUCCESS on success */
647 /* UNUR_FAILURE when matrix is not positive definite */
648 /* other error code, otherwise */
649 /*----------------------------------------------------------------------*/
650 {
651 #define idx(a,b) ((a)*dim+(b))
652
653 int i,j,k;
654 double sum1,sum2;
655
656 /* check arguments */
657 CHECK_NULL(S,UNUR_ERR_NULL);
658 if (dim<1) {
659 _unur_error("matrix",UNUR_ERR_GENERIC,"dimension < 1");
660 return UNUR_ERR_GENERIC;
661 }
662
663 L[idx(0,0)] = sqrt( S[idx(0,0)] );
664
665 for(j=1; j<dim; j++) {
666 L[idx(j,0)] = S[idx(j,0)] / L[idx(0,0)];
667
668 sum1 = L[idx(j,0)] * L[idx(j,0)];
669 for(k=1; k<j; k++) {
670 sum2 = 0.;
671 for(i=0; i<k; i++)
672 sum2 += L[idx(j,i)] * L[idx(k,i)];
673
674 L[idx(j,k)] = (S[idx(j,k)] - sum2) / L[idx(k,k)];
675 sum1 += L[idx(j,k)] * L[idx(j,k)];
676 }
677
678 if (!(S[idx(j,j)] > sum1)) {
679 /* covariance matrix not positive definite */
680 return UNUR_FAILURE;
681 }
682
683 L[idx(j,j)] = sqrt( S[idx(j,j)] - sum1 );
684 }
685
686 /* although not necessary upper triangular of L - matrix is set to 0 */
687 for(j=0; j<dim; j++)
688 for(k=j+1; k<dim; k++)
689 L[idx(j,k)]=0.;
690
691 /* o.k. */
692 return UNUR_SUCCESS;
693
694 #undef idx
695 } /* end of _unur_matrix_cholesky_decomposition() */
696
697 /*--------------------------------------------------------------------------*/
698
699 void
_unur_matrix_print_vector(int dim,const double * vec,const char * info,FILE * LOG,const char * genid,const char * indent)700 _unur_matrix_print_vector ( int dim, const double *vec, const char *info,
701 FILE *LOG, const char *genid, const char *indent )
702 /*----------------------------------------------------------------------*/
703 /* Print elements of vector in a single row enclosed by parenthesis */
704 /* into LOG file. */
705 /* */
706 /* parameters: */
707 /* dim ... dimension */
708 /* vec ... vector with <dim> entries */
709 /* info ... additional info-string to be printed as first line */
710 /* LOG ... output stream */
711 /* genid ... id string */
712 /* indent ... left margin of printed vector */
713 /*----------------------------------------------------------------------*/
714 {
715 int i;
716
717 if (vec) {
718 fprintf(LOG,"%s: %s\n", genid, info );
719 fprintf(LOG,"%s: %s( %g", genid, indent, vec[0]);
720 for (i=1; i<dim; i++)
721 fprintf(LOG,", %g", vec[i]);
722 fprintf(LOG," )\n");
723 }
724 else {
725 fprintf(LOG,"%s: %s [unknown]\n", genid, info );
726 }
727
728 fprintf(LOG,"%s:\n",genid);
729 } /* end of _unur_matrix_print_vector() */
730
731 /*--------------------------------------------------------------------------*/
732
733 void
_unur_matrix_print_matrix(int dim,const double * mat,const char * info,FILE * LOG,const char * genid,const char * indent)734 _unur_matrix_print_matrix ( int dim, const double *mat, const char *info,
735 FILE *LOG, const char *genid, const char *indent )
736 /*----------------------------------------------------------------------*/
737 /* Print elements of the given <dim>x<dim> square matrix into LOG file. */
738 /* The matrix is stored row-wise in <mat>. */
739 /* */
740 /* parameters: */
741 /* dim ... dimension */
742 /* mat ... square matrix with <dim> rows and columns */
743 /* info ... additional info-string to be printed as first line */
744 /* LOG ... output stream */
745 /* genid ... id string */
746 /* indent ... left margin of printed vector */
747 /*----------------------------------------------------------------------*/
748 {
749 #define idx(a,b) ((a)*dim+(b))
750 int i,j;
751
752 if (mat) {
753 fprintf(LOG,"%s: %s\n", genid, info);
754 for (i=0; i<dim; i++) {
755 fprintf(LOG,"%s: %s(% e", genid, indent, mat[idx(i,0)]);
756 for (j=1; j<dim; j++)
757 fprintf(LOG,",% e",mat[idx(i,j)]);
758 fprintf(LOG," )\n");
759 }
760 }
761 else {
762 fprintf(LOG,"%s: %s [unknown]\n", genid, info );
763 }
764
765 fprintf(LOG,"%s:\n",genid);
766
767 #undef idx
768 } /* end of _unur_matrix_print_matrix() */
769
770 /*--------------------------------------------------------------------------*/
771