LAPACK: SRC/sorbdb.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
00002      $                   X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
00003      $                   TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
00004       IMPLICIT NONE
00005 *
00006 *  -- LAPACK routine (version 3.3.0) --
00007 *
00008 *  -- Contributed by Brian Sutton of the Randolph-Macon College --
00009 *  -- November 2010
00010 *
00011 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00012 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--     
00013 *
00014 *     .. Scalar Arguments ..
00015       CHARACTER          SIGNS, TRANS
00016       INTEGER            INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
00017      $                   Q
00018 *     ..
00019 *     .. Array Arguments ..
00020       REAL               PHI( * ), THETA( * )
00021       REAL               TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
00022      $                   WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
00023      $                   X21( LDX21, * ), X22( LDX22, * )
00024 *     ..
00025 *
00026 *  Purpose
00027 *  =======
00028 *
00029 *  SORBDB simultaneously bidiagonalizes the blocks of an M-by-M
00030 *  partitioned orthogonal matrix X:
00031 *
00032 *                                  [ B11 | B12 0  0 ]
00033 *      [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**T
00034 *  X = [-----------] = [---------] [----------------] [---------]   .
00035 *      [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
00036 *                                  [  0  |  0  0  I ]
00037 *
00038 *  X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
00039 *  not the case, then X must be transposed and/or permuted. This can be
00040 *  done in constant time using the TRANS and SIGNS options. See SORCSD
00041 *  for details.)
00042 *
00043 *  The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
00044 *  (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
00045 *  represented implicitly by Householder vectors.
00046 *
00047 *  B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
00048 *  implicitly by angles THETA, PHI.
00049 *
00050 *  Arguments
00051 *  =========
00052 *
00053 *  TRANS   (input) CHARACTER
00054 *          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
00055 *                      order;
00056 *          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
00057 *                      major order.
00058 *
00059 *  SIGNS   (input) CHARACTER
00060 *          = 'O':      The lower-left block is made nonpositive (the
00061 *                      "other" convention);
00062 *          otherwise:  The upper-right block is made nonpositive (the
00063 *                      "default" convention).
00064 *
00065 *  M       (input) INTEGER
00066 *          The number of rows and columns in X.
00067 *
00068 *  P       (input) INTEGER
00069 *          The number of rows in X11 and X12. 0 <= P <= M.
00070 *
00071 *  Q       (input) INTEGER
00072 *          The number of columns in X11 and X21. 0 <= Q <=
00073 *          MIN(P,M-P,M-Q).
00074 *
00075 *  X11     (input/output) REAL array, dimension (LDX11,Q)
00076 *          On entry, the top-left block of the orthogonal matrix to be
00077 *          reduced. On exit, the form depends on TRANS:
00078 *          If TRANS = 'N', then
00079 *             the columns of tril(X11) specify reflectors for P1,
00080 *             the rows of triu(X11,1) specify reflectors for Q1;
00081 *          else TRANS = 'T', and
00082 *             the rows of triu(X11) specify reflectors for P1,
00083 *             the columns of tril(X11,-1) specify reflectors for Q1.
00084 *
00085 *  LDX11   (input) INTEGER
00086 *          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
00087 *          P; else LDX11 >= Q.
00088 *
00089 *  X12     (input/output) REAL array, dimension (LDX12,M-Q)
00090 *          On entry, the top-right block of the orthogonal matrix to
00091 *          be reduced. On exit, the form depends on TRANS:
00092 *          If TRANS = 'N', then
00093 *             the rows of triu(X12) specify the first P reflectors for
00094 *             Q2;
00095 *          else TRANS = 'T', and
00096 *             the columns of tril(X12) specify the first P reflectors
00097 *             for Q2.
00098 *
00099 *  LDX12   (input) INTEGER
00100 *          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
00101 *          P; else LDX11 >= M-Q.
00102 *
00103 *  X21     (input/output) REAL array, dimension (LDX21,Q)
00104 *          On entry, the bottom-left block of the orthogonal matrix to
00105 *          be reduced. On exit, the form depends on TRANS:
00106 *          If TRANS = 'N', then
00107 *             the columns of tril(X21) specify reflectors for P2;
00108 *          else TRANS = 'T', and
00109 *             the rows of triu(X21) specify reflectors for P2.
00110 *
00111 *  LDX21   (input) INTEGER
00112 *          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
00113 *          M-P; else LDX21 >= Q.
00114 *
00115 *  X22     (input/output) REAL array, dimension (LDX22,M-Q)
00116 *          On entry, the bottom-right block of the orthogonal matrix to
00117 *          be reduced. On exit, the form depends on TRANS:
00118 *          If TRANS = 'N', then
00119 *             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
00120 *             M-P-Q reflectors for Q2,
00121 *          else TRANS = 'T', and
00122 *             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
00123 *             M-P-Q reflectors for P2.
00124 *
00125 *  LDX22   (input) INTEGER
00126 *          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
00127 *          M-P; else LDX22 >= M-Q.
00128 *
00129 *  THETA   (output) REAL array, dimension (Q)
00130 *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
00131 *          be computed from the angles THETA and PHI. See Further
00132 *          Details.
00133 *
00134 *  PHI     (output) REAL array, dimension (Q-1)
00135 *          The entries of the bidiagonal blocks B11, B12, B21, B22 can
00136 *          be computed from the angles THETA and PHI. See Further
00137 *          Details.
00138 *
00139 *  TAUP1   (output) REAL array, dimension (P)
00140 *          The scalar factors of the elementary reflectors that define
00141 *          P1.
00142 *
00143 *  TAUP2   (output) REAL array, dimension (M-P)
00144 *          The scalar factors of the elementary reflectors that define
00145 *          P2.
00146 *
00147 *  TAUQ1   (output) REAL array, dimension (Q)
00148 *          The scalar factors of the elementary reflectors that define
00149 *          Q1.
00150 *
00151 *  TAUQ2   (output) REAL array, dimension (M-Q)
00152 *          The scalar factors of the elementary reflectors that define
00153 *          Q2.
00154 *
00155 *  WORK    (workspace) REAL array, dimension (LWORK)
00156 *
00157 *  LWORK   (input) INTEGER
00158 *          The dimension of the array WORK. LWORK >= M-Q.
00159 *
00160 *          If LWORK = -1, then a workspace query is assumed; the routine
00161 *          only calculates the optimal size of the WORK array, returns
00162 *          this value as the first entry of the WORK array, and no error
00163 *          message related to LWORK is issued by XERBLA.
00164 *
00165 *  INFO    (output) INTEGER
00166 *          = 0:  successful exit.
00167 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00168 *
00169 *  Further Details
00170 *  ===============
00171 *
00172 *  The bidiagonal blocks B11, B12, B21, and B22 are represented
00173 *  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
00174 *  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
00175 *  lower bidiagonal. Every entry in each bidiagonal band is a product
00176 *  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
00177 *  [1] or SORCSD for details.
00178 *
00179 *  P1, P2, Q1, and Q2 are represented as products of elementary
00180 *  reflectors. See SORCSD for details on generating P1, P2, Q1, and Q2
00181 *  using SORGQR and SORGLQ.
00182 *
00183 *  Reference
00184 *  =========
00185 *
00186 *  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
00187 *      Algorithms, 50(1):33-65, 2009.
00188 *
00189 *  ====================================================================
00190 *
00191 *     .. Parameters ..
00192       REAL               REALONE
00193       PARAMETER          ( REALONE = 1.0E0 )
00194       REAL               NEGONE, ONE
00195       PARAMETER          ( NEGONE = -1.0E0, ONE = 1.0E0 )
00196 *     ..
00197 *     .. Local Scalars ..
00198       LOGICAL            COLMAJOR, LQUERY
00199       INTEGER            I, LWORKMIN, LWORKOPT
00200       REAL               Z1, Z2, Z3, Z4
00201 *     ..
00202 *     .. External Subroutines ..
00203       EXTERNAL           SAXPY, SLARF, SLARFGP, SSCAL, XERBLA
00204 *     ..
00205 *     .. External Functions ..
00206       REAL               SNRM2
00207       LOGICAL            LSAME
00208       EXTERNAL           SNRM2, LSAME
00209 *     ..
00210 *     .. Intrinsic Functions
00211       INTRINSIC          ATAN2, COS, MAX, MIN, SIN
00212 *     ..
00213 *     .. Executable Statements ..
00214 *
00215 *     Test input arguments
00216 *
00217       INFO = 0
00218       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
00219       IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
00220          Z1 = REALONE
00221          Z2 = REALONE
00222          Z3 = REALONE
00223          Z4 = REALONE
00224       ELSE
00225          Z1 = REALONE
00226          Z2 = -REALONE
00227          Z3 = REALONE
00228          Z4 = -REALONE
00229       END IF
00230       LQUERY = LWORK .EQ. -1
00231 *
00232       IF( M .LT. 0 ) THEN
00233          INFO = -3
00234       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
00235          INFO = -4
00236       ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
00237      $         Q .GT. M-Q ) THEN
00238          INFO = -5
00239       ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
00240          INFO = -7
00241       ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
00242          INFO = -7
00243       ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
00244          INFO = -9
00245       ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
00246          INFO = -9
00247       ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
00248          INFO = -11
00249       ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
00250          INFO = -11
00251       ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
00252          INFO = -13
00253       ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
00254          INFO = -13
00255       END IF
00256 *
00257 *     Compute workspace
00258 *
00259       IF( INFO .EQ. 0 ) THEN
00260          LWORKOPT = M - Q
00261          LWORKMIN = M - Q
00262          WORK(1) = LWORKOPT
00263          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
00264             INFO = -21
00265          END IF
00266       END IF
00267       IF( INFO .NE. 0 ) THEN
00268          CALL XERBLA( 'xORBDB', -INFO )
00269          RETURN
00270       ELSE IF( LQUERY ) THEN
00271          RETURN
00272       END IF
00273 *
00274 *     Handle column-major and row-major separately
00275 *
00276       IF( COLMAJOR ) THEN
00277 *
00278 *        Reduce columns 1, ..., Q of X11, X12, X21, and X22 
00279 *
00280          DO I = 1, Q
00281 *
00282             IF( I .EQ. 1 ) THEN
00283                CALL SSCAL( P-I+1, Z1, X11(I,I), 1 )
00284             ELSE
00285                CALL SSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 )
00286                CALL SAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1),
00287      $                     1, X11(I,I), 1 )
00288             END IF
00289             IF( I .EQ. 1 ) THEN
00290                CALL SSCAL( M-P-I+1, Z2, X21(I,I), 1 )
00291             ELSE
00292                CALL SSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 )
00293                CALL SAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1),
00294      $                     1, X21(I,I), 1 )
00295             END IF
00296 *
00297             THETA(I) = ATAN2( SNRM2( M-P-I+1, X21(I,I), 1 ),
00298      $                 SNRM2( P-I+1, X11(I,I), 1 ) )
00299 *
00300             CALL SLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
00301             X11(I,I) = ONE
00302             CALL SLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
00303             X21(I,I) = ONE
00304 *
00305             CALL SLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
00306      $                  X11(I,I+1), LDX11, WORK )
00307             CALL SLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
00308      $                  X12(I,I), LDX12, WORK )
00309             CALL SLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
00310      $                  X21(I,I+1), LDX21, WORK )
00311             CALL SLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
00312      $                  X22(I,I), LDX22, WORK )
00313 *
00314             IF( I .LT. Q ) THEN
00315                CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
00316      $                     LDX11 )
00317                CALL SAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21,
00318      $                     X11(I,I+1), LDX11 )
00319             END IF
00320             CALL SSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 )
00321             CALL SAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22,
00322      $                  X12(I,I), LDX12 )
00323 *
00324             IF( I .LT. Q )
00325      $         PHI(I) = ATAN2( SNRM2( Q-I, X11(I,I+1), LDX11 ),
00326      $                  SNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
00327 *
00328             IF( I .LT. Q ) THEN
00329                CALL SLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
00330      $                       TAUQ1(I) )
00331                X11(I,I+1) = ONE
00332             END IF
00333             CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
00334      $                    TAUQ2(I) )
00335             X12(I,I) = ONE
00336 *
00337             IF( I .LT. Q ) THEN
00338                CALL SLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
00339      $                     X11(I+1,I+1), LDX11, WORK )
00340                CALL SLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
00341      $                     X21(I+1,I+1), LDX21, WORK )
00342             END IF
00343             CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00344      $                  X12(I+1,I), LDX12, WORK )
00345             CALL SLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00346      $                  X22(I+1,I), LDX22, WORK )
00347 *
00348          END DO
00349 *
00350 *        Reduce columns Q + 1, ..., P of X12, X22
00351 *
00352          DO I = Q + 1, P
00353 *
00354             CALL SSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
00355             CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
00356      $                    TAUQ2(I) )
00357             X12(I,I) = ONE
00358 *
00359             CALL SLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
00360      $                  X12(I+1,I), LDX12, WORK )
00361             IF( M-P-Q .GE. 1 )
00362      $         CALL SLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
00363      $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
00364 *
00365          END DO
00366 *
00367 *        Reduce columns P + 1, ..., M - Q of X12, X22
00368 *
00369          DO I = 1, M - P - Q
00370 *
00371             CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
00372             CALL SLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
00373      $                    LDX22, TAUQ2(P+I) )
00374             X22(Q+I,P+I) = ONE
00375             CALL SLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
00376      $                  TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
00377 *
00378          END DO
00379 *
00380       ELSE
00381 *
00382 *        Reduce columns 1, ..., Q of X11, X12, X21, X22
00383 *
00384          DO I = 1, Q
00385 *
00386             IF( I .EQ. 1 ) THEN
00387                CALL SSCAL( P-I+1, Z1, X11(I,I), LDX11 )
00388             ELSE
00389                CALL SSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 )
00390                CALL SAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I),
00391      $                     LDX12, X11(I,I), LDX11 )
00392             END IF
00393             IF( I .EQ. 1 ) THEN
00394                CALL SSCAL( M-P-I+1, Z2, X21(I,I), LDX21 )
00395             ELSE
00396                CALL SSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 )
00397                CALL SAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I),
00398      $                     LDX22, X21(I,I), LDX21 )
00399             END IF
00400 *
00401             THETA(I) = ATAN2( SNRM2( M-P-I+1, X21(I,I), LDX21 ),
00402      $                 SNRM2( P-I+1, X11(I,I), LDX11 ) )
00403 *
00404             CALL SLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
00405             X11(I,I) = ONE
00406             CALL SLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
00407      $                    TAUP2(I) )
00408             X21(I,I) = ONE
00409 *
00410             CALL SLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
00411      $                  X11(I+1,I), LDX11, WORK )
00412             CALL SLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
00413      $                  X12(I,I), LDX12, WORK )
00414             CALL SLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
00415      $                  X21(I+1,I), LDX21, WORK )
00416             CALL SLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
00417      $                  TAUP2(I), X22(I,I), LDX22, WORK )
00418 *
00419             IF( I .LT. Q ) THEN
00420                CALL SSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
00421                CALL SAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1,
00422      $                     X11(I+1,I), 1 )
00423             END IF
00424             CALL SSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 )
00425             CALL SAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1,
00426      $                  X12(I,I), 1 )
00427 *
00428             IF( I .LT. Q )
00429      $         PHI(I) = ATAN2( SNRM2( Q-I, X11(I+1,I), 1 ),
00430      $                  SNRM2( M-Q-I+1, X12(I,I), 1 ) )
00431 *
00432             IF( I .LT. Q ) THEN
00433                CALL SLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
00434                X11(I+1,I) = ONE
00435             END IF
00436             CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
00437             X12(I,I) = ONE
00438 *
00439             IF( I .LT. Q ) THEN
00440                CALL SLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I),
00441      $                     X11(I+1,I+1), LDX11, WORK )
00442                CALL SLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I),
00443      $                     X21(I+1,I+1), LDX21, WORK )
00444             END IF
00445             CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
00446      $                  X12(I,I+1), LDX12, WORK )
00447             CALL SLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
00448      $                  X22(I,I+1), LDX22, WORK )
00449 *
00450          END DO
00451 *
00452 *        Reduce columns Q + 1, ..., P of X12, X22
00453 *
00454          DO I = Q + 1, P
00455 *
00456             CALL SSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 )
00457             CALL SLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
00458             X12(I,I) = ONE
00459 *
00460             CALL SLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
00461      $                  X12(I,I+1), LDX12, WORK )
00462             IF( M-P-Q .GE. 1 )
00463      $         CALL SLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
00464      $                     X22(I,Q+1), LDX22, WORK )
00465 *
00466          END DO
00467 *
00468 *        Reduce columns P + 1, ..., M - Q of X12, X22
00469 *
00470          DO I = 1, M - P - Q
00471 *
00472             CALL SSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
00473             CALL SLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
00474      $                    TAUQ2(P+I) )
00475             X22(P+I,Q+I) = ONE
00476 *
00477             CALL SLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
00478      $                  TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
00479 *
00480          END DO
00481 *
00482       END IF
00483 *
00484       RETURN
00485 *
00486 *     End of SORBDB
00487 *
00488       END
00489