Write an optimized program in C to multiply two linear matrices.  This program should take as input the dimensions of each matrix, as well as the components of each matrix.  State any assumptions you made when writing this program.  Use the following matrices as a test:

 

1          0                           

-2         3             0      6      1

5          4             3      8      -2

0          1                           

 

Assumptions:

1.      No checking for output errors OK.

2.      Terse messages on input errors OK.

 

#include <stdio.h>

#include <stdlib.h>

 

typedef int elem_type;

#define SCAN_ELEM "%d"

#define PRINT_ELEM "%12d"

#define ELEM_ZERO 0

 

/* Matrix in the heap. */

typedef struct

  {

    unsigned num_rows, num_columns;

    /* The size of this array is really num_rows * num_columns.  Elements are

    ** stored in row-major order. */

    elem_type elem[1];

  }

matrix;

 

void *safe_alloc(size_t sz)

  {

    void *p = malloc(sz);

    if (!p)

      {

        fprintf(stderr, "\nheap exhausted\n");

        exit(1);

      }

    return(p);

  }

 

matrix *alloc_matrix_in_heap(unsigned n_rows, unsigned n_cols)

  {

    matrix *p =

      (matrix *)

        safe_alloc(sizeof(matrix) +

                   sizeof(elem_type) * (n_rows * n_cols - 1));

    p->num_rows = n_rows;

    p->num_columns = n_cols;

  }

 

matrix *matrix_mul(const matrix *a, const matrix *b)

  {

    /* The product. */

    matrix *c = alloc_matrix_in_heap(a->num_rows, b->num_columns);

 

    register const elem_type *pa = a->elem, *pb;

    elem_type *pc = c->elem;

    unsigned a_cols = a->num_columns;

    register unsigned b_cols = b->num_columns;

    const elem_type *pb_col = b->elem, *pb_col_last = pb_col + b_cols - 1;

    elem_type *pc_end = c->elem + a->num_rows * b_cols;

    register unsigned dc;

    register elem_type tmp;

 

    for ( ; ; )

      {

        dc = a_cols;

        pb = pb_col;

        tmp = ELEM_ZERO;

        do

          {

            tmp += *(pa++) * *pb;

            pb += b_cols;

          }

        while (--dc);

        *(pc++) = tmp;

        if (pc == pc_end)

          break;

        if (pb_col == pb_col_last)

          pb_col = b->elem;

        else

          {

            pb_col++;

            pa -= a_cols;

          }

      }

    return(c);

  }

 

unsigned read_unsigned(void)

  {

    unsigned u;

    if (scanf(" %u", &u) != 1)

      {

        fprintf(stderr, "\nbad input\n");

        exit(1);

      }

    return(u);

  }

 

elem_type read_elem(void)

  {

    elem_type u;

    if (scanf(" " SCAN_ELEM, &u) != 1)

      {

        fprintf(stderr, "\nbad input\n");

        exit(1);

      }

    return(u);

  }

 

matrix *read_matrix(void)

  {

    register unsigned r = read_unsigned();

    unsigned c = read_unsigned();

    matrix *m;

    register elem_type *p;

 

    if ((r == 0) || (c == 0))

      {

        fprintf(stderr, "\nbad input\n");

        exit(1);

      }

 

    m = alloc_matrix_in_heap(r, c);

 

    r *= c;

    p = m->elem;

    do

      *(p++) = read_elem();

    while (--r);

 

    return(m);

  }

 

void print_matrix(const matrix *m)

  {

    unsigned r = m->num_rows;

    register unsigned c;

    register const elem_type *p = m->elem;

    do

      {

        c = m->num_columns;

        do

          printf(PRINT_ELEM, *(p++));

        while (--c);

        putchar('\n');

      }

    while (--r);

  }

 

int main()

  {

    matrix *a = read_matrix();

    matrix *b = read_matrix();

    matrix *c;

 

    if (a->num_columns != b->num_rows)

      {

        fprintf(stderr, "\nbad input\n");

        exit(1);

      }

 

    c = matrix_mul(a, b);

 

    free(a);

    free(b);

 

    print_matrix(c);

 

    free(c);

 

    return(0);

  }