PolynomialExpansion: 2x2caseWithoutLoop.cvl

#ifdef _CIVL
#include <civlc.cvh>
#endif
#include <stdio.h>
#define n 2

$input double a1,a2,a3;
$input double b[n];
double x[n];
double xcg[n];

void cg(double A[n][n], double b[n], double x[n], int steps) {
  double r[n];
  double p[n];
  double temp[n];
  double tempp[n];
  double rsold;
  double rsnew;
  double rsfrac;
  double alpha;

  for(int i=0; i<n; i++) x[i] = 0;
  for(int i=0; i<n; i++) {
    temp[i] = 0.0;
    for(int j=0; j<n; j++) {
      temp[i] += A[i][j]*x[j];
    }
  }  // Step 1
  for(int i=0; i<n; i++) {
      r[i] = b[i] -temp[i];
      printf("\n===Residual vector r[%d] is: %f\n", i, r[i] );
  }
  rsold = 0.0;
  for(int i=0; i<n; i++) {
      rsold += r[i]*r[i];
  }
  for(int i=0; i<n; i++) {
      p[i] = r[i];
  }
  
  //first loop, i == 0;
  for(int i=0; i<n; i++) {
    temp[i] = 0.0;
    for(int j=0; j<n; j++) {
      temp[i] += A[i][j]*p[j];
    }
  }
  alpha = 0.0;
  for(int i=0; i<n; i++) {
    alpha += p[i]*temp[i];
  }
  printf("\n===alphaDenom is: %f\n", alpha);
  
#ifdef _CIVL
  $assume(alpha !=0); 
#endif
  alpha = rsold/alpha; //Step 2
  printf("\n===Scalar alpha is: %f\n", alpha);
  for(int i=0; i<n; i++) {
    tempp[i] = r[i] -alpha*temp[i];  //Step 3
    printf("\n===tempp[%d] is: %f\n", i, tempp[i]);
  }

  for(int i=0; i<n; i++) {
    r[i] = tempp[i];
  }

  for(int i=0; i<n; i++) {
    tempp[i] = x[i] +alpha*p[i]; //Step 4
  }

  for(int i=0; i<n; i++) {
    x[i] = tempp[i];
    printf("\n===x[%d] is: %f\n", i, x[i]);
  }

  rsnew = 0.0;
  for(int i=0; i<n; i++) {
    rsnew += r[i]*r[i];
  }
  printf("\n===rsnew is: %f\n", rsnew);
  
#ifdef _CIVL
  $assume(rsold !=0); 
#endif        
  rsfrac = rsnew/rsold; //Step 5
  printf("\n===Scalar beta is: %f\n", rsfrac);
  for(int i=0; i<n; i++) {
    tempp[i] = r[i] +rsfrac*p[i]; //Step 6
    printf("\n===Next search direction is: %f\n", tempp[i]);
  }

  for(int i=0; i<n; i++) {
    p[i] = tempp[i];
  }

  rsold = rsnew;

  //second loop, i == 1;
  for(int i=0; i<n; i++) {
    temp[i] = 0.0;
    for(int j=0; j<n; j++) {
      temp[i] += A[i][j]*p[j];
    }
  }
  alpha = 0.0;
  for(int i=0; i<n; i++) {
    alpha += p[i]*temp[i];
  }
  printf("\n===alphaDenom2 is: %f\n", alpha);
  
#ifdef _CIVL
  $assume(alpha !=0); 
#endif
  alpha = rsold/alpha; //Step 7
  printf("\n===Scalar alpha is: %f\n", alpha);

  for(int i=0; i<n; i++) {
    tempp[i] = r[i] -alpha*temp[i]; //Step 8
    printf("\n===tempp[%d] is: %f\n", i, tempp[i]);
  }

  for(int i=0; i<n; i++) {
    r[i] = tempp[i];
  }

  for(int i=0; i<n; i++) {
    tempp[i] = x[i] +alpha*p[i]; //Step 9
  }

  for(int i=0; i<n; i++) {
    x[i] = tempp[i];
    printf("\n===x[%d] is: %f\n", i, x[i]);
  } 
}

void main() {
  double bncg[n];
  double A[n][n];
  A[0][0] = a1;
  A[0][1] = a2;
  A[1][0] = a2;
  A[1][1] = a3;
  $assume(b[0]!=0 || b[1]!=0);
  cg(A,b,xcg,n);
  printf("\n================Solution x:================\n");
  for(int i=0; i<n; i++) {
    printf("CG(%d): %f\n",i, xcg[i]);
  }
  for(int i=0; i<n; i++) {
    bncg[i] = 0;
    for(int j=0; j<n; j++) {
      bncg[i] += A[i][j]*xcg[j];
    }
    printf("b final   (%d)=%f\n", i, b[i]);
    printf("b finalnCG(%d)=%f\n", i, bncg[i]);
    $assert(bncg[i] == b[i]);
  }
}