 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
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

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

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*
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition: slassq.f90:137
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
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: