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