1 SUBROUTINE DDRVRFP( NOUT, NN, NVAL, NNS, NSVAL, NNT, NTVAL,
2 + THRESH, A, ASAV, AFAC, AINV, B,
3 + BSAV, XACT, X, ARF, ARFINV,
4 + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02,
5 + D_TEMP_DPOT03, D_WORK_DLANSY,
6 + D_WORK_DPOT02, D_WORK_DPOT03 )
7 *
8 * -- LAPACK test routine (version 3.2.0) --
9 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
10 * November 2008
11 *
12 * .. Scalar Arguments ..
13 INTEGER NN, NNS, NNT, NOUT
14 DOUBLE PRECISION THRESH
15 * ..
16 * .. Array Arguments ..
17 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
18 DOUBLE PRECISION A( * )
19 DOUBLE PRECISION AINV( * )
20 DOUBLE PRECISION ASAV( * )
21 DOUBLE PRECISION B( * )
22 DOUBLE PRECISION BSAV( * )
23 DOUBLE PRECISION AFAC( * )
24 DOUBLE PRECISION ARF( * )
25 DOUBLE PRECISION ARFINV( * )
26 DOUBLE PRECISION XACT( * )
27 DOUBLE PRECISION X( * )
28 DOUBLE PRECISION D_WORK_DLATMS( * )
29 DOUBLE PRECISION D_WORK_DPOT01( * )
30 DOUBLE PRECISION D_TEMP_DPOT02( * )
31 DOUBLE PRECISION D_TEMP_DPOT03( * )
32 DOUBLE PRECISION D_WORK_DLANSY( * )
33 DOUBLE PRECISION D_WORK_DPOT02( * )
34 DOUBLE PRECISION D_WORK_DPOT03( * )
35 * ..
36 *
37 * Purpose
38 * =======
39 *
40 * DDRVRFP tests the LAPACK RFP routines:
41 * DPFTRF, DPFTRS, and DPFTRI.
42 *
43 * This testing routine follow the same tests as DDRVPO (test for the full
44 * format Symmetric Positive Definite solver).
45 *
46 * The tests are performed in Full Format, convertion back and forth from
47 * full format to RFP format are performed using the routines DTRTTF and
48 * DTFTTR.
49 *
50 * First, a specific matrix A of size N is created. There is nine types of
51 * different matrixes possible.
52 * 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS)
53 * 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS
54 * *3. First row and column zero 8. Scaled near underflow
55 * *4. Last row and column zero 9. Scaled near overflow
56 * *5. Middle row and column zero
57 * (* - tests error exits from DPFTRF, no test ratios are computed)
58 * A solution XACT of size N-by-NRHS is created and the associated right
59 * hand side B as well. Then DPFTRF is called to compute L (or U), the
60 * Cholesky factor of A. Then L (or U) is used to solve the linear system
61 * of equations AX = B. This gives X. Then L (or U) is used to compute the
62 * inverse of A, AINV. The following four tests are then performed:
63 * (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
64 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
65 * (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
66 * (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
67 * (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
68 * where EPS is the machine precision, RCOND the condition number of A, and
69 * norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
70 * Errors occur when INFO parameter is not as expected. Failures occur when
71 * a test ratios is greater than THRES.
72 *
73 * Arguments
74 * =========
75 *
76 * NOUT (input) INTEGER
77 * The unit number for output.
78 *
79 * NN (input) INTEGER
80 * The number of values of N contained in the vector NVAL.
81 *
82 * NVAL (input) INTEGER array, dimension (NN)
83 * The values of the matrix dimension N.
84 *
85 * NNS (input) INTEGER
86 * The number of values of NRHS contained in the vector NSVAL.
87 *
88 * NSVAL (input) INTEGER array, dimension (NNS)
89 * The values of the number of right-hand sides NRHS.
90 *
91 * NNT (input) INTEGER
92 * The number of values of MATRIX TYPE contained in the vector NTVAL.
93 *
94 * NTVAL (input) INTEGER array, dimension (NNT)
95 * The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
96 *
97 * THRESH (input) DOUBLE PRECISION
98 * The threshold value for the test ratios. A result is
99 * included in the output file if RESULT >= THRESH. To have
100 * every test ratio printed, use THRESH = 0.
101 *
102 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
103 *
104 * ASAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
105 *
106 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
107 *
108 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
109 *
110 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
111 *
112 * BSAV (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
113 *
114 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
115 *
116 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
117 *
118 * ARF (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
119 *
120 * ARFINV (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
121 *
122 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( 3*NMAX )
123 *
124 * D_WORK_DPOT01 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
125 *
126 * D_TEMP_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX*MAXRHS )
127 *
128 * D_TEMP_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX*NMAX )
129 *
130 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( NMAX )
131 *
132 * D_WORK_DLANSY (workspace) DOUBLE PRECISION array, dimension ( NMAX )
133 *
134 * D_WORK_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
135 *
136 * D_WORK_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141 DOUBLE PRECISION ONE, ZERO
142 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
143 INTEGER NTESTS
144 PARAMETER ( NTESTS = 4 )
145 * ..
146 * .. Local Scalars ..
147 LOGICAL ZEROT
148 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL,
149 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN,
150 + IIT, IIS
151 CHARACTER DIST, CTYPE, UPLO, CFORM
152 INTEGER KL, KU, MODE
153 DOUBLE PRECISION ANORM, AINVNM, CNDNUM, RCONDC
154 * ..
155 * .. Local Arrays ..
156 CHARACTER UPLOS( 2 ), FORMS( 2 )
157 INTEGER ISEED( 4 ), ISEEDY( 4 )
158 DOUBLE PRECISION RESULT( NTESTS )
159 * ..
160 * .. External Functions ..
161 DOUBLE PRECISION DLANSY
162 EXTERNAL DLANSY
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL ALADHD, ALAERH, ALASVM, DGET04, DTFTTR, DLACPY,
166 + DLARHS, DLATB4, DLATMS, DPFTRI, DPFTRF, DPFTRS,
167 + DPOT01, DPOT02, DPOT03, DPOTRI, DPOTRF, DTRTTF
168 * ..
169 * .. Scalars in Common ..
170 CHARACTER*32 SRNAMT
171 * ..
172 * .. Common blocks ..
173 COMMON / SRNAMC / SRNAMT
174 * ..
175 * .. Data statements ..
176 DATA ISEEDY / 1988, 1989, 1990, 1991 /
177 DATA UPLOS / 'U', 'L' /
178 DATA FORMS / 'N', 'T' /
179 * ..
180 * .. Executable Statements ..
181 *
182 * Initialize constants and the random number seed.
183 *
184 NRUN = 0
185 NFAIL = 0
186 NERRS = 0
187 DO 10 I = 1, 4
188 ISEED( I ) = ISEEDY( I )
189 10 CONTINUE
190 *
191 DO 130 IIN = 1, NN
192 *
193 N = NVAL( IIN )
194 LDA = MAX( N, 1 )
195 LDB = MAX( N, 1 )
196 *
197 DO 980 IIS = 1, NNS
198 *
199 NRHS = NSVAL( IIS )
200 *
201 DO 120 IIT = 1, NNT
202 *
203 IMAT = NTVAL( IIT )
204 *
205 * If N.EQ.0, only consider the first type
206 *
207 IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120
208 *
209 * Skip types 3, 4, or 5 if the matrix size is too small.
210 *
211 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120
212 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120
213 *
214 * Do first for UPLO = 'U', then for UPLO = 'L'
215 *
216 DO 110 IUPLO = 1, 2
217 UPLO = UPLOS( IUPLO )
218 *
219 * Do first for CFORM = 'N', then for CFORM = 'C'
220 *
221 DO 100 IFORM = 1, 2
222 CFORM = FORMS( IFORM )
223 *
224 * Set up parameters with DLATB4 and generate a test
225 * matrix with DLATMS.
226 *
227 CALL DLATB4( 'DPO', IMAT, N, N, CTYPE, KL, KU,
228 + ANORM, MODE, CNDNUM, DIST )
229 *
230 SRNAMT = 'DLATMS'
231 CALL DLATMS( N, N, DIST, ISEED, CTYPE,
232 + D_WORK_DLATMS,
233 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A,
234 + LDA, D_WORK_DLATMS, INFO )
235 *
236 * Check error code from DLATMS.
237 *
238 IF( INFO.NE.0 ) THEN
239 CALL ALAERH( 'DPF', 'DLATMS', INFO, 0, UPLO, N,
240 + N, -1, -1, -1, IIT, NFAIL, NERRS,
241 + NOUT )
242 GO TO 100
243 END IF
244 *
245 * For types 3-5, zero one row and column of the matrix to
246 * test that INFO is returned correctly.
247 *
248 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
249 IF( ZEROT ) THEN
250 IF( IIT.EQ.3 ) THEN
251 IZERO = 1
252 ELSE IF( IIT.EQ.4 ) THEN
253 IZERO = N
254 ELSE
255 IZERO = N / 2 + 1
256 END IF
257 IOFF = ( IZERO-1 )*LDA
258 *
259 * Set row and column IZERO of A to 0.
260 *
261 IF( IUPLO.EQ.1 ) THEN
262 DO 20 I = 1, IZERO - 1
263 A( IOFF+I ) = ZERO
264 20 CONTINUE
265 IOFF = IOFF + IZERO
266 DO 30 I = IZERO, N
267 A( IOFF ) = ZERO
268 IOFF = IOFF + LDA
269 30 CONTINUE
270 ELSE
271 IOFF = IZERO
272 DO 40 I = 1, IZERO - 1
273 A( IOFF ) = ZERO
274 IOFF = IOFF + LDA
275 40 CONTINUE
276 IOFF = IOFF - IZERO
277 DO 50 I = IZERO, N
278 A( IOFF+I ) = ZERO
279 50 CONTINUE
280 END IF
281 ELSE
282 IZERO = 0
283 END IF
284 *
285 * Save a copy of the matrix A in ASAV.
286 *
287 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
288 *
289 * Compute the condition number of A (RCONDC).
290 *
291 IF( ZEROT ) THEN
292 RCONDC = ZERO
293 ELSE
294 *
295 * Compute the 1-norm of A.
296 *
297 ANORM = DLANSY( '1', UPLO, N, A, LDA,
298 + D_WORK_DLANSY )
299 *
300 * Factor the matrix A.
301 *
302 CALL DPOTRF( UPLO, N, A, LDA, INFO )
303 *
304 * Form the inverse of A.
305 *
306 CALL DPOTRI( UPLO, N, A, LDA, INFO )
307 *
308 * Compute the 1-norm condition number of A.
309 *
310 AINVNM = DLANSY( '1', UPLO, N, A, LDA,
311 + D_WORK_DLANSY )
312 RCONDC = ( ONE / ANORM ) / AINVNM
313 *
314 * Restore the matrix A.
315 *
316 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
317 *
318 END IF
319 *
320 * Form an exact solution and set the right hand side.
321 *
322 SRNAMT = 'DLARHS'
323 CALL DLARHS( 'DPO', 'N', UPLO, ' ', N, N, KL, KU,
324 + NRHS, A, LDA, XACT, LDA, B, LDA,
325 + ISEED, INFO )
326 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
327 *
328 * Compute the L*L' or U'*U factorization of the
329 * matrix and solve the system.
330 *
331 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
332 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDB )
333 *
334 SRNAMT = 'DTRTTF'
335 CALL DTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO )
336 SRNAMT = 'DPFTRF'
337 CALL DPFTRF( CFORM, UPLO, N, ARF, INFO )
338 *
339 * Check error code from DPFTRF.
340 *
341 IF( INFO.NE.IZERO ) THEN
342 *
343 * LANGOU: there is a small hick here: IZERO should
344 * always be INFO however if INFO is ZERO, ALAERH does not
345 * complain.
346 *
347 CALL ALAERH( 'DPF', 'DPFSV ', INFO, IZERO,
348 + UPLO, N, N, -1, -1, NRHS, IIT,
349 + NFAIL, NERRS, NOUT )
350 GO TO 100
351 END IF
352 *
353 * Skip the tests if INFO is not 0.
354 *
355 IF( INFO.NE.0 ) THEN
356 GO TO 100
357 END IF
358 *
359 SRNAMT = 'DPFTRS'
360 CALL DPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB,
361 + INFO )
362 *
363 SRNAMT = 'DTFTTR'
364 CALL DTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO )
365 *
366 * Reconstruct matrix from factors and compute
367 * residual.
368 *
369 CALL DLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA )
370 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA,
371 + D_WORK_DPOT01, RESULT( 1 ) )
372 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
373 *
374 * Form the inverse and compute the residual.
375 *
376 IF(MOD(N,2).EQ.0)THEN
377 CALL DLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV,
378 + N+1 )
379 ELSE
380 CALL DLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV,
381 + N )
382 END IF
383 *
384 SRNAMT = 'DPFTRI'
385 CALL DPFTRI( CFORM, UPLO, N, ARFINV , INFO )
386 *
387 SRNAMT = 'DTFTTR'
388 CALL DTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA,
389 + INFO )
390 *
391 * Check error code from DPFTRI.
392 *
393 IF( INFO.NE.0 )
394 + CALL ALAERH( 'DPO', 'DPFTRI', INFO, 0, UPLO, N,
395 + N, -1, -1, -1, IMAT, NFAIL, NERRS,
396 + NOUT )
397 *
398 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA,
399 + D_TEMP_DPOT03, LDA, D_WORK_DPOT03,
400 + RCONDC, RESULT( 2 ) )
401 *
402 * Compute residual of the computed solution.
403 *
404 CALL DLACPY( 'Full', N, NRHS, B, LDA,
405 + D_TEMP_DPOT02, LDA )
406 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
407 + D_TEMP_DPOT02, LDA, D_WORK_DPOT02,
408 + RESULT( 3 ) )
409 *
410 * Check solution from generated exact solution.
411
412 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
413 + RESULT( 4 ) )
414 NT = 4
415 *
416 * Print information about the tests that did not
417 * pass the threshold.
418 *
419 DO 60 K = 1, NT
420 IF( RESULT( K ).GE.THRESH ) THEN
421 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
422 + CALL ALADHD( NOUT, 'DPF' )
423 WRITE( NOUT, FMT = 9999 )'DPFSV ', UPLO,
424 + N, IIT, K, RESULT( K )
425 NFAIL = NFAIL + 1
426 END IF
427 60 CONTINUE
428 NRUN = NRUN + NT
429 100 CONTINUE
430 110 CONTINUE
431 120 CONTINUE
432 980 CONTINUE
433 130 CONTINUE
434 *
435 * Print a summary of the results.
436 *
437 CALL ALASVM( 'DPF', NOUT, NFAIL, NRUN, NERRS )
438 *
439 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
440 + ', test(', I1, ')=', G12.5 )
441 *
442 RETURN
443 *
444 * End of DDRVRFP
445 *
446 END
2 + THRESH, A, ASAV, AFAC, AINV, B,
3 + BSAV, XACT, X, ARF, ARFINV,
4 + D_WORK_DLATMS, D_WORK_DPOT01, D_TEMP_DPOT02,
5 + D_TEMP_DPOT03, D_WORK_DLANSY,
6 + D_WORK_DPOT02, D_WORK_DPOT03 )
7 *
8 * -- LAPACK test routine (version 3.2.0) --
9 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
10 * November 2008
11 *
12 * .. Scalar Arguments ..
13 INTEGER NN, NNS, NNT, NOUT
14 DOUBLE PRECISION THRESH
15 * ..
16 * .. Array Arguments ..
17 INTEGER NVAL( NN ), NSVAL( NNS ), NTVAL( NNT )
18 DOUBLE PRECISION A( * )
19 DOUBLE PRECISION AINV( * )
20 DOUBLE PRECISION ASAV( * )
21 DOUBLE PRECISION B( * )
22 DOUBLE PRECISION BSAV( * )
23 DOUBLE PRECISION AFAC( * )
24 DOUBLE PRECISION ARF( * )
25 DOUBLE PRECISION ARFINV( * )
26 DOUBLE PRECISION XACT( * )
27 DOUBLE PRECISION X( * )
28 DOUBLE PRECISION D_WORK_DLATMS( * )
29 DOUBLE PRECISION D_WORK_DPOT01( * )
30 DOUBLE PRECISION D_TEMP_DPOT02( * )
31 DOUBLE PRECISION D_TEMP_DPOT03( * )
32 DOUBLE PRECISION D_WORK_DLANSY( * )
33 DOUBLE PRECISION D_WORK_DPOT02( * )
34 DOUBLE PRECISION D_WORK_DPOT03( * )
35 * ..
36 *
37 * Purpose
38 * =======
39 *
40 * DDRVRFP tests the LAPACK RFP routines:
41 * DPFTRF, DPFTRS, and DPFTRI.
42 *
43 * This testing routine follow the same tests as DDRVPO (test for the full
44 * format Symmetric Positive Definite solver).
45 *
46 * The tests are performed in Full Format, convertion back and forth from
47 * full format to RFP format are performed using the routines DTRTTF and
48 * DTFTTR.
49 *
50 * First, a specific matrix A of size N is created. There is nine types of
51 * different matrixes possible.
52 * 1. Diagonal 6. Random, CNDNUM = sqrt(0.1/EPS)
53 * 2. Random, CNDNUM = 2 7. Random, CNDNUM = 0.1/EPS
54 * *3. First row and column zero 8. Scaled near underflow
55 * *4. Last row and column zero 9. Scaled near overflow
56 * *5. Middle row and column zero
57 * (* - tests error exits from DPFTRF, no test ratios are computed)
58 * A solution XACT of size N-by-NRHS is created and the associated right
59 * hand side B as well. Then DPFTRF is called to compute L (or U), the
60 * Cholesky factor of A. Then L (or U) is used to solve the linear system
61 * of equations AX = B. This gives X. Then L (or U) is used to compute the
62 * inverse of A, AINV. The following four tests are then performed:
63 * (1) norm( L*L' - A ) / ( N * norm(A) * EPS ) or
64 * norm( U'*U - A ) / ( N * norm(A) * EPS ),
65 * (2) norm(B - A*X) / ( norm(A) * norm(X) * EPS ),
66 * (3) norm( I - A*AINV ) / ( N * norm(A) * norm(AINV) * EPS ),
67 * (4) ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),
68 * where EPS is the machine precision, RCOND the condition number of A, and
69 * norm( . ) the 1-norm for (1,2,3) and the inf-norm for (4).
70 * Errors occur when INFO parameter is not as expected. Failures occur when
71 * a test ratios is greater than THRES.
72 *
73 * Arguments
74 * =========
75 *
76 * NOUT (input) INTEGER
77 * The unit number for output.
78 *
79 * NN (input) INTEGER
80 * The number of values of N contained in the vector NVAL.
81 *
82 * NVAL (input) INTEGER array, dimension (NN)
83 * The values of the matrix dimension N.
84 *
85 * NNS (input) INTEGER
86 * The number of values of NRHS contained in the vector NSVAL.
87 *
88 * NSVAL (input) INTEGER array, dimension (NNS)
89 * The values of the number of right-hand sides NRHS.
90 *
91 * NNT (input) INTEGER
92 * The number of values of MATRIX TYPE contained in the vector NTVAL.
93 *
94 * NTVAL (input) INTEGER array, dimension (NNT)
95 * The values of matrix type (between 0 and 9 for PO/PP/PF matrices).
96 *
97 * THRESH (input) DOUBLE PRECISION
98 * The threshold value for the test ratios. A result is
99 * included in the output file if RESULT >= THRESH. To have
100 * every test ratio printed, use THRESH = 0.
101 *
102 * A (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
103 *
104 * ASAV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
105 *
106 * AFAC (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
107 *
108 * AINV (workspace) DOUBLE PRECISION array, dimension (NMAX*NMAX)
109 *
110 * B (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
111 *
112 * BSAV (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
113 *
114 * XACT (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
115 *
116 * X (workspace) DOUBLE PRECISION array, dimension (NMAX*MAXRHS)
117 *
118 * ARF (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
119 *
120 * ARFINV (workspace) DOUBLE PRECISION array, dimension ((NMAX*(NMAX+1))/2)
121 *
122 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( 3*NMAX )
123 *
124 * D_WORK_DPOT01 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
125 *
126 * D_TEMP_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX*MAXRHS )
127 *
128 * D_TEMP_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX*NMAX )
129 *
130 * D_WORK_DLATMS (workspace) DOUBLE PRECISION array, dimension ( NMAX )
131 *
132 * D_WORK_DLANSY (workspace) DOUBLE PRECISION array, dimension ( NMAX )
133 *
134 * D_WORK_DPOT02 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
135 *
136 * D_WORK_DPOT03 (workspace) DOUBLE PRECISION array, dimension ( NMAX )
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141 DOUBLE PRECISION ONE, ZERO
142 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
143 INTEGER NTESTS
144 PARAMETER ( NTESTS = 4 )
145 * ..
146 * .. Local Scalars ..
147 LOGICAL ZEROT
148 INTEGER I, INFO, IUPLO, LDA, LDB, IMAT, NERRS, NFAIL,
149 + NRHS, NRUN, IZERO, IOFF, K, NT, N, IFORM, IIN,
150 + IIT, IIS
151 CHARACTER DIST, CTYPE, UPLO, CFORM
152 INTEGER KL, KU, MODE
153 DOUBLE PRECISION ANORM, AINVNM, CNDNUM, RCONDC
154 * ..
155 * .. Local Arrays ..
156 CHARACTER UPLOS( 2 ), FORMS( 2 )
157 INTEGER ISEED( 4 ), ISEEDY( 4 )
158 DOUBLE PRECISION RESULT( NTESTS )
159 * ..
160 * .. External Functions ..
161 DOUBLE PRECISION DLANSY
162 EXTERNAL DLANSY
163 * ..
164 * .. External Subroutines ..
165 EXTERNAL ALADHD, ALAERH, ALASVM, DGET04, DTFTTR, DLACPY,
166 + DLARHS, DLATB4, DLATMS, DPFTRI, DPFTRF, DPFTRS,
167 + DPOT01, DPOT02, DPOT03, DPOTRI, DPOTRF, DTRTTF
168 * ..
169 * .. Scalars in Common ..
170 CHARACTER*32 SRNAMT
171 * ..
172 * .. Common blocks ..
173 COMMON / SRNAMC / SRNAMT
174 * ..
175 * .. Data statements ..
176 DATA ISEEDY / 1988, 1989, 1990, 1991 /
177 DATA UPLOS / 'U', 'L' /
178 DATA FORMS / 'N', 'T' /
179 * ..
180 * .. Executable Statements ..
181 *
182 * Initialize constants and the random number seed.
183 *
184 NRUN = 0
185 NFAIL = 0
186 NERRS = 0
187 DO 10 I = 1, 4
188 ISEED( I ) = ISEEDY( I )
189 10 CONTINUE
190 *
191 DO 130 IIN = 1, NN
192 *
193 N = NVAL( IIN )
194 LDA = MAX( N, 1 )
195 LDB = MAX( N, 1 )
196 *
197 DO 980 IIS = 1, NNS
198 *
199 NRHS = NSVAL( IIS )
200 *
201 DO 120 IIT = 1, NNT
202 *
203 IMAT = NTVAL( IIT )
204 *
205 * If N.EQ.0, only consider the first type
206 *
207 IF( N.EQ.0 .AND. IIT.GT.1 ) GO TO 120
208 *
209 * Skip types 3, 4, or 5 if the matrix size is too small.
210 *
211 IF( IMAT.EQ.4 .AND. N.LE.1 ) GO TO 120
212 IF( IMAT.EQ.5 .AND. N.LE.2 ) GO TO 120
213 *
214 * Do first for UPLO = 'U', then for UPLO = 'L'
215 *
216 DO 110 IUPLO = 1, 2
217 UPLO = UPLOS( IUPLO )
218 *
219 * Do first for CFORM = 'N', then for CFORM = 'C'
220 *
221 DO 100 IFORM = 1, 2
222 CFORM = FORMS( IFORM )
223 *
224 * Set up parameters with DLATB4 and generate a test
225 * matrix with DLATMS.
226 *
227 CALL DLATB4( 'DPO', IMAT, N, N, CTYPE, KL, KU,
228 + ANORM, MODE, CNDNUM, DIST )
229 *
230 SRNAMT = 'DLATMS'
231 CALL DLATMS( N, N, DIST, ISEED, CTYPE,
232 + D_WORK_DLATMS,
233 + MODE, CNDNUM, ANORM, KL, KU, UPLO, A,
234 + LDA, D_WORK_DLATMS, INFO )
235 *
236 * Check error code from DLATMS.
237 *
238 IF( INFO.NE.0 ) THEN
239 CALL ALAERH( 'DPF', 'DLATMS', INFO, 0, UPLO, N,
240 + N, -1, -1, -1, IIT, NFAIL, NERRS,
241 + NOUT )
242 GO TO 100
243 END IF
244 *
245 * For types 3-5, zero one row and column of the matrix to
246 * test that INFO is returned correctly.
247 *
248 ZEROT = IMAT.GE.3 .AND. IMAT.LE.5
249 IF( ZEROT ) THEN
250 IF( IIT.EQ.3 ) THEN
251 IZERO = 1
252 ELSE IF( IIT.EQ.4 ) THEN
253 IZERO = N
254 ELSE
255 IZERO = N / 2 + 1
256 END IF
257 IOFF = ( IZERO-1 )*LDA
258 *
259 * Set row and column IZERO of A to 0.
260 *
261 IF( IUPLO.EQ.1 ) THEN
262 DO 20 I = 1, IZERO - 1
263 A( IOFF+I ) = ZERO
264 20 CONTINUE
265 IOFF = IOFF + IZERO
266 DO 30 I = IZERO, N
267 A( IOFF ) = ZERO
268 IOFF = IOFF + LDA
269 30 CONTINUE
270 ELSE
271 IOFF = IZERO
272 DO 40 I = 1, IZERO - 1
273 A( IOFF ) = ZERO
274 IOFF = IOFF + LDA
275 40 CONTINUE
276 IOFF = IOFF - IZERO
277 DO 50 I = IZERO, N
278 A( IOFF+I ) = ZERO
279 50 CONTINUE
280 END IF
281 ELSE
282 IZERO = 0
283 END IF
284 *
285 * Save a copy of the matrix A in ASAV.
286 *
287 CALL DLACPY( UPLO, N, N, A, LDA, ASAV, LDA )
288 *
289 * Compute the condition number of A (RCONDC).
290 *
291 IF( ZEROT ) THEN
292 RCONDC = ZERO
293 ELSE
294 *
295 * Compute the 1-norm of A.
296 *
297 ANORM = DLANSY( '1', UPLO, N, A, LDA,
298 + D_WORK_DLANSY )
299 *
300 * Factor the matrix A.
301 *
302 CALL DPOTRF( UPLO, N, A, LDA, INFO )
303 *
304 * Form the inverse of A.
305 *
306 CALL DPOTRI( UPLO, N, A, LDA, INFO )
307 *
308 * Compute the 1-norm condition number of A.
309 *
310 AINVNM = DLANSY( '1', UPLO, N, A, LDA,
311 + D_WORK_DLANSY )
312 RCONDC = ( ONE / ANORM ) / AINVNM
313 *
314 * Restore the matrix A.
315 *
316 CALL DLACPY( UPLO, N, N, ASAV, LDA, A, LDA )
317 *
318 END IF
319 *
320 * Form an exact solution and set the right hand side.
321 *
322 SRNAMT = 'DLARHS'
323 CALL DLARHS( 'DPO', 'N', UPLO, ' ', N, N, KL, KU,
324 + NRHS, A, LDA, XACT, LDA, B, LDA,
325 + ISEED, INFO )
326 CALL DLACPY( 'Full', N, NRHS, B, LDA, BSAV, LDA )
327 *
328 * Compute the L*L' or U'*U factorization of the
329 * matrix and solve the system.
330 *
331 CALL DLACPY( UPLO, N, N, A, LDA, AFAC, LDA )
332 CALL DLACPY( 'Full', N, NRHS, B, LDB, X, LDB )
333 *
334 SRNAMT = 'DTRTTF'
335 CALL DTRTTF( CFORM, UPLO, N, AFAC, LDA, ARF, INFO )
336 SRNAMT = 'DPFTRF'
337 CALL DPFTRF( CFORM, UPLO, N, ARF, INFO )
338 *
339 * Check error code from DPFTRF.
340 *
341 IF( INFO.NE.IZERO ) THEN
342 *
343 * LANGOU: there is a small hick here: IZERO should
344 * always be INFO however if INFO is ZERO, ALAERH does not
345 * complain.
346 *
347 CALL ALAERH( 'DPF', 'DPFSV ', INFO, IZERO,
348 + UPLO, N, N, -1, -1, NRHS, IIT,
349 + NFAIL, NERRS, NOUT )
350 GO TO 100
351 END IF
352 *
353 * Skip the tests if INFO is not 0.
354 *
355 IF( INFO.NE.0 ) THEN
356 GO TO 100
357 END IF
358 *
359 SRNAMT = 'DPFTRS'
360 CALL DPFTRS( CFORM, UPLO, N, NRHS, ARF, X, LDB,
361 + INFO )
362 *
363 SRNAMT = 'DTFTTR'
364 CALL DTFTTR( CFORM, UPLO, N, ARF, AFAC, LDA, INFO )
365 *
366 * Reconstruct matrix from factors and compute
367 * residual.
368 *
369 CALL DLACPY( UPLO, N, N, AFAC, LDA, ASAV, LDA )
370 CALL DPOT01( UPLO, N, A, LDA, AFAC, LDA,
371 + D_WORK_DPOT01, RESULT( 1 ) )
372 CALL DLACPY( UPLO, N, N, ASAV, LDA, AFAC, LDA )
373 *
374 * Form the inverse and compute the residual.
375 *
376 IF(MOD(N,2).EQ.0)THEN
377 CALL DLACPY( 'A', N+1, N/2, ARF, N+1, ARFINV,
378 + N+1 )
379 ELSE
380 CALL DLACPY( 'A', N, (N+1)/2, ARF, N, ARFINV,
381 + N )
382 END IF
383 *
384 SRNAMT = 'DPFTRI'
385 CALL DPFTRI( CFORM, UPLO, N, ARFINV , INFO )
386 *
387 SRNAMT = 'DTFTTR'
388 CALL DTFTTR( CFORM, UPLO, N, ARFINV, AINV, LDA,
389 + INFO )
390 *
391 * Check error code from DPFTRI.
392 *
393 IF( INFO.NE.0 )
394 + CALL ALAERH( 'DPO', 'DPFTRI', INFO, 0, UPLO, N,
395 + N, -1, -1, -1, IMAT, NFAIL, NERRS,
396 + NOUT )
397 *
398 CALL DPOT03( UPLO, N, A, LDA, AINV, LDA,
399 + D_TEMP_DPOT03, LDA, D_WORK_DPOT03,
400 + RCONDC, RESULT( 2 ) )
401 *
402 * Compute residual of the computed solution.
403 *
404 CALL DLACPY( 'Full', N, NRHS, B, LDA,
405 + D_TEMP_DPOT02, LDA )
406 CALL DPOT02( UPLO, N, NRHS, A, LDA, X, LDA,
407 + D_TEMP_DPOT02, LDA, D_WORK_DPOT02,
408 + RESULT( 3 ) )
409 *
410 * Check solution from generated exact solution.
411
412 CALL DGET04( N, NRHS, X, LDA, XACT, LDA, RCONDC,
413 + RESULT( 4 ) )
414 NT = 4
415 *
416 * Print information about the tests that did not
417 * pass the threshold.
418 *
419 DO 60 K = 1, NT
420 IF( RESULT( K ).GE.THRESH ) THEN
421 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
422 + CALL ALADHD( NOUT, 'DPF' )
423 WRITE( NOUT, FMT = 9999 )'DPFSV ', UPLO,
424 + N, IIT, K, RESULT( K )
425 NFAIL = NFAIL + 1
426 END IF
427 60 CONTINUE
428 NRUN = NRUN + NT
429 100 CONTINUE
430 110 CONTINUE
431 120 CONTINUE
432 980 CONTINUE
433 130 CONTINUE
434 *
435 * Print a summary of the results.
436 *
437 CALL ALASVM( 'DPF', NOUT, NFAIL, NRUN, NERRS )
438 *
439 9999 FORMAT( 1X, A6, ', UPLO=''', A1, ''', N =', I5, ', type ', I1,
440 + ', test(', I1, ')=', G12.5 )
441 *
442 RETURN
443 *
444 * End of DDRVRFP
445 *
446 END