Annotation of rpl/lapack/lapack/dlags2.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAGS2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAGS2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlags2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlags2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlags2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
! 22: * SNV, CSQ, SNQ )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * LOGICAL UPPER
! 26: * DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
! 27: * $ SNU, SNV
! 28: * ..
! 29: *
! 30: *
! 31: *> \par Purpose:
! 32: * =============
! 33: *>
! 34: *> \verbatim
! 35: *>
! 36: *> DLAGS2 computes 2-by-2 orthogonal matrices U, V and Q, such
! 37: *> that if ( UPPER ) then
! 38: *>
! 39: *> U**T *A*Q = U**T *( A1 A2 )*Q = ( x 0 )
! 40: *> ( 0 A3 ) ( x x )
! 41: *> and
! 42: *> V**T*B*Q = V**T *( B1 B2 )*Q = ( x 0 )
! 43: *> ( 0 B3 ) ( x x )
! 44: *>
! 45: *> or if ( .NOT.UPPER ) then
! 46: *>
! 47: *> U**T *A*Q = U**T *( A1 0 )*Q = ( x x )
! 48: *> ( A2 A3 ) ( 0 x )
! 49: *> and
! 50: *> V**T*B*Q = V**T*( B1 0 )*Q = ( x x )
! 51: *> ( B2 B3 ) ( 0 x )
! 52: *>
! 53: *> The rows of the transformed A and B are parallel, where
! 54: *>
! 55: *> U = ( CSU SNU ), V = ( CSV SNV ), Q = ( CSQ SNQ )
! 56: *> ( -SNU CSU ) ( -SNV CSV ) ( -SNQ CSQ )
! 57: *>
! 58: *> Z**T denotes the transpose of Z.
! 59: *>
! 60: *> \endverbatim
! 61: *
! 62: * Arguments:
! 63: * ==========
! 64: *
! 65: *> \param[in] UPPER
! 66: *> \verbatim
! 67: *> UPPER is LOGICAL
! 68: *> = .TRUE.: the input matrices A and B are upper triangular.
! 69: *> = .FALSE.: the input matrices A and B are lower triangular.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] A1
! 73: *> \verbatim
! 74: *> A1 is DOUBLE PRECISION
! 75: *> \endverbatim
! 76: *>
! 77: *> \param[in] A2
! 78: *> \verbatim
! 79: *> A2 is DOUBLE PRECISION
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] A3
! 83: *> \verbatim
! 84: *> A3 is DOUBLE PRECISION
! 85: *> On entry, A1, A2 and A3 are elements of the input 2-by-2
! 86: *> upper (lower) triangular matrix A.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] B1
! 90: *> \verbatim
! 91: *> B1 is DOUBLE PRECISION
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[in] B2
! 95: *> \verbatim
! 96: *> B2 is DOUBLE PRECISION
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] B3
! 100: *> \verbatim
! 101: *> B3 is DOUBLE PRECISION
! 102: *> On entry, B1, B2 and B3 are elements of the input 2-by-2
! 103: *> upper (lower) triangular matrix B.
! 104: *> \endverbatim
! 105: *>
! 106: *> \param[out] CSU
! 107: *> \verbatim
! 108: *> CSU is DOUBLE PRECISION
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[out] SNU
! 112: *> \verbatim
! 113: *> SNU is DOUBLE PRECISION
! 114: *> The desired orthogonal matrix U.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[out] CSV
! 118: *> \verbatim
! 119: *> CSV is DOUBLE PRECISION
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[out] SNV
! 123: *> \verbatim
! 124: *> SNV is DOUBLE PRECISION
! 125: *> The desired orthogonal matrix V.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[out] CSQ
! 129: *> \verbatim
! 130: *> CSQ is DOUBLE PRECISION
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[out] SNQ
! 134: *> \verbatim
! 135: *> SNQ is DOUBLE PRECISION
! 136: *> The desired orthogonal matrix Q.
! 137: *> \endverbatim
! 138: *
! 139: * Authors:
! 140: * ========
! 141: *
! 142: *> \author Univ. of Tennessee
! 143: *> \author Univ. of California Berkeley
! 144: *> \author Univ. of Colorado Denver
! 145: *> \author NAG Ltd.
! 146: *
! 147: *> \date November 2011
! 148: *
! 149: *> \ingroup doubleOTHERauxiliary
! 150: *
! 151: * =====================================================================
1.1 bertrand 152: SUBROUTINE DLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
153: $ SNV, CSQ, SNQ )
154: *
1.9 ! bertrand 155: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 156: * -- LAPACK is a software package provided by Univ. of Tennessee, --
157: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 158: * November 2011
1.1 bertrand 159: *
160: * .. Scalar Arguments ..
161: LOGICAL UPPER
162: DOUBLE PRECISION A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, SNQ,
163: $ SNU, SNV
164: * ..
165: *
166: * =====================================================================
167: *
168: * .. Parameters ..
169: DOUBLE PRECISION ZERO
170: PARAMETER ( ZERO = 0.0D+0 )
171: * ..
172: * .. Local Scalars ..
173: DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
174: $ AVB21, AVB22, B, C, CSL, CSR, D, R, S1, S2,
175: $ SNL, SNR, UA11, UA11R, UA12, UA21, UA22, UA22R,
176: $ VB11, VB11R, VB12, VB21, VB22, VB22R
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL DLARTG, DLASV2
180: * ..
181: * .. Intrinsic Functions ..
182: INTRINSIC ABS
183: * ..
184: * .. Executable Statements ..
185: *
186: IF( UPPER ) THEN
187: *
188: * Input matrices A and B are upper triangular matrices
189: *
190: * Form matrix C = A*adj(B) = ( a b )
191: * ( 0 d )
192: *
193: A = A1*B3
194: D = A3*B1
195: B = A2*B1 - A1*B2
196: *
197: * The SVD of real 2-by-2 triangular C
198: *
199: * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
200: * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
201: *
202: CALL DLASV2( A, B, D, S1, S2, SNR, CSR, SNL, CSL )
203: *
204: IF( ABS( CSL ).GE.ABS( SNL ) .OR. ABS( CSR ).GE.ABS( SNR ) )
205: $ THEN
206: *
1.8 bertrand 207: * Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
208: * and (1,2) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 209: *
210: UA11R = CSL*A1
211: UA12 = CSL*A2 + SNL*A3
212: *
213: VB11R = CSR*B1
214: VB12 = CSR*B2 + SNR*B3
215: *
216: AUA12 = ABS( CSL )*ABS( A2 ) + ABS( SNL )*ABS( A3 )
217: AVB12 = ABS( CSR )*ABS( B2 ) + ABS( SNR )*ABS( B3 )
218: *
1.8 bertrand 219: * zero (1,2) elements of U**T *A and V**T *B
1.1 bertrand 220: *
221: IF( ( ABS( UA11R )+ABS( UA12 ) ).NE.ZERO ) THEN
222: IF( AUA12 / ( ABS( UA11R )+ABS( UA12 ) ).LE.AVB12 /
223: $ ( ABS( VB11R )+ABS( VB12 ) ) ) THEN
224: CALL DLARTG( -UA11R, UA12, CSQ, SNQ, R )
225: ELSE
226: CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
227: END IF
228: ELSE
229: CALL DLARTG( -VB11R, VB12, CSQ, SNQ, R )
230: END IF
231: *
232: CSU = CSL
233: SNU = -SNL
234: CSV = CSR
235: SNV = -SNR
236: *
237: ELSE
238: *
1.8 bertrand 239: * Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
240: * and (2,2) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 241: *
242: UA21 = -SNL*A1
243: UA22 = -SNL*A2 + CSL*A3
244: *
245: VB21 = -SNR*B1
246: VB22 = -SNR*B2 + CSR*B3
247: *
248: AUA22 = ABS( SNL )*ABS( A2 ) + ABS( CSL )*ABS( A3 )
249: AVB22 = ABS( SNR )*ABS( B2 ) + ABS( CSR )*ABS( B3 )
250: *
1.8 bertrand 251: * zero (2,2) elements of U**T*A and V**T*B, and then swap.
1.1 bertrand 252: *
253: IF( ( ABS( UA21 )+ABS( UA22 ) ).NE.ZERO ) THEN
254: IF( AUA22 / ( ABS( UA21 )+ABS( UA22 ) ).LE.AVB22 /
255: $ ( ABS( VB21 )+ABS( VB22 ) ) ) THEN
256: CALL DLARTG( -UA21, UA22, CSQ, SNQ, R )
257: ELSE
258: CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
259: END IF
260: ELSE
261: CALL DLARTG( -VB21, VB22, CSQ, SNQ, R )
262: END IF
263: *
264: CSU = SNL
265: SNU = CSL
266: CSV = SNR
267: SNV = CSR
268: *
269: END IF
270: *
271: ELSE
272: *
273: * Input matrices A and B are lower triangular matrices
274: *
275: * Form matrix C = A*adj(B) = ( a 0 )
276: * ( c d )
277: *
278: A = A1*B3
279: D = A3*B1
280: C = A2*B3 - A3*B2
281: *
282: * The SVD of real 2-by-2 triangular C
283: *
284: * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
285: * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
286: *
287: CALL DLASV2( A, C, D, S1, S2, SNR, CSR, SNL, CSL )
288: *
289: IF( ABS( CSR ).GE.ABS( SNR ) .OR. ABS( CSL ).GE.ABS( SNL ) )
290: $ THEN
291: *
1.8 bertrand 292: * Compute the (2,1) and (2,2) elements of U**T *A and V**T *B,
293: * and (2,1) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 294: *
295: UA21 = -SNR*A1 + CSR*A2
296: UA22R = CSR*A3
297: *
298: VB21 = -SNL*B1 + CSL*B2
299: VB22R = CSL*B3
300: *
301: AUA21 = ABS( SNR )*ABS( A1 ) + ABS( CSR )*ABS( A2 )
302: AVB21 = ABS( SNL )*ABS( B1 ) + ABS( CSL )*ABS( B2 )
303: *
1.8 bertrand 304: * zero (2,1) elements of U**T *A and V**T *B.
1.1 bertrand 305: *
306: IF( ( ABS( UA21 )+ABS( UA22R ) ).NE.ZERO ) THEN
307: IF( AUA21 / ( ABS( UA21 )+ABS( UA22R ) ).LE.AVB21 /
308: $ ( ABS( VB21 )+ABS( VB22R ) ) ) THEN
309: CALL DLARTG( UA22R, UA21, CSQ, SNQ, R )
310: ELSE
311: CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
312: END IF
313: ELSE
314: CALL DLARTG( VB22R, VB21, CSQ, SNQ, R )
315: END IF
316: *
317: CSU = CSR
318: SNU = -SNR
319: CSV = CSL
320: SNV = -SNL
321: *
322: ELSE
323: *
1.8 bertrand 324: * Compute the (1,1) and (1,2) elements of U**T *A and V**T *B,
325: * and (1,1) element of |U|**T *|A| and |V|**T *|B|.
1.1 bertrand 326: *
327: UA11 = CSR*A1 + SNR*A2
328: UA12 = SNR*A3
329: *
330: VB11 = CSL*B1 + SNL*B2
331: VB12 = SNL*B3
332: *
333: AUA11 = ABS( CSR )*ABS( A1 ) + ABS( SNR )*ABS( A2 )
334: AVB11 = ABS( CSL )*ABS( B1 ) + ABS( SNL )*ABS( B2 )
335: *
1.8 bertrand 336: * zero (1,1) elements of U**T*A and V**T*B, and then swap.
1.1 bertrand 337: *
338: IF( ( ABS( UA11 )+ABS( UA12 ) ).NE.ZERO ) THEN
339: IF( AUA11 / ( ABS( UA11 )+ABS( UA12 ) ).LE.AVB11 /
340: $ ( ABS( VB11 )+ABS( VB12 ) ) ) THEN
341: CALL DLARTG( UA12, UA11, CSQ, SNQ, R )
342: ELSE
343: CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
344: END IF
345: ELSE
346: CALL DLARTG( VB12, VB11, CSQ, SNQ, R )
347: END IF
348: *
349: CSU = SNR
350: SNU = CSR
351: CSV = SNL
352: SNV = CSL
353: *
354: END IF
355: *
356: END IF
357: *
358: RETURN
359: *
360: * End of DLAGS2
361: *
362: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlags2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack