LAPACK 3.11.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 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
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 158 of file sorbdb6.f.

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