SRC/zunbdb6.f

   1 *> \brief \b ZUNBDB6
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZUNBDB6 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb6.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb6.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb6.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
  22 *                           LDQ2, WORK, LWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
  26 *      $                   N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       COMPLEX*16         Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *> =============
  35 *>
  36 *>\verbatim
  37 *>
  38 *> ZUNBDB6 orthogonalizes the column vector
  39 *>      X = [ X1 ]
  40 *>          [ X2 ]
  41 *> with respect to the columns of
  42 *>      Q = [ Q1 ] .
  43 *>          [ Q2 ]
  44 *> The columns of Q must be orthonormal.
  45 *>
  46 *> If the projection is zero according to Kahan's "twice is enough"
  47 *> criterion, then the zero vector is returned.
  48 *>
  49 *>\endverbatim
  50 *
  51 *  Arguments:
  52 *  ==========
  53 *
  54 *> \param[in] M1
  55 *> \verbatim
  56 *>          M1 is INTEGER
  57 *>           The dimension of X1 and the number of rows in Q1. 0 <= M1.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] M2
  61 *> \verbatim
  62 *>          M2 is INTEGER
  63 *>           The dimension of X2 and the number of rows in Q2. 0 <= M2.
  64 *> \endverbatim
  65 *>
  66 *> \param[in] N
  67 *> \verbatim
  68 *>          N is INTEGER
  69 *>           The number of columns in Q1 and Q2. 0 <= N.
  70 *> \endverbatim
  71 *>
  72 *> \param[in,out] X1
  73 *> \verbatim
  74 *>          X1 is COMPLEX*16 array, dimension (M1)
  75 *>           On entry, the top part of the vector to be orthogonalized.
  76 *>           On exit, the top part of the projected vector.
  77 *> \endverbatim
  78 *>
  79 *> \param[in] INCX1
  80 *> \verbatim
  81 *>          INCX1 is INTEGER
  82 *>           Increment for entries of X1.
  83 *> \endverbatim
  84 *>
  85 *> \param[in,out] X2
  86 *> \verbatim
  87 *>          X2 is COMPLEX*16 array, dimension (M2)
  88 *>           On entry, the bottom part of the vector to be
  89 *>           orthogonalized. On exit, the bottom part of the projected
  90 *>           vector.
  91 *> \endverbatim
  92 *>
  93 *> \param[in] INCX2
  94 *> \verbatim
  95 *>          INCX2 is INTEGER
  96 *>           Increment for entries of X2.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] Q1
 100 *> \verbatim
 101 *>          Q1 is COMPLEX*16 array, dimension (LDQ1, N)
 102 *>           The top part of the orthonormal basis matrix.
 103 *> \endverbatim
 104 *>
 105 *> \param[in] LDQ1
 106 *> \verbatim
 107 *>          LDQ1 is INTEGER
 108 *>           The leading dimension of Q1. LDQ1 >= M1.
 109 *> \endverbatim
 110 *>
 111 *> \param[in] Q2
 112 *> \verbatim
 113 *>          Q2 is COMPLEX*16 array, dimension (LDQ2, N)
 114 *>           The bottom part of the orthonormal basis matrix.
 115 *> \endverbatim
 116 *>
 117 *> \param[in] LDQ2
 118 *> \verbatim
 119 *>          LDQ2 is INTEGER
 120 *>           The leading dimension of Q2. LDQ2 >= M2.
 121 *> \endverbatim
 122 *>
 123 *> \param[out] WORK
 124 *> \verbatim
 125 *>          WORK is COMPLEX*16 array, dimension (LWORK)
 126 *> \endverbatim
 127 *>
 128 *> \param[in] LWORK
 129 *> \verbatim
 130 *>          LWORK is INTEGER
 131 *>           The dimension of the array WORK. LWORK >= N.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] INFO
 135 *> \verbatim
 136 *>          INFO is INTEGER
 137 *>           = 0:  successful exit.
 138 *>           < 0:  if INFO = -i, the i-th argument had an illegal value.
 139 *> \endverbatim
 140 *
 141 *  Authors:
 142 *  ========
 143 *
 144 *> \author Univ. of Tennessee
 145 *> \author Univ. of California Berkeley
 146 *> \author Univ. of Colorado Denver
 147 *> \author NAG Ltd.
 148 *
 149 *> \date July 2012
 150 *
 151 *> \ingroup complex16OTHERcomputational
 152 *
 153 *  =====================================================================
 154       SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
 155      $                    LDQ2, WORK, LWORK, INFO )
 156 *
 157 *  -- LAPACK computational routine (version 3.5.0) --
 158 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 159 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 160 *     July 2012
 161 *
 162 *     .. Scalar Arguments ..
 163       INTEGER            INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
 164      $                   N
 165 *     ..
 166 *     .. Array Arguments ..
 167       COMPLEX*16         Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
 168 *     ..
 169 *
 170 *  =====================================================================
 171 *
 172 *     .. Parameters ..
 173       DOUBLE PRECISION   ALPHASQ, REALONE, REALZERO
 174       PARAMETER          ( ALPHASQ = 0.01D0, REALONE = 1.0D0,
 175      $                     REALZERO = 0.0D0 )
 176       COMPLEX*16         NEGONE, ONE, ZERO
 177       PARAMETER          ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0),
 178      $                     ZERO = (0.0D0,0.0D0) )
 179 *     ..
 180 *     .. Local Scalars ..
 181       INTEGER            I
 182       DOUBLE PRECISION   NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
 183 *     ..
 184 *     .. External Subroutines ..
 185       EXTERNAL           ZGEMV, ZLASSQ, XERBLA
 186 *     ..
 187 *     .. Intrinsic Function ..
 188       INTRINSIC          MAX
 189 *     ..
 190 *     .. Executable Statements ..
 191 *
 192 *     Test input arguments
 193 *
 194       INFO = 0
 195       IF( M1 .LT. 0 ) THEN
 196          INFO = -1
 197       ELSE IF( M2 .LT. 0 ) THEN
 198          INFO = -2
 199       ELSE IF( N .LT. 0 ) THEN
 200          INFO = -3
 201       ELSE IF( INCX1 .LT. 1 ) THEN
 202          INFO = -5
 203       ELSE IF( INCX2 .LT. 1 ) THEN
 204          INFO = -7
 205       ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
 206          INFO = -9
 207       ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
 208          INFO = -11
 209       ELSE IF( LWORK .LT. N ) THEN
 210          INFO = -13
 211       END IF
 212 *
 213       IF( INFO .NE. 0 ) THEN
 214          CALL XERBLA( 'ZUNBDB6', -INFO )
 215          RETURN
 216       END IF
 217 *
 218 *     First, project X onto the orthogonal complement of Q's column
 219 *     space
 220 *
 221       SCL1 = REALZERO
 222       SSQ1 = REALONE
 223       CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 224       SCL2 = REALZERO
 225       SSQ2 = REALONE
 226       CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
 227       NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 228 *
 229       IF( M1 .EQ. 0 ) THEN
 230          DO I = 1, N
 231             WORK(I) = ZERO
 232          END DO
 233       ELSE
 234          CALL ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
 235      $               1 )
 236       END IF
 237 *
 238       CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
 239 *
 240       CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
 241      $            INCX1 )
 242       CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
 243      $            INCX2 )
 244 *
 245       SCL1 = REALZERO
 246       SSQ1 = REALONE
 247       CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 248       SCL2 = REALZERO
 249       SSQ2 = REALONE
 250       CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
 251       NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 252 *
 253 *     If projection is sufficiently large in norm, then stop.
 254 *     If projection is zero, then stop.
 255 *     Otherwise, project again.
 256 *
 257       IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
 258          RETURN
 259       END IF
 260 *
 261       IF( NORMSQ2 .EQ. ZERO ) THEN
 262          RETURN
 263       END IF
 264 *
 265       NORMSQ1 = NORMSQ2
 266 *
 267       DO I = 1, N
 268          WORK(I) = ZERO
 269       END DO
 270 *
 271       IF( M1 .EQ. 0 ) THEN
 272          DO I = 1, N
 273             WORK(I) = ZERO
 274          END DO
 275       ELSE
 276          CALL ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
 277      $               1 )
 278       END IF
 279 *
 280       CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
 281 *
 282       CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
 283      $            INCX1 )
 284       CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
 285      $            INCX2 )
 286 *
 287       SCL1 = REALZERO
 288       SSQ1 = REALONE
 289       CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 290       SCL2 = REALZERO
 291       SSQ2 = REALONE
 292       CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 293       NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 294 *
 295 *     If second projection is sufficiently large in norm, then do
 296 *     nothing more. Alternatively, if it shrunk significantly, then
 297 *     truncate it to zero.
 298 *
 299       IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
 300          DO I = 1, M1
 301             X1(I) = ZERO
 302          END DO
 303          DO I = 1, M2
 304             X2(I) = ZERO
 305          END DO
 306       END IF
 307 *
 308       RETURN
 309 *
 310 *     End of ZUNBDB6
 311 *
 312       END
 313