1 SUBROUTINE DGLMTS( 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 * DGLMTS tests DGGGLM - 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) DOUBLE PRECISION array, dimension (LDA,M)
33 * The N-by-M matrix A.
34 *
35 * AF (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDB,P)
41 * The N-by-P matrix A.
42 *
43 * BF (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension( N )
49 * On input, the left hand side of the GLM.
50 *
51 * DF (workspace) DOUBLE PRECISION array, dimension( N )
52 *
53 * X (output) DOUBLE PRECISION array, dimension( M )
54 * solution vector X in the GLM problem.
55 *
56 * U (output) DOUBLE PRECISION array, dimension( P )
57 * solution vector U in the GLM problem.
58 *
59 * WORK (workspace) DOUBLE PRECISION 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 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
75 $ BF( LDB, * ), D( * ), DF( * ), RWORK( * ),
76 $ U( * ), WORK( LWORK ), X( * )
77 * ..
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE
80 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER INFO
84 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
85 * ..
86 * .. External Functions ..
87 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
88 EXTERNAL DASUM, DLAMCH, DLANGE
89 * ..
90 * .. External Subroutines ..
91 *
92 EXTERNAL DCOPY, DGEMV, DGGGLM, DLACPY
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 EPS = DLAMCH( 'Epsilon' )
100 UNFL = DLAMCH( 'Safe minimum' )
101 ANORM = MAX( DLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
102 BNORM = MAX( DLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
103 *
104 * Copy the matrices A and B to the arrays AF and BF,
105 * and the vector D the array DF.
106 *
107 CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA )
108 CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB )
109 CALL DCOPY( N, D, 1, DF, 1 )
110 *
111 * Solve GLM problem
112 *
113 CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
114 $ INFO )
115 *
116 * Test the residual for the solution of LSE
117 *
118 * norm( d - A*x - B*u )
119 * RESULT = -----------------------------------------
120 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
121 *
122 CALL DCOPY( N, D, 1, DF, 1 )
123 CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, DF, 1 )
124 *
125 CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, DF, 1 )
126 *
127 DNORM = DASUM( N, DF, 1 )
128 XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 )
129 YNORM = ANORM + BNORM
130 *
131 IF( XNORM.LE.ZERO ) THEN
132 RESULT = ZERO
133 ELSE
134 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
135 END IF
136 *
137 RETURN
138 *
139 * End of DGLMTS
140 *
141 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 * DGLMTS tests DGGGLM - 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) DOUBLE PRECISION array, dimension (LDA,M)
33 * The N-by-M matrix A.
34 *
35 * AF (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (LDB,P)
41 * The N-by-P matrix A.
42 *
43 * BF (workspace) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension( N )
49 * On input, the left hand side of the GLM.
50 *
51 * DF (workspace) DOUBLE PRECISION array, dimension( N )
52 *
53 * X (output) DOUBLE PRECISION array, dimension( M )
54 * solution vector X in the GLM problem.
55 *
56 * U (output) DOUBLE PRECISION array, dimension( P )
57 * solution vector U in the GLM problem.
58 *
59 * WORK (workspace) DOUBLE PRECISION 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 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
75 $ BF( LDB, * ), D( * ), DF( * ), RWORK( * ),
76 $ U( * ), WORK( LWORK ), X( * )
77 * ..
78 * .. Parameters ..
79 DOUBLE PRECISION ZERO, ONE
80 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
81 * ..
82 * .. Local Scalars ..
83 INTEGER INFO
84 DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM
85 * ..
86 * .. External Functions ..
87 DOUBLE PRECISION DASUM, DLAMCH, DLANGE
88 EXTERNAL DASUM, DLAMCH, DLANGE
89 * ..
90 * .. External Subroutines ..
91 *
92 EXTERNAL DCOPY, DGEMV, DGGGLM, DLACPY
93 * ..
94 * .. Intrinsic Functions ..
95 INTRINSIC MAX
96 * ..
97 * .. Executable Statements ..
98 *
99 EPS = DLAMCH( 'Epsilon' )
100 UNFL = DLAMCH( 'Safe minimum' )
101 ANORM = MAX( DLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
102 BNORM = MAX( DLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
103 *
104 * Copy the matrices A and B to the arrays AF and BF,
105 * and the vector D the array DF.
106 *
107 CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA )
108 CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB )
109 CALL DCOPY( N, D, 1, DF, 1 )
110 *
111 * Solve GLM problem
112 *
113 CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
114 $ INFO )
115 *
116 * Test the residual for the solution of LSE
117 *
118 * norm( d - A*x - B*u )
119 * RESULT = -----------------------------------------
120 * (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
121 *
122 CALL DCOPY( N, D, 1, DF, 1 )
123 CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, DF, 1 )
124 *
125 CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, DF, 1 )
126 *
127 DNORM = DASUM( N, DF, 1 )
128 XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 )
129 YNORM = ANORM + BNORM
130 *
131 IF( XNORM.LE.ZERO ) THEN
132 RESULT = ZERO
133 ELSE
134 RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS
135 END IF
136 *
137 RETURN
138 *
139 * End of DGLMTS
140 *
141 END