Linpack

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

//Some pseudocode:

/*
  i := 1
  j := 1
  while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
  if abs(A[k,j]) > abs(A[maxi,j]) then
  maxi := k
  end if
  end for
  if A[maxi,j] ≠ 0 then
  swap rows i and maxi, but do not change the value of i
  Now A[i,j] will contain the old value of A[maxi,j].
  divide each entry in row i by A[i,j]
  Now A[i,j] will have the value 1.
  for u := i+1 to m do
  subtract A[u,j] * row i from row u
  Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
  end for
  i := i + 1
  end if
  j := j + 1
  end while

*/


//==============================================================================
// dim - matrix dimension
// pmat - matrix with coefficients
// pfvec - vector with free members
// pmat becames trianglular after function call
//==============================================================================
int gauss_elimination(int dim, double* pmat, double* pfvec)
{
  double fMaxElem;
  double fAcc;

  int i , j, k, m;


  for(k=0; k<(dim-1); k++) // base row of matrix
    {
      // search of line with max element
      fMaxElem = fabs( pmat[k*dim + k] );
      m = k;
      for(i=k+1; i<dim; i++)
	{
	  if(fMaxElem < fabs(pmat[i*dim + k]) )
	    {
	      fMaxElem = pmat[i*dim + k];
	      m = i;
	    }
	}
    
      // permutation of base line (index k) and max element line(index m)
      if(m != k)
	{
	  for(i=k; i<dim; i++)
	    {
	      fAcc               = pmat[k*dim + i];
	      pmat[k*dim + i] = pmat[m*dim + i];
	      pmat[m*dim + i] = fAcc;
	    }
	  fAcc = pfvec[k];
	  pfvec[k] = pfvec[m];
	  pfvec[m] = fAcc;
	}

      if( pmat[k*dim + k] == 0.) return 1; // needs improvement !!!

      // triangulation of matrix with coefficients
      for(j=(k+1); j<dim; j++) // current row of matrix
	{
	  fAcc = - pmat[j*dim + k] / pmat[k*dim + k];
	  for(i=k; i<dim; i++)
	    {
	      pmat[j*dim + i] = pmat[j*dim + i] + fAcc*pmat[k*dim + i];
	    }
	  pfvec[j] = pfvec[j] + fAcc*pfvec[k]; // free member recalculation
	}
    }


  return 0;
}
//==============================================================================
// testing of function
//==============================================================================

unsigned int get_matrix_sz(int ram, double coeff)
{
  int msz;

  msz = (int) sqrt( (double)ram * coeff / sizeof(double) );

  return msz;

}

//===============================================================================
//get total number of operations
//===============================================================================
double get_ops(int N){
  return 2.0*N*N*N/3;
}

int main(int nArgs, char** pArgs)
{
  int res;
  int i,j;
  char ch;
     
  int ram;
  double coef;

  int MSZ;

  double *pmat;
  double *pfcof;

  time_t t1,t2, dt;

  double coefarr[3] = {1. / 3, 1. / 2, 2. / 3};
  int cn;

  printf("Please enter your RAM [MB]: "); 
  scanf("%d", &ram);

  for (cn = 0; cn < 3; cn++) {

    MSZ = get_matrix_sz(ram*1024*1024, coefarr[cn]);

    pmat = (double *) malloc(MSZ*MSZ*sizeof(double));
    pfcof = (double *) malloc(MSZ*sizeof(double));

    printf("Going with RAM filled on %f\n\n", coefarr[cn]);

    printf("Allocated %ux%u matrix\n", MSZ, MSZ);
     
    printf("Filling it with data ... \n");
    fflush(stdout);

    for (i = 0; i < MSZ; i++)
      for (j = 0; j < MSZ; j++)
	*(pmat + j  + MSZ * i) = (double) (rand() % 1000000);


    printf("Performing gaussian elimination...\n");
    fflush(stdout);

    t1 = time(NULL);

    res = gauss_elimination(MSZ, pmat, pfcof); // !!!

    t2 = time(NULL);

    dt = t2 - t1;

    if (!dt) {
      printf("Cannot calc approx speed. Increase problem size!\n");
      return 0;
    }

    printf("~Number of operations: %.0f\nCalculation took %d sec.\n", get_ops(MSZ), dt);
    printf("\n---------------------------------------\n");
    printf("Approx performance is %f MFLOPS\n", get_ops(MSZ) / dt / 1000000.0);
    printf("---------------------------------------\n\n");

    if(res)
      {
	printf("No solution! Try again\n");
	return 1;
      }

  }
  return 0;
}