ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
(Generate patch)

Comparing ray/src/util/rmtxop.c (file contents):
Revision 2.12 by greg, Mon Aug 27 23:03:05 2018 UTC vs.
Revision 2.17 by greg, Sat Dec 28 18:05:14 2019 UTC

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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines