We develop a surface hopping algorithm based on frozen Gaussian approximation for semiclassical matrix Schr(o)dinger equations.The algorithm is asymptotically derived from the Schr(o)dinger equation with rigorous approximation error analysis.The resulting algorithm can be viewed as a path integral stochastic representation of the semiclassical matrix Schr(o)dinger equations.