← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/lapack/dlanv2.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/lapack/dlanv2.f
 
1
      SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )
 
2
*
 
3
*  -- LAPACK auxiliary routine (version 2.0) --
 
4
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
 
5
*     Courant Institute, Argonne National Lab, and Rice University
 
6
*     September 30, 1994 
 
7
*
 
8
*     .. Scalar Arguments ..
 
9
      DOUBLE PRECISION   A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
 
10
*     ..
 
11
*
 
12
*  Purpose
 
13
*  =======
 
14
*
 
15
*  DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric
 
16
*  matrix in standard form:
 
17
*
 
18
*       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ]
 
19
*       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ]
 
20
*
 
21
*  where either
 
22
*  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or
 
23
*  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex
 
24
*  conjugate eigenvalues.
 
25
*
 
26
*  Arguments
 
27
*  =========
 
28
*
 
29
*  A       (input/output) DOUBLE PRECISION
 
30
*  B       (input/output) DOUBLE PRECISION
 
31
*  C       (input/output) DOUBLE PRECISION
 
32
*  D       (input/output) DOUBLE PRECISION
 
33
*          On entry, the elements of the input matrix.
 
34
*          On exit, they are overwritten by the elements of the
 
35
*          standardised Schur form.
 
36
*
 
37
*  RT1R    (output) DOUBLE PRECISION
 
38
*  RT1I    (output) DOUBLE PRECISION
 
39
*  RT2R    (output) DOUBLE PRECISION
 
40
*  RT2I    (output) DOUBLE PRECISION
 
41
*          The real and imaginary parts of the eigenvalues. If the
 
42
*          eigenvalues are both real, abs(RT1R) >= abs(RT2R); if the
 
43
*          eigenvalues are a complex conjugate pair, RT1I > 0.
 
44
*
 
45
*  CS      (output) DOUBLE PRECISION
 
46
*  SN      (output) DOUBLE PRECISION
 
47
*          Parameters of the rotation matrix.
 
48
*
 
49
*  =====================================================================
 
50
*
 
51
*     .. Parameters ..
 
52
      DOUBLE PRECISION   ZERO, HALF, ONE
 
53
      PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
 
54
*     ..
 
55
*     .. Local Scalars ..
 
56
      DOUBLE PRECISION   AA, BB, CC, CS1, DD, P, SAB, SAC, SIGMA, SN1,
 
57
     $                   TAU, TEMP
 
58
*     ..
 
59
*     .. External Functions ..
 
60
      DOUBLE PRECISION   DLAPY2
 
61
      EXTERNAL           DLAPY2
 
62
*     ..
 
63
*     .. Intrinsic Functions ..
 
64
      INTRINSIC          ABS, SIGN, SQRT
 
65
*     ..
 
66
*     .. Executable Statements ..
 
67
*
 
68
*     Initialize CS and SN
 
69
*
 
70
      CS = ONE
 
71
      SN = ZERO
 
72
*
 
73
      IF( C.EQ.ZERO ) THEN
 
74
         GO TO 10
 
75
*
 
76
      ELSE IF( B.EQ.ZERO ) THEN
 
77
*
 
78
*        Swap rows and columns
 
79
*
 
80
         CS = ZERO
 
81
         SN = ONE
 
82
         TEMP = D
 
83
         D = A
 
84
         A = TEMP
 
85
         B = -C
 
86
         C = ZERO
 
87
         GO TO 10
 
88
      ELSE IF( (A-D).EQ.ZERO .AND. SIGN( ONE, B ).NE.
 
89
     $   SIGN( ONE, C ) ) THEN
 
90
         GO TO 10
 
91
      ELSE
 
92
*
 
93
*        Make diagonal elements equal
 
94
*
 
95
         TEMP = A - D
 
96
         P = HALF*TEMP
 
97
         SIGMA = B + C
 
98
         TAU = DLAPY2( SIGMA, TEMP )
 
99
         CS1 = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) )
 
100
         SN1 = -( P / ( TAU*CS1 ) )*SIGN( ONE, SIGMA )
 
101
*
 
102
*        Compute [ AA  BB ] = [ A  B ] [ CS1 -SN1 ]
 
103
*                [ CC  DD ]   [ C  D ] [ SN1  CS1 ]
 
104
*
 
105
         AA = A*CS1 + B*SN1
 
106
         BB = -A*SN1 + B*CS1
 
107
         CC = C*CS1 + D*SN1
 
108
         DD = -C*SN1 + D*CS1
 
109
*
 
110
*        Compute [ A  B ] = [ CS1  SN1 ] [ AA  BB ]
 
111
*                [ C  D ]   [-SN1  CS1 ] [ CC  DD ]
 
112
*
 
113
         A = AA*CS1 + CC*SN1
 
114
         B = BB*CS1 + DD*SN1
 
115
         C = -AA*SN1 + CC*CS1
 
116
         D = -BB*SN1 + DD*CS1
 
117
*
 
118
*        Accumulate transformation
 
119
*
 
120
         TEMP = CS*CS1 - SN*SN1
 
121
         SN = CS*SN1 + SN*CS1
 
122
         CS = TEMP
 
123
*
 
124
         TEMP = HALF*( A+D )
 
125
         A = TEMP
 
126
         D = TEMP
 
127
*
 
128
         IF( C.NE.ZERO ) THEN
 
129
            IF( B.NE.ZERO ) THEN
 
130
               IF( SIGN( ONE, B ).EQ.SIGN( ONE, C ) ) THEN
 
131
*
 
132
*                 Real eigenvalues: reduce to upper triangular form
 
133
*
 
134
                  SAB = SQRT( ABS( B ) )
 
135
                  SAC = SQRT( ABS( C ) )
 
136
                  P = SIGN( SAB*SAC, C )
 
137
                  TAU = ONE / SQRT( ABS( B+C ) )
 
138
                  A = TEMP + P
 
139
                  D = TEMP - P
 
140
                  B = B - C
 
141
                  C = ZERO
 
142
                  CS1 = SAB*TAU
 
143
                  SN1 = SAC*TAU
 
144
                  TEMP = CS*CS1 - SN*SN1
 
145
                  SN = CS*SN1 + SN*CS1
 
146
                  CS = TEMP
 
147
               END IF
 
148
            ELSE
 
149
               B = -C
 
150
               C = ZERO
 
151
               TEMP = CS
 
152
               CS = -SN
 
153
               SN = TEMP
 
154
            END IF
 
155
         END IF
 
156
      END IF
 
157
*
 
158
   10 CONTINUE
 
159
*
 
160
*     Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
 
161
*
 
162
      RT1R = A
 
163
      RT2R = D
 
164
      IF( C.EQ.ZERO ) THEN
 
165
         RT1I = ZERO
 
166
         RT2I = ZERO
 
167
      ELSE
 
168
         RT1I = SQRT( ABS( B ) )*SQRT( ABS( C ) )
 
169
         RT2I = -RT1I
 
170
      END IF
 
171
      RETURN
 
172
*
 
173
*     End of DLANV2
 
174
*
 
175
      END