LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dorbdb6 ( integer  M1,
integer  M2,
integer  N,
double precision, dimension(*)  X1,
integer  INCX1,
double precision, dimension(*)  X2,
integer  INCX2,
double precision, dimension(ldq1,*)  Q1,
integer  LDQ1,
double precision, dimension(ldq2,*)  Q2,
integer  LDQ2,
double precision, dimension(*)  WORK,
integer  LWORK,
integer  INFO 
)

DORBDB6

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

Purpose: 

 DORBDB6 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dorbdb6.f.

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  DOUBLE PRECISION 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  DOUBLE PRECISION negone, one, zero
177  parameter ( negone = -1.0d0, one = 1.0d0, zero = 0.0d0 )
178 * ..
179 * .. Local Scalars ..
180  INTEGER i
181  DOUBLE PRECISION normsq1, normsq2, scl1, scl2, ssq1, ssq2
182 * ..
183 * .. External Subroutines ..
184  EXTERNAL dgemv, dlassq, 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( 'DORBDB6', -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 dlassq( m1, x1, incx1, scl1, ssq1 )
223  scl2 = realzero
224  ssq2 = realone
225  CALL dlassq( 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 dgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
234  $ 1 )
235  END IF
236 *
237  CALL dgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
238 *
239  CALL dgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
240  $ incx1 )
241  CALL dgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
242  $ incx2 )
243 *
244  scl1 = realzero
245  ssq1 = realone
246  CALL dlassq( m1, x1, incx1, scl1, ssq1 )
247  scl2 = realzero
248  ssq2 = realone
249  CALL dlassq( 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 dgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
276  $ 1 )
277  END IF
278 *
279  CALL dgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
280 *
281  CALL dgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
282  $ incx1 )
283  CALL dgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
284  $ incx2 )
285 *
286  scl1 = realzero
287  ssq1 = realone
288  CALL dlassq( m1, x1, incx1, scl1, ssq1 )
289  scl2 = realzero
290  ssq2 = realone
291  CALL dlassq( 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 DORBDB6
310 *
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dlassq
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105

Here is the call graph for this function:

Here is the caller graph for this function: