I am sampling covariance matrix from a Inverse Wishart distribution. In one dimensional case, after doing sufficient iterations I am taking the mode value for variance (after removing the burn-in values). How to do the same in a multivariate case?
You shouldn't need a Gibbs sampler: the mode of an inverse-Wishart has a closed form.
Also, independent random samples from the Cholesky factor of a Wishart can be obtained from the Bartlett decomposition: as it is triangular, it can be inverted easily by forward subsitution to get the Cholesky factor of an inverse-Wishart.
If I understand your question correctly:
Covariance matrix for 1-dim case reduces to the variance. Wishart Distribution (or Inv wishart distribution depending on your formulation) is a prior of covariance matrices, which for dimensions $\geq$ 2 correspond to multivariate case.
However, I may have misunderstood you. Please correct me if that is the case.