LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
ddrgsx.f
Go to the documentation of this file.
00001 *> \brief \b DDRGSX
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
00012 *                          BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
00013 *                          WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
00017 *      $                   NOUT, NSIZE
00018 *       DOUBLE PRECISION   THRESH
00019 *       ..
00020 *       .. Array Arguments ..
00021 *       LOGICAL            BWORK( * )
00022 *       INTEGER            IWORK( * )
00023 *       DOUBLE PRECISION   A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
00024 *      $                   ALPHAR( * ), B( LDA, * ), BETA( * ),
00025 *      $                   BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
00026 *      $                   WORK( * ), Z( LDA, * )
00027 *       ..
00028 *  
00029 *
00030 *> \par Purpose:
00031 *  =============
00032 *>
00033 *> \verbatim
00034 *>
00035 *> DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
00036 *> problem expert driver DGGESX.
00037 *>
00038 *> DGGESX factors A and B as Q S Z' and Q T Z', where ' means
00039 *> transpose, T is upper triangular, S is in generalized Schur form
00040 *> (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
00041 *> the 2x2 blocks corresponding to complex conjugate pairs of
00042 *> generalized eigenvalues), and Q and Z are orthogonal.  It also
00043 *> computes the generalized eigenvalues (alpha(1),beta(1)), ...,
00044 *> (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
00045 *> characteristic equation
00046 *>
00047 *>     det( A - w(j) B ) = 0
00048 *>
00049 *> Optionally it also reorders the eigenvalues so that a selected
00050 *> cluster of eigenvalues appears in the leading diagonal block of the
00051 *> Schur forms; computes a reciprocal condition number for the average
00052 *> of the selected eigenvalues; and computes a reciprocal condition
00053 *> number for the right and left deflating subspaces corresponding to
00054 *> the selected eigenvalues.
00055 *>
00056 *> When DDRGSX is called with NSIZE > 0, five (5) types of built-in
00057 *> matrix pairs are used to test the routine DGGESX.
00058 *>
00059 *> When DDRGSX is called with NSIZE = 0, it reads in test matrix data
00060 *> to test DGGESX.
00061 *>
00062 *> For each matrix pair, the following tests will be performed and
00063 *> compared with the threshhold THRESH except for the tests (7) and (9):
00064 *>
00065 *> (1)   | A - Q S Z' | / ( |A| n ulp )
00066 *>
00067 *> (2)   | B - Q T Z' | / ( |B| n ulp )
00068 *>
00069 *> (3)   | I - QQ' | / ( n ulp )
00070 *>
00071 *> (4)   | I - ZZ' | / ( n ulp )
00072 *>
00073 *> (5)   if A is in Schur form (i.e. quasi-triangular form)
00074 *>
00075 *> (6)   maximum over j of D(j)  where:
00076 *>
00077 *>       if alpha(j) is real:
00078 *>                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
00079 *>           D(j) = ------------------------ + -----------------------
00080 *>                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
00081 *>
00082 *>       if alpha(j) is complex:
00083 *>                                 | det( s S - w T ) |
00084 *>           D(j) = ---------------------------------------------------
00085 *>                  ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
00086 *>
00087 *>           and S and T are here the 2 x 2 diagonal blocks of S and T
00088 *>           corresponding to the j-th and j+1-th eigenvalues.
00089 *>
00090 *> (7)   if sorting worked and SDIM is the number of eigenvalues
00091 *>       which were selected.
00092 *>
00093 *> (8)   the estimated value DIF does not differ from the true values of
00094 *>       Difu and Difl more than a factor 10*THRESH. If the estimate DIF
00095 *>       equals zero the corresponding true values of Difu and Difl
00096 *>       should be less than EPS*norm(A, B). If the true value of Difu
00097 *>       and Difl equal zero, the estimate DIF should be less than
00098 *>       EPS*norm(A, B).
00099 *>
00100 *> (9)   If INFO = N+3 is returned by DGGESX, the reordering "failed"
00101 *>       and we check that DIF = PL = PR = 0 and that the true value of
00102 *>       Difu and Difl is < EPS*norm(A, B). We count the events when
00103 *>       INFO=N+3.
00104 *>
00105 *> For read-in test matrices, the above tests are run except that the
00106 *> exact value for DIF (and PL) is input data.  Additionally, there is
00107 *> one more test run for read-in test matrices:
00108 *>
00109 *> (10)  the estimated value PL does not differ from the true value of
00110 *>       PLTRU more than a factor THRESH. If the estimate PL equals
00111 *>       zero the corresponding true value of PLTRU should be less than
00112 *>       EPS*norm(A, B). If the true value of PLTRU equal zero, the
00113 *>       estimate PL should be less than EPS*norm(A, B).
00114 *>
00115 *> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
00116 *> matrix pairs are generated and tested. NSIZE should be kept small.
00117 *>
00118 *> SVD (routine DGESVD) is used for computing the true value of DIF_u
00119 *> and DIF_l when testing the built-in test problems.
00120 *>
00121 *> Built-in Test Matrices
00122 *> ======================
00123 *>
00124 *> All built-in test matrices are the 2 by 2 block of triangular
00125 *> matrices
00126 *>
00127 *>          A = [ A11 A12 ]    and      B = [ B11 B12 ]
00128 *>              [     A22 ]                 [     B22 ]
00129 *>
00130 *> where for different type of A11 and A22 are given as the following.
00131 *> A12 and B12 are chosen so that the generalized Sylvester equation
00132 *>
00133 *>          A11*R - L*A22 = -A12
00134 *>          B11*R - L*B22 = -B12
00135 *>
00136 *> have prescribed solution R and L.
00137 *>
00138 *> Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
00139 *>          B11 = I_m, B22 = I_k
00140 *>          where J_k(a,b) is the k-by-k Jordan block with ``a'' on
00141 *>          diagonal and ``b'' on superdiagonal.
00142 *>
00143 *> Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
00144 *>          B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
00145 *>          A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
00146 *>          B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
00147 *>
00148 *> Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
00149 *>          second diagonal block in A_11 and each third diagonal block
00150 *>          in A_22 are made as 2 by 2 blocks.
00151 *>
00152 *> Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
00153 *>             for i=1,...,m,  j=1,...,m and
00154 *>          A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
00155 *>             for i=m+1,...,k,  j=m+1,...,k
00156 *>
00157 *> Type 5:  (A,B) and have potentially close or common eigenvalues and
00158 *>          very large departure from block diagonality A_11 is chosen
00159 *>          as the m x m leading submatrix of A_1:
00160 *>                  |  1  b                            |
00161 *>                  | -b  1                            |
00162 *>                  |        1+d  b                    |
00163 *>                  |         -b 1+d                   |
00164 *>           A_1 =  |                  d  1            |
00165 *>                  |                 -1  d            |
00166 *>                  |                        -d  1     |
00167 *>                  |                        -1 -d     |
00168 *>                  |                               1  |
00169 *>          and A_22 is chosen as the k x k leading submatrix of A_2:
00170 *>                  | -1  b                            |
00171 *>                  | -b -1                            |
00172 *>                  |       1-d  b                     |
00173 *>                  |       -b  1-d                    |
00174 *>           A_2 =  |                 d 1+b            |
00175 *>                  |               -1-b d             |
00176 *>                  |                       -d  1+b    |
00177 *>                  |                      -1+b  -d    |
00178 *>                  |                              1-d |
00179 *>          and matrix B are chosen as identity matrices (see DLATM5).
00180 *>
00181 *> \endverbatim
00182 *
00183 *  Arguments:
00184 *  ==========
00185 *
00186 *> \param[in] NSIZE
00187 *> \verbatim
00188 *>          NSIZE is INTEGER
00189 *>          The maximum size of the matrices to use. NSIZE >= 0.
00190 *>          If NSIZE = 0, no built-in tests matrices are used, but
00191 *>          read-in test matrices are used to test DGGESX.
00192 *> \endverbatim
00193 *>
00194 *> \param[in] NCMAX
00195 *> \verbatim
00196 *>          NCMAX is INTEGER
00197 *>          Maximum allowable NMAX for generating Kroneker matrix
00198 *>          in call to DLAKF2
00199 *> \endverbatim
00200 *>
00201 *> \param[in] THRESH
00202 *> \verbatim
00203 *>          THRESH is DOUBLE PRECISION
00204 *>          A test will count as "failed" if the "error", computed as
00205 *>          described above, exceeds THRESH.  Note that the error
00206 *>          is scaled to be O(1), so THRESH should be a reasonably
00207 *>          small multiple of 1, e.g., 10 or 100.  In particular,
00208 *>          it should not depend on the precision (single vs. double)
00209 *>          or the size of the matrix.  THRESH >= 0.
00210 *> \endverbatim
00211 *>
00212 *> \param[in] NIN
00213 *> \verbatim
00214 *>          NIN is INTEGER
00215 *>          The FORTRAN unit number for reading in the data file of
00216 *>          problems to solve.
00217 *> \endverbatim
00218 *>
00219 *> \param[in] NOUT
00220 *> \verbatim
00221 *>          NOUT is INTEGER
00222 *>          The FORTRAN unit number for printing out error messages
00223 *>          (e.g., if a routine returns IINFO not equal to 0.)
00224 *> \endverbatim
00225 *>
00226 *> \param[out] A
00227 *> \verbatim
00228 *>          A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00229 *>          Used to store the matrix whose eigenvalues are to be
00230 *>          computed.  On exit, A contains the last matrix actually used.
00231 *> \endverbatim
00232 *>
00233 *> \param[in] LDA
00234 *> \verbatim
00235 *>          LDA is INTEGER
00236 *>          The leading dimension of A, B, AI, BI, Z and Q,
00237 *>          LDA >= max( 1, NSIZE ). For the read-in test,
00238 *>          LDA >= max( 1, N ), N is the size of the test matrices.
00239 *> \endverbatim
00240 *>
00241 *> \param[out] B
00242 *> \verbatim
00243 *>          B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00244 *>          Used to store the matrix whose eigenvalues are to be
00245 *>          computed.  On exit, B contains the last matrix actually used.
00246 *> \endverbatim
00247 *>
00248 *> \param[out] AI
00249 *> \verbatim
00250 *>          AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00251 *>          Copy of A, modified by DGGESX.
00252 *> \endverbatim
00253 *>
00254 *> \param[out] BI
00255 *> \verbatim
00256 *>          BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00257 *>          Copy of B, modified by DGGESX.
00258 *> \endverbatim
00259 *>
00260 *> \param[out] Z
00261 *> \verbatim
00262 *>          Z is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00263 *>          Z holds the left Schur vectors computed by DGGESX.
00264 *> \endverbatim
00265 *>
00266 *> \param[out] Q
00267 *> \verbatim
00268 *>          Q is DOUBLE PRECISION array, dimension (LDA, NSIZE)
00269 *>          Q holds the right Schur vectors computed by DGGESX.
00270 *> \endverbatim
00271 *>
00272 *> \param[out] ALPHAR
00273 *> \verbatim
00274 *>          ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
00275 *> \endverbatim
00276 *>
00277 *> \param[out] ALPHAI
00278 *> \verbatim
00279 *>          ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
00280 *> \endverbatim
00281 *>
00282 *> \param[out] BETA
00283 *> \verbatim
00284 *>          BETA is DOUBLE PRECISION array, dimension (NSIZE)
00285 *>
00286 *>          On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
00287 *> \endverbatim
00288 *>
00289 *> \param[out] C
00290 *> \verbatim
00291 *>          C is DOUBLE PRECISION array, dimension (LDC, LDC)
00292 *>          Store the matrix generated by subroutine DLAKF2, this is the
00293 *>          matrix formed by Kronecker products used for estimating
00294 *>          DIF.
00295 *> \endverbatim
00296 *>
00297 *> \param[in] LDC
00298 *> \verbatim
00299 *>          LDC is INTEGER
00300 *>          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
00301 *> \endverbatim
00302 *>
00303 *> \param[out] S
00304 *> \verbatim
00305 *>          S is DOUBLE PRECISION array, dimension (LDC)
00306 *>          Singular values of C
00307 *> \endverbatim
00308 *>
00309 *> \param[out] WORK
00310 *> \verbatim
00311 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
00312 *> \endverbatim
00313 *>
00314 *> \param[in] LWORK
00315 *> \verbatim
00316 *>          LWORK is INTEGER
00317 *>          The dimension of the array WORK.
00318 *>          LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
00319 *> \endverbatim
00320 *>
00321 *> \param[out] IWORK
00322 *> \verbatim
00323 *>          IWORK is INTEGER array, dimension (LIWORK)
00324 *> \endverbatim
00325 *>
00326 *> \param[in] LIWORK
00327 *> \verbatim
00328 *>          LIWORK is INTEGER
00329 *>          The dimension of the array IWORK. LIWORK >= NSIZE + 6.
00330 *> \endverbatim
00331 *>
00332 *> \param[out] BWORK
00333 *> \verbatim
00334 *>          BWORK is LOGICAL array, dimension (LDA)
00335 *> \endverbatim
00336 *>
00337 *> \param[out] INFO
00338 *> \verbatim
00339 *>          INFO is INTEGER
00340 *>          = 0:  successful exit
00341 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00342 *>          > 0:  A routine returned an error code.
00343 *> \endverbatim
00344 *
00345 *  Authors:
00346 *  ========
00347 *
00348 *> \author Univ. of Tennessee 
00349 *> \author Univ. of California Berkeley 
00350 *> \author Univ. of Colorado Denver 
00351 *> \author NAG Ltd. 
00352 *
00353 *> \date November 2011
00354 *
00355 *> \ingroup double_eig
00356 *
00357 *  =====================================================================
00358       SUBROUTINE DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
00359      $                   BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
00360      $                   WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
00361 *
00362 *  -- LAPACK test routine (version 3.4.0) --
00363 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00364 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00365 *     November 2011
00366 *
00367 *     .. Scalar Arguments ..
00368       INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
00369      $                   NOUT, NSIZE
00370       DOUBLE PRECISION   THRESH
00371 *     ..
00372 *     .. Array Arguments ..
00373       LOGICAL            BWORK( * )
00374       INTEGER            IWORK( * )
00375       DOUBLE PRECISION   A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
00376      $                   ALPHAR( * ), B( LDA, * ), BETA( * ),
00377      $                   BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
00378      $                   WORK( * ), Z( LDA, * )
00379 *     ..
00380 *
00381 *  =====================================================================
00382 *
00383 *     .. Parameters ..
00384       DOUBLE PRECISION   ZERO, ONE, TEN
00385       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 )
00386 *     ..
00387 *     .. Local Scalars ..
00388       LOGICAL            ILABAD
00389       CHARACTER          SENSE
00390       INTEGER            BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
00391      $                   MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
00392      $                   PRTYPE, QBA, QBB
00393       DOUBLE PRECISION   ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
00394      $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
00395 *     ..
00396 *     .. Local Arrays ..
00397       DOUBLE PRECISION   DIFEST( 2 ), PL( 2 ), RESULT( 10 )
00398 *     ..
00399 *     .. External Functions ..
00400       LOGICAL            DLCTSX
00401       INTEGER            ILAENV
00402       DOUBLE PRECISION   DLAMCH, DLANGE
00403       EXTERNAL           DLCTSX, ILAENV, DLAMCH, DLANGE
00404 *     ..
00405 *     .. External Subroutines ..
00406       EXTERNAL           ALASVM, DGESVD, DGET51, DGET53, DGGESX, DLABAD,
00407      $                   DLACPY, DLAKF2, DLASET, DLATM5, XERBLA
00408 *     ..
00409 *     .. Intrinsic Functions ..
00410       INTRINSIC          ABS, MAX, SQRT
00411 *     ..
00412 *     .. Scalars in Common ..
00413       LOGICAL            FS
00414       INTEGER            K, M, MPLUSN, N
00415 *     ..
00416 *     .. Common blocks ..
00417       COMMON             / MN / M, N, MPLUSN, K, FS
00418 *     ..
00419 *     .. Executable Statements ..
00420 *
00421 *     Check for errors
00422 *
00423       IF( NSIZE.LT.0 ) THEN
00424          INFO = -1
00425       ELSE IF( THRESH.LT.ZERO ) THEN
00426          INFO = -2
00427       ELSE IF( NIN.LE.0 ) THEN
00428          INFO = -3
00429       ELSE IF( NOUT.LE.0 ) THEN
00430          INFO = -4
00431       ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
00432          INFO = -6
00433       ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
00434          INFO = -17
00435       ELSE IF( LIWORK.LT.NSIZE+6 ) THEN
00436          INFO = -21
00437       END IF
00438 *
00439 *     Compute workspace
00440 *      (Note: Comments in the code beginning "Workspace:" describe the
00441 *       minimal amount of workspace needed at that point in the code,
00442 *       as well as the preferred amount for good performance.
00443 *       NB refers to the optimal block size for the immediately
00444 *       following subroutine, as returned by ILAENV.)
00445 *
00446       MINWRK = 1
00447       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00448          MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
00449 *
00450 *        workspace for sggesx
00451 *
00452          MAXWRK = 9*( NSIZE+1 ) + NSIZE*
00453      $            ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
00454          MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
00455      $            ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
00456 *
00457 *        workspace for dgesvd
00458 *
00459          BDSPAC = 5*NSIZE*NSIZE / 2
00460          MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
00461      $            ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2,
00462      $            NSIZE*NSIZE / 2, -1, -1 ) )
00463          MAXWRK = MAX( MAXWRK, BDSPAC )
00464 *
00465          MAXWRK = MAX( MAXWRK, MINWRK )
00466 *
00467          WORK( 1 ) = MAXWRK
00468       END IF
00469 *
00470       IF( LWORK.LT.MINWRK )
00471      $   INFO = -19
00472 *
00473       IF( INFO.NE.0 ) THEN
00474          CALL XERBLA( 'DDRGSX', -INFO )
00475          RETURN
00476       END IF
00477 *
00478 *     Important constants
00479 *
00480       ULP = DLAMCH( 'P' )
00481       ULPINV = ONE / ULP
00482       SMLNUM = DLAMCH( 'S' ) / ULP
00483       BIGNUM = ONE / SMLNUM
00484       CALL DLABAD( SMLNUM, BIGNUM )
00485       THRSH2 = TEN*THRESH
00486       NTESTT = 0
00487       NERRS = 0
00488 *
00489 *     Go to the tests for read-in matrix pairs
00490 *
00491       IFUNC = 0
00492       IF( NSIZE.EQ.0 )
00493      $   GO TO 70
00494 *
00495 *     Test the built-in matrix pairs.
00496 *     Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
00497 *     of test matrices, different size (M+N)
00498 *
00499       PRTYPE = 0
00500       QBA = 3
00501       QBB = 4
00502       WEIGHT = SQRT( ULP )
00503 *
00504       DO 60 IFUNC = 0, 3
00505          DO 50 PRTYPE = 1, 5
00506             DO 40 M = 1, NSIZE - 1
00507                DO 30 N = 1, NSIZE - M
00508 *
00509                   WEIGHT = ONE / WEIGHT
00510                   MPLUSN = M + N
00511 *
00512 *                 Generate test matrices
00513 *
00514                   FS = .TRUE.
00515                   K = 0
00516 *
00517                   CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
00518      $                         LDA )
00519                   CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
00520      $                         LDA )
00521 *
00522                   CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
00523      $                         LDA, AI( 1, M+1 ), LDA, BI, LDA,
00524      $                         BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
00525      $                         Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
00526 *
00527 *                 Compute the Schur factorization and swapping the
00528 *                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00529 *                 Swapping is accomplished via the function DLCTSX
00530 *                 which is supplied below.
00531 *
00532                   IF( IFUNC.EQ.0 ) THEN
00533                      SENSE = 'N'
00534                   ELSE IF( IFUNC.EQ.1 ) THEN
00535                      SENSE = 'E'
00536                   ELSE IF( IFUNC.EQ.2 ) THEN
00537                      SENSE = 'V'
00538                   ELSE IF( IFUNC.EQ.3 ) THEN
00539                      SENSE = 'B'
00540                   END IF
00541 *
00542                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00543                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00544 *
00545                   CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI,
00546      $                         LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
00547      $                         Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
00548      $                         IWORK, LIWORK, BWORK, LINFO )
00549 *
00550                   IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00551                      RESULT( 1 ) = ULPINV
00552                      WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN,
00553      $                  PRTYPE
00554                      INFO = LINFO
00555                      GO TO 30
00556                   END IF
00557 *
00558 *                 Compute the norm(A, B)
00559 *
00560                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
00561      $                         MPLUSN )
00562                   CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00563      $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00564                   ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
00565      $                    WORK )
00566 *
00567 *                 Do tests (1) to (4)
00568 *
00569                   CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
00570      $                         LDA, WORK, RESULT( 1 ) )
00571                   CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
00572      $                         LDA, WORK, RESULT( 2 ) )
00573                   CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
00574      $                         LDA, WORK, RESULT( 3 ) )
00575                   CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
00576      $                         LDA, WORK, RESULT( 4 ) )
00577                   NTEST = 4
00578 *
00579 *                 Do tests (5) and (6): check Schur form of A and
00580 *                 compare eigenvalues with diagonals.
00581 *
00582                   TEMP1 = ZERO
00583                   RESULT( 5 ) = ZERO
00584                   RESULT( 6 ) = ZERO
00585 *
00586                   DO 10 J = 1, MPLUSN
00587                      ILABAD = .FALSE.
00588                      IF( ALPHAI( J ).EQ.ZERO ) THEN
00589                         TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00590      $                          MAX( SMLNUM, ABS( ALPHAR( J ) ),
00591      $                          ABS( AI( J, J ) ) )+
00592      $                          ABS( BETA( J )-BI( J, J ) ) /
00593      $                          MAX( SMLNUM, ABS( BETA( J ) ),
00594      $                          ABS( BI( J, J ) ) ) ) / ULP
00595                         IF( J.LT.MPLUSN ) THEN
00596                            IF( AI( J+1, J ).NE.ZERO ) THEN
00597                               ILABAD = .TRUE.
00598                               RESULT( 5 ) = ULPINV
00599                            END IF
00600                         END IF
00601                         IF( J.GT.1 ) THEN
00602                            IF( AI( J, J-1 ).NE.ZERO ) THEN
00603                               ILABAD = .TRUE.
00604                               RESULT( 5 ) = ULPINV
00605                            END IF
00606                         END IF
00607                      ELSE
00608                         IF( ALPHAI( J ).GT.ZERO ) THEN
00609                            I1 = J
00610                         ELSE
00611                            I1 = J - 1
00612                         END IF
00613                         IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00614                            ILABAD = .TRUE.
00615                         ELSE IF( I1.LT.MPLUSN-1 ) THEN
00616                            IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00617                               ILABAD = .TRUE.
00618                               RESULT( 5 ) = ULPINV
00619                            END IF
00620                         ELSE IF( I1.GT.1 ) THEN
00621                            IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00622                               ILABAD = .TRUE.
00623                               RESULT( 5 ) = ULPINV
00624                            END IF
00625                         END IF
00626                         IF( .NOT.ILABAD ) THEN
00627                            CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
00628      $                                  LDA, BETA( J ), ALPHAR( J ),
00629      $                                  ALPHAI( J ), TEMP2, IINFO )
00630                            IF( IINFO.GE.3 ) THEN
00631                               WRITE( NOUT, FMT = 9997 )IINFO, J,
00632      $                           MPLUSN, PRTYPE
00633                               INFO = ABS( IINFO )
00634                            END IF
00635                         ELSE
00636                            TEMP2 = ULPINV
00637                         END IF
00638                      END IF
00639                      TEMP1 = MAX( TEMP1, TEMP2 )
00640                      IF( ILABAD ) THEN
00641                         WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
00642                      END IF
00643    10             CONTINUE
00644                   RESULT( 6 ) = TEMP1
00645                   NTEST = NTEST + 2
00646 *
00647 *                 Test (7) (if sorting worked)
00648 *
00649                   RESULT( 7 ) = ZERO
00650                   IF( LINFO.EQ.MPLUSN+3 ) THEN
00651                      RESULT( 7 ) = ULPINV
00652                   ELSE IF( MM.NE.N ) THEN
00653                      RESULT( 7 ) = ULPINV
00654                   END IF
00655                   NTEST = NTEST + 1
00656 *
00657 *                 Test (8): compare the estimated value DIF and its
00658 *                 value. first, compute the exact DIF.
00659 *
00660                   RESULT( 8 ) = ZERO
00661                   MN2 = MM*( MPLUSN-MM )*2
00662                   IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
00663 *
00664 *                    Note: for either following two causes, there are
00665 *                    almost same number of test cases fail the test.
00666 *
00667                      CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
00668      $                            AI( MM+1, MM+1 ), BI,
00669      $                            BI( MM+1, MM+1 ), C, LDC )
00670 *
00671                      CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
00672      $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
00673      $                            INFO )
00674                      DIFTRU = S( MN2 )
00675 *
00676                      IF( DIFEST( 2 ).EQ.ZERO ) THEN
00677                         IF( DIFTRU.GT.ABNRM*ULP )
00678      $                     RESULT( 8 ) = ULPINV
00679                      ELSE IF( DIFTRU.EQ.ZERO ) THEN
00680                         IF( DIFEST( 2 ).GT.ABNRM*ULP )
00681      $                     RESULT( 8 ) = ULPINV
00682                      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00683      $                        ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00684                         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
00685      $                                DIFEST( 2 ) / DIFTRU )
00686                      END IF
00687                      NTEST = NTEST + 1
00688                   END IF
00689 *
00690 *                 Test (9)
00691 *
00692                   RESULT( 9 ) = ZERO
00693                   IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00694                      IF( DIFTRU.GT.ABNRM*ULP )
00695      $                  RESULT( 9 ) = ULPINV
00696                      IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00697      $                  RESULT( 9 ) = ULPINV
00698                      IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00699      $                  RESULT( 9 ) = ULPINV
00700                      NTEST = NTEST + 1
00701                   END IF
00702 *
00703                   NTESTT = NTESTT + NTEST
00704 *
00705 *                 Print out tests which fail.
00706 *
00707                   DO 20 J = 1, 9
00708                      IF( RESULT( J ).GE.THRESH ) THEN
00709 *
00710 *                       If this is the first test to fail,
00711 *                       print a header to the data file.
00712 *
00713                         IF( NERRS.EQ.0 ) THEN
00714                            WRITE( NOUT, FMT = 9995 )'DGX'
00715 *
00716 *                          Matrix types
00717 *
00718                            WRITE( NOUT, FMT = 9993 )
00719 *
00720 *                          Tests performed
00721 *
00722                            WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00723      $                        'transpose', ( '''', I = 1, 4 )
00724 *
00725                         END IF
00726                         NERRS = NERRS + 1
00727                         IF( RESULT( J ).LT.10000.0D0 ) THEN
00728                            WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
00729      $                        WEIGHT, M, J, RESULT( J )
00730                         ELSE
00731                            WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
00732      $                        WEIGHT, M, J, RESULT( J )
00733                         END IF
00734                      END IF
00735    20             CONTINUE
00736 *
00737    30          CONTINUE
00738    40       CONTINUE
00739    50    CONTINUE
00740    60 CONTINUE
00741 *
00742       GO TO 150
00743 *
00744    70 CONTINUE
00745 *
00746 *     Read in data from file to check accuracy of condition estimation
00747 *     Read input data until N=0
00748 *
00749       NPTKNT = 0
00750 *
00751    80 CONTINUE
00752       READ( NIN, FMT = *, END = 140 )MPLUSN
00753       IF( MPLUSN.EQ.0 )
00754      $   GO TO 140
00755       READ( NIN, FMT = *, END = 140 )N
00756       DO 90 I = 1, MPLUSN
00757          READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
00758    90 CONTINUE
00759       DO 100 I = 1, MPLUSN
00760          READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
00761   100 CONTINUE
00762       READ( NIN, FMT = * )PLTRU, DIFTRU
00763 *
00764       NPTKNT = NPTKNT + 1
00765       FS = .TRUE.
00766       K = 0
00767       M = MPLUSN - N
00768 *
00769       CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00770       CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00771 *
00772 *     Compute the Schur factorization while swaping the
00773 *     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00774 *
00775       CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
00776      $             MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
00777      $             WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
00778 *
00779       IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00780          RESULT( 1 ) = ULPINV
00781          WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
00782          GO TO 130
00783       END IF
00784 *
00785 *     Compute the norm(A, B)
00786 *        (should this be norm of (A,B) or (AI,BI)?)
00787 *
00788       CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
00789       CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00790      $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00791       ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
00792 *
00793 *     Do tests (1) to (4)
00794 *
00795       CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
00796      $             RESULT( 1 ) )
00797       CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
00798      $             RESULT( 2 ) )
00799       CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
00800      $             RESULT( 3 ) )
00801       CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
00802      $             RESULT( 4 ) )
00803 *
00804 *     Do tests (5) and (6): check Schur form of A and compare
00805 *     eigenvalues with diagonals.
00806 *
00807       NTEST = 6
00808       TEMP1 = ZERO
00809       RESULT( 5 ) = ZERO
00810       RESULT( 6 ) = ZERO
00811 *
00812       DO 110 J = 1, MPLUSN
00813          ILABAD = .FALSE.
00814          IF( ALPHAI( J ).EQ.ZERO ) THEN
00815             TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
00816      $              MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
00817      $              J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
00818      $              MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
00819      $               / ULP
00820             IF( J.LT.MPLUSN ) THEN
00821                IF( AI( J+1, J ).NE.ZERO ) THEN
00822                   ILABAD = .TRUE.
00823                   RESULT( 5 ) = ULPINV
00824                END IF
00825             END IF
00826             IF( J.GT.1 ) THEN
00827                IF( AI( J, J-1 ).NE.ZERO ) THEN
00828                   ILABAD = .TRUE.
00829                   RESULT( 5 ) = ULPINV
00830                END IF
00831             END IF
00832          ELSE
00833             IF( ALPHAI( J ).GT.ZERO ) THEN
00834                I1 = J
00835             ELSE
00836                I1 = J - 1
00837             END IF
00838             IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
00839                ILABAD = .TRUE.
00840             ELSE IF( I1.LT.MPLUSN-1 ) THEN
00841                IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
00842                   ILABAD = .TRUE.
00843                   RESULT( 5 ) = ULPINV
00844                END IF
00845             ELSE IF( I1.GT.1 ) THEN
00846                IF( AI( I1, I1-1 ).NE.ZERO ) THEN
00847                   ILABAD = .TRUE.
00848                   RESULT( 5 ) = ULPINV
00849                END IF
00850             END IF
00851             IF( .NOT.ILABAD ) THEN
00852                CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
00853      $                      BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
00854      $                      IINFO )
00855                IF( IINFO.GE.3 ) THEN
00856                   WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
00857                   INFO = ABS( IINFO )
00858                END IF
00859             ELSE
00860                TEMP2 = ULPINV
00861             END IF
00862          END IF
00863          TEMP1 = MAX( TEMP1, TEMP2 )
00864          IF( ILABAD ) THEN
00865             WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
00866          END IF
00867   110 CONTINUE
00868       RESULT( 6 ) = TEMP1
00869 *
00870 *     Test (7) (if sorting worked)  <--------- need to be checked.
00871 *
00872       NTEST = 7
00873       RESULT( 7 ) = ZERO
00874       IF( LINFO.EQ.MPLUSN+3 )
00875      $   RESULT( 7 ) = ULPINV
00876 *
00877 *     Test (8): compare the estimated value of DIF and its true value.
00878 *
00879       NTEST = 8
00880       RESULT( 8 ) = ZERO
00881       IF( DIFEST( 2 ).EQ.ZERO ) THEN
00882          IF( DIFTRU.GT.ABNRM*ULP )
00883      $      RESULT( 8 ) = ULPINV
00884       ELSE IF( DIFTRU.EQ.ZERO ) THEN
00885          IF( DIFEST( 2 ).GT.ABNRM*ULP )
00886      $      RESULT( 8 ) = ULPINV
00887       ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00888      $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00889          RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
00890       END IF
00891 *
00892 *     Test (9)
00893 *
00894       NTEST = 9
00895       RESULT( 9 ) = ZERO
00896       IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00897          IF( DIFTRU.GT.ABNRM*ULP )
00898      $      RESULT( 9 ) = ULPINV
00899          IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00900      $      RESULT( 9 ) = ULPINV
00901          IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00902      $      RESULT( 9 ) = ULPINV
00903       END IF
00904 *
00905 *     Test (10): compare the estimated value of PL and it true value.
00906 *
00907       NTEST = 10
00908       RESULT( 10 ) = ZERO
00909       IF( PL( 1 ).EQ.ZERO ) THEN
00910          IF( PLTRU.GT.ABNRM*ULP )
00911      $      RESULT( 10 ) = ULPINV
00912       ELSE IF( PLTRU.EQ.ZERO ) THEN
00913          IF( PL( 1 ).GT.ABNRM*ULP )
00914      $      RESULT( 10 ) = ULPINV
00915       ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
00916      $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
00917          RESULT( 10 ) = ULPINV
00918       END IF
00919 *
00920       NTESTT = NTESTT + NTEST
00921 *
00922 *     Print out tests which fail.
00923 *
00924       DO 120 J = 1, NTEST
00925          IF( RESULT( J ).GE.THRESH ) THEN
00926 *
00927 *           If this is the first test to fail,
00928 *           print a header to the data file.
00929 *
00930             IF( NERRS.EQ.0 ) THEN
00931                WRITE( NOUT, FMT = 9995 )'DGX'
00932 *
00933 *              Matrix types
00934 *
00935                WRITE( NOUT, FMT = 9994 )
00936 *
00937 *              Tests performed
00938 *
00939                WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
00940      $            'transpose', ( '''', I = 1, 4 )
00941 *
00942             END IF
00943             NERRS = NERRS + 1
00944             IF( RESULT( J ).LT.10000.0D0 ) THEN
00945                WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
00946             ELSE
00947                WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
00948             END IF
00949          END IF
00950 *
00951   120 CONTINUE
00952 *
00953   130 CONTINUE
00954       GO TO 80
00955   140 CONTINUE
00956 *
00957   150 CONTINUE
00958 *
00959 *     Summary
00960 *
00961       CALL ALASVM( 'DGX', NOUT, NERRS, NTESTT, 0 )
00962 *
00963       WORK( 1 ) = MAXWRK
00964 *
00965       RETURN
00966 *
00967  9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00968      $      I6, ', JTYPE=', I6, ')' )
00969 *
00970  9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00971      $      I6, ', Input Example #', I2, ')' )
00972 *
00973  9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
00974      $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00975 *
00976  9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
00977      $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00978 *
00979  9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
00980      $      ' problem driver' )
00981 *
00982  9994 FORMAT( 'Input Example' )
00983 *
00984  9993 FORMAT( ' Matrix types: ', /
00985      $      '  1:  A is a block diagonal matrix of Jordan blocks ',
00986      $      'and B is the identity ', / '      matrix, ',
00987      $      / '  2:  A and B are upper triangular matrices, ',
00988      $      / '  3:  A and B are as type 2, but each second diagonal ',
00989      $      'block in A_11 and ', /
00990      $      '      each third diaongal block in A_22 are 2x2 blocks,',
00991      $      / '  4:  A and B are block diagonal matrices, ',
00992      $      / '  5:  (A,B) has potentially close or common ',
00993      $      'eigenvalues.', / )
00994 *
00995  9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00996      $      'Q and Z are ', A, ',', / 19X,
00997      $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
00998      $      / '  1 = | A - Q S Z', A,
00999      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
01000      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
01001      $      ' | / ( n ulp )             4 = | I - ZZ', A,
01002      $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
01003      $      'Schur form S', / '  6 = difference between (alpha,beta)',
01004      $      ' and diagonals of (S,T)', /
01005      $      '  7 = 1/ULP  if SDIM is not the correct number of ',
01006      $      'selected eigenvalues', /
01007      $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
01008      $      'DIFTRU/DIFEST > 10*THRESH',
01009      $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
01010      $      'when reordering fails', /
01011      $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
01012      $      'PLTRU/PLEST > THRESH', /
01013      $      '    ( Test 10 is only for input examples )', / )
01014  9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3,
01015      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
01016  9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3,
01017      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.3 )
01018  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
01019      $      ' result ', I2, ' is', 0P, F8.2 )
01020  9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
01021      $      ' result ', I2, ' is', 1P, D10.3 )
01022 *
01023 *     End of DDRGSX
01024 *
01025       END
 All Files Functions