docs: def cholesky(A, ordering_method='default', return_type=RETURN_P_L, use_long=False):
      '''
      P A P' = L L'
      '''
      logger.debug('Calculating cholesky decomposition for matrix {!r} with ordering method {}, return type {} and use_long {}.'.format(A, ordering_method, return_type, use_long))
  
      ## check input
      return_types = (RETURN_L, RETURN_L_D, RETURN_P_L, RETURN_P_L_D)
      if ordering_method not in CHOLMOD_ORDERING_METHODS:
          raise ValueError('Unknown ordering method {}. Only values in {} are supported.'.format(ordering_method, CHOLMOD_ORDERING_METHODS))
      if return_type not in return_types:
          raise ValueError('Unknown return type {}. Only values in {} are supported.'.format(return_type, return_types))
          if ordering_method != 'natural' and return_type in (RETURN_L, RETURN_L_D):
              raise ValueError('Return type {} is only supported for "natural" ordering method.'.format(return_type))
  
      #TODO symmetry check
      A = util.math.sparse.check.sorted_squared_csc(A)
  
      ## calculate cholesky decomposition
      try:
          try:
              f = sksparse.cholmod.cholesky(A, ordering_method=ordering_method, use_long=use_long)
          except sksparse.cholmod.CholmodTooLargeError as e:
              if not use_long:
                  warnings.warn('Problem to large for int, switching to long.')
                  return cholesky(A, ordering_method=ordering_method, return_type=return_type, use_long=True)
              else:
                  raise
      except sksparse.cholmod.CholmodNotPositiveDefiniteError as e:
          raise util.math.matrix.NoPositiveDefiniteMatrixError(A, 'Row/column {} makes matrix not positive definite.'.format(e.column))
      del A
  
      ## calculate permutation matrix
      p = f.P()
      n = len(p)
      if return_type in (RETURN_P_L, RETURN_P_L_D):
          P = scipy.sparse.dok_matrix((n,n), dtype=np.int8)
          for i in range(n):
              P[i,p[i]] = 1
          P = P.tocsr()
          P.astype(np.int8)
  
      ## return P, L
      if return_type in (RETURN_L, RETURN_P_L):
          L = f.L().tocsr()
          if return_type == RETURN_L:
              assert np.all(p == np.arange(n))
              logger.debug('Returning lower triangular matrix {!r}.'.format(L))
              return (L,)
          else:
              logger.debug('Returning permutation matrix {!r} and lower triangular matrix {!r}.'.format(P, L))
              return (P, L)
  
      ## return P, L, D
      if return_type in (RETURN_L_D, RETURN_P_L_D):
          L, D = f.L_D()
          # Do not use f.L_D() -> higher memory consumption
          # LD = f.LD()
          if return_type == RETURN_L_D:
              logger.debug('Returning lower triangular matrix {!r} and diagonal matrix {!r}.'.format(P, L, D))
              return (L, D)
          else:
              logger.debug('Returning permutation matrix {!r}, lower triangular matrix {!r} and diagonal matrix {!r}.'.format(P, L, D))
              return (P, L, D)