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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dsytf2_(char *uplo, integer *n, doublereal *a, integer *
   4         lda, integer *ipiv, integer *info)
   5 {
   6 /*  -- LAPACK routine (version 3.1) --
   7        Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   8        November 2006
   9
  10
  11     Purpose
  12     =======
  13
  14     DSYTF2 computes the factorization of a real symmetric matrix A using
  15     the Bunch-Kaufman diagonal pivoting method:
  16
  17        A = U*D*U'  or  A = L*D*L'
  18
  19     where U (or L) is a product of permutation and unit upper (lower)
  20     triangular matrices, U' is the transpose of U, and D is symmetric and
  21     block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
  22
  23     This is the unblocked version of the algorithm, calling Level 2 BLAS.
  24
  25     Arguments
  26     =========
  27
  28     UPLO    (input) CHARACTER*1
  29             Specifies whether the upper or lower triangular part of the
  30             symmetric matrix A is stored:
  31             = 'U':  Upper triangular
  32             = 'L':  Lower triangular
  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     INFO    (output) INTEGER
  63             = 0: successful exit
  64             < 0: if INFO = -k, the k-th argument had an illegal value
  65             > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
  66                  has been completed, but the block diagonal matrix D is
  67                  exactly singular, and division by zero will occur if it
  68                  is used to solve a system of equations.
  69
  70     Further Details
  71     ===============
  72
  73     09-29-06 - patch from
  74       Bobby Cheng, MathWorks
  75
  76       Replace l.204 and l.372
  77            IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
  78       by
  79            IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
  80
  81     01-01-96 - Based on modifications by
  82       J. Lewis, Boeing Computer Services Company
  83       A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
  84     1-96 - Based on modifications by J. Lewis, Boeing Computer Services
  85            Company
  86
  87     If UPLO = 'U', then A = U*D*U', where
  88        U = P(n)*U(n)* ... *P(k)U(k)* ...,
  89     i.e., U is a product of terms P(k)*U(k), where k decreases from n to
  90     1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
  91     and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
  92     defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
  93     that if the diagonal block D(k) is of order s (s = 1 or 2), then
  94
  95                (   I    v    0   )   k-s
  96        U(k) =  (   0    I    0   )   s
  97                (   0    0    I   )   n-k
  98                   k-s   s   n-k
  99
 100     If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
 101     If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
 102     and A(k,k), and v overwrites A(1:k-2,k-1:k).
 103
 104     If UPLO = 'L', then A = L*D*L', where
 105        L = P(1)*L(1)* ... *P(k)*L(k)* ...,
 106     i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
 107     n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
 108     and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
 109     defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
 110     that if the diagonal block D(k) is of order s (s = 1 or 2), then
 111
 112                (   I    0     0   )  k-1
 113        L(k) =  (   0    I     0   )  s
 114                (   0    v     I   )  n-k-s+1
 115                   k-1   s  n-k-s+1
 116
 117     If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
 118     If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
 119     and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
 120
 121     =====================================================================
 122
 123
 124        Test the input parameters.
 125
 126        Parameter adjustments */
 127     /* Table of constant values */
 128     static integer c__1 = 1;
 129
 130     /* System generated locals */
 131     integer a_dim1, a_offset, i__1, i__2;
 132     doublereal d__1, d__2, d__3;
 133     /* Builtin functions */
 134     double sqrt(doublereal);
 135     /* Local variables */
 136     static integer i__, j, k;
 137     static doublereal t, r1, d11, d12, d21, d22;
 138     static integer kk, kp;
 139     static doublereal wk, wkm1, wkp1;
 140     static integer imax, jmax;
 141     extern /* Subroutine */ int dsyr_(char *, integer *, doublereal *,
 142             doublereal *, integer *, doublereal *, integer *);
 143     static doublereal alpha;
 144     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
 145             integer *);
 146     extern logical lsame_(char *, char *);
 147     extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *,
 148             doublereal *, integer *);
 149     static integer kstep;
 150     static logical upper;
 151     static doublereal absakk;
 152     extern integer idamax_(integer *, doublereal *, integer *);
 153     extern logical disnan_(doublereal *);
 154     extern /* Subroutine */ int xerbla_(char *, integer *);
 155     static doublereal colmax, rowmax;
 156
 157
 158     a_dim1 = *lda;
 159     a_offset = 1 + a_dim1;
 160     a -= a_offset;
 161     --ipiv;
 162
 163     /* Function Body */
 164     *info = 0;
 165     upper = lsame_(uplo, "U");
 166     if (! upper && ! lsame_(uplo, "L")) {
 167         *info = -1;
 168     } else if (*n < 0) {
 169         *info = -2;
 170     } else if (*lda < max(1,*n)) {
 171         *info = -4;
 172     }
 173     if (*info != 0) {
 174         i__1 = -(*info);
 175         xerbla_("DSYTF2", &i__1);
 176         return 0;
 177     }
 178
 179 /*     Initialize ALPHA for use in choosing pivot block size. */
 180
 181     alpha = (sqrt(17.) + 1.) / 8.;
 182
 183     if (upper) {
 184
 185 /*        Factorize A as U*D*U' using the upper triangle of A
 186
 187           K is the main loop index, decreasing from N to 1 in steps of
 188           1 or 2 */
 189
 190         k = *n;
 191 L10:
 192
 193 /*        If K < 1, exit from loop */
 194
 195         if (k < 1) {
 196             goto L70;
 197         }
 198         kstep = 1;
 199
 200 /*        Determine rows and columns to be interchanged and whether
 201           a 1-by-1 or 2-by-2 pivot block will be used */
 202
 203         absakk = (d__1 = a[k + k * a_dim1], abs(d__1));
 204
 205 /*        IMAX is the row-index of the largest off-diagonal element in
 206           column K, and COLMAX is its absolute value */
 207
 208         if (k > 1) {
 209             i__1 = k - 1;
 210             imax = idamax_(&i__1, &a[k * a_dim1 + 1], &c__1);
 211             colmax = (d__1 = a[imax + k * a_dim1], abs(d__1));
 212         } else {
 213             colmax = 0.;
 214         }
 215
 216         if (max(absakk,colmax) == 0. || disnan_(&absakk)) {
 217
 218 /*           Column K is zero or contains a NaN: set INFO and continue */
 219
 220             if (*info == 0) {
 221                 *info = k;
 222             }
 223             kp = k;
 224         } else {
 225             if (absakk >= alpha * colmax) {
 226
 227 /*              no interchange, use 1-by-1 pivot block */
 228
 229                 kp = k;
 230             } else {
 231
 232 /*              JMAX is the column-index of the largest off-diagonal
 233                 element in row IMAX, and ROWMAX is its absolute value */
 234
 235                 i__1 = k - imax;
 236                 jmax = imax + idamax_(&i__1, &a[imax + (imax + 1) * a_dim1],
 237                         lda);
 238                 rowmax = (d__1 = a[imax + jmax * a_dim1], abs(d__1));
 239                 if (imax > 1) {
 240                     i__1 = imax - 1;
 241                     jmax = idamax_(&i__1, &a[imax * a_dim1 + 1], &c__1);
 242 /* Computing MAX */
 243                     d__2 = rowmax, d__3 = (d__1 = a[jmax + imax * a_dim1],
 244                             abs(d__1));
 245                     rowmax = max(d__2,d__3);
 246                 }
 247
 248                 if (absakk >= alpha * colmax * (colmax / rowmax)) {
 249
 250 /*                 no interchange, use 1-by-1 pivot block */
 251
 252                     kp = k;
 253                 } else if ((d__1 = a[imax + imax * a_dim1], abs(d__1)) >=
 254                         alpha * rowmax) {
 255
 256 /*                 interchange rows and columns K and IMAX, use 1-by-1
 257                    pivot block */
 258
 259                     kp = imax;
 260                 } else {
 261
 262 /*                 interchange rows and columns K-1 and IMAX, use 2-by-2
 263                    pivot block */
 264
 265                     kp = imax;
 266                     kstep = 2;
 267                 }
 268             }
 269
 270             kk = k - kstep + 1;
 271             if (kp != kk) {
 272
 273 /*              Interchange rows and columns KK and KP in the leading
 274                 submatrix A(1:k,1:k) */
 275
 276                 i__1 = kp - 1;
 277                 dswap_(&i__1, &a[kk * a_dim1 + 1], &c__1, &a[kp * a_dim1 + 1],
 278                          &c__1);
 279                 i__1 = kk - kp - 1;
 280                 dswap_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + (kp +
 281                         1) * a_dim1], lda);
 282                 t = a[kk + kk * a_dim1];
 283                 a[kk + kk * a_dim1] = a[kp + kp * a_dim1];
 284                 a[kp + kp * a_dim1] = t;
 285                 if (kstep == 2) {
 286                     t = a[k - 1 + k * a_dim1];
 287                     a[k - 1 + k * a_dim1] = a[kp + k * a_dim1];
 288                     a[kp + k * a_dim1] = t;
 289                 }
 290             }
 291
 292 /*           Update the leading submatrix */
 293
 294             if (kstep == 1) {
 295
 296 /*              1-by-1 pivot block D(k): column k now holds
 297
 298                 W(k) = U(k)*D(k)
 299
 300                 where U(k) is the k-th column of U
 301
 302                 Perform a rank-1 update of A(1:k-1,1:k-1) as
 303
 304                 A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)' */
 305
 306                 r1 = 1. / a[k + k * a_dim1];
 307                 i__1 = k - 1;
 308                 d__1 = -r1;
 309                 dsyr_(uplo, &i__1, &d__1, &a[k * a_dim1 + 1], &c__1, &a[
 310                         a_offset], lda);
 311
 312 /*              Store U(k) in column k */
 313
 314                 i__1 = k - 1;
 315                 dscal_(&i__1, &r1, &a[k * a_dim1 + 1], &c__1);
 316             } else {
 317
 318 /*              2-by-2 pivot block D(k): columns k and k-1 now hold
 319
 320                 ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 321
 322                 where U(k) and U(k-1) are the k-th and (k-1)-th columns
 323                 of U
 324
 325                 Perform a rank-2 update of A(1:k-2,1:k-2) as
 326
 327                 A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
 328                    = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )' */
 329
 330                 if (k > 2) {
 331
 332                     d12 = a[k - 1 + k * a_dim1];
 333                     d22 = a[k - 1 + (k - 1) * a_dim1] / d12;
 334                     d11 = a[k + k * a_dim1] / d12;
 335                     t = 1. / (d11 * d22 - 1.);
 336                     d12 = t / d12;
 337
 338                     for (j = k - 2; j >= 1; --j) {
 339                         wkm1 = d12 * (d11 * a[j + (k - 1) * a_dim1] - a[j + k
 340                                 * a_dim1]);
 341                         wk = d12 * (d22 * a[j + k * a_dim1] - a[j + (k - 1) *
 342                                 a_dim1]);
 343                         for (i__ = j; i__ >= 1; --i__) {
 344                             a[i__ + j * a_dim1] = a[i__ + j * a_dim1] - a[i__
 345                                     + k * a_dim1] * wk - a[i__ + (k - 1) *
 346                                     a_dim1] * wkm1;
 347 /* L20: */
 348                         }
 349                         a[j + k * a_dim1] = wk;
 350                         a[j + (k - 1) * a_dim1] = wkm1;
 351 /* L30: */
 352                     }
 353
 354                 }
 355
 356             }
 357         }
 358
 359 /*        Store details of the interchanges in IPIV */
 360
 361         if (kstep == 1) {
 362             ipiv[k] = kp;
 363         } else {
 364             ipiv[k] = -kp;
 365             ipiv[k - 1] = -kp;
 366         }
 367
 368 /*        Decrease K and return to the start of the main loop */
 369
 370         k -= kstep;
 371         goto L10;
 372
 373     } else {
 374
 375 /*        Factorize A as L*D*L' using the lower triangle of A
 376
 377           K is the main loop index, increasing from 1 to N in steps of
 378           1 or 2 */
 379
 380         k = 1;
 381 L40:
 382
 383 /*        If K > N, exit from loop */
 384
 385         if (k > *n) {
 386             goto L70;
 387         }
 388         kstep = 1;
 389
 390 /*        Determine rows and columns to be interchanged and whether
 391           a 1-by-1 or 2-by-2 pivot block will be used */
 392
 393         absakk = (d__1 = a[k + k * a_dim1], abs(d__1));
 394
 395 /*        IMAX is the row-index of the largest off-diagonal element in
 396           column K, and COLMAX is its absolute value */
 397
 398         if (k < *n) {
 399             i__1 = *n - k;
 400             imax = k + idamax_(&i__1, &a[k + 1 + k * a_dim1], &c__1);
 401             colmax = (d__1 = a[imax + k * a_dim1], abs(d__1));
 402         } else {
 403             colmax = 0.;
 404         }
 405
 406         if (max(absakk,colmax) == 0. || disnan_(&absakk)) {
 407
 408 /*           Column K is zero or contains a NaN: set INFO and continue */
 409
 410             if (*info == 0) {
 411                 *info = k;
 412             }
 413             kp = k;
 414         } else {
 415             if (absakk >= alpha * colmax) {
 416
 417 /*              no interchange, use 1-by-1 pivot block */
 418
 419                 kp = k;
 420             } else {
 421
 422 /*              JMAX is the column-index of the largest off-diagonal
 423                 element in row IMAX, and ROWMAX is its absolute value */
 424
 425                 i__1 = imax - k;
 426                 jmax = k - 1 + idamax_(&i__1, &a[imax + k * a_dim1], lda);
 427                 rowmax = (d__1 = a[imax + jmax * a_dim1], abs(d__1));
 428                 if (imax < *n) {
 429                     i__1 = *n - imax;
 430                     jmax = imax + idamax_(&i__1, &a[imax + 1 + imax * a_dim1],
 431                              &c__1);
 432 /* Computing MAX */
 433                     d__2 = rowmax, d__3 = (d__1 = a[jmax + imax * a_dim1],
 434                             abs(d__1));
 435                     rowmax = max(d__2,d__3);
 436                 }
 437
 438                 if (absakk >= alpha * colmax * (colmax / rowmax)) {
 439
 440 /*                 no interchange, use 1-by-1 pivot block */
 441
 442                     kp = k;
 443                 } else if ((d__1 = a[imax + imax * a_dim1], abs(d__1)) >=
 444                         alpha * rowmax) {
 445
 446 /*                 interchange rows and columns K and IMAX, use 1-by-1
 447                    pivot block */
 448
 449                     kp = imax;
 450                 } else {
 451
 452 /*                 interchange rows and columns K+1 and IMAX, use 2-by-2
 453                    pivot block */
 454
 455                     kp = imax;
 456                     kstep = 2;
 457                 }
 458             }
 459
 460             kk = k + kstep - 1;
 461             if (kp != kk) {
 462
 463 /*              Interchange rows and columns KK and KP in the trailing
 464                 submatrix A(k:n,k:n) */
 465
 466                 if (kp < *n) {
 467                     i__1 = *n - kp;
 468                     dswap_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + 1
 469                             + kp * a_dim1], &c__1);
 470                 }
 471                 i__1 = kp - kk - 1;
 472                 dswap_(&i__1, &a[kk + 1 + kk * a_dim1], &c__1, &a[kp + (kk +
 473                         1) * a_dim1], lda);
 474                 t = a[kk + kk * a_dim1];
 475                 a[kk + kk * a_dim1] = a[kp + kp * a_dim1];
 476                 a[kp + kp * a_dim1] = t;
 477                 if (kstep == 2) {
 478                     t = a[k + 1 + k * a_dim1];
 479                     a[k + 1 + k * a_dim1] = a[kp + k * a_dim1];
 480                     a[kp + k * a_dim1] = t;
 481                 }
 482             }
 483
 484 /*           Update the trailing submatrix */
 485
 486             if (kstep == 1) {
 487
 488 /*              1-by-1 pivot block D(k): column k now holds
 489
 490                 W(k) = L(k)*D(k)
 491
 492                 where L(k) is the k-th column of L */
 493
 494                 if (k < *n) {
 495
 496 /*                 Perform a rank-1 update of A(k+1:n,k+1:n) as
 497
 498                    A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)' */
 499
 500                     d11 = 1. / a[k + k * a_dim1];
 501                     i__1 = *n - k;
 502                     d__1 = -d11;
 503                     dsyr_(uplo, &i__1, &d__1, &a[k + 1 + k * a_dim1], &c__1, &
 504                             a[k + 1 + (k + 1) * a_dim1], lda);
 505
 506 /*                 Store L(k) in column K */
 507
 508                     i__1 = *n - k;
 509                     dscal_(&i__1, &d11, &a[k + 1 + k * a_dim1], &c__1);
 510                 }
 511             } else {
 512
 513 /*              2-by-2 pivot block D(k) */
 514
 515                 if (k < *n - 1) {
 516
 517 /*                 Perform a rank-2 update of A(k+2:n,k+2:n) as
 518
 519                    A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))'
 520
 521                    where L(k) and L(k+1) are the k-th and (k+1)-th
 522                    columns of L */
 523
 524                     d21 = a[k + 1 + k * a_dim1];
 525                     d11 = a[k + 1 + (k + 1) * a_dim1] / d21;
 526                     d22 = a[k + k * a_dim1] / d21;
 527                     t = 1. / (d11 * d22 - 1.);
 528                     d21 = t / d21;
 529
 530                     i__1 = *n;
 531                     for (j = k + 2; j <= i__1; ++j) {
 532
 533                         wk = d21 * (d11 * a[j + k * a_dim1] - a[j + (k + 1) *
 534                                 a_dim1]);
 535                         wkp1 = d21 * (d22 * a[j + (k + 1) * a_dim1] - a[j + k
 536                                 * a_dim1]);
 537
 538                         i__2 = *n;
 539                         for (i__ = j; i__ <= i__2; ++i__) {
 540                             a[i__ + j * a_dim1] = a[i__ + j * a_dim1] - a[i__
 541                                     + k * a_dim1] * wk - a[i__ + (k + 1) *
 542                                     a_dim1] * wkp1;
 543 /* L50: */
 544                         }
 545
 546                         a[j + k * a_dim1] = wk;
 547                         a[j + (k + 1) * a_dim1] = wkp1;
 548
 549 /* L60: */
 550                     }
 551                 }
 552             }
 553         }
 554
 555 /*        Store details of the interchanges in IPIV */
 556
 557         if (kstep == 1) {
 558             ipiv[k] = kp;
 559         } else {
 560             ipiv[k] = -kp;
 561             ipiv[k + 1] = -kp;
 562         }
 563
 564 /*        Increase K and return to the start of the main loop */
 565
 566         k += kstep;
 567         goto L40;
 568
 569     }
 570
 571 L70:
 572
 573     return 0;
 574
 575 /*     End of DSYTF2 */
 576
 577 } /* dsytf2_ */