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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7 static real c_b13 = -1.f;
   8 static real c_b14 = 1.f;
   9
  10 /* Subroutine */ int spotrf_(char *uplo, integer *n, real *a, integer *lda,
  11         integer *info)
  12 {
  13     /* System generated locals */
  14     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
  15
  16     /* Local variables */
  17     integer j, jb, nb;
  18     extern logical lsame_(char *, char *);
  19     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
  20             integer *, real *, real *, integer *, real *, integer *, real *,
  21             real *, integer *);
  22     logical upper;
  23     extern /* Subroutine */ int strsm_(char *, char *, char *, char *,
  24             integer *, integer *, real *, real *, integer *, real *, integer *
  25 ), ssyrk_(char *, char *, integer
  26             *, integer *, real *, real *, integer *, real *, real *, integer *
  27 ), spotf2_(char *, integer *, real *, integer *,
  28             integer *), xerbla_(char *, integer *);
  29     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  30             integer *, integer *);
  31
  32
  33 /*  -- LAPACK routine (version 3.1) -- */
  34 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  35 /*     November 2006 */
  36
  37 /*     .. Scalar Arguments .. */
  38 /*     .. */
  39 /*     .. Array Arguments .. */
  40 /*     .. */
  41
  42 /*  Purpose */
  43 /*  ======= */
  44
  45 /*  SPOTRF computes the Cholesky factorization of a real symmetric */
  46 /*  positive definite matrix A. */
  47
  48 /*  The factorization has the form */
  49 /*     A = U**T * U,  if UPLO = 'U', or */
  50 /*     A = L  * L**T,  if UPLO = 'L', */
  51 /*  where U is an upper triangular matrix and L is lower triangular. */
  52
  53 /*  This is the block version of the algorithm, calling Level 3 BLAS. */
  54
  55 /*  Arguments */
  56 /*  ========= */
  57
  58 /*  UPLO    (input) CHARACTER*1 */
  59 /*          = 'U':  Upper triangle of A is stored; */
  60 /*          = 'L':  Lower triangle of A is stored. */
  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**T*U or A = L*L**T. */
  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 = -i, the i-th argument had an illegal value */
  83 /*          > 0:  if INFO = i, the leading minor of order i 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_("SPOTRF", &i__1);
 121         return 0;
 122     }
 123
 124 /*     Quick return if possible */
 125
 126     if (*n == 0) {
 127         return 0;
 128     }
 129
 130 /*     Determine the block size for this environment. */
 131
 132     nb = ilaenv_(&c__1, "SPOTRF", uplo, n, &c_n1, &c_n1, &c_n1);
 133     if (nb <= 1 || nb >= *n) {
 134
 135 /*        Use unblocked code. */
 136
 137         spotf2_(uplo, n, &a[a_offset], lda, info);
 138     } else {
 139
 140 /*        Use blocked code. */
 141
 142         if (upper) {
 143
 144 /*           Compute the Cholesky factorization A = U'*U. */
 145
 146             i__1 = *n;
 147             i__2 = nb;
 148             for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
 149
 150 /*              Update and factorize the current diagonal block and test */
 151 /*              for non-positive-definiteness. */
 152
 153 /* Computing MIN */
 154                 i__3 = nb, i__4 = *n - j + 1;
 155                 jb = min(i__3,i__4);
 156                 i__3 = j - 1;
 157                 ssyrk_("Upper", "Transpose", &jb, &i__3, &c_b13, &a[j *
 158                         a_dim1 + 1], lda, &c_b14, &a[j + j * a_dim1], lda);
 159                 spotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info);
 160                 if (*info != 0) {
 161                     goto L30;
 162                 }
 163                 if (j + jb <= *n) {
 164
 165 /*                 Compute the current block row. */
 166
 167                     i__3 = *n - j - jb + 1;
 168                     i__4 = j - 1;
 169                     sgemm_("Transpose", "No transpose", &jb, &i__3, &i__4, &
 170                             c_b13, &a[j * a_dim1 + 1], lda, &a[(j + jb) *
 171                             a_dim1 + 1], lda, &c_b14, &a[j + (j + jb) *
 172                             a_dim1], lda);
 173                     i__3 = *n - j - jb + 1;
 174                     strsm_("Left", "Upper", "Transpose", "Non-unit", &jb, &
 175                             i__3, &c_b14, &a[j + j * a_dim1], lda, &a[j + (j
 176                             + jb) * a_dim1], lda);
 177                 }
 178 /* L10: */
 179             }
 180
 181         } else {
 182
 183 /*           Compute the Cholesky factorization A = L*L'. */
 184
 185             i__2 = *n;
 186             i__1 = nb;
 187             for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) {
 188
 189 /*              Update and factorize the current diagonal block and test */
 190 /*              for non-positive-definiteness. */
 191
 192 /* Computing MIN */
 193                 i__3 = nb, i__4 = *n - j + 1;
 194                 jb = min(i__3,i__4);
 195                 i__3 = j - 1;
 196                 ssyrk_("Lower", "No transpose", &jb, &i__3, &c_b13, &a[j +
 197                         a_dim1], lda, &c_b14, &a[j + j * a_dim1], lda);
 198                 spotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info);
 199                 if (*info != 0) {
 200                     goto L30;
 201                 }
 202                 if (j + jb <= *n) {
 203
 204 /*                 Compute the current block column. */
 205
 206                     i__3 = *n - j - jb + 1;
 207                     i__4 = j - 1;
 208                     sgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, &
 209                             c_b13, &a[j + jb + a_dim1], lda, &a[j + a_dim1],
 210                             lda, &c_b14, &a[j + jb + j * a_dim1], lda);
 211                     i__3 = *n - j - jb + 1;
 212                     strsm_("Right", "Lower", "Transpose", "Non-unit", &i__3, &
 213                             jb, &c_b14, &a[j + j * a_dim1], lda, &a[j + jb +
 214                             j * a_dim1], lda);
 215                 }
 216 /* L20: */
 217             }
 218         }
 219     }
 220     goto L40;
 221
 222 L30:
 223     *info = *info + j - 1;
 224
 225 L40:
 226     return 0;
 227
 228 /*     End of SPOTRF */
 229
 230 } /* spotrf_ */