Basic operators and adjoints

Coding chains and arrays

With a collection of operators, we can build more elaborate operators. One method is chaining. For example, the operator product $\bold A = \bold B \bold C$ is represented in the subroutine chain2( op1, op2, ...). Likewise the operator product $\bold A = \bold B \bold C \bold D$ is represented in the in the subroutine chain3( op1, op2, op3,...).

Another amazing thing about matrices is that their elements may themselves be matrices. A very common example is a row matrix containing two elements, each element itself a matrix. This is done by subroutine row also in module chain.

api/c/chain.c

void sf_chain( sf_operator oper1     /* outer operator */, 
	       sf_operator oper2     /* inner operator */, 
	       bool adj              /* adjoint flag */, 
	       bool add              /* addition flag */, 
	       int nm                /* model size */, 
	       int nd                /* data size */, 
	       int nt                /* intermediate size */, 
	       /*@out@*/ float* mod  /* [nm] model */, 
	       /*@out@*/ float* dat  /* [nd] data */, 
	       float* tmp            /* [nt] intermediate */) 
/*< Chains two operators, computing oper1{oper2{mod}} 
  or its adjoint. The tmp array is used for temporary storage. >*/
{
    if (adj) {
	oper1 (true, false, nt, nd, tmp, dat);
	oper2 (true, add, nm, nt, mod, tmp);
    } else {
	oper2 (false, false, nm, nt, mod, tmp);
	oper1 (false, add, nt, nd, tmp, dat);
    }
}

api/c/chain.c

void sf_array( sf_operator oper1     /* top operator */, 
	       sf_operator oper2     /* bottom operator */, 
	       bool adj              /* adjoint flag */, 
	       bool add              /* addition flag */, 
	       int nm                /* model size */, 
	       int nd1               /* top data size */, 
	       int nd2               /* bottom data size */, 
	       /*@out@*/ float* mod  /* [nm] model */, 
	       /*@out@*/ float* dat1 /* [nd1] top data */, 
	       /*@out@*/ float* dat2 /* [nd2] bottom data */) 
/*< Constructs an array of two operators, 
  computing {oper1{mod},oper2{mod}} or its adjoint. >*/
{
    if (adj) {
	oper1 (true, add,  nm, nd1, mod, dat1);
	oper2 (true, true, nm, nd2, mod, dat2);
    } else {
	oper1 (false, add, nm, nd1, mod, dat1);
	oper2 (false, add, nm, nd2, mod, dat2);
    }
}

Basic operators and adjoints