code: public static void solveL_diagReal(double L[], double[] b, int n)
      {
  //        for( int i = 0; i < n; i++ ) {
  //            double sum = b[i];
  //            for( int k=0; k<i; k++ ) {
  //                sum -= L[i*n+k]* b[k];
  //            }
  //            b[i] = sum / L[i*n+i];
  //        }
          int stride = n*2;
  
          for( int i = 0; i < n; i++ ) {
              double realSum = b[i*2];
              double imagSum = b[i*2+1];
  
              int indexL = i*stride;
              int indexB = 0;
              for( int k=0; k<i; k++ ) {
                  double realL = L[indexL++];
                  double imagL = L[indexL++];
  
                  double realB = b[indexB++];
                  double imagB = b[indexB++];
  
                  realSum -= realL*realB - imagL*imagB;
                  imagSum -= realL*imagB + imagL*realB;
              }
  
              double realL = L[indexL];
  
              b[i*2    ] = realSum/realL;
              b[i*2 + 1] = imagSum/realL;
          }
      }
docs: <p>
  Solves for non-singular lower triangular matrices with real valued diagonal elements
  using forward substitution.
  <br>
  b = L<sup>-1</sup>b<br>
  <br>
  where b is a vector, L is an n by n matrix.<br>
  </p>
  
  @param L An n by n non-singular lower triangular matrix. Not modified.
  @param b A vector of length n. Modified.
  @param n The size of the matrices.