Added cholesky decomp routines.
This commit is contained in:
Harry Moffat
parent
701dd1e1cd
commit
ac660e53c0
6 changed files with 746 additions and 0 deletions
|
|
@ -100,6 +100,8 @@ dormbr.o \
|
|||
dorml2.o \
|
||||
dormlq.o \
|
||||
dormqr.o \
|
||||
dpotrf.o \
|
||||
dpotrs.o \
|
||||
drscl.o \
|
||||
dtrcon.o \
|
||||
dtrtri.o \
|
||||
|
|
|
|||
254
ext/f2c_lapack/dpotrf.c
Normal file
254
ext/f2c_lapack/dpotrf.c
Normal file
|
|
@ -0,0 +1,254 @@
|
|||
/* dpotrf.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#include "f2c.h"
|
||||
|
||||
/* Table of constant values */
|
||||
|
||||
static integer c__1 = 1;
|
||||
static integer c_n1 = -1;
|
||||
static doublereal c_b13 = -1.;
|
||||
static doublereal c_b14 = 1.;
|
||||
|
||||
/* Subroutine */ int dpotrf_(char *uplo, integer *n, doublereal *a, integer *
|
||||
lda, integer *info, ftnlen uplo_len)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
|
||||
|
||||
/* Local variables */
|
||||
static integer j, jb, nb;
|
||||
extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
|
||||
integer *, doublereal *, doublereal *, integer *, doublereal *,
|
||||
integer *, doublereal *, doublereal *, integer *, ftnlen, ftnlen);
|
||||
extern logical lsame_(char *, char *, ftnlen, ftnlen);
|
||||
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublereal *, doublereal *, integer *,
|
||||
doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen);
|
||||
static logical upper;
|
||||
extern /* Subroutine */ int dsyrk_(char *, char *, integer *, integer *,
|
||||
doublereal *, doublereal *, integer *, doublereal *, doublereal *,
|
||||
integer *, ftnlen, ftnlen), dpotf2_(char *, integer *,
|
||||
doublereal *, integer *, integer *, ftnlen), xerbla_(char *,
|
||||
integer *, ftnlen);
|
||||
extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
|
||||
integer *, integer *, ftnlen, ftnlen);
|
||||
|
||||
|
||||
/* -- LAPACK routine (version 3.0) -- */
|
||||
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
|
||||
/* Courant Institute, Argonne National Lab, and Rice University */
|
||||
/* March 31, 1993 */
|
||||
|
||||
/* .. Scalar Arguments .. */
|
||||
/* .. */
|
||||
/* .. Array Arguments .. */
|
||||
/* .. */
|
||||
|
||||
/* Purpose */
|
||||
/* ======= */
|
||||
|
||||
/* DPOTRF computes the Cholesky factorization of a real symmetric */
|
||||
/* positive definite matrix A. */
|
||||
|
||||
/* The factorization has the form */
|
||||
/* A = U**T * U, if UPLO = 'U', or */
|
||||
/* A = L * L**T, if UPLO = 'L', */
|
||||
/* where U is an upper triangular matrix and L is lower triangular. */
|
||||
|
||||
/* This is the block version of the algorithm, calling Level 3 BLAS. */
|
||||
|
||||
/* 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. */
|
||||
|
||||
/* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
|
||||
/* On entry, the symmetric matrix A. If UPLO = 'U', the leading */
|
||||
/* N-by-N upper triangular part of A contains the upper */
|
||||
/* triangular part of the matrix A, and the strictly lower */
|
||||
/* triangular part of A is not referenced. If UPLO = 'L', the */
|
||||
/* leading N-by-N lower triangular part of A contains the lower */
|
||||
/* triangular part of the matrix A, and the strictly upper */
|
||||
/* triangular part of A is not referenced. */
|
||||
|
||||
/* On exit, if INFO = 0, the factor U or L from the Cholesky */
|
||||
/* factorization A = U**T*U or A = L*L**T. */
|
||||
|
||||
/* LDA (input) INTEGER */
|
||||
/* The leading dimension of the array A. LDA >= max(1,N). */
|
||||
|
||||
/* INFO (output) INTEGER */
|
||||
/* = 0: successful exit */
|
||||
/* < 0: if INFO = -i, the i-th argument had an illegal value */
|
||||
/* > 0: if INFO = i, the leading minor of order i is not */
|
||||
/* positive definite, and the factorization could not be */
|
||||
/* completed. */
|
||||
|
||||
/* ===================================================================== */
|
||||
|
||||
/* .. 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;
|
||||
|
||||
/* Function Body */
|
||||
*info = 0;
|
||||
upper = lsame_(uplo, "U", (ftnlen)1, (ftnlen)1);
|
||||
if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
|
||||
*info = -1;
|
||||
} else if (*n < 0) {
|
||||
*info = -2;
|
||||
} else if (*lda < max(1,*n)) {
|
||||
*info = -4;
|
||||
}
|
||||
if (*info != 0) {
|
||||
i__1 = -(*info);
|
||||
xerbla_("DPOTRF", &i__1, (ftnlen)6);
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Quick return if possible */
|
||||
|
||||
if (*n == 0) {
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* Determine the block size for this environment. */
|
||||
|
||||
nb = ilaenv_(&c__1, "DPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, (
|
||||
ftnlen)1);
|
||||
if (nb <= 1 || nb >= *n) {
|
||||
|
||||
/* Use unblocked code. */
|
||||
|
||||
dpotf2_(uplo, n, &a[a_offset], lda, info, (ftnlen)1);
|
||||
} else {
|
||||
|
||||
/* Use blocked code. */
|
||||
|
||||
if (upper) {
|
||||
|
||||
/* Compute the Cholesky factorization A = U'*U. */
|
||||
|
||||
i__1 = *n;
|
||||
i__2 = nb;
|
||||
for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
|
||||
|
||||
/* Update and factorize the current diagonal block and test */
|
||||
/* for non-positive-definiteness. */
|
||||
|
||||
/* Computing MIN */
|
||||
i__3 = nb, i__4 = *n - j + 1;
|
||||
jb = min(i__3,i__4);
|
||||
i__3 = j - 1;
|
||||
dsyrk_("Upper", "Transpose", &jb, &i__3, &c_b13, &a[j *
|
||||
a_dim1 + 1], lda, &c_b14, &a[j + j * a_dim1], lda, (
|
||||
ftnlen)5, (ftnlen)9);
|
||||
dpotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info, (ftnlen)
|
||||
5);
|
||||
if (*info != 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (j + jb <= *n) {
|
||||
|
||||
/* Compute the current block row. */
|
||||
|
||||
i__3 = *n - j - jb + 1;
|
||||
i__4 = j - 1;
|
||||
dgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, &
|
||||
c_b13, &a[j * a_dim1 + 1], lda, &a[(j + jb) *
|
||||
a_dim1 + 1], lda, &c_b14, &a[j + (j + jb) *
|
||||
a_dim1], lda, (ftnlen)9, (ftnlen)12);
|
||||
i__3 = *n - j - jb + 1;
|
||||
dtrsm_("Left", "Upper", "Transpose", "Non-unit", &jb, &
|
||||
i__3, &c_b14, &a[j + j * a_dim1], lda, &a[j + (j
|
||||
+ jb) * a_dim1], lda, (ftnlen)4, (ftnlen)5, (
|
||||
ftnlen)9, (ftnlen)8);
|
||||
}
|
||||
/* L10: */
|
||||
}
|
||||
|
||||
} else {
|
||||
|
||||
/* Compute the Cholesky factorization A = L*L'. */
|
||||
|
||||
i__2 = *n;
|
||||
i__1 = nb;
|
||||
for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
|
||||
|
||||
/* Update and factorize the current diagonal block and test */
|
||||
/* for non-positive-definiteness. */
|
||||
|
||||
/* Computing MIN */
|
||||
i__3 = nb, i__4 = *n - j + 1;
|
||||
jb = min(i__3,i__4);
|
||||
i__3 = j - 1;
|
||||
dsyrk_("Lower", "No transpose", &jb, &i__3, &c_b13, &a[j +
|
||||
a_dim1], lda, &c_b14, &a[j + j * a_dim1], lda, (
|
||||
ftnlen)5, (ftnlen)12);
|
||||
dpotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info, (ftnlen)
|
||||
5);
|
||||
if (*info != 0) {
|
||||
goto L30;
|
||||
}
|
||||
if (j + jb <= *n) {
|
||||
|
||||
/* Compute the current block column. */
|
||||
|
||||
i__3 = *n - j - jb + 1;
|
||||
i__4 = j - 1;
|
||||
dgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, &
|
||||
c_b13, &a[j + jb + a_dim1], lda, &a[j + a_dim1],
|
||||
lda, &c_b14, &a[j + jb + j * a_dim1], lda, (
|
||||
ftnlen)12, (ftnlen)9);
|
||||
i__3 = *n - j - jb + 1;
|
||||
dtrsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, &
|
||||
jb, &c_b14, &a[j + j * a_dim1], lda, &a[j + jb +
|
||||
j * a_dim1], lda, (ftnlen)5, (ftnlen)5, (ftnlen)9,
|
||||
(ftnlen)8);
|
||||
}
|
||||
/* L20: */
|
||||
}
|
||||
}
|
||||
}
|
||||
goto L40;
|
||||
|
||||
L30:
|
||||
*info = *info + j - 1;
|
||||
|
||||
L40:
|
||||
return 0;
|
||||
|
||||
/* End of DPOTRF */
|
||||
|
||||
} /* dpotrf_ */
|
||||
|
||||
171
ext/f2c_lapack/dpotrs.c
Normal file
171
ext/f2c_lapack/dpotrs.c
Normal file
|
|
@ -0,0 +1,171 @@
|
|||
/* dpotrs.f -- translated by f2c (version 20031025).
|
||||
You must link the resulting object file with libf2c:
|
||||
on Microsoft Windows system, link with libf2c.lib;
|
||||
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
|
||||
or, if you install libf2c.a in a standard place, with -lf2c -lm
|
||||
-- in that order, at the end of the command line, as in
|
||||
cc *.o -lf2c -lm
|
||||
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
|
||||
|
||||
http://www.netlib.org/f2c/libf2c.zip
|
||||
*/
|
||||
|
||||
#include "f2c.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, ftnlen uplo_len)
|
||||
{
|
||||
/* System generated locals */
|
||||
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
|
||||
|
||||
/* Local variables */
|
||||
extern logical lsame_(char *, char *, ftnlen, ftnlen);
|
||||
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
|
||||
integer *, integer *, doublereal *, doublereal *, integer *,
|
||||
doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen);
|
||||
static logical upper;
|
||||
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
|
||||
|
||||
|
||||
/* -- LAPACK routine (version 3.0) -- */
|
||||
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
|
||||
/* Courant Institute, Argonne National Lab, and Rice University */
|
||||
/* March 31, 1993 */
|
||||
|
||||
/* .. 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", (ftnlen)1, (ftnlen)1);
|
||||
if (! upper && ! lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
|
||||
*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, (ftnlen)6);
|
||||
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, (ftnlen)4, (ftnlen)5, (
|
||||
ftnlen)9, (ftnlen)8);
|
||||
|
||||
/* 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, (ftnlen)4, (ftnlen)5, (
|
||||
ftnlen)12, (ftnlen)8);
|
||||
} 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, (ftnlen)4, (ftnlen)5, (
|
||||
ftnlen)12, (ftnlen)8);
|
||||
|
||||
/* 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, (ftnlen)4, (ftnlen)5, (
|
||||
ftnlen)9, (ftnlen)8);
|
||||
}
|
||||
|
||||
return 0;
|
||||
|
||||
/* End of DPOTRS */
|
||||
|
||||
} /* dpotrs_ */
|
||||
|
||||
|
|
@ -69,6 +69,8 @@ dormbr.o \
|
|||
dorml2.o \
|
||||
dormlq.o \
|
||||
dormqr.o \
|
||||
dpotrf.o \
|
||||
dpotrs.o \
|
||||
drscl.o \
|
||||
dtrcon.o \
|
||||
dtrtrs.o \
|
||||
|
|
|
|||
184
ext/lapack/dpotrf.f
Normal file
184
ext/lapack/dpotrf.f
Normal file
|
|
@ -0,0 +1,184 @@
|
|||
SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
|
||||
*
|
||||
* -- LAPACK routine (version 3.0) --
|
||||
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
||||
* Courant Institute, Argonne National Lab, and Rice University
|
||||
* March 31, 1993
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER UPLO
|
||||
INTEGER INFO, LDA, N
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION A( LDA, * )
|
||||
* ..
|
||||
*
|
||||
* Purpose
|
||||
* =======
|
||||
*
|
||||
* DPOTRF computes the Cholesky factorization of a real symmetric
|
||||
* positive definite matrix A.
|
||||
*
|
||||
* The factorization has the form
|
||||
* A = U**T * U, if UPLO = 'U', or
|
||||
* A = L * L**T, if UPLO = 'L',
|
||||
* where U is an upper triangular matrix and L is lower triangular.
|
||||
*
|
||||
* This is the block version of the algorithm, calling Level 3 BLAS.
|
||||
*
|
||||
* 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.
|
||||
*
|
||||
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
|
||||
* On entry, the symmetric matrix A. If UPLO = 'U', the leading
|
||||
* N-by-N upper triangular part of A contains the upper
|
||||
* triangular part of the matrix A, and the strictly lower
|
||||
* triangular part of A is not referenced. If UPLO = 'L', the
|
||||
* leading N-by-N lower triangular part of A contains the lower
|
||||
* triangular part of the matrix A, and the strictly upper
|
||||
* triangular part of A is not referenced.
|
||||
*
|
||||
* On exit, if INFO = 0, the factor U or L from the Cholesky
|
||||
* factorization A = U**T*U or A = L*L**T.
|
||||
*
|
||||
* LDA (input) INTEGER
|
||||
* The leading dimension of the array A. LDA >= max(1,N).
|
||||
*
|
||||
* INFO (output) INTEGER
|
||||
* = 0: successful exit
|
||||
* < 0: if INFO = -i, the i-th argument had an illegal value
|
||||
* > 0: if INFO = i, the leading minor of order i is not
|
||||
* positive definite, and the factorization could not be
|
||||
* completed.
|
||||
*
|
||||
* =====================================================================
|
||||
*
|
||||
* .. Parameters ..
|
||||
DOUBLE PRECISION ONE
|
||||
PARAMETER ( ONE = 1.0D+0 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL UPPER
|
||||
INTEGER J, JB, NB
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
INTEGER ILAENV
|
||||
EXTERNAL LSAME, ILAENV
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX, MIN
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
UPPER = LSAME( UPLO, 'U' )
|
||||
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
||||
INFO = -1
|
||||
ELSE IF( N.LT.0 ) THEN
|
||||
INFO = -2
|
||||
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
||||
INFO = -4
|
||||
END IF
|
||||
IF( INFO.NE.0 ) THEN
|
||||
CALL XERBLA( 'DPOTRF', -INFO )
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
IF( N.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
* Determine the block size for this environment.
|
||||
*
|
||||
NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
|
||||
IF( NB.LE.1 .OR. NB.GE.N ) THEN
|
||||
*
|
||||
* Use unblocked code.
|
||||
*
|
||||
CALL DPOTF2( UPLO, N, A, LDA, INFO )
|
||||
ELSE
|
||||
*
|
||||
* Use blocked code.
|
||||
*
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* Compute the Cholesky factorization A = U'*U.
|
||||
*
|
||||
DO 10 J = 1, N, NB
|
||||
*
|
||||
* Update and factorize the current diagonal block and test
|
||||
* for non-positive-definiteness.
|
||||
*
|
||||
JB = MIN( NB, N-J+1 )
|
||||
CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE,
|
||||
$ A( 1, J ), LDA, ONE, A( J, J ), LDA )
|
||||
CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
|
||||
IF( INFO.NE.0 )
|
||||
$ GO TO 30
|
||||
IF( J+JB.LE.N ) THEN
|
||||
*
|
||||
* Compute the current block row.
|
||||
*
|
||||
CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1,
|
||||
$ J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ),
|
||||
$ LDA, ONE, A( J, J+JB ), LDA )
|
||||
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit',
|
||||
$ JB, N-J-JB+1, ONE, A( J, J ), LDA,
|
||||
$ A( J, J+JB ), LDA )
|
||||
END IF
|
||||
10 CONTINUE
|
||||
*
|
||||
ELSE
|
||||
*
|
||||
* Compute the Cholesky factorization A = L*L'.
|
||||
*
|
||||
DO 20 J = 1, N, NB
|
||||
*
|
||||
* Update and factorize the current diagonal block and test
|
||||
* for non-positive-definiteness.
|
||||
*
|
||||
JB = MIN( NB, N-J+1 )
|
||||
CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE,
|
||||
$ A( J, 1 ), LDA, ONE, A( J, J ), LDA )
|
||||
CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
|
||||
IF( INFO.NE.0 )
|
||||
$ GO TO 30
|
||||
IF( J+JB.LE.N ) THEN
|
||||
*
|
||||
* Compute the current block column.
|
||||
*
|
||||
CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
|
||||
$ J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ),
|
||||
$ LDA, ONE, A( J+JB, J ), LDA )
|
||||
CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit',
|
||||
$ N-J-JB+1, JB, ONE, A( J, J ), LDA,
|
||||
$ A( J+JB, J ), LDA )
|
||||
END IF
|
||||
20 CONTINUE
|
||||
END IF
|
||||
END IF
|
||||
GO TO 40
|
||||
*
|
||||
30 CONTINUE
|
||||
INFO = INFO + J - 1
|
||||
*
|
||||
40 CONTINUE
|
||||
RETURN
|
||||
*
|
||||
* End of DPOTRF
|
||||
*
|
||||
END
|
||||
133
ext/lapack/dpotrs.f
Normal file
133
ext/lapack/dpotrs.f
Normal file
|
|
@ -0,0 +1,133 @@
|
|||
SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
|
||||
*
|
||||
* -- LAPACK routine (version 3.0) --
|
||||
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
|
||||
* Courant Institute, Argonne National Lab, and Rice University
|
||||
* March 31, 1993
|
||||
*
|
||||
* .. Scalar Arguments ..
|
||||
CHARACTER UPLO
|
||||
INTEGER INFO, LDA, LDB, N, NRHS
|
||||
* ..
|
||||
* .. Array Arguments ..
|
||||
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
|
||||
* ..
|
||||
*
|
||||
* 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 ..
|
||||
DOUBLE PRECISION ONE
|
||||
PARAMETER ( ONE = 1.0D+0 )
|
||||
* ..
|
||||
* .. Local Scalars ..
|
||||
LOGICAL UPPER
|
||||
* ..
|
||||
* .. External Functions ..
|
||||
LOGICAL LSAME
|
||||
EXTERNAL LSAME
|
||||
* ..
|
||||
* .. External Subroutines ..
|
||||
EXTERNAL DTRSM, XERBLA
|
||||
* ..
|
||||
* .. Intrinsic Functions ..
|
||||
INTRINSIC MAX
|
||||
* ..
|
||||
* .. Executable Statements ..
|
||||
*
|
||||
* Test the input parameters.
|
||||
*
|
||||
INFO = 0
|
||||
UPPER = LSAME( UPLO, 'U' )
|
||||
IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
|
||||
INFO = -1
|
||||
ELSE IF( N.LT.0 ) THEN
|
||||
INFO = -2
|
||||
ELSE IF( NRHS.LT.0 ) THEN
|
||||
INFO = -3
|
||||
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
|
||||
INFO = -5
|
||||
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
|
||||
INFO = -7
|
||||
END IF
|
||||
IF( INFO.NE.0 ) THEN
|
||||
CALL XERBLA( 'DPOTRS', -INFO )
|
||||
RETURN
|
||||
END IF
|
||||
*
|
||||
* Quick return if possible
|
||||
*
|
||||
IF( N.EQ.0 .OR. NRHS.EQ.0 )
|
||||
$ RETURN
|
||||
*
|
||||
IF( UPPER ) THEN
|
||||
*
|
||||
* Solve A*X = B where A = U'*U.
|
||||
*
|
||||
* Solve U'*X = B, overwriting B with X.
|
||||
*
|
||||
CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
|
||||
$ ONE, A, LDA, B, LDB )
|
||||
*
|
||||
* Solve U*X = B, overwriting B with X.
|
||||
*
|
||||
CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
|
||||
$ NRHS, ONE, A, LDA, B, LDB )
|
||||
ELSE
|
||||
*
|
||||
* Solve A*X = B where A = L*L'.
|
||||
*
|
||||
* Solve L*X = B, overwriting B with X.
|
||||
*
|
||||
CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
|
||||
$ NRHS, ONE, A, LDA, B, LDB )
|
||||
*
|
||||
* Solve L'*X = B, overwriting B with X.
|
||||
*
|
||||
CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS,
|
||||
$ ONE, A, LDA, B, LDB )
|
||||
END IF
|
||||
*
|
||||
RETURN
|
||||
*
|
||||
* End of DPOTRS
|
||||
*
|
||||
END
|
||||
Loading…
Add table
Reference in a new issue