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

int main() {
   double *a, *b, *c, *d, *x;
   int   i, j, l, n, p, q;

   n = 12;  /* dimension of matrix A */
   p = 3;  /* number of blocks in one direction */

   q = n/p; 

   a = (double *) malloc( n * sizeof(double)  );
   b = (double *) malloc( n * sizeof(double)  );
   c = (double *) malloc( n * sizeof(double)  );
   d = (double *) malloc( n * sizeof(double)  );
   x = (double *) malloc( n * sizeof(double)  );

   srand( (unsigned)time( (long *)0 ));
   for(i = 0; i < n ; i++) {
/*
      a[i] = rand() /2147483647.0;  
      b[i] = rand() /2147483647.0;      
      c[i] = rand() /2147483647.0;      
      d[i] = rand() /2147483647.0;   
*/
      a[i] = 1.0; 
      b[i] = 4.0;
      c[i] = 5.0;
      d[i] = 3.0;
   }
   a[0]   = 0.0;
   c[n-1] = 0.0; 

   /* step 1: */
   for(l = 0; l < p ; l++) {
      for(j = 1; j < q ; j++) {    
         b[l*q+j] += -a[l*q+j] / b[l*q+j-1] * c[l*q+j-1];
         d[l*q+j] += -a[l*q+j] / b[l*q+j-1] * d[l*q+j-1];
         if(l > 0) {
            a[l*q+j] = -a[l*q+j] / b[l*q+j-1] * a[l*q+j-1];
         }
      }
   }

   /* step 2: */
   for(l = 1; l < p ; l++)
      for(j = q-3; j >= 0 ; j--)
         a[l*q+j] +=  -c[l*q+j] / b[l*q+j+1] * a[l*q+j+1];

   for(l = 0; l < p ; l++)
      for(j = q-3; j >= 0 ; j--) {
         d[l*q+j] +=  -c[l*q+j] / b[l*q+j+1] * d[l*q+j+1];
         c[l*q+j]  =  -c[l*q+j] / b[l*q+j+1] * c[l*q+j+1];
   }
   for(l = 1; l < p ; l++) {
      b[l*q-1] += -c[l*q-1] / b[l*q] * a[l*q];
      d[l*q-1] += -c[l*q-1] / b[l*q] * d[l*q];
      c[l*q-1]  = -c[l*q-1] / b[l*q] * c[l*q];
   }

   /* step 3: */
   for(l = 1; l < p ; l++) {
      for(j = 0; j <= q-2 ; j++) {
         c[l*q+j] += -a[l*q+j] / b[l*q-1] * c[l*q-1];
         d[l*q+j] += -a[l*q+j] / b[l*q-1] * d[l*q-1];
      }
      b[(l+1)*q-1] += -a[(l+1)*q-1] / b[l*q-1] * c[l*q-1];
      d[(l+1)*q-1] += -a[(l+1)*q-1] / b[l*q-1] * d[l*q-1];
   }

   /* step 4: */
   for(l = p-1; l > 0  ; l--)
      for(j = -1; j <= q-2 ; j++)
         d[l*q+j] += -c[l*q+j] / b[(l+1)*q-1] * d[(l+1)*q-1];
   l = 0; 
   for(j = 0; j <= q-2 ; j++)
      d[l*q+j] += -c[l*q+j] / b[(l+1)*q-1] * d[(l+1)*q-1];

   /* calculate the solution: */ 
   for(l = 0; l < p ; l++)
      for(j = 0; j < q ; j++)
         x[l*q+j] = d[l*q+j] / b[l*q+j];

   /* print the solution: */
   for(j = 0; j < n; j++)
      printf(" x[%d] = %6.2f\n", j, x[j] );
}