LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
debchvxx.f
Go to the documentation of this file.
1*> \brief \b DEBCHVXX
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 DEBCHVXX( THRESH, PATH )
12*
13* .. Scalar Arguments ..
14* DOUBLE PRECISION THRESH
15* CHARACTER*3 PATH
16* ..
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> DEBCHVXX will run D**SVXX on a series of Hilbert matrices and then
25*> compare the error bounds returned by D**SVXX to see if the returned
26*> answer indeed falls within those bounds.
27*>
28*> Eight test ratios will be computed. The tests will pass if they are .LT.
29*> THRESH. There are two cases that are determined by 1 / (SQRT( N ) * EPS).
30*> If that value is .LE. to the component wise reciprocal condition number,
31*> it uses the guaranteed case, other wise it uses the unguaranteed case.
32*>
33*> Test ratios:
34*> Let Xc be X_computed and Xt be X_truth.
35*> The norm used is the infinity norm.
36*>
37*> Let A be the guaranteed case and B be the unguaranteed case.
38*>
39*> 1. Normwise guaranteed forward error bound.
40*> A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
41*> ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
42*> If these conditions are met, the test ratio is set to be
43*> ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
44*> B: For this case, CGESVXX should just return 1. If it is less than
45*> one, treat it the same as in 1A. Otherwise it fails. (Set test
46*> ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
47*>
48*> 2. Componentwise guaranteed forward error bound.
49*> A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
50*> for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
51*> If these conditions are met, the test ratio is set to be
52*> ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10). Otherwise it is 1/EPS.
53*> B: Same as normwise test ratio.
54*>
55*> 3. Backwards error.
56*> A: The test ratio is set to BERR/EPS.
57*> B: Same test ratio.
58*>
59*> 4. Reciprocal condition number.
60*> A: A condition number is computed with Xt and compared with the one
61*> returned from CGESVXX. Let RCONDc be the RCOND returned by D**SVXX
62*> and RCONDt be the RCOND from the truth value. Test ratio is set to
63*> MAX(RCONDc/RCONDt, RCONDt/RCONDc).
64*> B: Test ratio is set to 1 / (EPS * RCONDc).
65*>
66*> 5. Reciprocal normwise condition number.
67*> A: The test ratio is set to
68*> MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
69*> B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
70*>
71*> 6. Reciprocal componentwise condition number.
72*> A: Test ratio is set to
73*> MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
74*> B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
75*>
76*> .. Parameters ..
77*> NMAX is determined by the largest number in the inverse of the hilbert
78*> matrix. Precision is exhausted when the largest entry in it is greater
79*> than 2 to the power of the number of bits in the fraction of the data
80*> type used plus one, which is 24 for single precision.
81*> NMAX should be 6 for single and 11 for double.
82*> \endverbatim
83*
84* Authors:
85* ========
86*
87*> \author Univ. of Tennessee
88*> \author Univ. of California Berkeley
89*> \author Univ. of Colorado Denver
90*> \author NAG Ltd.
91*
92*> \ingroup double_lin
93*
94* =====================================================================
95 SUBROUTINE debchvxx( THRESH, PATH )
96 IMPLICIT NONE
97* .. Scalar Arguments ..
98 DOUBLE PRECISION THRESH
99 CHARACTER*3 PATH
100
101 INTEGER NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
102 parameter(nmax = 10, nparams = 2, nerrbnd = 3,
103 $ ntests = 6)
104
105* .. Local Scalars ..
106 INTEGER N, NRHS, INFO, I ,J, k, NFAIL, LDA,
107 $ N_AUX_TESTS, LDAB, LDAFB
108 CHARACTER FACT, TRANS, UPLO, EQUED
109 CHARACTER*2 C2
110 CHARACTER(3) NGUAR, CGUAR
111 LOGICAL printed_guide
112 DOUBLE PRECISION NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
113 $ RNORM, RINORM, SUMR, SUMRI, EPS,
114 $ BERR(NMAX), RPVGRW, ORCOND,
115 $ CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
116 $ CWISE_RCOND, NWISE_RCOND,
117 $ CONDTHRESH, ERRTHRESH
118
119* .. Local Arrays ..
120 DOUBLE PRECISION TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
121 $ S(NMAX),R(NMAX),C(NMAX), DIFF(NMAX, NMAX),
122 $ ERRBND_N(NMAX*3), ERRBND_C(NMAX*3),
123 $ A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
124 $ AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
125 $ ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
126 $ AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX ),
127 $ WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
128 $ ACOPY(NMAX, NMAX)
129 INTEGER IPIV(NMAX), IWORK(3*NMAX)
130
131* .. External Functions ..
132 DOUBLE PRECISION DLAMCH
133
134* .. External Subroutines ..
135 EXTERNAL dlahilb, dgesvxx, dposvxx, dsysvxx,
137 LOGICAL LSAMEN
138
139* .. Intrinsic Functions ..
140 INTRINSIC sqrt, max, abs, dble
141
142* .. Parameters ..
143 INTEGER NWISE_I, CWISE_I
144 parameter(nwise_i = 1, cwise_i = 1)
145 INTEGER BND_I, COND_I
146 parameter(bnd_i = 2, cond_i = 3)
147
148* Create the loop to test out the Hilbert matrices
149
150 fact = 'E'
151 uplo = 'U'
152 trans = 'N'
153 equed = 'N'
154 eps = dlamch('Epsilon')
155 nfail = 0
156 n_aux_tests = 0
157 lda = nmax
158 ldab = (nmax-1)+(nmax-1)+1
159 ldafb = 2*(nmax-1)+(nmax-1)+1
160 c2 = path( 2: 3 )
161
162* Main loop to test the different Hilbert Matrices.
163
164 printed_guide = .false.
165
166 DO n = 1 , nmax
167 params(1) = -1
168 params(2) = -1
169
170 kl = n-1
171 ku = n-1
172 nrhs = n
173 m = max(sqrt(dble(n)), 10.0d+0)
174
175* Generate the Hilbert matrix, its inverse, and the
176* right hand side, all scaled by the LCM(1,..,2N-1).
177 CALL dlahilb(n, n, a, lda, invhilb, lda, b, lda, work, info)
178
179* Copy A into ACOPY.
180 CALL dlacpy('ALL', n, n, a, nmax, acopy, nmax)
181
182* Store A in band format for GB tests
183 DO j = 1, n
184 DO i = 1, kl+ku+1
185 ab( i, j ) = 0.0d+0
186 END DO
187 END DO
188 DO j = 1, n
189 DO i = max( 1, j-ku ), min( n, j+kl )
190 ab( ku+1+i-j, j ) = a( i, j )
191 END DO
192 END DO
193
194* Copy AB into ABCOPY.
195 DO j = 1, n
196 DO i = 1, kl+ku+1
197 abcopy( i, j ) = 0.0d+0
198 END DO
199 END DO
200 CALL dlacpy('ALL', kl+ku+1, n, ab, ldab, abcopy, ldab)
201
202* Call D**SVXX with default PARAMS and N_ERR_BND = 3.
203 IF ( lsamen( 2, c2, 'SY' ) ) THEN
204 CALL dsysvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
205 $ ipiv, equed, s, b, lda, x, lda, orcond,
206 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
207 $ params, work, iwork, info)
208 ELSE IF ( lsamen( 2, c2, 'PO' ) ) THEN
209 CALL dposvxx(fact, uplo, n, nrhs, acopy, lda, af, lda,
210 $ equed, s, b, lda, x, lda, orcond,
211 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
212 $ params, work, iwork, info)
213 ELSE IF ( lsamen( 2, c2, 'GB' ) ) THEN
214 CALL dgbsvxx(fact, trans, n, kl, ku, nrhs, abcopy,
215 $ ldab, afb, ldafb, ipiv, equed, r, c, b,
216 $ lda, x, lda, orcond, rpvgrw, berr, nerrbnd,
217 $ errbnd_n, errbnd_c, nparams, params, work, iwork,
218 $ info)
219 ELSE
220 CALL dgesvxx(fact, trans, n, nrhs, acopy, lda, af, lda,
221 $ ipiv, equed, r, c, b, lda, x, lda, orcond,
222 $ rpvgrw, berr, nerrbnd, errbnd_n, errbnd_c, nparams,
223 $ params, work, iwork, info)
224 END IF
225
226 n_aux_tests = n_aux_tests + 1
227 IF (orcond .LT. eps) THEN
228! Either factorization failed or the matrix is flagged, and 1 <=
229! INFO <= N+1. We don't decide based on rcond anymore.
230! IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
231! NFAIL = NFAIL + 1
232! WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
233! END IF
234 ELSE
235! Either everything succeeded (INFO == 0) or some solution failed
236! to converge (INFO > N+1).
237 IF (info .GT. 0 .AND. info .LE. n+1) THEN
238 nfail = nfail + 1
239 WRITE (*, fmt=8000) c2, n, info, orcond, rcond
240 END IF
241 END IF
242
243* Calculating the difference between D**SVXX's X and the true X.
244 DO i = 1,n
245 DO j =1,nrhs
246 diff(i,j) = x(i,j) - invhilb(i,j)
247 END DO
248 END DO
249
250* Calculating the RCOND
251 rnorm = 0.0d+0
252 rinorm = 0.0d+0
253 IF ( lsamen( 2, c2, 'PO' ) .OR. lsamen( 2, c2, 'SY' ) ) THEN
254 DO i = 1, n
255 sumr = 0.0d+0
256 sumri = 0.0d+0
257 DO j = 1, n
258 sumr = sumr + s(i) * abs(a(i,j)) * s(j)
259 sumri = sumri + abs(invhilb(i, j)) / (s(j) * s(i))
260
261 END DO
262 rnorm = max(rnorm,sumr)
263 rinorm = max(rinorm,sumri)
264 END DO
265 ELSE IF ( lsamen( 2, c2, 'GE' ) .OR. lsamen( 2, c2, 'GB' ) )
266 $ THEN
267 DO i = 1, n
268 sumr = 0.0d+0
269 sumri = 0.0d+0
270 DO j = 1, n
271 sumr = sumr + r(i) * abs(a(i,j)) * c(j)
272 sumri = sumri + abs(invhilb(i, j)) / (r(j) * c(i))
273 END DO
274 rnorm = max(rnorm,sumr)
275 rinorm = max(rinorm,sumri)
276 END DO
277 END IF
278
279 rnorm = rnorm / abs(a(1, 1))
280 rcond = 1.0d+0/(rnorm * rinorm)
281
282* Calculating the R for normwise rcond.
283 DO i = 1, n
284 rinv(i) = 0.0d+0
285 END DO
286 DO j = 1, n
287 DO i = 1, n
288 rinv(i) = rinv(i) + abs(a(i,j))
289 END DO
290 END DO
291
292* Calculating the Normwise rcond.
293 rinorm = 0.0d+0
294 DO i = 1, n
295 sumri = 0.0d+0
296 DO j = 1, n
297 sumri = sumri + abs(invhilb(i,j) * rinv(j))
298 END DO
299 rinorm = max(rinorm, sumri)
300 END DO
301
302! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
303! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
304 ncond = abs(a(1,1)) / rinorm
305
306 condthresh = m * eps
307 errthresh = m * eps
308
309 DO k = 1, nrhs
310 normt = 0.0d+0
311 normdif = 0.0d+0
312 cwise_err = 0.0d+0
313 DO i = 1, n
314 normt = max(abs(invhilb(i, k)), normt)
315 normdif = max(abs(x(i,k) - invhilb(i,k)), normdif)
316 IF (invhilb(i,k) .NE. 0.0d+0) THEN
317 cwise_err = max(abs(x(i,k) - invhilb(i,k))
318 $ /abs(invhilb(i,k)), cwise_err)
319 ELSE IF (x(i, k) .NE. 0.0d+0) THEN
320 cwise_err = dlamch('OVERFLOW')
321 END IF
322 END DO
323 IF (normt .NE. 0.0d+0) THEN
324 nwise_err = normdif / normt
325 ELSE IF (normdif .NE. 0.0d+0) THEN
326 nwise_err = dlamch('OVERFLOW')
327 ELSE
328 nwise_err = 0.0d+0
329 ENDIF
330
331 DO i = 1, n
332 rinv(i) = 0.0d+0
333 END DO
334 DO j = 1, n
335 DO i = 1, n
336 rinv(i) = rinv(i) + abs(a(i, j) * invhilb(j, k))
337 END DO
338 END DO
339 rinorm = 0.0d+0
340 DO i = 1, n
341 sumri = 0.0d+0
342 DO j = 1, n
343 sumri = sumri
344 $ + abs(invhilb(i, j) * rinv(j) / invhilb(i, k))
345 END DO
346 rinorm = max(rinorm, sumri)
347 END DO
348! invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
349! by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
350 ccond = abs(a(1,1))/rinorm
351
352! Forward error bound tests
353 nwise_bnd = errbnd_n(k + (bnd_i-1)*nrhs)
354 cwise_bnd = errbnd_c(k + (bnd_i-1)*nrhs)
355 nwise_rcond = errbnd_n(k + (cond_i-1)*nrhs)
356 cwise_rcond = errbnd_c(k + (cond_i-1)*nrhs)
357! write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
358! $ condthresh, ncond.ge.condthresh
359! write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
360 IF (ncond .GE. condthresh) THEN
361 nguar = 'YES'
362 IF (nwise_bnd .GT. errthresh) THEN
363 tstrat(1) = 1/(2.0d+0*eps)
364 ELSE
365 IF (nwise_bnd .NE. 0.0d+0) THEN
366 tstrat(1) = nwise_err / nwise_bnd
367 ELSE IF (nwise_err .NE. 0.0d+0) THEN
368 tstrat(1) = 1/(16.0*eps)
369 ELSE
370 tstrat(1) = 0.0d+0
371 END IF
372 IF (tstrat(1) .GT. 1.0d+0) THEN
373 tstrat(1) = 1/(4.0d+0*eps)
374 END IF
375 END IF
376 ELSE
377 nguar = 'NO'
378 IF (nwise_bnd .LT. 1.0d+0) THEN
379 tstrat(1) = 1/(8.0d+0*eps)
380 ELSE
381 tstrat(1) = 1.0d+0
382 END IF
383 END IF
384! write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
385! $ condthresh, ccond.ge.condthresh
386! write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
387 IF (ccond .GE. condthresh) THEN
388 cguar = 'YES'
389 IF (cwise_bnd .GT. errthresh) THEN
390 tstrat(2) = 1/(2.0d+0*eps)
391 ELSE
392 IF (cwise_bnd .NE. 0.0d+0) THEN
393 tstrat(2) = cwise_err / cwise_bnd
394 ELSE IF (cwise_err .NE. 0.0d+0) THEN
395 tstrat(2) = 1/(16.0d+0*eps)
396 ELSE
397 tstrat(2) = 0.0d+0
398 END IF
399 IF (tstrat(2) .GT. 1.0d+0) tstrat(2) = 1/(4.0d+0*eps)
400 END IF
401 ELSE
402 cguar = 'NO'
403 IF (cwise_bnd .LT. 1.0d+0) THEN
404 tstrat(2) = 1/(8.0d+0*eps)
405 ELSE
406 tstrat(2) = 1.0d+0
407 END IF
408 END IF
409
410! Backwards error test
411 tstrat(3) = berr(k)/eps
412
413! Condition number tests
414 tstrat(4) = rcond / orcond
415 IF (rcond .GE. condthresh .AND. tstrat(4) .LT. 1.0d+0)
416 $ tstrat(4) = 1.0d+0 / tstrat(4)
417
418 tstrat(5) = ncond / nwise_rcond
419 IF (ncond .GE. condthresh .AND. tstrat(5) .LT. 1.0d+0)
420 $ tstrat(5) = 1.0d+0 / tstrat(5)
421
422 tstrat(6) = ccond / nwise_rcond
423 IF (ccond .GE. condthresh .AND. tstrat(6) .LT. 1.0d+0)
424 $ tstrat(6) = 1.0d+0 / tstrat(6)
425
426 DO i = 1, ntests
427 IF (tstrat(i) .GT. thresh) THEN
428 IF (.NOT.printed_guide) THEN
429 WRITE(*,*)
430 WRITE( *, 9996) 1
431 WRITE( *, 9995) 2
432 WRITE( *, 9994) 3
433 WRITE( *, 9993) 4
434 WRITE( *, 9992) 5
435 WRITE( *, 9991) 6
436 WRITE( *, 9990) 7
437 WRITE( *, 9989) 8
438 WRITE(*,*)
439 printed_guide = .true.
440 END IF
441 WRITE( *, 9999) c2, n, k, nguar, cguar, i, tstrat(i)
442 nfail = nfail + 1
443 END IF
444 END DO
445 END DO
446
447c$$$ WRITE(*,*)
448c$$$ WRITE(*,*) 'Normwise Error Bounds'
449c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
450c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
451c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
452c$$$ WRITE(*,*)
453c$$$ WRITE(*,*) 'Componentwise Error Bounds'
454c$$$ WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
455c$$$ WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
456c$$$ WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
457c$$$ print *, 'Info: ', info
458c$$$ WRITE(*,*)
459* WRITE(*,*) 'TSTRAT: ',TSTRAT
460
461 END DO
462
463 WRITE(*,*)
464 IF( nfail .GT. 0 ) THEN
465 WRITE(*,9998) c2, nfail, ntests*n+n_aux_tests
466 ELSE
467 WRITE(*,9997) c2
468 END IF
469 9999 FORMAT( ' D', a2, 'SVXX: N =', i2, ', RHS = ', i2,
470 $ ', NWISE GUAR. = ', a, ', CWISE GUAR. = ', a,
471 $ ' test(',i1,') =', g12.5 )
472 9998 FORMAT( ' D', a2, 'SVXX: ', i6, ' out of ', i6,
473 $ ' tests failed to pass the threshold' )
474 9997 FORMAT( ' D', a2, 'SVXX passed the tests of error bounds' )
475* Test ratios.
476 9996 FORMAT( 3x, i2, ': Normwise guaranteed forward error', / 5x,
477 $ 'Guaranteed case: if norm ( abs( Xc - Xt )',
478 $ .LE.' / norm ( Xt ) ERRBND( *, nwise_i, bnd_i ), then',
479 $ / 5x,
480 $ .LE.'ERRBND( *, nwise_i, bnd_i ) MAX(SQRT(N), 10) * EPS')
481 9995 FORMAT( 3x, i2, ': Componentwise guaranteed forward error' )
482 9994 FORMAT( 3x, i2, ': Backwards error' )
483 9993 FORMAT( 3x, i2, ': Reciprocal condition number' )
484 9992 FORMAT( 3x, i2, ': Reciprocal normwise condition number' )
485 9991 FORMAT( 3x, i2, ': Raw normwise error estimate' )
486 9990 FORMAT( 3x, i2, ': Reciprocal componentwise condition number' )
487 9989 FORMAT( 3x, i2, ': Raw componentwise error estimate' )
488
489 8000 FORMAT( ' D', a2, 'SVXX: N =', i2, ', INFO = ', i3,
490 $ ', ORCOND = ', g12.5, ', real RCOND = ', g12.5 )
491*
492* End of DEBCHVXX
493*
494 END
subroutine dlahilb(n, nrhs, a, lda, x, ldx, b, ldb, work, info)
DLAHILB
Definition dlahilb.f:124
subroutine debchvxx(thresh, path)
DEBCHVXX
Definition debchvxx.f:96
subroutine dgbsvxx(fact, trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv, equed, r, c, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DGBSVXX computes the solution to system of linear equations A * X = B for GB matrices
Definition dgbsvxx.f:560
subroutine dgesvxx(fact, trans, n, nrhs, a, lda, af, ldaf, ipiv, equed, r, c, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DGESVXX computes the solution to system of linear equations A * X = B for GE matrices
Definition dgesvxx.f:540
subroutine dsysvxx(fact, uplo, n, nrhs, a, lda, af, ldaf, ipiv, equed, s, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DSYSVXX
Definition dsysvxx.f:505
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
logical function lsamen(n, ca, cb)
LSAMEN
Definition lsamen.f:74
subroutine dposvxx(fact, uplo, n, nrhs, a, lda, af, ldaf, equed, s, b, ldb, x, ldx, rcond, rpvgrw, berr, n_err_bnds, err_bnds_norm, err_bnds_comp, nparams, params, work, iwork, info)
DPOSVXX computes the solution to system of linear equations A * X = B for PO matrices
Definition dposvxx.f:494