#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 partitions (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)  );

  for(i=0; i<n ;i++) {
     a[i] = 1.0; 
     b[i] = 4.0;
     c[i] = 5.0;
     d[i] = 3.0;
  }

  a[0]   = 0.;
  c[n-1] = 0.;

  /* step 1: Elimination of the subdiagonal elements */ 
  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: Elimination of the upper side diagonal */
  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: transformation to an upper triangular matrix: eliminate
             a's by means of the last line of the upper neighbour */
  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: transform matrix to diagonal form*/
  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];

  /* 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] );
}