1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
|
SUBROUTINE MC03NY( NBLCKS, NRA, NCA, A, LDA, E, LDE, IMUK, INUK,
$ VEPS, LDVEPS, INFO )
C
C SLICOT RELEASE 5.0.
C
C Copyright (c) 2002-2009 NICONET e.V.
C
C This program is free software: you can redistribute it and/or
C modify it under the terms of the GNU General Public License as
C published by the Free Software Foundation, either version 2 of
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program. If not, see
C <http://www.gnu.org/licenses/>.
C
C PURPOSE
C
C To determine a minimal basis of the right nullspace of the
C subpencil s*E(eps)-A(eps) using the method given in [1] (see
C Eqs.(4.6.8), (4.6.9)).
C This pencil only contains Kronecker column indices, and it must be
C in staircase form as supplied by SLICOT Library Routine MB04VD.
C The basis vectors are represented by matrix V(s) having the form
C
C | V11(s) V12(s) V13(s) . . V1n(s) |
C | V22(s) V23(s) V2n(s) |
C | V33(s) . |
C V(s) = | . . |
C | . . |
C | . . |
C | Vnn(s) |
C
C where n is the number of full row rank blocks in matrix A(eps) and
C
C k j-i
C Vij(s) = Vij,0 + Vij,1*s +...+ Vij,k*s +...+ Vij,j-i*s . (1)
C
C In other words, Vij,k is the coefficient corresponding to degree k
C in the matrix polynomial Vij(s).
C Vij,k has dimensions mu(i)-by-(mu(j)-nu(j)).
C The coefficients Vij,k are stored in the matrix VEPS as follows
C (for the case n = 3):
C
C sizes m1-n1 m2-n2 m2-n2 m3-n3 m3-n3 m3-n3
C
C m1 { | V11,0 || V12,0 | V12,1 || V13,0 | V13,1 | V13,2 ||
C | || | || | | ||
C VEPS = m2 { | || V22,0 | || V23,0 | V23,1 | ||
C | || | || | | ||
C m3 { | || | || V33,0 | | ||
C
C where mi = mu(i), ni = nu(i).
C Matrix VEPS has dimensions nrv-by-ncv where
C nrv = Sum(i=1,...,n) mu(i)
C ncv = Sum(i=1,...,n) i*(mu(i)-nu(i))
C
C ==================================================================
C REMARK: This routine is intended to be called only from the SLICOT
C routine MC03ND.
C ==================================================================
C
C ARGUMENTS
C
C Input/Output Parameters
C
C NBLCKS (input) INTEGER
C Number of full row rank blocks in subpencil
C s*E(eps)-A(eps) that contains all Kronecker column indices
C of s*E-A. NBLCKS >= 0.
C
C NRA (input) INTEGER
C Number of rows of the subpencil s*E(eps)-A(eps) in s*E-A.
C NRA = nu(1) + nu(2) + ... + nu(NBLCKS). NRA >= 0.
C
C NCA (input) INTEGER
C Number of columns of the subpencil s*E(eps)-A(eps) in
C s*E-A.
C NCA = mu(1) + mu(2) + ... + mu(NBLCKS). NCA >= 0.
C
C A (input/output) DOUBLE PRECISION array, dimension (LDA,NCA)
C E (input/output) DOUBLE PRECISION array, dimension (LDE,NCA)
C On entry, the leading NRA-by-NCA part of these arrays must
C contain the matrices A and E, where s*E-A is the
C transformed pencil s*E0-A0 which is the pencil associated
C with P(s) as described in [1] Section 4.6. The pencil
C s*E-A is assumed to be in generalized Schur form.
C On exit, these arrays contain no useful information.
C
C LDA INTEGER
C The leading dimension of array A. LDA >= MAX(1,NRA).
C
C LDE INTEGER
C The leading dimension of array E. LDE >= MAX(1,NRA).
C
C IMUK (input) INTEGER array, dimension (NBLCKS)
C This array must contain the column dimensions mu(k) of the
C full column rank blocks in the subpencil s*E(eps)-A(eps)
C of s*E-A. The content of IMUK is modified by the routine
C but restored on exit.
C
C INUK (input) INTEGER array, dimension (NBLCKS)
C This array must contain the row dimensions nu(k) of the
C full row rank blocks in the subpencil s*E(eps)-A(eps) of
C s*E-A.
C
C VEPS (output) DOUBLE PRECISION array, dimension (LDVEPS,ncv)
C Let nrv = Sum(i=1,...,NBLCKS) mu(i) = NCA,
C ncv = Sum(i=1,...,NBLCKS) i*(mu(i)-nu(i)).
C The leading nrv-by-ncv part of this array contains the
C column vectors of a minimal polynomial basis for the right
C nullspace of the subpencil s*E(eps)-A(eps). (See [1]
C Section 4.6.4.) An upper bound for ncv is (NRA+1)*NCA.
C
C LDVEPS INTEGER
C The leading dimension of array VEPS.
C LDVEPS >= MAX(1,NCA).
C
C Error Indicator
C
C INFO INTEGER
C = 0: successful exit;
C < 0: if INFO = -i, the i-th argument had an illegal
C value;
C > 0: if INFO = k, the k-th diagonal block of A had not a
C full row rank.
C
C REFERENCES
C
C [1] Th.G.J. Beelen, New Algorithms for Computing the Kronecker
C structure of a Pencil with Applications to Systems and
C Control Theory.
C Ph.D.Thesis, Eindhoven University of Technology, 1987.
C
C NUMERICAL ASPECTS
C
C None.
C
C CONTRIBUTORS
C
C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Mar. 1997.
C Supersedes Release 2.0 routine MC03BY by Th.G.J. Beelen,
C A.J. Geurts, and G.J.H.H. van den Hurk.
C
C REVISIONS
C
C Dec. 1997.
C
C KEYWORDS
C
C Elementary polynomial operations, Kronecker form, polynomial
C matrix, polynomial operations, staircase form.
C
C ******************************************************************
C
C .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
C .. Scalar Arguments ..
INTEGER INFO, LDA, LDE, LDVEPS, NBLCKS, NCA, NRA
C .. Array Arguments ..
INTEGER IMUK(*), INUK(*)
DOUBLE PRECISION A(LDA,*), E(LDE,*), VEPS(LDVEPS,*)
C .. Local Scalars ..
INTEGER AC1, AC2, AR1, ARI, ARK, DIF, EC1, ER1, I, J, K,
$ MUI, NCV, NRV, NUI, SMUI, SMUI1, VC1, VC2, VR1,
$ VR2, WC1, WR1
C .. Local Arrays ..
DOUBLE PRECISION DUMMY(1)
C .. External Subroutines ..
EXTERNAL DCOPY, DGEMM, DLASET, DSCAL, DTRTRS, XERBLA
C .. Executable Statements ..
C
INFO = 0
IF( NBLCKS.LT.0 ) THEN
INFO = -1
ELSE IF( NRA.LT.0 ) THEN
INFO = -2
ELSE IF( NCA.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, NRA ) ) THEN
INFO = -5
ELSE IF( LDE.LT.MAX( 1, NRA ) ) THEN
INFO = -7
ELSE IF( LDVEPS.LT.MAX( 1, NCA ) ) THEN
INFO = -11
END IF
C
IF ( INFO.NE.0 ) THEN
C
C Error return.
C
CALL XERBLA( 'MC03NY', -INFO )
RETURN
END IF
C
C Quick return if possible.
C
IF ( NBLCKS.EQ.0 .OR. NRA.EQ.0 .OR. NCA.EQ.0 )
$ RETURN
C
C Computation of the nonzero parts of W1 and W2:
C
C | AH11 AH12 ... AH1n | | EH11 EH12 ... EH1n |
C | AH22 AH2n | | EH22 EH2n |
C W1 = | . . |, W2 = | . . |
C | . . | | . . |
C | AHnn | | EHnn |
C
C with AHij = -pinv(Aii) * Aij, EHij = pinv(Aii) * Eij and EHii = 0,
C AHij and EHij have dimensions mu(i)-by-mu(j), Aii = [ Oi | Ri ],
C and
C Ri is a regular nu(i)-by-nu(i) upper triangular matrix;
C Oi is a not necessarily square null matrix.
C Note that the first mu(i)-nu(i) rows in AHij and EHij are zero.
C For memory savings, the nonzero parts of W1 and W2 are constructed
C over A and E, respectively.
C
C (AR1,AC1) denotes the position of the first element of the
C submatrix Ri in matrix Aii.
C EC1 is the index of the first column of Ai,i+1/Ei,i+1.
C
EC1 = 1
AR1 = 1
C
DO 40 I = 1, NBLCKS - 1
NUI = INUK(I)
IF ( NUI.EQ.0 ) GO TO 60
MUI = IMUK(I)
EC1 = EC1 + MUI
AC1 = EC1 - NUI
CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', NUI,
$ NCA-EC1+1, A(AR1,AC1), LDA, E(AR1,EC1), LDE,
$ INFO )
IF ( INFO.GT.0 ) THEN
INFO = I
RETURN
END IF
C
DO 20 J = 1, NUI
CALL DSCAL( J, -ONE, A(AR1,AC1+J-1), 1 )
20 CONTINUE
C
CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', NUI,
$ NCA-EC1+1, A(AR1,AC1), LDA, A(AR1,EC1), LDA,
$ INFO )
AR1 = AR1 + NUI
40 CONTINUE
C
60 CONTINUE
C
C The contents of the array IMUK is changed for temporary use in
C this routine as follows:
C
C IMUK(i) = Sum(j=1,...,i) mu(j).
C
C On return, the original contents of IMUK is restored.
C In the same loop the actual number of columns of VEPS is computed.
C The number of rows of VEPS is NCA.
C
C NRV = Sum(i=1,...,NBLCKS) mu(i) = NCA,
C NCV = Sum(i=1,...,NBLCKS) i*(mu(i)-nu(i)).
C
SMUI = 0
NCV = 0
C
DO 80 I = 1, NBLCKS
MUI = IMUK(I)
SMUI = SMUI + MUI
IMUK(I) = SMUI
NCV = NCV + I*( MUI - INUK(I) )
80 CONTINUE
C
NRV = NCA
C
C Computation of the matrix VEPS.
C
C Initialisation of VEPS to zero.
C
CALL DLASET( 'Full', NRV, NCV, ZERO, ZERO, VEPS, LDVEPS )
C | I |
C Set Vii,0 = Kii in VEPS , i=1,...,NBLCKS, where Kii = |---|
C | O |
C and I is an identity matrix of size mu(i)-nu(i),
C O is a null matrix, dimensions nu(i)-by-(mu(i)-nu(i)).
C
C WR1 := Sum(j=1,...,i-1) mu(j) + 1
C is the index of the first row in Vii,0 in VEPS.
C WC1 := Sum(j=1,...,i-1) j*(mu(j)-nu(j)) + 1
C is the index of the first column in Vii,0 in VEPS.
C
DUMMY(1) = ONE
NUI = IMUK(1) - INUK(1)
CALL DCOPY( NUI, DUMMY, 0, VEPS, LDVEPS+1 )
WR1 = IMUK(1) + 1
WC1 = NUI + 1
C
DO 100 I = 2, NBLCKS
NUI = IMUK(I) - IMUK(I-1) - INUK(I)
CALL DCOPY( NUI, DUMMY, 0, VEPS(WR1,WC1), LDVEPS+1 )
WR1 = IMUK(I) + 1
WC1 = WC1 + I*NUI
100 CONTINUE
C
C Determination of the remaining nontrivial matrices in Vij,k
C block column by block column with decreasing block row index.
C
C The computation starts with the second block column since V11,0
C has already been determined.
C The coefficients Vij,k satisfy the recurrence relation:
C
C Vij,k = Sum(r=i+1,...,j-k) AHir*Vrj,k +
C + Sum(r=i+1,...,j-k+1) EHir*Vrj,k-1, i + k < j,
C
C = EHi,i+1 * Vi+1,j,k-1 i + k = j.
C
C This recurrence relation can be derived from [1], (4.6.8)
C and formula (1) in Section PURPOSE.
C
VC1 = IMUK(1) - INUK(1) + 1
ARI = 1
C
DO 180 J = 2, NBLCKS
DIF = IMUK(J) - IMUK(J-1) - INUK(J)
ARI = ARI + INUK(J-1)
ARK = ARI
C
C Computation of the matrices Vij,k where i + k < j.
C Each matrix Vij,k has dimension mu(i)-by-(mu(j) - nu(j)).
C
DO 160 K = 0, J - 2
C
C VC1, VC2 are the first and last column index of Vij,k.
C
VC2 = VC1 + DIF - 1
AC2 = IMUK(J-K)
AR1 = ARK
ARK = ARK - INUK(J-K-1)
C
DO 120 I = J - K - 1, 1, -1
C
C Compute the first part of Vij,k in decreasing order:
C Vij,k := Vij,k + Sum(r=i+1,..,j-k) AHir*Vrj,k.
C The non-zero parts of AHir are stored in
C A(AR1:AR1+nu(i)-1,AC1:AC2) and Vrj,k are stored in
C VEPS(AC1:AC2,VC1:VC2).
C The non-zero part of the result is stored in
C VEPS(VR1:VR2,VC1:VC2).
C
VR2 = IMUK(I)
AC1 = VR2 + 1
VR1 = AC1 - INUK(I)
AR1 = AR1 - INUK(I)
CALL DGEMM( 'No transpose', 'No transpose', INUK(I),
$ DIF, AC2-VR2, ONE, A(AR1,AC1), LDA,
$ VEPS(AC1,VC1), LDVEPS, ONE, VEPS(VR1,VC1),
$ LDVEPS )
120 CONTINUE
C
ER1 = 1
C
DO 140 I = 1, J - K - 1
C
C Compute the second part of Vij,k+1 in normal order:
C Vij,k+1 := Sum(r=i+1,..,j-k) EHir*Vrj,k.
C The non-zero parts of EHir are stored in
C E(ER1:ER1+nu(i)-1,EC1:AC2) and Vrj,k are stored in
C VEPS(EC1:AC2,VC1:VC2).
C The non-zero part of the result is stored in
C VEPS(VR1:VR2,VC2+1:VC2+DIF), where
C DIF = VC2 - VC1 + 1 = mu(j) - nu(j).
C This code portion also computes Vij,k+1 for i + k = j.
C
VR2 = IMUK(I)
EC1 = VR2 + 1
VR1 = EC1 - INUK(I)
CALL DGEMM( 'No transpose', 'No transpose', INUK(I),
$ DIF, AC2-VR2, ONE, E(ER1,EC1), LDE,
$ VEPS(EC1,VC1), LDVEPS, ZERO, VEPS(VR1,VC2+1),
$ LDVEPS )
ER1 = ER1 + INUK(I)
140 CONTINUE
C
VC1 = VC2 + 1
160 CONTINUE
C
VC1 = VC1 + DIF
180 CONTINUE
C
C Restore original contents of the array IMUK.
C
C Since, at the moment:
C IMUK(i) = Sum(j=1,...,i) mu(j), (i=1,...,NBLCKS),
C the original values are:
C mu(i) = IMUK(i) - IMUK(i-1) with IMUK(0 ) = 0.
C
SMUI1 = 0
C
DO 200 I = 1, NBLCKS
SMUI = IMUK(I)
IMUK(I) = SMUI - SMUI1
SMUI1 = SMUI
200 CONTINUE
C
RETURN
C *** Last line of MC03NY ***
END
|
