1*> \brief \b DDRVRF4
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 DDRVRF4( NOUT, NN, NVAL, THRESH, C1, C2, LDC, CRF, A,
12*      +                    LDA, D_WORK_DLANGE )
13*
14*       .. Scalar Arguments ..
15*       INTEGER            LDA, LDC, NN, NOUT
16*       DOUBLE PRECISION   THRESH
17*       ..
18*       .. Array Arguments ..
19*       INTEGER            NVAL( NN )
20*       DOUBLE PRECISION   A( LDA, * ), C1( LDC, * ), C2( LDC, *),
21*      +                   CRF( * ), D_WORK_DLANGE( * )
22*       ..
23*
24*
25*> \par Purpose:
26*  =============
27*>
28*> \verbatim
29*>
30*> DDRVRF4 tests the LAPACK RFP routines:
31*>     DSFRK
32*> \endverbatim
33*
34*  Arguments:
35*  ==========
36*
37*> \param[in] NOUT
38*> \verbatim
39*>          NOUT is INTEGER
40*>                The unit number for output.
41*> \endverbatim
42*>
43*> \param[in] NN
44*> \verbatim
45*>          NN is INTEGER
46*>                The number of values of N contained in the vector NVAL.
47*> \endverbatim
48*>
49*> \param[in] NVAL
50*> \verbatim
51*>          NVAL is INTEGER array, dimension (NN)
52*>                The values of the matrix dimension N.
53*> \endverbatim
54*>
55*> \param[in] THRESH
56*> \verbatim
57*>          THRESH is DOUBLE PRECISION
58*>                The threshold value for the test ratios.  A result is
59*>                included in the output file if RESULT >= THRESH.  To
60*>                have every test ratio printed, use THRESH = 0.
61*> \endverbatim
62*>
63*> \param[out] C1
64*> \verbatim
65*>          C1 is DOUBLE PRECISION array,
66*>                dimension (LDC,NMAX)
67*> \endverbatim
68*>
69*> \param[out] C2
70*> \verbatim
71*>          C2 is DOUBLE PRECISION array,
72*>                dimension (LDC,NMAX)
73*> \endverbatim
74*>
75*> \param[in] LDC
76*> \verbatim
77*>          LDC is INTEGER
78*>                The leading dimension of the array A.
79*>                LDA >= max(1,NMAX).
80*> \endverbatim
81*>
82*> \param[out] CRF
83*> \verbatim
84*>          CRF is DOUBLE PRECISION array,
85*>                dimension ((NMAX*(NMAX+1))/2).
86*> \endverbatim
87*>
88*> \param[out] A
89*> \verbatim
90*>          A is DOUBLE PRECISION array,
91*>                dimension (LDA,NMAX)
92*> \endverbatim
93*>
94*> \param[in] LDA
95*> \verbatim
96*>          LDA is INTEGER
97*>                The leading dimension of the array A.  LDA >= max(1,NMAX).
98*> \endverbatim
99*>
100*> \param[out] D_WORK_DLANGE
101*> \verbatim
102*>          D_WORK_DLANGE is DOUBLE PRECISION array, dimension (NMAX)
103*> \endverbatim
104*
105*  Authors:
106*  ========
107*
108*> \author Univ. of Tennessee
109*> \author Univ. of California Berkeley
110*> \author Univ. of Colorado Denver
111*> \author NAG Ltd.
112*
113*> \ingroup double_lin
114*
115*  =====================================================================
116      SUBROUTINE DDRVRF4( NOUT, NN, NVAL, THRESH, C1, C2, LDC, CRF, A,
117     +                    LDA, D_WORK_DLANGE )
118*
119*  -- LAPACK test routine --
120*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
121*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
122*
123*     .. Scalar Arguments ..
124      INTEGER            LDA, LDC, NN, NOUT
125      DOUBLE PRECISION   THRESH
126*     ..
127*     .. Array Arguments ..
128      INTEGER            NVAL( NN )
129      DOUBLE PRECISION   A( LDA, * ), C1( LDC, * ), C2( LDC, *),
130     +                   CRF( * ), D_WORK_DLANGE( * )
131*     ..
132*
133*  =====================================================================
134*     ..
135*     .. Parameters ..
136      DOUBLE PRECISION   ZERO, ONE
137      PARAMETER          ( ZERO = 0.0D+0, ONE  = 1.0D+0 )
138      INTEGER            NTESTS
139      PARAMETER          ( NTESTS = 1 )
140*     ..
141*     .. Local Scalars ..
142      CHARACTER          UPLO, CFORM, TRANS
143      INTEGER            I, IFORM, IIK, IIN, INFO, IUPLO, J, K, N,
144     +                   NFAIL, NRUN, IALPHA, ITRANS
145      DOUBLE PRECISION   ALPHA, BETA, EPS, NORMA, NORMC
146*     ..
147*     .. Local Arrays ..
148      CHARACTER          UPLOS( 2 ), FORMS( 2 ), TRANSS( 2 )
149      INTEGER            ISEED( 4 ), ISEEDY( 4 )
150      DOUBLE PRECISION   RESULT( NTESTS )
151*     ..
152*     .. External Functions ..
153      DOUBLE PRECISION   DLAMCH, DLARND, DLANGE
154      EXTERNAL           DLAMCH, DLARND, DLANGE
155*     ..
156*     .. External Subroutines ..
157      EXTERNAL           DSYRK, DSFRK, DTFTTR, DTRTTF
158*     ..
159*     .. Intrinsic Functions ..
160      INTRINSIC          ABS, MAX
161*     ..
162*     .. Scalars in Common ..
163      CHARACTER*32       SRNAMT
164*     ..
165*     .. Common blocks ..
166      COMMON             / SRNAMC / SRNAMT
167*     ..
168*     .. Data statements ..
169      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
170      DATA               UPLOS  / 'U', 'L' /
171      DATA               FORMS  / 'N', 'T' /
172      DATA               TRANSS / 'N', 'T' /
173*     ..
174*     .. Executable Statements ..
175*
176*     Initialize constants and the random number seed.
177*
178      NRUN = 0
179      NFAIL = 0
180      INFO = 0
181      DO 10 I = 1, 4
182         ISEED( I ) = ISEEDY( I )
183   10 CONTINUE
184      EPS = DLAMCH( 'Precision' )
185*
186      DO 150 IIN = 1, NN
187*
188         N = NVAL( IIN )
189*
190         DO 140 IIK = 1, NN
191*
192            K = NVAL( IIN )
193*
194            DO 130 IFORM = 1, 2
195*
196               CFORM = FORMS( IFORM )
197*
198               DO 120 IUPLO = 1, 2
199*
200                  UPLO = UPLOS( IUPLO )
201*
202                  DO 110 ITRANS = 1, 2
203*
204                     TRANS = TRANSS( ITRANS )
205*
206                     DO 100 IALPHA = 1, 4
207*
208                        IF ( IALPHA.EQ. 1) THEN
209                           ALPHA = ZERO
210                           BETA = ZERO
211                        ELSE IF ( IALPHA.EQ. 2) THEN
212                           ALPHA = ONE
213                           BETA = ZERO
214                        ELSE IF ( IALPHA.EQ. 3) THEN
215                           ALPHA = ZERO
216                           BETA = ONE
217                        ELSE
218                           ALPHA = DLARND( 2, ISEED )
219                           BETA = DLARND( 2, ISEED )
220                        END IF
221*
222*                       All the parameters are set:
223*                          CFORM, UPLO, TRANS, M, N,
224*                          ALPHA, and BETA
225*                       READY TO TEST!
226*
227                        NRUN = NRUN + 1
228*
229                        IF ( ITRANS.EQ.1 ) THEN
230*
231*                          In this case we are NOTRANS, so A is N-by-K
232*
233                           DO J = 1, K
234                              DO I = 1, N
235                                 A( I, J) = DLARND( 2, ISEED )
236                              END DO
237                           END DO
238*
239                           NORMA = DLANGE( 'I', N, K, A, LDA,
240     +                                      D_WORK_DLANGE )
241*
242
243                        ELSE
244*
245*                          In this case we are TRANS, so A is K-by-N
246*
247                           DO J = 1,N
248                              DO I = 1, K
249                                 A( I, J) = DLARND( 2, ISEED )
250                              END DO
251                           END DO
252*
253                           NORMA = DLANGE( 'I', K, N, A, LDA,
254     +                                      D_WORK_DLANGE )
255*
256                        END IF
257*
258*                       Generate C1 our N--by--N symmetric matrix.
259*                       Make sure C2 has the same upper/lower part,
260*                       (the one that we do not touch), so
261*                       copy the initial C1 in C2 in it.
262*
263                        DO J = 1, N
264                           DO I = 1, N
265                              C1( I, J) = DLARND( 2, ISEED )
266                              C2(I,J) = C1(I,J)
267                           END DO
268                        END DO
269*
270*                       (See comment later on for why we use DLANGE and
271*                       not DLANSY for C1.)
272*
273                        NORMC = DLANGE( 'I', N, N, C1, LDC,
274     +                                      D_WORK_DLANGE )
275*
276                        SRNAMT = 'DTRTTF'
277                        CALL DTRTTF( CFORM, UPLO, N, C1, LDC, CRF,
278     +                               INFO )
279*
280*                       call dsyrk the BLAS routine -> gives C1
281*
282                        SRNAMT = 'DSYRK '
283                        CALL DSYRK( UPLO, TRANS, N, K, ALPHA, A, LDA,
284     +                              BETA, C1, LDC )
285*
286*                       call dsfrk the RFP routine -> gives CRF
287*
288                        SRNAMT = 'DSFRK '
289                        CALL DSFRK( CFORM, UPLO, TRANS, N, K, ALPHA, A,
290     +                              LDA, BETA, CRF )
291*
292*                       convert CRF in full format -> gives C2
293*
294                        SRNAMT = 'DTFTTR'
295                        CALL DTFTTR( CFORM, UPLO, N, CRF, C2, LDC,
296     +                               INFO )
297*
298*                       compare C1 and C2
299*
300                        DO J = 1, N
301                           DO I = 1, N
302                              C1(I,J) = C1(I,J)-C2(I,J)
303                           END DO
304                        END DO
305*
306*                       Yes, C1 is symmetric so we could call DLANSY,
307*                       but we want to check the upper part that is
308*                       supposed to be unchanged and the diagonal that
309*                       is supposed to be real -> DLANGE
310*
311                        RESULT(1) = DLANGE( 'I', N, N, C1, LDC,
312     +                                      D_WORK_DLANGE )
313                        RESULT(1) = RESULT(1)
314     +                              / MAX( ABS( ALPHA ) * NORMA
315     +                                   + ABS( BETA ) , ONE )
316     +                              / MAX( N , 1 ) / EPS
317*
318                        IF( RESULT(1).GE.THRESH ) THEN
319                           IF( NFAIL.EQ.0 ) THEN
320                              WRITE( NOUT, * )
321                              WRITE( NOUT, FMT = 9999 )
322                           END IF
323                           WRITE( NOUT, FMT = 9997 ) 'DSFRK',
324     +                        CFORM, UPLO, TRANS, N, K, RESULT(1)
325                           NFAIL = NFAIL + 1
326                        END IF
327*
328  100                CONTINUE
329  110             CONTINUE
330  120          CONTINUE
331  130       CONTINUE
332  140    CONTINUE
333  150 CONTINUE
334*
335*     Print a summary of the results.
336*
337      IF ( NFAIL.EQ.0 ) THEN
338         WRITE( NOUT, FMT = 9996 ) 'DSFRK', NRUN
339      ELSE
340         WRITE( NOUT, FMT = 9995 ) 'DSFRK', NFAIL, NRUN
341      END IF
342*
343 9999 FORMAT( 1X, ' *** Error(s) or Failure(s) while testing DSFRK
344     +         ***')
345 9997 FORMAT( 1X, '     Failure in ',A5,', CFORM=''',A1,''',',
346     + ' UPLO=''',A1,''',',' TRANS=''',A1,''',', ' N=',I3,', K =', I3,
347     + ', test=',G12.5)
348 9996 FORMAT( 1X, 'All tests for ',A5,' auxiliary routine passed the ',
349     +        'threshold ( ',I5,' tests run)')
350 9995 FORMAT( 1X, A6, ' auxiliary routine: ',I5,' out of ',I5,
351     +        ' tests failed to pass the threshold')
352*
353      RETURN
354*
355*     End of DDRVRF4
356*
357      END
358