LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
schkqp3rk.f
Go to the documentation of this file.
1*> \brief \b SCHKQP3RK
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE SCHKQP3RK( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
12* $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
13* $ B, COPYB, S, TAU,
14* $ WORK, IWORK, NOUT )
15* IMPLICIT NONE
16*
17* -- LAPACK test routine --
18* -- LAPACK is a software package provided by Univ. of Tennessee, --
19* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
20*
21* .. Scalar Arguments ..
22* INTEGER NM, NN, NNS, NNB, NOUT
23* REAL THRESH
24* ..
25* .. Array Arguments ..
26* LOGICAL DOTYPE( * )
27* INTEGER IWORK( * ), MVAL( * ), NBVAL( * ), NSVAL( * ),
28* $ NVAL( * ), NXVAL( * )
29* REAL A( * ), COPYA( * ), B( * ), COPYB( * ),
30* $ S( * ), TAU( * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SCHKQP3RK tests SGEQP3RK.
40*> \endverbatim
41*
42* Arguments:
43* ==========
44*
45*> \param[in] DOTYPE
46*> \verbatim
47*> DOTYPE is LOGICAL array, dimension (NTYPES)
48*> The matrix types to be used for testing. Matrices of type j
49*> (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
50*> .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
51*> \endverbatim
52*>
53*> \param[in] NM
54*> \verbatim
55*> NM is INTEGER
56*> The number of values of M contained in the vector MVAL.
57*> \endverbatim
58*>
59*> \param[in] MVAL
60*> \verbatim
61*> MVAL is INTEGER array, dimension (NM)
62*> The values of the matrix row dimension M.
63*> \endverbatim
64*>
65*> \param[in] NN
66*> \verbatim
67*> NN is INTEGER
68*> The number of values of N contained in the vector NVAL.
69*> \endverbatim
70*>
71*> \param[in] NVAL
72*> \verbatim
73*> NVAL is INTEGER array, dimension (NN)
74*> The values of the matrix column dimension N.
75*> \endverbatim
76*>
77*> \param[in] NNS
78*> \verbatim
79*> NNS is INTEGER
80*> The number of values of NRHS contained in the vector NSVAL.
81*> \endverbatim
82*>
83*> \param[in] NSVAL
84*> \verbatim
85*> NSVAL is INTEGER array, dimension (NNS)
86*> The values of the number of right hand sides NRHS.
87*> \endverbatim
88*>
89*> \param[in] NNB
90*> \verbatim
91*> NNB is INTEGER
92*> The number of values of NB and NX contained in the
93*> vectors NBVAL and NXVAL. The blocking parameters are used
94*> in pairs (NB,NX).
95*> \endverbatim
96*>
97*> \param[in] NBVAL
98*> \verbatim
99*> NBVAL is INTEGER array, dimension (NNB)
100*> The values of the blocksize NB.
101*> \endverbatim
102*>
103*> \param[in] NXVAL
104*> \verbatim
105*> NXVAL is INTEGER array, dimension (NNB)
106*> The values of the crossover point NX.
107*> \endverbatim
108*>
109*> \param[in] THRESH
110*> \verbatim
111*> THRESH is REAL
112*> The threshold value for the test ratios. A result is
113*> included in the output file if RESULT >= THRESH. To have
114*> every test ratio printed, use THRESH = 0.
115*> \endverbatim
116*>
117*> \param[out] A
118*> \verbatim
119*> A is REAL array, dimension (MMAX*NMAX)
120*> where MMAX is the maximum value of M in MVAL and NMAX is the
121*> maximum value of N in NVAL.
122*> \endverbatim
123*>
124*> \param[out] COPYA
125*> \verbatim
126*> COPYA is REAL array, dimension (MMAX*NMAX)
127*> \endverbatim
128*>
129*> \param[out] B
130*> \verbatim
131*> B is REAL array, dimension (MMAX*NSMAX)
132*> where MMAX is the maximum value of M in MVAL and NSMAX is the
133*> maximum value of NRHS in NSVAL.
134*> \endverbatim
135*>
136*> \param[out] COPYB
137*> \verbatim
138*> COPYB is REAL array, dimension (MMAX*NSMAX)
139*> \endverbatim
140*>
141*> \param[out] S
142*> \verbatim
143*> S is REAL array, dimension
144*> (min(MMAX,NMAX))
145*> \endverbatim
146*>
147*> \param[out] TAU
148*> \verbatim
149*> TAU is REAL array, dimension (MMAX)
150*> \endverbatim
151*>
152*> \param[out] WORK
153*> \verbatim
154*> WORK is REAL array, dimension
155*> (MMAX*NMAX + 4*NMAX + MMAX)
156*> \endverbatim
157*>
158*> \param[out] IWORK
159*> \verbatim
160*> IWORK is INTEGER array, dimension (2*NMAX)
161*> \endverbatim
162*>
163*> \param[in] NOUT
164*> \verbatim
165*> NOUT is INTEGER
166*> The unit number for output.
167*> \endverbatim
168*
169* Authors:
170* ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup single_lin
178*
179* =====================================================================
180 SUBROUTINE schkqp3rk( DOTYPE, NM, MVAL, NN, NVAL, NNS, NSVAL,
181 $ NNB, NBVAL, NXVAL, THRESH, A, COPYA,
182 $ B, COPYB, S, TAU,
183 $ WORK, IWORK, NOUT )
184 IMPLICIT NONE
185*
186* -- LAPACK test routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 INTEGER NM, NN, NNB, NNS, NOUT
192 REAL THRESH
193* ..
194* .. Array Arguments ..
195 LOGICAL DOTYPE( * )
196 INTEGER IWORK( * ), NBVAL( * ), MVAL( * ), NVAL( * ),
197 $ NSVAL( * ), NXVAL( * )
198 REAL A( * ), COPYA( * ), B( * ), COPYB( * ),
199 $ s( * ), tau( * ), work( * )
200* ..
201*
202* =====================================================================
203*
204* .. Parameters ..
205 INTEGER NTYPES
206 PARAMETER ( NTYPES = 19 )
207 integer ntests
208 parameter( ntests = 5 )
209 REAL ONE, ZERO, BIGNUM
210 parameter( one = 1.0e+0, zero = 0.0e+0,
211 $ bignum = 1.0e+38 )
212* ..
213* .. Local Scalars ..
214 CHARACTER DIST, TYPE
215 CHARACTER*3 PATH
216 INTEGER I, IHIGH, ILOW, IM, IMAT, IN, INC_ZERO,
217 $ inb, ind_offset_gen,
218 $ ind_in, ind_out, ins, info,
219 $ istep, j, j_inc, j_first_nz, jb_zero,
220 $ kfact, kl, kmax, ku, lda, lw, lwork,
221 $ lwork_mqr, m, minmn, minmnb_gen, mode, n,
222 $ nb, nb_zero, nerrs, nfail, nb_gen, nrhs,
223 $ nrun, nx, t
224 REAL ANORM, CNDNUM, EPS, ABSTOL, RELTOL,
225 $ DTEMP, MAXC2NRMK, RELMAXC2NRMK
226* ..
227* .. Local Arrays ..
228 INTEGER ISEED( 4 ), ISEEDY( 4 )
229 REAL RESULT( NTESTS ), RDUMMY( 1 )
230* ..
231* .. External Functions ..
232 REAL SLAMCH, SQPT01, SQRT11, SQRT12, SLANGE
233 EXTERNAL SLAMCH, SQPT01, SQRT11, SQRT12, SLANGE
234* ..
235* .. External Subroutines ..
236 EXTERNAL alaerh, alahd, alasum, saxpy, sgeqp3rk,
237 $ slacpy, slaord, slaset, slatb4, slatms,
238 $ sormqr, sswap, icopy, xlaenv
239* ..
240* .. Intrinsic Functions ..
241 INTRINSIC abs, max, min, mod, real
242* ..
243* .. Scalars in Common ..
244 LOGICAL LERR, OK
245 CHARACTER*32 SRNAMT
246 INTEGER INFOT, IOUNIT
247* ..
248* .. Common blocks ..
249 COMMON / infoc / infot, iounit, ok, lerr
250 COMMON / srnamc / srnamt
251* ..
252* .. Data statements ..
253 DATA iseedy / 1988, 1989, 1990, 1991 /
254* ..
255* .. Executable Statements ..
256*
257* Initialize constants and the random number seed.
258*
259 path( 1: 1 ) = 'Single precision'
260 path( 2: 3 ) = 'QK'
261 nrun = 0
262 nfail = 0
263 nerrs = 0
264 DO i = 1, 4
265 iseed( i ) = iseedy( i )
266 END DO
267 eps = slamch( 'Epsilon' )
268 infot = 0
269*
270 DO im = 1, nm
271*
272* Do for each value of M in MVAL.
273*
274 m = mval( im )
275 lda = max( 1, m )
276*
277 DO in = 1, nn
278*
279* Do for each value of N in NVAL.
280*
281 n = nval( in )
282 minmn = min( m, n )
283 lwork = max( 1, m*max( m, n )+4*minmn+max( m, n ),
284 $ m*n + 2*minmn + 4*n )
285*
286 DO ins = 1, nns
287 nrhs = nsval( ins )
288*
289* Set up parameters with SLATB4 and generate
290* M-by-NRHS B matrix with SLATMS.
291* IMAT = 14:
292* Random matrix, CNDNUM = 2, NORM = ONE,
293* MODE = 3 (geometric distribution of singular values).
294*
295 CALL slatb4( path, 14, m, nrhs, TYPE, kl, ku, anorm,
296 $ mode, cndnum, dist )
297*
298 srnamt = 'SLATMS'
299 CALL slatms( m, nrhs, dist, iseed, TYPE, s, mode,
300 $ cndnum, anorm, kl, ku, 'No packing',
301 $ copyb, lda, work, info )
302
303
304*
305* Check error code from SLATMS.
306*
307 IF( info.NE.0 ) THEN
308 CALL alaerh( path, 'SLATMS', info, 0, ' ', m,
309 $ nrhs, -1, -1, -1, 6, nfail, nerrs,
310 $ nout )
311 cycle
312 END IF
313*
314 DO imat = 1, ntypes
315*
316* Do the tests only if DOTYPE( IMAT ) is true.
317*
318 IF( .NOT.dotype( imat ) )
319 $ cycle
320*
321* The type of distribution used to generate the random
322* eigen-/singular values:
323* ( 'S' for symmetric distribution ) => UNIFORM( -1, 1 )
324*
325* Do for each type of NON-SYMMETRIC matrix: CNDNUM NORM MODE
326* 1. Zero matrix
327* 2. Random, Diagonal, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
328* 3. Random, Upper triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
329* 4. Random, Lower triangular, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
330* 5. Random, First column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
331* 6. Random, Last MINMN column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
332* 7. Random, Last N column is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
333* 8. Random, Middle column in MINMN is zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
334* 9. Random, First half of MINMN columns are zero, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
335* 10. Random, Last columns are zero starting from MINMN/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
336* 11. Random, Half MINMN columns in the middle are zero starting
337* from MINMN/2-(MINMN/2)/2+1, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
338* 12. Random, Odd columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
339* 13. Random, Even columns are ZERO, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
340* 14. Random, CNDNUM = 2 CNDNUM = 2 ONE 3 ( geometric distribution of singular values )
341* 15. Random, CNDNUM = sqrt(0.1/EPS) CNDNUM = BADC1 = sqrt(0.1/EPS) ONE 3 ( geometric distribution of singular values )
342* 16. Random, CNDNUM = 0.1/EPS CNDNUM = BADC2 = 0.1/EPS ONE 3 ( geometric distribution of singular values )
343* 17. Random, CNDNUM = 0.1/EPS, CNDNUM = BADC2 = 0.1/EPS ONE 2 ( one small singular value, S(N)=1/CNDNUM )
344* one small singular value S(N)=1/CNDNUM
345* 18. Random, CNDNUM = 2, scaled near underflow CNDNUM = 2 SMALL = SAFMIN
346* 19. Random, CNDNUM = 2, scaled near overflow CNDNUM = 2 LARGE = 1.0/( 0.25 * ( SAFMIN / EPS ) ) 3 ( geometric distribution of singular values )
347*
348 IF( imat.EQ.1 ) THEN
349*
350* Matrix 1: Zero matrix
351*
352 CALL slaset( 'Full', m, n, zero, zero, copya, lda )
353 DO i = 1, minmn
354 s( i ) = zero
355 END DO
356*
357 ELSE IF( (imat.GE.2 .AND. imat.LE.4 )
358 $ .OR. (imat.GE.14 .AND. imat.LE.19 ) ) THEN
359*
360* Matrices 2-5.
361*
362* Set up parameters with SLATB4 and generate a test
363* matrix with SLATMS.
364*
365 CALL slatb4( path, imat, m, n, TYPE, kl, ku, anorm,
366 $ mode, cndnum, dist )
367*
368 srnamt = 'SLATMS'
369 CALL slatms( m, n, dist, iseed, TYPE, s, mode,
370 $ cndnum, anorm, kl, ku, 'No packing',
371 $ copya, lda, work, info )
372*
373* Check error code from SLATMS.
374*
375 IF( info.NE.0 ) THEN
376 CALL alaerh( path, 'SLATMS', info, 0, ' ', m, n,
377 $ -1, -1, -1, imat, nfail, nerrs,
378 $ nout )
379 cycle
380 END IF
381*
382 CALL slaord( 'Decreasing', minmn, s, 1 )
383*
384 ELSE IF( minmn.GE.2
385 $ .AND. imat.GE.5 .AND. imat.LE.13 ) THEN
386*
387* Rectangular matrices 5-13 that contain zero columns,
388* only for matrices MINMN >=2.
389*
390* JB_ZERO is the column index of ZERO block.
391* NB_ZERO is the column block size of ZERO block.
392* NB_GEN is the column blcok size of the
393* generated block.
394* J_INC in the non_zero column index increment
395* for matrix 12 and 13.
396* J_FIRS_NZ is the index of the first non-zero
397* column.
398*
399 IF( imat.EQ.5 ) THEN
400*
401* First column is zero.
402*
403 jb_zero = 1
404 nb_zero = 1
405 nb_gen = n - nb_zero
406*
407 ELSE IF( imat.EQ.6 ) THEN
408*
409* Last column MINMN is zero.
410*
411 jb_zero = minmn
412 nb_zero = 1
413 nb_gen = n - nb_zero
414*
415 ELSE IF( imat.EQ.7 ) THEN
416*
417* Last column N is zero.
418*
419 jb_zero = n
420 nb_zero = 1
421 nb_gen = n - nb_zero
422*
423 ELSE IF( imat.EQ.8 ) THEN
424*
425* Middle column in MINMN is zero.
426*
427 jb_zero = minmn / 2 + 1
428 nb_zero = 1
429 nb_gen = n - nb_zero
430*
431 ELSE IF( imat.EQ.9 ) THEN
432*
433* First half of MINMN columns is zero.
434*
435 jb_zero = 1
436 nb_zero = minmn / 2
437 nb_gen = n - nb_zero
438*
439 ELSE IF( imat.EQ.10 ) THEN
440*
441* Last columns are zero columns,
442* starting from (MINMN / 2 + 1) column.
443*
444 jb_zero = minmn / 2 + 1
445 nb_zero = n - jb_zero + 1
446 nb_gen = n - nb_zero
447*
448 ELSE IF( imat.EQ.11 ) THEN
449*
450* Half of the columns in the middle of MINMN
451* columns is zero, starting from
452* MINMN/2 - (MINMN/2)/2 + 1 column.
453*
454 jb_zero = minmn / 2 - (minmn / 2) / 2 + 1
455 nb_zero = minmn / 2
456 nb_gen = n - nb_zero
457*
458 ELSE IF( imat.EQ.12 ) THEN
459*
460* Odd-numbered columns are zero,
461*
462 nb_gen = n / 2
463 nb_zero = n - nb_gen
464 j_inc = 2
465 j_first_nz = 2
466*
467 ELSE IF( imat.EQ.13 ) THEN
468*
469* Even-numbered columns are zero.
470*
471 nb_zero = n / 2
472 nb_gen = n - nb_zero
473 j_inc = 2
474 j_first_nz = 1
475*
476 END IF
477*
478*
479* 1) Set the first NB_ZERO columns in COPYA(1:M,1:N)
480* to zero.
481*
482 CALL slaset( 'Full', m, nb_zero, zero, zero,
483 $ copya, lda )
484*
485* 2) Generate an M-by-(N-NB_ZERO) matrix with the
486* chosen singular value distribution
487* in COPYA(1:M,NB_ZERO+1:N).
488*
489 CALL slatb4( path, imat, m, nb_gen, TYPE, kl, ku,
490 $ anorm, mode, cndnum, dist )
491*
492 srnamt = 'SLATMS'
493*
494 ind_offset_gen = nb_zero * lda
495*
496 CALL slatms( m, nb_gen, dist, iseed, TYPE, s, mode,
497 $ cndnum, anorm, kl, ku, 'No packing',
498 $ copya( ind_offset_gen + 1 ), lda,
499 $ work, info )
500*
501* Check error code from SLATMS.
502*
503 IF( info.NE.0 ) THEN
504 CALL alaerh( path, 'SLATMS', info, 0, ' ', m,
505 $ nb_gen, -1, -1, -1, imat, nfail,
506 $ nerrs, nout )
507 cycle
508 END IF
509*
510* 3) Swap the gererated colums from the right side
511* NB_GEN-size block in COPYA into correct column
512* positions.
513*
514 IF( imat.EQ.6
515 $ .OR. imat.EQ.7
516 $ .OR. imat.EQ.8
517 $ .OR. imat.EQ.10
518 $ .OR. imat.EQ.11 ) THEN
519*
520* Move by swapping the generated columns
521* from the right NB_GEN-size block from
522* (NB_ZERO+1:NB_ZERO+JB_ZERO)
523* into columns (1:JB_ZERO-1).
524*
525 DO j = 1, jb_zero-1, 1
526 CALL sswap( m,
527 $ copya( ( nb_zero+j-1)*lda+1), 1,
528 $ copya( (j-1)*lda + 1 ), 1 )
529 END DO
530*
531 ELSE IF( imat.EQ.12 .OR. imat.EQ.13 ) THEN
532*
533* ( IMAT = 12, Odd-numbered ZERO columns. )
534* Swap the generated columns from the right
535* NB_GEN-size block into the even zero colums in the
536* left NB_ZERO-size block.
537*
538* ( IMAT = 13, Even-numbered ZERO columns. )
539* Swap the generated columns from the right
540* NB_GEN-size block into the odd zero colums in the
541* left NB_ZERO-size block.
542*
543 DO j = 1, nb_gen, 1
544 ind_out = ( nb_zero+j-1 )*lda + 1
545 ind_in = ( j_inc*(j-1)+(j_first_nz-1) )*lda
546 $ + 1
547 CALL sswap( m,
548 $ copya( ind_out ), 1,
549 $ copya( ind_in), 1 )
550 END DO
551*
552 END IF
553*
554* 5) Order the singular values generated by
555* DLAMTS in decreasing order and add trailing zeros
556* that correspond to zero columns.
557* The total number of singular values is MINMN.
558*
559 minmnb_gen = min( m, nb_gen )
560*
561 DO i = minmnb_gen+1, minmn
562 s( i ) = zero
563 END DO
564*
565 ELSE
566*
567* IF(MINMN.LT.2) skip this size for this matrix type.
568*
569 cycle
570 END IF
571*
572* Initialize a copy array for a pivot array for SGEQP3RK.
573*
574 DO i = 1, n
575 iwork( i ) = 0
576 END DO
577*
578 DO inb = 1, nnb
579*
580* Do for each pair of values (NB,NX) in NBVAL and NXVAL.
581*
582 nb = nbval( inb )
583 CALL xlaenv( 1, nb )
584 nx = nxval( inb )
585 CALL xlaenv( 3, nx )
586*
587* We do MIN(M,N)+1 because we need a test for KMAX > N,
588* when KMAX is larger than MIN(M,N), KMAX should be
589* KMAX = MIN(M,N)
590*
591 DO kmax = 0, min(m,n)+1
592*
593* Get a working copy of COPYA into A( 1:M,1:N ).
594* Get a working copy of COPYB into A( 1:M, (N+1):NRHS ).
595* Get a working copy of COPYB into into B( 1:M, 1:NRHS ).
596* Get a working copy of IWORK(1:N) awith zeroes into
597* which is going to be used as pivot array IWORK( N+1:2N ).
598* NOTE: IWORK(2N+1:3N) is going to be used as a WORK array
599* for the routine.
600*
601 CALL slacpy( 'All', m, n, copya, lda, a, lda )
602 CALL slacpy( 'All', m, nrhs, copyb, lda,
603 $ a( lda*n + 1 ), lda )
604 CALL slacpy( 'All', m, nrhs, copyb, lda,
605 $ b, lda )
606 CALL icopy( n, iwork( 1 ), 1, iwork( n+1 ), 1 )
607*
608 abstol = -1.0
609 reltol = -1.0
610*
611* Compute the QR factorization with pivoting of A
612*
613 lw = max( 1, max( 2*n + nb*( n+nrhs+1 ),
614 $ 3*n + nrhs - 1 ) )
615*
616* Compute SGEQP3RK factorization of A.
617*
618 srnamt = 'SGEQP3RK'
619 CALL sgeqp3rk( m, n, nrhs, kmax, abstol, reltol,
620 $ a, lda, kfact, maxc2nrmk,
621 $ relmaxc2nrmk, iwork( n+1 ), tau,
622 $ work, lw, iwork( 2*n+1 ), info )
623*
624* Check error code from SGEQP3RK.
625*
626 IF( info.LT.0 )
627 $ CALL alaerh( path, 'SGEQP3RK', info, 0, ' ',
628 $ m, n, nx, -1, nb, imat,
629 $ nfail, nerrs, nout )
630*
631* Compute test 1:
632*
633* This test in only for the full rank factorization of
634* the matrix A.
635*
636* Array S(1:min(M,N)) contains svd(A) the sigular values
637* of the original matrix A in decreasing absolute value
638* order. The test computes svd(R), the vector sigular
639* values of the upper trapezoid of A(1:M,1:N) that
640* contains the factor R, in decreasing order. The test
641* returns the ratio:
642*
643* 2-norm(svd(R) - svd(A)) / ( max(M,N) * 2-norm(svd(A)) * EPS )
644*
645 IF( kfact.EQ.minmn ) THEN
646*
647 result( 1 ) = sqrt12( m, n, a, lda, s, work,
648 $ lwork )
649*
650 DO t = 1, 1
651 IF( result( t ).GE.thresh ) THEN
652 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
653 $ CALL alahd( nout, path )
654 WRITE( nout, fmt = 9999 ) 'SGEQP3RK', m, n,
655 $ nrhs, kmax, abstol, reltol, nb, nx,
656 $ imat, t, result( t )
657 nfail = nfail + 1
658 END IF
659 END DO
660 nrun = nrun + 1
661*
662* End test 1
663*
664 END IF
665*
666* Compute test 2:
667*
668* The test returns the ratio:
669*
670* 1-norm( A*P - Q*R ) / ( max(M,N) * 1-norm(A) * EPS )
671*
672 result( 2 ) = sqpt01( m, n, kfact, copya, a, lda, tau,
673 $ iwork( n+1 ), work, lwork )
674*
675* Compute test 3:
676*
677* The test returns the ratio:
678*
679* 1-norm( Q**T * Q - I ) / ( M * EPS )
680*
681 result( 3 ) = sqrt11( m, kfact, a, lda, tau, work,
682 $ lwork )
683*
684* Print information about the tests that did not pass
685* the threshold.
686*
687 DO t = 2, 3
688 IF( result( t ).GE.thresh ) THEN
689 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
690 $ CALL alahd( nout, path )
691 WRITE( nout, fmt = 9999 ) 'SGEQP3RK', m, n,
692 $ nrhs, kmax, abstol, reltol,
693 $ nb, nx, imat, t, result( t )
694 nfail = nfail + 1
695 END IF
696 END DO
697 nrun = nrun + 2
698*
699* Compute test 4:
700*
701* This test is only for the factorizations with the
702* rank greater than 2.
703* The elements on the diagonal of R should be non-
704* increasing.
705*
706* The test returns the ratio:
707*
708* Returns 1.0D+100 if abs(R(K+1,K+1)) > abs(R(K,K)),
709* K=1:KFACT-1
710*
711 IF( min(kfact, minmn).GE.2 ) THEN
712*
713 DO j = 1, kfact-1, 1
714
715 dtemp = (( abs( a( (j-1)*m+j ) ) -
716 $ abs( a( (j)*m+j+1 ) ) ) /
717 $ abs( a(1) ) )
718*
719 IF( dtemp.LT.zero ) THEN
720 result( 4 ) = bignum
721 END IF
722*
723 END DO
724*
725* Print information about the tests that did not
726* pass the threshold.
727*
728 DO t = 4, 4
729 IF( result( t ).GE.thresh ) THEN
730 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
731 $ CALL alahd( nout, path )
732 WRITE( nout, fmt = 9999 ) 'SGEQP3RK',
733 $ m, n, nrhs, kmax, abstol, reltol,
734 $ nb, nx, imat, t,
735 $ result( t )
736 nfail = nfail + 1
737 END IF
738 END DO
739 nrun = nrun + 1
740*
741* End test 4.
742*
743 END IF
744*
745* Compute test 5:
746*
747* This test in only for matrix A with min(M,N) > 0.
748*
749* The test returns the ratio:
750*
751* 1-norm(Q**T * B - Q**T * B ) /
752* ( M * EPS )
753*
754* (1) Compute B:=Q**T * B in the matrix B.
755*
756 IF( minmn.GT.0 ) THEN
757*
758 lwork_mqr = max(1, nrhs)
759 CALL sormqr( 'Left', 'Transpose',
760 $ m, nrhs, kfact, a, lda, tau, b, lda,
761 $ work, lwork_mqr, info )
762*
763 DO i = 1, nrhs
764*
765* Compare N+J-th column of A and J-column of B.
766*
767 CALL saxpy( m, -one, a( ( n+i-1 )*lda+1 ), 1,
768 $ b( ( i-1 )*lda+1 ), 1 )
769 END DO
770*
771 result( 5 ) =
772 $ abs(
773 $ slange( 'One-norm', m, nrhs, b, lda, rdummy ) /
774 $ ( real( m )*slamch( 'Epsilon' ) )
775 $ )
776*
777* Print information about the tests that did not pass
778* the threshold.
779*
780 DO t = 5, 5
781 IF( result( t ).GE.thresh ) THEN
782 IF( nfail.EQ.0 .AND. nerrs.EQ.0 )
783 $ CALL alahd( nout, path )
784 WRITE( nout, fmt = 9999 ) 'SGEQP3RK', m, n,
785 $ nrhs, kmax, abstol, reltol,
786 $ nb, nx, imat, t, result( t )
787 nfail = nfail + 1
788 END IF
789 END DO
790 nrun = nrun + 1
791*
792* End compute test 5.
793*
794 END IF
795*
796* END DO KMAX = 1, MIN(M,N)+1
797*
798 END DO
799*
800* END DO for INB = 1, NNB
801*
802 END DO
803*
804* END DO for IMAT = 1, NTYPES
805*
806 END DO
807*
808* END DO for INS = 1, NNS
809*
810 END DO
811*
812* END DO for IN = 1, NN
813*
814 END DO
815*
816* END DO for IM = 1, NM
817*
818 END DO
819*
820* Print a summary of the results.
821*
822 CALL alasum( path, nout, nfail, nrun, nerrs )
823*
824 9999 FORMAT( 1x, a, ' M =', i5, ', N =', i5, ', NRHS =', i5,
825 $ ', KMAX =', i5, ', ABSTOL =', g12.5,
826 $ ', RELTOL =', g12.5, ', NB =', i4, ', NX =', i4,
827 $ ', type ', i2, ', test ', i2, ', ratio =', g12.5 )
828*
829* End of SCHKQP3RK
830*
831 END
alasum
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.f:73
xlaenv
subroutine xlaenv(ispec, nvalue)
XLAENV
Definition xlaenv.f:81
alaerh
subroutine alaerh(path, subnam, info, infoe, opts, m, n, kl, ku, n5, imat, nfail, nerrs, nout)
ALAERH
Definition alaerh.f:147
alahd
subroutine alahd(iounit, path)
ALAHD
Definition alahd.f:107
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:110
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:168
icopy
subroutine icopy(n, sx, incx, sy, incy)
ICOPY
Definition icopy.f:75
schkqp3rk
subroutine schkqp3rk(dotype, nm, mval, nn, nval, nns, nsval, nnb, nbval, nxval, thresh, a, copya, b, copyb, s, tau, work, iwork, nout)
SCHKQP3RK
Definition schkqp3rk.f:184
sgeqp3rk
subroutine sgeqp3rk(m, n, nrhs, kmax, abstol, reltol, a, lda, k, maxc2nrmk, relmaxc2nrmk, jpiv, tau, work, lwork, iwork, info)
SGEQP3RK computes a truncated Householder QR factorization with column pivoting of a real m-by-n matr...
Definition sgeqp3rk.f:585
slaord
subroutine slaord(job, n, x, incx)
SLAORD
Definition slaord.f:73
slatb4
subroutine slatb4(path, imat, m, n, type, kl, ku, anorm, mode, cndnum, dist)
SLATB4
Definition slatb4.f:120
slatms
subroutine slatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
SLATMS
Definition slatms.f:321