The Cholesky Method is our primary method of solving a GUNNS network’s [A] {p} = {w} system of equations. It is a direct method as opposed to our alternate iterative method, the SOR_Method. Although Cholesky is much slower than SOR, it is much more stable, can handle a more poorly-conditioned [A] matrix, and requires less tuning. For these reasons we have historically preferred Cholesky in GUNNS and its predecessors (PFN, EPSIM, etc), and it is a more mature and widely-used implementation than our SOR method.

We don’t use the classical Cholesky A = LL* decomposition method. Instead, we use the closely related A = LDL* (which we call Cholesky LDU) decomposition. The main reason for this is that when we first were prototyping GUNNS, we compared the LL* and LDL* methods and found them to run at very close the same speed, so the LDL* was preferred because it avoids taking the square root of the diagonal in [A], which we figured would be more robust overall. As it turns out, all GUNNS implementations to-date ensure the diagonals are positive, so LL* would be a perfectly valid option for us.

We haven’t compared the performance of the LL* and LDL* methods since GUNNS was a prototype. Two factors that might lead us to switch to LL* would be if it is significantly faster on today’s equipment, or if it exhibits less rounding errors for poorly conditioned [A].

This is a quick reference of the Cholesky LDU algorithms that GUNNS currently uses. Our implementation has extra modifications to improve speed and protect against arithmetic faults, not detailed here.

Input A is the pre-decomposed Admittance Matrix. The decomposed LDU output is stored in A. Input n is the number of rows in the matrix.

Decompose(double *A, int n):

  for i=1 to n-1:
    for j=0 to i-1:
      for k=0 to j-1:
        Aij -= Aik * Ajk

    for j=0 to i-1:
      Aii -= Aij * Aij / Ajj
      Aij /= Ajj
      Aji  = Aij

Input LDU is the decomposed Admittance Matrix from above. Input B is the network Source Vector. Output x is the Potential Vector solution. Input n is the number of rows in LDU and the size of B and x.

Solve(double *LDU, double B[], double x[], int n):

  SolveUnitLowerTriangular(LDU, B, x, n);

  for i=0 to n-1:
    x[i] /= Aii

  SolveUnitUpperTriangular(LDU, x, x, n);

SolveUnitLowerTriangular(double *L, double B[], double x[], int n):

  x[0] = B[0]
  for i=1 to n-1:
    x[i] = B[i] - (L[i][0] * x[0] + ... + L[i][i-1] * x[i-1])

SolveUnitUpperTriangular(double *U, double B[], double x[], int n):

  x[n-1] = B[n-1]
  for i=n-2 to 0:
    x[i] = B[i] - (U[i][i+1] * x[i+1] + ... + U[i][n-1] * x[n-1])