How can I multiply an NxN matrix A in Fortran x times to get its power without amplifying rounding errors?
If A can be diagonalized as
A P = P D,
where P is some NxN matrix (each column is called 'eigenvector'), and D is an NxN diagonal matrix (the diagonal elements are called 'eigenvalues'), then
A = P D P^{-1},
where P^{-1} is the inverse matrix of P. Therefore the second power of A results in
A A= P D P^{-1} P D P^{-1} = P D D P^{-1}.
Repeating multiplication of A for x times yields
A^x = P D^x P^{-1}.
Note here that D^x is still a diagonal matrix. Let the i-th diagonal element of D be D_{ii}. Then, the i-th diagonal element of D^x is
[D^x]_{ii} = (D_{ii})^x.
That is, the elements of D^x is simply x-th power of the elements of D and can be computed without much rounding error, I guess. Now, you multiply P and P^{-1} from left and right, respectively, to this D^x to obtain A^x. The error in A^x depends on the error of P and P^{-1}, which can be calculated by some subroutines in numerical packages such as LAPACK.