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.16 by greg, Mon Aug 12 20:38:19 2019 UTC

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

Diff Legend

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