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

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