← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/slicot/sb04qr.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/slicot/sb04qr.f
 
1
      SUBROUTINE SB04QR( M, D, IPR, INFO )
 
2
C
 
3
C     RELEASE 4.0, WGS COPYRIGHT 2000.
 
4
C
 
5
C     PURPOSE
 
6
C
 
7
C     To solve a linear algebraic system of order M whose coefficient
 
8
C     matrix has zeros below the third subdiagonal and zero elements on
 
9
C     the third subdiagonal with even column indices. The matrix is 
 
10
C     stored compactly, row-wise. 
 
11
C
 
12
C     ARGUMENTS
 
13
C
 
14
C     Input/Output Parameters
 
15
C
 
16
C     M       (input) INTEGER
 
17
C             The order of the system.  M >= 0, M even.
 
18
C             Note that parameter M should have twice the value in the
 
19
C             original problem (see SLICOT Library routine SB04QU).
 
20
C
 
21
C     D       (input/output) DOUBLE PRECISION array, dimension 
 
22
C             (M*M/2+4*M)
 
23
C             On entry, the first M*M/2 + 3*M elements of this array
 
24
C             must contain the coefficient matrix, stored compactly,
 
25
C             row-wise, and the next M elements must contain the right
 
26
C             hand side of the linear system, as set by SLICOT Library
 
27
C             routine SB04QU.
 
28
C             On exit, the content of this array is updated, the last M
 
29
C             elements containing the solution with components 
 
30
C             interchanged (see IPR). 
 
31
C
 
32
C     IPR     (output) INTEGER array, dimension (2*M)
 
33
C             The leading M elements contain information about the
 
34
C             row interchanges performed for solving the system.
 
35
C             Specifically, the i-th component of the solution is
 
36
C             specified by IPR(i).
 
37
C
 
38
C     Error Indicator
 
39
C
 
40
C     INFO    INTEGER
 
41
C             = 0:  successful exit;
 
42
C             = 1:  if a singular matrix was encountered.
 
43
C
 
44
C     METHOD
 
45
C
 
46
C     Gaussian elimination with partial pivoting is used. The rows of
 
47
C     the matrix are not actually permuted, only their indices are
 
48
C     interchanged in array IPR.
 
49
C
 
50
C     REFERENCES
 
51
C
 
52
C     [1] Golub, G.H., Nash, S. and Van Loan, C.F.
 
53
C         A Hessenberg-Schur method for the problem AX + XB = C.
 
54
C         IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979.
 
55
C
 
56
C     [2] Sima, V.
 
57
C         Algorithms for Linear-quadratic Optimization.
 
58
C         Marcel Dekker, Inc., New York, 1996.
 
59
C
 
60
C     NUMERICAL ASPECTS
 
61
C     
 
62
C     None.
 
63
C
 
64
C     CONTRIBUTORS
 
65
C
 
66
C     D. Sima, University of Bucharest, May 2000.
 
67
C
 
68
C     REVISIONS
 
69
C
 
70
C     -
 
71
C
 
72
C     KEYWORDS
 
73
C
 
74
C     Hessenberg form, orthogonal transformation, real Schur form,
 
75
C     Sylvester equation.
 
76
C
 
77
C     ******************************************************************
 
78
C
 
79
      DOUBLE PRECISION  ZERO
 
80
      PARAMETER         ( ZERO = 0.0D0 )
 
81
C     .. Scalar Arguments ..
 
82
      INTEGER           INFO, M
 
83
C     .. Array Arguments ..
 
84
      INTEGER           IPR(*)
 
85
      DOUBLE PRECISION  D(*)
 
86
C     .. Local Scalars ..
 
87
      INTEGER           I, I1, I2, IPRM, IPRM1, J, K, L, M1, MPI, MPI1,
 
88
     $                  MPI2
 
89
      DOUBLE PRECISION  D1, D2, D3, DMAX
 
90
C     .. External Subroutines ..
 
91
      EXTERNAL          DAXPY
 
92
C     .. Intrinsic Functions ..
 
93
      INTRINSIC         ABS, MOD
 
94
C     .. Executable Statements ..
 
95
C
 
96
      INFO = 0
 
97
      I2   = M*M/2 + 3*M
 
98
      MPI  = M
 
99
      IPRM = I2
 
100
      M1 = M
 
101
      I1 = 1
 
102
C
 
103
      DO 20 I = 1, M
 
104
         MPI = MPI + 1
 
105
         IPRM = IPRM + 1
 
106
         IPR(MPI) = I1
 
107
         IPR(I) = IPRM
 
108
         I1 = I1 + M1
 
109
         IF ( I.GE.4 .AND. MOD( I, 2 ).EQ.0 ) M1 = M1 - 2 
 
110
   20 CONTINUE
 
111
C
 
112
      M1   = M - 1
 
113
      MPI1 = M + 1 
 
114
C
 
115
C     Reduce to upper triangular form.
 
116
C
 
117
      DO 80 I = 1, M1
 
118
         MPI  = MPI1
 
119
         MPI1 = MPI1 + 1
 
120
         IPRM = IPR(MPI)
 
121
         D1 = D(IPRM)
 
122
         I1 = 3
 
123
         IF ( MOD( I, 2 ).EQ.0 ) I1 = 2
 
124
         IF ( I.EQ.M1 ) I1 = 1
 
125
         MPI2 = MPI + I1
 
126
         L = 0
 
127
         DMAX = ABS( D1 )
 
128
C
 
129
         DO 40 J = MPI1, MPI2
 
130
            D2 = D(IPR(J))
 
131
            D3 = ABS( D2 )
 
132
            IF ( D3.GT.DMAX ) THEN
 
133
               DMAX = D3
 
134
               D1 = D2
 
135
               L = J - MPI
 
136
            END IF 
 
137
   40    CONTINUE
 
138
C
 
139
C        Check singularity.
 
140
C
 
141
         IF ( DMAX.EQ.ZERO ) THEN
 
142
            INFO = 1
 
143
            RETURN
 
144
         END IF
 
145
C
 
146
         IF ( L.GT.0 ) THEN
 
147
C
 
148
C           Permute the row indices.
 
149
C
 
150
            K = IPRM
 
151
            J = MPI + L
 
152
            IPRM   = IPR(J)
 
153
            IPR(J) = K
 
154
            IPR(MPI) = IPRM
 
155
            K  = IPR(I)
 
156
            I2 = I + L
 
157
            IPR(I)  = IPR(I2)
 
158
            IPR(I2) = K
 
159
         END IF
 
160
         IPRM = IPRM + 1
 
161
C 
 
162
C        Annihilate the subdiagonal elements of the matrix.
 
163
C
 
164
         I2 = I
 
165
         D3 = D(IPR(I))
 
166
C
 
167
         DO 60 J = MPI1, MPI2
 
168
            I2 = I2 + 1
 
169
            IPRM1 = IPR(J)
 
170
            DMAX = -D(IPRM1)/D1
 
171
            D(IPR(I2)) = D(IPR(I2)) + DMAX*D3
 
172
            CALL DAXPY( M-I, DMAX, D(IPRM), 1, D(IPRM1+1), 1 )
 
173
            IPR(J) = IPR(J) + 1
 
174
   60    CONTINUE
 
175
C
 
176
   80 CONTINUE
 
177
C
 
178
      MPI  = M + M
 
179
      IPRM = IPR(MPI)
 
180
C
 
181
C     Check singularity.
 
182
C
 
183
      IF ( D(IPRM).EQ.ZERO ) THEN
 
184
         INFO = 1
 
185
         RETURN
 
186
      END IF
 
187
C
 
188
C     Back substitution.
 
189
C
 
190
      D(IPR(M)) = D(IPR(M))/D(IPRM)
 
191
C
 
192
      DO 120 I = M1, 1, -1
 
193
         MPI   = MPI - 1
 
194
         IPRM  = IPR(MPI)
 
195
         IPRM1 = IPRM
 
196
         DMAX  = ZERO
 
197
C
 
198
         DO 100 K = I+1, M
 
199
            IPRM1 = IPRM1 + 1
 
200
            DMAX  = DMAX + D(IPR(K))*D(IPRM1)
 
201
  100    CONTINUE
 
202
C
 
203
         D(IPR(I)) = ( D(IPR(I)) - DMAX )/D(IPRM)
 
204
  120 CONTINUE
 
205
C
 
206
      RETURN
 
207
C *** Last line of SB04QR ***
 
208
      END