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