#include <stdio.h>
#include <math.h>
#define n 3
#define tol 1.e-10
/* Fortran (f): use parameter declaration */
/* Computational Linear Algebra (CLA) Problem
Using Virtual Row Pivoting and Scaling */
main()
{
int i, j, k, kp, pt, p[n]; /* CAUTION: int p[n] declares p[0] .. p[n-1] */
/* f: use integer for int */
double relres, bnorm, relnorm, Tm, Tp, Ts, Tx; /* f: use real*8 for double */
double b[n], r[n], s[n], x[n], a[n][n+1], as[n][n]; /* f: use b(n), a(n,n+1) */
/* Input;  f:  use b(1), .. ,b(n) and a(1,1), .. ,a(n,n) */
b[0] = 96.79;
b[1] = 74.85;
b[2] = 51.22;
a[0][0] = 2.178;
a[0][1] = 3.145;
a[0][2] = 9.138;
a[1][0] = 1.456;
a[1][1] = 2.805;
a[1][2] = 6.731;     
a[2][0] = 1.688;
a[2][1] = 3.763;
a[2][2] = 2.341;
printf("Forward Gaussian Elimination Using VRP and VRS:");
printf("\nOriginal [ a | b | p ] ="); /* f: use write or print */
for(i = 0; i < n; i++){ p[i] = i; a[i][n]=b[i]; /* ij save loops i*/
/* f: use do 10 i=1,n */
   for(j = 0; j < n; j++){ as[i][j] = a[i][j]; } 
      printf("\n [%11.3e%11.3e%11.3e |%11.3e | %d ]",
      a[i][0],a[i][1],a[i][2],b[i],p[i]);
   } /* end i save loops */
for(k = 0; k < n-1; k++){ /* fge loop nest */
   kp = k; /* current pivot (change) */
   for(i = k; i < n; i++){ s[p[i]] = 0; /* ij pivot and scale search loops */
      for(j = k; j < n; j++){ /* search for max scale in remainder ith row */
	 Ts = a[p[i]][j];
	 if(fabs(Ts) > s[p[i]]) {s[p[i]] = Ts;} 
	 /* f: use if(abs(Ts).gt.s(p(i))) s(p(i))=Ts */
      } /* end j scale search loop */  
      if(Ts < tol ) {
         printf("\n Scale too small warning in %dth row:%11.3e",i,Ts);
         return;}
      if(fabs(a[p[i]][k])/Ts > fabs(a[p[kp]][k]/s[p[kp]])) {kp = i;}
   } /* end  i pivot search loop */
   Tp = a[p[kp]][k];
   printf("\n %dth Pivot: a[p[%d]][%d] =%11.3e;", k, kp, k, Tp);
   if(fabs(Tp) < tol){ 
      printf("\n Pivot too small warning in %dth row:%11.3e",kp,Tp); return;}
      if(kp != k){pt = p[kp]; p[kp] = p[k]; p[k] = pt;} /* Pivot interchange */
      /* f: use if(kp.ne.k) ... */
   for(i = k+1; i < n; i++){ /* multiplier-elim. loop */
      Tm = a[p[i]][k]/Tp; a[p[i]][k] = Tm;  /* Save Multiplier Over Zeros */
      for(j = k+1; j < n+1; j++){ /* elimination Loop */
	 a[p[i]][j] = a[p[i]][j] - Tm*a[p[k]][j];} /* end j elim. loop */ 
   }   /* end i multiplier/elim. loop */
printf("\nk=%d Virtual Pivoted Form [ a | b | p | s ] =",k);
for(i = 0; i < n; i++){
  printf("\n [");
  for(j = 0; j < n; j++){ printf("%11.3e", a[i][j]); }
  printf(" |%11.3e | %d |%11.3e ]",a[i][n], p[i], s[i]); }
}   /* end k pivot loop and FGE */
printf("\nFinal Actual Pivoted Form [ Pa | Pb | p | s ] =");
for(i = 0; i < n; i++){
  printf("\n [");
  for(j = 0; j < n; j++){ printf("%11.3e", a[p[i]][j]); }
  printf(" |%11.3e | %d |%11.3e ]",a[p[i]][n], i, s[p[i]]); }
/* Begin Back Substitution Steps */
Tp = a[p[n-1]][n-1]; 
if(fabs(Tp) < tol){
   printf("\n Pivot too small warning in %dth row:%11.3e",n-1,Tp); return;}
Tx = a[p[n-1]][n]/Tp; x[n-1]=Tx; /* get x[n-1] */
printf("\nBack Substituted Results in Order of Solution:");
printf("\n x[%1d] =%11.3e;",n-1, Tx);
for(k = n-2; k >= 0; k--) { /* Begin Get x[k] Loop */
   for(i = k; i >= 0; i--) { /* Back Substitution of x[k+1] Loop */
      a[p[i]][n] = a[p[i]][n] - a[p[i]][k+1]*Tx; } /* end i loop */
   Tx = a[p[k]][n]/a[p[k]][k]; x[k] = Tx;  
   printf("\n x[%1d] =%11.3e;", k, Tx); }  /* End Get x[k] Loop */
printf("\nAnswer: FGE with BackSubstition, VRP & VRS Results: x(approx) =");  
for(k = 0; k < n; k++) { printf("\n[%11.3e ]", x[k]); } 
/* Begin Residual Norm Calculation */
bnorm = 0;
relnorm = 0;
for(i = 0; i < n; i++) {
   r[i] = b[i];
   for(k = 0; k < n; k++) { r[i] = r[i] - as[i][k]*x[k];  }
   if(fabs(b[i])>bnorm) { bnorm = fabs(b[i]); }
   if(fabs(r[i])>relnorm) { relnorm = fabs(r[i]); } } 
printf("\nEquation substitution error of residual =  b-a*x(approx) = ");
for(k = 0; k < n; k++){ printf("\n[%11.3e ]", r[k]); } 
printf("\nResidual Norm =%11.3e; b Norm =%11.3e;",relnorm,bnorm);
relres = relnorm/bnorm;
printf("\nRelative Residual = ||residual||/||b||=%11.3e", relres);
} /* end main procedure */