Listing 1: Matrix Add-Subtract routine

/* ------------ */
/* madsub.c     */
/* ------------ */
#include <assert.h>
#include "matdefs.h"
123456789012345678901234567890123456789012345
/* ================================================ */
/* madsub - performs matrix add/subtract operations */
/* ================================================ */
/* ------------------------------------------------- */
/* void madsub(a,nra,nca,b,nrb,ncb,c,nrc,ncc,m1,m2)  */
/* ------------------------------------------------- */
/*  Descriptions of Calling Sequence Parameters      */
/*                                                   */
/*    a   - augend or minuend matrix                 */
/*    nra - declared number of rows in a             */
/*    nca - declared number of columns in a          */
/*    b   - addend or subtrahend matrix              */
/*    nrb - declared number of rows in b             */
/*    ncb - declared number of columns in b          */
/*    c   - result matrix                            */
/*    nrc - declared number of rows in c             */
/*    ncc - declared number of columns in c          */
/*                                                   */
/* Remarks                                           */
/*                                                   */
/*    The dimensions on c determine the number of    */
/*    elements taken from matrices a and b.          */
/*                                                   */
/*    Code combinations for optional additions or    */
/*    subtractions:                                  */
/*                  elements       stored            */
/*             m1    in sum   m2   result (c)        */
/*             --    ------   --   ----------        */
/*             AB    a, b    APB     a + b           */
/*            ATB    a',b    AMB     a - b           */
/*            ABT    a, b'  MAMB   -(a + b)          */
/*           ATBT    a',b'                           */
/* ------------------------------------------------- */
# if defined(__STDC__) || defined(__PROTO__)
void
madsub( double *a, int nra, int nca,
        double *b, int nrb, int ncb,
        double *c, int nrc, int ncc, int m1, int m2)
# else
void    madsub(a, nra, nca, b, nrb, ncb, c, nrc, ncc, m1, m2)
double *a;
int     nra, nca;
double *b;
int     nrb, ncb;
double *c;
int     nrc, ncc, m1, m2;
# endif
{
    int     i, j;
    double  s;
    int     incra,              /* Increment starting a   */
            incrb,              /* Increment starting a   */
            indxa,              /* Increment element of a */
            indxb;              /* Increment element of b */

    int     nsa,                /* Next start on a   */
            nsb,                /* Next start on b   */
            nwa,                /* Next element in a */
            nwb,                /* Next element in b */
            nwc;                /* Next element in c */

    double  sb;                 /* Temp for next b */

# ifndef NDEBUG
    int     NumRowsA, NumColsA, NumRowsB, NumColsB;
# else
    nra = nra, nrb = nrb;       /* Keeps Compiler Quiet */
# endif

    if (m1 == AB || m1 == ABT)  /* Augend is a */
    {
        indxa = 1;              /* Increment on each element */
        incra = nca;            /* Increment for next start  */
# ifndef NDEBUG
        NumRowsA = nra;
        NumColsA = nca;
# endif
    }
    else                        /* Augend is a' */
    {
        indxa = nca;            /* Increment on each element */
        incra = 1;              /* Increment for next start  */
# ifndef NDEBUG
        NumRowsA = nca;
        NumColsA = nra;
# endif
    }
    if (m1 == AB || m1 == ATB)  /* Addend is b */
    {
        indxb = 1;              /* Increment on each element */
        incrb = ncb;            /* Increment for next start  */
# ifndef NDEBUG
        NumRowsB = nrb;
        NumColsB = ncb;
# endif
    }
    else                        /* Addend is b' */
    {
        indxb = ncb;            /* Increment on each element */
        incrb = 1;              /* Increment for next start  */
# ifndef NDEBUG
        NumRowsB = ncb;
        NumColsB = nrb;
# endif
    }

    assert(nrc <= NumRowsA && nrc <= NumRowsB);
    assert(ncc <= NumColsA && ncc <= NumColsB);

    /* ------------------------ */
    /* Initialize Loop Controls */
    /* ------------------------ */
    nsa = nsb = 1;
    nwa = nwb = nwc = 1;

    for (i = 1; i <= nrc; ++i)
    {
        for (j = 1; j <= ncc; ++j)
        {

            sb = b[nwb - 1];

            if (m2 == AMB)      /* If a-b */
            {
                sb = -sb;
            }
            s = a[nwa - 1] + sb;

            if (m2 == MAMB)     /* If -(a+b) */
            {
                s = -s;
            }

            c[nwc - 1] = s;     /* Store result */

            ++nwc;
            nwa += indxa;
            nwb += indxb;
        }
        nsa += incra;           /* Advance controls */
        nsb += incrb;
        nwa = nsa;
        nwb = nsb;
    }
}                               /* madsub */