For some families of classical orthogonal polynomials defined on appropriate intervals, it is shown that the corresponding Jacobi matrices are totally positive and their bidiagonal factorizations can be accurately computed. By exploiting these facts, an algorithm to compute with high relative accuracy the eigenvalues of those Jacobi matrices, and consequently the nodes of Gaussian quadrature formulae for those families of orthogonal polynomials, is presented. An algorithm is also presented for the computation of the eigenvectors of these Jacobi matrices, and hence the weights of Gaussian quadrature formulae. Although in this case high relative accuracy is not theoretically guaranteed, the numerical experiments with our algorithm provide very accurate results.