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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dsytrf_(char *uplo, integer *n, doublereal *a, integer *
   4         lda, integer *ipiv, doublereal *work, integer *lwork, integer *info)
   5 {
   6 /*  -- LAPACK routine (version 3.0) --
   7        Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
   8        Courant Institute, Argonne National Lab, and Rice University
   9        June 30, 1999
  10
  11
  12     Purpose
  13     =======
  14
  15     DSYTRF computes the factorization of a real symmetric matrix A using
  16     the Bunch-Kaufman diagonal pivoting method.  The form of the
  17     factorization is
  18
  19        A = U*D*U**T  or  A = L*D*L**T
  20
  21     where U (or L) is a product of permutation and unit upper (lower)
  22     triangular matrices, and D is symmetric and block diagonal with
  23     1-by-1 and 2-by-2 diagonal blocks.
  24
  25     This is the blocked version of the algorithm, calling Level 3 BLAS.
  26
  27     Arguments
  28     =========
  29
  30     UPLO    (input) CHARACTER*1
  31             = 'U':  Upper triangle of A is stored;
  32             = 'L':  Lower triangle of A is stored.
  33
  34     N       (input) INTEGER
  35             The order of the matrix A.  N >= 0.
  36
  37     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  38             On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  39             N-by-N upper triangular part of A contains the upper
  40             triangular part of the matrix A, and the strictly lower
  41             triangular part of A is not referenced.  If UPLO = 'L', the
  42             leading N-by-N lower triangular part of A contains the lower
  43             triangular part of the matrix A, and the strictly upper
  44             triangular part of A is not referenced.
  45
  46             On exit, the block diagonal matrix D and the multipliers used
  47             to obtain the factor U or L (see below for further details).
  48
  49     LDA     (input) INTEGER
  50             The leading dimension of the array A.  LDA >= max(1,N).
  51
  52     IPIV    (output) INTEGER array, dimension (N)
  53             Details of the interchanges and the block structure of D.
  54             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  55             interchanged and D(k,k) is a 1-by-1 diagonal block.
  56             If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
  57             columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  58             is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
  59             IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
  60             interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  61
  62     WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  63             On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  64
  65     LWORK   (input) INTEGER
  66             The length of WORK.  LWORK >=1.  For best performance
  67             LWORK >= N*NB, where NB is the block size returned by ILAENV.
  68
  69             If LWORK = -1, then a workspace query is assumed; the routine
  70             only calculates the optimal size of the WORK array, returns
  71             this value as the first entry of the WORK array, and no error
  72             message related to LWORK is issued by XERBLA.
  73
  74     INFO    (output) INTEGER
  75             = 0:  successful exit
  76             < 0:  if INFO = -i, the i-th argument had an illegal value
  77             > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
  78                   has been completed, but the block diagonal matrix D is
  79                   exactly singular, and division by zero will occur if it
  80                   is used to solve a system of equations.
  81
  82     Further Details
  83     ===============
  84
  85     If UPLO = 'U', then A = U*D*U', where
  86        U = P(n)*U(n)* ... *P(k)U(k)* ...,
  87     i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  88     1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  89     and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  90     defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  91     that if the diagonal block D(k) is of order s (s = 1 or 2), then
  92
  93                (   I    v    0   )   k-s
  94        U(k) =  (   0    I    0   )   s
  95                (   0    0    I   )   n-k
  96                   k-s   s   n-k
  97
  98     If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
  99     If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
 100     and A(k,k), and v overwrites A(1:k-2,k-1:k).
 101
 102     If UPLO = 'L', then A = L*D*L', where
 103        L = P(1)*L(1)* ... *P(k)*L(k)* ...,
 104     i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
 105     n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 106     and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 107     defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
 108     that if the diagonal block D(k) is of order s (s = 1 or 2), then
 109
 110                (   I    0     0   )  k-1
 111        L(k) =  (   0    I     0   )  s
 112                (   0    v     I   )  n-k-s+1
 113                   k-1   s  n-k-s+1
 114
 115     If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
 116     If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
 117     and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
 118
 119     =====================================================================
 120
 121
 122        Test the input parameters.
 123
 124        Parameter adjustments */
 125     /* Table of constant values */
 126     static integer c__1 = 1;
 127     static integer c_n1 = -1;
 128     static integer c__2 = 2;
 129
 130     /* System generated locals */
 131     integer a_dim1, a_offset, i__1, i__2;
 132     /* Local variables */
 133     static integer j, k;
 134     extern logical lsame_(char *, char *);
 135     static integer nbmin, iinfo;
 136     static logical upper;
 137     extern /* Subroutine */ int dsytf2_(char *, integer *, doublereal *,
 138             integer *, integer *, integer *);
 139     static integer kb, nb;
 140     extern /* Subroutine */ int xerbla_(char *, integer *);
 141     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
 142             integer *, integer *, ftnlen, ftnlen);
 143     extern /* Subroutine */ int dlasyf_(char *, integer *, integer *, integer
 144             *, doublereal *, integer *, integer *, doublereal *, integer *,
 145             integer *);
 146     static integer ldwork, lwkopt;
 147     static logical lquery;
 148     static integer iws;
 149 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
 150
 151
 152     a_dim1 = *lda;
 153     a_offset = 1 + a_dim1 * 1;
 154     a -= a_offset;
 155     --ipiv;
 156     --work;
 157
 158     /* Function Body */
 159     *info = 0;
 160     upper = lsame_(uplo, "U");
 161     lquery = *lwork == -1;
 162     if (! upper && ! lsame_(uplo, "L")) {
 163         *info = -1;
 164     } else if (*n < 0) {
 165         *info = -2;
 166     } else if (*lda < max(1,*n)) {
 167         *info = -4;
 168     } else if (*lwork < 1 && ! lquery) {
 169         *info = -7;
 170     }
 171
 172     if (*info == 0) {
 173
 174 /*        Determine the block size */
 175
 176         nb = ilaenv_(&c__1, "DSYTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
 177                  (ftnlen)1);
 178         lwkopt = *n * nb;
 179         work[1] = (doublereal) lwkopt;
 180     }
 181
 182     if (*info != 0) {
 183         i__1 = -(*info);
 184         xerbla_("DSYTRF", &i__1);
 185         return 0;
 186     } else if (lquery) {
 187         return 0;
 188     }
 189
 190     nbmin = 2;
 191     ldwork = *n;
 192     if (nb > 1 && nb < *n) {
 193         iws = ldwork * nb;
 194         if (*lwork < iws) {
 195 /* Computing MAX */
 196             i__1 = *lwork / ldwork;
 197             nb = max(i__1,1);
 198 /* Computing MAX */
 199             i__1 = 2, i__2 = ilaenv_(&c__2, "DSYTRF", uplo, n, &c_n1, &c_n1, &
 200                     c_n1, (ftnlen)6, (ftnlen)1);
 201             nbmin = max(i__1,i__2);
 202         }
 203     } else {
 204         iws = 1;
 205     }
 206     if (nb < nbmin) {
 207         nb = *n;
 208     }
 209
 210     if (upper) {
 211
 212 /*        Factorize A as U*D*U' using the upper triangle of A
 213
 214           K is the main loop index, decreasing from N to 1 in steps of
 215           KB, where KB is the number of columns factorized by DLASYF;
 216           KB is either NB or NB-1, or K for the last block */
 217
 218         k = *n;
 219 L10:
 220
 221 /*        If K < 1, exit from loop */
 222
 223         if (k < 1) {
 224             goto L40;
 225         }
 226
 227         if (k > nb) {
 228
 229 /*           Factorize columns k-kb+1:k of A and use blocked code to
 230              update columns 1:k-kb */
 231
 232             dlasyf_(uplo, &k, &nb, &kb, &a[a_offset], lda, &ipiv[1], &work[1],
 233                      &ldwork, &iinfo);
 234         } else {
 235
 236 /*           Use unblocked code to factorize columns 1:k of A */
 237
 238             dsytf2_(uplo, &k, &a[a_offset], lda, &ipiv[1], &iinfo);
 239             kb = k;
 240         }
 241
 242 /*        Set INFO on the first occurrence of a zero pivot */
 243
 244         if (*info == 0 && iinfo > 0) {
 245             *info = iinfo;
 246         }
 247
 248 /*        Decrease K and return to the start of the main loop */
 249
 250         k -= kb;
 251         goto L10;
 252
 253     } else {
 254
 255 /*        Factorize A as L*D*L' using the lower triangle of A
 256
 257           K is the main loop index, increasing from 1 to N in steps of
 258           KB, where KB is the number of columns factorized by DLASYF;
 259           KB is either NB or NB-1, or N-K+1 for the last block */
 260
 261         k = 1;
 262 L20:
 263
 264 /*        If K > N, exit from loop */
 265
 266         if (k > *n) {
 267             goto L40;
 268         }
 269
 270         if (k <= *n - nb) {
 271
 272 /*           Factorize columns k:k+kb-1 of A and use blocked code to
 273              update columns k+kb:n */
 274
 275             i__1 = *n - k + 1;
 276             dlasyf_(uplo, &i__1, &nb, &kb, &a_ref(k, k), lda, &ipiv[k], &work[
 277                     1], &ldwork, &iinfo);
 278         } else {
 279
 280 /*           Use unblocked code to factorize columns k:n of A */
 281
 282             i__1 = *n - k + 1;
 283             dsytf2_(uplo, &i__1, &a_ref(k, k), lda, &ipiv[k], &iinfo);
 284             kb = *n - k + 1;
 285         }
 286
 287 /*        Set INFO on the first occurrence of a zero pivot */
 288
 289         if (*info == 0 && iinfo > 0) {
 290             *info = iinfo + k - 1;
 291         }
 292
 293 /*        Adjust IPIV */
 294
 295         i__1 = k + kb - 1;
 296         for (j = k; j <= i__1; ++j) {
 297             if (ipiv[j] > 0) {
 298                 ipiv[j] = ipiv[j] + k - 1;
 299             } else {
 300                 ipiv[j] = ipiv[j] - k + 1;
 301             }
 302 /* L30: */
 303         }
 304
 305 /*        Increase K and return to the start of the main loop */
 306
 307         k += kb;
 308         goto L20;
 309
 310     }
 311
 312 L40:
 313     work[1] = (doublereal) lwkopt;
 314     return 0;
 315
 316 /*     End of DSYTRF */
 317
 318 } /* dsytrf_ */
 319
 320 #undef a_ref
 321
 322