Editing Guide to madagascar programs (section)

==sfadd==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Add, multiply, or divide  RSF datasets.
|-
! colspan="4" | sfadd > out.rsf scale= add= sqrt= abs= log= exp= mode= [< file0.rsf] file1.rsf file2.rsf ...
|-
|  colspan="4" | The various operations, if selected, occur in the following order:<br><br>(1) Take absolute value, abs=<br>(2) Add a scalar, add=<br>(3) Take the natural logarithm, log=<br>(4) Take the square root, sqrt=<br>(5) Multiply by a scalar, scale=<br>(6) Compute the base-e exponential, exp=<br>(7) Add, multiply, or divide the data sets, mode=<br><br>sfadd operates on integer, float, or complex data, but all the input<br>and output files must be of the same data type.<br><br>An alternative to sfadd is sfmath, which is more versatile, but may be<br>less efficient.
|-
| ''bools  '' || '''abs=''' ||   || 	If true take absolute value  [nin]
|-
| ''floats '' || '''add=''' ||   || 	Scalar values to add to each dataset  [nin]
|-
| ''bools  '' || '''exp=''' ||   || 	If true compute exponential  [nin]
|-
| ''bools  '' || '''log=''' ||   || 	If true take logarithm  [nin]
|-
| ''string '' || '''mode=''' ||   || 	'a' means add (default), <br>       'p' or 'm' means multiply, <br>       'd' means divide
|-
| ''floats '' || '''scale=''' ||   || 	Scalar values to multiply each dataset with  [nin]
|-
| ''bools  '' || '''sqrt=''' ||   || 	If true take square root  [nin]
|}

<tt>sfadd</tt> is useful for combining (adding, dividing, or
multiplying) several datasets. What if you want to subtract two
datasets? Easy. Use the <tt>scale</tt> parameter as follows:
<pre>
bash&#36; sfadd data1.rsf data2.rsf scale=1,-1 > diff.rsf
</pre>
or
<pre>
bash&#36; sfadd < data1.rsf data2.rsf scale=1,-1 > diff.rsf
</pre>
The same task can be accomplished with the more general <tt>sfmath</tt> program:
<pre>
bash&#36; sfmath one=data1.rsf two=data2.rsf output='one-two' > diff.rsf
</pre>
or
<pre>
bash&#36; sfmath < data1.rsf two=data2.rsf output='input-two' > diff.rsf
</pre>
In both cases, the size and shape of <tt>data1.rsf</tt> and
<tt>data2.rsf</tt> hypercubes should be the same, and a warning message is printed out if the axis sampling parameters (such as
<tt>o1</tt> or <tt>d1</tt>) in these files are different.
====Implementation: [https://github.com/ahay/src/blob/master/system/main/add.c system/main/add.c]====
The first input file is either in the list or in the standard input.
<syntaxhighlight lang="c">
    /* find number of input files */
    if (isatty(fileno(stdin))) { 
        /* no input file in stdin */
	nin=0;
    } else {
        filename[0] = "in";
	nin=1;
    }
</syntaxhighlight>

Collect input files in the <tt>in</tt> array from all command-line
parameters that don't contain an "<tt>=</tt>" sign. The total number
of input files is <tt>nin</tt>.
<syntaxhighlight lang="c">
    for (i=1; i< argc; i++) { /* collect inputs */
	if (NULL != strchr(argv[i],'=')) continue; /* not a file */
	filename[nin] = argv[i];
	nin++;
    }
    if (0==nin) sf_error ("no input");
    /* nin = no of input files*/
</syntaxhighlight>

A helper function <tt>check_compat</tt> checks the compatibility of
input files.
<syntaxhighlight lang="c">
static void 
check_compat (sf_datatype type /* data type */, 
	      size_t      nin  /* number of files */, 
	      sf_file*    in   /* input files [nin] */, 
	      int         dim  /* file dimensionality */, 
	      const int*  n    /* dimensions [dim] */)
/* Check that the input files are compatible. 
   Issue error for type mismatch or size mismatch.
   Issue warning for grid parameters mismatch. */
{
    int ni, id;
    size_t i;
    float d, di, o, oi;
    char key[3];
    const float tol=1.e-5; /* tolerance for comparison */
    
    for (i=1; i < nin; i++) {
	if (sf_gettype(in[i]) != type) 
	    sf_error ("type mismatch: need %d",type);
	for (id=1; id <= dim; id++) {
	    (void) snprintf(key,3,"n%d",id);
	    if (!sf_histint(in[i],key,&ni) || ni != n[id-1])
		sf_error("%s mismatch: need %d",key,n[id-1]);
	    (void) snprintf(key,3,"d%d",id);
	    if (sf_histfloat(in[0],key,&d)) {
		if (!sf_histfloat(in[i],key,&di) || 
		    (fabsf(di-d) > tol*fabsf(d)))
		    sf_warning("%s mismatch: need %g",key,d);
	    } else {
		d = 1.;
	    }
	    (void) snprintf(key,3,"o%d",id);
	    if (sf_histfloat(in[0],key,&o) && 
		(!sf_histfloat(in[i],key,&oi) || 
		 (fabsf(oi-o) > tol*fabsf(d))))
		sf_warning("%s mismatch: need %g",key,o);
	}
    }
}
</syntaxhighlight>

Finally, we enter the main loop, where the input data are getting read
buffer by buffer and combined in the total product depending on the
data type.
<syntaxhighlight lang="c">
    for (nbuf /= sf_esize(in[0]); nsiz > 0; nsiz -= nbuf) {
	if (nbuf > nsiz) nbuf=nsiz;

	for (j=0; j < nin; j++) {
	    collect = (bool) (j != 0);
	    switch(type) {
		case SF_FLOAT:
		    sf_floatread((float*) bufi,
				 nbuf,
				 in[j]);	    
		    add_float(collect, 
			      nbuf,
			      (float*) buf,
			      (const float*) bufi, 
			      cmode, 
			      scale[j], 
			      add[j], 
			      abs_flag[j], 
			      log_flag[j], 
			      sqrt_flag[j], 
			      exp_flag[j]);
		    break;
</syntaxhighlight>

The data combination program for floating point numbers is
<tt>add_float</tt>.  
<syntaxhighlight lang="c">
static void add_float (bool   collect,    /* if collect */
		       size_t nbuf,       /* buffer size */
		       float* buf,        /* output [nbuf] */
		       const float* bufi, /* input  [nbuf] */  
		       char   cmode,      /* operation */
		       float  scale,      /* scale factor */
		       float  add,        /* add factor */
		       bool   abs_flag,   /* if abs */
		       bool   log_flag,   /* if log */
		       bool   sqrt_flag,  /* if sqrt */
		       bool   exp_flag    /* if exp */)
/* Add floating point numbers */
{
    size_t j;
    float f;

    for (j=0; j < nbuf; j++) {
	f = bufi[j];
	if (abs_flag)    f = fabsf(f);
	f += add;
	if (log_flag)    f = logf(f);
	if (sqrt_flag)   f = sqrtf(f);
	if (1. != scale) f *= scale;
	if (exp_flag)    f = expf(f);
	if (collect) {
	    switch (cmode) {
		case 'p': /* product */
		case 'm': /* multiply */
		    buf[j] *= f;
		    break;
		case 'd': /* delete */
		    if (f != 0.) buf[j] /= f;
		    break;
		default:  /* add */
		    buf[j] += f;
		    break;
	    }
	} else {
	    buf[j] = f;
	}
    }
}
</syntaxhighlight>