1       SUBROUTINE DCHKQP( DOTYPE, NM, MVAL, NN, NVAL, THRESH, TSTERR, A,
  2      $                   COPYA, S, COPYS, TAU, WORK, IWORK, NOUT )
  3 *
  4 *  -- LAPACK test routine (version 3.1.1) --
  5 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  6 *     January 2007
  7 *
  8 *     .. Scalar Arguments ..
  9       LOGICAL            TSTERR
 10       INTEGER            NM, NN, NOUT
 11       DOUBLE PRECISION   THRESH
 12 *     ..
 13 *     .. Array Arguments ..
 14       LOGICAL            DOTYPE( * )
 15       INTEGER            IWORK( * ), MVAL( * ), NVAL( * )
 16       DOUBLE PRECISION   A( * ), COPYA( * ), COPYS( * ), S( * ),
 17      $                   TAU( * ), WORK( * )
 18 *     ..
 19 *
 20 *  Purpose
 21 *  =======
 22 *
 23 *  DCHKQP tests DGEQPF.
 24 *
 25 *  Arguments
 26 *  =========
 27 *
 28 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
 29 *          The matrix types to be used for testing.  Matrices of type j
 30 *          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
 31 *          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
 32 *
 33 *  NM      (input) INTEGER
 34 *          The number of values of M contained in the vector MVAL.
 35 *
 36 *  MVAL    (input) INTEGER array, dimension (NM)
 37 *          The values of the matrix row dimension M.
 38 *
 39 *  NN      (input) INTEGER
 40 *          The number of values of N contained in the vector NVAL.
 41 *
 42 *  NVAL    (input) INTEGER array, dimension (NN)
 43 *          The values of the matrix column dimension N.
 44 *
 45 *  THRESH  (input) DOUBLE PRECISION
 46 *          The threshold value for the test ratios.  A result is
 47 *          included in the output file if RESULT >= THRESH.  To have
 48 *          every test ratio printed, use THRESH = 0.
 49 *
 50 *  TSTERR  (input) LOGICAL
 51 *          Flag that indicates whether error exits are to be tested.
 52 *
 53 *  A       (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
 54 *          where MMAX is the maximum value of M in MVAL and NMAX is the
 55 *          maximum value of N in NVAL.
 56 *
 57 *  COPYA   (workspace) DOUBLE PRECISION array, dimension (MMAX*NMAX)
 58 *
 59 *  S       (workspace) DOUBLE PRECISION array, dimension
 60 *                      (min(MMAX,NMAX))
 61 *
 62 *  COPYS   (workspace) DOUBLE PRECISION array, dimension
 63 *                      (min(MMAX,NMAX))
 64 *
 65 *  TAU     (workspace) DOUBLE PRECISION array, dimension (MMAX)
 66 *
 67 *  WORK    (workspace) DOUBLE PRECISION array, dimension
 68 *                      (MMAX*NMAX + 4*NMAX + MMAX)
 69 *
 70 *  IWORK   (workspace) INTEGER array, dimension (NMAX)
 71 *
 72 *  NOUT    (input) INTEGER
 73 *          The unit number for output.
 74 *
 75 *  =====================================================================
 76 *
 77 *     .. Parameters ..
 78       INTEGER            NTYPES
 79       PARAMETER          ( NTYPES = 6 )
 80       INTEGER            NTESTS
 81       PARAMETER          ( NTESTS = 3 )
 82       DOUBLE PRECISION   ONE, ZERO
 83       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
 84 *     ..
 85 *     .. Local Scalars ..
 86       CHARACTER*3        PATH
 87       INTEGER            I, IHIGH, ILOW, IM, IMODE, IN, INFO, ISTEP, K,
 88      $                   LDA, LWORK, M, MNMIN, MODE, N, NERRS, NFAIL,
 89      $                   NRUN
 90       DOUBLE PRECISION   EPS
 91 *     ..
 92 *     .. Local Arrays ..
 93       INTEGER            ISEED( 4 ), ISEEDY( 4 )
 94       DOUBLE PRECISION   RESULT( NTESTS )
 95 *     ..
 96 *     .. External Functions ..
 97       DOUBLE PRECISION   DLAMCH, DQPT01, DQRT11, DQRT12
 98       EXTERNAL           DLAMCH, DQPT01, DQRT11, DQRT12
 99 *     ..
100 *     .. External Subroutines ..
101       EXTERNAL           ALAHD, ALASUM, DERRQP, DGEQPF, DLACPY, DLAORD,
102      $                   DLASET, DLATMS
103 *     ..
104 *     .. Intrinsic Functions ..
105       INTRINSIC          MAXMIN
106 *     ..
107 *     .. Scalars in Common ..
108       LOGICAL            LERR, OK
109       CHARACTER*32       SRNAMT
110       INTEGER            INFOT, IOUNIT
111 *     ..
112 *     .. Common blocks ..
113       COMMON             / INFOC / INFOT, IOUNIT, OK, LERR
114       COMMON             / SRNAMC / SRNAMT
115 *     ..
116 *     .. Data statements ..
117       DATA               ISEEDY / 1988198919901991 /
118 *     ..
119 *     .. Executable Statements ..
120 *
121 *     Initialize constants and the random number seed.
122 *
123       PATH( 11 ) = 'Double precision'
124       PATH( 23 ) = 'QP'
125       NRUN = 0
126       NFAIL = 0
127       NERRS = 0
128       DO 10 I = 14
129          ISEED( I ) = ISEEDY( I )
130    10 CONTINUE
131       EPS = DLAMCH( 'Epsilon' )
132 *
133 *     Test the error exits
134 *
135       IF( TSTERR )
136      $   CALL DERRQP( PATH, NOUT )
137       INFOT = 0
138 *
139       DO 80 IM = 1, NM
140 *
141 *        Do for each value of M in MVAL.
142 *
143          M = MVAL( IM )
144          LDA = MAX1, M )
145 *
146          DO 70 IN = 1, NN
147 *
148 *           Do for each value of N in NVAL.
149 *
150             N = NVAL( IN )
151             MNMIN = MIN( M, N )
152             LWORK = MAX1, M*MAX( M, N ) + 4*MNMIN + MAX( M, N ),
153      $                   M*+ 2*MNMIN + 4*N )
154 *
155             DO 60 IMODE = 1, NTYPES
156                IF.NOT.DOTYPE( IMODE ) )
157      $            GO TO 60
158 *
159 *              Do for each type of matrix
160 *                 1:  zero matrix
161 *                 2:  one small singular value
162 *                 3:  geometric distribution of singular values
163 *                 4:  first n/2 columns fixed
164 *                 5:  last n/2 columns fixed
165 *                 6:  every second column fixed
166 *
167                MODE = IMODE
168                IF( IMODE.GT.3 )
169      $            MODE = 1
170 *
171 *              Generate test matrix of size m by n using
172 *              singular value distribution indicated by `mode'.
173 *
174                DO 20 I = 1, N
175                   IWORK( I ) = 0
176    20          CONTINUE
177                IF( IMODE.EQ.1 ) THEN
178                   CALL DLASET( 'Full', M, N, ZERO, ZERO, COPYA, LDA )
179                   DO 30 I = 1, MNMIN
180                      COPYS( I ) = ZERO
181    30             CONTINUE
182                ELSE
183                   CALL DLATMS( M, N, 'Uniform', ISEED, 'Nonsymm', COPYS,
184      $                         MODE, ONE / EPS, ONE, M, N, 'No packing',
185      $                         COPYA, LDA, WORK, INFO )
186                   IF( IMODE.GE.4 ) THEN
187                      IF( IMODE.EQ.4 ) THEN
188                         ILOW = 1
189                         ISTEP = 1
190                         IHIGH = MAX1, N / 2 )
191                      ELSE IF( IMODE.EQ.5 ) THEN
192                         ILOW = MAX1, N / 2 )
193                         ISTEP = 1
194                         IHIGH = N
195                      ELSE IF( IMODE.EQ.6 ) THEN
196                         ILOW = 1
197                         ISTEP = 2
198                         IHIGH = N
199                      END IF
200                      DO 40 I = ILOW, IHIGH, ISTEP
201                         IWORK( I ) = 1
202    40                CONTINUE
203                   END IF
204                   CALL DLAORD( 'Decreasing', MNMIN, COPYS, 1 )
205                END IF
206 *
207 *              Save A and its singular values
208 *
209                CALL DLACPY( 'All', M, N, COPYA, LDA, A, LDA )
210 *
211 *              Compute the QR factorization with pivoting of A
212 *
213                SRNAMT = 'DGEQPF'
214                CALL DGEQPF( M, N, A, LDA, IWORK, TAU, WORK, INFO )
215 *
216 *              Compute norm(svd(a) - svd(r))
217 *
218                RESULT1 ) = DQRT12( M, N, A, LDA, COPYS, WORK, LWORK )
219 *
220 *              Compute norm( A*P - Q*R )
221 *
222                RESULT2 ) = DQPT01( M, N, MNMIN, COPYA, A, LDA, TAU,
223      $                       IWORK, WORK, LWORK )
224 *
225 *              Compute Q'*Q
226 *
227                RESULT3 ) = DQRT11( M, MNMIN, A, LDA, TAU, WORK,
228      $                       LWORK )
229 *
230 *              Print information about the tests that did not pass
231 *              the threshold.
232 *
233                DO 50 K = 13
234                   IFRESULT( K ).GE.THRESH ) THEN
235                      IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
236      $                  CALL ALAHD( NOUT, PATH )
237                      WRITE( NOUT, FMT = 9999 )M, N, IMODE, K,
238      $                  RESULT( K )
239                      NFAIL = NFAIL + 1
240                   END IF
241    50          CONTINUE
242                NRUN = NRUN + 3
243    60       CONTINUE
244    70    CONTINUE
245    80 CONTINUE
246 *
247 *     Print a summary of the results.
248 *
249       CALL ALASUM( PATH, NOUT, NFAIL, NRUN, NERRS )
250 *
251  9999 FORMAT' M =', I5, ', N =', I5, ', type ', I2, ', test ', I2,
252      $      ', ratio ='G12.5 )
253 *
254 *     End of DCHKQP
255 *
256       END