SPRAL_SCALINGv1.0.0
Sparse Matrix Scalings
Fortran User Guide
This package generates various scalings
(and matchings) of real sparse matrices.
Given a symmetric matrix ,
it finds a diagonal matrix
such that the scaled matrix
has specific numerical properties.
Given a unsymmetric or rectangular matrix ,
it finds diagonal matrices
and
such that the scaled matrix
has specific numerical properties.
The specific numerical properties delivered depends on the algorithm used:
-
Matching-based
- algorithms scale
such that the maximum (absolute) value in each row and column of
is exactly 1.0, where the entries of maximum value form a maximum cardinality matching. The
Hungarian algorithm delivers an optimal matching slowly, whereas the auction algorithm delivers
an approximate matching quickly.
-
Norm-equilibration
- algorithms scale
such that the infinity norm of each row and column of
is
(for some user specified tolerance ).
Jonathan Hogg (STFC Rutherford Appleton Laboratory)
Major version history
-
2014-12-17 Version 1.0.0
- Initial public release
4.1 Installation
Please see the SPRAL install documentation.
4.2 Usage overview
4.2.1 Calling sequences
Access to the package requires a USE statement
use spral_scaling
The following procedures are available to the user:
- auction_scale_sym() and auction_scale_unsym() generate approximate matching-based scalings
for symmetric and unsymmetric/rectangular matrices respectively using an auction algorithm.
- equilib_scale_sym() and equilib_scale_unsym() generate norm-equilibration scalings for
symmetric and unsymmetric/rectangular matrices respectively.
- hungarian_scale_sym() and hungarian_scale_unsym() generate matching-based scalings for a
symmetric and unsymmetric/rectangular matrices respectively using the Hungarian algorithm.
4.2.2 Optional arguments
We use square brackets [ ] to indicate optional arguments. In each call, optional arguments appear last in
argument list. Since we reserve the right to add additional optional arguments in future releases of the
code, we strongly recommend that all optional arguments be called by keyword, not by
position.
4.2.3 Integer, real and package types
INTEGER denotes default INTEGER and INTEGER(long) denotes INTEGER(kind=selected_int_kind(18)).
REAL denotes double precision real. We also use the term package type to mean the same.
4.2.4 Data formats
Compressed Sparse Column (CSC) Format
This standard data format consists of the following data:
integer :: m ! number of rows (unsymmetric only)
integer :: n ! number of columns
integer, size(n+1) :: ptr ! column pointers
integer, size(ptr(n+1)-1) :: row ! row indices
real, size(ptr(n+1)-1) :: val ! numerical values
Non-zero matrix entries are ordered by increasing column index and stored in the arrays row(:) and val(:)
such that row(k) holds the row number and val(k) holds the value of the k-th entry. The ptr(:)
array stores column pointers such that ptr(i) is the position in row(:) and val(:) of the first entry
in the i-th column, and ptr(n+1) is one more than the total number of entries. There must be no
duplicate or out of range entries. Entries that are zero, including those on the diagonal, need not be
specified.
For symmetric matrices, only the lower triangular entries of
should
be supplied. For unsymmetric matrices, all entries in the matrix should be supplied.
Note that these routines offer no checking of user data, and the behaviour of these routines with misformatted
data is undefined.
To illustrate the CSC format, the following arrays describe the symmetric matrix shown in Figure 4.1.
n = 5
ptr(1:6) = (/ 1, 4, 5, 7, 9, 10 /)
row(1:9) = (/ 1, 2, 4, 3, 3, 5, 4, 5, 5 /)
val(1:9) = (/ 1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9 /)
4.3 Auction Algorithm
4.3.1 auction_scale_sym()
To generate a scaling for a symmetric matrix using an auction algorithm such that the entry of
maximum absolute value in each row and column is approximately 1.0,
call auction_scale_sym(n, ptr, row, val, scaling, options, inform[, match])
-
n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold the lower triangular part of
in compressed sparse column format as described in Section 4.2.4.
-
scaling(n)
- is an INTENT(OUT) array of package type. On exit, scaling(i) holds ,
the scaling corresponding to the i-th row and column.
-
options
- is an INTENT(IN) scalar of type auction_options. Its components specify the algorithmic options
used by the subroutine, as explained in Section 4.3.3.
-
inform
- is an INTENT(OUT) scalar of type auction_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.3.4.
-
match(n)
- is an optional INTENT(OUT) array of type INTEGER. If present, then on exit it specifies the matching
of rows to columns. Row i is matched to column match(i), or is unmatched if match(i)0.
4.3.2 auction_scale_unsym()
To generate a scaling for an unsymmetric or rectangular matrix using an auction algorithm such
that the entry of maximum absolute value in each row and column is approximately 1.0,
call auction_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform[,
match])
-
m, n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold
in compressed sparse column format as described in Section 4.2.4.
-
rscaling(m)
- is an INTENT(OUT) array of package type. On exit, rscaling(i) holds ,
the scaling corresponding to the i-th row.
-
cscaling(n)
- is an INTENT(OUT) array of package type. On exit, scaling(j) holds ,
the scaling corresponding to the j-th column.
-
options
- is an INTENT(IN) scalar of type auction_options. Its components specify the algorithmic options
used by the subroutine, as explained in Section 4.3.3.
-
inform
- is an INTENT(OUT) scalar of type auction_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.3.4.
-
match(m)
- is an optional INTENT(OUT) array of type INTEGER. If present, then on exit it specifies the matching
of rows to columns. Row i is matched to column match(i), or is unmatched if match(i)0.
4.3.3 type(auction_options)
The derived data type auction_options is used to specify the options used by the routines auction_scale_sym()
and auction_scale_unsym(). The components, that are automatically given default values in the definition of the
type, are:
-
eps_initial
- is a scalar of type default REAL that specifies the initial value of the minimum improvement
parameter
as described in Section 4.3.6.
-
max_iterations
- is a scalar of type INTEGER that specifies the maximum number of iterations the algorithm
may perform. The default is max_iterations=30000.
-
max_unchanged(3)
- is an array of type INTEGER that, together with min_proportion(:) specifies termination
conditions for the algorithm, as described in Section 4.3.6. The default is max_unchanged(:) = (/ 10,
100, 100 /).
-
min_proportion(3)
- is an array of type default REAL that, together with max_unchanged(:) specifies
termination conditions for the algorithm, as described in Section 4.3.6. The default is
min_proportion(:) = (/ 0.90, 0.0, 0.0 /).
4.3.4 type(auction_inform)
The derived data type auction_inform is used to hold parameters that give information about the progress of the
routines auction_scale_sym() and auction_scale_unsym(). The components are:
-
flag
- gives the exit status of the algorithm (details in Section 4.3.5).
-
iterations
- is a scalar of type INTEGER that holds the number of iterations performed.
-
matched
- is a scalar of type integer that holds the number of rows and columns that have been matched.
As the algorithm may terminate before a full matching is obtained, this only provides a lower bound
on the structural rank.
-
stat
- is a scalar of type INTEGER. In the event of an allocation error, it holds the Fortran stat parameter if
it is available (and is set to 0 otherwise).
-
unmatchable
- is a scalar of type integer that holds the number of columns designated as unmatchable. A
column is designated as unmatchable if there is no way to match it that improves the quality of the
matching. It provides an approximate lower bound on the structural rank deficiency.
4.3.5 Error Flags
A successful return from a routine is indicated by inform%flag having the value zero. A negative value is associated
with an error message.
Possible negative (error) values are:
-
-1
- Allocation error. If available, the Fortran stat parameter is returned in inform%stat.
4.3.6 Algorithm description
This algorithm finds a fast approximation to the matching and scaling produced by the HSL package MC64. If an
optimal matching is required, use the Hungarian algorithm instead. The algorithm works by solving the following
maximum product optimization problem using an auction algorithm. The scaling is derived from the dual variables
associated with the solution.
The array
gives a matching of rows to columns.
By using the transformation
where , the maximum product
problem in is replaced by a
minimum sum problem in
where all entries are positive. By standard optimization theory, there exist dual variables
and
corresponding to the constraints that satisfy the first order optimality conditions
To obtain a scaling we define scaling matrices
and
as
If a symmetric scaling is required, we average these as
By the first order optimality conditions, these scaling matrices guarantee that
To solve the minimum sum problem an auction algorithm is used. The algorithm is not guaranteed to find an optimal
matching. However it can find an approximate matching very quickly. A matching is maintained along with the row pricing
vector .
In each major iteration, we loop over each column in turn. If the column
is unmatched, we
calculate the value for
each entry and find the maximum across the column. If this maximum is positive, the current matching can be improved by matching
column with row
. This may mean that the previous
match of row now becomes unmatched.
We update the price of row ,
that is ,
to reflect this new benefit and continue to the next column.
To prevent incremental shuffling, we insist that the value of adding a new column is at least a threshold value
above zero, where
is based on the last iteration
in which row changed its
match. This is done by adding
to the price ,
where ,
where itr is the current iteration number.
The algorithm terminates if any of the following are satsified:
- All entries are matched.
- The number of major iterations exceeds options%max_iterations.
- At least options%max_unchanged(1) iterations have passed without the cardinality of the matching
increasing, and the proportion of matched columns is options%min_proportion(1).
- At least options%max_unchanged(2) iterations have passed without the cardinality of the matching
increasing, and the proportion of matched columns is options%min_proportion(2).
- At least options%max_unchanged(3) iterations have passed without the cardinality of the matching
increasing, and the proportion of matched columns is options%min_proportion(3).
The different combinations given by options%max_unchanged(1:3) and options%min_proportion(1:3)
allow a wide range of termination heuristics to be specified by the user depending on their
particular needs. Note that the matching and scaling produced will always be approximate as
is
non-zero.
Further details are given in the following paper:
- J.D. Hogg and J.A. Scott. (2014). On the efficient scaling of sparse symmetric matrices using an auction
algorithm. RAL Technical Report RAL-P-2014-002.
4.3.7 Example of auction_scale_sym()
The following code shows an example usage of auction_scale_sym().
! examples/Fortran/scaling/auction_sym.f90 - Example code for SPRAL_SCALING
program auction_scale_sym_example
use spral_scaling
use spral_matrix_util, only : print_matrix, &
SPRAL_MATRIX_REAL_SYM_INDEF
implicit none
! Derived types
type (auction_options) :: options
type (auction_inform) :: inform
! Parameters
integer, parameter :: wp = kind(0.0d0)
! Matrix data
integer :: n, ptr(6), row(8)
real(wp) :: val(8)
! Other variables
integer :: match(5), i, j
real(wp) :: scaling(5)
! Data for symmetric matrix:
! ( 2 1 )
! ( 1 4 1 8 )
! ( 1 3 2 )
! ( 2 )
! ( 8 2 )
n = 5
ptr(1:n+1) = (/ 1, 3, 6, 8,8, 9 /)
row(1:ptr(n+1)-1) = (/ 1, 2, 2, 3, 5, 3, 4, 5 /)
val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 /)
write(*, "(a)") "Initial matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)
! Perform symmetric scaling
call auction_scale_sym(n, ptr, row, val, scaling, options, inform, &
match=match)
if(inform%flag<0) then
write(*, "(a, i5)") "auction_scale_sym() returned with error ", &
inform%flag
stop
endif
! Print scaling and matching
write(*,"(a,10i10)") ’Matching:’, match(1:n)
write(*,"(a,10es10.2)") ’Scaling: ’, scaling(1:n)
! Calculate scaled matrix and print it
do i = 1, n
do j = ptr(i), ptr(i+1)-1
val(j) = scaling(i) * val(j) * scaling(row(j))
end do
end do
write(*, "(a)") "Scaled matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)
end program auction_scale_sym_example
The above code produces the following output.
Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1: 2.0000E+00 1.0000E+00
2: 1.0000E+00 4.0000E+00 1.0000E+00 8.0000E+00
3: 1.0000E+00 3.0000E+00 2.0000E+00
4: 2.0000E+00
5: 8.0000E+00 2.0000E+00
Matching: 1 5 4 3 2
Scaling: 7.07E-01 1.62E-01 2.78E-01 1.80E+00 7.72E-01
Scaled matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1: 1.0000E+00 1.1443E-01
2: 1.1443E-01 1.0476E-01 4.5008E-02 1.0000E+00
3: 4.5008E-02 2.3204E-01 1.0000E+00
4: 1.0000E+00
5: 1.0000E+00 1.1932E+00
4.4 Norm-equilibration algorithm
4.4.1 equilib_scale_sym()
To generate a scaling for a symmetric matrix using a norm equilibration
algorithm such that the infinity norm of each row and column is equal to
,
call equilib_scale_sym(n, ptr, row, val, scaling, options, inform)
-
n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold the lower triangular part of
in compressed sparse column format as described in Section 4.2.4.
-
scaling(n)
- is an INTENT(OUT) array of package type. On exit, scaling(i) holds ,
the scaling corresponding to the i-th row and column.
-
options
- is an INTENT(IN) scalar of type equilib_options. Its components specify the algorithmic options
used by the subroutine, as explained in Section 4.4.3.
-
inform
- is an INTENT(OUT) scalar of type equilib_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.4.4.
4.4.2 equilib_scale_unsym()
To generate a scaling for an unsymmetric or rectangular matrix using a norm
equilibration algorithm such that the infinity norm of each row and column is equal to
,
call equilib_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform)
-
m, n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold
in compressed sparse column format as described in Section 4.2.4.
-
rscaling(m)
- is an INTENT(OUT) array of package type. On exit, rscaling(i) holds ,
the scaling corresponding to the i-th row.
-
cscaling(n)
- is an INTENT(OUT) array of package type. On exit, cscaling(j) holds ,
the scaling corresponding to the j-th column.
-
options
- is an INTENT(IN) scalar of type equilib_options. Its components specify the algorithmic options
used by the subroutine, as explained in Section 4.4.3.
-
inform
- is an INTENT(OUT) scalar of type equilib_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.4.4.
4.4.3 type(equilib_options)
The derived data type equilib_options is used to specify the options used by the routine equilib_scale_sym().
The components, that are automatically given default values in the definition of the type, are:
-
max_iterations
- is a scalar of type INTEGER that specifies the maximum number of iterations the algorithm
may perform. The default is max_iterations=10.
-
tol
- is a scalar of type default REAL that specifies the convergence tolerance
for the algorithm (though often termination is based on max_iterations). The default is tol = 1e-8.
4.4.4 type(equilib_inform)
The derived data type equilib_inform is used to hold parameters that give information about the progress of the
routines equilib_scale_sym() and equilib_scale_unsym(). The components are:
-
flag
- gives the exit status of the algorithm (details in Section 4.4.5).
-
iterations
- is a scalar of type INTEGER. On exit, it holds the number of iterations performed.
-
stat
- is a scalar of type INTEGER. In the event of an allocation error, it holds the Fortran stat parameter if
it is available (and is set to 0 otherwise).
4.4.5 Error Flags
A successful return from a routine is indicated by inform%flag having the value zero. A negative value is associated
with an error message.
Possible negative (error) values are:
-
-1
- Allocation error. If available, the Fortran stat parameter is returned in inform%stat.
4.4.6 Algorithm description
This algorithm is very similar to that used by the HSL routine MC77. An iterative method is used to
scale the infinity norm of both rows and columns to 1 with an asymptotic linear rate of convergence of
,
preserving symmetry if the matrix is symmetric.
For unsymmetric matrices, the algorithm outline is as follows:
for do
if()
exit
end for
For symmetric matrices,
is symmetric, so
,
and some operations can be skipped.
Further details are given in the following paper:
- P. Knight, D. Ruiz and B. Ucar. (2012). A symmetry preserving algorithm for matrix scaling. INRIA
Research Report 7552.
4.4.7 Example of equilib_scale_sym()
The following code shows an example usage of equilib_scale_sym().
! examples/Fortran/scaling/equilib_sym.f90 - Example code for SPRAL_SCALING
program equilib_scale_sym_example
use spral_scaling
use spral_matrix_util, only : print_matrix, &
SPRAL_MATRIX_REAL_SYM_INDEF
implicit none
! Derived types
type (equilib_options) :: options
type (equilib_inform) :: inform
! Parameters
integer, parameter :: wp = kind(0.0d0)
! Matrix data
integer :: n, ptr(6), row(8)
real(wp) :: val(8)
! Other variables
integer :: i, j
real(wp) :: scaling(5)
! Data for symmetric matrix:
! ( 2 1 )
! ( 1 4 1 8 )
! ( 1 3 2 )
! ( 2 )
! ( 8 2 )
n = 5
ptr(1:n+1) = (/ 1, 3, 6, 8,8, 9 /)
row(1:ptr(n+1)-1) = (/ 1, 2, 2, 3, 5, 3, 4, 5 /)
val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 /)
write(*, "(a)") "Initial matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)
! Perform symmetric scaling
call equilib_scale_sym(n, ptr, row, val, scaling, options, inform)
if(inform%flag<0) then
write(*, "(a, i5)") "equilib_scale_sym() returned with error ", &
inform%flag
stop
endif
! Print scaling and matching
write(*,"(a,10es10.2)") ’Scaling: ’, scaling(1:n)
! Calculate scaled matrix and print it
do i = 1, n
do j = ptr(i), ptr(i+1)-1
val(j) = scaling(i) * val(j) * scaling(row(j))
end do
end do
write(*, "(a)") "Scaled matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val)
end program equilib_scale_sym_example
The above code produces the following output.
Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1: 2.0000E+00 1.0000E+00
2: 1.0000E+00 4.0000E+00 1.0000E+00 8.0000E+00
3: 1.0000E+00 3.0000E+00 2.0000E+00
4: 2.0000E+00
5: 8.0000E+00 2.0000E+00
Scaling: 7.07E-01 3.54E-01 5.77E-01 8.66E-01 3.54E-01
Scaled matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
1: 1.0000E+00 2.5000E-01
2: 2.5000E-01 5.0000E-01 2.0412E-01 1.0000E+00
3: 2.0412E-01 1.0000E+00 9.9960E-01
4: 9.9960E-01
5: 1.0000E+00 2.5000E-01
4.5 Hungarian algorithm
4.5.1 hungarian_scale_sym()
To generate a scaling for a symmetric matrix using the Hungarian algorithm such that the entry of
maximum absolute value in each row and column is 1.0,
call hungarian_scale_sym(n, ptr, row, val, scaling, options, inform[, match])
-
n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold the lower triangular part of
in compressed sparse column format as described in Section 4.2.4.
-
scaling(n)
- is an INTENT(OUT) array of package type. On exit, scaling(i) holds ,
the scaling corresponding to the i-th row and column.
-
options
- is an INTENT(IN) scalar of type hungarian_options. Its components specify the algorithmic
options used by the subroutine, as explained in Section 4.5.3.
-
inform
- is an INTENT(OUT) scalar of type hungarian_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.5.4.
-
match(m)
- is an optional INTENT(OUT) array of type INTEGER. If present, then on exit it specifies the matching
of rows to columns. Row i is matched to column match(i), or is unmatched if match(i)0.
4.5.2 hungarian_scale_unsym()
To generate a scaling for an unsymmetric or rectangular matrix using the Hungarian
algorithm such that the entry of maximum absolute value in each row and column is 1.0,
call hungarian_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, inform[, match])
-
m, n, ptr(:), row(:), val(:)
- are INTENT(IN) variables that must hold the lower triangular part of
in compressed sparse column format as described in Section 4.2.4.
-
rscaling(m)
- is an INTENT(OUT) array of package type. On exit, rscaling(i) holds ,
the scaling corresponding to the i-th row.
-
rscaling(n)
- is an INTENT(OUT) array of package type. On exit, rscaling(j) holds ,
the scaling corresponding to the j-th column.
-
options
- is an INTENT(IN) scalar of type hungarian_options. Its components specify the algorithmic
options used by the subroutine, as explained in Section 4.5.3.
-
inform
- is an INTENT(OUT) scalar of type hungarian_inform. On exit, its components provide information
about the execution of the subroutine, as explained in Section 4.5.4.
-
match(m)
- is an optional INTENT(OUT) array of type INTEGER. If present, then on exit it specifies the matching
of rows to columns. Row i is matched to column match(i), or is unmatched if match(i)0.
4.5.3 type(hungarian_options)
The derived data type hungarian_options is used to specify the options used by the routines
hungarian_scale_sym() and hungarian_scale_unsym(). The components, that are automatically given default
values in the definition of the type, are:
-
scale_if_singular
- is a scalar of type default LOGICAL that specifies whether scaling should continue if the
matrix
is found to be structurally singular. If scale_if_singular.true.,
and the
is structurally singular, a partial scaling corresponding to a maximum cardinality matching will be
returned and a warning issued. Otherwise, an identity scaling will be returned and an error issued.
4.5.4 type(hungarian_inform)
The derived data type hungarian_inform is used to hold parameters that give information about the
progress of the routine hungarian_scale_sym() and hungarian_scale_unsym(). The components
are:
-
flag
- gives the exit status of the algorithm (details in Section 4.5.5).
-
matched
- is a scalar of type INTEGER that holds the number of rows and columns that have been matched
(i.e. the structural rank).
-
stat
- is a scalar of type INTEGER. In the event of an allocation error, it holds the Fortran stat parameter if
it is available (and is set to 0 otherwise).
4.5.5 Error Flags
A successful return from a routine is indicated by inform%flag having the value zero. A negative value is associated
with an error message and a positive value with a warning.
Possible negative (error) values are:
-
-1
- Allocation error. If available, the Fortran stat parameter is returned in inform%stat.
-
-2
- Matrix
is structurally rank-deficient. This error is returned only if options%scale_if_singular.false..
The scaling vector is set to 1.0 and a matching of maximum cardinality returned in the optional
argument match(:), if present.
Possible positive (warning) values are:
-
+1
- Matrix
is structurally rank-deficient. This warning is returned only if options%scale_if_singular.true..
4.5.6 Algorithm description
This algorithm is the same as used by the HSL package MC64. A scaling is derived from dual variables
found during the solution of the below maximum product optimization problem using the Hungarian
algorithm.
The array
gives a matching of rows to columns.
By using the transformation
where , the maximum product
problem in is replaced by a
minimum sum problem in
where all entries are positive. By standard optimization theory, there exist dual variables
and
corresponding to the constraints that satisfy the first order optimality conditions
To obtain a scaling we define scaling matrices
and
as
If a symmetric scaling is required, we average these as
By the first order optimality conditions, these scaling matrices guarantee that
To solve the minimum sum problem, the Hungarian algorithm maintains an optimal matching on a subset of the
rows and columns. It proceeds to grow this set by finding augmenting paths from an unmatched row to an
unmatched column. The algorithm is guaranteed to find the optimal solution in a fixed number of steps, but can be
very slow as it may need to explore the full matrix a number of times equal to the dimension of the matrix. To
minimize the solution time, a warmstarting heuristic is used to construct an initial optimal subset
matching.
Further details are given in the following paper:
- I.S. Duff and J. Koster. (1997). The design and use of algorithms for permuting large entries to the
diagonal of sparse matrices. SIAM J. Matrix Anal. Applics. 20(4), pp 889–901.
4.5.7 Example usage of hungarian_scale_unsym()
The following code shows an example usage of hungarian_scale_unsym().
! examples/Fortran/scaling/hungarian_unsym.f90 - Example code for SPRAL_SCALING
program hungarian_scale_unsym_example
use spral_scaling
use spral_matrix_util, only : print_matrix, &
SPRAL_MATRIX_REAL_UNSYM
implicit none
! Derived types
type (hungarian_options) :: options
type (hungarian_inform) :: inform
! Parameters
integer, parameter :: wp = kind(0.0d0)
! Matrix data
integer :: m, n, ptr(6), row(10)
real(wp) :: val(10)
! Other variables
integer :: match(5), i, j
real(wp) :: rscaling(5), cscaling(5)
! Data for unsymmetric matrix:
! ( 2 5 )
! ( 1 4 7 )
! ( 1 2 )
! ( 3 )
! ( 8 2 )
m = 5; n = 5
ptr(1:n+1) = (/ 1, 3, 7, 8, 9, 11 /)
row(1:ptr(n+1)-1) = (/ 1, 2, 1, 2, 3, 5, 4, 3, 2, 5 /)
val(1:ptr(n+1)-1) = (/ 2.0, 1.0, 5.0, 4.0, 1.0, 8.0, 3.0, 2.0, 7.0, 2.0 /)
write(*, "(a)") "Initial matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_UNSYM, m, n, ptr, row, val)
! Perform unsymmetric scaling
call hungarian_scale_unsym(m, n, ptr, row, val, rscaling, cscaling, options, &
inform, match=match)
if(inform%flag<0) then
write(*, "(a, i5)") "hungarian_scale_unsym() returned with error ", &
inform%flag
stop
endif
! Print scaling and matching
write(*,"(a,10i10)") ’Matching:’, match(1:m)
write(*,"(a,10es10.2)") ’Row Scaling: ’, rscaling(1:m)
write(*,"(a,10es10.2)") ’Col Scaling: ’, cscaling(1:n)
! Calculate scaled matrix and print it
do i = 1, n
do j = ptr(i), ptr(i+1)-1
val(j) = rscaling(row(j)) * val(j) * cscaling(i)
end do
end do
write(*, "(a)") "Scaled matrix:"
call print_matrix(6, -1, SPRAL_MATRIX_REAL_UNSYM, m, n, ptr, row, val)
end program hungarian_scale_unsym_example
The above code produces the following output.
Initial matrix:
Real unsymmetric matrix, dimension 5x5 with 10 entries.
1: 2.0000E+00 5.0000E+00
2: 1.0000E+00 4.0000E+00 7.0000E+00
3: 1.0000E+00 2.0000E+00
4: 3.0000E+00
5: 8.0000E+00 2.0000E+00
Matching: 1 5 4 3 2
Row Scaling: 5.22E-01 5.22E-01 5.22E-01 5.22E-01 5.22E-01
Col Scaling: 9.59E-01 2.40E-01 6.39E-01 9.59E-01 2.74E-01
Scaled matrix:
Real unsymmetric matrix, dimension 5x5 with 10 entries.
1: 1.0000E+00 6.2500E-01
2: 5.0000E-01 5.0000E-01 1.0000E+00
3: 1.2500E-01 1.0000E+00
4: 1.0000E+00
5: 1.0000E+00 2.8571E-01