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 columns of Q must be orthonormal.
45: *>
46: *> If the projection is zero according to Kahan's "twice is enough"
47: *> criterion, then the zero vector is returned.
48: *>
49: *>\endverbatim
50: *
51: * Arguments:
52: * ==========
53: *
54: *> \param[in] M1
55: *> \verbatim
56: *> M1 is INTEGER
57: *> The dimension of X1 and the number of rows in Q1. 0 <= M1.
58: *> \endverbatim
59: *>
60: *> \param[in] M2
61: *> \verbatim
62: *> M2 is INTEGER
63: *> The dimension of X2 and the number of rows in Q2. 0 <= M2.
64: *> \endverbatim
65: *>
66: *> \param[in] N
67: *> \verbatim
68: *> N is INTEGER
69: *> The number of columns in Q1 and Q2. 0 <= N.
70: *> \endverbatim
71: *>
72: *> \param[in,out] X1
73: *> \verbatim
74: *> X1 is COMPLEX*16 array, dimension (M1)
75: *> On entry, the top part of the vector to be orthogonalized.
76: *> On exit, the top part of the projected vector.
77: *> \endverbatim
78: *>
79: *> \param[in] INCX1
80: *> \verbatim
81: *> INCX1 is INTEGER
82: *> Increment for entries of X1.
83: *> \endverbatim
84: *>
85: *> \param[in,out] X2
86: *> \verbatim
87: *> X2 is COMPLEX*16 array, dimension (M2)
88: *> On entry, the bottom part of the vector to be
89: *> orthogonalized. On exit, the bottom part of the projected
90: *> vector.
91: *> \endverbatim
92: *>
93: *> \param[in] INCX2
94: *> \verbatim
95: *> INCX2 is INTEGER
96: *> Increment for entries of X2.
97: *> \endverbatim
98: *>
99: *> \param[in] Q1
100: *> \verbatim
101: *> Q1 is COMPLEX*16 array, dimension (LDQ1, N)
102: *> The top part of the orthonormal basis matrix.
103: *> \endverbatim
104: *>
105: *> \param[in] LDQ1
106: *> \verbatim
107: *> LDQ1 is INTEGER
108: *> The leading dimension of Q1. LDQ1 >= M1.
109: *> \endverbatim
110: *>
111: *> \param[in] Q2
112: *> \verbatim
113: *> Q2 is COMPLEX*16 array, dimension (LDQ2, N)
114: *> The bottom part of the orthonormal basis matrix.
115: *> \endverbatim
116: *>
117: *> \param[in] LDQ2
118: *> \verbatim
119: *> LDQ2 is INTEGER
120: *> The leading dimension of Q2. LDQ2 >= M2.
121: *> \endverbatim
122: *>
123: *> \param[out] WORK
124: *> \verbatim
125: *> WORK is COMPLEX*16 array, dimension (LWORK)
126: *> \endverbatim
127: *>
128: *> \param[in] LWORK
129: *> \verbatim
130: *> LWORK is INTEGER
131: *> The dimension of the array WORK. LWORK >= N.
132: *> \endverbatim
133: *>
134: *> \param[out] INFO
135: *> \verbatim
136: *> INFO is INTEGER
137: *> = 0: successful exit.
138: *> < 0: if INFO = -i, the i-th argument had an illegal value.
139: *> \endverbatim
140: *
141: * Authors:
142: * ========
143: *
144: *> \author Univ. of Tennessee
145: *> \author Univ. of California Berkeley
146: *> \author Univ. of Colorado Denver
147: *> \author NAG Ltd.
148: *
149: *> \date July 2012
150: *
151: *> \ingroup complex16OTHERcomputational
152: *
153: * =====================================================================
154: SUBROUTINE ZUNBDB6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
155: $ LDQ2, WORK, LWORK, INFO )
156: *
157: * -- LAPACK computational routine (version 3.5.0) --
158: * -- LAPACK is a software package provided by Univ. of Tennessee, --
159: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
160: * July 2012
161: *
162: * .. Scalar Arguments ..
163: INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
164: $ N
165: * ..
166: * .. Array Arguments ..
167: COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
168: * ..
169: *
170: * =====================================================================
171: *
172: * .. Parameters ..
173: DOUBLE PRECISION ALPHASQ, REALONE, REALZERO
174: PARAMETER ( ALPHASQ = 0.01D0, REALONE = 1.0D0,
175: $ REALZERO = 0.0D0 )
176: COMPLEX*16 NEGONE, ONE, ZERO
177: PARAMETER ( NEGONE = (-1.0D0,0.0D0), ONE = (1.0D0,0.0D0),
178: $ ZERO = (0.0D0,0.0D0) )
179: * ..
180: * .. Local Scalars ..
181: INTEGER I
182: DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
183: * ..
184: * .. External Subroutines ..
185: EXTERNAL ZGEMV, ZLASSQ, XERBLA
186: * ..
187: * .. Intrinsic Function ..
188: INTRINSIC MAX
189: * ..
190: * .. Executable Statements ..
191: *
192: * Test input arguments
193: *
194: INFO = 0
195: IF( M1 .LT. 0 ) THEN
196: INFO = -1
197: ELSE IF( M2 .LT. 0 ) THEN
198: INFO = -2
199: ELSE IF( N .LT. 0 ) THEN
200: INFO = -3
201: ELSE IF( INCX1 .LT. 1 ) THEN
202: INFO = -5
203: ELSE IF( INCX2 .LT. 1 ) THEN
204: INFO = -7
205: ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
206: INFO = -9
207: ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
208: INFO = -11
209: ELSE IF( LWORK .LT. N ) THEN
210: INFO = -13
211: END IF
212: *
213: IF( INFO .NE. 0 ) THEN
214: CALL XERBLA( 'ZUNBDB6', -INFO )
215: RETURN
216: END IF
217: *
218: * First, project X onto the orthogonal complement of Q's column
219: * space
220: *
221: SCL1 = REALZERO
222: SSQ1 = REALONE
223: CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
224: SCL2 = REALZERO
225: SSQ2 = REALONE
226: CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
227: NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
228: *
229: IF( M1 .EQ. 0 ) THEN
230: DO I = 1, N
231: WORK(I) = ZERO
232: END DO
233: ELSE
234: CALL ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
235: $ 1 )
236: END IF
237: *
238: CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
239: *
240: CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
241: $ INCX1 )
242: CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
243: $ INCX2 )
244: *
245: SCL1 = REALZERO
246: SSQ1 = REALONE
247: CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
248: SCL2 = REALZERO
249: SSQ2 = REALONE
250: CALL ZLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
251: NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
252: *
253: * If projection is sufficiently large in norm, then stop.
254: * If projection is zero, then stop.
255: * Otherwise, project again.
256: *
257: IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
258: RETURN
259: END IF
260: *
261: IF( NORMSQ2 .EQ. ZERO ) THEN
262: RETURN
263: END IF
264: *
265: NORMSQ1 = NORMSQ2
266: *
267: DO I = 1, N
268: WORK(I) = ZERO
269: END DO
270: *
271: IF( M1 .EQ. 0 ) THEN
272: DO I = 1, N
273: WORK(I) = ZERO
274: END DO
275: ELSE
276: CALL ZGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
277: $ 1 )
278: END IF
279: *
280: CALL ZGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
281: *
282: CALL ZGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
283: $ INCX1 )
284: CALL ZGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
285: $ INCX2 )
286: *
287: SCL1 = REALZERO
288: SSQ1 = REALONE
289: CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
290: SCL2 = REALZERO
291: SSQ2 = REALONE
292: CALL ZLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
293: NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
294: *
295: * If second projection is sufficiently large in norm, then do
296: * nothing more. Alternatively, if it shrunk significantly, then
297: * truncate it to zero.
298: *
299: IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
300: DO I = 1, M1
301: X1(I) = ZERO
302: END DO
303: DO I = 1, M2
304: X2(I) = ZERO
305: END DO
306: END IF
307: *
308: RETURN
309: *
310: * End of ZUNBDB6
311: *
312: END
313:
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