ahay/src/authoring-sf-programs-c

When to use

Load this skill whenever you are writing a new sf<name> main program in C. C is the reference implementation for Madagascar: the Fortran 77, Fortran 90, C++, Python, Julia, Java, Matlab, Octave, and Chapel APIs all wrap or mirror the C API. Every claim in this skill is grounded in api/c/rsf.h (the amalgamated public header at build/api/c/rsf.h) and the test programs under api/c/. (CUDA programs exist in user directories but Madagascar has no dedicated api/cuda/ binding; they use the C API for RSF I/O and CUDA for compute.)

This skill is C-specific. For language-agnostic conventions (file naming, self-doc format, parameter style, build integration) see the companion:

• ../authoring-sf-programs/SKILL.md — shared conventions (load this too)

Typical placement: user/<youruser>/M<name>.c. The build system finds every M*.c in a user directory and compiles it into sf<name>.

Skeleton

This is the minimal correct skeleton for a C main program. Copy it verbatim and extend from here.

#include <rsf.h>
 
int main(int argc, char* argv[])
{
    sf_file in, out;
    int n1, n2;
    float d1, o1;
    float *trace;
 
    sf_init(argc, argv);
    in  = sf_input("in");
    out = sf_output("out");
 
    if (!sf_histint  (in, "n1", &n1)) sf_error("No n1= in input");
    if (!sf_histint  (in, "n2", &n2)) n2 = 1;
    if (!sf_histfloat(in, "d1", &d1)) d1 = 1.0f;
    if (!sf_histfloat(in, "o1", &o1)) o1 = 0.0f;
 
    trace = sf_floatalloc(n1);
 
    for (int i2 = 0; i2 < n2; i2++) {
        sf_floatread(trace, n1, in);
        /* process trace here */
        sf_floatwrite(trace, n1, out);
    }
 
    free(trace);
    exit(0);
}

Key invariants:

• sf_init must be the first RSF call.
• Open all sf_input files before sf_output — the first output inherits dimensions from the first input by default.
• Always call exit(0) (not return 0) at the end; this is the Madagascar convention.

Self-doc header

Every M<name>.c must begin with a comment block that the build system extracts to produce sfdoc sf<name> output. The format (sourced from user/fomels/Mpick.c in this tree):

/* Automatic picking from semblance-like panels.
 
Takes: rect1=1 rect2=1 ...
 
rectN defines the size of the smoothing stencil in N-th dimension.
 
Theory in Appendix B of:
S. Fomel, 2009,
Velocity analysis using AB semblance: Geophysical Prospecting, v. 57, 311-321.
Reproducible version in RSFSRC/book/tccs/avo
http://ahay.org/RSF/book/tccs/avo/paper_html/
 
August 2012 program of the month:
http://ahay.org/blog/2012/08/01/program-of-the-month-sfpick/
*/
 
/*
  Copyright (C) 2004 University of Texas at Austin
 
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

Rules extracted from real programs in this tree:

1. First comment — brief one-line description on the opening line, optional extended prose below. This entire block becomes the program's documentation.
2. Parameter inline comments — each sf_getint/sf_getfloat/sf_getbool call is followed immediately by a /* description */ comment on the next line. The doc extractor captures these and lists them as named parameters in sfdoc.
3. Second comment — the GPL copyright block. Keep it verbatim.
4. #include <rsf.h> comes after both comment blocks, before any other code.

Example parameter documentation from user/fomels/Mpick.c:

if (!sf_getint("niter",&niter)) niter=100;
/* number of iterations */
if (!sf_getfloat("an",&an)) an=1.;
/* axes anisotropy */
if (!sf_getint("gate",&gate)) gate=3;
/* picking gate */
if (!sf_getbool("smooth",&smooth)) smooth=true;
/* if apply smoothing */

Each /* ... */ directly following a sf_get* call is parsed as the parameter description; it appears under "parameters:" in sfdoc sf<name>.

I/O API

Signatures from build/api/c/rsf.h.

Initialization

void sf_init(int argc, char *argv[]);

Must be the first RSF call. Parses the command line into an internal parameter table; opens stdin/stdout for piped RSF data.

Opening files

sf_file sf_input  (const char* tag);
sf_file sf_output (const char* tag);

tag is a logical name used on the command line: sfmyprog in=data.rsf out=result.rsf. Use "in" and "out" for the primary pipes (stdin / stdout). Additional named files ("den", "mask", etc.) are opened the same way. Output files inherit dimension headers from the first input by default — override with explicit sf_put* calls before the first write.

sf_datatype sf_gettype(sf_file file);

Returns the data type (SF_FLOAT, SF_COMPLEX, SF_INT, etc.). Use to validate input type before reading.

void sf_fileclose(sf_file file);

Closes an sf_file handle and flushes metadata. For secondary outputs (not stdin/stdout) it is good practice to call this explicitly; for primary stdin/stdout it is called implicitly on exit.

Reading and writing float data

void sf_floatread  (float* arr, size_t size, sf_file file);
void sf_floatwrite (float* arr, size_t size, sf_file file);

Read/write size consecutive floats. arr must be pre-allocated. These are the most common I/O calls.

Reading and writing complex data

void sf_complexread  (sf_complex* arr, size_t size, sf_file file);
void sf_complexwrite (sf_complex* arr, size_t size, sf_file file);

sf_complex is float _Complex (or the platform equivalent). Use when the input type is SF_COMPLEX.

Header metadata API

RSF files carry key=value pairs in a text header. There are two directions:

• Read from input header — sf_hist* family. Returns true if the key exists.
• Write to output header — sf_put* family. Must be called before the first data write.

Signatures from build/api/c/rsf.h:

Reading header values

bool sf_histint    (sf_file file, const char* key, int* par);
bool sf_histfloat  (sf_file file, const char* key, float* par);
char* sf_histstring(sf_file file, const char* key);

sf_histstring returns a heap-allocated string (or NULL); test with if (NULL != (...)). Example:

char *label;
if (NULL != (label = sf_histstring(scn, "label2")))
    sf_putstring(pik, "label", label);

Writing header values

void sf_putint    (sf_file file, const char* key, int par);
void sf_putfloat  (sf_file file, const char* key, float par);
void sf_putstring (sf_file file, const char* key, const char* par);

Axis objects

Madagascar uses sf_axis structs to bundle n, d, o, label, and unit for each dimension. These simplify copying or transforming axes between files.

sf_axis sf_iaxa(sf_file FF, int i);    /* read axis i (1-based) */
void    sf_oaxa(sf_file FF, const sf_axis AA, int i); /* write axis i */

Example — copy axis 1 from input to output unchanged, then set a new axis 2:

sf_axis ax1, ax2;
ax1 = sf_iaxa(in, 1);
sf_oaxa(out, ax1, 1);
sf_putint(out, "n2", new_n2);
sf_putfloat(out, "d2", new_d2);

sf_filedims is a convenience wrapper that reads all n values at once:

int sf_filedims(sf_file file, int *n);  /* returns number of dimensions */

Parameter API

Command-line parameters are set by the caller as key=value pairs. Signatures from build/api/c/rsf.h:

Scalar getters

bool sf_getint      (const char* key, int* par);
bool sf_getfloat    (const char* key, float* par);
bool sf_getbool     (const char* key, bool* par);
char* sf_getstring  (const char* key);           /* returns NULL if absent */
bool sf_getlargeint (const char* key, off_t* par);

All scalar getters return true if the key was found (false otherwise). sf_getstring returns a pointer or NULL.

Array getters

bool sf_getints    (const char* key, int* par,   size_t n);
bool sf_getfloats  (const char* key, float* par, size_t n);
bool sf_getbools   (const char* key, bool* par,  size_t n);

Array getters read up to n comma-separated values (e.g. rect=3,5,1).

Default-value pattern

The canonical Madagascar idiom — always provide a fallback default when the parameter is optional:

if (!sf_getint("n", &n))     n = 100;
/* number of samples */
if (!sf_getfloat("eps", &eps)) eps = 0.01f;
/* regularization parameter */
if (!sf_getbool("verb", &verb)) verb = true;
/* verbosity flag */

The /* ... */ on the line after each sf_get* call is parsed by the doc extractor and shown in sfdoc sf<name>.

For required parameters that have no sane default, call sf_error when absent:

if (!sf_histint(in, "n1", &n1)) sf_error("No n1= in input");

Memory API

Signatures from build/api/c/rsf.h. All allocators abort with a meaningful error message on allocation failure — never return NULL.

1-D allocators

float*      sf_floatalloc  (size_t n);
sf_complex* sf_complexalloc(size_t n);
int*        sf_intalloc    (size_t n);

2-D and 3-D allocators

float**      sf_floatalloc2 (size_t n1, size_t n2);
float***     sf_floatalloc3 (size_t n1, size_t n2, size_t n3);
sf_complex** sf_complexalloc2(size_t n1, size_t n2);
int**        sf_intalloc2   (size_t n1, size_t n2);
int***       sf_intalloc3   (size_t n1, size_t n2, size_t n3);

The multi-dimensional allocators return a pointer-to-pointer (standard C array-of-arrays), but the underlying data is one contiguous block — out[0] (for 2-D) or out[0][0] (for 3-D) points to it. This makes it safe to pass arr[0] to sf_floatread/sf_floatwrite as a flat buffer:

float **scan = sf_floatalloc2(n1, n2);
sf_floatread(scan[0], n1 * n2, in);   /* reads entire 2-D block at once */

Deallocation

Madagascar does not provide a custom deallocator. Use standard free() directly:

free(trace);          /* 1-D */
free(scan[0]);        /* 2-D: free the contiguous data block */
free(scan);           /* then free the pointer array */

Most production Madagascar programs skip explicit free calls before exit(0) because the OS reclaims all memory on process exit. This is acceptable; add free calls when running under Valgrind or in regression test programs that check for leaks.

Error API

Signatures from build/api/c/rsf.h:

void sf_error  (const char *format, ...);
void sf_warning(const char *format, ...);

sf_error prints the formatted message to stderr, then calls exit(1). Use it for unrecoverable conditions (bad input, missing required header key, wrong data type):

if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
if (!sf_histint(in, "n1", &n1)) sf_error("No n1= in input");
sf_error("Size mismatch [n%d]: %d != %d", j+1, m[j], n[j]);

sf_warning prints to stderr but does not exit — execution continues. Use it for progress messages and non-fatal conditions:

sf_warning("cmp %d of %d;", i3+1, n3);  /* semicolon suppresses newline */
sf_warning(".");                          /* trailing dot flushes the line */

Both functions accept printf-style format strings and variadic arguments.

Extra libraries

Most programs need only #include <rsf.h> and link against librsf. When a program needs helpers from api/c/ (FFT, solvers, splines, etc.) or a user-directory helper module, the SConstruct in the user directory must declare the dependency.