struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
    dim::Int                    # matrix dimension
    mat::S                      # input matrix
    chol::Cholesky{T,S}         # Cholesky factorization of mat
end

# Constructors

PDMat(mat, chol)    # with both the input matrix and a Cholesky factorization

PDMat(mat)          # with the input matrix, of type Matrix or Symmetric
                    # Remarks: the Cholesky factorization will be computed
                    # upon construction.

PDMat(chol)         # with the Cholesky factorization
                    # Remarks: the full matrix will be computed upon
                    # construction.

struct PDiagMat{T<:Real,V<:AbstractVector{T}} <: AbstractPDMat{T}
    dim::Int                    # matrix dimension
    diag::V                     # the vector of diagonal elements
end

# Constructors

PDiagMat(v)         # with the vector of diagonal elements

struct ScalMat{T<:Real} <: AbstractPDMat{T}
    dim::Int         # matrix dimension
    value::T         # diagonal value (shared by all diagonal elements)
end


# Constructors

ScalMat(d, v)        # with dimension d and diagonal value v

struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T}
    dim::Int                       # matrix dimension
    mat::SparseMatrixCSC           # input matrix
    chol::CholTypeSparse           # Cholesky factorization of mat
end

# Constructors

PDSparseMat(mat, chol)    # with both the input matrix and a Cholesky factorization

PDSparseMat(mat)          # with the sparse input matrix, of type SparseMatrixCSC
                          # Remarks: the Cholesky factorization will be computed
                          # upon construction.

PDSparseMat(chol)         # with the Cholesky factorization
                          # Remarks: the sparse matrix 'mat' will be computed upon
                          # construction.

size(a)     # return the size of `a`.

size(a, i)  # return the i-th dimension of `a`.

ndims(a)    # the number of dimensions, which is always 2.

eltype(a)   # the element type

Matrix(a)   # return a copy of the matrix in full form.

diag(a)     # return a vector of diagonal elements.

inv(a)      # inverse of `a`, of a proper subtype of `AbstractPDMat`.
            # Note: when `a` is an instance of either `PDMat`, `PDiagMat`,
            # and `ScalMat`, `inv(a)` is of the same type of `a`.
            # This needs not be required for customized subtypes -- the
            # inverse does not always has the same pattern as `a`.

eigmax(a)   # maximum eigenvalue of `a`.

eigmin(a)   # minimum eigenvalue of `a`.

logdet(a)   # log-determinant of `a`, computed in a numerically stable way.

a * x       # multiple `a` with `x` (forward transform)

a \ x       # multiply `inv(a)` with `x` (backward transform).
            # The internal implementation may not explicitly instantiate
            # the inverse of `a`.

a * c       # scale `a` by a positive scale `c`.
            # The result is in general of the same type of `a`.

c * a       # equivalent to a * c.

a + b       # add two positive definite matrices

pdadd(a, b, c)      # add `a` with `b * c`, where both `a` and `b` are
                    # instances of `AbstractPDMat`.

pdadd(m, a)         # add `a` to a dense matrix `m` of the same size.

pdadd(m, a, c)      # add `a * c` to a dense matrix `m` of the same size.

pdadd!(m, a)        # add `a` to a dense matrix `m` of the same size inplace.

pdadd!(m, a, c)     # add `a * c` to a dense matrix `m` of the same size,
                    # inplace.

pdadd!(r, m, a)     # add `a` to a dense matrix `m` of the same size, and write
                    # the result to `r`.

pdadd!(r, m, a, c)  # add `a * c` to a dense matrix `m` of the same size, and
                    # write the result to `r`.

quad(a, x)          # compute x' * a * x when `x` is a vector.
                    # perform such computation in a column-wise fashion, when
                    # `x` is a matrix, and return a vector of length `n`,
                    # where `n` is the number of columns in `x`.

quad!(r, a, x)      # compute x' * a * x in a column-wise fashion, and write
                    # the results to `r`.

invquad(a, x)       # compute x' * inv(a) * x when `x` is a vector.
                    # perform such computation in a column-wise fashion, when
                    # `x` is a matrix, and return a vector of length `n`.

invquad!(r, a, x)   # compute x' * inv(a) * x in a column-wise fashion, and
                    # write the results to `r`.

X_A_Xt(a, x)        # compute `x * a * x'` for a matrix `x`.

Xt_A_X(a, x)        # compute `x' * a * x` for a matrix `x`.

X_invA_Xt(a, x)     # compute `x * inv(a) * x'` for a matrix `x`.

Xt_invA_X(a, x)     # compute `x' * inv(a) * x` for a matrix `x`.

whiten(a, x)        # whitening transform. `x` can be a vector or a matrix.
                    #
                    # Note: If the covariance of `x` is `a`, then the
                    # covariance of the transformed result is an identity
                    # matrix.

whiten!(a, x)       # whitening transform inplace, directly updating `x`.

whiten!(r, a, x)    # write the transformed result to `r`.

unwhiten(a, x)      # inverse of whitening transform. `x` can be a vector or
                    # a matrix.
                    #
                    # Note: If the covariance of `x` is an identity matrix,
                    # then the covariance of the transformed result is `a`.
                    # Note: the un-whitening transform is useful for
                    # generating Gaussian samples.

unwhiten!(a, x)     # un-whitening transform inplace, updating `x`.

unwhiten!(r, a, x)  # write the transformed result to `r`.

# Let `M` be the name of the subtype, then the following methods need
# to be implemented for `M`:

Matrix(a::M)    # return a copy of the matrix in full form, of type
                # `Matrix{eltype(M)}`.

diag(a::M)      # return the vector of diagonal elements, of type
                # `Vector{eltype(M)}`.

pdadd!(m, a, c)     # add `a * c` to a dense matrix `m` of the same size
                    # inplace.

* (a::M, c::Real)        # return a scaled version of `a`.

* (a::M, x::DenseVecOrMat)  # transform `x`, i.e. compute `a * x`.

\ (a::M, x::DenseVecOrMat)  # inverse transform `x`, i.e. compute `inv(a) * x`.

inv(a::M)       # compute the inverse of `a`.

logdet(a::M)    # compute the log-determinant of `a`.

eigmax(a::M)    # compute the maximum eigenvalue of `a`.

eigmin(a::M)    # compute the minimum eigenvalue of `a`.

whiten!(r::DenseVecOrMat, a::M, x::DenseVecOrMat)  # whitening transform,
                                                   # write result to `r`.

unwhiten!(r::DenseVecOrMat, a::M, x::DenseVecOrMat)  # un-whitening transform,
                                                     # write result to `r`.

quad(a::M, x::DenseVector)      # compute `x' * a * x`

quad!(r::AbstractArray, a::M, x::DenseMatrix)   # compute `x' * a * x` in
                                                # a column-wise manner

invquad(a::M, x::DenseVector)   # compute `x' * inv(a) * x`

invquad!(r::AbstractArray, a::M, x::DenseMatrix) # compute `x' * inv(a) * x`
                                                 # in a column-wise manner

X_A_Xt(a::M, x::DenseMatrix)        # compute `x * a * x'`

Xt_A_X(a::M, x::DenseMatrix)        # compute `x' * a * x`

X_invA_Xt(a::M, x::DenseMatrix)     # compute `x * inv(a) * x'`

Xt_invA_X(a::M, x::DenseMatrix)     # compute `x' * inv(a) * x`