code: def multiple_pdinv(A):
      
      N = A.shape[-1]
      chols = [jitchol(A[:, :, i]) for i in range(N)]
      halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols]
      invs = [dpotri(L[0], True)[0] for L in chols]
      invs = [np.triu(I) + np.triu(I, 1).T for I in invs]
      return np.dstack(invs), np.array(halflogdets)
docs: :param A: A DxDxN numpy array (each A[:,:,i] is pd)
  
      :rval invs: the inverses of A
      :rtype invs: np.ndarray
      :rval hld: 0.5* the log of the determinants of A
      :rtype hld: np.array