MAGMA  1.2.0 MatrixAlgebraonGPUandMulticoreArchitectures
zget01.f
Go to the documentation of this file.
1  SUBROUTINE zget01( M, N, A, LDA, AFAC, LDAFAC, IPIV, RWORK,
2  \$ resid )
3 *
4 * -- LAPACK test routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9  INTEGER lda, ldafac, m, n
10  DOUBLE PRECISION resid
11 * ..
12 * .. Array Arguments ..
13  INTEGER ipiv( * )
14  DOUBLE PRECISION rwork( * )
15  COMPLEX*16 a( lda, * ), afac( ldafac, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * ZGET01 reconstructs a matrix A from its L*U factorization and
22 * computes the residual
23 * norm(L*U - A) / ( N * norm(A) * EPS ),
24 * where EPS is the machine epsilon.
25 *
26 * Arguments
27 * ==========
28 *
29 * M (input) INTEGER
30 * The number of rows of the matrix A. M >= 0.
31 *
32 * N (input) INTEGER
33 * The number of columns of the matrix A. N >= 0.
34 *
35 * A (input) COMPLEX*16 array, dimension (LDA,N)
36 * The original M x N matrix A.
37 *
38 * LDA (input) INTEGER
39 * The leading dimension of the array A. LDA >= max(1,M).
40 *
41 * AFAC (input/output) COMPLEX*16 array, dimension (LDAFAC,N)
42 * The factored form of the matrix A. AFAC contains the factors
43 * L and U from the L*U factorization as computed by ZGETRF.
44 * Overwritten with the reconstructed matrix, and then with the
45 * difference L*U - A.
46 *
47 * LDAFAC (input) INTEGER
48 * The leading dimension of the array AFAC. LDAFAC >= max(1,M).
49 *
50 * IPIV (input) INTEGER array, dimension (N)
51 * The pivot indices from ZGETRF.
52 *
53 * RWORK (workspace) DOUBLE PRECISION array, dimension (M)
54 *
55 * RESID (output) DOUBLE PRECISION
56 * norm(L*U - A) / ( N * norm(A) * EPS )
57 *
58 * =====================================================================
59 *
60 * .. Parameters ..
61  DOUBLE PRECISION zero, one
62  parameter( zero = 0.0d+0, one = 1.0d+0 )
63  COMPLEX*16 cone
64  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
65 * ..
66 * .. Local Scalars ..
67  INTEGER i, j, k
68  DOUBLE PRECISION anorm, eps
69  COMPLEX*16 t
70 * ..
71 * .. External Functions ..
72  DOUBLE PRECISION dlamch, zlange
73  COMPLEX*16 zdotu
74  EXTERNAL dlamch, zlange, zdotu
75 * ..
76 * .. External Subroutines ..
77  EXTERNAL zgemv, zlaswp, zscal, ztrmv
78 * ..
79 * .. Intrinsic Functions ..
80  INTRINSIC dble, min
81 * ..
82 * .. Executable Statements ..
83 *
84 * Quick exit if M = 0 or N = 0.
85 *
86  IF( m.LE.0 .OR. n.LE.0 ) THEN
87  resid = zero
88  return
89  END IF
90 *
91 * Determine EPS and the norm of A.
92 *
93  eps = dlamch( 'Epsilon' )
94  anorm = zlange( '1', m, n, a, lda, rwork )
95 *
96 * Compute the product L*U and overwrite AFAC with the result.
97 * A column at a time of the product is obtained, starting with
98 * column N.
99 *
100  DO 10 k = n, 1, -1
101  IF( k.GT.m ) THEN
102  CALL ztrmv( 'Lower', 'No transpose', 'Unit', m, afac,
103  \$ ldafac, afac( 1, k ), 1 )
104  ELSE
105 *
106 * Compute elements (K+1:M,K)
107 *
108  t = afac( k, k )
109  IF( k+1.LE.m ) THEN
110  CALL zscal( m-k, t, afac( k+1, k ), 1 )
111  CALL zgemv( 'No transpose', m-k, k-1, cone,
112  \$ afac( k+1, 1 ), ldafac, afac( 1, k ), 1,
113  \$ cone, afac( k+1, k ), 1 )
114  END IF
115 *
116 * Compute the (K,K) element
117 *
118  afac( k, k ) = t + zdotu( k-1, afac( k, 1 ), ldafac,
119  \$ afac( 1, k ), 1 )
120 *
121 * Compute elements (1:K-1,K)
122 *
123  CALL ztrmv( 'Lower', 'No transpose', 'Unit', k-1, afac,
124  \$ ldafac, afac( 1, k ), 1 )
125  END IF
126  10 continue
127  CALL zlaswp( n, afac, ldafac, 1, min( m, n ), ipiv, -1 )
128 *
129 * Compute the difference L*U - A and store in AFAC.
130 *
131  DO 30 j = 1, n
132  DO 20 i = 1, m
133  afac( i, j ) = afac( i, j ) - a( i, j )
134  20 continue
135  30 continue
136 *
137 * Compute norm( L*U - A ) / ( N * norm(A) * EPS )
138 *
139  resid = zlange( '1', m, n, afac, ldafac, rwork )
140 *
141  IF( anorm.LE.zero ) THEN
142  IF( resid.NE.zero )
143  \$ resid = one / eps
144  ELSE
145  resid = ( ( resid / dble( n ) ) / anorm ) / eps
146  END IF
147 *
148  return
149 *
150 * End of ZGET01
151 *
152  END