LAPACK 3.3.0

# sglmts.f

Go to the documentation of this file.
```00001       SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF,
00002      \$                   X, U, WORK, LWORK, RWORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            LDA, LDB, LWORK, M, P, N
00010       REAL               RESULT
00011 *     ..
00012 *     .. Array Arguments ..
00013       REAL               A( LDA, * ), AF( LDA, * ), B( LDB, * ),
00014      \$                   BF( LDB, * ), RWORK( * ), D( * ), DF( * ),
00015      \$                   U( * ), WORK( LWORK ), X( * )
00016 *
00017 *  Purpose
00018 *  =======
00019 *
00020 *  SGLMTS tests SGGGLM - a subroutine for solving the generalized
00021 *  linear model problem.
00022 *
00023 *  Arguments
00024 *  =========
00025 *
00026 *  N       (input) INTEGER
00027 *          The number of rows of the matrices A and B.  N >= 0.
00028 *
00029 *  M       (input) INTEGER
00030 *          The number of columns of the matrix A.  M >= 0.
00031 *
00032 *  P       (input) INTEGER
00033 *          The number of columns of the matrix B.  P >= 0.
00034 *
00035 *  A       (input) REAL array, dimension (LDA,M)
00036 *          The N-by-M matrix A.
00037 *
00038 *  AF      (workspace) REAL array, dimension (LDA,M)
00039 *
00040 *  LDA     (input) INTEGER
00041 *          The leading dimension of the arrays A, AF. LDA >= max(M,N).
00042 *
00043 *  B       (input) REAL array, dimension (LDB,P)
00044 *          The N-by-P matrix A.
00045 *
00046 *  BF      (workspace) REAL array, dimension (LDB,P)
00047 *
00048 *  LDB     (input) INTEGER
00049 *          The leading dimension of the arrays B, BF. LDB >= max(P,N).
00050 *
00051 *  D       (input) REAL array, dimension( N )
00052 *          On input, the left hand side of the GLM.
00053 *
00054 *  DF      (workspace) REAL array, dimension( N )
00055 *
00056 *  X       (output) REAL array, dimension( M )
00057 *          solution vector X in the GLM problem.
00058 *
00059 *  U       (output) REAL array, dimension( P )
00060 *          solution vector U in the GLM problem.
00061 *
00062 *  WORK    (workspace) REAL array, dimension (LWORK)
00063 *
00064 *  LWORK   (input) INTEGER
00065 *          The dimension of the array WORK.
00066 *
00067 *  RWORK   (workspace) REAL array, dimension (M)
00068 *
00069 *  RESULT   (output) REAL
00070 *          The test ratio:
00071 *                           norm( d - A*x - B*u )
00072 *            RESULT = -----------------------------------------
00073 *                     (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
00074 *
00075 *  ====================================================================
00076 *
00077 *     .. Parameters ..
00078       REAL               ZERO, ONE
00079       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00080 *     ..
00081 *     .. Local Scalars ..
00082       INTEGER            INFO
00083       REAL               ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL
00084 *     ..
00085 *     .. External Functions ..
00086       REAL               SASUM, SLAMCH, SLANGE
00087       EXTERNAL           SASUM, SLAMCH, SLANGE
00088 *     ..
00089 *     .. External Subroutines ..
00090       EXTERNAL           SLACPY
00091 *
00092 *     .. Intrinsic Functions ..
00093       INTRINSIC          MAX
00094 *     ..
00095 *     .. Executable Statements ..
00096 *
00097       EPS = SLAMCH( 'Epsilon' )
00098       UNFL = SLAMCH( 'Safe minimum' )
00099       ANORM = MAX( SLANGE( '1', N, M, A, LDA, RWORK ), UNFL )
00100       BNORM = MAX( SLANGE( '1', N, P, B, LDB, RWORK ), UNFL )
00101 *
00102 *     Copy the matrices A and B to the arrays AF and BF,
00103 *     and the vector D the array DF.
00104 *
00105       CALL SLACPY( 'Full', N, M, A, LDA, AF, LDA )
00106       CALL SLACPY( 'Full', N, P, B, LDB, BF, LDB )
00107       CALL SCOPY( N, D, 1, DF, 1 )
00108 *
00109 *     Solve GLM problem
00110 *
00111       CALL SGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, LWORK,
00112      \$             INFO )
00113 *
00114 *     Test the residual for the solution of LSE
00115 *
00116 *                       norm( d - A*x - B*u )
00117 *       RESULT = -----------------------------------------
00118 *                (norm(A)+norm(B))*(norm(x)+norm(u))*EPS
00119 *
00120       CALL SCOPY( N, D, 1, DF, 1 )
00121       CALL SGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1,
00122      \$             ONE, DF, 1 )
00123 *
00124       CALL SGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1,
00125      \$             ONE, DF, 1 )
00126 *
00127       DNORM = SASUM( N, DF, 1 )
00128       XNORM = SASUM( M, X, 1 ) + SASUM( P, U, 1 )
00129       YNORM = ANORM + BNORM
00130 *
00131       IF( XNORM.LE.ZERO ) THEN
00132          RESULT = ZERO
00133       ELSE
00134          RESULT =  ( ( DNORM / YNORM ) / XNORM ) /EPS
00135       END IF
00136 *
00137       RETURN
00138 *
00139 *     End of SGLMTS
00140 *
00141       END
```