On the PARDISO overview webpage under the section "Direct-Iterative Preconditioning for Nonsymmetric Linear Systems" a technique to sometimes avoid factorization of the LHS is described. The relevant parameter iparm(4) controls this technique.
The description of this technique implies it is targeted towards non-symmetric linear systems by using the conjugate-gradient squared method. However, referring to the iparm(4) parameter description, I see some conflicting information. It states at the top it uses the conjugate-gradient squared method with LU preconditioning for non-symmetric matrices and regular conjugate-gradient with Cholesky preconditioning for symmetric positive-definite matrices. This is the ideal configuration. But setting K=2 states it would use conjugate-gradient squared for symmetric positive definite matrices. This is not ideal.
What is really being used for positive-definite matrices?
For default setting, regular CG for spd or Hpd; and BiCG which applied CGS for non-symmetric A. However, if you use iparm(4) and set K=2, that means you specify to use CGS for spd. As I know, CG can only apply for symmetric, but (bi)CGS can be applied for both symmetric or non-symmetric. Because, BiCG requires two MVs to extend residual and shadow residual that one multiplication by A and one by AT. And if Pi(A) reduces r0 to smaller matrix r1, will use CGS to operate contraction process twice (compute Pi2(A)r0 ) for convergence. If your are solving problem for very large, and very sparse coefficient matrices, CGS may could provide a quick way for finding x. But if you are performing with small spd matrix, CG would be enough.