OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
zunbdb6.f
Go to the documentation of this file.
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*> \ingroup complex16OTHERcomputational
150*
151* =====================================================================
152 SUBROUTINE zunbdb6( M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2,
153 $ LDQ2, WORK, LWORK, INFO )
154*
155* -- LAPACK computational routine --
156* -- LAPACK is a software package provided by Univ. of Tennessee, --
157* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
158*
159* .. Scalar Arguments ..
160 INTEGER INCX1, INCX2, INFO, LDQ1, LDQ2, LWORK, M1, M2,
161 $ n
162* ..
163* .. Array Arguments ..
164 COMPLEX*16 Q1(LDQ1,*), Q2(LDQ2,*), WORK(*), X1(*), X2(*)
165* ..
166*
167* =====================================================================
168*
169* .. Parameters ..
170 DOUBLE PRECISION ALPHASQ, REALONE, REALZERO
171 parameter( alphasq = 0.01d0, realone = 1.0d0,
172 $ realzero = 0.0d0 )
173 COMPLEX*16 NEGONE, ONE, ZERO
174 parameter( negone = (-1.0d0,0.0d0), one = (1.0d0,0.0d0),
175 $ zero = (0.0d0,0.0d0) )
176* ..
177* .. Local Scalars ..
178 INTEGER I
179 DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
180* ..
181* .. External Subroutines ..
182 EXTERNAL zgemv, zlassq, xerbla
183* ..
184* .. Intrinsic Function ..
185 INTRINSIC max
186* ..
187* .. Executable Statements ..
188*
189* Test input arguments
190*
191 info = 0
192 IF( m1 .LT. 0 ) THEN
193 info = -1
194 ELSE IF( m2 .LT. 0 ) THEN
195 info = -2
196 ELSE IF( n .LT. 0 ) THEN
197 info = -3
198 ELSE IF( incx1 .LT. 1 ) THEN
199 info = -5
200 ELSE IF( incx2 .LT. 1 ) THEN
201 info = -7
202 ELSE IF( ldq1 .LT. max( 1, m1 ) ) THEN
203 info = -9
204 ELSE IF( ldq2 .LT. max( 1, m2 ) ) THEN
205 info = -11
206 ELSE IF( lwork .LT. n ) THEN
207 info = -13
208 END IF
209*
210 IF( info .NE. 0 ) THEN
211 CALL xerbla( 'ZUNBDB6', -info )
212 RETURN
213 END IF
214*
215* First, project X onto the orthogonal complement of Q's column
216* space
217*
218 scl1 = realzero
219 ssq1 = realone
220 CALL zlassq( m1, x1, incx1, scl1, ssq1 )
221 scl2 = realzero
222 ssq2 = realone
223 CALL zlassq( m2, x2, incx2, scl2, ssq2 )
224 normsq1 = scl1**2*ssq1 + scl2**2*ssq2
225*
226 IF( m1 .EQ. 0 ) THEN
227 DO i = 1, n
228 work(i) = zero
229 END DO
230 ELSE
231 CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
232 $ 1 )
233 END IF
234*
235 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
236*
237 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
238 $ incx1 )
239 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
240 $ incx2 )
241*
242 scl1 = realzero
243 ssq1 = realone
244 CALL zlassq( m1, x1, incx1, scl1, ssq1 )
245 scl2 = realzero
246 ssq2 = realone
247 CALL zlassq( m2, x2, incx2, scl2, ssq2 )
248 normsq2 = scl1**2*ssq1 + scl2**2*ssq2
249*
250* If projection is sufficiently large in norm, then stop.
251* If projection is zero, then stop.
252* Otherwise, project again.
253*
254 IF( normsq2 .GE. alphasq*normsq1 ) THEN
255 RETURN
256 END IF
257*
258 IF( normsq2 .EQ. zero ) THEN
259 RETURN
260 END IF
261*
262 normsq1 = normsq2
263*
264 DO i = 1, n
265 work(i) = zero
266 END DO
267*
268 IF( m1 .EQ. 0 ) THEN
269 DO i = 1, n
270 work(i) = zero
271 END DO
272 ELSE
273 CALL zgemv( 'C', m1, n, one, q1, ldq1, x1, incx1, zero, work,
274 $ 1 )
275 END IF
276*
277 CALL zgemv( 'C', m2, n, one, q2, ldq2, x2, incx2, one, work, 1 )
278*
279 CALL zgemv( 'N', m1, n, negone, q1, ldq1, work, 1, one, x1,
280 $ incx1 )
281 CALL zgemv( 'N', m2, n, negone, q2, ldq2, work, 1, one, x2,
282 $ incx2 )
283*
284 scl1 = realzero
285 ssq1 = realone
286 CALL zlassq( m1, x1, incx1, scl1, ssq1 )
287 scl2 = realzero
288 ssq2 = realone
289 CALL zlassq( m1, x1, incx1, scl1, ssq1 )
290 normsq2 = scl1**2*ssq1 + scl2**2*ssq2
291*
292* If second projection is sufficiently large in norm, then do
293* nothing more. Alternatively, if it shrunk significantly, then
294* truncate it to zero.
295*
296 IF( normsq2 .LT. alphasq*normsq1 ) THEN
297 DO i = 1, m1
298 x1(i) = zero
299 END DO
300 DO i = 1, m2
301 x2(i) = zero
302 END DO
303 END IF
304*
305 RETURN
306*
307* End of ZUNBDB6
308*
309 END
310
zlassq
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.f90:137
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zunbdb6
subroutine zunbdb6(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
ZUNBDB6
Definition zunbdb6.f:154
zgemv
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
Definition zgemv.f:158
max
#define max(a, b)
Definition macros.h:21