File:
[local] /
rpl /
lapack /
lapack /
dorbdb6.f
Revision
1.2:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:24 2014 UTC (10 years, 3 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
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: *> \date July 2012
150: *
151: *> \ingroup doubleOTHERcomputational
152: *
153: * =====================================================================
154: SUBROUTINE DORBDB6( 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: DOUBLE PRECISION 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: DOUBLE PRECISION NEGONE, ONE, ZERO
177: PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0, ZERO = 0.0D0 )
178: * ..
179: * .. Local Scalars ..
180: INTEGER I
181: DOUBLE PRECISION NORMSQ1, NORMSQ2, SCL1, SCL2, SSQ1, SSQ2
182: * ..
183: * .. External Subroutines ..
184: EXTERNAL DGEMV, DLASSQ, XERBLA
185: * ..
186: * .. Intrinsic Function ..
187: INTRINSIC MAX
188: * ..
189: * .. Executable Statements ..
190: *
191: * Test input arguments
192: *
193: INFO = 0
194: IF( M1 .LT. 0 ) THEN
195: INFO = -1
196: ELSE IF( M2 .LT. 0 ) THEN
197: INFO = -2
198: ELSE IF( N .LT. 0 ) THEN
199: INFO = -3
200: ELSE IF( INCX1 .LT. 1 ) THEN
201: INFO = -5
202: ELSE IF( INCX2 .LT. 1 ) THEN
203: INFO = -7
204: ELSE IF( LDQ1 .LT. MAX( 1, M1 ) ) THEN
205: INFO = -9
206: ELSE IF( LDQ2 .LT. MAX( 1, M2 ) ) THEN
207: INFO = -11
208: ELSE IF( LWORK .LT. N ) THEN
209: INFO = -13
210: END IF
211: *
212: IF( INFO .NE. 0 ) THEN
213: CALL XERBLA( 'DORBDB6', -INFO )
214: RETURN
215: END IF
216: *
217: * First, project X onto the orthogonal complement of Q's column
218: * space
219: *
220: SCL1 = REALZERO
221: SSQ1 = REALONE
222: CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
223: SCL2 = REALZERO
224: SSQ2 = REALONE
225: CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
226: NORMSQ1 = SCL1**2*SSQ1 + SCL2**2*SSQ2
227: *
228: IF( M1 .EQ. 0 ) THEN
229: DO I = 1, N
230: WORK(I) = ZERO
231: END DO
232: ELSE
233: CALL DGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
234: $ 1 )
235: END IF
236: *
237: CALL DGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
238: *
239: CALL DGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
240: $ INCX1 )
241: CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
242: $ INCX2 )
243: *
244: SCL1 = REALZERO
245: SSQ1 = REALONE
246: CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
247: SCL2 = REALZERO
248: SSQ2 = REALONE
249: CALL DLASSQ( M2, X2, INCX2, SCL2, SSQ2 )
250: NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
251: *
252: * If projection is sufficiently large in norm, then stop.
253: * If projection is zero, then stop.
254: * Otherwise, project again.
255: *
256: IF( NORMSQ2 .GE. ALPHASQ*NORMSQ1 ) THEN
257: RETURN
258: END IF
259: *
260: IF( NORMSQ2 .EQ. ZERO ) THEN
261: RETURN
262: END IF
263: *
264: NORMSQ1 = NORMSQ2
265: *
266: DO I = 1, N
267: WORK(I) = ZERO
268: END DO
269: *
270: IF( M1 .EQ. 0 ) THEN
271: DO I = 1, N
272: WORK(I) = ZERO
273: END DO
274: ELSE
275: CALL DGEMV( 'C', M1, N, ONE, Q1, LDQ1, X1, INCX1, ZERO, WORK,
276: $ 1 )
277: END IF
278: *
279: CALL DGEMV( 'C', M2, N, ONE, Q2, LDQ2, X2, INCX2, ONE, WORK, 1 )
280: *
281: CALL DGEMV( 'N', M1, N, NEGONE, Q1, LDQ1, WORK, 1, ONE, X1,
282: $ INCX1 )
283: CALL DGEMV( 'N', M2, N, NEGONE, Q2, LDQ2, WORK, 1, ONE, X2,
284: $ INCX2 )
285: *
286: SCL1 = REALZERO
287: SSQ1 = REALONE
288: CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
289: SCL2 = REALZERO
290: SSQ2 = REALONE
291: CALL DLASSQ( M1, X1, INCX1, SCL1, SSQ1 )
292: NORMSQ2 = SCL1**2*SSQ1 + SCL2**2*SSQ2
293: *
294: * If second projection is sufficiently large in norm, then do
295: * nothing more. Alternatively, if it shrunk significantly, then
296: * truncate it to zero.
297: *
298: IF( NORMSQ2 .LT. ALPHASQ*NORMSQ1 ) THEN
299: DO I = 1, M1
300: X1(I) = ZERO
301: END DO
302: DO I = 1, M2
303: X2(I) = ZERO
304: END DO
305: END IF
306: *
307: RETURN
308: *
309: * End of DORBDB6
310: *
311: END
312:
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