git.maemo.org Git - opencv/blob - 3rdparty/lapack/dpotf2.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static doublereal c_b10 = -1.;
   7 static doublereal c_b12 = 1.;
   8
   9 /* Subroutine */ int dpotf2_(char *uplo, integer *n, doublereal *a, integer *
  10         lda, integer *info)
  11 {
  12     /* System generated locals */
  13     integer a_dim1, a_offset, i__1, i__2, i__3;
  14     doublereal d__1;
  15
  16     /* Builtin functions */
  17     double sqrt(doublereal);
  18
  19     /* Local variables */
  20     integer j;
  21     doublereal ajj;
  22     extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
  23             integer *);
  24     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
  25             integer *);
  26     extern logical lsame_(char *, char *);
  27     extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
  28             doublereal *, doublereal *, integer *, doublereal *, integer *,
  29             doublereal *, doublereal *, integer *);
  30     logical upper;
  31     extern /* Subroutine */ int xerbla_(char *, integer *);
  32
  33
  34 /*  -- LAPACK routine (version 3.1) -- */
  35 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  36 /*     November 2006 */
  37
  38 /*     .. Scalar Arguments .. */
  39 /*     .. */
  40 /*     .. Array Arguments .. */
  41 /*     .. */
  42
  43 /*  Purpose */
  44 /*  ======= */
  45
  46 /*  DPOTF2 computes the Cholesky factorization of a real symmetric */
  47 /*  positive definite matrix A. */
  48
  49 /*  The factorization has the form */
  50 /*     A = U' * U ,  if UPLO = 'U', or */
  51 /*     A = L  * L',  if UPLO = 'L', */
  52 /*  where U is an upper triangular matrix and L is lower triangular. */
  53
  54 /*  This is the unblocked version of the algorithm, calling Level 2 BLAS. */
  55
  56 /*  Arguments */
  57 /*  ========= */
  58
  59 /*  UPLO    (input) CHARACTER*1 */
  60 /*          Specifies whether the upper or lower triangular part of the */
  61 /*          symmetric matrix A is stored. */
  62 /*          = 'U':  Upper triangular */
  63 /*          = 'L':  Lower triangular */
  64
  65 /*  N       (input) INTEGER */
  66 /*          The order of the matrix A.  N >= 0. */
  67
  68 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
  69 /*          On entry, the symmetric matrix A.  If UPLO = 'U', the leading */
  70 /*          n by n upper triangular part of A contains the upper */
  71 /*          triangular part of the matrix A, and the strictly lower */
  72 /*          triangular part of A is not referenced.  If UPLO = 'L', the */
  73 /*          leading n by n lower triangular part of A contains the lower */
  74 /*          triangular part of the matrix A, and the strictly upper */
  75 /*          triangular part of A is not referenced. */
  76
  77 /*          On exit, if INFO = 0, the factor U or L from the Cholesky */
  78 /*          factorization A = U'*U  or A = L*L'. */
  79
  80 /*  LDA     (input) INTEGER */
  81 /*          The leading dimension of the array A.  LDA >= max(1,N). */
  82
  83 /*  INFO    (output) INTEGER */
  84 /*          = 0: successful exit */
  85 /*          < 0: if INFO = -k, the k-th argument had an illegal value */
  86 /*          > 0: if INFO = k, the leading minor of order k is not */
  87 /*               positive definite, and the factorization could not be */
  88 /*               completed. */
  89
  90 /*  ===================================================================== */
  91
  92 /*     .. Parameters .. */
  93 /*     .. */
  94 /*     .. Local Scalars .. */
  95 /*     .. */
  96 /*     .. External Functions .. */
  97 /*     .. */
  98 /*     .. External Subroutines .. */
  99 /*     .. */
 100 /*     .. Intrinsic Functions .. */
 101 /*     .. */
 102 /*     .. Executable Statements .. */
 103
 104 /*     Test the input parameters. */
 105
 106     /* Parameter adjustments */
 107     a_dim1 = *lda;
 108     a_offset = 1 + a_dim1;
 109     a -= a_offset;
 110
 111     /* Function Body */
 112     *info = 0;
 113     upper = lsame_(uplo, "U");
 114     if (! upper && ! lsame_(uplo, "L")) {
 115         *info = -1;
 116     } else if (*n < 0) {
 117         *info = -2;
 118     } else if (*lda < max(1,*n)) {
 119         *info = -4;
 120     }
 121     if (*info != 0) {
 122         i__1 = -(*info);
 123         xerbla_("DPOTF2", &i__1);
 124         return 0;
 125     }
 126
 127 /*     Quick return if possible */
 128
 129     if (*n == 0) {
 130         return 0;
 131     }
 132
 133     if (upper) {
 134
 135 /*        Compute the Cholesky factorization A = U'*U. */
 136
 137         i__1 = *n;
 138         for (j = 1; j <= i__1; ++j) {
 139
 140 /*           Compute U(J,J) and test for non-positive-definiteness. */
 141
 142             i__2 = j - 1;
 143             ajj = a[j + j * a_dim1] - ddot_(&i__2, &a[j * a_dim1 + 1], &c__1,
 144                     &a[j * a_dim1 + 1], &c__1);
 145             if (ajj <= 0.) {
 146                 a[j + j * a_dim1] = ajj;
 147                 goto L30;
 148             }
 149             ajj = sqrt(ajj);
 150             a[j + j * a_dim1] = ajj;
 151
 152 /*           Compute elements J+1:N of row J. */
 153
 154             if (j < *n) {
 155                 i__2 = j - 1;
 156                 i__3 = *n - j;
 157                 dgemv_("Transpose", &i__2, &i__3, &c_b10, &a[(j + 1) * a_dim1
 158                         + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b12, &a[j + (
 159                         j + 1) * a_dim1], lda);
 160                 i__2 = *n - j;
 161                 d__1 = 1. / ajj;
 162                 dscal_(&i__2, &d__1, &a[j + (j + 1) * a_dim1], lda);
 163             }
 164 /* L10: */
 165         }
 166     } else {
 167
 168 /*        Compute the Cholesky factorization A = L*L'. */
 169
 170         i__1 = *n;
 171         for (j = 1; j <= i__1; ++j) {
 172
 173 /*           Compute L(J,J) and test for non-positive-definiteness. */
 174
 175             i__2 = j - 1;
 176             ajj = a[j + j * a_dim1] - ddot_(&i__2, &a[j + a_dim1], lda, &a[j
 177                     + a_dim1], lda);
 178             if (ajj <= 0.) {
 179                 a[j + j * a_dim1] = ajj;
 180                 goto L30;
 181             }
 182             ajj = sqrt(ajj);
 183             a[j + j * a_dim1] = ajj;
 184
 185 /*           Compute elements J+1:N of column J. */
 186
 187             if (j < *n) {
 188                 i__2 = *n - j;
 189                 i__3 = j - 1;
 190                 dgemv_("No transpose", &i__2, &i__3, &c_b10, &a[j + 1 +
 191                         a_dim1], lda, &a[j + a_dim1], lda, &c_b12, &a[j + 1 +
 192                         j * a_dim1], &c__1);
 193                 i__2 = *n - j;
 194                 d__1 = 1. / ajj;
 195                 dscal_(&i__2, &d__1, &a[j + 1 + j * a_dim1], &c__1);
 196             }
 197 /* L20: */
 198         }
 199     }
 200     goto L40;
 201
 202 L30:
 203     *info = j;
 204
 205 L40:
 206     return 0;
 207
 208 /*     End of DPOTF2 */
 209
 210 } /* dpotf2_ */