1: *> \brief \b DORBDB6
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DORBDB6 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorbdb6.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorbdb6.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorbdb6.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
22: * LDQ2, WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
26: * $ N
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *>\verbatim
37: *>
38: *> DORBDB6 orthogonalizes the column vector
39: *> X = [ X1 ]
40: *> [ X2 ]
41: *> with respect to the columns of
42: *> Q = [ Q1 ] .
43: *> [ Q2 ]
44: *> The Euclidean norm of X must be one and the columns of Q must be
45: *> orthonormal. The orthogonalized vector will be zero if and only if it
46: *> lies entirely in the range of Q.
47: *>
48: *> The projection is computed with at most two iterations of the
49: *> classical Gram-Schmidt algorithm, see
50: *> * L. Giraud, J. Langou, M. Rozložník. "On the round-off error
51: *> analysis of the Gram-Schmidt algorithm with reorthogonalization."
52: *> 2002. CERFACS Technical Report No. TR/PA/02/33. URL:
53: *> https://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf
54: *>
55: *>\endverbatim
56: *
57: * Arguments:
58: * ==========
59: *
60: *> \param[in] M1
61: *> \verbatim
62: *> M1 is INTEGER
63: *> The dimension of X1 and the number of rows in Q1. 0 <= M1.
64: *> \endverbatim
65: *>
66: *> \param[in] M2
67: *> \verbatim
68: *> M2 is INTEGER
69: *> The dimension of X2 and the number of rows in Q2. 0 <= M2.
70: *> \endverbatim
71: *>
72: *> \param[in] N
73: *> \verbatim
74: *> N is INTEGER
75: *> The number of columns in Q1 and Q2. 0 <= N.
76: *> \endverbatim
77: *>
78: *> \param[in,out] X1
79: *> \verbatim
80: *> X1 is DOUBLE PRECISION array, dimension (M1)
81: *> On entry, the top part of the vector to be orthogonalized.
82: *> On exit, the top part of the projected vector.
83: *> \endverbatim
84: *>
85: *> \param[in] INCX1
86: *> \verbatim
87: *> INCX1 is INTEGER
88: *> Increment for entries of X1.
89: *> \endverbatim
90: *>
91: *> \param[in,out] X2
92: *> \verbatim
93: *> X2 is DOUBLE PRECISION array, dimension (M2)
94: *> On entry, the bottom part of the vector to be
95: *> orthogonalized. On exit, the bottom part of the projected
96: *> vector.
97: *> \endverbatim
98: *>
99: *> \param[in] INCX2
100: *> \verbatim
101: *> INCX2 is INTEGER
102: *> Increment for entries of X2.
103: *> \endverbatim
104: *>
105: *> \param[in] Q1
106: *> \verbatim
107: *> Q1 is DOUBLE PRECISION array, dimension (LDQ1, N)
108: *> The top part of the orthonormal basis matrix.
109: *> \endverbatim
110: *>
111: *> \param[in] LDQ1
112: *> \verbatim
113: *> LDQ1 is INTEGER
114: *> The leading dimension of Q1. LDQ1 >= M1.
115: *> \endverbatim
116: *>
117: *> \param[in] Q2
118: *> \verbatim
119: *> Q2 is DOUBLE PRECISION array, dimension (LDQ2, N)
120: *> The bottom part of the orthonormal basis matrix.
121: *> \endverbatim
122: *>
123: *> \param[in] LDQ2
124: *> \verbatim
125: *> LDQ2 is INTEGER
126: *> The leading dimension of Q2. LDQ2 >= M2.
127: *> \endverbatim
128: *>
129: *> \param[out] WORK
130: *> \verbatim
131: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
132: *> \endverbatim
133: *>
134: *> \param[in] LWORK
135: *> \verbatim
136: *> LWORK is INTEGER
137: *> The dimension of the array WORK. LWORK >= N.
138: *> \endverbatim
139: *>
140: *> \param[out] INFO
141: *> \verbatim
142: *> INFO is INTEGER
143: *> = 0: successful exit.
144: *> < 0: if INFO = -i, the i-th argument had an illegal value.
145: *> \endverbatim
146: *
147: * Authors:
148: * ========
149: *
150: *> \author Univ. of Tennessee
151: *> \author Univ. of California Berkeley
152: *> \author Univ. of Colorado Denver
153: *> \author NAG Ltd.
154: *
155: *> \ingroup doubleOTHERcomputational
156: *
157: * =====================================================================
158: SUBROUTINE DORBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
159: $ LDQ2, WORK, LWORK, INFO )
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: DOUBLE PRECISION Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
171: * ..
172: *
173: * =====================================================================
174: *
175: * .. Parameters ..
176: DOUBLE PRECISION ALPHA, REALONE, REALZERO
177: PARAMETER ( ALPHA = 0.01D0, REALONE = 1.0D0,
178: $ REALZERO = 0.0D0 )
179: DOUBLE PRECISION NEGONE, ONE, ZERO
180: PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
181: * ..
182: * .. Local Scalars ..
183: INTEGER I, IX
184: DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
185: * ..
186: * .. External Functions ..
187: DOUBLE PRECISION DLAMCH
188: * ..
189: * .. External Subroutines ..
190: EXTERNAL DGEMV, DLASSQ, 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( 'DORBDB6', -INFO )
220: RETURN
221: END IF
222: *
223: EPS = DLAMCH( '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 DGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
240: $ 1 )
241: END IF
242: *
243: CALL DGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
244: *
245: CALL DGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
246: $ INCX1 )
247: CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
248: $ INCX2 )
249: *
250: SCL = REALZERO
251: SSQ = REALZERO
252: CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
253: CALL DLASSQ( 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 DGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
286: $ 1 )
287: END IF
288: *
289: CALL DGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
290: *
291: CALL DGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
292: $ INCX1 )
293: CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
294: $ INCX2 )
295: *
296: SCL = REALZERO
297: SSQ = REALZERO
298: CALL DLASSQ( M1, X1, INCX1, SCL, SSQ )
299: CALL DLASSQ( 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 DORBDB6
318: *
319: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dorbdb6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack