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