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

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlasyf_(char *uplo, integer *n, integer *nb, integer *kb,
   4          doublereal *a, integer *lda, integer *ipiv, doublereal *w, integer *
   5         ldw, integer *info)
   6 {
   7 /*  -- LAPACK routine (version 3.1) --
   8        Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   9        November 2006
  10
  11
  12     Purpose
  13     =======
  14
  15     DLASYF computes a partial factorization of a real symmetric matrix A
  16     using the Bunch-Kaufman diagonal pivoting method. The partial
  17     factorization has the form:
  18
  19     A  =  ( I  U12 ) ( A11  0  ) (  I    0   )  if UPLO = 'U', or:
  20           ( 0  U22 ) (  0   D  ) ( U12' U22' )
  21
  22     A  =  ( L11  0 ) (  D   0  ) ( L11' L21' )  if UPLO = 'L'
  23           ( L21  I ) (  0  A22 ) (  0    I   )
  24
  25     where the order of D is at most NB. The actual order is returned in
  26     the argument KB, and is either NB or NB-1, or N if N <= NB.
  27
  28     DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
  29     (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
  30     A22 (if UPLO = 'L').
  31
  32     Arguments
  33     =========
  34
  35     UPLO    (input) CHARACTER*1
  36             Specifies whether the upper or lower triangular part of the
  37             symmetric matrix A is stored:
  38             = 'U':  Upper triangular
  39             = 'L':  Lower triangular
  40
  41     N       (input) INTEGER
  42             The order of the matrix A.  N >= 0.
  43
  44     NB      (input) INTEGER
  45             The maximum number of columns of the matrix A that should be
  46             factored.  NB should be at least 2 to allow for 2-by-2 pivot
  47             blocks.
  48
  49     KB      (output) INTEGER
  50             The number of columns of A that were actually factored.
  51             KB is either NB-1 or NB, or N if N <= NB.
  52
  53     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  54             On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  55             n-by-n upper triangular part of A contains the upper
  56             triangular part of the matrix A, and the strictly lower
  57             triangular part of A is not referenced.  If UPLO = 'L', the
  58             leading n-by-n lower triangular part of A contains the lower
  59             triangular part of the matrix A, and the strictly upper
  60             triangular part of A is not referenced.
  61             On exit, A contains details of the partial factorization.
  62
  63     LDA     (input) INTEGER
  64             The leading dimension of the array A.  LDA >= max(1,N).
  65
  66     IPIV    (output) INTEGER array, dimension (N)
  67             Details of the interchanges and the block structure of D.
  68             If UPLO = 'U', only the last KB elements of IPIV are set;
  69             if UPLO = 'L', only the first KB elements are set.
  70
  71             If IPIV(k) > 0, then rows and columns k and IPIV(k) were
  72             interchanged and D(k,k) is a 1-by-1 diagonal block.
  73             If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
  74             columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
  75             is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
  76             IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
  77             interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
  78
  79     W       (workspace) DOUBLE PRECISION array, dimension (LDW,NB)
  80
  81     LDW     (input) INTEGER
  82             The leading dimension of the array W.  LDW >= max(1,N).
  83
  84     INFO    (output) INTEGER
  85             = 0: successful exit
  86             > 0: if INFO = k, D(k,k) is exactly zero.  The factorization
  87                  has been completed, but the block diagonal matrix D is
  88                  exactly singular.
  89
  90     =====================================================================
  91
  92
  93        Parameter adjustments */
  94     /* Table of constant values */
  95     static integer c__1 = 1;
  96     static doublereal c_b8 = -1.;
  97     static doublereal c_b9 = 1.;
  98
  99     /* System generated locals */
 100     integer a_dim1, a_offset, w_dim1, w_offset, i__1, i__2, i__3, i__4, i__5;
 101     doublereal d__1, d__2, d__3;
 102     /* Builtin functions */
 103     double sqrt(doublereal);
 104     /* Local variables */
 105     static integer j, k;
 106     static doublereal t, r1, d11, d21, d22;
 107     static integer jb, jj, kk, jp, kp, kw, kkw, imax, jmax;
 108     static doublereal alpha;
 109     extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
 110             integer *), dgemm_(char *, char *, integer *, integer *, integer *
 111 , doublereal *, doublereal *, integer *, doublereal *, integer *,
 112             doublereal *, doublereal *, integer *);
 113     extern logical lsame_(char *, char *);
 114     extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
 115             doublereal *, doublereal *, integer *, doublereal *, integer *,
 116             doublereal *, doublereal *, integer *), dcopy_(integer *,
 117             doublereal *, integer *, doublereal *, integer *), dswap_(integer
 118             *, doublereal *, integer *, doublereal *, integer *);
 119     static integer kstep;
 120     static doublereal absakk;
 121     extern integer idamax_(integer *, doublereal *, integer *);
 122     static doublereal colmax, rowmax;
 123
 124
 125     a_dim1 = *lda;
 126     a_offset = 1 + a_dim1;
 127     a -= a_offset;
 128     --ipiv;
 129     w_dim1 = *ldw;
 130     w_offset = 1 + w_dim1;
 131     w -= w_offset;
 132
 133     /* Function Body */
 134     *info = 0;
 135
 136 /*     Initialize ALPHA for use in choosing pivot block size. */
 137
 138     alpha = (sqrt(17.) + 1.) / 8.;
 139
 140     if (lsame_(uplo, "U")) {
 141
 142 /*        Factorize the trailing columns of A using the upper triangle
 143           of A and working backwards, and compute the matrix W = U12*D
 144           for use in updating A11
 145
 146           K is the main loop index, decreasing from N in steps of 1 or 2
 147
 148           KW is the column of W which corresponds to column K of A */
 149
 150         k = *n;
 151 L10:
 152         kw = *nb + k - *n;
 153
 154 /*        Exit from loop */
 155
 156         if (k <= *n - *nb + 1 && *nb < *n || k < 1) {
 157             goto L30;
 158         }
 159
 160 /*        Copy column K of A to column KW of W and update it */
 161
 162         dcopy_(&k, &a[k * a_dim1 + 1], &c__1, &w[kw * w_dim1 + 1], &c__1);
 163         if (k < *n) {
 164             i__1 = *n - k;
 165             dgemv_("No transpose", &k, &i__1, &c_b8, &a[(k + 1) * a_dim1 + 1],
 166                      lda, &w[k + (kw + 1) * w_dim1], ldw, &c_b9, &w[kw *
 167                     w_dim1 + 1], &c__1);
 168         }
 169
 170         kstep = 1;
 171
 172 /*        Determine rows and columns to be interchanged and whether
 173           a 1-by-1 or 2-by-2 pivot block will be used */
 174
 175         absakk = (d__1 = w[k + kw * w_dim1], abs(d__1));
 176
 177 /*        IMAX is the row-index of the largest off-diagonal element in
 178           column K, and COLMAX is its absolute value */
 179
 180         if (k > 1) {
 181             i__1 = k - 1;
 182             imax = idamax_(&i__1, &w[kw * w_dim1 + 1], &c__1);
 183             colmax = (d__1 = w[imax + kw * w_dim1], abs(d__1));
 184         } else {
 185             colmax = 0.;
 186         }
 187
 188         if (max(absakk,colmax) == 0.) {
 189
 190 /*           Column K is zero: set INFO and continue */
 191
 192             if (*info == 0) {
 193                 *info = k;
 194             }
 195             kp = k;
 196         } else {
 197             if (absakk >= alpha * colmax) {
 198
 199 /*              no interchange, use 1-by-1 pivot block */
 200
 201                 kp = k;
 202             } else {
 203
 204 /*              Copy column IMAX to column KW-1 of W and update it */
 205
 206                 dcopy_(&imax, &a[imax * a_dim1 + 1], &c__1, &w[(kw - 1) *
 207                         w_dim1 + 1], &c__1);
 208                 i__1 = k - imax;
 209                 dcopy_(&i__1, &a[imax + (imax + 1) * a_dim1], lda, &w[imax +
 210                         1 + (kw - 1) * w_dim1], &c__1);
 211                 if (k < *n) {
 212                     i__1 = *n - k;
 213                     dgemv_("No transpose", &k, &i__1, &c_b8, &a[(k + 1) *
 214                             a_dim1 + 1], lda, &w[imax + (kw + 1) * w_dim1],
 215                             ldw, &c_b9, &w[(kw - 1) * w_dim1 + 1], &c__1);
 216                 }
 217
 218 /*              JMAX is the column-index of the largest off-diagonal
 219                 element in row IMAX, and ROWMAX is its absolute value */
 220
 221                 i__1 = k - imax;
 222                 jmax = imax + idamax_(&i__1, &w[imax + 1 + (kw - 1) * w_dim1],
 223                          &c__1);
 224                 rowmax = (d__1 = w[jmax + (kw - 1) * w_dim1], abs(d__1));
 225                 if (imax > 1) {
 226                     i__1 = imax - 1;
 227                     jmax = idamax_(&i__1, &w[(kw - 1) * w_dim1 + 1], &c__1);
 228 /* Computing MAX */
 229                     d__2 = rowmax, d__3 = (d__1 = w[jmax + (kw - 1) * w_dim1],
 230                              abs(d__1));
 231                     rowmax = max(d__2,d__3);
 232                 }
 233
 234                 if (absakk >= alpha * colmax * (colmax / rowmax)) {
 235
 236 /*                 no interchange, use 1-by-1 pivot block */
 237
 238                     kp = k;
 239                 } else if ((d__1 = w[imax + (kw - 1) * w_dim1], abs(d__1)) >=
 240                         alpha * rowmax) {
 241
 242 /*                 interchange rows and columns K and IMAX, use 1-by-1
 243                    pivot block */
 244
 245                     kp = imax;
 246
 247 /*                 copy column KW-1 of W to column KW */
 248
 249                     dcopy_(&k, &w[(kw - 1) * w_dim1 + 1], &c__1, &w[kw *
 250                             w_dim1 + 1], &c__1);
 251                 } else {
 252
 253 /*                 interchange rows and columns K-1 and IMAX, use 2-by-2
 254                    pivot block */
 255
 256                     kp = imax;
 257                     kstep = 2;
 258                 }
 259             }
 260
 261             kk = k - kstep + 1;
 262             kkw = *nb + kk - *n;
 263
 264 /*           Updated column KP is already stored in column KKW of W */
 265
 266             if (kp != kk) {
 267
 268 /*              Copy non-updated column KK to column KP */
 269
 270                 a[kp + k * a_dim1] = a[kk + k * a_dim1];
 271                 i__1 = k - 1 - kp;
 272                 dcopy_(&i__1, &a[kp + 1 + kk * a_dim1], &c__1, &a[kp + (kp +
 273                         1) * a_dim1], lda);
 274                 dcopy_(&kp, &a[kk * a_dim1 + 1], &c__1, &a[kp * a_dim1 + 1], &
 275                         c__1);
 276
 277 /*              Interchange rows KK and KP in last KK columns of A and W */
 278
 279                 i__1 = *n - kk + 1;
 280                 dswap_(&i__1, &a[kk + kk * a_dim1], lda, &a[kp + kk * a_dim1],
 281                          lda);
 282                 i__1 = *n - kk + 1;
 283                 dswap_(&i__1, &w[kk + kkw * w_dim1], ldw, &w[kp + kkw *
 284                         w_dim1], ldw);
 285             }
 286
 287             if (kstep == 1) {
 288
 289 /*              1-by-1 pivot block D(k): column KW of W now holds
 290
 291                 W(k) = U(k)*D(k)
 292
 293                 where U(k) is the k-th column of U
 294
 295                 Store U(k) in column k of A */
 296
 297                 dcopy_(&k, &w[kw * w_dim1 + 1], &c__1, &a[k * a_dim1 + 1], &
 298                         c__1);
 299                 r1 = 1. / a[k + k * a_dim1];
 300                 i__1 = k - 1;
 301                 dscal_(&i__1, &r1, &a[k * a_dim1 + 1], &c__1);
 302             } else {
 303
 304 /*              2-by-2 pivot block D(k): columns KW and KW-1 of W now
 305                 hold
 306
 307                 ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
 308
 309                 where U(k) and U(k-1) are the k-th and (k-1)-th columns
 310                 of U */
 311
 312                 if (k > 2) {
 313
 314 /*                 Store U(k) and U(k-1) in columns k and k-1 of A */
 315
 316                     d21 = w[k - 1 + kw * w_dim1];
 317                     d11 = w[k + kw * w_dim1] / d21;
 318                     d22 = w[k - 1 + (kw - 1) * w_dim1] / d21;
 319                     t = 1. / (d11 * d22 - 1.);
 320                     d21 = t / d21;
 321                     i__1 = k - 2;
 322                     for (j = 1; j <= i__1; ++j) {
 323                         a[j + (k - 1) * a_dim1] = d21 * (d11 * w[j + (kw - 1)
 324                                 * w_dim1] - w[j + kw * w_dim1]);
 325                         a[j + k * a_dim1] = d21 * (d22 * w[j + kw * w_dim1] -
 326                                 w[j + (kw - 1) * w_dim1]);
 327 /* L20: */
 328                     }
 329                 }
 330
 331 /*              Copy D(k) to A */
 332
 333                 a[k - 1 + (k - 1) * a_dim1] = w[k - 1 + (kw - 1) * w_dim1];
 334                 a[k - 1 + k * a_dim1] = w[k - 1 + kw * w_dim1];
 335                 a[k + k * a_dim1] = w[k + kw * w_dim1];
 336             }
 337         }
 338
 339 /*        Store details of the interchanges in IPIV */
 340
 341         if (kstep == 1) {
 342             ipiv[k] = kp;
 343         } else {
 344             ipiv[k] = -kp;
 345             ipiv[k - 1] = -kp;
 346         }
 347
 348 /*        Decrease K and return to the start of the main loop */
 349
 350         k -= kstep;
 351         goto L10;
 352
 353 L30:
 354
 355 /*        Update the upper triangle of A11 (= A(1:k,1:k)) as
 356
 357           A11 := A11 - U12*D*U12' = A11 - U12*W'
 358
 359           computing blocks of NB columns at a time */
 360
 361         i__1 = -(*nb);
 362         for (j = (k - 1) / *nb * *nb + 1; i__1 < 0 ? j >= 1 : j <= 1; j +=
 363                 i__1) {
 364 /* Computing MIN */
 365             i__2 = *nb, i__3 = k - j + 1;
 366             jb = min(i__2,i__3);
 367
 368 /*           Update the upper triangle of the diagonal block */
 369
 370             i__2 = j + jb - 1;
 371             for (jj = j; jj <= i__2; ++jj) {
 372                 i__3 = jj - j + 1;
 373                 i__4 = *n - k;
 374                 dgemv_("No transpose", &i__3, &i__4, &c_b8, &a[j + (k + 1) *
 375                         a_dim1], lda, &w[jj + (kw + 1) * w_dim1], ldw, &c_b9,
 376                         &a[j + jj * a_dim1], &c__1);
 377 /* L40: */
 378             }
 379
 380 /*           Update the rectangular superdiagonal block */
 381
 382             i__2 = j - 1;
 383             i__3 = *n - k;
 384             dgemm_("No transpose", "Transpose", &i__2, &jb, &i__3, &c_b8, &a[(
 385                     k + 1) * a_dim1 + 1], lda, &w[j + (kw + 1) * w_dim1], ldw,
 386                      &c_b9, &a[j * a_dim1 + 1], lda);
 387 /* L50: */
 388         }
 389
 390 /*        Put U12 in standard form by partially undoing the interchanges
 391           in columns k+1:n */
 392
 393         j = k + 1;
 394 L60:
 395         jj = j;
 396         jp = ipiv[j];
 397         if (jp < 0) {
 398             jp = -jp;
 399             ++j;
 400         }
 401         ++j;
 402         if (jp != jj && j <= *n) {
 403             i__1 = *n - j + 1;
 404             dswap_(&i__1, &a[jp + j * a_dim1], lda, &a[jj + j * a_dim1], lda);
 405         }
 406         if (j <= *n) {
 407             goto L60;
 408         }
 409
 410 /*        Set KB to the number of columns factorized */
 411
 412         *kb = *n - k;
 413
 414     } else {
 415
 416 /*        Factorize the leading columns of A using the lower triangle
 417           of A and working forwards, and compute the matrix W = L21*D
 418           for use in updating A22
 419
 420           K is the main loop index, increasing from 1 in steps of 1 or 2 */
 421
 422         k = 1;
 423 L70:
 424
 425 /*        Exit from loop */
 426
 427         if (k >= *nb && *nb < *n || k > *n) {
 428             goto L90;
 429         }
 430
 431 /*        Copy column K of A to column K of W and update it */
 432
 433         i__1 = *n - k + 1;
 434         dcopy_(&i__1, &a[k + k * a_dim1], &c__1, &w[k + k * w_dim1], &c__1);
 435         i__1 = *n - k + 1;
 436         i__2 = k - 1;
 437         dgemv_("No transpose", &i__1, &i__2, &c_b8, &a[k + a_dim1], lda, &w[k
 438                 + w_dim1], ldw, &c_b9, &w[k + k * w_dim1], &c__1);
 439
 440         kstep = 1;
 441
 442 /*        Determine rows and columns to be interchanged and whether
 443           a 1-by-1 or 2-by-2 pivot block will be used */
 444
 445         absakk = (d__1 = w[k + k * w_dim1], abs(d__1));
 446
 447 /*        IMAX is the row-index of the largest off-diagonal element in
 448           column K, and COLMAX is its absolute value */
 449
 450         if (k < *n) {
 451             i__1 = *n - k;
 452             imax = k + idamax_(&i__1, &w[k + 1 + k * w_dim1], &c__1);
 453             colmax = (d__1 = w[imax + k * w_dim1], abs(d__1));
 454         } else {
 455             colmax = 0.;
 456         }
 457
 458         if (max(absakk,colmax) == 0.) {
 459
 460 /*           Column K is zero: set INFO and continue */
 461
 462             if (*info == 0) {
 463                 *info = k;
 464             }
 465             kp = k;
 466         } else {
 467             if (absakk >= alpha * colmax) {
 468
 469 /*              no interchange, use 1-by-1 pivot block */
 470
 471                 kp = k;
 472             } else {
 473
 474 /*              Copy column IMAX to column K+1 of W and update it */
 475
 476                 i__1 = imax - k;
 477                 dcopy_(&i__1, &a[imax + k * a_dim1], lda, &w[k + (k + 1) *
 478                         w_dim1], &c__1);
 479                 i__1 = *n - imax + 1;
 480                 dcopy_(&i__1, &a[imax + imax * a_dim1], &c__1, &w[imax + (k +
 481                         1) * w_dim1], &c__1);
 482                 i__1 = *n - k + 1;
 483                 i__2 = k - 1;
 484                 dgemv_("No transpose", &i__1, &i__2, &c_b8, &a[k + a_dim1],
 485                         lda, &w[imax + w_dim1], ldw, &c_b9, &w[k + (k + 1) *
 486                         w_dim1], &c__1);
 487
 488 /*              JMAX is the column-index of the largest off-diagonal
 489                 element in row IMAX, and ROWMAX is its absolute value */
 490
 491                 i__1 = imax - k;
 492                 jmax = k - 1 + idamax_(&i__1, &w[k + (k + 1) * w_dim1], &c__1)
 493                         ;
 494                 rowmax = (d__1 = w[jmax + (k + 1) * w_dim1], abs(d__1));
 495                 if (imax < *n) {
 496                     i__1 = *n - imax;
 497                     jmax = imax + idamax_(&i__1, &w[imax + 1 + (k + 1) *
 498                             w_dim1], &c__1);
 499 /* Computing MAX */
 500                     d__2 = rowmax, d__3 = (d__1 = w[jmax + (k + 1) * w_dim1],
 501                             abs(d__1));
 502                     rowmax = max(d__2,d__3);
 503                 }
 504
 505                 if (absakk >= alpha * colmax * (colmax / rowmax)) {
 506
 507 /*                 no interchange, use 1-by-1 pivot block */
 508
 509                     kp = k;
 510                 } else if ((d__1 = w[imax + (k + 1) * w_dim1], abs(d__1)) >=
 511                         alpha * rowmax) {
 512
 513 /*                 interchange rows and columns K and IMAX, use 1-by-1
 514                    pivot block */
 515
 516                     kp = imax;
 517
 518 /*                 copy column K+1 of W to column K */
 519
 520                     i__1 = *n - k + 1;
 521                     dcopy_(&i__1, &w[k + (k + 1) * w_dim1], &c__1, &w[k + k *
 522                             w_dim1], &c__1);
 523                 } else {
 524
 525 /*                 interchange rows and columns K+1 and IMAX, use 2-by-2
 526                    pivot block */
 527
 528                     kp = imax;
 529                     kstep = 2;
 530                 }
 531             }
 532
 533             kk = k + kstep - 1;
 534
 535 /*           Updated column KP is already stored in column KK of W */
 536
 537             if (kp != kk) {
 538
 539 /*              Copy non-updated column KK to column KP */
 540
 541                 a[kp + k * a_dim1] = a[kk + k * a_dim1];
 542                 i__1 = kp - k - 1;
 543                 dcopy_(&i__1, &a[k + 1 + kk * a_dim1], &c__1, &a[kp + (k + 1)
 544                         * a_dim1], lda);
 545                 i__1 = *n - kp + 1;
 546                 dcopy_(&i__1, &a[kp + kk * a_dim1], &c__1, &a[kp + kp *
 547                         a_dim1], &c__1);
 548
 549 /*              Interchange rows KK and KP in first KK columns of A and W */
 550
 551                 dswap_(&kk, &a[kk + a_dim1], lda, &a[kp + a_dim1], lda);
 552                 dswap_(&kk, &w[kk + w_dim1], ldw, &w[kp + w_dim1], ldw);
 553             }
 554
 555             if (kstep == 1) {
 556
 557 /*              1-by-1 pivot block D(k): column k of W now holds
 558
 559                 W(k) = L(k)*D(k)
 560
 561                 where L(k) is the k-th column of L
 562
 563                 Store L(k) in column k of A */
 564
 565                 i__1 = *n - k + 1;
 566                 dcopy_(&i__1, &w[k + k * w_dim1], &c__1, &a[k + k * a_dim1], &
 567                         c__1);
 568                 if (k < *n) {
 569                     r1 = 1. / a[k + k * a_dim1];
 570                     i__1 = *n - k;
 571                     dscal_(&i__1, &r1, &a[k + 1 + k * a_dim1], &c__1);
 572                 }
 573             } else {
 574
 575 /*              2-by-2 pivot block D(k): columns k and k+1 of W now hold
 576
 577                 ( W(k) W(k+1) ) = ( L(k) L(k+1) )*D(k)
 578
 579                 where L(k) and L(k+1) are the k-th and (k+1)-th columns
 580                 of L */
 581
 582                 if (k < *n - 1) {
 583
 584 /*                 Store L(k) and L(k+1) in columns k and k+1 of A */
 585
 586                     d21 = w[k + 1 + k * w_dim1];
 587                     d11 = w[k + 1 + (k + 1) * w_dim1] / d21;
 588                     d22 = w[k + k * w_dim1] / d21;
 589                     t = 1. / (d11 * d22 - 1.);
 590                     d21 = t / d21;
 591                     i__1 = *n;
 592                     for (j = k + 2; j <= i__1; ++j) {
 593                         a[j + k * a_dim1] = d21 * (d11 * w[j + k * w_dim1] -
 594                                 w[j + (k + 1) * w_dim1]);
 595                         a[j + (k + 1) * a_dim1] = d21 * (d22 * w[j + (k + 1) *
 596                                  w_dim1] - w[j + k * w_dim1]);
 597 /* L80: */
 598                     }
 599                 }
 600
 601 /*              Copy D(k) to A */
 602
 603                 a[k + k * a_dim1] = w[k + k * w_dim1];
 604                 a[k + 1 + k * a_dim1] = w[k + 1 + k * w_dim1];
 605                 a[k + 1 + (k + 1) * a_dim1] = w[k + 1 + (k + 1) * w_dim1];
 606             }
 607         }
 608
 609 /*        Store details of the interchanges in IPIV */
 610
 611         if (kstep == 1) {
 612             ipiv[k] = kp;
 613         } else {
 614             ipiv[k] = -kp;
 615             ipiv[k + 1] = -kp;
 616         }
 617
 618 /*        Increase K and return to the start of the main loop */
 619
 620         k += kstep;
 621         goto L70;
 622
 623 L90:
 624
 625 /*        Update the lower triangle of A22 (= A(k:n,k:n)) as
 626
 627           A22 := A22 - L21*D*L21' = A22 - L21*W'
 628
 629           computing blocks of NB columns at a time */
 630
 631         i__1 = *n;
 632         i__2 = *nb;
 633         for (j = k; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
 634 /* Computing MIN */
 635             i__3 = *nb, i__4 = *n - j + 1;
 636             jb = min(i__3,i__4);
 637
 638 /*           Update the lower triangle of the diagonal block */
 639
 640             i__3 = j + jb - 1;
 641             for (jj = j; jj <= i__3; ++jj) {
 642                 i__4 = j + jb - jj;
 643                 i__5 = k - 1;
 644                 dgemv_("No transpose", &i__4, &i__5, &c_b8, &a[jj + a_dim1],
 645                         lda, &w[jj + w_dim1], ldw, &c_b9, &a[jj + jj * a_dim1]
 646 , &c__1);
 647 /* L100: */
 648             }
 649
 650 /*           Update the rectangular subdiagonal block */
 651
 652             if (j + jb <= *n) {
 653                 i__3 = *n - j - jb + 1;
 654                 i__4 = k - 1;
 655                 dgemm_("No transpose", "Transpose", &i__3, &jb, &i__4, &c_b8,
 656                         &a[j + jb + a_dim1], lda, &w[j + w_dim1], ldw, &c_b9,
 657                         &a[j + jb + j * a_dim1], lda);
 658             }
 659 /* L110: */
 660         }
 661
 662 /*        Put L21 in standard form by partially undoing the interchanges
 663           in columns 1:k-1 */
 664
 665         j = k - 1;
 666 L120:
 667         jj = j;
 668         jp = ipiv[j];
 669         if (jp < 0) {
 670             jp = -jp;
 671             --j;
 672         }
 673         --j;
 674         if (jp != jj && j >= 1) {
 675             dswap_(&j, &a[jp + a_dim1], lda, &a[jj + a_dim1], lda);
 676         }
 677         if (j >= 1) {
 678             goto L120;
 679         }
 680
 681 /*        Set KB to the number of columns factorized */
 682
 683         *kb = k - 1;
 684
 685     }
 686     return 0;
 687
 688 /*     End of DLASYF */
 689
 690 } /* dlasyf_ */