SRC/sorbdb6.f

   1 *> \brief \b SORBDB6
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SORBDB6 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorbdb6.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorbdb6.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorbdb6.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SORBDB6( 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 *       REAL               Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
  30 *       ..
  31 *
  32 *
  33 *> \par Purpose:
  34 *> =============
  35 *>
  36 *>\verbatim
  37 *>
  38 *> SORBDB6 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 realOTHERcomputational
 152 *
 153 *  =====================================================================
 154       SUBROUTINE SORBDB6( 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       REAL               Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
 168 *     ..
 169 *
 170 *  =====================================================================
 171 *
 172 *     .. Parameters ..
 173       REAL               ALPHASQ, REALONE, REALZERO
 174       PARAMETER          ( ALPHASQ = 0.01E0, REALONE = 1.0E0,
 175      $                     REALZERO = 0.0E0 )
 176       REAL               NEGONE, ONE, ZERO
 177       PARAMETER          ( NEGONE = -1.0E0, ONE = 1.0E0, ZERO = 0.0E0 )
 178 *     ..
 179 *     .. Local Scalars ..
 180       INTEGER            I
 181       REAL               NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
 182 *     ..
 183 *     .. External Subroutines ..
 184       EXTERNAL           SGEMV, SLASSQ, XERBLA
 185 *     ..
 186 *     .. Intrinsic Function ..
 187       INTRINSIC          MAX
 188 *     ..
 189 *     .. Executable Statements ..
 190 *
 191 *     Test input arguments
 192 *
 193       INFO = 0
 194       IF( M1 .LT. 0 ) THEN
 195          INFO = -1
 196       ELSE IF( M2 .LT. 0 ) THEN
 197          INFO = -2
 198       ELSE IF( N .LT. 0 ) THEN
 199          INFO = -3
 200       ELSE IF( INCX1 .LT. 1 ) THEN
 201          INFO = -5
 202       ELSE IF( INCX2 .LT. 1 ) THEN
 203          INFO = -7
 204       ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
 205          INFO = -9
 206       ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
 207          INFO = -11
 208       ELSE IF( LWORK .LT. N ) THEN
 209          INFO = -13
 210       END IF
 211 *
 212       IF( INFO .NE. 0 ) THEN
 213          CALL XERBLA( 'SORBDB6', -INFO )
 214          RETURN
 215       END IF
 216 *
 217 *     First, project X onto the orthogonal complement of Q's column
 218 *     space
 219 *
 220       SCL1 = REALZERO
 221       SSQ1 = REALONE
 222       CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 223       SCL2 = REALZERO
 224       SSQ2 = REALONE
 225       CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
 226       NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 227 *
 228       IF( M1 .EQ. 0 ) THEN
 229          DO I = 1, N
 230             WORK(I) = ZERO
 231          END DO
 232       ELSE
 233          CALL SGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
 234      $               1 )
 235       END IF
 236 *
 237       CALL SGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
 238 *
 239       CALL SGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
 240      $            INCX1 )
 241       CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
 242      $            INCX2 )
 243 *
 244       SCL1 = REALZERO
 245       SSQ1 = REALONE
 246       CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 247       SCL2 = REALZERO
 248       SSQ2 = REALONE
 249       CALL SLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
 250       NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 251 *
 252 *     If projection is sufficiently large in norm, then stop.
 253 *     If projection is zero, then stop.
 254 *     Otherwise, project again.
 255 *
 256       IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
 257          RETURN
 258       END IF
 259 *
 260       IF( NORMSQ2 .EQ. ZERO ) THEN
 261          RETURN
 262       END IF
 263 *
 264       NORMSQ1 = NORMSQ2
 265 *
 266       DO I = 1, N
 267          WORK(I) = ZERO
 268       END DO
 269 *
 270       IF( M1 .EQ. 0 ) THEN
 271          DO I = 1, N
 272             WORK(I) = ZERO
 273          END DO
 274       ELSE
 275          CALL SGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
 276      $               1 )
 277       END IF
 278 *
 279       CALL SGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
 280 *
 281       CALL SGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
 282      $            INCX1 )
 283       CALL SGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
 284      $            INCX2 )
 285 *
 286       SCL1 = REALZERO
 287       SSQ1 = REALONE
 288       CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 289       SCL2 = REALZERO
 290       SSQ2 = REALONE
 291       CALL SLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
 292       NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
 293 *
 294 *     If second projection is sufficiently large in norm, then do
 295 *     nothing more. Alternatively, if it shrunk significantly, then
 296 *     truncate it to zero.
 297 *
 298       IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
 299          DO I = 1, M1
 300             X1(I) = ZERO
 301          END DO
 302          DO I = 1, M2
 303             X2(I) = ZERO
 304          END DO
 305       END IF
 306 *
 307       RETURN
 308 *
 309 *     End of SORBDB6
 310 *
 311       END
 312