ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.20
Committed: Wed Mar 23 00:58:35 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4
Changes since 2.19: +3 -3 lines
Log Message:
fix(rmtxop): Made error reporting consistent regardless of concatination direction

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmtxop.c,v 2.19 2021/01/21 17:47:07 greg Exp $";
3 #endif
4 /*
5 * General component matrix operations.
6 */
7
8 #include <stdlib.h>
9 #include <errno.h>
10 #include "rtio.h"
11 #include "resolu.h"
12 #include "rmatrix.h"
13 #include "platform.h"
14
15 #define MAXCOMP 16 /* #components we support */
16
17 /* unary matrix operation(s) */
18 typedef struct {
19 double sca[MAXCOMP]; /* scalar coefficients */
20 double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
21 short nsf; /* number of scalars */
22 short clen; /* number of coefficients */
23 short transpose; /* do transpose? */
24 } RUNARYOP;
25
26 /* matrix input source and requested operation(s) */
27 typedef struct {
28 const char *inspec; /* input specification */
29 RMPref rmp; /* matrix preference */
30 RUNARYOP preop; /* unary operation(s) */
31 RMATRIX *mtx; /* original matrix if loaded */
32 int binop; /* binary op with next (or 0) */
33 } ROPMAT;
34
35 int verbose = 0; /* verbose reporting? */
36
37 /* Load matrix */
38 static int
39 loadmatrix(ROPMAT *rop)
40 {
41 if (rop->mtx != NULL) /* already loaded? */
42 return(0);
43
44 rop->mtx = rmx_load(rop->inspec, rop->rmp);
45 if (rop->mtx == NULL) {
46 fputs(rop->inspec, stderr);
47 fputs(": cannot load matrix\n", stderr);
48 return(-1);
49 }
50 return(1);
51 }
52
53 /* Get matrix and perform unary operations */
54 static RMATRIX *
55 loadop(ROPMAT *rop)
56 {
57 RMATRIX *mres;
58 int i;
59
60 if (loadmatrix(rop) < 0) /* make sure we're loaded */
61 return(NULL);
62
63 if (rop->preop.nsf > 0) { /* apply scalar(s) */
64 if (rop->preop.clen > 0) {
65 fputs("Options -s and -c are exclusive\n", stderr);
66 goto failure;
67 }
68 if (rop->preop.nsf == 1) {
69 for (i = rop->mtx->ncomp; --i; )
70 rop->preop.sca[i] = rop->preop.sca[0];
71 } else if (rop->preop.nsf != rop->mtx->ncomp) {
72 fprintf(stderr, "%s: -s must have one or %d factors\n",
73 rop->inspec, rop->mtx->ncomp);
74 goto failure;
75 }
76 if (!rmx_scale(rop->mtx, rop->preop.sca)) {
77 fputs(rop->inspec, stderr);
78 fputs(": scalar operation failed\n", stderr);
79 goto failure;
80 }
81 if (verbose) {
82 fputs(rop->inspec, stderr);
83 fputs(": applied scalar (", stderr);
84 for (i = 0; i < rop->preop.nsf; i++)
85 fprintf(stderr, " %f", rop->preop.sca[i]);
86 fputs(" )\n", stderr);
87 }
88 }
89 if (rop->preop.clen > 0) { /* apply transform */
90 if (rop->preop.clen % rop->mtx->ncomp) {
91 fprintf(stderr, "%s: -c must have N x %d coefficients\n",
92 rop->inspec, rop->mtx->ncomp);
93 goto failure;
94 }
95 mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
96 rop->preop.cmat);
97 if (mres == NULL) {
98 fprintf(stderr, "%s: matrix transform failed\n",
99 rop->inspec);
100 goto failure;
101 }
102 if (verbose)
103 fprintf(stderr, "%s: applied %d x %d transform\n",
104 rop->inspec, mres->ncomp,
105 rop->mtx->ncomp);
106 rmx_free(rop->mtx);
107 rop->mtx = mres;
108 }
109 if (rop->preop.transpose) { /* transpose matrix? */
110 mres = rmx_transpose(rop->mtx);
111 if (mres == NULL) {
112 fputs(rop->inspec, stderr);
113 fputs(": transpose failed\n", stderr);
114 goto failure;
115 }
116 if (verbose) {
117 fputs(rop->inspec, stderr);
118 fputs(": transposed rows and columns\n", stderr);
119 }
120 rmx_free(rop->mtx);
121 rop->mtx = mres;
122 }
123 mres = rop->mtx;
124 rop->mtx = NULL;
125 return(mres);
126 failure:
127 rmx_free(rop->mtx);
128 return(rop->mtx = NULL);
129 }
130
131 /* Execute binary operation, free matrix arguments and return new result */
132 static RMATRIX *
133 binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
134 {
135 RMATRIX *mres = NULL;
136 int i;
137
138 if ((mleft == NULL) | (mright == NULL))
139 return(NULL);
140 switch (op) {
141 case '.': /* concatenate */
142 if (mleft->ncomp != mright->ncomp) {
143 fputs(inspec, stderr);
144 fputs(": # components do not match\n", stderr);
145 } else if (mleft->ncols != mright->nrows) {
146 fputs(inspec, stderr);
147 fputs(": mismatched dimensions\n",
148 stderr);
149 } else
150 mres = rmx_multiply(mleft, mright);
151 rmx_free(mleft);
152 rmx_free(mright);
153 if (mres == NULL) {
154 fputs(inspec, stderr);
155 fputs(": concatenation failed\n", stderr);
156 return(NULL);
157 }
158 if (verbose) {
159 fputs(inspec, stderr);
160 fputs(": concatenated matrix\n", stderr);
161 }
162 break;
163 case '+':
164 if (!rmx_sum(mleft, mright, NULL)) {
165 fputs(inspec, stderr);
166 fputs(": matrix sum failed\n", stderr);
167 rmx_free(mleft);
168 rmx_free(mright);
169 return(NULL);
170 }
171 if (verbose) {
172 fputs(inspec, stderr);
173 fputs(": added in matrix\n", stderr);
174 }
175 rmx_free(mright);
176 mres = mleft;
177 break;
178 case '*':
179 case '/': {
180 const char * tnam = (op == '/') ?
181 "division" : "multiplication";
182 errno = 0;
183 if (!rmx_elemult(mleft, mright, (op == '/'))) {
184 fprintf(stderr, "%s: element-wise %s failed\n",
185 inspec, tnam);
186 rmx_free(mleft);
187 rmx_free(mright);
188 return(NULL);
189 }
190 if (errno)
191 fprintf(stderr,
192 "%s: warning - error during element-wise %s\n",
193 inspec, tnam);
194 else if (verbose)
195 fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
196 rmx_free(mright);
197 mres = mleft;
198 } break;
199 default:
200 fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
201 rmx_free(mleft);
202 rmx_free(mright);
203 return(NULL);
204 }
205 return(mres);
206 }
207
208 /* Perform matrix operations from left to right */
209 static RMATRIX *
210 op_left2right(ROPMAT *mop)
211 {
212 RMATRIX *mleft = loadop(mop);
213
214 while (mop->binop) {
215 if (mleft == NULL)
216 break;
217 mleft = binaryop(mop[1].inspec,
218 mleft, mop->binop, loadop(mop+1));
219 mop++;
220 }
221 return(mleft);
222 }
223
224 /* Perform matrix operations from right to left */
225 static RMATRIX *
226 op_right2left(ROPMAT *mop)
227 {
228 RMATRIX *mright;
229 int rpos = 0;
230 /* find end of list */
231 while (mop[rpos].binop)
232 if (mop[rpos++].binop != '.') {
233 fputs(
234 "Right-to-left evaluation only for matrix multiplication!\n",
235 stderr);
236 return(NULL);
237 }
238 mright = loadop(mop+rpos);
239 while (rpos-- > 0) {
240 if (mright == NULL)
241 break;
242 mright = binaryop(mop[rpos+1].inspec,
243 loadop(mop+rpos), mop[rpos].binop, mright);
244 }
245 return(mright);
246 }
247
248 #define t_nrows(mop) ((mop)->preop.transpose ? (mop)->mtx->ncols \
249 : (mop)->mtx->nrows)
250 #define t_ncols(mop) ((mop)->preop.transpose ? (mop)->mtx->nrows \
251 : (mop)->mtx->ncols)
252
253 /* Should we prefer concatenating from rightmost matrix towards left? */
254 static int
255 prefer_right2left(ROPMAT *mop)
256 {
257 int mri = 0;
258
259 while (mop[mri].binop) /* find rightmost matrix */
260 if (mop[mri++].binop != '.')
261 return(0); /* pre-empt reversal for other ops */
262
263 if (mri <= 1)
264 return(0); /* won't matter */
265
266 if (loadmatrix(mop+mri) < 0) /* load rightmost cat */
267 return(1); /* fail will bail in a moment */
268
269 if (t_ncols(mop+mri) == 1)
270 return(1); /* definitely better R->L */
271
272 if (t_ncols(mop+mri) >= t_nrows(mop+mri))
273 return(0); /* ...probably worse */
274
275 if (loadmatrix(mop) < 0) /* load leftmost */
276 return(0); /* fail will bail in a moment */
277
278 return(t_ncols(mop+mri) < t_nrows(mop));
279 }
280
281 static int
282 get_factors(double da[], int n, char *av[])
283 {
284 int ac;
285
286 for (ac = 0; ac < n && isflt(av[ac]); ac++)
287 da[ac] = atof(av[ac]);
288 return(ac);
289 }
290
291 static ROPMAT *
292 grow_moparray(ROPMAT *mop, int n2alloc)
293 {
294 int nmats = 0;
295
296 while (mop[nmats++].binop)
297 ;
298 mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
299 if (mop == NULL) {
300 fputs("Out of memory in grow_moparray()\n", stderr);
301 exit(1);
302 }
303 if (n2alloc > nmats)
304 memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
305 return(mop);
306 }
307
308 /* Load one or more matrices and operate on them, sending results to stdout */
309 int
310 main(int argc, char *argv[])
311 {
312 int outfmt = DTfromHeader;
313 int nall = 2;
314 ROPMAT *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
315 int nmats = 0;
316 RMATRIX *mres = NULL;
317 int stdin_used = 0;
318 int i;
319 /* get options and arguments */
320 for (i = 1; i < argc; i++) {
321 if (argv[i][0] && !argv[i][1] &&
322 strchr(".+*/", argv[i][0]) != NULL) {
323 if (!nmats || mop[nmats-1].binop) {
324 fprintf(stderr,
325 "%s: missing matrix argument before '%c' operation\n",
326 argv[0], argv[i][0]);
327 return(1);
328 }
329 mop[nmats-1].binop = argv[i][0];
330 } else if (argv[i][0] != '-' || !argv[i][1]) {
331 if (argv[i][0] == '-') {
332 if (stdin_used++) {
333 fprintf(stderr,
334 "%s: standard input used for more than one matrix\n",
335 argv[0]);
336 return(1);
337 }
338 mop[nmats].inspec = stdin_name;
339 } else
340 mop[nmats].inspec = argv[i];
341 if (nmats > 0 && !mop[nmats-1].binop)
342 mop[nmats-1].binop = '.';
343 nmats++;
344 } else {
345 int n = argc-1 - i;
346 switch (argv[i][1]) { /* get option */
347 case 'v':
348 verbose++;
349 break;
350 case 'f':
351 switch (argv[i][2]) {
352 case 'd':
353 outfmt = DTdouble;
354 break;
355 case 'f':
356 outfmt = DTfloat;
357 break;
358 case 'a':
359 outfmt = DTascii;
360 break;
361 case 'c':
362 outfmt = DTrgbe;
363 break;
364 default:
365 goto userr;
366 }
367 break;
368 case 't':
369 mop[nmats].preop.transpose = 1;
370 break;
371 case 's':
372 if (n > MAXCOMP) n = MAXCOMP;
373 i += mop[nmats].preop.nsf =
374 get_factors(mop[nmats].preop.sca,
375 n, argv+i+1);
376 break;
377 case 'c':
378 if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
379 i += mop[nmats].preop.clen =
380 get_factors(mop[nmats].preop.cmat,
381 n, argv+i+1);
382 break;
383 case 'r':
384 if (argv[i][2] == 'f')
385 mop[nmats].rmp = RMPreflF;
386 else if (argv[i][2] == 'b')
387 mop[nmats].rmp = RMPreflB;
388 else
389 goto userr;
390 break;
391 default:
392 fprintf(stderr, "%s: unknown operation '%s'\n",
393 argv[0], argv[i]);
394 goto userr;
395 }
396 }
397 if (nmats >= nall)
398 mop = grow_moparray(mop, nall += 2);
399 }
400 if (mop[0].inspec == NULL) /* nothing to do? */
401 goto userr;
402 if (mop[nmats-1].binop) {
403 fprintf(stderr,
404 "%s: missing matrix argument after '%c' operation\n",
405 argv[0], mop[nmats-1].binop);
406 return(1);
407 }
408 /* favor quicker concatenation */
409 mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
410 : op_left2right(mop);
411 if (mop[nmats].mtx == NULL)
412 return(1);
413 /* apply trailing unary operations */
414 mop[nmats].inspec = "trailing_ops";
415 mres = loadop(mop+nmats);
416 if (mres == NULL)
417 return(1);
418 /* write result to stdout */
419 if (outfmt == DTfromHeader)
420 outfmt = mres->dtype;
421 if (outfmt != DTascii)
422 SET_FILE_BINARY(stdout);
423 newheader("RADIANCE", stdout);
424 printargs(argc, argv, stdout);
425 if (!rmx_write(mres, outfmt, stdout)) {
426 fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
427 return(1);
428 }
429 /* rmx_free(mres); free(mop); */
430 return(0);
431 userr:
432 fprintf(stderr,
433 "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
434 argv[0]);
435 return(1);
436 }