double dot (double v1[], double v2[], int N)
    {
            double dot_product=0.0;
            int i;
            for (i=0; i<N, ++i)
            {
                dot_product+=v1[i]*v2[i];
            }
    return dot_product;
    }

Here's how we would call it:
int main(void){
        double a[3], b[3], a_dot_b;
        int N;
        N=Ndim;
                ….
        a_dot_b = dot(a,b,N);
                …
        return 0;}
 

Last time we also discussed that there are two (straightforward) ways to get output out of a C function back to the calling routine:

1) Return a scalar, e.g. double bounded_newton(…) returns a double precision number to the calling routine:
double Xbest;
Xbest = bounded_newton(….);

2) Give a memory address in the argument list; then the function can write to that memory address (and to adjacent memory blocks).
int N_evals;
Xbest = bounded_newton(…., &N_evals,…);

for this to work we have to define bounded_newton() so that it expects
a memory address as an argument:

bounded_newton(…,int *pN_evals,…)
{
        …
        pN_evals += 1; /* each time a function is called */
}

Arrays are specified by the memory address of the first element, so you can use this to output a whole array of numbers:

int scale(double Vin[], double Vout[], double a, int N)
{
        int j;
        for(j=0;j<N;++j)
        {
                Vout[j] = a*Vin[j];
        }
return N;
}

In the calling routine both arrays have to be declared before you pass them down:
double X[10], Y[10], multiplier=7.0;
int N=5;
X[0] =….
X[1] = …
scale(X,Y,multiplier,N);

This works because the array names X, Y are actually memory addresses

int I=3;
double x=4.0;
double a[3] = {4,5,6};

if you ask the computer what I is, it will say 3.
it will say x is 4.0
it will say a is 0x43219 A MEMORY ADDRESS!

To get the numbers stored in a, you need indices: a[0], a[1], …

Of course the computer knows the memory location of I and x as well, but
they are called &I and &x. If we pass those memory addresses down, a function can stuff numbers in them too.

So far we have discussed:
Scalar Functions of Scalar Arguments
Scalar Functions of Vector arguments (e.g. dot_product() )
Scalar Functionals (e.g. root-finding, integration, numerical derivative)
Vector Functions (e.g. scale(), cross_product)

Now let's do
VECTOR FUNCTIONALS: e.g. Numerical Gradient

int grad(double f(double v[]), double x[], double gradf[], double ATOL, double RTOL, int Ndim); /*returns gradf[], the gradient of f at x. */

first define a function which computes the derivative along one dimension, with index direction (i.e. if direction=0, we want d/dx; if direction=1 we want d/dy):

double derivN(double f(double v[]), double x[], int Ndim, int direction, double RTOL, double ATOL)
{
        …
        h = …. /* depends on tolerances etc */
        for(j=0;j<Ndim;++j)
        {
                x_plus_h[j] = x[j];
                if(j==direction) x_plus_h[j] += h;
        }
        return ((f(x_plus_h)-f(x))/h);
}

int grad(double f(double v[]), double x[], double gradf[], double ATOL, double RTOL, int Ndim)
{
        for(direction=0;direction<Ndim;++I)
        {
                gradf[direction]=derivN(f,x,Ndim,direction,RTOL,ATOL);
          }
return 0;}

You might actually call grad() from main() like this:

double X0[3],grad_phi[3],ATOL,RTOL;
…
grad(phi,X0,grad_phi,ATOL,RTOL,3);
…

where phi() was declared and defined as a scalar function of a 3-vector.

Note how phi replaces f, X0 replaces x, grad_phi replaces gradf

Multi-d arrays

M[i][j], T[l][m][n], …

for 2-d arrays, we use I to index the row, and j to index the column.
Inside the computer, the whole array is actually stored sequentially;
the convention is that the last index increases most rapidly, so

M[2][2] is stored this way: M00, M01, M10, M11.

say we wanted to store a whole bunch of ordered pairs (x,y), say 20 pairs,
and then we want to send these pairs down one at a time to a function which
computed the distance between (x,y) and the origin (0,0).

Here is one way:

double distance_from_origin(double v[]);

double xydataset[20][2];
double *xy[20];
double vector[2];

xy[j] = &xydataset[j][0]; /* xy[j] is the memory location where the jth
row of xydataset starts */

vector = xy[j];

r = distance_from_origin(vector);
 

Note that this works if we just use a 1-d array for the data, e.g.

double data[40];
…
xy[j] = &data[2*j];

xy has all the information needed to find whatever data values you
want in data[] or xydataset[][], so you can write programs that just
take xy as their argument, without any need to specify how the original
data set was dimensioned: this is exactly what the Numerical Recipes
routines do.

The amazing thing is that in C

xy[m][n] == xydataset[m][n]

since xy[m][n] = *(xy[m]+n) = *(*(xy+m)+n)

Note that xy is the memory address of a memory address of a double, so
you can see why matrices are often handled using pointers-to-pointers:

double **M;

this is discussed more in the ABC text, and in the beginning of the NRC text. We will talk more about it next lecture.

Using the same memory address for both input and output:

Note that you can use the same memory addresses for both input and output
(but be careful!! you are overwriting the input!):

        int miles_to_km(double distances[], int N)
{
        int I;
        for(I=0;I<N;I+=1)
        {
                distances[I]=distances[I]*1.609344;
        }
        return N;
}

usually it is better to use different arguments for input and output, to keep your functions modular, encapsulated, and general. The user could actually pass the same memory addresses down as both the input and the output argument, e.g.:

int scale(double Vin[], double Vout[], double a, int N);

can be invoked using
double X[3]={3,4,5};
scale(X,X,7.0,3);

and X[] will be changed to {21,28,35}