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

◆ sorbdb6()

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

Definition at line 157 of file sorbdb6.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 REAL Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
170* ..
171*
172* =====================================================================
173*
174* .. Parameters ..
175 REAL ALPHA, REALONE, REALZERO
176 parameter( alpha = 0.83e0, realone = 1.0e0,
177 $ realzero = 0.0e0 )
178 REAL NEGONE, ONE, ZERO
179 parameter( negone = -1.0e0, one = 1.0e0, zero = 0.0e0 )
180* ..
181* .. Local Scalars ..
182 INTEGER I, IX
183 REAL EPS, NORM, NORM_NEW, SCL, SSQ
184* ..
185* .. External Functions ..
186 REAL SLAMCH
187* ..
188* .. External Subroutines ..
189 EXTERNAL sgemv, slassq, 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( 'SORBDB6', -info )
219 RETURN
220 END IF
221*
222 eps = slamch( 'Precision' )
223*
224* Compute the Euclidean norm of X
225*
226 scl = realzero
227 ssq = realzero
228 CALL slassq( m1, x1, incx1, scl, ssq )
229 CALL slassq( 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 sgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
241 $ 1 )
242 END IF
243*
244 CALL sgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
245*
246 CALL sgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
247 $ incx1 )
248 CALL sgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
249 $ incx2 )
250*
251 scl = realzero
252 ssq = realzero
253 CALL slassq( m1, x1, incx1, scl, ssq )
254 CALL slassq( 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 sgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
287 $ 1 )
288 END IF
289*
290 CALL sgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
291*
292 CALL sgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
293 $ incx1 )
294 CALL sgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
295 $ incx2 )
296*
297 scl = realzero
298 ssq = realzero
299 CALL slassq( m1, x1, incx1, scl, ssq )
300 CALL slassq( 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 SORBDB6
319*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slassq
subroutine slassq(n, x, incx, scale, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:124
Here is the call graph for this function:
Here is the caller graph for this function: