LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dorbdb6()

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. 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
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.

Definition at line 157 of file dorbdb6.f.

159*
160* -- LAPACK computational routine --
161* -- LAPACK is a software package provided by Univ. of Tennessee, --
162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
163*
164* .. Scalar Arguments ..
165 INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
166 $ N
167* ..
168* .. Array Arguments ..
169 DOUBLE PRECISION Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 DOUBLE PRECISION ALPHA, REALONE, REALZERO
176 parameter( alpha = 0.83d0, realone = 1.0d0,
177 $ realzero = 0.0d0 )
178 DOUBLE PRECISION NEGONE, ONE, ZERO
179 parameter( negone = -1.0d0, one = 1.0d0, zero = 0.0d0 )
180* ..
181* .. Local Scalars ..
182 INTEGER I, IX
183 DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
184* ..
185* .. External Functions ..
186 DOUBLE PRECISION DLAMCH
187* ..
188* .. External Subroutines ..
189 EXTERNAL dgemv, dlassq, xerbla
190* ..
191* .. Intrinsic Function ..
192 INTRINSIC max
193* ..
194* .. Executable Statements ..
195*
196* Test input arguments
197*
198 info = 0
199 IF( m1 .LT. 0 ) THEN
200 info = -1
201 ELSE IF( m2 .LT. 0 ) THEN
202 info = -2
203 ELSE IF( n .LT. 0 ) THEN
204 info = -3
205 ELSE IF( incx1 .LT. 1 ) THEN
206 info = -5
207 ELSE IF( incx2 .LT. 1 ) THEN
208 info = -7
209 ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
210 info = -9
211 ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
212 info = -11
213 ELSE IF( lwork .LT. n ) THEN
214 info = -13
215 END IF
216*
217 IF( info .NE. 0 ) THEN
218 CALL xerbla( 'DORBDB6', -info )
219 RETURN
220 END IF
221*
222 eps = dlamch( 'Precision' )
223*
224* Compute the Euclidean norm of X
225*
226 scl = realzero
227 ssq = realzero
228 CALL dlassq( m1, x1, incx1, scl, ssq )
229 CALL dlassq( m2, x2, incx2, scl, ssq )
230 norm = scl * sqrt( ssq )
231*
232* First, project X onto the orthogonal complement of Q's column
233* space
234*
235 IF( m1 .EQ. 0 ) THEN
236 DO i = 1, n
237 work(i) = zero
238 END DO
239 ELSE
240 CALL dgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
241 $ 1 )
242 END IF
243*
244 CALL dgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
245*
246 CALL dgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
247 $ incx1 )
248 CALL dgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
249 $ incx2 )
250*
251 scl = realzero
252 ssq = realzero
253 CALL dlassq( m1, x1, incx1, scl, ssq )
254 CALL dlassq( m2, x2, incx2, scl, ssq )
255 norm_new = scl * sqrt(ssq)
256*
257* If projection is sufficiently large in norm, then stop.
258* If projection is zero, then stop.
259* Otherwise, project again.
260*
261 IF( norm_new .GE. alpha * norm ) THEN
262 RETURN
263 END IF
264*
265 IF( norm_new .LE. n * eps * norm ) THEN
266 DO ix = 1, 1 + (m1-1)*incx1, incx1
267 x1( ix ) = zero
268 END DO
269 DO ix = 1, 1 + (m2-1)*incx2, incx2
270 x2( ix ) = zero
271 END DO
272 RETURN
273 END IF
274*
275 norm = norm_new
276*
277 DO i = 1, n
278 work(i) = zero
279 END DO
280*
281 IF( m1 .EQ. 0 ) THEN
282 DO i = 1, n
283 work(i) = zero
284 END DO
285 ELSE
286 CALL dgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
287 $ 1 )
288 END IF
289*
290 CALL dgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
291*
292 CALL dgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
293 $ incx1 )
294 CALL dgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
295 $ incx2 )
296*
297 scl = realzero
298 ssq = realzero
299 CALL dlassq( m1, x1, incx1, scl, ssq )
300 CALL dlassq( m2, x2, incx2, scl, ssq )
301 norm_new = scl * sqrt(ssq)
302*
303* If second projection is sufficiently large in norm, then do
304* nothing more. Alternatively, if it shrunk significantly, then
305* truncate it to zero.
306*
307 IF( norm_new .LT. alpha * norm ) THEN
308 DO ix = 1, 1 + (m1-1)*incx1, incx1
309 x1(ix) = zero
310 END DO
311 DO ix = 1, 1 + (m2-1)*incx2, incx2
312 x2(ix) = zero
313 END DO
314 END IF
315*
316 RETURN
317*
318* End of DORBDB6
319*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:124
Here is the call graph for this function:
Here is the caller graph for this function: