1: *> \brief \b ZUNBDB6
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZUNBDB6 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunbdb6.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunbdb6.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunbdb6.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZUNBDB6( 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: * COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *>\verbatim
37: *>
38: *> ZUNBDB6 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 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 complex16OTHERcomputational
156: *
157: * =====================================================================
158: SUBROUTINE ZUNBDB6( 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: COMPLEX*16 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: COMPLEX*16 NEGONE, ONE, ZERO
180: PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0),
181: $ ZERO = (0.0D0,0.0D0) )
182: * ..
183: * .. Local Scalars ..
184: INTEGER I, IX
185: DOUBLE PRECISION EPS, NORM, NORM_NEW, SCL, SSQ
186: * ..
187: * .. External Functions ..
188: DOUBLE PRECISION DLAMCH
189: * ..
190: * .. External Subroutines ..
191: EXTERNAL ZGEMV, ZLASSQ, XERBLA
192: * ..
193: * .. Intrinsic Function ..
194: INTRINSIC MAX
195: * ..
196: * .. Executable Statements ..
197: *
198: * Test input arguments
199: *
200: INFO = 0
201: IF( M1 .LT. 0 ) THEN
202: INFO = -1
203: ELSE IF( M2 .LT. 0 ) THEN
204: INFO = -2
205: ELSE IF( N .LT. 0 ) THEN
206: INFO = -3
207: ELSE IF( INCX1 .LT. 1 ) THEN
208: INFO = -5
209: ELSE IF( INCX2 .LT. 1 ) THEN
210: INFO = -7
211: ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
212: INFO = -9
213: ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
214: INFO = -11
215: ELSE IF( LWORK .LT. N ) THEN
216: INFO = -13
217: END IF
218: *
219: IF( INFO .NE. 0 ) THEN
220: CALL XERBLA( 'ZUNBDB6', -INFO )
221: RETURN
222: END IF
223: *
224: EPS = DLAMCH( 'Precision' )
225: *
226: * First, project X onto the orthogonal complement of Q's column
227: * space
228: *
229: * Christoph Conrads: In debugging mode the norm should be computed
230: * and an assertion added comparing the norm with one. Alas, Fortran
231: * never made it into 1989 when assert() was introduced into the C
232: * programming language.
233: NORM = REALONE
234: *
235: IF( M1 .EQ. 0 ) THEN
236: DO I = 1, N
237: WORK(I) = ZERO
238: END DO
239: ELSE
240: CALL ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
241: $ 1 )
242: END IF
243: *
244: CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
245: *
246: CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
247: $ INCX1 )
248: CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
249: $ INCX2 )
250: *
251: SCL = REALZERO
252: SSQ = REALZERO
253: CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ )
254: CALL ZLASSQ( 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 ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
287: $ 1 )
288: END IF
289: *
290: CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
291: *
292: CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
293: $ INCX1 )
294: CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
295: $ INCX2 )
296: *
297: SCL = REALZERO
298: SSQ = REALZERO
299: CALL ZLASSQ( M1, X1, INCX1, SCL, SSQ )
300: CALL ZLASSQ( 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 ZUNBDB6
319: *
320: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zunbdb6.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack