1 SUBROUTINE DPPT03( UPLO, N, A, AINV, WORK, LDWORK, RWORK, RCOND,
2 $ RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( * ), AINV( * ), RWORK( * ),
15 $ WORK( LDWORK, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DPPT03 computes the residual for a symmetric packed matrix times its
22 * inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * symmetric matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
39 * The original symmetric matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
43 * The (symmetric) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) DOUBLE PRECISION
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, J, JJ
68 DOUBLE PRECISION AINVNM, ANORM, EPS
69 * ..
70 * .. External Functions ..
71 LOGICAL LSAME
72 DOUBLE PRECISION DLAMCH, DLANGE, DLANSP
73 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSP
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL DCOPY, DSPMV
80 * ..
81 * .. Executable Statements ..
82 *
83 * Quick exit if N = 0.
84 *
85 IF( N.LE.0 ) THEN
86 RCOND = ONE
87 RESID = ZERO
88 RETURN
89 END IF
90 *
91 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
92 *
93 EPS = DLAMCH( 'Epsilon' )
94 ANORM = DLANSP( '1', UPLO, N, A, RWORK )
95 AINVNM = DLANSP( '1', UPLO, N, AINV, RWORK )
96 IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN
97 RCOND = ZERO
98 RESID = ONE / EPS
99 RETURN
100 END IF
101 RCOND = ( ONE / ANORM ) / AINVNM
102 *
103 * UPLO = 'U':
104 * Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
105 * expand it to a full matrix, then multiply by A one column at a
106 * time, moving the result one column to the left.
107 *
108 IF( LSAME( UPLO, 'U' ) ) THEN
109 *
110 * Copy AINV
111 *
112 JJ = 1
113 DO 10 J = 1, N - 1
114 CALL DCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
115 CALL DCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK )
116 JJ = JJ + J
117 10 CONTINUE
118 JJ = ( ( N-1 )*N ) / 2 + 1
119 CALL DCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK )
120 *
121 * Multiply by A
122 *
123 DO 20 J = 1, N - 1
124 CALL DSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO,
125 $ WORK( 1, J ), 1 )
126 20 CONTINUE
127 CALL DSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO,
128 $ WORK( 1, N ), 1 )
129 *
130 * UPLO = 'L':
131 * Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
132 * and multiply by A, moving each column to the right.
133 *
134 ELSE
135 *
136 * Copy AINV
137 *
138 CALL DCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK )
139 JJ = N + 1
140 DO 30 J = 2, N
141 CALL DCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
142 CALL DCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK )
143 JJ = JJ + N - J + 1
144 30 CONTINUE
145 *
146 * Multiply by A
147 *
148 DO 40 J = N, 2, -1
149 CALL DSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO,
150 $ WORK( 1, J ), 1 )
151 40 CONTINUE
152 CALL DSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO,
153 $ WORK( 1, 1 ), 1 )
154 *
155 END IF
156 *
157 * Add the identity matrix to WORK .
158 *
159 DO 50 I = 1, N
160 WORK( I, I ) = WORK( I, I ) + ONE
161 50 CONTINUE
162 *
163 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
164 *
165 RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK )
166 *
167 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
168 *
169 RETURN
170 *
171 * End of DPPT03
172 *
173 END
2 $ RESID )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER UPLO
10 INTEGER LDWORK, N
11 DOUBLE PRECISION RCOND, RESID
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( * ), AINV( * ), RWORK( * ),
15 $ WORK( LDWORK, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DPPT03 computes the residual for a symmetric packed matrix times its
22 * inverse:
23 * norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * UPLO (input) CHARACTER*1
30 * Specifies whether the upper or lower triangular part of the
31 * symmetric matrix A is stored:
32 * = 'U': Upper triangular
33 * = 'L': Lower triangular
34 *
35 * N (input) INTEGER
36 * The number of rows and columns of the matrix A. N >= 0.
37 *
38 * A (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
39 * The original symmetric matrix A, stored as a packed
40 * triangular matrix.
41 *
42 * AINV (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
43 * The (symmetric) inverse of the matrix A, stored as a packed
44 * triangular matrix.
45 *
46 * WORK (workspace) DOUBLE PRECISION array, dimension (LDWORK,N)
47 *
48 * LDWORK (input) INTEGER
49 * The leading dimension of the array WORK. LDWORK >= max(1,N).
50 *
51 * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
52 *
53 * RCOND (output) DOUBLE PRECISION
54 * The reciprocal of the condition number of A, computed as
55 * ( 1/norm(A) ) / norm(AINV).
56 *
57 * RESID (output) DOUBLE PRECISION
58 * norm(I - A*AINV) / ( N * norm(A) * norm(AINV) * EPS )
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63 DOUBLE PRECISION ZERO, ONE
64 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
65 * ..
66 * .. Local Scalars ..
67 INTEGER I, J, JJ
68 DOUBLE PRECISION AINVNM, ANORM, EPS
69 * ..
70 * .. External Functions ..
71 LOGICAL LSAME
72 DOUBLE PRECISION DLAMCH, DLANGE, DLANSP
73 EXTERNAL LSAME, DLAMCH, DLANGE, DLANSP
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC DBLE
77 * ..
78 * .. External Subroutines ..
79 EXTERNAL DCOPY, DSPMV
80 * ..
81 * .. Executable Statements ..
82 *
83 * Quick exit if N = 0.
84 *
85 IF( N.LE.0 ) THEN
86 RCOND = ONE
87 RESID = ZERO
88 RETURN
89 END IF
90 *
91 * Exit with RESID = 1/EPS if ANORM = 0 or AINVNM = 0.
92 *
93 EPS = DLAMCH( 'Epsilon' )
94 ANORM = DLANSP( '1', UPLO, N, A, RWORK )
95 AINVNM = DLANSP( '1', UPLO, N, AINV, RWORK )
96 IF( ANORM.LE.ZERO .OR. AINVNM.EQ.ZERO ) THEN
97 RCOND = ZERO
98 RESID = ONE / EPS
99 RETURN
100 END IF
101 RCOND = ( ONE / ANORM ) / AINVNM
102 *
103 * UPLO = 'U':
104 * Copy the leading N-1 x N-1 submatrix of AINV to WORK(1:N,2:N) and
105 * expand it to a full matrix, then multiply by A one column at a
106 * time, moving the result one column to the left.
107 *
108 IF( LSAME( UPLO, 'U' ) ) THEN
109 *
110 * Copy AINV
111 *
112 JJ = 1
113 DO 10 J = 1, N - 1
114 CALL DCOPY( J, AINV( JJ ), 1, WORK( 1, J+1 ), 1 )
115 CALL DCOPY( J-1, AINV( JJ ), 1, WORK( J, 2 ), LDWORK )
116 JJ = JJ + J
117 10 CONTINUE
118 JJ = ( ( N-1 )*N ) / 2 + 1
119 CALL DCOPY( N-1, AINV( JJ ), 1, WORK( N, 2 ), LDWORK )
120 *
121 * Multiply by A
122 *
123 DO 20 J = 1, N - 1
124 CALL DSPMV( 'Upper', N, -ONE, A, WORK( 1, J+1 ), 1, ZERO,
125 $ WORK( 1, J ), 1 )
126 20 CONTINUE
127 CALL DSPMV( 'Upper', N, -ONE, A, AINV( JJ ), 1, ZERO,
128 $ WORK( 1, N ), 1 )
129 *
130 * UPLO = 'L':
131 * Copy the trailing N-1 x N-1 submatrix of AINV to WORK(1:N,1:N-1)
132 * and multiply by A, moving each column to the right.
133 *
134 ELSE
135 *
136 * Copy AINV
137 *
138 CALL DCOPY( N-1, AINV( 2 ), 1, WORK( 1, 1 ), LDWORK )
139 JJ = N + 1
140 DO 30 J = 2, N
141 CALL DCOPY( N-J+1, AINV( JJ ), 1, WORK( J, J-1 ), 1 )
142 CALL DCOPY( N-J, AINV( JJ+1 ), 1, WORK( J, J ), LDWORK )
143 JJ = JJ + N - J + 1
144 30 CONTINUE
145 *
146 * Multiply by A
147 *
148 DO 40 J = N, 2, -1
149 CALL DSPMV( 'Lower', N, -ONE, A, WORK( 1, J-1 ), 1, ZERO,
150 $ WORK( 1, J ), 1 )
151 40 CONTINUE
152 CALL DSPMV( 'Lower', N, -ONE, A, AINV( 1 ), 1, ZERO,
153 $ WORK( 1, 1 ), 1 )
154 *
155 END IF
156 *
157 * Add the identity matrix to WORK .
158 *
159 DO 50 I = 1, N
160 WORK( I, I ) = WORK( I, I ) + ONE
161 50 CONTINUE
162 *
163 * Compute norm(I - A*AINV) / (N * norm(A) * norm(AINV) * EPS)
164 *
165 RESID = DLANGE( '1', N, N, WORK, LDWORK, RWORK )
166 *
167 RESID = ( ( RESID*RCOND ) / EPS ) / DBLE( N )
168 *
169 RETURN
170 *
171 * End of DPPT03
172 *
173 END