1 /* Simple matrix operations.
2 *
3 * Contents:
4 * 1. The matrixops API.
5 * 2. Debugging and development tools
6 * 3. Unit tests
7 * 4. Test driver
8 *
9 * Compare:
10 * esl_vectorops: simple vector operations
11 * esl_dmatrix: matrix algebra, double precision
12 *
13 * Note:
14 * Don't use <const> qualifier on input matrices. C's rules
15 * for const qualifier on nested pointers (int **foo) are
16 * problematic.
17 */
18 #include "esl_config.h"
19
20 #include <stdlib.h>
21
22 #include "easel.h"
23 #include "esl_matrixops.h"
24 #include "esl_vectorops.h"
25
26 /* Function: esl_mat_{DFIC}Create()
27 * Synopsis: Create an MxN matrix of doubles (floats, ints, chars).
28 */
29 double **
esl_mat_DCreate(int M,int N)30 esl_mat_DCreate(int M, int N)
31 {
32 double **A = NULL;
33 int i;
34 int status;
35
36 ESL_DASSERT1(( M > 0 ));
37 ESL_DASSERT1(( N > 0 ));
38
39 ESL_ALLOC(A, sizeof(double *) * M);
40 A[0] = NULL;
41
42 ESL_ALLOC(A[0], sizeof(double) * M * N);
43 for (i = 1; i < M; i++)
44 A[i] = A[0] + i * N;
45 return A;
46
47 ERROR:
48 esl_mat_DDestroy(A);
49 return NULL;
50 }
51 float **
esl_mat_FCreate(int M,int N)52 esl_mat_FCreate(int M, int N)
53 {
54 float **A = NULL;
55 int i;
56 int status;
57
58 ESL_DASSERT1(( M > 0 ));
59 ESL_DASSERT1(( N > 0 ));
60
61 ESL_ALLOC(A, sizeof(float *) * M);
62 A[0] = NULL;
63
64 ESL_ALLOC(A[0], sizeof(float) * M * N);
65 for (i = 1; i < M; i++)
66 A[i] = A[0] + i * N;
67 return A;
68
69 ERROR:
70 esl_mat_FDestroy(A);
71 return NULL;
72 }
73 int **
esl_mat_ICreate(int M,int N)74 esl_mat_ICreate(int M, int N)
75 {
76 int **A = NULL;
77 int i;
78 int status;
79
80 ESL_DASSERT1(( M > 0 ));
81 ESL_DASSERT1(( N > 0 ));
82
83 ESL_ALLOC(A, sizeof(int *) * M);
84 A[0] = NULL;
85
86 ESL_ALLOC(A[0], sizeof(int) * M * N);
87 for (i = 1; i < M; i++)
88 A[i] = A[0] + i * N;
89
90 return A;
91
92 ERROR:
93 esl_mat_IDestroy(A);
94 return NULL;
95 }
96 char **
esl_mat_CCreate(int M,int N)97 esl_mat_CCreate(int M, int N)
98 {
99 char **A = NULL;
100 int i;
101 int status;
102
103 ESL_DASSERT1(( M > 0 ));
104 ESL_DASSERT1(( N > 0 ));
105
106 ESL_ALLOC(A, sizeof(char *) * M);
107 A[0] = NULL;
108
109 ESL_ALLOC(A[0], sizeof(char) * M * N);
110 for (i = 1; i < M; i++)
111 A[i] = A[0] + i * N;
112
113 return A;
114
115 ERROR:
116 esl_mat_CDestroy(A);
117 return NULL;
118 }
119
120
121 /* Function: esl_mat_{DFI}Clone()
122 * Synopsis: Duplicate a 2D matrix
123 * Incept: SRE, Fri 05 Apr 2019
124 *
125 * Purpose: Duplicate a matrix, into freshly allocated space.
126 * Return a ptr to the duplicate.
127 *
128 * Throws: <NULL> on allocation failure.
129 */
130 double **
esl_mat_DClone(double ** A,int M,int N)131 esl_mat_DClone(double **A, int M, int N)
132 {
133 double **B;
134
135 if ((B = esl_mat_DCreate(M, N)) == NULL) return NULL;
136 esl_mat_DCopy(A, M, N, B);
137 return B;
138 }
139 float **
esl_mat_FClone(float ** A,int M,int N)140 esl_mat_FClone(float **A, int M, int N)
141 {
142 float **B;
143
144 if ((B = esl_mat_FCreate(M, N)) == NULL) return NULL;
145 esl_mat_FCopy(A, M, N, B);
146 return B;
147 }
148 int **
esl_mat_IClone(int ** A,int M,int N)149 esl_mat_IClone(int **A, int M, int N)
150 {
151 int **B;
152
153 if ((B = esl_mat_ICreate(M, N)) == NULL) return NULL;
154 esl_mat_ICopy(A, M, N, B);
155 return B;
156 }
157
158
159
160
161 /* Function: esl_mat_{DFIC}GrowTo()
162 * Synopsis: Increase the allocation for a matrix
163 * Incept: SRE, Fri 06 Jul 2018 [World Cup, France v. Uruguay]
164 *
165 * Purpose: Reallocate an existing matrix to <M> rows and <N> columns.
166 *
167 * If only <M> has increased, and <N> is the same as it
168 * was, then the existing contents of <*ret_A> remain valid
169 * and unchanged. Newly allocated values are undefined, and
170 * caller needs to initialize them.
171 *
172 * If <N> is changed, then the layout of the matrix
173 * <*ret_A> is necessarily changed too, and it will need to
174 * be reset to new values. Caller should treat the entire
175 * matrix as undefined, and reinitialize all of it.
176 *
177 * Return: <eslOK> on success, and <*ret_A> has been reallocated.
178 *
179 * Throws: <eslEMEM> on allocation failure; now <*ret_A> remains
180 * valid with its previous size and contents.
181 */
182 int
esl_mat_DGrowTo(double *** ret_A,int M,int N)183 esl_mat_DGrowTo(double ***ret_A, int M, int N)
184 {
185 double **A = *ret_A;
186 int i;
187 int status;
188
189 ESL_REALLOC(A[0], sizeof(double) * (M*N)); // must reallocate contents first
190 ESL_REALLOC(A, sizeof(double *) * M); // ... then the pointers
191 for (i = 1; i < M; i++) // ... then reset row pointers.
192 A[i] = A[0] + i * N;
193
194 *ret_A = A;
195 return eslOK;
196
197 ERROR:
198 *ret_A = A; // 1st realloc could succeed, moving A, before 2nd realloc fails.
199 return status;
200 }
201 int
esl_mat_FGrowTo(float *** ret_A,int M,int N)202 esl_mat_FGrowTo(float ***ret_A, int M, int N)
203 {
204 float **A = *ret_A;
205 int i;
206 int status;
207
208 ESL_REALLOC(A[0], sizeof(float) * (M*N));
209 ESL_REALLOC(A, sizeof(float *) * M);
210 for (i = 1; i < M; i++) A[i] = A[0] + i * N;
211 *ret_A = A;
212 return eslOK;
213
214 ERROR:
215 *ret_A = A;
216 return status;
217 }
218 int
esl_mat_IGrowTo(int *** ret_A,int M,int N)219 esl_mat_IGrowTo(int ***ret_A, int M, int N)
220 {
221 int **A = *ret_A;
222 int i;
223 int status;
224
225 ESL_REALLOC(A[0], sizeof(int) * (M*N));
226 ESL_REALLOC(A, sizeof(int *) * M);
227 for (i = 1; i < M; i++) A[i] = A[0] + i * N;
228 *ret_A = A;
229 return eslOK;
230
231 ERROR:
232 *ret_A = A;
233 return status;
234 }
235 int
esl_mat_CGrowTo(char *** ret_A,int M,int N)236 esl_mat_CGrowTo(char ***ret_A, int M, int N)
237 {
238 char **A = *ret_A;
239 int i;
240 int status;
241
242 ESL_REALLOC(A[0], sizeof(char) * (M*N));
243 ESL_REALLOC(A, sizeof(char *) * M);
244 for (i = 1; i < M; i++) A[i] = A[0] + i * N;
245 *ret_A = A;
246 return eslOK;
247
248 ERROR:
249 *ret_A = A;
250 return status;
251 }
252
253
254
255
256
257
258
259 /* Function: esl_mat_{DFIC}Sizeof()
260 * Synopsis: Returns size of a matrix, in bytes.
261 *
262 * Note: Doesn't need a particular matrix to calculate this;
263 */
264 size_t
esl_mat_DSizeof(int M,int N)265 esl_mat_DSizeof(int M, int N)
266 {
267 size_t n = 0;
268 n += sizeof(double) * M * N;
269 n += sizeof(double *) * M;
270 return n;
271 }
272 size_t
esl_mat_FSizeof(int M,int N)273 esl_mat_FSizeof(int M, int N)
274 {
275 size_t n = 0;
276 n += sizeof(float) * M * N;
277 n += sizeof(float *) * M;
278 return n;
279 }
280 size_t
esl_mat_ISizeof(int M,int N)281 esl_mat_ISizeof(int M, int N)
282 {
283 size_t n = 0;
284 n += sizeof(int) * M * N;
285 n += sizeof(int *) * M;
286 return n;
287 }
288 size_t
esl_mat_CSizeof(int M,int N)289 esl_mat_CSizeof(int M, int N)
290 {
291 size_t n = 0;
292 n += sizeof(char) * M * N;
293 n += sizeof(char *) * M;
294 return n;
295 }
296
297
298 /* Function: esl_mat_{DFI}Set()
299 * Synopsis: Set all values in a matrix to one number.
300 */
301 void
esl_mat_DSet(double ** A,int M,int N,double value)302 esl_mat_DSet(double **A, int M, int N, double value)
303 {
304 esl_vec_DSet(A[0], M*N, value);
305 }
306 void
esl_mat_FSet(float ** A,int M,int N,float value)307 esl_mat_FSet(float **A, int M, int N, float value)
308 {
309 esl_vec_FSet(A[0], M*N, value);
310 }
311 void
esl_mat_ISet(int ** A,int M,int N,int value)312 esl_mat_ISet(int **A, int M, int N, int value)
313 {
314 esl_vec_ISet(A[0], M*N, value);
315 }
316
317
318 /* Function: esl_mat_{DFI}Scale()
319 * Synopsis: Scale all values in a matrix by one scalar multiplier
320 * Incept: SRE, Wed 03 Apr 2019
321 */
322 void
esl_mat_DScale(double ** A,int M,int N,double x)323 esl_mat_DScale(double **A, int M, int N, double x)
324 {
325 esl_vec_DScale(A[0], M*N, x);
326 }
327 void
esl_mat_FScale(float ** A,int M,int N,float x)328 esl_mat_FScale(float **A, int M, int N, float x)
329 {
330 esl_vec_FScale(A[0], M*N, x);
331 }
332 void
esl_mat_IScale(int ** A,int M,int N,int x)333 esl_mat_IScale(int **A, int M, int N, int x)
334 {
335 esl_vec_IScale(A[0], M*N, x);
336 }
337
338
339
340 /* Function: esl_mat_{DFIWB}Copy()
341 * Synopsis: Copy <src> matrix to <dest>.
342 */
343 void
esl_mat_DCopy(double ** src,int M,int N,double ** dest)344 esl_mat_DCopy(double **src, int M, int N, double **dest)
345 {
346 esl_vec_DCopy(src[0], M*N, dest[0]);
347 }
348 void
esl_mat_FCopy(float ** src,int M,int N,float ** dest)349 esl_mat_FCopy(float **src, int M, int N, float **dest)
350 {
351 esl_vec_FCopy(src[0], M*N, dest[0]);
352 }
353 void
esl_mat_ICopy(int ** src,int M,int N,int ** dest)354 esl_mat_ICopy(int **src, int M, int N, int **dest)
355 {
356 esl_vec_ICopy(src[0], M*N, dest[0]);
357 }
358 void
esl_mat_WCopy(int16_t ** src,int M,int N,int16_t ** dest)359 esl_mat_WCopy(int16_t **src, int M, int N, int16_t **dest)
360 {
361 esl_vec_WCopy(src[0], M*N, dest[0]);
362 }
363 void
esl_mat_BCopy(int8_t ** src,int M,int N,int8_t ** dest)364 esl_mat_BCopy(int8_t **src, int M, int N, int8_t **dest)
365 {
366 esl_vec_BCopy(src[0], M*N, dest[0]);
367 }
368
369
370 /* Function: esl_mat_{DFI}Max()
371 * Synopsis: Return max value in a matrix.
372 */
373 double
esl_mat_DMax(double ** A,int M,int N)374 esl_mat_DMax(double **A, int M, int N)
375 {
376 return esl_vec_DMax(A[0], M*N);
377 }
378 float
esl_mat_FMax(float ** A,int M,int N)379 esl_mat_FMax(float **A, int M, int N)
380 {
381 return esl_vec_FMax(A[0], M*N);
382 }
383 int
esl_mat_IMax(int ** A,int M,int N)384 esl_mat_IMax(int **A, int M, int N)
385 {
386 return esl_vec_IMax(A[0], M*N);
387 }
388
389
390 /* Function: esl_mat_{DFI}Compare()
391 * Synopsis: Compare two matrices for equality.
392 * Incept: SRE, Tue 26 Jun 2018
393 *
394 * Purpose: Returns <eslOK> if two matrices <A>, <B> are
395 * identical (integer version) or identical
396 * within tolerance <tol> (double, float
397 * version; using <esl_{DF}Compare()>);
398 * Otherwise return <eslFAIL>.
399 */
400 int
esl_mat_DCompare(double ** A,double ** B,int M,int N,double tol)401 esl_mat_DCompare(double **A, double **B, int M, int N, double tol)
402 {
403 return esl_vec_DCompare(A[0], B[0], M*N, tol);
404 }
405 int
esl_mat_FCompare(float ** A,float ** B,int M,int N,float tol)406 esl_mat_FCompare(float **A, float **B, int M, int N, float tol)
407 {
408 return esl_vec_FCompare(A[0], B[0], M*N, tol);
409 }
410 int
esl_mat_ICompare(int ** A,int ** B,int M,int N)411 esl_mat_ICompare(int **A, int **B, int M, int N)
412 {
413 return esl_vec_ICompare(A[0], B[0], M*N);
414 }
415
416
417 /* Function: esl_mat_{DFIC}Destroy()
418 * Synopsis: Free a matrix.
419 * Incept: SRE, Tue 26 Jun 2018
420 *
421 * Note: If you try to use a single esl_mat_Destroy(void **A...),
422 * compiler will complain about incompatible pointer types;
423 * caller would have to explicitly cast (void **) A.
424 */
425 void
esl_mat_DDestroy(double ** A)426 esl_mat_DDestroy(double **A)
427 {
428 if (A) {
429 free(A[0]);
430 free(A);
431 }
432 }
433 void
esl_mat_FDestroy(float ** A)434 esl_mat_FDestroy(float **A)
435 {
436 if (A) {
437 free(A[0]);
438 free(A);
439 }
440 }
441 void
esl_mat_IDestroy(int ** A)442 esl_mat_IDestroy(int **A)
443 {
444 if (A) {
445 free(A[0]);
446 free(A);
447 }
448 }
449 void
esl_mat_CDestroy(char ** A)450 esl_mat_CDestroy(char **A)
451 {
452 if (A) {
453 free(A[0]);
454 free(A);
455 }
456 }
457
458
459
460 /*****************************************************************
461 * 2. Debugging and development tools
462 *****************************************************************/
463
464 int
esl_mat_DDump(double ** A,int M,int N)465 esl_mat_DDump(double **A, int M, int N)
466 {
467 int i,j;
468 for (i = 0; i < M; i++)
469 for (j = 0; j < N; j++)
470 printf("%10.4g%c", A[i][j], j==N-1? '\n' : ' ');
471 return eslOK;
472 }
473 int
esl_mat_FDump(float ** A,int M,int N)474 esl_mat_FDump(float **A, int M, int N)
475 {
476 int i,j;
477 for (i = 0; i < M; i++)
478 for (j = 0; j < N; j++)
479 printf("%10.4g%c", A[i][j], j==N-1? '\n' : ' ');
480 return eslOK;
481 }
482 int
esl_mat_IDump(int ** A,int M,int N)483 esl_mat_IDump(int **A, int M, int N)
484 {
485 int i,j;
486 for (i = 0; i < M; i++)
487 for (j = 0; j < N; j++)
488 printf("%8d%c", A[i][j], j==N-1? '\n' : ' ');
489 return eslOK;
490 }
491
492
493
494
495 /*****************************************************************
496 * 3. Unit tests
497 *****************************************************************/
498 #ifdef eslMATRIXOPS_TESTDRIVE
499
500 #include "esl_random.h"
501
502 /* utest_idiocy()
503 *
504 * Basically just a valgrind and compile warning test. Inconceivable
505 * that it could fail. (Why yes, I do know what that word means.)
506 */
507 static void
utest_idiocy(ESL_RANDOMNESS * rng)508 utest_idiocy(ESL_RANDOMNESS *rng)
509 {
510 char msg[] = "esl_matrixops utest_idiocy() test failed";
511 int m = 1 + esl_rnd_Roll(rng, 10); // 1..10
512 int n = 1 + esl_rnd_Roll(rng, 10); // 1..10
513 double **D = esl_mat_DCreate(m, n);
514 float **F = esl_mat_FCreate(m, n);
515 int **A = esl_mat_ICreate(m, n);
516 char **S = esl_mat_CCreate(m, n);
517 double **D2 = esl_mat_DCreate(m, n);
518 float **F2 = esl_mat_FCreate(m, n);
519 int **A2 = esl_mat_ICreate(m, n);
520 char **S2 = esl_mat_CCreate(m, n);
521 int testval = 42;
522
523 esl_mat_DSet(D, m, n, (double) testval);
524 esl_mat_FSet(F, m, n, (float) testval);
525 esl_mat_ISet(A, m, n, testval);
526
527 if (esl_mat_DMax(D, m, n) != (double) testval) esl_fatal(msg);
528 if (esl_mat_FMax(F, m, n) != (float) testval) esl_fatal(msg);
529 if (esl_mat_IMax(A, m, n) != testval) esl_fatal(msg);
530
531 esl_mat_DCopy(D, m, n, D2);
532 esl_mat_FCopy(F, m, n, F2);
533 esl_mat_ICopy(A, m, n, A2);
534
535 if (esl_mat_DCompare(D, D2, m, n, 1e-8) != eslOK) esl_fatal(msg);
536 if (esl_mat_FCompare(F, F2, m, n, 1e-4) != eslOK) esl_fatal(msg);
537 if (esl_mat_ICompare(A, A2, m, n) != eslOK) esl_fatal(msg);
538
539 esl_mat_DDestroy(D); esl_mat_DDestroy(D2);
540 esl_mat_FDestroy(F); esl_mat_FDestroy(F2);
541 esl_mat_IDestroy(A); esl_mat_IDestroy(A2);
542 esl_mat_CDestroy(S); esl_mat_CDestroy(S2);
543 }
544
545 static void
utest_grow(void)546 utest_grow(void)
547 {
548 char msg[] = "esl_matrixops utest_grow() test failed";
549 double **D1 = esl_mat_DCreate(5, 3);
550 double **D2 = esl_mat_DCreate(10, 3);
551 int i,j;
552
553 for (i = 0; i < 5; i++)
554 for (j = 0; j < 3; j++)
555 D1[i][j] = (double) (i*3 + j);
556 esl_mat_DGrowTo(&D1, 10, 3);
557 for (i = 5; i < 10; i++)
558 for (j = 0; j < 3; j++)
559 D1[i][j] = (double) (i*3 + j);
560
561 for (i = 0; i < 10; i++)
562 for (j = 0; j < 3; j++)
563 D2[i][j] = (double) (i*3 + j);
564
565 if (esl_mat_DCompare(D1, D2, 10, 3, 1e-5) != eslOK) esl_fatal(msg);
566
567 esl_mat_DDestroy(D1);
568 esl_mat_DDestroy(D2);
569 }
570 #endif // eslMATRIXOPS_TESTDRIVE
571
572
573 /*****************************************************************
574 * 4. Test driver
575 *****************************************************************/
576 #ifdef eslMATRIXOPS_TESTDRIVE
577
578 #include "easel.h"
579 #include "esl_getopts.h"
580 #include "esl_random.h"
581
582 static ESL_OPTIONS options[] = {
583 /* name type default env range togs reqs incomp help docgrp */
584 {"-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show help and usage", 0},
585 {"-s", eslARG_INT, "0", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0},
586 { 0,0,0,0,0,0,0,0,0,0},
587 };
588 static char usage[] = "[-options]";
589 static char banner[] = "test driver for matrixops module";
590
591 int
main(int argc,char ** argv)592 main(int argc, char **argv)
593 {
594 ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
595 ESL_RANDOMNESS *rng = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
596
597 fprintf(stderr, "## %s\n", argv[0]);
598 fprintf(stderr, "# rng seed = %" PRIu32 "\n", esl_randomness_GetSeed(rng));
599
600 utest_idiocy(rng);
601 utest_grow();
602
603 fprintf(stderr, "# status = ok\n");
604 esl_randomness_Destroy(rng);
605 esl_getopts_Destroy(go);
606 return 0;
607 }
608 #endif // eslMATRIXOPS_TESTDRIVE
609