/
zppsvx_example.f90
executable file
·96 lines (82 loc) · 3.23 KB
/
zppsvx_example.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
Program zppsvx_example
! ZPPSVX Example Program Text
! Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
! For licence see
! https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md
! .. Use Statements ..
Use lapack_example_aux, Only: nagf_file_print_matrix_complex_gen_comp
Use lapack_interfaces, Only: zppsvx
Use lapack_precision, Only: dp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: nin = 5, nout = 6
Character (1), Parameter :: uplo = 'U'
! .. Local Scalars ..
Real (Kind=dp) :: rcond
Integer :: i, ifail, info, j, ldb, ldx, n, nrhs
Character (1) :: equed
! .. Local Arrays ..
Complex (Kind=dp), Allocatable :: afp(:), ap(:), b(:, :), work(:), &
x(:, :)
Real (Kind=dp), Allocatable :: berr(:), ferr(:), rwork(:), s(:)
Character (1) :: clabs(1), rlabs(1)
! .. Executable Statements ..
Write (nout, *) 'ZPPSVX Example Program Results'
Write (nout, *)
Flush (nout)
! Skip heading in data file
Read (nin, *)
Read (nin, *) n, nrhs
ldb = n
ldx = n
Allocate (afp((n*(n+1))/2), ap((n*(n+1))/2), b(ldb,nrhs), work(2*n), x( &
ldx,nrhs), berr(nrhs), ferr(nrhs), rwork(n), s(n))
! Read the upper or lower triangular part of the matrix A from
! data file
If (uplo=='U') Then
Read (nin, *)((ap(i+(j*(j-1))/2),j=i,n), i=1, n)
Else If (uplo=='L') Then
Read (nin, *)((ap(i+((2*n-j)*(j-1))/2),j=1,i), i=1, n)
End If
! Read B from data file
Read (nin, *)(b(i,1:nrhs), i=1, n)
! Solve the equations AX = B for X
Call zppsvx('Equilibration', uplo, n, nrhs, ap, afp, equed, s, b, ldb, &
x, ldx, rcond, ferr, berr, work, rwork, info)
If ((info==0) .Or. (info==n+1)) Then
! Print solution, error bounds, condition number and the form
! of equilibration
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
Call nagf_file_print_matrix_complex_gen_comp('General', ' ', n, nrhs, &
x, ldx, 'Bracketed', 'F7.4', 'Solution(s)', 'Integer', rlabs, &
'Integer', clabs, 80, 0, ifail)
Write (nout, *)
Write (nout, *) 'Backward errors (machine-dependent)'
Write (nout, 100) berr(1:nrhs)
Write (nout, *)
Write (nout, *) 'Estimated forward error bounds (machine-dependent)'
Write (nout, 100) ferr(1:nrhs)
Write (nout, *)
Write (nout, *) 'Estimate of reciprocal condition number'
Write (nout, 100) rcond
Write (nout, *)
If (equed=='N') Then
Write (nout, *) 'A has not been equilibrated'
Else If (equed=='Y') Then
Write (nout, *) &
'A has been row and column scaled as diag(S)*A*diag(S)'
End If
If (info==n+1) Then
Write (nout, *)
Write (nout, *) 'The matrix A is singular to working precision'
End If
Else
Write (nout, 110) 'The leading minor of order ', info, &
' is not positive definite'
End If
100 Format ((3X,1P,7E11.1))
110 Format (1X, A, I3, A)
End Program