src/util/rmtxop.c

#ifndef lint
static const char RCSid[] = "$Id: rmtxop.c,v 2.15 2019/08/12 16:55:24 greg Exp $";
#endif
/*
 * General component matrix operations.
 */

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include "rtio.h"
#include "resolu.h"
#include "rmatrix.h"
#include "platform.h"

#define MAXCOMP         16              /* #components we support */

static const char       stdin_name[] = "<stdin>";

/* unary matrix operation(s) */
typedef struct {
        double          sca[MAXCOMP];           /* scalar coefficients */
        double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
        short           nsf;                    /* number of scalars */
        short           clen;                   /* number of coefficients */
        short           transpose;              /* do transpose? */
} RUNARYOP;

/* matrix input source and requested operation(s) */
typedef struct {
        const char      *inspec;                /* input specification */
        RUNARYOP        preop;                  /* unary operation(s) */
        RMATRIX         *mtx;                   /* original matrix if loaded */
        int             binop;                  /* binary op with next (or 0) */
} ROPMAT;

int     verbose = 0;                    /* verbose reporting? */

/* Load matrix */
static int
loadmatrix(ROPMAT *rop)
{
        if (rop->mtx != NULL)
                return(0);

        rop->mtx = rmx_load(rop->inspec == stdin_name ?
                                (const char *)NULL : rop->inspec);
        if (rop->mtx == NULL) {
                fputs(rop->inspec, stderr);
                fputs(": cannot load matrix\n", stderr);
                return(-1);
        }
        return(1);
}

/* Get matrix and perform unary operations */
static RMATRIX *
loadop(ROPMAT *rop)
{
        RMATRIX *mres;
        int     i;

        if (loadmatrix(rop) < 0)                /* make sure we're loaded */
                return(NULL);

        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
                if (rop->preop.clen > 0) {
                        fputs("Options -s and -c are exclusive\n", stderr);
                        goto failure;
                }
                if (rop->preop.nsf == 1) {
                        for (i = rop->mtx->ncomp; --i; )
                                rop->preop.sca[i] = rop->preop.sca[0];
                } else if (rop->preop.nsf != rop->mtx->ncomp) {
                        fprintf(stderr, "%s: -s must have one or %d factors\n",
                                        rop->inspec, rop->mtx->ncomp);
                        goto failure;
                }
                if (!rmx_scale(rop->mtx, rop->preop.sca)) {
                        fputs(rop->inspec, stderr);
                        fputs(": scalar operation failed\n", stderr);
                        goto failure;
                }
                if (verbose) {
                        fputs(rop->inspec, stderr);
                        fputs(": applied scalar (", stderr);
                        for (i = 0; i < rop->preop.nsf; i++)
                                fprintf(stderr, " %f", rop->preop.sca[i]);
                        fputs(" )\n", stderr);
                }
        }
        if (rop->preop.clen > 0) {      /* apply transform */
                if (rop->preop.clen % rop->mtx->ncomp) {
                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
                                        rop->inspec, rop->mtx->ncomp);
                        goto failure;
                }
                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
                                        rop->preop.cmat);
                if (mres == NULL) {
                        fprintf(stderr, "%s: matrix transform failed\n",
                                                rop->inspec);
                        goto failure;
                }
                if (verbose)
                        fprintf(stderr, "%s: applied %d x %d transform\n",
                                        rop->inspec, mres->ncomp,
                                        rop->mtx->ncomp);
                rmx_free(rop->mtx);
                rop->mtx = mres;
        }
        if (rop->preop.transpose) {     /* transpose matrix? */
                mres = rmx_transpose(rop->mtx);
                if (mres == NULL) {
                        fputs(rop->inspec, stderr);
                        fputs(": transpose failed\n", stderr);
                        goto failure;
                }
                if (verbose) {
                        fputs(rop->inspec, stderr);
                        fputs(": transposed rows and columns\n", stderr);
                }
                rmx_free(rop->mtx);
                rop->mtx = mres;
        }
        mres = rop->mtx;
        rop->mtx = NULL;
        return(mres);
failure:
        rmx_free(rop->mtx);
        return(rop->mtx = NULL);
}

/* Execute binary operation, free matrix arguments and return new result */
static RMATRIX *
binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
{
        RMATRIX *mres = NULL;
        int     i;

        if ((mleft == NULL) | (mright == NULL))
                return(NULL);
        switch (op) {
        case '.':                       /* concatenate */
                if (mleft->ncomp != mright->ncomp) {
                        fputs(inspec, stderr);
                        fputs(": # components do not match\n", stderr);
                } else if (mleft->ncols != mright->nrows) {
                        fputs(inspec, stderr);
                        fputs(": mismatched dimensions\n",
                                        stderr);
                } else
                        mres = rmx_multiply(mleft, mright);
                rmx_free(mleft);
                rmx_free(mright);
                if (mres == NULL) {
                        fputs(inspec, stderr);
                        fputs(": concatenation failed\n", stderr);
                        return(NULL);
                }
                if (verbose) {
                        fputs(inspec, stderr);
                        fputs(": concatenated matrix\n", stderr);
                }
                break;
        case '+':
                if (!rmx_sum(mleft, mright, NULL)) {
                        fputs(inspec, stderr);
                        fputs(": matrix sum failed\n", stderr);
                        rmx_free(mleft);
                        rmx_free(mright);
                        return(NULL);
                }
                if (verbose) {
                        fputs(inspec, stderr);
                        fputs(": added in matrix\n", stderr);
                }
                rmx_free(mright);
                mres = mleft;
                break;
        case '*':
        case '/': {
                const char *    tnam = (op == '/') ?
                                        "division" : "multiplication";
                errno = 0;
                if (!rmx_elemult(mleft, mright, (op == '/'))) {
                        fprintf(stderr, "%s: element-wise %s failed\n",
                                        inspec, tnam);
                        rmx_free(mleft);
                        rmx_free(mright);
                        return(NULL);
                }
                if (errno)
                        fprintf(stderr,
                                "%s: warning - error during element-wise %s\n",
                                        inspec, tnam);
                else if (verbose)
                        fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
                rmx_free(mright);
                mres = mleft;
                } break;
        default:
                fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
                rmx_free(mleft);
                rmx_free(mright);
                return(NULL);
        }
        return(mres);
}

/* Perform matrix operations from left to right */
static RMATRIX *
op_left2right(ROPMAT *mop)
{
        RMATRIX *mleft = loadop(mop);

        while (mop->binop) {
                if (mleft == NULL)
                        break;
                mleft = binaryop(mop[1].inspec,
                                mleft, mop->binop, loadop(mop+1));
                mop++;
        }
        return(mleft);
}

/* Perform matrix operations from right to left */
static RMATRIX *
op_right2left(ROPMAT *mop)
{
        RMATRIX *mright;
        int     rpos = 0;
                                        /* find end of list */
        while (mop[rpos].binop)
                if (mop[rpos++].binop != '.') {
                        fputs(
                "Right-to-left evaluation only for matrix multiplication!\n",
                                        stderr);
                        return(NULL);
                }
        mright = loadop(mop+rpos);
        while (rpos-- > 0) {
                if (mright == NULL)
                        break;
                mright = binaryop(mop[rpos].inspec,
                                loadop(mop+rpos), mop[rpos].binop, mright);
        }
        return(mright);
}

#define t_nrows(mop)    ((mop)->preop.transpose ? (mop)->mtx->ncols \
                                                : (mop)->mtx->nrows)
#define t_ncols(mop)    ((mop)->preop.transpose ? (mop)->mtx->nrows \
                                                : (mop)->mtx->ncols)

/* Should we prefer concatenating from rightmost matrix towards left? */
static int
prefer_right2left(ROPMAT *mop)
{
        int     mri = 0;

        while (mop[mri].binop)          /* find rightmost matrix */
                if (mop[mri++].binop != '.')
                        return(0);      /* pre-empt reversal for other ops */

        if (mri <= 1)
                return(0);              /* won't matter */

        if (loadmatrix(mop+mri) < 0)    /* load rightmost cat */
                return(1);              /* fail will bail in a moment */

        if (t_ncols(mop+mri) == 1)
                return(1);              /* definitely better R->L */

        if (t_ncols(mop+mri) >= t_nrows(mop+mri))
                return(0);              /* ...probably worse */

        if (loadmatrix(mop) < 0)        /* load leftmost */
                return(0);              /* fail will bail in a moment */

        return(t_ncols(mop+mri) < t_nrows(mop));
}

static int
get_factors(double da[], int n, char *av[])
{
        int     ac;

        for (ac = 0; ac < n && isflt(av[ac]); ac++)
                da[ac] = atof(av[ac]);
        return(ac);
}

static ROPMAT *
grow_moparray(ROPMAT *mop, int n2alloc)
{
        int     nmats = 0;

        while (mop[nmats++].binop)
                ;
        mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
        if (mop == NULL) {
                fputs("Out of memory in grow_moparray()\n", stderr);
                exit(1);
        }
        if (n2alloc > nmats)
                memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
        return(mop);
}

/* Load one or more matrices and operate on them, sending results to stdout */
int
main(int argc, char *argv[])
{
        int     outfmt = DTfromHeader;
        int     nall = 2;
        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
        int     nmats = 0;
        RMATRIX *mres = NULL;
        int     stdin_used = 0;
        int     i;
                                        /* get options and arguments */
        for (i = 1; i < argc; i++) {
                if (argv[i][0] && !argv[i][1] &&
                                strchr(".+*/", argv[i][0]) != NULL) {
                        if (!nmats || mop[nmats-1].binop) {
                                fprintf(stderr,
                        "%s: missing matrix argument before '%c' operation\n",
                                                argv[0], argv[i][0]);
                                return(1);
                        }
                        mop[nmats-1].binop = argv[i][0];
                } else if (argv[i][0] != '-' || !argv[i][1]) {
                        if (argv[i][0] == '-') {
                                if (stdin_used++) {
                                        fprintf(stderr,
                        "%s: standard input used for more than one matrix\n",
                                                argv[0]);
                                        return(1);
                                }
                                mop[nmats].inspec = stdin_name;
                        } else
                                mop[nmats].inspec = argv[i];
                        if (nmats > 0 && !mop[nmats-1].binop)
                                mop[nmats-1].binop = '.';
                        nmats++;
                } else {
                        int     n = argc-1 - i;
                        switch (argv[i][1]) {   /* get option */
                        case 'v':
                                verbose++;
                                break;
                        case 'f':
                                switch (argv[i][2]) {
                                case 'd':
                                        outfmt = DTdouble;
                                        break;
                                case 'f':
                                        outfmt = DTfloat;
                                        break;
                                case 'a':
                                        outfmt = DTascii;
                                        break;
                                case 'c':
                                        outfmt = DTrgbe;
                                        break;
                                default:
                                        goto userr;
                                }
                                break;
                        case 't':
                                mop[nmats].preop.transpose = 1;
                                break;
                        case 's':
                                if (n > MAXCOMP) n = MAXCOMP;
                                i += mop[nmats].preop.nsf =
                                        get_factors(mop[nmats].preop.sca,
                                                        n, argv+i+1);
                                break;
                        case 'c':
                                if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
                                i += mop[nmats].preop.clen =
                                        get_factors(mop[nmats].preop.cmat,
                                                        n, argv+i+1);
                                break;
                        default:
                                fprintf(stderr, "%s: unknown operation '%s'\n",
                                                argv[0], argv[i]);
                                goto userr;
                        }
                }
                if (nmats >= nall)
                        mop = grow_moparray(mop, nall += 2);
        }
        if (mop[0].inspec == NULL)      /* nothing to do? */
                goto userr;
        if (mop[nmats-1].binop) {
                fprintf(stderr,
                        "%s: missing matrix argument after '%c' operation\n",
                                argv[0], mop[nmats-1].binop);
                return(1);
        }
                                        /* favor quicker concatenation */
        mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
                                                : op_left2right(mop);
        if (mop[nmats].mtx == NULL)
                return(1);
                                        /* apply trailing unary operations */
        mop[nmats].inspec = "trailing_ops";
        mres = loadop(mop+nmats);
        if (mres == NULL)
                return(1);
                                        /* write result to stdout */
        if (outfmt == DTfromHeader)
                outfmt = mres->dtype;
        if (outfmt != DTascii)
                SET_FILE_BINARY(stdout);
        newheader("RADIANCE", stdout);
        printargs(argc, argv, stdout);
        if (!rmx_write(mres, outfmt, stdout)) {
                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
                return(1);
        }
        /* rmx_free(mres); free(mop); */
        return(0);
userr:
        fprintf(stderr,
        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..] m1 [.+*/] .. > mres\n",
                        argv[0]);
        return(1);
}
Revision:	2.16
Committed:	Mon Aug 12 20:38:19 2019 UTC (6 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.15:	+11 -8 lines
Log Message:	Added warning about divide-by-zero to matrix information header
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.16	static const char RCSid[] = "$Id: rmtxop.c,v 2.15 2019/08/12 16:55:24 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[MAXCOMPMAXCOMP]; / 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			switch (op) {
144			case '.': /* concatenate */
145	greg	2.16	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	greg	2.13	rmx_free(mleft);
155	greg	2.12	rmx_free(mright);
156	greg	2.1	if (mres == NULL) {
157	greg	2.13	fputs(inspec, stderr);
158	greg	2.16	fputs(": concatenation failed\n", stderr);
159	greg	2.1	return(NULL);
160			}
161			if (verbose) {
162	greg	2.13	fputs(inspec, stderr);
163	greg	2.1	fputs(": concatenated matrix\n", stderr);
164			}
165	greg	2.13	break;
166			case '+':
167			if (!rmx_sum(mleft, mright, NULL)) {
168			fputs(inspec, stderr);
169	greg	2.1	fputs(": matrix sum failed\n", stderr);
170	greg	2.13	rmx_free(mleft);
171	greg	2.1	rmx_free(mright);
172			return(NULL);
173			}
174			if (verbose) {
175	greg	2.13	fputs(inspec, stderr);
176	greg	2.1	fputs(": added in matrix\n", stderr);
177			}
178			rmx_free(mright);
179	greg	2.13	mres = mleft;
180			break;
181			case '*':
182			case '/': {
183			const char * tnam = (op == '/') ?
184	greg	2.11	"division" : "multiplication";
185			errno = 0;
186	greg	2.13	if (!rmx_elemult(mleft, mright, (op == '/'))) {
187	greg	2.11	fprintf(stderr, "%s: element-wise %s failed\n",
188	greg	2.13	inspec, tnam);
189			rmx_free(mleft);
190	greg	2.11	rmx_free(mright);
191			return(NULL);
192			}
193			if (errno)
194			fprintf(stderr,
195			"%s: warning - error during element-wise %s\n",
196	greg	2.13	inspec, tnam);
197	greg	2.11	else if (verbose)
198	greg	2.13	fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
199	greg	2.11	rmx_free(mright);
200	greg	2.13	mres = mleft;
201			} break;
202			default:
203			fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
204			rmx_free(mleft);
205	greg	2.1	rmx_free(mright);
206			return(NULL);
207			}
208	greg	2.13	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	greg	2.1	return(mleft);
225			}
226
227	greg	2.13	/* 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	greg	2.14	if (mop[rpos++].binop != '.') {
236			fputs(
237			"Right-to-left evaluation only for matrix multiplication!\n",
238			stderr);
239			return(NULL);
240			}
241	greg	2.13	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	greg	2.1	static int
285			get_factors(double da[], int n, char *av[])
286			{
287			int ac;
288
289			for (ac = 0; ac < n && isflt(av[ac]); ac++)
290			da[ac] = atof(av[ac]);
291			return(ac);
292			}
293
294	greg	2.13	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, n2allocsizeof(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	greg	2.1	/* Load one or more matrices and operate on them, sending results to stdout */
312			int
313			main(int argc, char *argv[])
314			{
315	greg	2.7	int outfmt = DTfromHeader;
316	greg	2.13	int nall = 2;
317			ROPMAT mop = (ROPMAT )calloc(nall, sizeof(ROPMAT));
318			int nmats = 0;
319	greg	2.1	RMATRIX *mres = NULL;
320	greg	2.13	int stdin_used = 0;
321	greg	2.1	int i;
322			/* get options and arguments */
323	greg	2.13	for (i = 1; i < argc; i++) {
324	greg	2.11	if (argv[i][0] && !argv[i][1] &&
325	greg	2.13	strchr(".+*/", argv[i][0]) != NULL) {
326	greg	2.15	if (!nmats \|\| mop[nmats-1].binop) {
327	greg	2.14	fprintf(stderr,
328	greg	2.15	"%s: missing matrix argument before '%c' operation\n",
329	greg	2.14	argv[0], argv[i][0]);
330	greg	2.13	return(1);
331			}
332	greg	2.15	mop[nmats-1].binop = argv[i][0];
333	greg	2.1	} else if (argv[i][0] != '-' \|\| !argv[i][1]) {
334	greg	2.13	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	greg	2.1	} else {
348			int n = argc-1 - i;
349			switch (argv[i][1]) { /* get option */
350			case 'v':
351	greg	2.14	verbose++;
352	greg	2.1	break;
353			case 'f':
354			switch (argv[i][2]) {
355			case 'd':
356			outfmt = DTdouble;
357			break;
358			case 'f':
359			outfmt = DTfloat;
360			break;
361			case 'a':
362			outfmt = DTascii;
363			break;
364			case 'c':
365			outfmt = DTrgbe;
366			break;
367			default:
368			goto userr;
369			}
370			break;
371			case 't':
372	greg	2.13	mop[nmats].preop.transpose = 1;
373	greg	2.1	break;
374			case 's':
375			if (n > MAXCOMP) n = MAXCOMP;
376	greg	2.13	i += mop[nmats].preop.nsf =
377			get_factors(mop[nmats].preop.sca,
378			n, argv+i+1);
379	greg	2.1	break;
380			case 'c':
381			if (n > MAXCOMPMAXCOMP) n = MAXCOMPMAXCOMP;
382	greg	2.13	i += mop[nmats].preop.clen =
383			get_factors(mop[nmats].preop.cmat,
384			n, argv+i+1);
385	greg	2.1	break;
386			default:
387			fprintf(stderr, "%s: unknown operation '%s'\n",
388			argv[0], argv[i]);
389			goto userr;
390			}
391			}
392	greg	2.14	if (nmats >= nall)
393			mop = grow_moparray(mop, nall += 2);
394	greg	2.13	}
395			if (mop[0].inspec == NULL) /* nothing to do? */
396	greg	2.1	goto userr;
397	greg	2.15	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	greg	2.13	/* favor quicker concatenation */
404	greg	2.14	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	greg	2.13	return(1);
413	greg	2.1	/* write result to stdout */
414	greg	2.6	if (outfmt == DTfromHeader)
415			outfmt = mres->dtype;
416	greg	2.4	if (outfmt != DTascii)
417	greg	2.10	SET_FILE_BINARY(stdout);
418	greg	2.1	newheader("RADIANCE", stdout);
419			printargs(argc, argv, stdout);
420	greg	2.5	if (!rmx_write(mres, outfmt, stdout)) {
421	greg	2.1	fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
422			return(1);
423			}
424	greg	2.13	/* rmx_free(mres); free(mop); */
425	greg	2.1	return(0);
426			userr:
427			fprintf(stderr,
428	greg	2.13	"Usage: %s [-v][-f[adfc][-t][-s sf .. \| -c ce ..] m1 [.+*/] .. > mres\n",
429	greg	2.1	argv[0]);
430			return(1);
431			}