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.20 by greg, Wed Mar 23 00:58:35 2022 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 + /* unary matrix operation(s) */
18   typedef struct {
19          double          sca[MAXCOMP];           /* scalar coefficients */
20 –        int             nsf;                    /* number of scalars */
20          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
21 <        int             clen;                   /* number of coefficients */
22 <        int             transpose;              /* do transpose? */
23 <        int             op;                     /* '*' or '+' */
24 < } ROPERAT;                              /* matrix operation */
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 < static void
38 < op_default(ROPERAT *op)
37 > /* Load matrix */
38 > static int
39 > loadmatrix(ROPMAT *rop)
40   {
41 <        memset(op, 0, sizeof(ROPERAT));
42 <        op->op = '.';
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 < operate(RMATRIX *mleft, ROPERAT *op, const char *fname)
55 > loadop(ROPMAT *rop)
56   {
57 <        RMATRIX *mright = rmx_load(fname);
40 <        RMATRIX *mtmp;
57 >        RMATRIX *mres;
58          int     i;
59  
60 <        if (fname == NULL)
44 <                fname = "<stdin>";
45 <        if (mright == NULL) {
46 <                fputs(fname, stderr);
47 <                fputs(": cannot load matrix\n", stderr);
60 >        if (loadmatrix(rop) < 0)                /* make sure we're loaded */
61                  return(NULL);
62 <        }
63 <        if (op->nsf > 0) {              /* apply scalar(s) */
64 <                if (op->clen > 0) {
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 <                        rmx_free(mright);
54 <                        return(NULL);
66 >                        goto failure;
67                  }
68 <                if (op->nsf == 1) {
69 <                        for (i = mright->ncomp; --i; )
70 <                                op->sca[i] = op->sca[0];
71 <                } else if (op->nsf != mright->ncomp) {
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 <                                        fname, mright->ncomp);
74 <                        rmx_free(mright);
63 <                        return(NULL);
73 >                                        rop->inspec, rop->mtx->ncomp);
74 >                        goto failure;
75                  }
76 <                if ((mleft == NULL) | (op->op != '+') &&
77 <                                !rmx_scale(mright, op->sca)) {
67 <                        fputs(fname, stderr);
76 >                if (!rmx_scale(rop->mtx, rop->preop.sca)) {
77 >                        fputs(rop->inspec, stderr);
78                          fputs(": scalar operation failed\n", stderr);
79 <                        rmx_free(mright);
70 <                        return(NULL);
79 >                        goto failure;
80                  }
81                  if (verbose) {
82 <                        fputs(fname, stderr);
82 >                        fputs(rop->inspec, stderr);
83                          fputs(": applied scalar (", stderr);
84 <                        for (i = 0; i < op->nsf; i++)
85 <                                fprintf(stderr, " %f", op->sca[i]);
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 (op->clen > 0) {             /* apply transform */
90 <                if (op->clen % mright->ncomp) {
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 <                                        fname, mright->ncomp);
93 <                        rmx_free(mright);
85 <                        return(NULL);
92 >                                        rop->inspec, rop->mtx->ncomp);
93 >                        goto failure;
94                  }
95 <                mtmp = rmx_transform(mright, op->clen/mright->ncomp, op->cmat);
96 <                if (mtmp == NULL) {
97 <                        fprintf(stderr, "%s: matrix transform failed\n", fname);
98 <                        rmx_free(mright);
99 <                        return(NULL);
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 <                                        fname, mtmp->ncomp, mright->ncomp);
105 <                rmx_free(mright);
106 <                mright = mtmp;
104 >                                        rop->inspec, mres->ncomp,
105 >                                        rop->mtx->ncomp);
106 >                rmx_free(rop->mtx);
107 >                rop->mtx = mres;
108          }
109 <        if (op->transpose) {            /* transpose matrix? */
110 <                mtmp = rmx_transpose(mright);
111 <                if (mtmp == NULL) {
112 <                        fputs(fname, stderr);
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 <                        rmx_free(mright);
105 <                        return(NULL);
114 >                        goto failure;
115                  }
116                  if (verbose) {
117 <                        fputs(fname, stderr);
117 >                        fputs(rop->inspec, stderr);
118                          fputs(": transposed rows and columns\n", stderr);
119                  }
120 <                rmx_free(mright);
121 <                mright = mtmp;
120 >                rmx_free(rop->mtx);
121 >                rop->mtx = mres;
122          }
123 <        if (mleft == NULL)              /* just one matrix */
124 <                return(mright);
125 <        if (op->op == '.') {            /* concatenate */
126 <                RMATRIX *mres = rmx_multiply(mleft, mright);
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(fname, stderr);
155 <                        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);
154 >                        fputs(inspec, stderr);
155 >                        fputs(": concatenation failed\n", stderr);
156                          return(NULL);
157                  }
158                  if (verbose) {
159 <                        fputs(fname, stderr);
159 >                        fputs(inspec, stderr);
160                          fputs(": concatenated matrix\n", stderr);
161                  }
162 <                rmx_free(mright);
163 <                rmx_free(mleft);
164 <                mleft = mres;
165 <        } else if (op->op == '+') {
136 <                if (!rmx_sum(mleft, mright, op->nsf ? op->sca : (double *)NULL)) {
137 <                        fputs(fname, stderr);
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(fname, stderr);
172 >                        fputs(inspec, stderr);
173                          fputs(": added in matrix\n", stderr);
174                  }
175                  rmx_free(mright);
176 <        } else if ((op->op == '*') | (op->op == '/')) {
177 <                const char *    tnam = (op->op == '/') ?
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->op == '/'))) {
183 >                if (!rmx_elemult(mleft, mright, (op == '/'))) {
184                          fprintf(stderr, "%s: element-wise %s failed\n",
185 <                                        fname, tnam);
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 <                                        fname, tnam);
193 >                                        inspec, tnam);
194                  else if (verbose)
195 <                        fprintf(stderr, "%s: element-wise %s\n", fname, tnam);
195 >                        fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
196                  rmx_free(mright);
197 <        } else {
198 <                fprintf(stderr, "%s: unknown operation '%c'\n", fname, op->op);
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;
# Line 179 | Line 288 | get_factors(double da[], int n, char *av[])
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 <        ROPERAT op;
317 >        int     stdin_used = 0;
318          int     i;
190 –                                        /* initialize */
191 –        op_default(&op);
319                                          /* get options and arguments */
320 <        for (i = 1; i < argc; i++)
320 >        for (i = 1; i < argc; i++) {
321                  if (argv[i][0] && !argv[i][1] &&
322 <                                strchr("+*/", argv[i][0]) != NULL) {
323 <                        op.op = argv[i][0];
324 <                } else if (argv[i][0] != '-' || !argv[i][1]) {
325 <                        char    *fname = NULL;  /* load matrix */
326 <                        if (argv[i][0] != '-')
327 <                                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);
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 <                        op_default(&op);        /* reset operator */
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 = !verbose;
348 >                                verbose++;
349                                  break;
350                          case 'f':
351                                  switch (argv[i][2]) {
# Line 230 | Line 366 | main(int argc, char *argv[])
366                                  }
367                                  break;
368                          case 't':
369 <                                op.transpose = 1;
369 >                                mop[nmats].preop.transpose = 1;
370                                  break;
371                          case 's':
372                                  if (n > MAXCOMP) n = MAXCOMP;
373 <                                op.nsf = get_factors(op.sca, n, argv+i+1);
374 <                                i += op.nsf;
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 <                                op.clen = get_factors(op.cmat, n, argv+i+1);
380 <                                i += op.clen;
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 (mres == NULL)               /* check that we got something */
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;
# Line 261 | Line 426 | main(int argc, char *argv[])
426                  fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
427                  return(1);
428          }
429 <        /* rmx_free(mres); mres = NULL; */
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 ..] m1 [+*/] .. > mres\n",
433 >        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
434                          argv[0]);
435          return(1);
436   }

Diff Legend

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