SCALING - Sparse matrix scalings

#include <spral_scaling.h> /* or <spral.h> for all packages */

Purpose

This package generates various scalings (and matchings) of real sparse matrices.

Given a symmetric matrix \(A\), it finds a diagonal matrix \(D\) such that the scaled matrix

\[\hat{A} = DAD\]

has specific numerical properties.

Given an unsymmetric or rectangular matrix \(A\), it finds diagonal matrices \(D_r\) and \(D_c\) such that the scaled matrix

\[\hat{A} = D_r A D_c\]

has specific numerical properties.

The specific numerical properties delivered depends on the algorithm used:

Matching-based

algorithms scale \(A\) such that the maximum (absolute) value in each row and column of \(\hat{A}\) 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 \(A\) such that the infinity norm of each row and column of \(\hat{A}\) is \(1.0\pm \tau\) (for some user specified tolerance \(\tau\)).

Auction Algorithm

Routines

void spral_scaling_auction_default_options(struct spral_scaling_auction_options *options)

Intialises members of options structure to default values.

Parameters:

• options – Structure to be initialised.

void spral_scaling_auction_sym(int n, const int *ptr, const int *row, const double *val, double *scaling, int *match, const struct spral_scaling_auction_options *options, struct spral_scaling_auction_inform *inform)

Find a matching-based symmetric scaling using the auction algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is (approximately) \(1.0\).

Parameters:

• n – number of columns in \(A\).

• ptr[n+1] – columns pointers for \(A\) (see CSC format).

• row[ptr[n]] – row indices for \(A\) (see CSC format).

• val[ptr[n]] – non-zero values for \(A\) (see CSC format).

• scaling[n] – returns scaling found by routine.

• match – may be NULL; otherwise, an array of size n to output the matching found by routine. Row i is matched to column match[i], or is unmatched if match[i]=0.

• options – controls behaviour of routine.

• inform – returns information on execution of routine.

void spral_scaling_auction_sym_long(int n, const int64_t *ptr, const int *row, const double *val, double *scaling, int *match, const struct spral_scaling_auction_options *options, struct spral_scaling_auction_inform *inform)

As spral_scaling_auction_sym(), except ptr has type int64_t.

void spral_scaling_auction_unsym(int m, int n, const int *ptr, const int *row, const double *val, double *rscaling, double *cscaling, int *match, const struct spral_scaling_auction_options *options, struct spral_scaling_auction_inform *inform)

Find a matching-based unsymmetric scaling using the auction algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is (approximately) \(1.0\).

Parameters:

• m – number of rows in \(A\)

• n – number of columns in \(A\)

• ptr[n+1] – columns pointers for \(A\) (see CSC format)

• row[ptr[n]] – row indices for \(A\) (see CSC format)

• val[ptr[n]] – non-zero values for \(A\) (see CSC format)

• rscaling[m] – returns row scaling found by routine

• cscaling[n] – returns column scaling found by routine

• match – may be NULL; otherwise, an array of size m to output the matching found by routine. Row i is matched to column match[i], or is unmatched if match[i]=0.

• options – controls behaviour of routine

• inform – returns information on execution of routine

void spral_scaling_auction_unsym_long(int m, int n, const int64_t *ptr, const int *row, const double *val, double *rscaling, double *cscaling, int *match, const struct spral_scaling_auction_options *options, struct spral_scaling_auction_inform *inform)

As spral_scaling_auction_unsym(), except ptr has type int64_t.

Data-types

struct spral_scaling_auction_options

Used to specify options to the routines spral_scaling_auction_sym() and spral_scaling_auction_unsym(). The routine spral_scaling_auction_default_options() may be used to intialise with default values.

Please refer to the method section for details on how these parameters are used.

int array_base

Indexing base for arrays. Either 0 (C indexing) or 1 (Fortran indexing). Default is 0.

float eps_initial

Initial value of improvement parameter \(\epsilon\). Default is 0.01.

int max_iterations

Maximum number of iterations. Default is 30000.

int max_unchanged[3]

Together with min_proportion[], specifies termination conditions. Default is {10, 100, 100}.

float min_proportion[3]

Together with max_unchanged[], specifies termination conditions. Default is {0.9, 0.0, 0.0}.

struct spral_scaling_auction_inform

Used to return information about the execution of the algorithm.

int flag

Gives the exit status of the algorithm (see table below)

int iterations

Number of iterations performed.

int matched

Number of rows and columns that have been matched.

int stat

Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

int unmatchable

Number of columns designated as unmatchable (there is no way to match it that improves the quality of the matching).

Note: As the algorithm may terminate before a full matching is obtained, spral_scaling_auction_inform.matched provides only a lower bound on the structural rank. However, spral_scaling_auction_inform.unmatchable provides an approximate lower bound on the structural rank deficiency.

inform.flag	Return status
0	Success.
-1	Allocation error. Fortran stat value is returned in spral_scaling_auction_inform.stat.

Example

The following code shows an example usage of auction_scale_sym().

/* examples/C/scaling/auction_sym.c - Example code for SPRAL_SCALING */
#include <stdlib.h>
#include <stdio.h>
#include "spral.h"

int main(void) {
   /* Derived types */
   struct spral_scaling_auction_options options;
   struct spral_scaling_auction_inform inform;

   /* Other variables */
   int match[5];
   double scaling[5];

   /* Data for symmetric matrix:
    * ( 2  1         )
    * ( 1  4  1    8 )
    * (    1  3  2   )
    * (       2      )
    * (    8       2 ) */
   int n = 5;
   int ptr[]    = { 0,        2,             5,      7,7,   8 };
   int row[]    = { 0,   1,   1,   2,   4,   2,   3,   4   };
   double val[] = { 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 };
   printf("Initial matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Perform symmetric scaling */
   spral_scaling_auction_default_options(&options);
   spral_scaling_auction_sym(n, ptr, row, val, scaling, match, &options, &inform);
   if(inform.flag<0) {
      printf("spral_scaling_auction_sym() returned with error %5d", inform.flag);
      exit(1);
   }

   /* Print scaling and matching */
   printf("Matching:");
   for(int i=0; i<n; i++) printf(" %10d", match[i]);
   printf("\nScaling: ");
   for(int i=0; i<n; i++) printf(" %10.2le", scaling[i]);
   printf("\n");

   /* Calculate scaled matrix and print it */
   for(int i=0; i<n; i++) {
      for(int j=ptr[i]; j<ptr[i+1]; j++)
         val[j] = scaling[row[j]] * val[j] * scaling[i];
   }
   printf("Scaled matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Success */
   return 0;
}

The above code produces the following output:

Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
0:   2.0000E+00   1.0000E+00
1:   1.0000E+00   4.0000E+00   1.0000E+00                8.0000E+00
2:                1.0000E+00   3.0000E+00   2.0000E+00
3:                             2.0000E+00
4:                8.0000E+00                             2.0000E+00
Matching:          0         4         3         2         1
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.
0:   1.0000E+00   1.1443E-01
1:   1.1443E-01   1.0476E-01   4.5008E-02                1.0000E+00
2:                4.5008E-02   2.3204E-01   1.0000E+00
3:                             1.0000E+00
4:                1.0000E+00                             1.1932E+00

Method

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.

\[\begin{split}\max_{\sigma} & \prod_{i=1}^m\prod_{j=1}^n |a_{ij}|\sigma_{ij} & \\ \mathrm{s.t.} & \sum_{i=1}^m\sigma_{ij} = 1, & \forall j=1,n \\ & \sum_{j=1}^n\sigma_{ij} = 1, & \forall i=1,m \\ & \sigma_{ij} \in \{0,1\}.\end{split}\]

The array \(\sigma\) gives a matching of rows to columns.

By using the transformation

\[w_{ij} = \log c_j - \log |a_{ij}|,\]

where \(c_j = \max_i |a_{ij}|\), the maximum product problem in \(a_{ij}\) is replaced by a minimum sum problem in \(w_{ij}\) where all entries are positive. By standard optimization theory, there exist dual variables \(u\) and \(v\) corresponding to the constraints that satisfy the first order optimality conditions

\[\begin{split}w_{ij} - u_i - v_j = 0, & \mbox{ if } \sigma_{ij }=1, \\ w_{ij} - u_i - v_j \ge 0, & \mbox{ if } \sigma_{ij }=0.\end{split}\]

To obtain a scaling we define scaling matrices \(D_r\) and \(D_c\) as

\[ \begin{align}\begin{aligned}d^r_i = e^{u_i},\\d^c_i = e^{v_i}.\end{aligned}\end{align} \]

If a symmetric scaling is required, we average these as

\[d_i = \sqrt{d^r_id^c_i}.\]

By the first order optimality conditions, these scaling matrices guarantee that

\[\begin{split}d^r_i|a_{ij}|d^c_j = 1, && \mbox{if } \sigma_{ij}=1, \\ d^r_i|a_{ij}|d^c_j \le 1, && \mbox{if } \sigma_{ij}=0.\end{split}\]

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 \(u\). In each major iteration, we loop over each column in turn. If the column \(j\) is unmatched, we calculate the value \(p_i = w_{ij} - u_i\) for each entry and find the maximum across the column. If this maximum is positive, the current matching can be improved by matching column \(j\) with row \(i\). This may mean that the previous match of row \(i\) now becomes unmatched. We update the price of row \(i\), that is \(u_i\), 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 \(\epsilon\) above zero, where \(\epsilon\) is based on the last iteration in which row \(i\) changed its match. This is done by adding \(\epsilon\) to the price \(u_i\), where \(\epsilon = \mathrm{options.eps_initial} + \mathrm{itr} / (n+1)\), 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[0] iterations have passed without the cardinality of the matching increasing, and the proportion of matched columns is options.min_proportion[0].

• 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].

The different combinations given by options.max_unchanged[] and options.min_proportion[] 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 \(\epsilon\) is non-zero.

Further details are given in the following paper:

[1]

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. [STFC TR]

Norm-equilibration Algorithm

Routines

void spral_scaling_equilib_default_options(struct spral_scaling_equilib_options *options)

Intialises members of options structure to default values.

Parameters:

• options – Structure to be initialised.

void spral_scaling_equilib_sym(int n, const int *ptr, const int *row, const double *val, double *scaling, const struct spral_scaling_equilib_options *options, struct spral_scaling_equilib_inform *inform)

Find a matching-based symmetric scaling using the norm-equilibration algorithm.

The scaled matrix is such that the infinity norm of each row and column are equal to \(1.0\).

Parameters:

• n – number of columns in \(A\).

• ptr[n] – columns pointers for \(A\) (see CSC format).

• row[ptr[n]] – row indices for \(A\) (see CSC format).

• val[ptr[n]] – non-zero values for \(A\) (see CSC format).

• scaling[n] – returns scaling found by routine.

• options – controls behaviour of routine.

• inform – returns information on execution of routine.

void spral_scaling_equilib_sym_long(int n, const int64_t *ptr, const int *row, const double *val, double *scaling, const struct spral_scaling_equilib_options *options, struct spral_scaling_equilib_inform *inform)

As spral_scaling_equilib_sym(), except ptr has type int64_t.

void spral_scaling_equilib_unsym(int m, int n, const int *ptr, const int *row, const double *val, double *rscaling, double *cscaling, const struct spral_scaling_equilib_options *options, struct spral_scaling_equilib_inform *inform)

Find a matching-based unsymmetric scaling using the norm-equilibration algorithm.

The scaled matrix is such that the infinity norm of each row and column are equal to \(1.0\).

Parameters:

• m – number of rows in \(A\).

• n – number of columns in \(A\).

• ptr[n+1] – columns pointers for \(A\) (see CSC format).

• row[ptr[n]] – row indices for \(A\) (see CSC format).

• val[ptr[n]] – non-zero values for \(A\) (see CSC format).

• rscaling[m] – returns row scaling found by routine.

• cscaling[n] – returns column scaling found by routine.

• options – controls behaviour of routine.

• inform – returns information on execution of routine.

void spral_scaling_equilib_unsym_long(int m, int n, const int64_t *ptr, const int *row, const double *val, double *rscaling, double *cscaling, const struct spral_scaling_equilib_options *options, struct spral_scaling_equilib_inform *inform)

As spral_scaling_equilib_unsym(), except ptr has type int64_t.

Data-types

struct spral_scaling_equilib_options

Used to specify options to the routines spral_scaling_equilib_sym() and spral_scaling_equilib_unsym(). The routine spral_scaling_equilib_default_options() may be used to intialise with default values.

Please refer to the method section for details on how these parameters are used.

int array_base

Indexing base for arrays. Either 0 (C indexing) or 1 (Fortran indexing). Default is 0.

int max_iterations

Maximum number of iterations. Default is 10.

float tol

Convergence tolerance \(\tau\). Default is 1e-8.

struct spral_scaling_equilib_inform

Used to return information about the execution of the algorithm.

int flag

Gives the exit status of the algorithm (see table below).

int iterations

Number of iteration performed.

int stat

Holds the Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

inform.flag	Return status
0	Success.
-1	Allocation error. Fortran stat value is returned in inform.stat.

Example

The following code shows an example usage of spral_scaling_equilib_sym().

/* examples/C/scaling/equilib_sym.c - Example code for SPRAL_SCALING */
#include <stdlib.h>
#include <stdio.h>
#include "spral.h"

int main(void) {
   /* Derived types */
   struct spral_scaling_equilib_options options;
   struct spral_scaling_equilib_inform inform;

   /* Other variables */
   double scaling[5];

   /* Data for symmetric matrix:
    * ( 2  1         )
    * ( 1  4  1    8 )
    * (    1  3  2   )
    * (       2      )
    * (    8       2 ) */
   int n = 5;
   int ptr[]    = { 0,        2,             5,      7,7,   8 };
   int row[]    = { 0,   1,   1,   2,   4,   2,   3,   4   };
   double val[] = { 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 };
   printf("Initial matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Perform symmetric scaling */
   spral_scaling_equilib_default_options(&options);
   spral_scaling_equilib_sym(n, ptr, row, val, scaling, &options, &inform);
   if(inform.flag<0) {
      printf("spral_scaling_equilib_sym() returned with error %5d", inform.flag);
      exit(1);
   }

   /* Print scaling */
   printf("Scaling: ");
   for(int i=0; i<n; i++) printf(" %10.2le", scaling[i]);
   printf("\n");

   /* Calculate scaled matrix and print it */
   for(int i=0; i<n; i++) {
      for(int j=ptr[i]; j<ptr[i+1]; j++)
         val[j] = scaling[row[j]] * val[j] * scaling[i];
   }
   printf("Scaled matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Success */
   return 0;
}

The above code produces the following output:

Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
0:   2.0000E+00   1.0000E+00
1:   1.0000E+00   4.0000E+00   1.0000E+00                8.0000E+00
2:                1.0000E+00   3.0000E+00   2.0000E+00
3:                             2.0000E+00
4:                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.
0:   1.0000E+00   2.5000E-01
1:   2.5000E-01   5.0000E-01   2.0412E-01                1.0000E+00
2:                2.0412E-01   1.0000E+00   9.9960E-01
3:                             9.9960E-01
4:                1.0000E+00                             2.5000E-01

Method

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.0\) with an asymptotic linear rate of convergence of \(\frac{1}{2}\), preserving symmetry if the matrix is symmetric.

For unsymmetric matrices, the algorithm outline is as follows:

• \(D_r^{(0)} = I, D_c^{(0)}=I\)

• for (\(k=1,\) options.max_iterations)

  • \(A^{(k-1)} = D_r^{(k-1)} A D_c^{(k-1)}\)

  • \((D_r^{(k)})_{ii} = (D_r^{(k-1)})_{ii}\; /\; \sqrt{\max_j(A^{(k-1)})_{ij}}\)

  • \((D_c^{(k)})_{jj} = (D_c^{(k-1)})_{jj}\; /\; \sqrt{\max_i(A^{(k-1)})_{ij}}\)

  • if ( \(|1-\|A^{(k-1)}\|_{\max}|\le\) options.tol ) exit

For symmetric matrices, \(A^{(k-1)}\) is symmetric, so \(D_r^{(k)} = D_c^{(k)}\), and some operations can be skipped.

Further details are given in the following paper:

[2]

16. Knight, D. Ruiz and B. Ucar. (2012). A symmetry preserving algorithm for matrix scaling. INRIA Research Report 7552. [INRIA TR]

Hungarian Algorithm

Routines

void spral_scaling_hungarian_default_options(struct spral_scaling_hungarian_options *options)

Intialises members of options structure to default values.

Parameters:

• options – Structure to be initialised.

void spral_scaling_hungarian_sym(int n, const int *ptr, const int *row, const double *val, double *scaling, int *match, const struct spral_scaling_hungarian_options *options, struct spral_scaling_hungarian_inform *inform)

Find a matching-based symmetric scaling using the Hungarian algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is \(1.0\).

Parameters:

• n – number of columns in \(A\).

• ptr[n+1] – columns pointers for \(A\) (see CSC format).

• row[ptr[n]] – row indices for \(A\) (see CSC format).

• val[ptr[n]] – non-zero values for \(A\) (see CSC format).

• scaling[n] – returns scaling found by routine.

• match – may be NULL; otherwise, an array of size n to output the matching found by routine. Row i is matched to column match[i], or is unmatched if match[i]=0.

• options – controls behaviour of routine.

• inform – returns information on execution of routine.

void spral_scaling_hungarian_sym_long(int n, const int64_t *ptr, const int *row, const double *val, double *scaling, int *match, const struct spral_scaling_hungarian_options *options, struct spral_scaling_hungarian_inform *inform)

As spral_scaling_hungarian_sym(), except ptr has type int64_t.

void spral_scaling_hungarian_unsym(int m, int n, const int *ptr, const int *row, const double *val, double *rscaling, double *cscaling, int *match, const struct spral_scaling_hungarian_options *options, struct spral_scaling_hungarian_inform *inform)

Find a matching-based symmetric scaling using the Hungarian algorithm.

The scaled matrix is such that the entry of maximum absolute value in each row and column is \(1.0\).

Parameters:

• m – number of rows in \(A\).

• n – number of columns in \(A\).

• ptr[n+1] – columns pointers for \(A\) (see CSC format).

• row[ptr[n]] – row indices for \(A\) (see CSC format).

• val[ptr[n]] – non-zero values for \(A\) (see CSC format).

• rscaling[m] – returns row scaling found by routine.

• cscaling[n] – returns column scaling found by routine.

• match – may be NULL; otherwise, an array of size n to output the matching found by routine. Row i is matched to column match[i], or is unmatched if match[i]=0.

• options – controls behaviour of routine.

• inform – returns information on execution of routine.

void spral_scaling_hungarian_unsym_long(int m, int n, const int64_t *ptr, const int *row, const double *val, double *rscaling, double *cscaling, int *match, const struct spral_scaling_hungarian_options *options, struct spral_scaling_hungarian_inform *inform)

As spral_scaling_hungarian_unsym(), except ptr has type int64_t.

Data-types

struct spral_scaling_hungarian_options

Used to specify options to the routines spral_scaling_hungarian_sym() and spral_scaling_hungarian_unsym(). The routine spral_scaling_hungarian_default_options() may be used to intialise with default values.

Please refer to the method section for details on how these parameters are used.

int array_base

Indexing base for arrays. Either 0 (C indexing) or 1 (Fortran indexing). Default is 0.

bool scale_if_singular

Behaviour for structurally singular matrices. If true, a partial scaling corresponding to a maximum cardinality matching will be returned. If false, an identity scaling is returned with an error code. Default is false.

Warning

If options.scale_if_singular=true, the resulting scaling will only be maximal for the matched rows/columns, and extreme care should be taken to ensure its use is meaningful!

type  hungarian_inform

Used to return information about the execution of the algorithm.

int flag

Exit status of the algorithm (see table below)

int matched

Number of rows and columns that have been matched.

int stat

Holds the Fortran stat parameter in the event of an allocation failure (set to 0 otherwise).

Note: The number matched gives the structural rank of the matrix.

inform.flag	Return status
0	Success.
+1	Warning: Matrix is structurally rank-deficient. Only returned if options.scale_if_singular=true.
-1	Error: Allocation failed. Fortran stat value is returned in inform.stat.
-2	Error: Matrix is structurally rank-deficient. Only returned if options.scale_if_singular=false. Scaling vector(s) will be set to the identity, and a maximum cardinality matching will be returned in match[] (if present).

Example

The following code shows an example usage of spral_scaling_hungarian_sym().

/* examples/C/scaling/hungarian_sym.c - Example code for SPRAL_SCALING */
#include <stdlib.h>
#include <stdio.h>
#include "spral.h"

int main(void) {
   /* Derived types */
   struct spral_scaling_hungarian_options options;
   struct spral_scaling_hungarian_inform inform;

   /* Other variables */
   int match[5];
   double scaling[5];

   /* Data for symmetric matrix:
    * ( 2  1         )
    * ( 1  4  1    8 )
    * (    1  3  2   )
    * (       2      )
    * (    8       2 ) */
   int n = 5;
   int ptr[]    = { 0,        2,             5,      7,7,   8 };
   int row[]    = { 0,   1,   1,   2,   4,   2,   3,   4   };
   double val[] = { 2.0, 1.0, 4.0, 1.0, 8.0, 3.0, 2.0, 2.0 };
   printf("Initial matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Perform symmetric scaling */
   spral_scaling_hungarian_default_options(&options);
   spral_scaling_hungarian_sym(n, ptr, row, val, scaling, match, &options, &inform);
   if(inform.flag<0) {
      printf("spral_scaling_hungarian_sym() returned with error %5d", inform.flag);
      exit(1);
   }

   /* Print scaling and matching */
   printf("Matching:");
   for(int i=0; i<n; i++) printf(" %10d", match[i]);
   printf("\nScaling: ");
   for(int i=0; i<n; i++) printf(" %10.2le", scaling[i]);
   printf("\n");

   /* Calculate scaled matrix and print it */
   for(int i=0; i<n; i++) {
      for(int j=ptr[i]; j<ptr[i+1]; j++)
         val[j] = scaling[row[j]] * val[j] * scaling[i];
   }
   printf("Scaled matrix:\n");
   spral_print_matrix(-1, SPRAL_MATRIX_REAL_SYM_INDEF, n, n, ptr, row, val, 0);

   /* Success */
   return 0;
}

The above code produces the following output:

Initial matrix:
Real symmetric indefinite matrix, dimension 5x5 with 8 entries.
0:   2.0000E+00   1.0000E+00
1:   1.0000E+00   4.0000E+00   1.0000E+00                8.0000E+00
2:                1.0000E+00   3.0000E+00   2.0000E+00
3:                             2.0000E+00
4:                8.0000E+00                             2.0000E+00
Matching:          0          4          3          2          1
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.
0:   1.0000E+00   2.5000E-01
1:   2.5000E-01   5.0000E-01   2.0412E-01                1.0000E+00
2:                2.0412E-01   1.0000E+00   1.0000E+00
3:                             1.0000E+00
4:                1.0000E+00                             2.5000E-01

Method

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.

\[\begin{split}\max_{\sigma} & \prod_{i=1}^m\prod_{j=1}^n |a_{ij}|\sigma_{ij} & \\ \mathrm{s.t.} & \sum_{i=1}^m\sigma_{ij} = 1, & \forall j=1,n \\ & \sum_{j=1}^n\sigma_{ij} = 1, & \forall i=1,m \\ & \sigma_{ij} \in \{0,1\}.\end{split}\]

The array \(\sigma\) gives a matching of rows to columns.

By using the transformation

\[w_{ij} = \log c_j - \log |a_{ij}|,\]

where \(c_j = \max_i |a_{ij}|\), the maximum product problem in \(a_{ij}\) is replaced by a minimum sum problem in \(w_{ij}\) where all entries are positive. By standard optimization theory, there exist dual variables \(u\) and \(v\) corresponding to the constraints that satisfy the first order optimality conditions

\[\begin{split}w_{ij} - u_i - v_j = 0, & \mbox{if } \sigma_{ij }=1, \\ w_{ij} - u_i - v_j \ge 0, & \mbox{if } \sigma_{ij }=0.\end{split}\]

To obtain a scaling we define scaling matrices \(D_r\) and \(D_c\) as

\[ \begin{align}\begin{aligned}d^r_i = e^{u_i},\\d^c_i = e^{v_i}.\end{aligned}\end{align} \]

If a symmetric scaling is required, we average these as

\[d_i = \sqrt{d^r_id^c_i}.\]

By the first order optimality conditions, these scaling matrices guarantee that

\[\begin{split}d^r_i|a_{ij}|d^c_j = 1, && \mbox{if } \sigma_{ij}=1, \\ d^r_i|a_{ij}|d^c_j \le 1, && \mbox{if } \sigma_{ij}=0.\end{split}\]

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:

[3]

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. [Journal] [Preprint]