Annotation of rpl/lapack/lapack/dorbdb.f, revision 1.2

1.1       bertrand    1:       SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
                      2:      $                   X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
                      3:      $                   TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
                      4:       IMPLICIT NONE
                      5: *
                      6: *  -- LAPACK routine (version 3.3.0) --
                      7: *
                      8: *  -- Contributed by Brian Sutton of the Randolph-Macon College --
                      9: *  -- November 2010
                     10: *
                     11: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                     12: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
                     13: *
                     14: *     .. Scalar Arguments ..
                     15:       CHARACTER          SIGNS, TRANS
                     16:       INTEGER            INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
                     17:      $                   Q
                     18: *     ..
                     19: *     .. Array Arguments ..
                     20:       DOUBLE PRECISION   PHI( * ), THETA( * )
                     21:       DOUBLE PRECISION   TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
                     22:      $                   WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
                     23:      $                   X21( LDX21, * ), X22( LDX22, * )
                     24: *     ..
                     25: *
                     26: *  Purpose
                     27: *  =======
                     28: *
                     29: *  DORBDB simultaneously bidiagonalizes the blocks of an M-by-M
                     30: *  partitioned orthogonal matrix X:
                     31: *
                     32: *                                  [ B11 | B12 0  0 ]
                     33: *      [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**T
                     34: *  X = [-----------] = [---------] [----------------] [---------]   .
                     35: *      [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
                     36: *                                  [  0  |  0  0  I ]
                     37: *
                     38: *  X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
                     39: *  not the case, then X must be transposed and/or permuted. This can be
                     40: *  done in constant time using the TRANS and SIGNS options. See DORCSD
                     41: *  for details.)
                     42: *
                     43: *  The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
                     44: *  (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
                     45: *  represented implicitly by Householder vectors.
                     46: *
                     47: *  B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
                     48: *  implicitly by angles THETA, PHI.
                     49: *
                     50: *  Arguments
                     51: *  =========
                     52: *
                     53: *  TRANS   (input) CHARACTER
                     54: *          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                     55: *                      order;
                     56: *          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                     57: *                      major order.
                     58: *
                     59: *  SIGNS   (input) CHARACTER
                     60: *          = 'O':      The lower-left block is made nonpositive (the
                     61: *                      "other" convention);
                     62: *          otherwise:  The upper-right block is made nonpositive (the
                     63: *                      "default" convention).
                     64: *
                     65: *  M       (input) INTEGER
                     66: *          The number of rows and columns in X.
                     67: *
                     68: *  P       (input) INTEGER
                     69: *          The number of rows in X11 and X12. 0 <= P <= M.
                     70: *
                     71: *  Q       (input) INTEGER
                     72: *          The number of columns in X11 and X21. 0 <= Q <=
                     73: *          MIN(P,M-P,M-Q).
                     74: *
                     75: *  X11     (input/output) DOUBLE PRECISION array, dimension (LDX11,Q)
                     76: *          On entry, the top-left block of the orthogonal matrix to be
                     77: *          reduced. On exit, the form depends on TRANS:
                     78: *          If TRANS = 'N', then
                     79: *             the columns of tril(X11) specify reflectors for P1,
                     80: *             the rows of triu(X11,1) specify reflectors for Q1;
                     81: *          else TRANS = 'T', and
                     82: *             the rows of triu(X11) specify reflectors for P1,
                     83: *             the columns of tril(X11,-1) specify reflectors for Q1.
                     84: *
                     85: *  LDX11   (input) INTEGER
                     86: *          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
                     87: *          P; else LDX11 >= Q.
                     88: *
                     89: *  X12     (input/output) DOUBLE PRECISION array, dimension (LDX12,M-Q)
                     90: *          On entry, the top-right block of the orthogonal matrix to
                     91: *          be reduced. On exit, the form depends on TRANS:
                     92: *          If TRANS = 'N', then
                     93: *             the rows of triu(X12) specify the first P reflectors for
                     94: *             Q2;
                     95: *          else TRANS = 'T', and
                     96: *             the columns of tril(X12) specify the first P reflectors
                     97: *             for Q2.
                     98: *
                     99: *  LDX12   (input) INTEGER
                    100: *          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
                    101: *          P; else LDX11 >= M-Q.
                    102: *
                    103: *  X21     (input/output) DOUBLE PRECISION array, dimension (LDX21,Q)
                    104: *          On entry, the bottom-left block of the orthogonal matrix to
                    105: *          be reduced. On exit, the form depends on TRANS:
                    106: *          If TRANS = 'N', then
                    107: *             the columns of tril(X21) specify reflectors for P2;
                    108: *          else TRANS = 'T', and
                    109: *             the rows of triu(X21) specify reflectors for P2.
                    110: *
                    111: *  LDX21   (input) INTEGER
                    112: *          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
                    113: *          M-P; else LDX21 >= Q.
                    114: *
                    115: *  X22     (input/output) DOUBLE PRECISION array, dimension (LDX22,M-Q)
                    116: *          On entry, the bottom-right block of the orthogonal matrix to
                    117: *          be reduced. On exit, the form depends on TRANS:
                    118: *          If TRANS = 'N', then
                    119: *             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
                    120: *             M-P-Q reflectors for Q2,
                    121: *          else TRANS = 'T', and
                    122: *             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
                    123: *             M-P-Q reflectors for P2.
                    124: *
                    125: *  LDX22   (input) INTEGER
                    126: *          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
                    127: *          M-P; else LDX22 >= M-Q.
                    128: *
                    129: *  THETA   (output) DOUBLE PRECISION array, dimension (Q)
                    130: *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
                    131: *          be computed from the angles THETA and PHI. See Further
                    132: *          Details.
                    133: *
                    134: *  PHI     (output) DOUBLE PRECISION array, dimension (Q-1)
                    135: *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
                    136: *          be computed from the angles THETA and PHI. See Further
                    137: *          Details.
                    138: *
                    139: *  TAUP1   (output) DOUBLE PRECISION array, dimension (P)
                    140: *          The scalar factors of the elementary reflectors that define
                    141: *          P1.
                    142: *
                    143: *  TAUP2   (output) DOUBLE PRECISION array, dimension (M-P)
                    144: *          The scalar factors of the elementary reflectors that define
                    145: *          P2.
                    146: *
                    147: *  TAUQ1   (output) DOUBLE PRECISION array, dimension (Q)
                    148: *          The scalar factors of the elementary reflectors that define
                    149: *          Q1.
                    150: *
                    151: *  TAUQ2   (output) DOUBLE PRECISION array, dimension (M-Q)
                    152: *          The scalar factors of the elementary reflectors that define
                    153: *          Q2.
                    154: *
                    155: *  WORK    (workspace) DOUBLE PRECISION array, dimension (LWORK)
                    156: *
                    157: *  LWORK   (input) INTEGER
                    158: *          The dimension of the array WORK. LWORK >= M-Q.
                    159: *
                    160: *          If LWORK = -1, then a workspace query is assumed; the routine
                    161: *          only calculates the optimal size of the WORK array, returns
                    162: *          this value as the first entry of the WORK array, and no error
                    163: *          message related to LWORK is issued by XERBLA.
                    164: *
                    165: *  INFO    (output) INTEGER
                    166: *          = 0:  successful exit.
                    167: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    168: *
                    169: *  Further Details
                    170: *  ===============
                    171: *
                    172: *  The bidiagonal blocks B11, B12, B21, and B22 are represented
                    173: *  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
                    174: *  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
                    175: *  lower bidiagonal. Every entry in each bidiagonal band is a product
                    176: *  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
                    177: *  [1] or DORCSD for details.
                    178: *
                    179: *  P1, P2, Q1, and Q2 are represented as products of elementary
                    180: *  reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2
                    181: *  using DORGQR and DORGLQ.
                    182: *
                    183: *  Reference
                    184: *  =========
                    185: *
                    186: *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
                    187: *      Algorithms, 50(1):33-65, 2009.
                    188: *
                    189: *  ====================================================================
                    190: *
                    191: *     .. Parameters ..
                    192:       DOUBLE PRECISION   REALONE
                    193:       PARAMETER          ( REALONE = 1.0D0 )
                    194:       DOUBLE PRECISION   NEGONE, ONE
                    195:       PARAMETER          ( NEGONE = -1.0D0, ONE = 1.0D0 )
                    196: *     ..
                    197: *     .. Local Scalars ..
                    198:       LOGICAL            COLMAJOR, LQUERY
                    199:       INTEGER            I, LWORKMIN, LWORKOPT
                    200:       DOUBLE PRECISION   Z1, Z2, Z3, Z4
                    201: *     ..
                    202: *     .. External Subroutines ..
                    203:       EXTERNAL           DAXPY, DLARF, DLARFGP, DSCAL, XERBLA
                    204: *     ..
                    205: *     .. External Functions ..
                    206:       DOUBLE PRECISION   DNRM2
                    207:       LOGICAL            LSAME
                    208:       EXTERNAL           DNRM2, LSAME
                    209: *     ..
                    210: *     .. Intrinsic Functions
                    211:       INTRINSIC          ATAN2, COS, MAX, MIN, SIN
                    212: *     ..
                    213: *     .. Executable Statements ..
                    214: *
                    215: *     Test input arguments
                    216: *
                    217:       INFO = 0
                    218:       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
                    219:       IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
                    220:          Z1 = REALONE
                    221:          Z2 = REALONE
                    222:          Z3 = REALONE
                    223:          Z4 = REALONE
                    224:       ELSE
                    225:          Z1 = REALONE
                    226:          Z2 = -REALONE
                    227:          Z3 = REALONE
                    228:          Z4 = -REALONE
                    229:       END IF
                    230:       LQUERY = LWORK .EQ. -1
                    231: *
                    232:       IF( M .LT. 0 ) THEN
                    233:          INFO = -3
                    234:       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
                    235:          INFO = -4
                    236:       ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
                    237:      $         Q .GT. M-Q ) THEN
                    238:          INFO = -5
                    239:       ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
                    240:          INFO = -7
                    241:       ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
                    242:          INFO = -7
                    243:       ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
                    244:          INFO = -9
                    245:       ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
                    246:          INFO = -9
                    247:       ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
                    248:          INFO = -11
                    249:       ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
                    250:          INFO = -11
                    251:       ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
                    252:          INFO = -13
                    253:       ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
                    254:          INFO = -13
                    255:       END IF
                    256: *
                    257: *     Compute workspace
                    258: *
                    259:       IF( INFO .EQ. 0 ) THEN
                    260:          LWORKOPT = M - Q
                    261:          LWORKMIN = M - Q
                    262:          WORK(1) = LWORKOPT
                    263:          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
                    264:             INFO = -21
                    265:          END IF
                    266:       END IF
                    267:       IF( INFO .NE. 0 ) THEN
                    268:          CALL XERBLA( 'xORBDB', -INFO )
                    269:          RETURN
                    270:       ELSE IF( LQUERY ) THEN
                    271:          RETURN
                    272:       END IF
                    273: *
                    274: *     Handle column-major and row-major separately
                    275: *
                    276:       IF( COLMAJOR ) THEN
                    277: *
                    278: *        Reduce columns 1, ..., Q of X11, X12, X21, and X22 
                    279: *
                    280:          DO I = 1, Q
                    281: *
                    282:             IF( I .EQ. 1 ) THEN
                    283:                CALL DSCAL( P-I+1, Z1, X11(I,I), 1 )
                    284:             ELSE
                    285:                CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 )
                    286:                CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1),
                    287:      $                     1, X11(I,I), 1 )
                    288:             END IF
                    289:             IF( I .EQ. 1 ) THEN
                    290:                CALL DSCAL( M-P-I+1, Z2, X21(I,I), 1 )
                    291:             ELSE
                    292:                CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 )
                    293:                CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1),
                    294:      $                     1, X21(I,I), 1 )
                    295:             END IF
                    296: *
                    297:             THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ),
                    298:      $                 DNRM2( P-I+1, X11(I,I), 1 ) )
                    299: *
                    300:             CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
                    301:             X11(I,I) = ONE
                    302:             CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
                    303:             X21(I,I) = ONE
                    304: *
                    305:             CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
                    306:      $                  X11(I,I+1), LDX11, WORK )
                    307:             CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
                    308:      $                  X12(I,I), LDX12, WORK )
                    309:             CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
                    310:      $                  X21(I,I+1), LDX21, WORK )
                    311:             CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
                    312:      $                  X22(I,I), LDX22, WORK )
                    313: *
                    314:             IF( I .LT. Q ) THEN
                    315:                CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
                    316:      $                     LDX11 )
                    317:                CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21,
                    318:      $                     X11(I,I+1), LDX11 )
                    319:             END IF
                    320:             CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 )
                    321:             CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22,
                    322:      $                  X12(I,I), LDX12 )
                    323: *
                    324:             IF( I .LT. Q )
                    325:      $         PHI(I) = ATAN2( DNRM2( Q-I, X11(I,I+1), LDX11 ),
                    326:      $                  DNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
                    327: *
                    328:             IF( I .LT. Q ) THEN
                    329:                CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
                    330:      $                       TAUQ1(I) )
                    331:                X11(I,I+1) = ONE
                    332:             END IF
                    333:             CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
                    334:      $                    TAUQ2(I) )
                    335:             X12(I,I) = ONE
                    336: *
                    337:             IF( I .LT. Q ) THEN
                    338:                CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
                    339:      $                     X11(I+1,I+1), LDX11, WORK )
                    340:                CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
                    341:      $                     X21(I+1,I+1), LDX21, WORK )
                    342:             END IF
                    343:             CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
                    344:      $                  X12(I+1,I), LDX12, WORK )
                    345:             CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
                    346:      $                  X22(I+1,I), LDX22, WORK )
                    347: *
                    348:          END DO
                    349: *
                    350: *        Reduce columns Q + 1, ..., P of X12, X22
                    351: *
                    352:          DO I = Q + 1, P
                    353: *
                    354:             CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
                    355:             CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
                    356:      $                    TAUQ2(I) )
                    357:             X12(I,I) = ONE
                    358: *
                    359:             CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
                    360:      $                  X12(I+1,I), LDX12, WORK )
                    361:             IF( M-P-Q .GE. 1 )
                    362:      $         CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
                    363:      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
                    364: *
                    365:          END DO
                    366: *
                    367: *        Reduce columns P + 1, ..., M - Q of X12, X22
                    368: *
                    369:          DO I = 1, M - P - Q
                    370: *
                    371:             CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
                    372:             CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
                    373:      $                    LDX22, TAUQ2(P+I) )
                    374:             X22(Q+I,P+I) = ONE
                    375:             CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
                    376:      $                  TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
                    377: *
                    378:          END DO
                    379: *
                    380:       ELSE
                    381: *
                    382: *        Reduce columns 1, ..., Q of X11, X12, X21, X22
                    383: *
                    384:          DO I = 1, Q
                    385: *
                    386:             IF( I .EQ. 1 ) THEN
                    387:                CALL DSCAL( P-I+1, Z1, X11(I,I), LDX11 )
                    388:             ELSE
                    389:                CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 )
                    390:                CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I),
                    391:      $                     LDX12, X11(I,I), LDX11 )
                    392:             END IF
                    393:             IF( I .EQ. 1 ) THEN
                    394:                CALL DSCAL( M-P-I+1, Z2, X21(I,I), LDX21 )
                    395:             ELSE
                    396:                CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 )
                    397:                CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I),
                    398:      $                     LDX22, X21(I,I), LDX21 )
                    399:             END IF
                    400: *
                    401:             THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), LDX21 ),
                    402:      $                 DNRM2( P-I+1, X11(I,I), LDX11 ) )
                    403: *
                    404:             CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
                    405:             X11(I,I) = ONE
                    406:             CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
                    407:      $                    TAUP2(I) )
                    408:             X21(I,I) = ONE
                    409: *
                    410:             CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
                    411:      $                  X11(I+1,I), LDX11, WORK )
                    412:             CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
                    413:      $                  X12(I,I), LDX12, WORK )
                    414:             CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
                    415:      $                  X21(I+1,I), LDX21, WORK )
                    416:             CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
                    417:      $                  TAUP2(I), X22(I,I), LDX22, WORK )
                    418: *
                    419:             IF( I .LT. Q ) THEN
                    420:                CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
                    421:                CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1,
                    422:      $                     X11(I+1,I), 1 )
                    423:             END IF
                    424:             CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 )
                    425:             CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1,
                    426:      $                  X12(I,I), 1 )
                    427: *
                    428:             IF( I .LT. Q )
                    429:      $         PHI(I) = ATAN2( DNRM2( Q-I, X11(I+1,I), 1 ),
                    430:      $                  DNRM2( M-Q-I+1, X12(I,I), 1 ) )
                    431: *
                    432:             IF( I .LT. Q ) THEN
                    433:                CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
                    434:                X11(I+1,I) = ONE
                    435:             END IF
                    436:             CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
                    437:             X12(I,I) = ONE
                    438: *
                    439:             IF( I .LT. Q ) THEN
                    440:                CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I),
                    441:      $                     X11(I+1,I+1), LDX11, WORK )
                    442:                CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I),
                    443:      $                     X21(I+1,I+1), LDX21, WORK )
                    444:             END IF
                    445:             CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
                    446:      $                  X12(I,I+1), LDX12, WORK )
                    447:             CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
                    448:      $                  X22(I,I+1), LDX22, WORK )
                    449: *
                    450:          END DO
                    451: *
                    452: *        Reduce columns Q + 1, ..., P of X12, X22
                    453: *
                    454:          DO I = Q + 1, P
                    455: *
                    456:             CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 )
                    457:             CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
                    458:             X12(I,I) = ONE
                    459: *
                    460:             CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
                    461:      $                  X12(I,I+1), LDX12, WORK )
                    462:             IF( M-P-Q .GE. 1 )
                    463:      $         CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
                    464:      $                     X22(I,Q+1), LDX22, WORK )
                    465: *
                    466:          END DO
                    467: *
                    468: *        Reduce columns P + 1, ..., M - Q of X12, X22
                    469: *
                    470:          DO I = 1, M - P - Q
                    471: *
                    472:             CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
                    473:             CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
                    474:      $                    TAUQ2(P+I) )
                    475:             X22(P+I,Q+I) = ONE
                    476: *
                    477:             CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
                    478:      $                  TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
                    479: *
                    480:          END DO
                    481: *
                    482:       END IF
                    483: *
                    484:       RETURN
                    485: *
                    486: *     End of DORBDB
                    487: *
                    488:       END
                    489: 

CVSweb interface <joel.bertrand@systella.fr>