Solving large linear system AQ^{-1}A'X = B

Suppose that A is a m by n full row rank sparse matrix, and Q is an n by n symmetric positive definite sparse matrix with m<n. Besides, m is about 10^5, and n is about 10^6. There is no other special structure for A and Q (i.e., not circulant, not Toeplitz, etc.). Besides, Q is always changing when it has different values for its parameters. I need to solve the linear system AQ^{-1}A'X=B for X on my Mac many times, where B is an m by m matrix.

Here is what I tried:

 1.  [R1, p1, s1] = chol(Q, 'vector');
 2.  R1A(s1, :) = R1'\A(:,s1)';  % this step most of the time fails for some Q, with my Mac frozen.
 3.  X1 = factorize(R1A'*R1A)\B; % factorize is  a function in SuiteSparse package.

Besides, I also tried to run the second line in a for loop or a parfor loop but it runs very slow, and it also make my Mac died. I tried to use UMFPACK package, and it can not solve this problem on my MAC. I preferred the direct solver. Any suggestions will be appreciated.

Thanks, Pulong