Fix bug in potrf wrapping: return a triangular matrix

The cholesky decomposition is supposed to be upper triangular. The wrapper for
Lapack's `potrf`, instead, filled the resulting matrix to make it symmetric.
Since this behaviour was introduced in commit 971973d0, which dealt with C89
and added such lines for several wrappers, I wonder if this was not a
copy/paste overlook.

This patch also adds a test for the cholesky decomposition.