--- /dev/null
+#include "clapack.h"
+
+/* Table of constant values */
+
+static doublereal c_b9 = 1.;
+
+/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs,
+ doublereal *a, integer *lda, doublereal *b, integer *ldb, integer *
+ info)
+{
+ /* System generated locals */
+ integer a_dim1, a_offset, b_dim1, b_offset, i__1;
+
+ /* Local variables */
+ extern logical lsame_(char *, char *);
+ extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
+ integer *, integer *, doublereal *, doublereal *, integer *,
+ doublereal *, integer *);
+ logical upper;
+ extern /* Subroutine */ int xerbla_(char *, integer *);
+
+
+/* -- LAPACK routine (version 3.1) -- */
+/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
+/* November 2006 */
+
+/* .. Scalar Arguments .. */
+/* .. */
+/* .. Array Arguments .. */
+/* .. */
+
+/* Purpose */
+/* ======= */
+
+/* DPOTRS solves a system of linear equations A*X = B with a symmetric */
+/* positive definite matrix A using the Cholesky factorization */
+/* A = U**T*U or A = L*L**T computed by DPOTRF. */
+
+/* Arguments */
+/* ========= */
+
+/* UPLO (input) CHARACTER*1 */
+/* = 'U': Upper triangle of A is stored; */
+/* = 'L': Lower triangle of A is stored. */
+
+/* N (input) INTEGER */
+/* The order of the matrix A. N >= 0. */
+
+/* NRHS (input) INTEGER */
+/* The number of right hand sides, i.e., the number of columns */
+/* of the matrix B. NRHS >= 0. */
+
+/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
+/* The triangular factor U or L from the Cholesky factorization */
+/* A = U**T*U or A = L*L**T, as computed by DPOTRF. */
+
+/* LDA (input) INTEGER */
+/* The leading dimension of the array A. LDA >= max(1,N). */
+
+/* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
+/* On entry, the right hand side matrix B. */
+/* On exit, the solution matrix X. */
+
+/* LDB (input) INTEGER */
+/* The leading dimension of the array B. LDB >= max(1,N). */
+
+/* INFO (output) INTEGER */
+/* = 0: successful exit */
+/* < 0: if INFO = -i, the i-th argument had an illegal value */
+
+/* ===================================================================== */
+
+/* .. Parameters .. */
+/* .. */
+/* .. Local Scalars .. */
+/* .. */
+/* .. External Functions .. */
+/* .. */
+/* .. External Subroutines .. */
+/* .. */
+/* .. Intrinsic Functions .. */
+/* .. */
+/* .. Executable Statements .. */
+
+/* Test the input parameters. */
+
+ /* Parameter adjustments */
+ a_dim1 = *lda;
+ a_offset = 1 + a_dim1;
+ a -= a_offset;
+ b_dim1 = *ldb;
+ b_offset = 1 + b_dim1;
+ b -= b_offset;
+
+ /* Function Body */
+ *info = 0;
+ upper = lsame_(uplo, "U");
+ if (! upper && ! lsame_(uplo, "L")) {
+ *info = -1;
+ } else if (*n < 0) {
+ *info = -2;
+ } else if (*nrhs < 0) {
+ *info = -3;
+ } else if (*lda < max(1,*n)) {
+ *info = -5;
+ } else if (*ldb < max(1,*n)) {
+ *info = -7;
+ }
+ if (*info != 0) {
+ i__1 = -(*info);
+ xerbla_("DPOTRS", &i__1);
+ return 0;
+ }
+
+/* Quick return if possible */
+
+ if (*n == 0 || *nrhs == 0) {
+ return 0;
+ }
+
+ if (upper) {
+
+/* Solve A*X = B where A = U'*U. */
+
+/* Solve U'*X = B, overwriting B with X. */
+
+ dtrsm_("Left", "Upper", "Transpose", "Non-unit", n, nrhs, &c_b9, &a[
+ a_offset], lda, &b[b_offset], ldb);
+
+/* Solve U*X = B, overwriting B with X. */
+
+ dtrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b9, &
+ a[a_offset], lda, &b[b_offset], ldb);
+ } else {
+
+/* Solve A*X = B where A = L*L'. */
+
+/* Solve L*X = B, overwriting B with X. */
+
+ dtrsm_("Left", "Lower", "No transpose", "Non-unit", n, nrhs, &c_b9, &
+ a[a_offset], lda, &b[b_offset], ldb);
+
+/* Solve L'*X = B, overwriting B with X. */
+
+ dtrsm_("Left", "Lower", "Transpose", "Non-unit", n, nrhs, &c_b9, &a[
+ a_offset], lda, &b[b_offset], ldb);
+ }
+
+ return 0;
+
+/* End of DPOTRS */
+
+} /* dpotrs_ */