subroutine sorbdb6	(	integer	M1,
		integer	M2,
		integer	N,
		real, dimension(*)	X1,
		integer	INCX1,
		real, dimension(*)	X2,
		integer	INCX2,
		real, dimension(ldq1,*)	Q1,
		integer	LDQ1,
		real, dimension(ldq2,*)	Q2,
		integer	LDQ2,
		real, dimension(*)	WORK,
		integer	LWORK,
		integer	INFO
	)

SORBDB6

Download SORBDB6 + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 SORBDB6 orthogonalizes the column vector
      X = [ X1 ]
          [ X2 ]
 with respect to the columns of
      Q = [ Q1 ] .
          [ Q2 ]
 The columns of Q must be orthonormal.

 If the projection is zero according to Kahan's "twice is enough"
 criterion, then the zero vector is returned.

Parameters

[in]	M1	M1 is INTEGER The dimension of X1 and the number of rows in Q1. 0 <= M1.
[in]	M2	M2 is INTEGER The dimension of X2 and the number of rows in Q2. 0 <= M2.
[in]	N	N is INTEGER The number of columns in Q1 and Q2. 0 <= N.
[in,out]	X1	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.
[in]	INCX1	INCX1 is INTEGER Increment for entries of X1.
[in,out]	X2	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.
[in]	INCX2	INCX2 is INTEGER Increment for entries of X2.
[in]	Q1	Q1 is REAL array, dimension (LDQ1, N) The top part of the orthonormal basis matrix.
[in]	LDQ1	LDQ1 is INTEGER The leading dimension of Q1. LDQ1 >= M1.
[in]	Q2	Q2 is REAL array, dimension (LDQ2, N) The bottom part of the orthonormal basis matrix.
[in]	LDQ2	LDQ2 is INTEGER The leading dimension of Q2. LDQ2 >= M2.
[out]	WORK	WORK is REAL array, dimension (LWORK)
[in]	LWORK	LWORK is INTEGER The dimension of the array WORK. LWORK >= N.
[out]	INFO	INFO is INTEGER = 0: successful exit. < 0: if INFO = -i, the i-th argument had an illegal value.

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Date

July 2012

Definition at line 156 of file sorbdb6.f.

*
*  -- LAPACK computational routine (version 3.5.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     July 2012
*
*     .. Scalar Arguments ..
      INTEGER            incx1, incx2, info, ldq1, ldq2, lwork, m1, m2,
     $                   n
*     ..
*     .. Array Arguments ..
      REAL               q1(ldq1,*), q2(ldq2,*), work(*), x1(*), x2(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               alphasq, realone, realzero
      parameter                ( alphasq = 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
      REAL               normsq1, normsq2, scl1, scl2, ssq1, ssq2
*     ..
*     .. 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
*
*     First, project X onto the orthogonal complement of Q's column
*     space
*
      scl1 = realzero
      ssq1 = realone
      CALL slassq( m1, x1, incx1, scl1, ssq1 )
      scl2 = realzero
      ssq2 = realone
      CALL slassq( m2, x2, incx2, scl2, ssq2 )
      normsq1 = scl1**2*ssq1 + scl2**2*ssq2
*
      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 )
*
      scl1 = realzero
      ssq1 = realone
      CALL slassq( m1, x1, incx1, scl1, ssq1 )
      scl2 = realzero
      ssq2 = realone
      CALL slassq( m2, x2, incx2, scl2, ssq2 )
      normsq2 = scl1**2*ssq1 + scl2**2*ssq2
*
*     If projection is sufficiently large in norm, then stop.
*     If projection is zero, then stop.
*     Otherwise, project again.
*
      IF( normsq2 .GE. alphasq*normsq1 ) THEN
         RETURN
      END IF
*
      IF( normsq2 .EQ. zero ) THEN
         RETURN
      END IF
*      
      normsq1 = normsq2
*
      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 )
*
      scl1 = realzero
      ssq1 = realone
      CALL slassq( m1, x1, incx1, scl1, ssq1 )
      scl2 = realzero
      ssq2 = realone
      CALL slassq( m1, x1, incx1, scl1, ssq1 )
      normsq2 = scl1**2*ssq1 + scl2**2*ssq2
*
*     If second projection is sufficiently large in norm, then do
*     nothing more. Alternatively, if it shrunk significantly, then
*     truncate it to zero.
*
      IF( normsq2 .LT. alphasq*normsq1 ) THEN
         DO i = 1, m1
            x1(i) = zero
         END DO
         DO i = 1, m2
            x2(i) = zero
         END DO
      END IF
*
      RETURN
*      
*     End of SORBDB6
*

Here is the call graph for this function:

Here is the caller graph for this function: