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