sorbdb6.f

*> \brief \b SORBDB6
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SORBDB6 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorbdb6.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorbdb6.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorbdb6.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE SORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
*                           LDQ2, WORK, LWORK, INFO )
*
*       .. Scalar Arguments ..
*       INTEGER            INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
*      $                   N
*       ..
*       .. Array Arguments ..
*       REAL               Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*>\verbatim
*>
*> SORBDB6 orthogonalizes the column vector
*>      X = [ X1 ]
*>          [ X2 ]
*> with respect to the columns of
*>      Q = [ Q1 ] .
*>          [ Q2 ]
*> The Euclidean norm of X must be one and the columns of Q must be
*> orthonormal. The orthogonalized vector will be zero if and only if it
*> lies entirely in the range of Q.
*>
*> The projection is computed with at most two iterations of the
*> classical Gram-Schmidt algorithm, see
*> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
*>   analysis of the Gram-Schmidt algorithm with reorthogonalization."
*>   2002. CERFACS Technical Report No. TR/PA/02/33. URL:
*>   https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
*>
*>\endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M1
*> \verbatim
*>          M1 is INTEGER
*>           The dimension of X1 and the number of rows in Q1. 0 <= M1.
*> \endverbatim
*>
*> \param[in] M2
*> \verbatim
*>          M2 is INTEGER
*>           The dimension of X2 and the number of rows in Q2. 0 <= M2.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>           The number of columns in Q1 and Q2. 0 <= N.
*> \endverbatim
*>
*> \param[in,out] X1
*> \verbatim
*>          X1 is REAL array, dimension (M1)
*>           On entry, the top part of the vector to be orthogonalized.
*>           On exit, the top part of the projected vector.
*> \endverbatim
*>
*> \param[in] INCX1
*> \verbatim
*>          INCX1 is INTEGER
*>           Increment for entries of X1.
*> \endverbatim
*>
*> \param[in,out] X2
*> \verbatim
*>          X2 is REAL array, dimension (M2)
*>           On entry, the bottom part of the vector to be
*>           orthogonalized. On exit, the bottom part of the projected
*>           vector.
*> \endverbatim
*>
*> \param[in] INCX2
*> \verbatim
*>          INCX2 is INTEGER
*>           Increment for entries of X2.
*> \endverbatim
*>
*> \param[in] Q1
*> \verbatim
*>          Q1 is REAL array, dimension (LDQ1, N)
*>           The top part of the orthonormal basis matrix.
*> \endverbatim
*>
*> \param[in] LDQ1
*> \verbatim
*>          LDQ1 is INTEGER
*>           The leading dimension of Q1. LDQ1 >= M1.
*> \endverbatim
*>
*> \param[in] Q2
*> \verbatim
*>          Q2 is REAL array, dimension (LDQ2, N)
*>           The bottom part of the orthonormal basis matrix.
*> \endverbatim
*>
*> \param[in] LDQ2
*> \verbatim
*>          LDQ2 is INTEGER
*>           The leading dimension of Q2. LDQ2 >= M2.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>           The dimension of the array WORK. LWORK >= N.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          INFO is INTEGER
*>           = 0:  successful exit.
*>           < 0:  if INFO = -i, the i-th argument had an illegal value.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup unbdb6
*
*  =====================================================================
      SUBROUTINE sorbdb6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
     $                    LDQ2, WORK, LWORK, INFO )
*
*  -- LAPACK computational routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
     $                   n
*     ..
*     .. Array Arguments ..
      REAL               Q1(ldq1,*), Q2(ldq2,*), WORK(*), X1(*), X2(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ALPHA, REALONE, REALZERO
      parameter( alpha = 0.01e0, realone = 1.0e0,
     $                     realzero = 0.0e0 )
      REAL               NEGONE, ONE, ZERO
      parameter( negone = -1.0e0, one = 1.0e0, zero = 0.0e0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IX
      REAL               EPS, NORM, NORM_NEW, SCL, SSQ
*     ..
*     .. External Functions ..
      REAL               SLAMCH
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgemv, slassq, xerbla
*     ..
*     .. Intrinsic Function ..
      INTRINSIC          max
*     ..
*     .. Executable Statements ..
*
*     Test input arguments
*
      info = 0
      IF( m1 .LT. 0 ) THEN
         info = -1
      ELSE IF( m2 .LT. 0 ) THEN
         info = -2
      ELSE IF( n .LT. 0 ) THEN
         info = -3
      ELSE IF( incx1 .LT. 1 ) THEN
         info = -5
      ELSE IF( incx2 .LT. 1 ) THEN
         info = -7
      ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
         info = -9
      ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
         info = -11
      ELSE IF( lwork .LT. n ) THEN
         info = -13
      END IF
*
      IF( info .NE. 0 ) THEN
         CALL xerbla( 'SORBDB6', -info )
         RETURN
      END IF
*
      eps = slamch( 'Precision' )
*
*     First, project X onto the orthogonal complement of Q's column
*     space
*
*     Christoph Conrads: In debugging mode the norm should be computed
*     and an assertion added comparing the norm with one. Alas, Fortran
*     never made it into 1989 when assert() was introduced into the C
*     programming language.
      norm = realone
*
      IF( m1 .EQ. 0 ) THEN
         DO i = 1, n
            work(i) = zero
         END DO
      ELSE
         CALL sgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
     $               1 )
      END IF
*
      CALL sgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
*
      CALL sgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
     $            incx1 )
      CALL sgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
     $            incx2 )
*
      scl = realzero
      ssq = realzero
      CALL slassq( m1, x1, incx1, scl, ssq )
      CALL slassq( m2, x2, incx2, scl, ssq )
      norm_new = scl * sqrt(ssq)
*
*     If projection is sufficiently large in norm, then stop.
*     If projection is zero, then stop.
*     Otherwise, project again.
*
      IF( norm_new .GE. alpha * norm ) THEN
         RETURN
      END IF
*
      IF( norm_new .LE. n * eps * norm ) THEN
         DO ix = 1, 1 + (m1-1)*incx1, incx1
           x1( ix ) = zero
         END DO
         DO ix = 1, 1 + (m2-1)*incx2, incx2
           x2( ix ) = zero
         END DO
         RETURN
      END IF
*
      norm = norm_new
*
      DO i = 1, n
         work(i) = zero
      END DO
*
      IF( m1 .EQ. 0 ) THEN
         DO i = 1, n
            work(i) = zero
         END DO
      ELSE
         CALL sgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
     $               1 )
      END IF
*
      CALL sgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
*
      CALL sgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
     $            incx1 )
      CALL sgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
     $            incx2 )
*
      scl = realzero
      ssq = realzero
      CALL slassq( m1, x1, incx1, scl, ssq )
      CALL slassq( m2, x2, incx2, scl, ssq )
      norm_new = scl * sqrt(ssq)
*
*     If second projection is sufficiently large in norm, then do
*     nothing more. Alternatively, if it shrunk significantly, then
*     truncate it to zero.
*
      IF( norm_new .LT. alpha * norm ) THEN
         DO ix = 1, 1 + (m1-1)*incx1, incx1
            x1(ix) = zero
         END DO
         DO ix = 1, 1 + (m2-1)*incx2, incx2
            x2(ix) = zero
         END DO
      END IF
*
      RETURN
*
*     End of SORBDB6
*
      END