ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.13
Committed: Mon Aug 12 01:20:26 2019 UTC (4 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +246 -107 lines
Log Message:
Take advantage of matrix multiplication associativity

File Contents

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