LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
zpbsvx.f
Go to the documentation of this file.
1 *> \brief <b> ZPBSVX computes the solution to system of linear equations A * X = B for OTHER matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpbsvx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpbsvx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpbsvx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZPBSVX( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB,
22 * EQUED, S, B, LDB, X, LDX, RCOND, FERR, BERR,
23 * WORK, RWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER EQUED, FACT, UPLO
27 * INTEGER INFO, KD, LDAB, LDAFB, LDB, LDX, N, NRHS
28 * DOUBLE PRECISION RCOND
29 * ..
30 * .. Array Arguments ..
31 * DOUBLE PRECISION BERR( * ), FERR( * ), RWORK( * ), S( * )
32 * COMPLEX*16 AB( LDAB, * ), AFB( LDAFB, * ), B( LDB, * ),
33 * \$ WORK( * ), X( LDX, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
43 *> compute the solution to a complex system of linear equations
44 *> A * X = B,
45 *> where A is an N-by-N Hermitian positive definite band matrix and X
46 *> and B are N-by-NRHS matrices.
47 *>
48 *> Error bounds on the solution and a condition estimate are also
49 *> provided.
50 *> \endverbatim
51 *
52 *> \par Description:
53 * =================
54 *>
55 *> \verbatim
56 *>
57 *> The following steps are performed:
58 *>
59 *> 1. If FACT = 'E', real scaling factors are computed to equilibrate
60 *> the system:
61 *> diag(S) * A * diag(S) * inv(diag(S)) * X = diag(S) * B
62 *> Whether or not the system will be equilibrated depends on the
63 *> scaling of the matrix A, but if equilibration is used, A is
64 *> overwritten by diag(S)*A*diag(S) and B by diag(S)*B.
65 *>
66 *> 2. If FACT = 'N' or 'E', the Cholesky decomposition is used to
67 *> factor the matrix A (after equilibration if FACT = 'E') as
68 *> A = U**H * U, if UPLO = 'U', or
69 *> A = L * L**H, if UPLO = 'L',
70 *> where U is an upper triangular band matrix, and L is a lower
71 *> triangular band matrix.
72 *>
73 *> 3. If the leading i-by-i principal minor is not positive definite,
74 *> then the routine returns with INFO = i. Otherwise, the factored
75 *> form of A is used to estimate the condition number of the matrix
76 *> A. If the reciprocal of the condition number is less than machine
77 *> precision, INFO = N+1 is returned as a warning, but the routine
78 *> still goes on to solve for X and compute error bounds as
79 *> described below.
80 *>
81 *> 4. The system of equations is solved for X using the factored form
82 *> of A.
83 *>
84 *> 5. Iterative refinement is applied to improve the computed solution
85 *> matrix and calculate error bounds and backward error estimates
86 *> for it.
87 *>
88 *> 6. If equilibration was used, the matrix X is premultiplied by
89 *> diag(S) so that it solves the original system before
90 *> equilibration.
91 *> \endverbatim
92 *
93 * Arguments:
94 * ==========
95 *
96 *> \param[in] FACT
97 *> \verbatim
98 *> FACT is CHARACTER*1
99 *> Specifies whether or not the factored form of the matrix A is
100 *> supplied on entry, and if not, whether the matrix A should be
101 *> equilibrated before it is factored.
102 *> = 'F': On entry, AFB contains the factored form of A.
103 *> If EQUED = 'Y', the matrix A has been equilibrated
104 *> with scaling factors given by S. AB and AFB will not
105 *> be modified.
106 *> = 'N': The matrix A will be copied to AFB and factored.
107 *> = 'E': The matrix A will be equilibrated if necessary, then
108 *> copied to AFB and factored.
109 *> \endverbatim
110 *>
111 *> \param[in] UPLO
112 *> \verbatim
113 *> UPLO is CHARACTER*1
114 *> = 'U': Upper triangle of A is stored;
115 *> = 'L': Lower triangle of A is stored.
116 *> \endverbatim
117 *>
118 *> \param[in] N
119 *> \verbatim
120 *> N is INTEGER
121 *> The number of linear equations, i.e., the order of the
122 *> matrix A. N >= 0.
123 *> \endverbatim
124 *>
125 *> \param[in] KD
126 *> \verbatim
127 *> KD is INTEGER
128 *> The number of superdiagonals of the matrix A if UPLO = 'U',
129 *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
130 *> \endverbatim
131 *>
132 *> \param[in] NRHS
133 *> \verbatim
134 *> NRHS is INTEGER
135 *> The number of right-hand sides, i.e., the number of columns
136 *> of the matrices B and X. NRHS >= 0.
137 *> \endverbatim
138 *>
139 *> \param[in,out] AB
140 *> \verbatim
141 *> AB is COMPLEX*16 array, dimension (LDAB,N)
142 *> On entry, the upper or lower triangle of the Hermitian band
143 *> matrix A, stored in the first KD+1 rows of the array, except
144 *> if FACT = 'F' and EQUED = 'Y', then A must contain the
145 *> equilibrated matrix diag(S)*A*diag(S). The j-th column of A
146 *> is stored in the j-th column of the array AB as follows:
147 *> if UPLO = 'U', AB(KD+1+i-j,j) = A(i,j) for max(1,j-KD)<=i<=j;
148 *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(N,j+KD).
149 *> See below for further details.
150 *>
151 *> On exit, if FACT = 'E' and EQUED = 'Y', A is overwritten by
152 *> diag(S)*A*diag(S).
153 *> \endverbatim
154 *>
155 *> \param[in] LDAB
156 *> \verbatim
157 *> LDAB is INTEGER
158 *> The leading dimension of the array A. LDAB >= KD+1.
159 *> \endverbatim
160 *>
161 *> \param[in,out] AFB
162 *> \verbatim
163 *> AFB is COMPLEX*16 array, dimension (LDAFB,N)
164 *> If FACT = 'F', then AFB is an input argument and on entry
165 *> contains the triangular factor U or L from the Cholesky
166 *> factorization A = U**H *U or A = L*L**H of the band matrix
167 *> A, in the same storage format as A (see AB). If EQUED = 'Y',
168 *> then AFB is the factored form of the equilibrated matrix A.
169 *>
170 *> If FACT = 'N', then AFB is an output argument and on exit
171 *> returns the triangular factor U or L from the Cholesky
172 *> factorization A = U**H *U or A = L*L**H.
173 *>
174 *> If FACT = 'E', then AFB is an output argument and on exit
175 *> returns the triangular factor U or L from the Cholesky
176 *> factorization A = U**H *U or A = L*L**H of the equilibrated
177 *> matrix A (see the description of A for the form of the
178 *> equilibrated matrix).
179 *> \endverbatim
180 *>
181 *> \param[in] LDAFB
182 *> \verbatim
183 *> LDAFB is INTEGER
184 *> The leading dimension of the array AFB. LDAFB >= KD+1.
185 *> \endverbatim
186 *>
187 *> \param[in,out] EQUED
188 *> \verbatim
189 *> EQUED is CHARACTER*1
190 *> Specifies the form of equilibration that was done.
191 *> = 'N': No equilibration (always true if FACT = 'N').
192 *> = 'Y': Equilibration was done, i.e., A has been replaced by
193 *> diag(S) * A * diag(S).
194 *> EQUED is an input argument if FACT = 'F'; otherwise, it is an
195 *> output argument.
196 *> \endverbatim
197 *>
198 *> \param[in,out] S
199 *> \verbatim
200 *> S is DOUBLE PRECISION array, dimension (N)
201 *> The scale factors for A; not accessed if EQUED = 'N'. S is
202 *> an input argument if FACT = 'F'; otherwise, S is an output
203 *> argument. If FACT = 'F' and EQUED = 'Y', each element of S
204 *> must be positive.
205 *> \endverbatim
206 *>
207 *> \param[in,out] B
208 *> \verbatim
209 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
210 *> On entry, the N-by-NRHS right hand side matrix B.
211 *> On exit, if EQUED = 'N', B is not modified; if EQUED = 'Y',
212 *> B is overwritten by diag(S) * B.
213 *> \endverbatim
214 *>
215 *> \param[in] LDB
216 *> \verbatim
217 *> LDB is INTEGER
218 *> The leading dimension of the array B. LDB >= max(1,N).
219 *> \endverbatim
220 *>
221 *> \param[out] X
222 *> \verbatim
223 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
224 *> If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X to
225 *> the original system of equations. Note that if EQUED = 'Y',
226 *> A and B are modified on exit, and the solution to the
227 *> equilibrated system is inv(diag(S))*X.
228 *> \endverbatim
229 *>
230 *> \param[in] LDX
231 *> \verbatim
232 *> LDX is INTEGER
233 *> The leading dimension of the array X. LDX >= max(1,N).
234 *> \endverbatim
235 *>
236 *> \param[out] RCOND
237 *> \verbatim
238 *> RCOND is DOUBLE PRECISION
239 *> The estimate of the reciprocal condition number of the matrix
240 *> A after equilibration (if done). If RCOND is less than the
241 *> machine precision (in particular, if RCOND = 0), the matrix
242 *> is singular to working precision. This condition is
243 *> indicated by a return code of INFO > 0.
244 *> \endverbatim
245 *>
246 *> \param[out] FERR
247 *> \verbatim
248 *> FERR is DOUBLE PRECISION array, dimension (NRHS)
249 *> The estimated forward error bound for each solution vector
250 *> X(j) (the j-th column of the solution matrix X).
251 *> If XTRUE is the true solution corresponding to X(j), FERR(j)
252 *> is an estimated upper bound for the magnitude of the largest
253 *> element in (X(j) - XTRUE) divided by the magnitude of the
254 *> largest element in X(j). The estimate is as reliable as
255 *> the estimate for RCOND, and is almost always a slight
256 *> overestimate of the true error.
257 *> \endverbatim
258 *>
259 *> \param[out] BERR
260 *> \verbatim
261 *> BERR is DOUBLE PRECISION array, dimension (NRHS)
262 *> The componentwise relative backward error of each solution
263 *> vector X(j) (i.e., the smallest relative change in
264 *> any element of A or B that makes X(j) an exact solution).
265 *> \endverbatim
266 *>
267 *> \param[out] WORK
268 *> \verbatim
269 *> WORK is COMPLEX*16 array, dimension (2*N)
270 *> \endverbatim
271 *>
272 *> \param[out] RWORK
273 *> \verbatim
274 *> RWORK is DOUBLE PRECISION array, dimension (N)
275 *> \endverbatim
276 *>
277 *> \param[out] INFO
278 *> \verbatim
279 *> INFO is INTEGER
280 *> = 0: successful exit
281 *> < 0: if INFO = -i, the i-th argument had an illegal value
282 *> > 0: if INFO = i, and i is
283 *> <= N: the leading minor of order i of A is
284 *> not positive definite, so the factorization
285 *> could not be completed, and the solution has not
286 *> been computed. RCOND = 0 is returned.
287 *> = N+1: U is nonsingular, but RCOND is less than machine
288 *> precision, meaning that the matrix is singular
289 *> to working precision. Nevertheless, the
290 *> solution and error bounds are computed because
291 *> there are a number of situations where the
292 *> computed solution can be more accurate than the
293 *> value of RCOND would suggest.
294 *> \endverbatim
295 *
296 * Authors:
297 * ========
298 *
299 *> \author Univ. of Tennessee
300 *> \author Univ. of California Berkeley
301 *> \author Univ. of Colorado Denver
302 *> \author NAG Ltd.
303 *
304 *> \date April 2012
305 *
306 *> \ingroup complex16OTHERsolve
307 *
308 *> \par Further Details:
309 * =====================
310 *>
311 *> \verbatim
312 *>
313 *> The band storage scheme is illustrated by the following example, when
314 *> N = 6, KD = 2, and UPLO = 'U':
315 *>
316 *> Two-dimensional storage of the Hermitian matrix A:
317 *>
318 *> a11 a12 a13
319 *> a22 a23 a24
320 *> a33 a34 a35
321 *> a44 a45 a46
322 *> a55 a56
323 *> (aij=conjg(aji)) a66
324 *>
325 *> Band storage of the upper triangle of A:
326 *>
327 *> * * a13 a24 a35 a46
328 *> * a12 a23 a34 a45 a56
329 *> a11 a22 a33 a44 a55 a66
330 *>
331 *> Similarly, if UPLO = 'L' the format of A is as follows:
332 *>
333 *> a11 a22 a33 a44 a55 a66
334 *> a21 a32 a43 a54 a65 *
335 *> a31 a42 a53 a64 * *
336 *>
337 *> Array elements marked * are not used by the routine.
338 *> \endverbatim
339 *>
340 * =====================================================================
341  SUBROUTINE zpbsvx( FACT, UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB,
342  \$ equed, s, b, ldb, x, ldx, rcond, ferr, berr,
343  \$ work, rwork, info )
344 *
345 * -- LAPACK driver routine (version 3.4.1) --
346 * -- LAPACK is a software package provided by Univ. of Tennessee, --
347 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
348 * April 2012
349 *
350 * .. Scalar Arguments ..
351  CHARACTER equed, fact, uplo
352  INTEGER info, kd, ldab, ldafb, ldb, ldx, n, nrhs
353  DOUBLE PRECISION rcond
354 * ..
355 * .. Array Arguments ..
356  DOUBLE PRECISION berr( * ), ferr( * ), rwork( * ), s( * )
357  COMPLEX*16 ab( ldab, * ), afb( ldafb, * ), b( ldb, * ),
358  \$ work( * ), x( ldx, * )
359 * ..
360 *
361 * =====================================================================
362 *
363 * .. Parameters ..
364  DOUBLE PRECISION zero, one
365  parameter( zero = 0.0d+0, one = 1.0d+0 )
366 * ..
367 * .. Local Scalars ..
368  LOGICAL equil, nofact, rcequ, upper
369  INTEGER i, infequ, j, j1, j2
370  DOUBLE PRECISION amax, anorm, bignum, scond, smax, smin, smlnum
371 * ..
372 * .. External Functions ..
373  LOGICAL lsame
374  DOUBLE PRECISION dlamch, zlanhb
375  EXTERNAL lsame, dlamch, zlanhb
376 * ..
377 * .. External Subroutines ..
378  EXTERNAL xerbla, zcopy, zlacpy, zlaqhb, zpbcon, zpbequ,
379  \$ zpbrfs, zpbtrf, zpbtrs
380 * ..
381 * .. Intrinsic Functions ..
382  INTRINSIC max, min
383 * ..
384 * .. Executable Statements ..
385 *
386  info = 0
387  nofact = lsame( fact, 'N' )
388  equil = lsame( fact, 'E' )
389  upper = lsame( uplo, 'U' )
390  IF( nofact .OR. equil ) THEN
391  equed = 'N'
392  rcequ = .false.
393  ELSE
394  rcequ = lsame( equed, 'Y' )
395  smlnum = dlamch( 'Safe minimum' )
396  bignum = one / smlnum
397  END IF
398 *
399 * Test the input parameters.
400 *
401  IF( .NOT.nofact .AND. .NOT.equil .AND. .NOT.lsame( fact, 'F' ) )
402  \$ THEN
403  info = -1
404  ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
405  info = -2
406  ELSE IF( n.LT.0 ) THEN
407  info = -3
408  ELSE IF( kd.LT.0 ) THEN
409  info = -4
410  ELSE IF( nrhs.LT.0 ) THEN
411  info = -5
412  ELSE IF( ldab.LT.kd+1 ) THEN
413  info = -7
414  ELSE IF( ldafb.LT.kd+1 ) THEN
415  info = -9
416  ELSE IF( lsame( fact, 'F' ) .AND. .NOT.
417  \$ ( rcequ .OR. lsame( equed, 'N' ) ) ) THEN
418  info = -10
419  ELSE
420  IF( rcequ ) THEN
421  smin = bignum
422  smax = zero
423  DO 10 j = 1, n
424  smin = min( smin, s( j ) )
425  smax = max( smax, s( j ) )
426  10 continue
427  IF( smin.LE.zero ) THEN
428  info = -11
429  ELSE IF( n.GT.0 ) THEN
430  scond = max( smin, smlnum ) / min( smax, bignum )
431  ELSE
432  scond = one
433  END IF
434  END IF
435  IF( info.EQ.0 ) THEN
436  IF( ldb.LT.max( 1, n ) ) THEN
437  info = -13
438  ELSE IF( ldx.LT.max( 1, n ) ) THEN
439  info = -15
440  END IF
441  END IF
442  END IF
443 *
444  IF( info.NE.0 ) THEN
445  CALL xerbla( 'ZPBSVX', -info )
446  return
447  END IF
448 *
449  IF( equil ) THEN
450 *
451 * Compute row and column scalings to equilibrate the matrix A.
452 *
453  CALL zpbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
454  IF( infequ.EQ.0 ) THEN
455 *
456 * Equilibrate the matrix.
457 *
458  CALL zlaqhb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
459  rcequ = lsame( equed, 'Y' )
460  END IF
461  END IF
462 *
463 * Scale the right-hand side.
464 *
465  IF( rcequ ) THEN
466  DO 30 j = 1, nrhs
467  DO 20 i = 1, n
468  b( i, j ) = s( i )*b( i, j )
469  20 continue
470  30 continue
471  END IF
472 *
473  IF( nofact .OR. equil ) THEN
474 *
475 * Compute the Cholesky factorization A = U**H *U or A = L*L**H.
476 *
477  IF( upper ) THEN
478  DO 40 j = 1, n
479  j1 = max( j-kd, 1 )
480  CALL zcopy( j-j1+1, ab( kd+1-j+j1, j ), 1,
481  \$ afb( kd+1-j+j1, j ), 1 )
482  40 continue
483  ELSE
484  DO 50 j = 1, n
485  j2 = min( j+kd, n )
486  CALL zcopy( j2-j+1, ab( 1, j ), 1, afb( 1, j ), 1 )
487  50 continue
488  END IF
489 *
490  CALL zpbtrf( uplo, n, kd, afb, ldafb, info )
491 *
492 * Return if INFO is non-zero.
493 *
494  IF( info.GT.0 )THEN
495  rcond = zero
496  return
497  END IF
498  END IF
499 *
500 * Compute the norm of the matrix A.
501 *
502  anorm = zlanhb( '1', uplo, n, kd, ab, ldab, rwork )
503 *
504 * Compute the reciprocal of the condition number of A.
505 *
506  CALL zpbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, rwork,
507  \$ info )
508 *
509 * Compute the solution matrix X.
510 *
511  CALL zlacpy( 'Full', n, nrhs, b, ldb, x, ldx )
512  CALL zpbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
513 *
514 * Use iterative refinement to improve the computed solution and
515 * compute error bounds and backward error estimates for it.
516 *
517  CALL zpbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,
518  \$ ldx, ferr, berr, work, rwork, info )
519 *
520 * Transform the solution matrix X to a solution of the original
521 * system.
522 *
523  IF( rcequ ) THEN
524  DO 70 j = 1, nrhs
525  DO 60 i = 1, n
526  x( i, j ) = s( i )*x( i, j )
527  60 continue
528  70 continue
529  DO 80 j = 1, nrhs
530  ferr( j ) = ferr( j ) / scond
531  80 continue
532  END IF
533 *
534 * Set INFO = N+1 if the matrix is singular to working precision.
535 *
536  IF( rcond.LT.dlamch( 'Epsilon' ) )
537  \$ info = n + 1
538 *
539  return
540 *
541 * End of ZPBSVX
542 *
543  END