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
|
SUBROUTINE MB01TD( N, A, LDA, B, LDB, DWORK, INFO )
C
C RELEASE 4.0, WGS COPYRIGHT 1999.
C
C PURPOSE
C
C To compute the matrix product A * B, where A and B are upper
C quasi-triangular matrices (that is, block upper triangular with
C 1-by-1 or 2-by-2 diagonal blocks) with the same structure.
C The result is returned in the array B.
C
C ARGUMENTS
C
C Input/Output Parameters
C
C N (input) INTEGER
C The order of the matrices A and B. N >= 0.
C
C A (input) DOUBLE PRECISION array, dimension (LDA,N)
C The leading N-by-N part of this array must contain the
C upper quasi-triangular matrix A. The elements below the
C subdiagonal are not referenced.
C
C LDA INTEGER
C The leading dimension of the array A. LDA >= max(1,N).
C
C B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
C On entry, the leading N-by-N part of this array must
C contain the upper quasi-triangular matrix B, with the same
C structure as matrix A.
C On exit, the leading N-by-N part of this array contains
C the computed product A * B, with the same structure as
C on entry.
C The elements below the subdiagonal are not referenced.
C
C LDB INTEGER
C The leading dimension of the array B. LDB >= max(1,N).
C
C Workspace
C
C DWORK DOUBLE PRECISION array, dimension (N-1)
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 = 1: if the matrices A and B have not the same structure,
C and/or A and B are not upper quasi-triangular.
C
C METHOD
C
C The matrix product A * B is computed column by column, using
C BLAS 2 and BLAS 1 operations.
C
C FURTHER COMMENTS
C
C This routine can be used, for instance, for computing powers of
C a real Schur form matrix.
C
C CONTRIBUTOR
C
C V. Sima, Katholieke Univ. Leuven, Belgium, June 1998.
C
C REVISIONS
C
C V. Sima, Feb. 2000.
C
C KEYWORDS
C
C Elementary matrix operations, matrix operations.
C
C ******************************************************************
C
C .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
C .. Scalar Arguments ..
INTEGER INFO, LDA, LDB, N
C .. Array Arguments ..
DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*)
C .. Local Scalars ..
INTEGER I, J, JMIN, JMNM
C .. External Subroutines ..
EXTERNAL DAXPY, DTRMV, XERBLA
C .. Intrinsic Functions ..
INTRINSIC MAX, MIN
C
C .. Executable Statements ..
C
C Test the input scalar arguments.
C
INFO = 0
IF( N.LT.0 ) THEN
INFO = -1
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -3
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -5
END IF
C
IF ( INFO.NE.0 ) THEN
C
C Error return.
C
CALL XERBLA( 'MB01TD', -INFO )
RETURN
END IF
C
C Quick return, if possible.
C
IF ( N.EQ.0 ) THEN
RETURN
ELSE IF ( N.EQ.1 ) THEN
B(1,1) = A(1,1)*B(1,1)
RETURN
END IF
C
C Test the upper quasi-triangular structure of A and B for identity.
C
DO 10 I = 1, N - 1
IF ( A(I+1,I).EQ.ZERO ) THEN
IF ( B(I+1,I).NE.ZERO ) THEN
INFO = 1
RETURN
END IF
ELSE IF ( I.LT.N-1 ) THEN
IF ( A(I+2,I+1).NE.ZERO ) THEN
INFO = 1
RETURN
END IF
END IF
10 CONTINUE
C
DO 30 J = 1, N
JMIN = MIN( J+1, N )
JMNM = MIN( JMIN, N-1 )
C
C Compute the contribution of the subdiagonal of A to the
C j-th column of the product.
C
DO 20 I = 1, JMNM
DWORK(I) = A(I+1,I)*B(I,J)
20 CONTINUE
C
C Multiply the upper triangle of A by the j-th column of B,
C and add to the above result.
C
CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', JMIN, A, LDA,
$ B(1,J), 1 )
CALL DAXPY( JMNM, ONE, DWORK, 1, B(2,J), 1 )
30 CONTINUE
C
RETURN
C *** Last line of MB01TD ***
END
|
