1 SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U,
2 $ WORK, LWORK, RWORK, RESULT )
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 INTEGER LDA, LDB, LWORK, M, N, P
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 *
14 * Purpose
15 * =======
16 *
17 * ZGLMTS tests ZGGGLM - a subroutine for solving the generalized
18 * linear model problem.
19 *
20 * Arguments
21 * =========
22 *
23 * N (input) INTEGER
24 * The number of rows of the matrices A and B. N >= 0.
25 *
26 * M (input) INTEGER
27 * The number of columns of the matrix A. M >= 0.
28 *
29 * P (input) INTEGER
30 * The number of columns of the matrix B. P >= 0.
31 *
32 * A (input) COMPLEX*16 array, dimension (LDA,M)
33 * The N-by-M matrix A.
34 *
35 * AF (workspace) COMPLEX*16 array, dimension (LDA,M)
36 *
37 * LDA (input) INTEGER
38 * The leading dimension of the arrays A, AF. LDA >= max(M,N).
39 *
40 * B (input) COMPLEX*16 array, dimension (LDB,P)
41 * The N-by-P matrix A.
42 *
43 * BF (workspace) COMPLEX*16 array, dimension (LDB,P)
44 *
45 * LDB (input) INTEGER
46 * The leading dimension of the arrays B, BF. LDB >= max(P,N).
47 *
48 * D (input) COMPLEX*16 array, dimension( N )
49 * On input, the left hand side of the GLM.
50 *
51 * DF (workspace) COMPLEX*16 array, dimension( N )
52 *
53 * X (output) COMPLEX*16 array, dimension( M )
54 * solution vector X in the GLM problem.
55 *
56 * U (output) COMPLEX*16 array, dimension( P )
57 * solution vector U in the GLM problem.
58 *
59 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
60 *
61 * LWORK (input) INTEGER
62 * The dimension of the array WORK.
63 *
64 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
65 *
66 * RESULT (output) DOUBLE PRECISION
67 * The test ratio:
68 * norm( d - A*x - B*u )
69 * RESULT = -----------------------------------------
70 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
71 *
72 * ====================================================================
73 *
74 DOUBLE PRECISION RWORK( * )
75 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
76 $ BF( LDB, * ), D( * ), DF( * ), U( * ),
77 $ WORK( LWORK ), X( * )
78 * ..
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO
81 PARAMETER ( ZERO = 0.0D+0 )
82 COMPLEX*16 CONE
83 PARAMETER ( CONE = 1.0D+0 )
84 * ..
85 * .. Local Scalars ..
86 INTEGER INFO
87 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
88 * ..
89 * .. External Functions ..
90 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
91 EXTERNAL DLAMCH, DZASUM, ZLANGE
92 * ..
93 * .. External Subroutines ..
94 *
95 EXTERNAL ZCOPY, ZGEMV, ZGGGLM, ZLACPY
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX
99 * ..
100 * .. Executable Statements ..
101 *
102 EPS = DLAMCH( 'Epsilon' )
103 UNFL = DLAMCH( 'Safe minimum' )
104 ANORM = MAX( ZLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
105 BNORM = MAX( ZLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
106 *
107 * Copy the matrices A and B to the arrays AF and BF,
108 * and the vector D the array DF.
109 *
110 CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA )
111 CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB )
112 CALL ZCOPY( N, D, 1, DF, 1 )
113 *
114 * Solve GLM problem
115 *
116 CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
117 $ INFO )
118 *
119 * Test the residual for the solution of LSE
120 *
121 * norm( d - A*x - B*u )
122 * RESULT = -----------------------------------------
123 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
124 *
125 CALL ZCOPY( N, D, 1, DF, 1 )
126 CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, DF,
127 $ 1 )
128 *
129 CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, DF,
130 $ 1 )
131 *
132 DNORM = DZASUM( N, DF, 1 )
133 XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 )
134 YNORM = ANORM + BNORM
135 *
136 IF( XNORM.LE.ZERO ) THEN
137 RESULT = ZERO
138 ELSE
139 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
140 END IF
141 *
142 RETURN
143 *
144 * End of ZGLMTS
145 *
146 END
2 $ WORK, LWORK, RWORK, RESULT )
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 INTEGER LDA, LDB, LWORK, M, N, P
10 DOUBLE PRECISION RESULT
11 * ..
12 * .. Array Arguments ..
13 *
14 * Purpose
15 * =======
16 *
17 * ZGLMTS tests ZGGGLM - a subroutine for solving the generalized
18 * linear model problem.
19 *
20 * Arguments
21 * =========
22 *
23 * N (input) INTEGER
24 * The number of rows of the matrices A and B. N >= 0.
25 *
26 * M (input) INTEGER
27 * The number of columns of the matrix A. M >= 0.
28 *
29 * P (input) INTEGER
30 * The number of columns of the matrix B. P >= 0.
31 *
32 * A (input) COMPLEX*16 array, dimension (LDA,M)
33 * The N-by-M matrix A.
34 *
35 * AF (workspace) COMPLEX*16 array, dimension (LDA,M)
36 *
37 * LDA (input) INTEGER
38 * The leading dimension of the arrays A, AF. LDA >= max(M,N).
39 *
40 * B (input) COMPLEX*16 array, dimension (LDB,P)
41 * The N-by-P matrix A.
42 *
43 * BF (workspace) COMPLEX*16 array, dimension (LDB,P)
44 *
45 * LDB (input) INTEGER
46 * The leading dimension of the arrays B, BF. LDB >= max(P,N).
47 *
48 * D (input) COMPLEX*16 array, dimension( N )
49 * On input, the left hand side of the GLM.
50 *
51 * DF (workspace) COMPLEX*16 array, dimension( N )
52 *
53 * X (output) COMPLEX*16 array, dimension( M )
54 * solution vector X in the GLM problem.
55 *
56 * U (output) COMPLEX*16 array, dimension( P )
57 * solution vector U in the GLM problem.
58 *
59 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
60 *
61 * LWORK (input) INTEGER
62 * The dimension of the array WORK.
63 *
64 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
65 *
66 * RESULT (output) DOUBLE PRECISION
67 * The test ratio:
68 * norm( d - A*x - B*u )
69 * RESULT = -----------------------------------------
70 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
71 *
72 * ====================================================================
73 *
74 DOUBLE PRECISION RWORK( * )
75 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
76 $ BF( LDB, * ), D( * ), DF( * ), U( * ),
77 $ WORK( LWORK ), X( * )
78 * ..
79 * .. Parameters ..
80 DOUBLE PRECISION ZERO
81 PARAMETER ( ZERO = 0.0D+0 )
82 COMPLEX*16 CONE
83 PARAMETER ( CONE = 1.0D+0 )
84 * ..
85 * .. Local Scalars ..
86 INTEGER INFO
87 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
88 * ..
89 * .. External Functions ..
90 DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE
91 EXTERNAL DLAMCH, DZASUM, ZLANGE
92 * ..
93 * .. External Subroutines ..
94 *
95 EXTERNAL ZCOPY, ZGEMV, ZGGGLM, ZLACPY
96 * ..
97 * .. Intrinsic Functions ..
98 INTRINSIC MAX
99 * ..
100 * .. Executable Statements ..
101 *
102 EPS = DLAMCH( 'Epsilon' )
103 UNFL = DLAMCH( 'Safe minimum' )
104 ANORM = MAX( ZLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
105 BNORM = MAX( ZLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
106 *
107 * Copy the matrices A and B to the arrays AF and BF,
108 * and the vector D the array DF.
109 *
110 CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA )
111 CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB )
112 CALL ZCOPY( N, D, 1, DF, 1 )
113 *
114 * Solve GLM problem
115 *
116 CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
117 $ INFO )
118 *
119 * Test the residual for the solution of LSE
120 *
121 * norm( d - A*x - B*u )
122 * RESULT = -----------------------------------------
123 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
124 *
125 CALL ZCOPY( N, D, 1, DF, 1 )
126 CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, DF,
127 $ 1 )
128 *
129 CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, DF,
130 $ 1 )
131 *
132 DNORM = DZASUM( N, DF, 1 )
133 XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 )
134 YNORM = ANORM + BNORM
135 *
136 IF( XNORM.LE.ZERO ) THEN
137 RESULT = ZERO
138 ELSE
139 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
140 END IF
141 *
142 RETURN
143 *
144 * End of ZGLMTS
145 *
146 END