ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.12
Committed: Mon Aug 27 23:03:05 2018 UTC (5 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2
Changes since 2.11: +16 -16 lines
Log Message:
Moved transpose operation later to save memory with -c (transform) reduction

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmtxop.c,v 2.11 2017/08/28 15:59:46 greg Exp $";
3 #endif
4 /*
5 * General component matrix operations.
6 */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <errno.h>
11 #include "rtio.h"
12 #include "resolu.h"
13 #include "rmatrix.h"
14 #include "platform.h"
15
16 #define MAXCOMP 50 /* #components we support */
17
18 typedef struct {
19 double sca[MAXCOMP]; /* scalar coefficients */
20 int nsf; /* number of scalars */
21 double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
22 int clen; /* number of coefficients */
23 int transpose; /* do transpose? */
24 int op; /* '*' or '+' */
25 } ROPERAT; /* matrix operation */
26
27 int verbose = 0; /* verbose reporting? */
28
29 static void
30 op_default(ROPERAT *op)
31 {
32 memset(op, 0, sizeof(ROPERAT));
33 op->op = '.';
34 }
35
36 static RMATRIX *
37 operate(RMATRIX *mleft, ROPERAT *op, const char *fname)
38 {
39 RMATRIX *mright = rmx_load(fname);
40 RMATRIX *mtmp;
41 int i;
42
43 if (fname == NULL)
44 fname = "<stdin>";
45 if (mright == NULL) {
46 fputs(fname, stderr);
47 fputs(": cannot load matrix\n", stderr);
48 return(NULL);
49 }
50 if (op->nsf > 0) { /* apply scalar(s) */
51 if (op->clen > 0) {
52 fputs("Options -s and -c are exclusive\n", stderr);
53 rmx_free(mright);
54 return(NULL);
55 }
56 if (op->nsf == 1) {
57 for (i = mright->ncomp; --i; )
58 op->sca[i] = op->sca[0];
59 } else if (op->nsf != mright->ncomp) {
60 fprintf(stderr, "%s: -s must have one or %d factors\n",
61 fname, mright->ncomp);
62 rmx_free(mright);
63 return(NULL);
64 }
65 if ((mleft == NULL) | (op->op != '+') &&
66 !rmx_scale(mright, op->sca)) {
67 fputs(fname, stderr);
68 fputs(": scalar operation failed\n", stderr);
69 rmx_free(mright);
70 return(NULL);
71 }
72 if (verbose) {
73 fputs(fname, stderr);
74 fputs(": applied scalar (", stderr);
75 for (i = 0; i < op->nsf; i++)
76 fprintf(stderr, " %f", op->sca[i]);
77 fputs(" )\n", stderr);
78 }
79 }
80 if (op->clen > 0) { /* apply transform */
81 if (op->clen % mright->ncomp) {
82 fprintf(stderr, "%s: -c must have N x %d coefficients\n",
83 fname, mright->ncomp);
84 rmx_free(mright);
85 return(NULL);
86 }
87 mtmp = rmx_transform(mright, op->clen/mright->ncomp, op->cmat);
88 if (mtmp == NULL) {
89 fprintf(stderr, "%s: matrix transform failed\n", fname);
90 rmx_free(mright);
91 return(NULL);
92 }
93 if (verbose)
94 fprintf(stderr, "%s: applied %d x %d transform\n",
95 fname, mtmp->ncomp, mright->ncomp);
96 rmx_free(mright);
97 mright = mtmp;
98 }
99 if (op->transpose) { /* transpose matrix? */
100 mtmp = rmx_transpose(mright);
101 if (mtmp == NULL) {
102 fputs(fname, stderr);
103 fputs(": transpose failed\n", stderr);
104 rmx_free(mright);
105 return(NULL);
106 }
107 if (verbose) {
108 fputs(fname, stderr);
109 fputs(": transposed rows and columns\n", stderr);
110 }
111 rmx_free(mright);
112 mright = mtmp;
113 }
114 if (mleft == NULL) /* just one matrix */
115 return(mright);
116 if (op->op == '.') { /* concatenate */
117 RMATRIX *mres = rmx_multiply(mleft, mright);
118 if (mres == NULL) {
119 fputs(fname, stderr);
120 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);
126 return(NULL);
127 }
128 if (verbose) {
129 fputs(fname, stderr);
130 fputs(": concatenated matrix\n", stderr);
131 }
132 rmx_free(mright);
133 rmx_free(mleft);
134 mleft = mres;
135 } else if (op->op == '+') {
136 if (!rmx_sum(mleft, mright, op->nsf ? op->sca : (double *)NULL)) {
137 fputs(fname, stderr);
138 fputs(": matrix sum failed\n", stderr);
139 rmx_free(mright);
140 return(NULL);
141 }
142 if (verbose) {
143 fputs(fname, stderr);
144 fputs(": added in matrix\n", stderr);
145 }
146 rmx_free(mright);
147 } else if ((op->op == '*') | (op->op == '/')) {
148 const char * tnam = (op->op == '/') ?
149 "division" : "multiplication";
150 errno = 0;
151 if (!rmx_elemult(mleft, mright, (op->op == '/'))) {
152 fprintf(stderr, "%s: element-wise %s failed\n",
153 fname, tnam);
154 rmx_free(mright);
155 return(NULL);
156 }
157 if (errno)
158 fprintf(stderr,
159 "%s: warning - error during element-wise %s\n",
160 fname, tnam);
161 else if (verbose)
162 fprintf(stderr, "%s: element-wise %s\n", fname, tnam);
163 rmx_free(mright);
164 } else {
165 fprintf(stderr, "%s: unknown operation '%c'\n", fname, op->op);
166 rmx_free(mright);
167 return(NULL);
168 }
169 return(mleft);
170 }
171
172 static int
173 get_factors(double da[], int n, char *av[])
174 {
175 int ac;
176
177 for (ac = 0; ac < n && isflt(av[ac]); ac++)
178 da[ac] = atof(av[ac]);
179 return(ac);
180 }
181
182 /* Load one or more matrices and operate on them, sending results to stdout */
183 int
184 main(int argc, char *argv[])
185 {
186 int outfmt = DTfromHeader;
187 RMATRIX *mres = NULL;
188 ROPERAT op;
189 int i;
190 /* initialize */
191 op_default(&op);
192 /* get options and arguments */
193 for (i = 1; i < argc; i++)
194 if (argv[i][0] && !argv[i][1] &&
195 strchr("+*/", argv[i][0]) != NULL) {
196 op.op = argv[i][0];
197 } else if (argv[i][0] != '-' || !argv[i][1]) {
198 char *fname = NULL; /* load matrix */
199 if (argv[i][0] != '-')
200 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);
206 }
207 op_default(&op); /* reset operator */
208 } else {
209 int n = argc-1 - i;
210 switch (argv[i][1]) { /* get option */
211 case 'v':
212 verbose = !verbose;
213 break;
214 case 'f':
215 switch (argv[i][2]) {
216 case 'd':
217 outfmt = DTdouble;
218 break;
219 case 'f':
220 outfmt = DTfloat;
221 break;
222 case 'a':
223 outfmt = DTascii;
224 break;
225 case 'c':
226 outfmt = DTrgbe;
227 break;
228 default:
229 goto userr;
230 }
231 break;
232 case 't':
233 op.transpose = 1;
234 break;
235 case 's':
236 if (n > MAXCOMP) n = MAXCOMP;
237 op.nsf = get_factors(op.sca, n, argv+i+1);
238 i += op.nsf;
239 break;
240 case 'c':
241 if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
242 op.clen = get_factors(op.cmat, n, argv+i+1);
243 i += op.clen;
244 break;
245 default:
246 fprintf(stderr, "%s: unknown operation '%s'\n",
247 argv[0], argv[i]);
248 goto userr;
249 }
250 }
251 if (mres == NULL) /* check that we got something */
252 goto userr;
253 /* write result to stdout */
254 if (outfmt == DTfromHeader)
255 outfmt = mres->dtype;
256 if (outfmt != DTascii)
257 SET_FILE_BINARY(stdout);
258 newheader("RADIANCE", stdout);
259 printargs(argc, argv, stdout);
260 if (!rmx_write(mres, outfmt, stdout)) {
261 fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
262 return(1);
263 }
264 /* rmx_free(mres); mres = NULL; */
265 return(0);
266 userr:
267 fprintf(stderr,
268 "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..] m1 [+*/] .. > mres\n",
269 argv[0]);
270 return(1);
271 }