I recently started using MKL_PARDISO. I noticed that phase 331 gives the wrong result if you want to solve for multiple right hand sides using the Schur complement feature.
Attached a code to reproduce the problem. I just copied the example you provide with mkl distribution for the Schur complement and added multiple rhs.
I'm using composer_xe_2015.3.187, with MKL 11.2.3
Just a clarification about the issue.
Since the 3 rhs I'm using are all equal to one, I'm expecting that the output of phase 331 will be the same for x[0...7],x[8...15] and x[16...23]. Instead, the numbers in x[8...15] and x[16...23] (that are the same) are not in the same order as in x[0..7].