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
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
|
<HTML>
<HEAD><TITLE>MB04VD - SLICOT Library Routine Documentation</TITLE>
</HEAD>
<BODY>
<H2><A Name="MB04VD">MB04VD</A></H2>
<H3>
Upper block triangular form for a rectangular pencil sE-A, with E in column echelon form (added functionality)
</H3>
<A HREF ="#Specification"><B>[Specification]</B></A>
<A HREF ="#Arguments"><B>[Arguments]</B></A>
<A HREF ="#Method"><B>[Method]</B></A>
<A HREF ="#References"><B>[References]</B></A>
<A HREF ="#Comments"><B>[Comments]</B></A>
<A HREF ="#Example"><B>[Example]</B></A>
<P>
<B><FONT SIZE="+1">Purpose</FONT></B>
<PRE>
To compute orthogonal transformations Q and Z such that the
transformed pencil Q'(sE-A)Z is in upper block triangular form,
where E is an M-by-N matrix in column echelon form (see SLICOT
Library routine MB04UD) and A is an M-by-N matrix.
If MODE = 'B', then the matrices A and E are transformed into the
following generalized Schur form by unitary transformations Q1
and Z1 :
| sE(eps,inf)-A(eps,inf) | X |
Q1'(sE-A)Z1 = |------------------------|------------|. (1)
| O | sE(r)-A(r) |
The pencil sE(eps,inf)-A(eps,inf) is in staircase form, and it
contains all Kronecker column indices and infinite elementary
divisors of the pencil sE-A. The pencil sE(r)-A(r) contains all
Kronecker row indices and elementary divisors of sE-A.
Note: X is a pencil.
If MODE = 'T', then the submatrices having full row and column
rank in the pencil sE(eps,inf)-A(eps,inf) in (1) are
triangularized by applying unitary transformations Q2 and Z2 to
Q1'*(sE-A)*Z1.
If MODE = 'S', then the pencil sE(eps,inf)-A(eps,inf) in (1) is
separated into sE(eps)-A(eps) and sE(inf)-A(inf) by applying
unitary transformations Q3 and Z3 to Q2'*Q1'*(sE-A)*Z1*Z2.
This gives
| sE(eps)-A(eps) | X | X |
|----------------|----------------|------------|
| O | sE(inf)-A(inf) | X |
Q'(sE-A)Z =|=================================|============| (2)
| | |
| O | sE(r)-A(r) |
where Q = Q1*Q2*Q3 and Z = Z1*Z2*Z3.
Note: the pencil sE(r)-A(r) is not reduced further.
</PRE>
<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
<PRE>
SUBROUTINE MB04VD( MODE, JOBQ, JOBZ, M, N, RANKE, A, LDA, E, LDE,
$ Q, LDQ, Z, LDZ, ISTAIR, NBLCKS, NBLCKI, IMUK,
$ INUK, IMUK0, MNEI, TOL, IWORK, INFO )
C .. Scalar Arguments ..
CHARACTER JOBQ, JOBZ, MODE
INTEGER INFO, LDA, LDE, LDQ, LDZ, M, N, NBLCKI, NBLCKS,
$ RANKE
DOUBLE PRECISION TOL
C .. Array Arguments ..
INTEGER IMUK(*), IMUK0(*), INUK(*), ISTAIR(*), IWORK(*),
$ MNEI(*)
DOUBLE PRECISION A(LDA,*), E(LDE,*), Q(LDQ,*), Z(LDZ,*)
</PRE>
<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
<P>
<B>Mode Parameters</B>
<PRE>
MODE CHARACTER*1
Specifies the desired structure of the transformed
pencil Q'(sE-A)Z to be computed as follows:
= 'B': Basic reduction given by (1);
= 'T': Further reduction of (1) to triangular form;
= 'S': Further separation of sE(eps,inf)-A(eps,inf)
in (1) into the two pencils in (2).
JOBQ CHARACTER*1
Indicates whether the user wishes to accumulate in a
matrix Q the orthogonal row transformations, as follows:
= 'N': Do not form Q;
= 'I': Q is initialized to the unit matrix and the
orthogonal transformation matrix Q is returned;
= 'U': The given matrix Q is updated by the orthogonal
row transformations used in the reduction.
JOBZ CHARACTER*1
Indicates whether the user wishes to accumulate in a
matrix Z the orthogonal column transformations, as
follows:
= 'N': Do not form Z;
= 'I': Z is initialized to the unit matrix and the
orthogonal transformation matrix Z is returned;
= 'U': The given matrix Z is updated by the orthogonal
transformations used in the reduction.
</PRE>
<B>Input/Output Parameters</B>
<PRE>
M (input) INTEGER
The number of rows in the matrices A, E and the order of
the matrix Q. M >= 0.
N (input) INTEGER
The number of columns in the matrices A, E and the order
of the matrix Z. N >= 0.
RANKE (input) INTEGER
The rank of the matrix E in column echelon form.
RANKE >= 0.
A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the leading M-by-N part of this array must
contain the matrix to be row compressed.
On exit, the leading M-by-N part of this array contains
the matrix that has been row compressed while keeping
matrix E in column echelon form.
LDA INTEGER
The leading dimension of array A. LDA >= MAX(1,M).
E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
On entry, the leading M-by-N part of this array must
contain the matrix in column echelon form to be
transformed equivalent to matrix A.
On exit, the leading M-by-N part of this array contains
the matrix that has been transformed equivalent to matrix
A.
LDE INTEGER
The leading dimension of array E. LDE >= MAX(1,M).
Q (input/output) DOUBLE PRECISION array, dimension (LDQ,*)
On entry, if JOBQ = 'U', then the leading M-by-M part of
this array must contain a given matrix Q (e.g. from a
previous call to another SLICOT routine), and on exit, the
leading M-by-M part of this array contains the product of
the input matrix Q and the row transformation matrix used
to transform the rows of matrices A and E.
On exit, if JOBQ = 'I', then the leading M-by-M part of
this array contains the matrix of accumulated orthogonal
row transformations performed.
If JOBQ = 'N', the array Q is not referenced and can be
supplied as a dummy array (i.e. set parameter LDQ = 1 and
declare this array to be Q(1,1) in the calling program).
LDQ INTEGER
The leading dimension of array Q. If JOBQ = 'U' or
JOBQ = 'I', LDQ >= MAX(1,M); if JOBQ = 'N', LDQ >= 1.
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,*)
On entry, if JOBZ = 'U', then the leading N-by-N part of
this array must contain a given matrix Z (e.g. from a
previous call to another SLICOT routine), and on exit, the
leading N-by-N part of this array contains the product of
the input matrix Z and the column transformation matrix
used to transform the columns of matrices A and E.
On exit, if JOBZ = 'I', then the leading N-by-N part of
this array contains the matrix of accumulated orthogonal
column transformations performed.
If JOBZ = 'N', the array Z is not referenced and can be
supplied as a dummy array (i.e. set parameter LDZ = 1 and
declare this array to be Z(1,1) in the calling program).
LDZ INTEGER
The leading dimension of array Z. If JOBZ = 'U' or
JOBZ = 'I', LDZ >= MAX(1,N); if JOBZ = 'N', LDZ >= 1.
ISTAIR (input/output) INTEGER array, dimension (M)
On entry, this array must contain information on the
column echelon form of the unitary transformed matrix E.
Specifically, ISTAIR(i) must be set to +j if the first
non-zero element E(i,j) is a corner point and -j
otherwise, for i = 1,2,...,M.
On exit, this array contains no useful information.
NBLCKS (output) INTEGER
The number of submatrices having full row rank greater
than or equal to 0 detected in matrix A in the pencil
sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
NBLCKI (output) INTEGER
If MODE = 'S', the number of diagonal submatrices in the
pencil sE(inf)-A(inf). If MODE = 'B' or 'T' then
NBLCKI = 0.
IMUK (output) INTEGER array, dimension (MAX(N,M+1))
The leading NBLCKS elements of this array contain the
column dimensions mu(1),...,mu(NBLCKS) of the submatrices
having full column rank in the pencil sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
INUK (output) INTEGER array, dimension (MAX(N,M+1))
The leading NBLCKS elements of this array contain the
row dimensions nu(1),...,nu(NBLCKS) of the submatrices
having full row rank in the pencil sE(x)-A(x),
where x = eps,inf if MODE = 'B' or 'T',
or x = eps if MODE = 'S'.
IMUK0 (output) INTEGER array, dimension (limuk0),
where limuk0 = N if MODE = 'S' and 1, otherwise.
If MODE = 'S', then the leading NBLCKI elements of this
array contain the dimensions mu0(1),...,mu0(NBLCKI)
of the square diagonal submatrices in the pencil
sE(inf)-A(inf).
Otherwise, IMUK0 is not referenced and can be supplied
as a dummy array.
MNEI (output) INTEGER array, dimension (3)
If MODE = 'B' or 'T' then
MNEI(1) contains the row dimension of
sE(eps,inf)-A(eps,inf);
MNEI(2) contains the column dimension of
sE(eps,inf)-A(eps,inf);
MNEI(3) = 0.
If MODE = 'S', then
MNEI(1) contains the row dimension of sE(eps)-A(eps);
MNEI(2) contains the column dimension of sE(eps)-A(eps);
MNEI(3) contains the order of the regular pencil
sE(inf)-A(inf).
</PRE>
<B>Tolerances</B>
<PRE>
TOL DOUBLE PRECISION
A tolerance below which matrix elements are considered
to be zero. If the user sets TOL to be less than (or
equal to) zero then the tolerance is taken as
EPS * MAX( ABS(A(I,J)), ABS(E(I,J)) ), where EPS is the
machine precision (see LAPACK Library routine DLAMCH),
I = 1,2,...,M and J = 1,2,...,N.
</PRE>
<B>Workspace</B>
<PRE>
IWORK INTEGER array, dimension (N)
</PRE>
<B>Error Indicator</B>
<PRE>
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if incorrect rank decisions were revealed during the
triangularization phase. This failure is not likely
to occur. The possible values are:
= 1: if incorrect dimensions of a full column rank
submatrix;
= 2: if incorrect dimensions of a full row rank
submatrix.
</PRE>
<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
<PRE>
Let sE - A be an arbitrary pencil. Prior to calling the routine,
this pencil must be transformed into a pencil with E in column
echelon form. This may be accomplished by calling the SLICOT
Library routine MB04UD. Depending on the value of MODE,
submatrices of A and E are then reduced to one of the forms
described above. Further details can be found in [1].
</PRE>
<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
<PRE>
[1] Beelen, Th. and Van Dooren, P.
An improved algorithm for the computation of Kronecker's
canonical form of a singular pencil.
Linear Algebra and Applications, 105, pp. 9-65, 1988.
</PRE>
<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
<PRE>
It is shown in [1] that the algorithm is numerically backward
stable. The operations count is proportional to (MAX(M,N))**3.
</PRE>
<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
<PRE>
The difference mu(k)-nu(k), for k = 1,2,...,NBLCKS, is the number
of elementary Kronecker blocks of size k x (k+1).
If MODE = 'B' or 'T' on entry, then the difference nu(k)-mu(k+1),
for k = 1,2,...,NBLCKS, is the number of infinite elementary
divisors of degree k (with mu(NBLCKS+1) = 0).
If MODE = 'S' on entry, then the difference mu0(k)-mu0(k+1),
for k = 1,2,...,NBLCKI, is the number of infinite elementary
divisors of degree k (with mu0(NBLCKI+1) = 0).
In the pencil sE(r)-A(r), the pencils sE(f)-A(f) and
sE(eta)-A(eta) can be separated by pertransposing the pencil
sE(r)-A(r) and calling the routine with MODE set to 'B'. The
result has got to be pertransposed again. (For more details see
[1]).
</PRE>
<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
<P>
<B>Program Text</B>
<PRE>
* MB04VD EXAMPLE PROGRAM TEXT
*
* .. Parameters ..
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER MMAX, NMAX
PARAMETER ( MMAX = 20, NMAX = 20 )
INTEGER LDA, LDE, LDQ, LDZ
PARAMETER ( LDA = MMAX, LDE = MMAX, LDQ = MMAX,
$ LDZ = NMAX )
INTEGER LINUK
PARAMETER ( LINUK = MAX( NMAX,MMAX+1 ) )
* PARAMETER ( LINUK = NMAX+MMAX+1 )
INTEGER LIWORK
PARAMETER ( LIWORK = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( NMAX,MMAX ) )
* PARAMETER ( LDWORK = NMAX+MMAX )
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
* .. Local Scalars ..
DOUBLE PRECISION TOL
INTEGER I, INFO, J, M, N, NBLCKI, NBLCKS, RANKE
CHARACTER*1 JOBQ, JOBZ, MODE
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), E(LDE,NMAX),
$ Q(LDQ,MMAX), Z(LDZ,NMAX)
INTEGER IMUK(LINUK), IMUK0(NMAX), INUK(LINUK),
$ ISTAIR(MMAX), IWORK(LIWORK), MNEI(3)
C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL MB04UD, MB04VD
* .. Intrinsic Functions ..
INTRINSIC MAX
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) M, N, TOL, MODE
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99984 ) M
ELSE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99983 ) N
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,M )
READ ( NIN, FMT = * ) ( ( E(I,J), J = 1,N ), I = 1,M )
JOBQ = 'I'
JOBZ = 'I'
* Reduce E to column echelon form and compute Q'*A*Z.
CALL MB04UD( JOBQ, JOBZ, M, N, A, LDA, E, LDE, Q, LDQ, Z, LDZ,
$ RANKE, ISTAIR, TOL, DWORK, INFO )
JOBQ = 'U'
JOBZ = 'U'
*
IF ( INFO.EQ.0 ) THEN
* Compute a unitary transformed pencil Q'*(s*E-A)*Z.
CALL MB04VD( MODE, JOBQ, JOBZ, M, N, RANKE, A, LDA, E, LDE,
$ Q, LDQ, Z, LDZ, ISTAIR, NBLCKS, NBLCKI, IMUK,
$ INUK, IMUK0, MNEI, TOL, IWORK, INFO )
*
IF ( INFO.EQ.0 ) THEN
WRITE ( NOUT, FMT = 99996 )
WRITE ( NOUT, FMT = 99995 )
DO 140 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( Q(I,J), J = 1,M )
140 CONTINUE
WRITE ( NOUT, FMT = 99993 )
DO 160 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( E(I,J), J = 1,N )
160 CONTINUE
WRITE ( NOUT, FMT = 99992 )
DO 180 I = 1, M
WRITE ( NOUT, FMT = 99994 ) ( A(I,J), J = 1,N )
180 CONTINUE
WRITE ( NOUT, FMT = 99991 )
DO 200 I = 1, N
WRITE ( NOUT, FMT = 99994 ) ( Z(I,J), J = 1,N )
200 CONTINUE
WRITE ( NOUT, FMT = 99990 ) NBLCKS
IF ( .NOT. LSAME( MODE, 'S' ) ) THEN
WRITE ( NOUT, FMT = 99989 ) ( IMUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99988 ) ( INUK(I), I = 1,NBLCKS )
ELSE
WRITE ( NOUT, FMT = 99987 ) ( IMUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99986 ) ( INUK(I), I = 1,NBLCKS )
WRITE ( NOUT, FMT = 99982 ) ( IMUK0(I), I = 1,NBLCKI )
WRITE ( NOUT, FMT = 99985 ) ( MNEI(I), I = 1,3 )
END IF
ELSE
WRITE ( NOUT, FMT = 99998 ) INFO
END IF
ELSE
WRITE ( NOUT, FMT = 99997 ) INFO
END IF
END IF
STOP
*
99999 FORMAT (' MB04VD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04VD = ',I2)
99997 FORMAT (' INFO on exit from MB04UD = ',I2)
99996 FORMAT (' The unitary transformed pencil is Q''*(s*E-A)*Z, where',
$ /)
99995 FORMAT (' Matrix Q',/)
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' Matrix E',/)
99992 FORMAT (/' Matrix A',/)
99991 FORMAT (/' Matrix Z',/)
99990 FORMAT (/' The number of submatrices having full row rank detect',
$ 'ed in matrix A = ',I3)
99989 FORMAT (/' The column dimensions of the submatrices having full ',
$ 'column rank in the pencil',/' sE(eps,inf) - A(eps,inf) a',
$ 're',/20(1X,I5))
99988 FORMAT (/' The row dimensions of the submatrices having full row',
$ ' rank in the pencil',/' sE(eps,inf) - A(eps,inf) are',
$ /20(1X,I5))
99987 FORMAT (/' The column dimensions of the submatrices having full ',
$ 'column rank in the pencil',/' sE(eps) - A(eps) are',
$ /20(1X,I5))
99986 FORMAT (/' The row dimensions of the submatrices having full row',
$ ' rank in the pencil',/' sE(eps) - A(eps) are',/20(1X,I5))
99985 FORMAT (/' MNEI is ',/20(1X,I5))
99984 FORMAT (/' M is out of range.',/' M = ',I5)
99983 FORMAT (/' N is out of range.',/' N = ',I5)
99982 FORMAT (/' The orders of the diagonal submatrices in the pencil ',
$ 'sE(inf) - A(inf) are',/20(1X,I5))
END
</PRE>
<B>Program Data</B>
<PRE>
MB04VD EXAMPLE PROGRAM DATA
2 4 0.0 S
1.0 0.0 -1.0 0.0
1.0 1.0 0.0 -1.0
0.0 -1.0 0.0 0.0
0.0 -1.0 0.0 0.0
</PRE>
<B>Program Results</B>
<PRE>
MB04VD EXAMPLE PROGRAM RESULTS
The unitary transformed pencil is Q'*(s*E-A)*Z, where
Matrix Q
0.7071 -0.7071
0.7071 0.7071
Matrix E
0.0000 0.0000 -1.1547 0.8165
0.0000 0.0000 0.0000 0.0000
Matrix A
0.0000 1.7321 0.5774 -0.4082
0.0000 0.0000 0.0000 -1.2247
Matrix Z
0.5774 0.8165 0.0000 0.0000
0.0000 0.0000 0.8165 -0.5774
0.5774 -0.4082 -0.4082 -0.5774
0.5774 -0.4082 0.4082 0.5774
The number of submatrices having full row rank detected in matrix A = 2
The column dimensions of the submatrices having full column rank in the pencil
sE(eps) - A(eps) are
2 1
The row dimensions of the submatrices having full row rank in the pencil
sE(eps) - A(eps) are
1 0
The orders of the diagonal submatrices in the pencil sE(inf) - A(inf) are
1
MNEI is
1 3 1
</PRE>
<HR>
<p>
<A HREF=..\libindex.html><B>Return to index</B></A></BODY>
</HTML>
|
