Developer Reference

Contents

p?posvx

Solves a symmetric or Hermitian positive definite system of linear equations.

Syntax

void
psposvx
(
char
*fact
,
char
*uplo
,
MKL_INT
*n
,
MKL_INT
*nrhs
,
float
*a
,
MKL_INT
*ia
,
MKL_INT
*ja
,
MKL_INT
*desca
,
float
*af
,
MKL_INT
*iaf
,
MKL_INT
*jaf
,
MKL_INT
*descaf
,
char
*equed
,
float
*sr
,
float
*sc
,
float
*b
,
MKL_INT
*ib
,
MKL_INT
*jb
,
MKL_INT
*descb
,
float
*x
,
MKL_INT
*ix
,
MKL_INT
*jx
,
MKL_INT
*descx
,
float
*rcond
,
float
*ferr
,
float
*berr
,
float
*work
,
MKL_INT
*lwork
,
MKL_INT
*iwork
,
MKL_INT
*liwork
,
MKL_INT
*info
);
void
pdposvx
(
char
*fact
,
char
*uplo
,
MKL_INT
*n
,
MKL_INT
*nrhs
,
double
*a
,
MKL_INT
*ia
,
MKL_INT
*ja
,
MKL_INT
*desca
,
double
*af
,
MKL_INT
*iaf
,
MKL_INT
*jaf
,
MKL_INT
*descaf
,
char
*equed
,
double
*sr
,
double
*sc
,
double
*b
,
MKL_INT
*ib
,
MKL_INT
*jb
,
MKL_INT
*descb
,
double
*x
,
MKL_INT
*ix
,
MKL_INT
*jx
,
MKL_INT
*descx
,
double
*rcond
,
double
*ferr
,
double
*berr
,
double
*work
,
MKL_INT
*lwork
,
MKL_INT
*iwork
,
MKL_INT
*liwork
,
MKL_INT
*info
);
void
pcposvx
(
char
*fact
,
char
*uplo
,
MKL_INT
*n
,
MKL_INT
*nrhs
,
MKL_Complex8
*a
,
MKL_INT
*ia
,
MKL_INT
*ja
,
MKL_INT
*desca
,
MKL_Complex8
*af
,
MKL_INT
*iaf
,
MKL_INT
*jaf
,
MKL_INT
*descaf
,
char
*equed
,
float
*sr
,
float
*sc
,
MKL_Complex8
*b
,
MKL_INT
*ib
,
MKL_INT
*jb
,
MKL_INT
*descb
,
MKL_Complex8
*x
,
MKL_INT
*ix
,
MKL_INT
*jx
,
MKL_INT
*descx
,
float
*rcond
,
float
*ferr
,
float
*berr
,
MKL_Complex8
*work
,
MKL_INT
*lwork
,
float
*rwork
,
MKL_INT
*lrwork
,
MKL_INT
*info
);
void
pzposvx
(
char
*fact
,
char
*uplo
,
MKL_INT
*n
,
MKL_INT
*nrhs
,
MKL_Complex16
*a
,
MKL_INT
*ia
,
MKL_INT
*ja
,
MKL_INT
*desca
,
MKL_Complex16
*af
,
MKL_INT
*iaf
,
MKL_INT
*jaf
,
MKL_INT
*descaf
,
char
*equed
,
double
*sr
,
double
*sc
,
MKL_Complex16
*b
,
MKL_INT
*ib
,
MKL_INT
*jb
,
MKL_INT
*descb
,
MKL_Complex16
*x
,
MKL_INT
*ix
,
MKL_INT
*jx
,
MKL_INT
*descx
,
double
*rcond
,
double
*ferr
,
double
*berr
,
MKL_Complex16
*work
,
MKL_INT
*lwork
,
double
*rwork
,
MKL_INT
*lrwork
,
MKL_INT
*info
);
Include Files
  • mkl_scalapack.h
Description
The
p?posvx
function
uses the Cholesky factorization
A
=
U
T
*U
or
A
=
L*L
T
to compute the solution to a real or complex system of linear equations
A
(
ia
:
ia
+
n
-1,
ja
:
ja
+
n
-1)*
X
=
B
(
ib
:
ib
+
n
-1,
jb
:
jb
+
nrhs
-1)
,
where
A
(
ia
:
ia
+
n
-1,
ja
:
ja
+
n
-1)
is a
n
-by-
n
matrix and
X
and
B
(
ib
:
ib
+
n
-1,
jb
:
jb
+
nrhs
-1)
are
n
-by-
nrhs
matrices.
Error bounds on the solution and a condition estimate are also provided.
In the following comments
y
denotes
Y
(
iy
:
iy
+
m
-1,
jy
:
jy
+
k
-1)
, an
m
-by-
k
matrix where
y
can be
a
,
af
,
b
and
x
.
The
function
p?posvx
performs the following steps:
  1. If
    fact
    =
    'E'
    , real scaling factors
    s
    are computed to equilibrate the system:
    diag(
    sr
    )*
    A
    *diag(
    sc
    )*inv(diag(
    sc
    ))*
    X
    = diag(
    sr
    )*B
    Whether or not the system will be equilibrated depends on the scaling of the matrix
    A
    , but if equilibration is used,
    A
    is overwritten by
    diag(
    sr
    )*
    A
    *diag(
    sc
    )
    and
    B
    by
    diag(
    sr
    )*
    B
    .
  2. If
    fact
    =
    'N'
    or
    'E'
    , the Cholesky decomposition is used to factor the matrix
    A
    (after equilibration if
    fact
    =
    'E'
    ) as
    A
    =
    U
    T
    *U
    , if
    uplo
    =
    'U'
    , or
    A
    =
    L*L
    T
    , if
    uplo
    =
    'L'
    ,
    where
    U
    is an upper triangular matrix and
    L
    is a lower triangular matrix.
  3. The factored form of
    A
    is used to estimate the condition number of the matrix
    A
    . If the reciprocal of the condition number is less than machine precision, steps 4-6 are skipped
  4. The system of equations is solved for
    X
    using the factored form of
    A
    .
  5. Iterative refinement is applied to improve the computed solution matrix and calculate error bounds and backward error estimates for it.
  6. If equilibration was used, the matrix
    X
    is premultiplied by diag(
    sr
    ) so that it solves the original system before equilibration.
Input Parameters
fact
(global) Must be
'F'
,
'N'
, or
'E'
.
Specifies whether or not the factored form of the matrix
A
is supplied on entry, and if not, whether the matrix
A
should be equilibrated before it is factored.
If
fact
=
'F'
: on entry,
af
contains the factored form of
A
. If
equed
=
'Y'
, the matrix
A
has been equilibrated with scaling factors given by
s
.
a
and
af
will not be modified.
If
fact
=
'N'
, the matrix
A
will be copied to
af
and factored.
If
fact
=
'E'
, the matrix
A
will be equilibrated if necessary, then copied to
af
and factored.
uplo
(global) Must be
'U'
or
'L'
.
Indicates whether the upper or lower triangular part of
A
is stored.
n
(global) The order of the distributed matrix sub(
A
)
(
n
0)
.
nrhs
(global) The number of right-hand sides; the number of columns of the distributed submatrices
B
and
X
. (
nrhs
0)
.
a
(local)
Pointer into the local memory to an array of local size
lld_a
*
LOCc
(
ja
+
n
-1)
. On entry, the symmetric/Hermitian matrix
A
, except if
fact
=
'F'
and
equed
=
'Y'
, then
A
must contain the equilibrated matrix
diag(
sr
)*
A
*diag(
sc
)
.
If
uplo
=
'U'
, the leading
n
-by-
n
upper triangular part of
A
contains the upper triangular part of the matrix
A
, and the strictly lower triangular part of
A
is not referenced.
If
uplo
=
'L'
, the leading
n
-by-
n
lower triangular part of
A
contains the lower triangular part of the matrix
A
, and the strictly upper triangular part of
A
is not referenced.
A
is not modified if
fact
=
'F'
or
'N'
, or if
fact
=
'E'
and
equed
=
'N'
on exit.
ia
,
ja
(global) The row and column indices in the global matrix
A
indicating the first row and the first column of the submatrix
A
, respectively.
desca
(global and local) array of size
dlen_
. The array descriptor for the distributed matrix
A
.
af
(local)
Pointer into the local memory to an array of local size
lld_af
*
LOCc
(
ja
+
n
-1)
.
If
fact
=
'F'
, then
af
is an input argument and on entry contains the triangular factor
U
or
L
from the Cholesky factorization
A
=
U
T
*
U
or
A
=
L
*
L
T
, in the same storage format as
A
. If
equed
'N'
, then
af
is the factored form of the equilibrated matrix
diag(
sr
)*
A
*diag(
sc
)
.
iaf
,
jaf
(global) The row and column indices in the global matrix
AF
indicating the first row and the first column of the submatrix
AF
, respectively.
descaf
(global and local) array of size
dlen_
. The array descriptor for the distributed matrix
AF
.
equed
(global) Must be
'N'
or
'Y'
.
equed
is an input argument if
fact
=
'F'
. It specifies the form of equilibration that was done:
If
equed
=
'N'
, no equilibration was done (always true if
fact
=
'N'
);
If
equed
=
'Y'
, equilibration was done and
A
has been replaced by
diag(
sr
)*
A
*diag(
sc
)
.
sr
(local)
Array of size
lld_a
.
The array
s
contains the scale factors for
A
. This array is an input argument if
fact
=
'F'
only; otherwise it is an output argument.
If
equed
=
'N'
,
s
is not accessed.
If
fact
=
'F'
and
equed
=
'Y'
, each element of
s
must be positive.
b
(local)
Pointer into the local memory to an array of local size
lld_b
*
LOCc
(
jb
+
nrhs
-1)
. On entry, the
n
-by-
nrhs
right-hand side matrix
B
.
ib
,
jb
(global) The row and column indices in the global matrix
B
indicating the first row and the first column of the submatrix
B
, respectively.
descb
(global and local) Array of size
dlen_
. The array descriptor for the distributed matrix
B
.
x
(local)
Pointer into the local memory to an array of local size
lld_x
*
LOCc
(
jx
+
nrhs
-1)
.
ix
,
jx
(global) The row and column indices in the global matrix
X
indicating the first row and the first column of the submatrix
X
, respectively.
descx
(global and local) array of size
dlen_
. The array descriptor for the distributed matrix
X
.
work
(local)
Workspace array of size
lwork
.
lwork
(local or global)
The size of the array
work
.
lwork
is local input and must be at least
lwork
=
max
(
p?pocon
(
lwork
),
p?porfs
(
lwork
)) +
LOCr
(
n_a
)
.
lwork
= 3*
desca
[
lld_
- 1]
.
If
lwork
= -1
, then
lwork
is global input and a workspace query is assumed; the
function
only calculates the minimum and optimal size for all work arrays. Each of these values is returned in the first entry of the corresponding work array, and no error message is issued by
pxerbla
.
iwork
(local) Workspace array of size
liwork
.
liwork
(local or global)
The size of the array
iwork
.
liwork
is local input and must be at least
liwork
=
desca
[
lld_
- 1]
liwork
=
LOCr
(
n_a
)
.
If
liwork
= -1, then
liwork
is global input and a workspace query is assumed; the
function
only calculates the minimum and optimal size for all work arrays. Each of these values is returned in the first entry of the corresponding work array, and no error message is issued by
pxerbla
.
Output Parameters
a
On exit, if
fact
=
'E'
and
equed
=
'Y'
,
a
is overwritten by
diag(
sr
)*
a
*diag(
sc
)
.
af
If
fact
=
'N'
, then
af
is an output argument and on exit returns the triangular factor
U
or
L
from the Cholesky factorization
A
=
U
T
*
U
or
A
=
L*L
T
of the original matrix
A
.
If
fact
=
'E'
, then
af
is an output argument and on exit returns the triangular factor
U
or
L
from the Cholesky factorization
A
=
U
T
*
U
or
A
=
L*L
T
of the equilibrated matrix
A
(see the description of
A
for the form of the equilibrated matrix).
equed
If
fact
'F'
, then
equed
is an output argument. It specifies the form of equilibration that was done (see the description of
equed
in
Input Arguments
section).
sr
This array is an output argument if
fact
'F'
.
See the description of
sr
in
Input Arguments
section.
sc
This array is an output argument if
fact
'F'
.
See the description of
sc
in
Input Arguments
section.
b
On exit, if
equed
=
'N'
,
b
is not modified; if
trans
=
'N'
and
equed
=
'R'
or
'B'
,
b
is overwritten by
diag(
r
)*
b
; if
trans
=
'T'
or
'C'
and
equed
=
'C'
or
'B'
,
b
is overwritten by
diag(
c
)*
b
.
x
(local)
If
info
= 0
the
n
-by-
nrhs
solution matrix
X
to the original system of equations.
Note that
A
and
B
are modified on exit if
equed
'N'
, and the solution to the equilibrated system is
inv(diag(
sc
))*
X
if
trans
=
'N'
and
equed
=
'C'
or
'B'
, or
inv(diag(
sr
))*
X
if
trans
=
'T'
or
'C'
and
equed
=
'R'
or
'B'
.
rcond
(global)
An estimate of the reciprocal condition number of the matrix
A
after equilibration (if done). If
rcond
is less than the machine precision (in particular, if
rcond
=0
), the matrix is singular to working precision. This condition is indicated by a return code of
info
> 0
.
ferr
Arrays of size at least
max(
LOC
,
n_b
)
. The estimated forward error bounds for each solution vector
X
(
j
) (the
j
-th column of the solution matrix
X
). If
xtrue
is the true solution,
ferr
[
j
- 1]
bounds the magnitude of the largest entry in
(
X
(
j
) -
xtrue
)
divided by the magnitude of the largest entry in
X
(
j
). The quality of the error bound depends on the quality of the estimate of
norm(inv(
A
))
computed in the code; if the estimate of
norm(inv(
A
))
is accurate, the error bound is guaranteed.
berr
(local)
Arrays of size at least
max(
LOC
,
n_b
)
. The componentwise relative backward error of each solution vector
X
(
j
) (the smallest relative change in any entry of
A
or
B
that makes
X
(
j
) an exact solution).
work
[0]
(local) On exit,
work
[0]
returns the minimal and optimal
liwork
.
info
(global)
If
info
=0
, the execution is successful.
< 0
: if
info
= -
i
, the
i
-th argument had an illegal value
> 0
: if
info
=
i
,