CS13002 PDS Lab, Spring 2003, Section 1/A, Assignment 6, Solution

CS13002 PDS Lab, Spring 2003, Section 1/A

Solution of Assignment 6

Read the in-line documentations in the program below to understand how the program works.

#include <stdio.h>
#include <stdlib.h>

#define MAX_DIMEN 20
typedef struct {
   int rowdim;                              /* Number of rows */
   int coldim;                              /* Number of columns */
   float element[MAX_DIMEN][MAX_DIMEN];     /* The elements */
} matrix;

typedef struct {
   int dim;                                 /* Dimension of a vector */
   float element[MAX_DIMEN];                /* The elements */
} vector;

int inverse ( matrix *B , matrix A )
/* Computes the inverse of a matrix.                                     */
/* A is the matrix to invert, B is the matrix to store the inverse of A. */
/* Return 1 on successful computation of inverse, 0 otherwise.           */
{
   int i,j,k,n;
   float t;

   n = A.rowdim;
   if (n != A.coldim) {
      fprintf(stderr, "Error in inverse: square matrix expected\n");
      return(0);
   }

   /* Initialize B to the identity matrix */
   (*B).rowdim = (*B).coldim = n;
   for (i=0;i<n;i++) {
      (*B).element[i][i] = 1;
      for (j=i+1;j<n;j++) (*B).element[i][j] = (*B).element[j][i] = 0;
   }

   for (i=0;i<n;i++) {
      /* Find non-zero entry in the i-th column */
      k = i;
      while ((k<n)&&(A.element[k][i]==0)) k++;
      if (k == n) return(0);
      if (k != i) {   /* Interchange i-th and k-th rows */
         for (j=0;j<i;j++) {
            /* There is no fun in interchanging 0's in A */

            t = (*B).element[i][j];
            (*B).element[i][j] = (*B).element[k][j];
            (*B).element[k][j] = t;
         }
         for (j=i;j<n;j++) {
            t = A.element[i][j];
            A.element[i][j] = A.element[k][j];
            A.element[k][j] = t;

            t = (*B).element[i][j];
            (*B).element[i][j] = (*B).element[k][j];
            (*B).element[k][j] = t;
         }
      }

      /* Multiply the i-th row by 1/a_ii, thereby forcing a_ii := 1 */
      t = 1.0 / A.element[i][i];
      A.element[i][i] = 1;
      for (j=0; j<=i; j++) (*B).element[i][j] *= t;
      for (; j<n; j++) {
         A.element[i][j] *= t;
         (*B).element[i][j] *= t;
      }

      /* Reduce other rows */
      for (k=0; k<n; k++) if (k != i) {
         t = A.element[k][i];
         for (j=0; j<=i; j++) (*B).element[k][j] -= t * (*B).element[i][j];
         A.element[k][i] = 0;
         for (j=i+1; j<n; j++) {
            A.element[k][j] -= t * A.element[i][j];
            (*B).element[k][j] -= t * (*B).element[i][j];
         }
      }
   }

   return(1);
}

void printMatrix ( matrix A )
/* Formatted printing of matrices */
{
   int i,j;

   for (i=0;i<A.rowdim;i++) {
      for (j=0;j<A.coldim;j++) printf("%9.5f ",A.element[i][j]);
      printf("\n");
   }
}

void printVector ( vector v )
/* Formatted printing of vectors. Well, our vectors are column vectors,
   but this printing makes vectors look like row vectors. Treat them as
   column vectors, whatever they look like... */
{
   int i;

   for (i=0;i<v.dim;i++) printf("%9.5f ",v.element[i]);
   printf("\n");
}

void randMatrix ( matrix *A , int m , int n )
/* Generate a random mxn matrix with integer entries between -10 and +10 */
{
   int i,j;

   (*A).rowdim = m;
   (*A).coldim = n;
   for (i=0;i<m;i++) for (j=0;j<n;j++) (*A).element[i][j] = (rand()%21)-10;
}

void randVector ( vector *v , int m )
/* Generate a random m-dimensional vector with integer entries
   between -10 and +10 */
{
   int i;

   (*v).dim = m;
   for (i=0;i<m;i++) (*v).element[i] = (rand()%21)-10;
}

void matVecMul ( vector *v , matrix A , vector u )
/* Multiplication of an mxn matrix A with an n-dimensional column vector u.
   The result (stored in *v) is an m-dimensional column vector. */
{
   int i,j,m,n;

   m = A.rowdim;
   n = A.coldim;
   if (n != u.dim) {
      fprintf(stderr, "Error in matVecMul: Incompatible dimensions...\n");
      return;
   }
   (*v).dim = m;
   for (i=0;i<m;i++) {
      (*v).element[i] = 0;
      for (j=0;j<n;j++) (*v).element[i] += A.element[i][j] * u.element[j];
   }
}

int main ()
{
   matrix A,Ainv;
   vector b,x;
   int n = 8;

   /* Seed the random number generator for getting different random systems */
   srand((unsigned int)time(NULL));

   /* Generate random system */
   randMatrix(&A,n,n);
   printf("A :\n"); printMatrix(A);
   randVector(&b,n);
   printf("b :\n"); printVector(b);

   if (inverse(&Ainv,A)) {
      printf("Inverse of A :\n"); printMatrix(Ainv);
      matVecMul(&x,Ainv,b);
      printf("x :\n"); printVector(x);
      printf("Original b :\n"); printVector(b);
      matVecMul(&b,A,x);
      printf("Reconstructed b :\n"); printVector(b);
   } else {
      printf("Error: unable to invert matrix... giving up!\n");
   }
}

Let me point out a couple of topics in connection with the above code.

The portion of A already processed is not altered by the loop of the routine inverse. We know that during the i-th step we need not tamper the first i-1 columns of A. The elementary operations reducing the i-th column do not affect the first i-1 columns. For example, interchanging zeros produces no effect. Of course, you can go for interchanging zeros and/or adding zeros, thereby sacrificing efficiency. That won't lead to errors, after all, but to some unnecessary wastage of computing resources. While solving big systems this wastage could be considerable. Note, however, that B grows from the identity matrix to a (possibly) dense matrix. We cannot skip certain columns of B during the elementary row operations.
Check how the functions `returning' matrices or vectors are declared. A typical example is:
```
int inverse ( matrix *B , matrix A );
```
This is called as:
```
inverse(&Ainv,A);
```
You may wonder what business the star in the prototype and the ampersand in the call are doing. These are necessary because of C's call-by-value parameter passing mechanism. We want that after the routine inverse returns, Ainv would store the inverse of A. If we had declarations and calls like the following:
```
int inverse ( matrix B , matrix A );
...
inverse(Ainv,A);
```
then the desired behavior is not exhibited. In this case the elements (including the entire array) in Ainv is copied to the corresponding elements in the parameter B. Any change in B inside the body of the function is local and does not affect Ainv in any way. So though the function computes and stores in B the inverse of A, the results of this computation is not visible outside (namely in main).
Whenever we want the value of a variable (here, Ainv) to be changed by a function, we must pass to the function a pointer to it (&Ainv). Then in the original inverse routine B receives a copy of this pointer. This copying does not matter, because copies of a pointer point to the same location. And it is that location that we want to modify.
Instead of passing the output matrix to the function, the function may actually return the matrix. For example, the following also works.
```
matrix inverse ( matrix A )
{
   matrix B;

   /* Compute in B the inverse of A */
   ...

   return(B);
}
...
Ainv = inverse(A);
```
In this case the local variable B is returned and assigned (i.e., copied) to Ainv and then finally destroyed. Thus Ainv now receives the changes we desire.
So far so good, but returning structures is an unsafe practice, if the structure contains pointers to dynamic arrays (instead of static arrays, as the case was in this exercise). To see why let us redefine the matrix data type as:
```
typedef struct {
   int rowdim;
   int coldim;
   float **elements;
} matrix2;
```
Now consider functions and calls involving return of matrix2 data types:
```
matrix2 inverse ( matrix2 A )
{
   matrix2 B;

   /* Allocate suitable amount of memory to B.element */
   ...
   /* Compute in B the inverse of A */
   ...
   return(B);
}
...
Ainv = inverse(A);
```
This call is dangerous and is almost certain to lead to segmentation faults. The reason is the following: B.element is allocated memory in the function. The return statement makes a copy of B to Ainv. Thus Ainv receives a copy of the element pointer, but not that of the locations pointed to by element. After the routine is returned, B is destroyed and the memory allocated to its element is returned to the system pool. This memory can later be reallocated to other data or even to other processes. But Ainv still knows that it has to read from these memory locations. This leads to erroneous data and/or run-time errors (segmentation faults).
While dealing with structures involving pointer(s), always pass pointer to the structure that you want to modify. That is, use the following:
```
int inverse ( matrix2 *B , matrix2 A );
...
inverse(&Ainv,A);
```
This is 100% safe. Also passing of a structure is more efficient with dynamic arrays than with static arrays. A static array leads to element-wise copying, whereas dynamic arrays lead to copying only the pointers. Faster as it is, there is another point to notice. A is passed as a read-only parameter and we don't want its fields to be modified inside the function. Under the call-by-value mechanism of C, any change in A.rowdim, A.coldim or A.element inside the function does not affect the A outside. However, modifying the array location A.element[i][j] inside the function produces the same changes in the array locations outside the function.
Be very very careful while handling pointers and take time and effort to understand what is passed, what is copied, what is changed and what changes are visible outside! Pointers give you immense control over memory, but this facility does not come for free. You have to be extremely watchful, so that you do not end up making unauthorized memory accesses. Understand and then use. Enjoy the power of C.

[Course home] [Home]

CS13002 PDS Lab, Spring 2003, Section 1/ASolution of Assignment 6

CS13002 PDS Lab, Spring 2003, Section 1/A

Solution of Assignment 6