--- ray/src/util/rcomb.c	2023/12/19 20:38:38	2.5
+++ ray/src/util/rcomb.c	2025/04/04 18:06:48	2.29
@@ -1,25 +1,26 @@
 #ifndef lint
-static const char RCSid[] = "$Id: rcomb.c,v 2.5 2023/12/19 20:38:38 greg Exp $";
+static const char RCSid[] = "$Id: rcomb.c,v 2.29 2025/04/04 18:06:48 greg Exp $";
 #endif
 /*
  * General component matrix combiner, operating on a row at a time.
+ *
+ * Multi-processing mode under Unix creates children that each work
+ * on one input row at a time, fed by the original process.  Final conversion
+ * and output to stdout is sorted by last child while its siblings send it
+ * their record calculations.
  */
 
-#include <errno.h>
 #include <math.h>
 #include "platform.h"
+#include "rtprocess.h"
 #include "rtio.h"
-#include "resolu.h"
 #include "rmatrix.h"
 #include "calcomp.h"
-#include "paths.h"
 
 #ifndef M_PI
 #define M_PI	3.14159265358979323846
 #endif
 
-#define MAXCOMP		MAXCSAMP	/* #components we support */
-
 /* Unary matrix operation(s) */
 typedef struct {
 	double		cmat[MAXCOMP*MAXCOMP];	/* component transformation */
@@ -55,9 +56,13 @@ int		cur_row;			/* current input/output row */
 int		cur_col;			/* current input/output column */
 int		cur_chan;			/* if we're looping channels */
 
-static int	checksymbolic(ROPMAT *rop);
+SUBPROC		*cproc = NULL;			/* child process array */
+int		nchildren = 0;			/* # of child processes */
+int		inchild = -1;			/* our child ID (-1: parent) */
 
-static int
+extern int	checksymbolic(ROPMAT *rop);
+
+int
 split_input(ROPMAT *rop)
 {
 	if (rop->rmp == &rop->imx && !(rop->rmp = rmx_copy(&rop->imx))) {
@@ -69,7 +74,7 @@ split_input(ROPMAT *rop)
 }
 
 /* Check/set transform based on a reference input file */
-static int
+int
 checkreffile(ROPMAT *rop)
 {
 	static const char	*curRF = NULL;
@@ -132,7 +137,7 @@ checkreffile(ROPMAT *rop)
 }
 
 /* Compute conversion row from spectrum to one channel of RGB */
-static void
+void
 rgbrow(ROPMAT *rop, int r, int p)
 {
 	const int	nc = rop->imx.ncomp;
@@ -149,7 +154,7 @@ rgbrow(ROPMAT *rop, int r, int p)
 }
 
 /* Compute conversion row from spectrum to one channel of XYZ */
-static void
+void
 xyzrow(ROPMAT *rop, int r, int p)
 {
 	const int	nc = rop->imx.ncomp;
@@ -166,8 +171,8 @@ xyzrow(ROPMAT *rop, int r, int p)
 }
 
 /* Use the spectral sensitivity function to compute matrix coefficients */
-static void
-sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
+void
+sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
 {
 	const int	nc = rop->imx.ncomp;
 	int		i;
@@ -181,7 +186,7 @@ sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, in
 }
 
 /* Check/set symbolic transform */
-static int
+int
 checksymbolic(ROPMAT *rop)
 {
 	const int	nc = rop->imx.ncomp;
@@ -295,7 +300,7 @@ checksymbolic(ROPMAT *rop)
 	return(1);
 }
 
-static int
+int
 get_component_xfm(ROPMAT *rop)
 {
 	int	i, j;
@@ -378,7 +383,7 @@ get_component_xfm(ROPMAT *rop)
 	return(1);
 }
 
-static int
+int
 apply_op(RMATRIX *dst, const RMATRIX *src, const RUNARYOP *ro)
 {
 	if (ro->clen > 0) {
@@ -391,14 +396,13 @@ apply_op(RMATRIX *dst, const RMATRIX *src, const RUNAR
 			return(0);
 		rmx_free(res);
 	} else if (dst != src)
-		memcpy(dst->mtx, src->mtx,
-				sizeof(double)*dst->ncomp*dst->ncols*dst->nrows);
+		memcpy(dst->mtx, src->mtx, rmx_array_size(dst));
 	if (ro->nsf == dst->ncomp)
 		rmx_scale(dst, ro->sca);
 	return(1);
 }
 
-static int
+int
 open_input(ROPMAT *rop)
 {
 	int	outtype;
@@ -420,7 +424,7 @@ open_input(ROPMAT *rop)
 }
 
 /* Return nominal wavelength associated with input component (return nm) */
-static double
+double
 l_wavelength(char *nam)
 {
 	double	comp = argument(1);
@@ -444,7 +448,7 @@ l_wavelength(char *nam)
 }
 
 /* Return ith input with optional channel selector */
-static double
+double
 l_chanin(char *nam)
 {
 	double	inp = argument(1);
@@ -469,7 +473,7 @@ l_chanin(char *nam)
 	return(mop[mi].rmp->mtx[cur_col*in_ncomp + chan]);
 }
 
-static int
+int
 initialize(RMATRIX *imp)
 {
 	int	i;
@@ -482,7 +486,7 @@ initialize(RMATRIX *imp)
 		restype = mop[i].rmp->dtype;
 		if (!imp->dtype || (restype = rmx_newtype(restype, imp->dtype)) > 0)
 			imp->dtype = restype;
-		else
+		else if (!nowarn)
 			fprintf(stderr, "%s: warning - data type mismatch\n",
 					mop[i].inspec);
 		if (!i) {
@@ -516,7 +520,7 @@ initialize(RMATRIX *imp)
 	return(1);
 }
 
-static void
+void
 output_headinfo(FILE *fp)
 {
 	int	i;
@@ -539,16 +543,23 @@ output_headinfo(FILE *fp)
 	}
 }
 
-static int
-combine_input(ROPMAT *res, FILE *fout)
+int
+spawned_children(int np)
 {
-	int	set_r, set_c;
-	RMATRIX	*tmp = NULL;
-	int	co_set;
-	int	i;
-					/* allocate input row buffers */
+	int	i, rv;
+
+#if defined(_WIN32) || defined(_WIN64)
+	if (np > 1) {
+		if (!nowarn)
+			fputs("Warning: only one process under Windows\n", stderr);
+		np = 1;
+	} else
+#endif
+	if ((in_nrows > 0) & (np*4 > in_nrows))
+		np = in_nrows/4;
+				/* we'll be doing a row at a time */
 	for (i = 0; i < nmats; i++) {
-		mop[i].imx.nrows = 1;	/* we'll be doing a row at a time */
+		mop[i].imx.nrows = 1;
 		if (!rmx_prepare(&mop[i].imx))
 			goto memerror;
 		if (mop[i].rmp != &mop[i].imx) {
@@ -557,22 +568,152 @@ combine_input(ROPMAT *res, FILE *fout)
 				goto memerror;
 		}
 	}
-					/* prep output row buffers */
-	if (mcat || res->preop.clen > 0) {
-		if (!split_input(res))	/* need separate buffer */
-			return(0);
-		if (res->preop.clen > 0)
-			res->rmp->ncomp = res->preop.clen / res->imx.ncomp;
-		res->rmp->nrows = 1;
-		if (!mcat | !mcat_last && !rmx_prepare(res->rmp))
+				/* prep output row buffer(s) */
+	if (mop[nmats].preop.clen > 0) {
+		if (!split_input(&mop[nmats]))	/* need separate buffer */
 			goto memerror;
+		mop[nmats].rmp->ncomp = mop[nmats].preop.clen /
+					mop[nmats].imx.ncomp;
 	}
-	if (mcat && mcat_last &&
-			!(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp)))
+	mop[nmats].imx.nrows = 1;
+	if (!rmx_prepare(&mop[nmats].imx))
 		goto memerror;
-	res->imx.nrows = 1;
-	if (!rmx_prepare(&res->imx))
+	if (mop[nmats].rmp != &mop[nmats].imx) {
+		mop[nmats].rmp->nrows = 1;
+		if (!rmx_prepare(mop[nmats].rmp))
+			goto memerror;
+	}
+	if (np <= 1) {		/* single process return */
+#ifdef getc_unlocked
+		for (i = 0; i < nmats; i++)
+			flockfile(mop[i].infp);
+		flockfile(stdout);
+#endif
+		return(0);
+	}
+	fflush(stdout);		/* flush header & spawn children */
+	nchildren = np + 1;	/* extra child to sequence output */
+	cproc = (SUBPROC *)malloc(sizeof(SUBPROC)*nchildren);
+	if (!cproc)
 		goto memerror;
+	for (i = nchildren; i--; ) cproc[i] = sp_inactive;
+	cproc[nchildren-1].flags |= PF_FILT_OUT;
+				/* start each child from parent */
+	for (i = 0; i < nchildren; i++)
+		if ((rv = open_process(&cproc[i], NULL)) <= 0)
+			break;	/* child breaks here */
+	if (rv < 0) {
+		perror("fork");	/* WTH? */
+		close_processes(cproc, i);
+		exit(1);
+	}
+	if (i != nchildren-1) {	/* last child is sole reader */
+		int	j = i;
+		while (j-- > 0) {
+			close(cproc[j].r);
+			cproc[j].r = -1;
+		}
+	}
+	if (rv > 0)
+		return(1);	/* parent return value */
+
+	inchild = i;		/* else set our child index */
+	while (i-- > 0)		/* only parent writes siblings */
+		close(cproc[i].w);
+
+	i = nmats;		/* close matrix streams (carefully) */
+	while (i-- > 0) {
+		if (mop[i].infp != stdin) {
+			close(fileno(mop[i].infp));	/* avoid lseek() */
+			fclose(mop[i].infp);		/* ! pclose() */
+		}
+		mop[i].infp = NULL;
+	}
+	fpurge(stdin);		/* discard previously buffered input */
+
+	if (inchild == nchildren-1)
+		return(-1);	/* output process return value */
+
+	i = nmats;		/* get matrix rows from parent */
+	while (i-- > 0) {
+		mop[i].infp = stdin;
+		mop[i].imx.dtype = DTrmx_native;
+		mop[i].imx.pflags &= ~RMF_SWAPIN;
+	}
+#ifdef getc_unlocked
+	flockfile(stdin);
+#endif
+	mop[nmats].rmp->dtype = DTrmx_native;
+	return(0);		/* worker child return value */
+memerror:
+	fputs("Out of memory in spawned_children()\n", stderr);
+	exit(1);
+}
+
+int
+parent_loop(void)
+{
+	int	i;
+
+	rmx_reset(&mop[nmats].imx);		/* not touching output side */
+	if (mop[nmats].rmp != &mop[nmats].imx) {
+		rmx_free(mop[nmats].rmp);
+		mop[nmats].rmp = &mop[nmats].imx;
+	}
+#ifdef getc_unlocked
+	for (i = 0; i < nmats; i++)		/* we handle matrix inputs */
+		flockfile(mop[i].infp);
+#endif
+						/* load & send rows to kids */
+	for (cur_row = 0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row++) {
+	    int		wfd = cproc[cur_row % (nchildren-1)].w;
+	    for (i = 0; i < nmats; i++)
+		if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
+			if (cur_row > in_nrows)	/* unknown #input rows? */
+				break;
+			fprintf(stderr, "%s: load error at row %d\n",
+					mop[i].inspec, cur_row);
+			return(0);
+		}
+	    if (i < nmats)
+		break;
+	    for (i = 0; i < nmats; i++)
+	    	if (writebuf(wfd, mop[i].imx.mtx, rmx_array_size(&mop[i].imx))
+	    				!= rmx_array_size(&mop[i].imx)) {
+			fprintf(stderr, "%s: write error at row %d\n",
+					mop[i].inspec, cur_row);
+			return(0);
+		}
+	}
+	i = close_processes(cproc, nchildren);	/* collect family */
+	free(cproc); cproc = NULL; nchildren = 0;
+	if (i < 0) {
+		if (!nowarn)
+			fputs("Warning: lost child process\n", stderr);
+		return(1);
+	}
+	if (i > 0) {
+		fprintf(stderr, "Child exited with status %d\n", i);
+		return(0);
+	}
+	return(1);				/* return success! */
+}
+
+int
+combine_input(void)
+{
+	const int	row0 = (inchild >= 0)*inchild;
+	const int	rstep = nchildren ? nchildren-1 : 1;
+	ROPMAT		*res = &mop[nmats];
+	int		set_r, set_c;
+	RMATRIX		*tmp = NULL;
+	int		co_set;
+	int		i;
+
+	if (mcat_last && !(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp))) {
+		fputs("Out of buffer space in combine_input()\n", stderr);
+		return(0);
+	}
 					/* figure out what the user set */
 	co_set = fundefined("co");
 	if (!co_set)
@@ -588,19 +729,21 @@ combine_input(ROPMAT *res, FILE *fout)
 	} else				/* save a little time */
 		set_r = set_c = 0;
 					/* read/process row-by-row */
-	for (cur_row = 0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row++) {
+	for (cur_row = row0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row += rstep) {
 	    RMATRIX	*mres = NULL;
-	    for (i = 0; i < nmats; i++) {
+	    for (i = 0; i < nmats; i++)
 		if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
 			if (cur_row > in_nrows)	/* unknown #input rows? */
-				goto loop_exit;
-			fprintf(stderr, "%s: read error at row %d\n",
+				break;
+			fprintf(stderr, "%s: load error at row %d\n",
 					mop[i].inspec, cur_row);
 			return(0);
 		}
+	    if (i < nmats)
+	    	break;
+	    for (i = 0; i < nmats; i++)
 		if (!apply_op(mop[i].rmp, &mop[i].imx, &mop[i].preop))
 			return(0);
-	    }
 	    if (set_r) varset("r", '=', cur_row);
 	    for (cur_col = 0; cur_col < in_ncols; cur_col++) {
 	    	if (set_c) varset("c", '=', cur_col);
@@ -638,34 +781,59 @@ combine_input(ROPMAT *res, FILE *fout)
 	    }
 	    rmx_free(mres); mres = NULL;
 	    if (!rmx_write_data(res->rmp->mtx, res->rmp->ncomp,
-	    			res->rmp->ncols, res->rmp->dtype, fout))
+	    			res->rmp->ncols, res->rmp->dtype, stdout) ||
+	    			 (inchild >= 0 && fflush(stdout) == EOF)) {
+		fprintf(stderr, "Conversion/write error at row %d\n",
+				cur_row);
 	    	return(0);
+	    }
 	}
-loop_exit:
-#if 0		/* we're about to exit, so who cares? */
-	rmx_free(tmp);			/* clean up */
-	rmx_reset(res->rmp);
-	rmx_reset(&res->imx);
-	for (i = 0; i < nmats; i++) {
-		rmx_reset(mop[i].rmp);
+	return(inchild >= 0 || fflush(stdout) != EOF);
+multerror:
+	fputs("Unexpected matrix multiply error\n", stderr);
+	return(0);
+}
+
+int
+output_loop(void)
+{
+	const size_t	row_size = rmx_array_size(mop[nmats].rmp);
+	int		cur_child = 0;
+	int		i = nmats;
+
+	while (i-- > 0) {				/* free input buffers */
 		rmx_reset(&mop[i].imx);
-		if (mop[i].inspec[0] == '!')
-			pclose(mop[i].infp);
-		else if (mop[i].inspec != stdin_name)
-			fclose(mop[i].infp);
-		mop[i].infp = NULL;
+		if (mop[i].rmp != &mop[i].imx) {
+			rmx_free(mop[i].rmp);
+			mop[i].rmp = &mop[i].imx;
+		}
 	}
+	if (mop[nmats].rmp != &mop[nmats].imx)		/* output is split? */
+		rmx_reset(&mop[nmats].imx);
+#ifdef getc_unlocked
+	flockfile(stdout);				/* we own this, now */
 #endif
-	return(fflush(fout) != EOF);
-memerror:
-	fputs("Out of buffer space in combine_input()\n", stderr);
-	return(0);
-multerror:
-	fputs("Unexpected matrix multiply error in combine_input()\n", stderr);
-	return(0);
+	for ( ; ; ) {					/* loop until no more */
+		ssize_t		rv;
+		rv = readbuf(cproc[cur_child].r, mop[nmats].rmp->mtx, row_size);
+		if (!rv)				/* out of rows? */
+			break;
+		if (rv != row_size) {
+			fputs("Read error\n", stderr);
+			return(0);
+		}					/* do final conversion */
+		if (!rmx_write_data(mop[nmats].rmp->mtx, mop[nmats].rmp->ncomp,
+	    			mop[nmats].rmp->ncols, mop[nmats].rmp->dtype, stdout)) {
+			fputs("Conversion/write error\n", stderr);
+			return(0);
+		}
+		cur_child++;
+		cur_child *= (cur_child < inchild);	/* loop over workers */
+	}
+	return(fflush(stdout) != EOF);
 }
 
-static int
+int
 get_factors(double da[], int n, char *av[])
 {
 	int	ac;
@@ -675,12 +843,14 @@ get_factors(double da[], int n, char *av[])
 	return(ac);
 }
 
-static void
+void
 resize_inparr(int n2alloc)
 {
 	int	i;
 
-	for (i = nmats; i > n2alloc; i--) {
+	if (n2alloc == nall)
+		return;
+	for (i = nall; i-- > n2alloc; ) {
 		rmx_reset(&mop[i].imx);
 		if (mop[i].rmp != &mop[i].imx)
 			rmx_free(mop[i].rmp);
@@ -690,8 +860,8 @@ resize_inparr(int n2alloc)
 		fputs("Out of memory in resize_inparr()\n", stderr);
 		exit(1);
 	}
-	if (n2alloc > nmats)
-		memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
+	if (n2alloc > nall)
+		memset(mop+nall, 0, (n2alloc-nall)*sizeof(ROPMAT));
 	nall = n2alloc;
 }
 
@@ -704,6 +874,7 @@ main(int argc, char *argv[])
 	const char	*defCsym = NULL;
 	int		echoheader = 1;
 	int		stdin_used = 0;
+	int		nproc = 1;
 	const char	*mcat_spec = NULL;
 	int		n2comp = 0;
 	uby8		comp_ndx[128];
@@ -731,6 +902,11 @@ main(int argc, char *argv[])
 			case 'h':
 				echoheader = !echoheader;
 				break;
+			case 'n':
+				nproc = atoi(argv[++i]);
+				if (nproc <= 0)
+					goto userr;
+				break;
 			case 'e':
 				if (!n) goto userr;
 				comp_ndx[n2comp++] = i++;
@@ -769,13 +945,14 @@ main(int argc, char *argv[])
 				}
 				break;
 			case 'C':
+				mcat_last = 0;
 				if (!n || isflt(argv[i+1]))
 					goto userr;
 				defCsym = mop[nmats].preop.csym = argv[++i];
 				mop[nmats].preop.clen = 0;
-				mcat_last = 0;
 				break;
 			case 'c':
+				mcat_last = 0;
 				if (n && !isflt(argv[i+1])) {
 					mop[nmats].preop.csym = argv[++i];
 					mop[nmats].preop.clen = 0;
@@ -791,16 +968,15 @@ main(int argc, char *argv[])
 					goto userr;
 				}
 				mop[nmats].preop.csym = NULL;
-				mcat_last = 0;
 				break;
 			case 'm':
+				mcat_last = 1;
 				if (!n) goto userr;
 				if (argv[++i][0] == '-' && !argv[i][1]) {
 					if (stdin_used++) goto stdin_error;
 					mcat_spec = stdin_name;
 				} else
 					mcat_spec = argv[i];
-				mcat_last = 1;
 				break;
 			default:
 				fprintf(stderr, "%s: unknown option '%s'\n",
@@ -868,15 +1044,21 @@ main(int argc, char *argv[])
 		fprintf(stderr, "%s: unsupported output format\n", argv[0]);
 		return(1);
 	}
-					/* process & write rows */
-	return(combine_input(&mop[nmats], stdout) ? 0 : 1);
+	doptimize(1);			/* optimize definitions */
+	i = spawned_children(nproc);	/* create multiple processes if requested */
+	if (i > 0)			/* running in parent process? */
+		return(parent_loop() ? 0 : 1);
+	if (i < 0)			/* running in output process? */
+		return(output_loop() ? 0 : 1);
+					/* else we are a worker process */
+	return(combine_input() ? 0 : 1);
 stdin_error:
 	fprintf(stderr, "%s: %s used for more than one input\n",
 			argv[0], stdin_name);
 	return(1);
 userr:
 	fprintf(stderr,
-	"Usage: %s [-h][-f{adfc}][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
+	"Usage: %s [-h][-f{adfc}][-n nproc][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
 			argv[0]);
 	return(1);
 }