source: CIVL/examples/cg/cg2_cholesky.cvl@ 897c8d7

/* 
 * CG 2x2 case with potitive definite assumption based on Cholesky decomposition.
 * From wiki:
 * The matrix M is positive definite if and only if there exists a unique lower 
 * triangular matrix L, with real and strictly positive diagonal elements, such 
 * that M = LL*. This factorization is called the Cholesky decomposition of M.
 */
#include <civlc.cvh>
#include <stdio.h>
#define n 2

$input double a1,a2,a3;
$assume(a1>0 && a3>0); //dignoal enties are positive
$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];
    }
  }
  for(int i=0; i<n; i++) {
    r[i] = b[i] -temp[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];
  }
  for(int i=0; i<steps; i++) {
    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];
    }
    
    $assume(alpha !=0);

    alpha = rsold/alpha;

    // tempp = r-alpha*temp
    for(int i=0; i<n; i++) {
      tempp[i] = r[i] -alpha*temp[i];
    }

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

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

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

    if(i<steps-1) {
      // rsnew = <r,r>
      rsnew = 0.0;
      for(int i=0; i<n; i++) {
        rsnew += r[i]*r[i];
      }

      $assume(rsold !=0);

      rsfrac = rsnew/rsold;
      for(int i=0; i<n; i++) {
        tempp[i] = r[i] +rsfrac*p[i];
      }

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

void main() {
  double bncg[n];
  double L[n][n]; //lower triangular matrix
  double LT[n][n];//its transpose
  double A[n][n];

  L[0][0] = a1;
  L[0][1] = 0;
  L[1][0] = a2;
  L[1][1] = a3;
  
  LT[0][0] = a1;
  LT[0][1] = a2;
  LT[1][0] = 0;
  LT[1][1] = a3;

  for(int i=0; i<n; i++) {
    for(int j=0; j<n;j++) {
      A[i][j] = 0.0;
      for(int k=0;k<n;k++) {
        A[i][j] += L[i][k] * LT[k][j]; //form the input matrix A
      }                                // to ensure A is Positive Definite
    }
  }
  
  $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("x[%d] = %f\n\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];
    }
    $assert(bncg[i] == b[i]);
    }
}