SRC/dorm22.f

   1 *> \brief \b DORM22 multiplies a general matrix by a banded orthogonal matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DORM22 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorm22.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorm22.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorm22.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *     SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
  22 *    $                   WORK, LWORK, INFO )
  23 *
  24 *     .. Scalar Arguments ..
  25 *     CHARACTER          SIDE, TRANS
  26 *     INTEGER            M, N, N1, N2, LDQ, LDC, LWORK, INFO
  27 *     ..
  28 *     .. Array Arguments ..
  29 *     DOUBLE PRECISION   Q( LDQ, * ), C( LDC, * ), WORK( * )
  30 *     ..
  31 *
  32 *> \par Purpose
  33 *  ============
  34 *>
  35 *> \verbatim
  36 *>
  37 *>
  38 *>  DORM22 overwrites the general real M-by-N matrix C with
  39 *>
  40 *>                  SIDE = 'L'     SIDE = 'R'
  41 *>  TRANS = 'N':      Q * C          C * Q
  42 *>  TRANS = 'T':      Q**T * C       C * Q**T
  43 *>
  44 *>  where Q is a real orthogonal matrix of order NQ, with NQ = M if
  45 *>  SIDE = 'L' and NQ = N if SIDE = 'R'.
  46 *>  The orthogonal matrix Q processes a 2-by-2 block structure
  47 *>
  48 *>         [  Q11  Q12  ]
  49 *>     Q = [            ]
  50 *>         [  Q21  Q22  ],
  51 *>
  52 *>  where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
  53 *>  N2-by-N2 upper triangular matrix.
  54 *> \endverbatim
  55 *
  56 *  Arguments
  57 *  =========
  58 *
  59 *> \param[in] SIDE
  60 *> \verbatim
  61 *>          SIDE is CHARACTER*1
  62 *>          = 'L': apply Q or Q**T from the Left;
  63 *>          = 'R': apply Q or Q**T from the Right.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] TRANS
  67 *> \verbatim
  68 *>          TRANS is CHARACTER*1
  69 *>          = 'N':  apply Q (No transpose);
  70 *>          = 'C':  apply Q**T (Conjugate transpose).
  71 *> \endverbatim
  72 *>
  73 *> \param[in] M
  74 *> \verbatim
  75 *>          M is INTEGER
  76 *>          The number of rows of the matrix C. M >= 0.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] N
  80 *> \verbatim
  81 *>          N is INTEGER
  82 *>          The number of columns of the matrix C. N >= 0.
  83 *> \endverbatim
  84 *>
  85 *> \param[in] N1
  86 *> \param[in] N2
  87 *> \verbatim
  88 *>          N1 is INTEGER
  89 *>          N2 is INTEGER
  90 *>          The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
  91 *>          The following requirement must be satisfied:
  92 *>          N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] Q
  96 *> \verbatim
  97 *>          Q is DOUBLE PRECISION array, dimension
  98 *>                                       (LDQ,M) if SIDE = 'L'
  99 *>                                       (LDQ,N) if SIDE = 'R'
 100 *> \endverbatim
 101 *>
 102 *> \param[in] LDQ
 103 *> \verbatim
 104 *>          LDQ is INTEGER
 105 *>          The leading dimension of the array Q.
 106 *>          LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
 107 *> \endverbatim
 108 *>
 109 *> \param[in,out] C
 110 *> \verbatim
 111 *>          C is DOUBLE PRECISION array, dimension (LDC,N)
 112 *>          On entry, the M-by-N matrix C.
 113 *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 114 *> \endverbatim
 115 *>
 116 *> \param[in] LDC
 117 *> \verbatim
 118 *>          LDC is INTEGER
 119 *>          The leading dimension of the array C. LDC >= max(1,M).
 120 *> \endverbatim
 121 *>
 122 *> \param[out] WORK
 123 *> \verbatim
 124 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 125 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LWORK
 129 *> \verbatim
 130 *>          LWORK is INTEGER
 131 *>          The dimension of the array WORK.
 132 *>          If SIDE = 'L', LWORK >= max(1,N);
 133 *>          if SIDE = 'R', LWORK >= max(1,M).
 134 *>          For optimum performance LWORK >= M*N.
 135 *>
 136 *>          If LWORK = -1, then a workspace query is assumed; the routine
 137 *>          only calculates the optimal size of the WORK array, returns
 138 *>          this value as the first entry of the WORK array, and no error
 139 *>          message related to LWORK is issued by XERBLA.
 140 *> \endverbatim
 141 *>
 142 *> \param[out] INFO
 143 *> \verbatim
 144 *>          INFO is INTEGER
 145 *>          = 0:  successful exit
 146 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 147 *> \endverbatim
 148 *
 149 *
 150 *  Authors:
 151 *  ========
 152 *
 153 *> \author Univ. of Tennessee
 154 *> \author Univ. of California Berkeley
 155 *> \author Univ. of Colorado Denver
 156 *> \author NAG Ltd.
 157 *
 158 *> \date January 2015
 159 *
 160 *> \ingroup complexOTHERcomputational
 161 *
 162 *  =====================================================================
 163       SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
 164      $                   WORK, LWORK, INFO )
 165 *
 166 *  -- LAPACK computational routine (version 3.6.0) --
 167 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 168 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 169 *     January 2015
 170 *
 171       IMPLICIT NONE
 172 *
 173 *     .. Scalar Arguments ..
 174       CHARACTER          SIDE, TRANS
 175       INTEGER            M, N, N1, N2, LDQ, LDC, LWORK, INFO
 176 *     ..
 177 *     .. Array Arguments ..
 178       DOUBLE PRECISION   Q( LDQ, * ), C( LDC, * ), WORK( * )
 179 *     ..
 180 *
 181 *  =====================================================================
 182 *
 183 *     .. Parameters ..
 184       DOUBLE PRECISION   ONE
 185       PARAMETER          ( ONE = 1.0D+0 )
 186 *
 187 *     .. Local Scalars ..
 188       LOGICAL            LEFT, LQUERY, NOTRAN
 189       INTEGER            I, LDWORK, LEN, LWKOPT, NB, NQ, NW
 190 *     ..
 191 *     .. External Functions ..
 192       LOGICAL            LSAME
 193       EXTERNAL           LSAME
 194 *     ..
 195 *     .. External Subroutines ..
 196       EXTERNAL           DGEMM, DLACPY, DTRMM, XERBLA
 197 *     ..
 198 *     .. Intrinsic Functions ..
 199       INTRINSIC          DBLE, MAX, MIN
 200 *     ..
 201 *     .. Executable Statements ..
 202 *
 203 *     Test the input arguments
 204 *
 205       INFO = 0
 206       LEFT = LSAME( SIDE, 'L' )
 207       NOTRAN = LSAME( TRANS, 'N' )
 208       LQUERY = ( LWORK.EQ.-1 )
 209 *
 210 *     NQ is the order of Q;
 211 *     NW is the minimum dimension of WORK.
 212 *
 213       IF( LEFT ) THEN
 214          NQ = M
 215       ELSE
 216          NQ = N
 217       END IF
 218       NW = NQ
 219       IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
 220       IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
 221          INFO = -1
 222       ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
 223      $          THEN
 224          INFO = -2
 225       ELSE IF( M.LT.0 ) THEN
 226          INFO = -3
 227       ELSE IF( N.LT.0 ) THEN
 228          INFO = -4
 229       ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
 230          INFO = -5
 231       ELSE IF( N2.LT.0 ) THEN
 232          INFO = -6
 233       ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
 234          INFO = -8
 235       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
 236          INFO = -10
 237       ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
 238          INFO = -12
 239       END IF
 240 *
 241       IF( INFO.EQ.0 ) THEN
 242          LWKOPT = M*N
 243          WORK( 1 ) = DBLE( LWKOPT )
 244       END IF
 245 *
 246       IF( INFO.NE.0 ) THEN
 247          CALL XERBLA( 'DORM22', -INFO )
 248          RETURN
 249       ELSE IF( LQUERY ) THEN
 250          RETURN
 251       END IF
 252 *
 253 *     Quick return if possible
 254 *
 255       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 256          WORK( 1 ) = 1
 257          RETURN
 258       END IF
 259 *
 260 *     Degenerate cases (N1 = 0 or N2 = 0) are handled using DTRMM.
 261 *
 262       IF( N1.EQ.0 ) THEN
 263          CALL DTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
 264      $               Q, LDQ, C, LDC )
 265          WORK( 1 ) = ONE
 266          RETURN
 267       ELSE IF( N2.EQ.0 ) THEN
 268          CALL DTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
 269      $               Q, LDQ, C, LDC )
 270          WORK( 1 ) = ONE
 271          RETURN
 272       END IF
 273 *
 274 *     Compute the largest chunk size available from the workspace.
 275 *
 276       NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
 277 *
 278       IF( LEFT ) THEN
 279          IF( NOTRAN ) THEN
 280             DO I = 1, N, NB
 281                LEN = MIN( NB, N-I+1 )
 282                LDWORK = M
 283 *
 284 *              Multiply bottom part of C by Q12.
 285 *
 286                CALL DLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
 287      $                      LDWORK )
 288                CALL DTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
 289      $                     N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
 290      $                     LDWORK )
 291 *
 292 *              Multiply top part of C by Q11.
 293 *
 294                CALL DGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
 295      $                     ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
 296      $                     LDWORK )
 297 *
 298 *              Multiply top part of C by Q21.
 299 *
 300                CALL DLACPY( 'All', N2, LEN, C( 1, I ), LDC,
 301      $                      WORK( N1+1 ), LDWORK )
 302                CALL DTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
 303      $                     N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
 304      $                     WORK( N1+1 ), LDWORK )
 305 *
 306 *              Multiply bottom part of C by Q22.
 307 *
 308                CALL DGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
 309      $                     ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
 310      $                     ONE, WORK( N1+1 ), LDWORK )
 311 *
 312 *              Copy everything back.
 313 *
 314                CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
 315      $                      LDC )
 316             END DO
 317          ELSE
 318             DO I = 1, N, NB
 319                LEN = MIN( NB, N-I+1 )
 320                LDWORK = M
 321 *
 322 *              Multiply bottom part of C by Q21**T.
 323 *
 324                CALL DLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
 325      $                      LDWORK )
 326                CALL DTRMM( 'Left', 'Upper', 'Transpose', 'Non-Unit',
 327      $                     N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
 328      $                     LDWORK )
 329 *
 330 *              Multiply top part of C by Q11**T.
 331 *
 332                CALL DGEMM( 'Transpose', 'No Transpose', N2, LEN, N1,
 333      $                     ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
 334      $                     LDWORK )
 335 *
 336 *              Multiply top part of C by Q12**T.
 337 *
 338                CALL DLACPY( 'All', N1, LEN, C( 1, I ), LDC,
 339      $                      WORK( N2+1 ), LDWORK )
 340                CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-Unit',
 341      $                     N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
 342      $                     WORK( N2+1 ), LDWORK )
 343 *
 344 *              Multiply bottom part of C by Q22**T.
 345 *
 346                CALL DGEMM( 'Transpose', 'No Transpose', N1, LEN, N2,
 347      $                     ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
 348      $                     ONE, WORK( N2+1 ), LDWORK )
 349 *
 350 *              Copy everything back.
 351 *
 352                CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
 353      $                      LDC )
 354             END DO
 355          END IF
 356       ELSE
 357          IF( NOTRAN ) THEN
 358             DO I = 1, M, NB
 359                LEN = MIN( NB, M-I+1 )
 360                LDWORK = LEN
 361 *
 362 *              Multiply right part of C by Q21.
 363 *
 364                CALL DLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
 365      $                      LDWORK )
 366                CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
 367      $                     LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
 368      $                     LDWORK )
 369 *
 370 *              Multiply left part of C by Q11.
 371 *
 372                CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
 373      $                     ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
 374      $                     LDWORK )
 375 *
 376 *              Multiply left part of C by Q12.
 377 *
 378                CALL DLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
 379      $                      WORK( 1 + N2*LDWORK ), LDWORK )
 380                CALL DTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
 381      $                     LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
 382      $                     WORK( 1 + N2*LDWORK ), LDWORK )
 383 *
 384 *              Multiply right part of C by Q22.
 385 *
 386                CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
 387      $                     ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
 388      $                     ONE, WORK( 1 + N2*LDWORK ), LDWORK )
 389 *
 390 *              Copy everything back.
 391 *
 392                CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
 393      $                      LDC )
 394             END DO
 395          ELSE
 396             DO I = 1, M, NB
 397                LEN = MIN( NB, M-I+1 )
 398                LDWORK = LEN
 399 *
 400 *              Multiply right part of C by Q12**T.
 401 *
 402                CALL DLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
 403      $                      LDWORK )
 404                CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Non-Unit',
 405      $                     LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
 406      $                     LDWORK )
 407 *
 408 *              Multiply left part of C by Q11**T.
 409 *
 410                CALL DGEMM( 'No Transpose', 'Transpose', LEN, N1, N2,
 411      $                     ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
 412      $                     LDWORK )
 413 *
 414 *              Multiply left part of C by Q21**T.
 415 *
 416                CALL DLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
 417      $                      WORK( 1 + N1*LDWORK ), LDWORK )
 418                CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-Unit',
 419      $                     LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
 420      $                     WORK( 1 + N1*LDWORK ), LDWORK )
 421 *
 422 *              Multiply right part of C by Q22**T.
 423 *
 424                CALL DGEMM( 'No Transpose', 'Transpose', LEN, N2, N1,
 425      $                     ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
 426      $                     ONE, WORK( 1 + N1*LDWORK ), LDWORK )
 427 *
 428 *              Copy everything back.
 429 *
 430                CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
 431      $                      LDC )
 432             END DO
 433          END IF
 434       END IF
 435 *
 436       WORK( 1 ) = DBLE( LWKOPT )
 437       RETURN
 438 *
 439 *     End of DORM22
 440 *
 441       END