In this paper, two efficient iterative methods are presented to solve the symmetric and skew symmetric solutions of a linear matrix equation AXB+CYD=E, respectively, with real pair matrices X and Y. By these two iterative methods, the solvability of the symmetric and skew symmetric solutions for the matrix equation can be determined automatically. When the matrix equation has symmetric and skew symmetric solutions, then, for any initial pair matrices X0 and Y0, symmetric and skew symmetric solutions can be obtained within finite iteration steps in the absence of roundoff errors, and the minimum norm of the symmetric and skew symmetric solutions can be obtained by choosing a special kind of initial pair matrices. In addition, the unique optimal approximation pair solution X̂ and Ŷ to the given matrices X¯ and Y¯ in Frobenius norm can be obtained by finding the minimum norm solution of a new matrix equation AX˜B+CY˜D=E˜, where E˜=E−AX¯B−CY¯D. The given numerical examples demonstrate that the iterative methods are quite efficient.